deal.II version GIT relicensing-2171-g98af52917d 2024-11-23 14:40: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# ifdef DEBUG
135 for (unsigned int i = 0; i < surrounding_points.size(); ++i)
136 Assert(closest_point(sh, surrounding_points[i], tolerance)
137 .distance(surrounding_points[i]) <
138 std::max(tolerance * surrounding_points[i].norm(), tolerance),
139 ExcPointNotOnManifold<spacedim>(surrounding_points[i]));
140# endif
141 return closest_point(sh, candidate, tolerance);
142 }
143
144
145 /*===================== DirectionalProjectionManifold ======================*/
146 template <int dim, int spacedim>
148 const TopoDS_Shape &sh,
149 const Tensor<1, spacedim> &direction,
150 const double tolerance)
151 : sh(sh)
152 , direction(direction)
153 , tolerance(tolerance)
154 {
155 Assert(spacedim == 3, ExcNotImplemented());
156 }
157
158
159
160 template <int dim, int spacedim>
161 std::unique_ptr<Manifold<dim, spacedim>>
163 {
164 return std::unique_ptr<Manifold<dim, spacedim>>(
165 new DirectionalProjectionManifold(sh, direction, tolerance));
166 }
167
168
169
170 template <int dim, int spacedim>
173 const ArrayView<const Point<spacedim>> &surrounding_points,
174 const Point<spacedim> &candidate) const
175 {
176 (void)surrounding_points;
177# ifdef DEBUG
178 for (unsigned int i = 0; i < surrounding_points.size(); ++i)
179 Assert(closest_point(sh, surrounding_points[i], tolerance)
180 .distance(surrounding_points[i]) <
181 std::max(tolerance * surrounding_points[i].norm(), tolerance),
182 ExcPointNotOnManifold<spacedim>(surrounding_points[i]));
183# endif
184 return line_intersection(sh, candidate, direction, tolerance);
185 }
186
187
188
189 /*===================== NormalToMeshProjectionManifold =====================*/
190 template <int dim, int spacedim>
192 const TopoDS_Shape &sh,
193 const double tolerance)
194 : sh(sh)
195 , tolerance(tolerance)
196 {
197 Assert(spacedim == 3, ExcNotImplemented());
198 Assert(
199 std::get<0>(count_elements(sh)) > 0,
201 "NormalToMeshProjectionManifold needs a shape containing faces to operate."));
202 }
203
204 template <int dim, int spacedim>
205 std::unique_ptr<Manifold<dim, spacedim>>
207 {
208 return std::unique_ptr<Manifold<dim, spacedim>>(
210 }
211
212
213 namespace
214 {
215 template <int spacedim>
217 internal_project_to_manifold(const TopoDS_Shape &,
218 const double,
219 const ArrayView<const Point<spacedim>> &,
220 const Point<spacedim> &)
221 {
223 return {};
224 }
225
226 template <>
228 internal_project_to_manifold(
229 const TopoDS_Shape &sh,
230 const double tolerance,
231 const ArrayView<const Point<3>> &surrounding_points,
232 const Point<3> &candidate)
233 {
234 constexpr int spacedim = 3;
235 TopoDS_Shape out_shape;
236 Tensor<1, spacedim> average_normal;
237# ifdef DEBUG
238 for (const auto &point : surrounding_points)
239 {
240 Assert(closest_point(sh, point, tolerance).distance(point) <
241 std::max(tolerance * point.norm(), tolerance),
242 ExcPointNotOnManifold<spacedim>(point));
243 }
244# endif
245
246 switch (surrounding_points.size())
247 {
248 case 2:
249 {
250 for (const auto &point : surrounding_points)
251 {
252 std::tuple<Point<3>, Tensor<1, 3>, double, double>
253 p_and_diff_forms =
255 point,
256 tolerance);
257 average_normal += std::get<1>(p_and_diff_forms);
258 }
259
260 average_normal /= 2.0;
261
262 Assert(
263 average_normal.norm() > 1e-4,
265 "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."));
266
267 Tensor<1, 3> T = surrounding_points[0] - surrounding_points[1];
268 T /= T.norm();
269 average_normal = average_normal - (average_normal * T) * T;
270 average_normal /= average_normal.norm();
271 break;
272 }
273 case 4:
274 {
275 Tensor<1, 3> u = surrounding_points[1] - surrounding_points[0];
276 Tensor<1, 3> v = surrounding_points[2] - surrounding_points[0];
277 const double n1_coords[3] = {u[1] * v[2] - u[2] * v[1],
278 u[2] * v[0] - u[0] * v[2],
279 u[0] * v[1] - u[1] * v[0]};
280 Tensor<1, 3> n1(n1_coords);
281 n1 = n1 / n1.norm();
282 u = surrounding_points[2] - surrounding_points[3];
283 v = surrounding_points[1] - surrounding_points[3];
284 const double n2_coords[3] = {u[1] * v[2] - u[2] * v[1],
285 u[2] * v[0] - u[0] * v[2],
286 u[0] * v[1] - u[1] * v[0]};
287 Tensor<1, 3> n2(n2_coords);
288 n2 = n2 / n2.norm();
289
290 average_normal = (n1 + n2) / 2.0;
291
292 Assert(
293 average_normal.norm() > tolerance,
295 "Failed to refine cell: the normal estimated via the surrounding points turns out to be a null vector, making the projection direction undetermined."));
296
297 average_normal /= average_normal.norm();
298 break;
299 }
300 case 8:
301 {
302 Tensor<1, 3> u = surrounding_points[1] - surrounding_points[0];
303 Tensor<1, 3> v = surrounding_points[2] - surrounding_points[0];
304 const double n1_coords[3] = {u[1] * v[2] - u[2] * v[1],
305 u[2] * v[0] - u[0] * v[2],
306 u[0] * v[1] - u[1] * v[0]};
307 Tensor<1, 3> n1(n1_coords);
308 n1 = n1 / n1.norm();
309 u = surrounding_points[2] - surrounding_points[3];
310 v = surrounding_points[1] - surrounding_points[3];
311 const double n2_coords[3] = {u[1] * v[2] - u[2] * v[1],
312 u[2] * v[0] - u[0] * v[2],
313 u[0] * v[1] - u[1] * v[0]};
314 Tensor<1, 3> n2(n2_coords);
315 n2 = n2 / n2.norm();
316 u = surrounding_points[4] - surrounding_points[7];
317 v = surrounding_points[6] - surrounding_points[7];
318 const double n3_coords[3] = {u[1] * v[2] - u[2] * v[1],
319 u[2] * v[0] - u[0] * v[2],
320 u[0] * v[1] - u[1] * v[0]};
321 Tensor<1, 3> n3(n3_coords);
322 n3 = n3 / n3.norm();
323 u = surrounding_points[6] - surrounding_points[7];
324 v = surrounding_points[5] - surrounding_points[7];
325 const double n4_coords[3] = {u[1] * v[2] - u[2] * v[1],
326 u[2] * v[0] - u[0] * v[2],
327 u[0] * v[1] - u[1] * v[0]};
328 Tensor<1, 3> n4(n4_coords);
329 n4 = n4 / n4.norm();
330
331 average_normal = (n1 + n2 + n3 + n4) / 4.0;
332
333 Assert(
334 average_normal.norm() > tolerance,
336 "Failed to refine cell: the normal estimated via the surrounding points turns out to be a null vector, making the projection direction undetermined."));
337
338 average_normal /= average_normal.norm();
339 break;
340 }
341 default:
342 {
343 // Given an arbitrary number of points we compute all the possible
344 // normal vectors
345 for (unsigned int i = 0; i < surrounding_points.size(); ++i)
346 for (unsigned int j = 0; j < surrounding_points.size(); ++j)
347 if (j != i)
348 for (unsigned int k = 0; k < surrounding_points.size(); ++k)
349 if (k != j && k != i)
350 {
351 Tensor<1, 3> u =
352 surrounding_points[i] - surrounding_points[j];
353 Tensor<1, 3> v =
354 surrounding_points[i] - surrounding_points[k];
355 const double n_coords[3] = {u[1] * v[2] - u[2] * v[1],
356 u[2] * v[0] - u[0] * v[2],
357 u[0] * v[1] -
358 u[1] * v[0]};
359 Tensor<1, 3> n1(n_coords);
360 if (n1.norm() > tolerance)
361 {
362 n1 = n1 / n1.norm();
363 if (average_normal.norm() < tolerance)
364 average_normal = n1;
365 else
366 {
367 auto dot_prod = n1 * average_normal;
368 // We check that the direction of the normal
369 // vector w.r.t the current average, and make
370 // sure we flip it if it is opposite
371 if (dot_prod > 0)
372 average_normal += n1;
373 else
374 average_normal -= n1;
375 }
376 }
377 }
378 Assert(
379 average_normal.norm() > tolerance,
381 "Failed to compute a normal: the normal estimated via the surrounding points turns out to be a null vector, making the projection direction undetermined."));
382 average_normal = average_normal / average_normal.norm();
383 break;
384 }
385 }
386
387 return line_intersection(sh, candidate, average_normal, tolerance);
388 }
389 } // namespace
390
391
392 template <int dim, int spacedim>
395 const ArrayView<const Point<spacedim>> &surrounding_points,
396 const Point<spacedim> &candidate) const
397 {
398 return internal_project_to_manifold(sh,
399 tolerance,
400 surrounding_points,
401 candidate);
402 }
403
404
405 /*==================== ArclengthProjectionLineManifold =====================*/
406 template <int dim, int spacedim>
408 ArclengthProjectionLineManifold(const TopoDS_Shape &sh,
409 const double tolerance)
410 :
411
412 ChartManifold<dim, spacedim, 1>(sh.Closed() ? Point<1>(shape_length(sh)) :
413 Point<1>())
414 , sh(sh)
415 , curve(curve_adaptor(sh))
416 , tolerance(tolerance)
417 , length(shape_length(sh))
418 {
419 Assert(spacedim >= 2, ExcImpossibleInDimSpacedim(dim, spacedim));
420 }
421
422
423
424 template <int dim, int spacedim>
425 std::unique_ptr<Manifold<dim, spacedim>>
427 {
428 return std::unique_ptr<Manifold<dim, spacedim>>(
429 new ArclengthProjectionLineManifold(sh, tolerance));
430 }
431
432
433
434 template <int dim, int spacedim>
437 const Point<spacedim> &space_point) const
438 {
439 double t(0.0);
440 ShapeAnalysis_Curve curve_analysis;
441 gp_Pnt proj;
442
443 const double dist = curve_analysis.Project(
445 *curve, point(space_point), tolerance, proj, t, true);
446# else
447 curve->GetCurve(), point(space_point), tolerance, proj, t, true);
448# endif
449
450 (void)dist;
451 Assert(dist < tolerance * length,
452 ExcPointNotOnManifold<spacedim>(space_point));
453
454 return Point<1>(GCPnts_AbscissaPoint::Length(
456 *curve, curve->FirstParameter(), t));
457# else
458 curve->GetCurve(), curve->GetCurve().FirstParameter(), t));
459# endif
460 }
461
462
463
464 template <int dim, int spacedim>
467 const Point<1> &chart_point) const
468 {
469# if DEAL_II_OPENCASCADE_VERSION_GTE(7, 6, 0)
470 GCPnts_AbscissaPoint AP(*curve, chart_point[0], curve->FirstParameter());
471 gp_Pnt P = curve->Value(AP.Parameter());
472# else
473 GCPnts_AbscissaPoint AP(curve->GetCurve(),
474 chart_point[0],
475 curve->GetCurve().FirstParameter());
476 gp_Pnt P = curve->GetCurve().Value(AP.Parameter());
477# endif
478
479 return point<spacedim>(P);
480 }
481
482 template <int dim, int spacedim>
484 const double tolerance)
485 : face(face)
486 , tolerance(tolerance)
487 {}
488
489
490
491 template <int dim, int spacedim>
492 std::unique_ptr<Manifold<dim, spacedim>>
494 {
495 return std::unique_ptr<Manifold<dim, spacedim>>(
496 new NURBSPatchManifold<dim, spacedim>(face, tolerance));
497 }
498
499
500
501 template <int dim, int spacedim>
504 const Point<spacedim> &space_point) const
505 {
506 Handle(Geom_Surface) SurfToProj = BRep_Tool::Surface(face);
507
508 ShapeAnalysis_Surface projector(SurfToProj);
509 gp_Pnt2d proj_params = projector.ValueOfUV(point(space_point), tolerance);
510
511 double u = proj_params.X();
512 double v = proj_params.Y();
513
514 return {u, v};
515 }
516
517 template <int dim, int spacedim>
520 const Point<2> &chart_point) const
521 {
522 return ::OpenCASCADE::push_forward<spacedim>(face,
523 chart_point[0],
524 chart_point[1]);
525 }
526
527 template <int dim, int spacedim>
530 const Point<2> &chart_point) const
531 {
533 Handle(Geom_Surface) surf = BRep_Tool::Surface(face);
534
535 gp_Pnt q;
536 gp_Vec Du, Dv;
537 surf->D1(chart_point[0], chart_point[1], q, Du, Dv);
538
539 DX[0][0] = Du.X();
540 DX[1][0] = Du.Y();
541 if (spacedim > 2)
542 DX[2][0] = Du.Z();
543 else
544 Assert(std::abs(Du.Z()) < tolerance,
546 "Expecting derivative along Z to be zero! Bailing out."));
547 DX[0][1] = Dv.X();
548 DX[1][1] = Dv.Y();
549 if (spacedim > 2)
550 DX[2][1] = Dv.Z();
551 else
552 Assert(std::abs(Dv.Z()) < tolerance,
554 "Expecting derivative along Z to be zero! Bailing out."));
555 return DX;
556 }
557
558 template <int dim, int spacedim>
559 std::tuple<double, double, double, double>
561 {
562 Standard_Real umin, umax, vmin, vmax;
563 BRepTools::UVBounds(face, umin, umax, vmin, vmax);
564 return std::make_tuple(umin, umax, vmin, vmax);
565 }
566
567// Explicit instantiations
568# include "manifold_lib.inst"
569} // end namespace OpenCASCADE
570
572
573#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:111
numbers::NumberTraits< Number >::real_type norm() const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:498
#define DEAL_II_OPENCASCADE_VERSION_GTE(major, minor, subminor)
Definition config.h:366
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
static ::ExceptionBase & ExcImpossibleInDimSpacedim(int arg1, int arg2)
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcUnsupportedShape()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define DEAL_II_ASSERT_UNREACHABLE()
#define DEAL_II_NOT_IMPLEMENTED()
std::tuple< unsigned int, unsigned int, unsigned int > count_elements(const TopoDS_Shape &shape)
Definition utilities.cc:91
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:798
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:191
Point< dim > closest_point(const TopoDS_Shape &in_shape, const Point< dim > &origin, const double tolerance=1e-7)
Definition utilities.cc:788
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:513
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)