deal.II version GIT relicensing-6960-g70e9c5b7be 2026-10-22 20:10:01+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
manifold_lib.cc
Go to the documentation of this file.
1// -----------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
4// Copyright (C) 2014 - 2026 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Detailed license information governing the source code and contributions
9// can be found in LICENSE.md and CONTRIBUTING.md at the top level directory.
10//
11// -----------------------------------------------------------------------------
12
13
14#include <deal.II/base/config.h>
15
17
18#ifdef DEAL_II_WITH_OPENCASCADE
19
20# include <BRepAdaptor_CompCurve.hxx>
21# include <BRepAdaptor_Curve.hxx>
22# if !DEAL_II_OPENCASCADE_VERSION_GTE(7, 6, 0)
23# include <BRepAdaptor_HCompCurve.hxx>
24# include <BRepAdaptor_HCurve.hxx>
25# endif
26# include <BRepTools.hxx>
27# include <BRep_Tool.hxx>
28# include <GCPnts_AbscissaPoint.hxx>
29# include <ShapeAnalysis_Curve.hxx>
30# include <ShapeAnalysis_Surface.hxx>
31# include <TopoDS.hxx>
32
33
34#endif
35
37
38#ifdef DEAL_II_WITH_OPENCASCADE
39
40
41namespace OpenCASCADE
42{
43 namespace
44 {
50# if DEAL_II_OPENCASCADE_VERSION_GTE(7, 6, 0)
51 opencascade::handle<Adaptor3d_Curve>
52 curve_adaptor(const TopoDS_Shape &shape)
53 {
54 Assert((shape.ShapeType() == TopAbs_WIRE) ||
55 (shape.ShapeType() == TopAbs_EDGE),
57 if (shape.ShapeType() == TopAbs_WIRE)
58 return opencascade::handle<BRepAdaptor_CompCurve>(
59 new BRepAdaptor_CompCurve(TopoDS::Wire(shape)));
60 else if (shape.ShapeType() == TopAbs_EDGE)
61 return opencascade::handle<BRepAdaptor_Curve>(
62 new BRepAdaptor_Curve(TopoDS::Edge(shape)));
63
65 return opencascade::handle<BRepAdaptor_Curve>(new BRepAdaptor_Curve());
66 }
67# else
68
69 opencascade::handle<Adaptor3d_HCurve>
70 curve_adaptor(const TopoDS_Shape &shape)
71 {
72 Assert((shape.ShapeType() == TopAbs_WIRE) ||
73 (shape.ShapeType() == TopAbs_EDGE),
75 if (shape.ShapeType() == TopAbs_WIRE)
76 return opencascade::handle<BRepAdaptor_HCompCurve>(
77 new BRepAdaptor_HCompCurve(TopoDS::Wire(shape)));
78
79
80 else if (shape.ShapeType() == TopAbs_EDGE)
81 return opencascade::handle<BRepAdaptor_HCurve>(
82 new BRepAdaptor_HCurve(TopoDS::Edge(shape)));
83
85 return opencascade::handle<BRepAdaptor_HCurve>(new BRepAdaptor_HCurve());
86 }
87# endif
88
89
90
91 // Helper internal functions.
92 double
93 shape_length(const TopoDS_Shape &sh)
94 {
95# if DEAL_II_OPENCASCADE_VERSION_GTE(7, 6, 0)
96 opencascade::handle<Adaptor3d_Curve> adapt = curve_adaptor(sh);
97 return GCPnts_AbscissaPoint::Length(*adapt);
98# else
99
100 opencascade::handle<Adaptor3d_HCurve> adapt = curve_adaptor(sh);
101
102 return GCPnts_AbscissaPoint::Length(adapt->GetCurve());
103# endif
104 }
105 } // namespace
106
107 /*======================= NormalProjectionManifold =========================*/
108 template <int dim, int spacedim>
110 const TopoDS_Shape &sh,
111 const double tolerance)
112 : sh(sh)
113 , tolerance(tolerance)
114 {
115 Assert(spacedim == 3, ExcNotImplemented());
116 }
117
118
119
120 template <int dim, int spacedim>
121 std::unique_ptr<Manifold<dim, spacedim>>
123 {
124 return std::unique_ptr<Manifold<dim, spacedim>>(
125 new NormalProjectionManifold(sh, tolerance));
126 }
127
128
129
130 template <int dim, int spacedim>
133 const ArrayView<const Point<spacedim>> &surrounding_points,
134 const Point<spacedim> &candidate) const
135 {
136 (void)surrounding_points;
137 if constexpr (running_in_debug_mode())
138 {
139 for (unsigned int i = 0; i < surrounding_points.size(); ++i)
140 Assert(closest_point(sh, surrounding_points[i], tolerance)
141 .distance(surrounding_points[i]) <
142 std::max(tolerance * surrounding_points[i].norm(),
143 tolerance),
144 ExcPointNotOnManifold<spacedim>(surrounding_points[i]));
145 }
146 return closest_point(sh, candidate, tolerance);
147 }
148
149
150 /*===================== DirectionalProjectionManifold ======================*/
151 template <int dim, int spacedim>
153 const TopoDS_Shape &sh,
154 const Tensor<1, spacedim> &direction,
155 const double tolerance)
156 : sh(sh)
157 , direction(direction)
158 , tolerance(tolerance)
159 {
160 Assert(spacedim == 3, ExcNotImplemented());
161 }
162
163
164
165 template <int dim, int spacedim>
166 std::unique_ptr<Manifold<dim, spacedim>>
168 {
169 return std::unique_ptr<Manifold<dim, spacedim>>(
170 new DirectionalProjectionManifold(sh, direction, tolerance));
171 }
172
173
174
175 template <int dim, int spacedim>
178 const ArrayView<const Point<spacedim>> &surrounding_points,
179 const Point<spacedim> &candidate) const
180 {
181 (void)surrounding_points;
182 if constexpr (running_in_debug_mode())
183 {
184 for (unsigned int i = 0; i < surrounding_points.size(); ++i)
185 Assert(closest_point(sh, surrounding_points[i], tolerance)
186 .distance(surrounding_points[i]) <
187 std::max(tolerance * surrounding_points[i].norm(),
188 tolerance),
189 ExcPointNotOnManifold<spacedim>(surrounding_points[i]));
190 }
191 return line_intersection(sh, candidate, direction, tolerance);
192 }
193
194
195
196 /*===================== NormalToMeshProjectionManifold =====================*/
197 template <int dim, int spacedim>
199 const TopoDS_Shape &sh,
200 const double tolerance)
201 : sh(sh)
202 , tolerance(tolerance)
203 {
204 Assert(spacedim == 3, ExcNotImplemented());
205 Assert(
206 std::get<0>(count_elements(sh)) > 0,
208 "NormalToMeshProjectionManifold needs a shape containing faces to operate."));
209 }
210
211 template <int dim, int spacedim>
212 std::unique_ptr<Manifold<dim, spacedim>>
214 {
215 return std::unique_ptr<Manifold<dim, spacedim>>(
217 }
218
219
220 namespace
221 {
222 template <int spacedim>
224 internal_project_to_manifold(const TopoDS_Shape &,
225 const double,
226 const ArrayView<const Point<spacedim>> &,
227 const Point<spacedim> &)
228 {
230 return {};
231 }
232
233 template <>
235 internal_project_to_manifold(
236 const TopoDS_Shape &sh,
237 const double tolerance,
238 const ArrayView<const Point<3>> &surrounding_points,
239 const Point<3> &candidate)
240 {
241 constexpr int spacedim = 3;
242 TopoDS_Shape out_shape;
243 Tensor<1, spacedim> average_normal;
244 if constexpr (running_in_debug_mode())
245 {
246 for (const auto &point : surrounding_points)
247 {
248 Assert(closest_point(sh, point, tolerance).distance(point) <
249 std::max(tolerance * point.norm(), tolerance),
250 ExcPointNotOnManifold<spacedim>(point));
251 }
252 }
253
254 switch (surrounding_points.size())
255 {
256 case 2:
257 {
258 for (const auto &point : surrounding_points)
259 {
260 std::tuple<Point<3>, Tensor<1, 3>, double, double>
261 p_and_diff_forms =
263 point,
264 tolerance);
265 average_normal += std::get<1>(p_and_diff_forms);
266 }
267
268 average_normal /= 2.0;
269
270 Assert(
271 average_normal.norm() > 1e-4,
273 "Failed to refine cell: the average of the surface normals at the surrounding edge turns out to be a null vector, making the projection direction undetermined."));
274
275 Tensor<1, 3> T = surrounding_points[0] - surrounding_points[1];
276 T /= T.norm();
277 average_normal = average_normal - (average_normal * T) * T;
278 average_normal /= average_normal.norm();
279 break;
280 }
281 case 4:
282 {
283 Tensor<1, 3> u = surrounding_points[1] - surrounding_points[0];
284 Tensor<1, 3> v = surrounding_points[2] - surrounding_points[0];
285 const double n1_coords[3] = {u[1] * v[2] - u[2] * v[1],
286 u[2] * v[0] - u[0] * v[2],
287 u[0] * v[1] - u[1] * v[0]};
288 Tensor<1, 3> n1(n1_coords);
289 n1 = n1 / n1.norm();
290 u = surrounding_points[2] - surrounding_points[3];
291 v = surrounding_points[1] - surrounding_points[3];
292 const double n2_coords[3] = {u[1] * v[2] - u[2] * v[1],
293 u[2] * v[0] - u[0] * v[2],
294 u[0] * v[1] - u[1] * v[0]};
295 Tensor<1, 3> n2(n2_coords);
296 n2 = n2 / n2.norm();
297
298 average_normal = (n1 + n2) / 2.0;
299
300 Assert(
301 average_normal.norm() > tolerance,
303 "Failed to refine cell: the normal estimated via the surrounding points turns out to be a null vector, making the projection direction undetermined."));
304
305 average_normal /= average_normal.norm();
306 break;
307 }
308 case 8:
309 {
310 Tensor<1, 3> u = surrounding_points[1] - surrounding_points[0];
311 Tensor<1, 3> v = surrounding_points[2] - surrounding_points[0];
312 const double n1_coords[3] = {u[1] * v[2] - u[2] * v[1],
313 u[2] * v[0] - u[0] * v[2],
314 u[0] * v[1] - u[1] * v[0]};
315 Tensor<1, 3> n1(n1_coords);
316 n1 = n1 / n1.norm();
317 u = surrounding_points[2] - surrounding_points[3];
318 v = surrounding_points[1] - surrounding_points[3];
319 const double n2_coords[3] = {u[1] * v[2] - u[2] * v[1],
320 u[2] * v[0] - u[0] * v[2],
321 u[0] * v[1] - u[1] * v[0]};
322 Tensor<1, 3> n2(n2_coords);
323 n2 = n2 / n2.norm();
324 u = surrounding_points[4] - surrounding_points[7];
325 v = surrounding_points[6] - surrounding_points[7];
326 const double n3_coords[3] = {u[1] * v[2] - u[2] * v[1],
327 u[2] * v[0] - u[0] * v[2],
328 u[0] * v[1] - u[1] * v[0]};
329 Tensor<1, 3> n3(n3_coords);
330 n3 = n3 / n3.norm();
331 u = surrounding_points[6] - surrounding_points[7];
332 v = surrounding_points[5] - surrounding_points[7];
333 const double n4_coords[3] = {u[1] * v[2] - u[2] * v[1],
334 u[2] * v[0] - u[0] * v[2],
335 u[0] * v[1] - u[1] * v[0]};
336 Tensor<1, 3> n4(n4_coords);
337 n4 = n4 / n4.norm();
338
339 average_normal = (n1 + n2 + n3 + n4) / 4.0;
340
341 Assert(
342 average_normal.norm() > tolerance,
344 "Failed to refine cell: the normal estimated via the surrounding points turns out to be a null vector, making the projection direction undetermined."));
345
346 average_normal /= average_normal.norm();
347 break;
348 }
349 default:
350 {
351 // Given an arbitrary number of points we compute all the possible
352 // normal vectors
353 for (unsigned int i = 0; i < surrounding_points.size(); ++i)
354 for (unsigned int j = 0; j < surrounding_points.size(); ++j)
355 if (j != i)
356 for (unsigned int k = 0; k < surrounding_points.size(); ++k)
357 if (k != j && k != i)
358 {
359 Tensor<1, 3> u =
360 surrounding_points[i] - surrounding_points[j];
361 Tensor<1, 3> v =
362 surrounding_points[i] - surrounding_points[k];
363 const double n_coords[3] = {u[1] * v[2] - u[2] * v[1],
364 u[2] * v[0] - u[0] * v[2],
365 u[0] * v[1] -
366 u[1] * v[0]};
367 Tensor<1, 3> n1(n_coords);
368 if (n1.norm() > tolerance)
369 {
370 n1 = n1 / n1.norm();
371 if (average_normal.norm() < tolerance)
372 average_normal = n1;
373 else
374 {
375 auto dot_prod = n1 * average_normal;
376 // We check that the direction of the normal
377 // vector w.r.t the current average, and make
378 // sure we flip it if it is opposite
379 if (dot_prod > 0)
380 average_normal += n1;
381 else
382 average_normal -= n1;
383 }
384 }
385 }
386 Assert(
387 average_normal.norm() > tolerance,
389 "Failed to compute a normal: the normal estimated via the surrounding points turns out to be a null vector, making the projection direction undetermined."));
390 average_normal = average_normal / average_normal.norm();
391 break;
392 }
393 }
394
395 return line_intersection(sh, candidate, average_normal, tolerance);
396 }
397 } // namespace
398
399
400 template <int dim, int spacedim>
403 const ArrayView<const Point<spacedim>> &surrounding_points,
404 const Point<spacedim> &candidate) const
405 {
406 return internal_project_to_manifold(sh,
407 tolerance,
408 surrounding_points,
409 candidate);
410 }
411
412
413 /*==================== ArclengthProjectionLineManifold =====================*/
414 template <int dim, int spacedim>
416 ArclengthProjectionLineManifold(const TopoDS_Shape &sh,
417 const double tolerance)
418 :
419
420 ChartManifold<dim, spacedim, 1>(sh.Closed() ? Point<1>(shape_length(sh)) :
421 Point<1>())
422 , sh(sh)
423 , curve(curve_adaptor(sh))
424 , tolerance(tolerance)
425 , length(shape_length(sh))
426 {
427 Assert(spacedim >= 2, ExcImpossibleInDimSpacedim(dim, spacedim));
428 }
429
430
431
432 template <int dim, int spacedim>
433 std::unique_ptr<Manifold<dim, spacedim>>
435 {
436 return std::unique_ptr<Manifold<dim, spacedim>>(
437 new ArclengthProjectionLineManifold(sh, tolerance));
438 }
439
440
441
442 template <int dim, int spacedim>
445 const Point<spacedim> &space_point) const
446 {
447 double t(0.0);
448 ShapeAnalysis_Curve curve_analysis;
449 gp_Pnt proj;
450
451 const double dist = curve_analysis.Project(
453 *curve, point(space_point), tolerance, proj, t, true);
454# else
455 curve->GetCurve(), point(space_point), tolerance, proj, t, true);
456# endif
457
458 (void)dist;
459 Assert(dist < tolerance * length,
460 ExcPointNotOnManifold<spacedim>(space_point));
461
462 return Point<1>(GCPnts_AbscissaPoint::Length(
464 *curve, curve->FirstParameter(), t));
465# else
466 curve->GetCurve(), curve->GetCurve().FirstParameter(), t));
467# endif
468 }
469
470
471
472 template <int dim, int spacedim>
475 const Point<1> &chart_point) const
476 {
477# if DEAL_II_OPENCASCADE_VERSION_GTE(7, 6, 0)
478 GCPnts_AbscissaPoint AP(*curve, chart_point[0], curve->FirstParameter());
479 gp_Pnt P = curve->Value(AP.Parameter());
480# else
481 GCPnts_AbscissaPoint AP(curve->GetCurve(),
482 chart_point[0],
483 curve->GetCurve().FirstParameter());
484 gp_Pnt P = curve->GetCurve().Value(AP.Parameter());
485# endif
486
487 return point<spacedim>(P);
488 }
489
490 template <int dim, int spacedim>
492 const double tolerance)
493 : face(face)
494 , tolerance(tolerance)
495 {}
496
497
498
499 template <int dim, int spacedim>
500 std::unique_ptr<Manifold<dim, spacedim>>
502 {
503 return std::unique_ptr<Manifold<dim, spacedim>>(
504 new NURBSPatchManifold<dim, spacedim>(face, tolerance));
505 }
506
507
508
509 template <int dim, int spacedim>
512 const Point<spacedim> &space_point) const
513 {
514 opencascade::handle<Geom_Surface> SurfToProj = BRep_Tool::Surface(face);
515
516 ShapeAnalysis_Surface projector(SurfToProj);
517 gp_Pnt2d proj_params = projector.ValueOfUV(point(space_point), tolerance);
518
519 double u = proj_params.X();
520 double v = proj_params.Y();
521
522 return {u, v};
523 }
524
525 template <int dim, int spacedim>
528 const Point<2> &chart_point) const
529 {
530 return ::OpenCASCADE::push_forward<spacedim>(face,
531 chart_point[0],
532 chart_point[1]);
533 }
534
535 template <int dim, int spacedim>
538 const Point<2> &chart_point) const
539 {
541
542 opencascade::handle<Geom_Surface> surf = BRep_Tool::Surface(face);
543
544 gp_Pnt q;
545 gp_Vec Du, Dv;
546 surf->D1(chart_point[0], chart_point[1], q, Du, Dv);
547
548 DX[0][0] = Du.X();
549 DX[1][0] = Du.Y();
550 if (spacedim > 2)
551 DX[2][0] = Du.Z();
552 else
553 Assert(std::abs(Du.Z()) < tolerance,
555 "Expecting derivative along Z to be zero! Bailing out."));
556 DX[0][1] = Dv.X();
557 DX[1][1] = Dv.Y();
558 if (spacedim > 2)
559 DX[2][1] = Dv.Z();
560 else
561 Assert(std::abs(Dv.Z()) < tolerance,
563 "Expecting derivative along Z to be zero! Bailing out."));
564 return DX;
565 }
566
567 template <int dim, int spacedim>
568 std::tuple<double, double, double, double>
570 {
571 Standard_Real umin, umax, vmin, vmax;
572 BRepTools::UVBounds(face, umin, umax, vmin, vmax);
573 return std::make_tuple(umin, umax, vmin, vmax);
574 }
575
576// We don't build the .inst file if deal.II isn't configured
577// with GMSH, but doxygen doesn't know that and tries to find that
578// file anyway for parsing -- which then of course it fails on. So
579// exclude the following from doxygen consideration.
580# ifndef DOXYGEN
581# include "opencascade/manifold_lib.inst"
582# endif
583} // end namespace OpenCASCADE
584
585
586#endif
587
virtual std::unique_ptr< Manifold< dim, spacedim > > clone() const override
virtual Point< spacedim > push_forward(const Point< 1 > &chart_point) const override
ArclengthProjectionLineManifold(const TopoDS_Shape &sh, const double tolerance=1e-7)
virtual Point< 1 > pull_back(const Point< spacedim > &space_point) const override
virtual Point< spacedim > project_to_manifold(const ArrayView< const Point< spacedim > > &surrounding_points, const Point< spacedim > &candidate) const override
DirectionalProjectionManifold(const TopoDS_Shape &sh, const Tensor< 1, spacedim > &direction, const double tolerance=1e-7)
virtual std::unique_ptr< Manifold< dim, spacedim > > clone() const override
virtual std::unique_ptr< Manifold< dim, spacedim > > clone() const override
std::tuple< double, double, double, double > get_uv_bounds() const
virtual Point< 2 > pull_back(const Point< spacedim > &space_point) const override
virtual Point< spacedim > push_forward(const Point< 2 > &chart_point) const override
NURBSPatchManifold(const TopoDS_Face &face, const double tolerance=1e-7)
virtual DerivativeForm< 1, 2, spacedim > push_forward_gradient(const Point< 2 > &chart_point) const override
virtual std::unique_ptr< Manifold< dim, spacedim > > clone() const override
virtual Point< spacedim > project_to_manifold(const ArrayView< const Point< spacedim > > &surrounding_points, const Point< spacedim > &candidate) const override
NormalProjectionManifold(const TopoDS_Shape &sh, const double tolerance=1e-7)
virtual Point< spacedim > project_to_manifold(const ArrayView< const Point< spacedim > > &surrounding_points, const Point< spacedim > &candidate) const override
virtual std::unique_ptr< Manifold< dim, spacedim > > clone() const override
NormalToMeshProjectionManifold(const TopoDS_Shape &sh, const double tolerance=1e-7)
Definition point.h:111
numbers::NumberTraits< Number >::real_type norm() const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:38
constexpr bool running_in_debug_mode()
Definition config.h:76
#define DEAL_II_OPENCASCADE_VERSION_GTE(major, minor, subminor)
Definition config.h:496
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:39
#define DEAL_II_ASSERT_UNREACHABLE()
#define DEAL_II_NOT_IMPLEMENTED()
static ::ExceptionBase & ExcImpossibleInDimSpacedim(int arg1, int arg2)
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcUnsupportedShape()
static ::ExceptionBase & ExcMessage(std::string arg1)
constexpr char T
std::tuple< unsigned int, unsigned int, unsigned int > count_elements(const TopoDS_Shape &shape)
Definition utilities.cc:90
std::tuple< Point< 3 >, Tensor< 1, 3 >, double, double > closest_point_and_differential_forms(const TopoDS_Shape &in_shape, const Point< 3 > &origin, const double tolerance=1e-7)
Definition utilities.cc:803
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:190
Point< dim > closest_point(const TopoDS_Shape &in_shape, const Point< dim > &origin, const double tolerance=1e-7)
Definition utilities.cc:793
Point< dim > line_intersection(const TopoDS_Shape &in_shape, const Point< dim > &origin, const Tensor< 1, dim > &direction, const double tolerance=1e-7)
Definition utilities.cc:518
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)