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