Reference documentation for deal.II version 9.5.0
\(\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// Copyright (C) 2014 - 2022 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
16
17#include <deal.II/base/config.h>
18
20
21#ifdef DEAL_II_WITH_OPENCASCADE
22
23
24# include <BRepAdaptor_CompCurve.hxx>
25# include <BRepAdaptor_Curve.hxx>
26# if !DEAL_II_OPENCASCADE_VERSION_GTE(7, 6, 0)
27# include <BRepAdaptor_HCompCurve.hxx>
28# include <BRepAdaptor_HCurve.hxx>
29# endif
30# include <BRepTools.hxx>
31# include <BRep_Tool.hxx>
32# include <GCPnts_AbscissaPoint.hxx>
33# include <ShapeAnalysis_Curve.hxx>
34# include <ShapeAnalysis_Surface.hxx>
35# include <TopoDS.hxx>
36# if !DEAL_II_OPENCASCADE_VERSION_GTE(7, 0, 0)
37# include <Handle_Adaptor3d_HCurve.hxx>
38# endif
39
40
42
43
44namespace OpenCASCADE
45{
46 namespace
47 {
53# if DEAL_II_OPENCASCADE_VERSION_GTE(7, 6, 0)
54 Handle_Adaptor3d_Curve
55 curve_adaptor(const TopoDS_Shape &shape)
56 {
57 Assert((shape.ShapeType() == TopAbs_WIRE) ||
58 (shape.ShapeType() == TopAbs_EDGE),
60 if (shape.ShapeType() == TopAbs_WIRE)
61 return Handle(BRepAdaptor_CompCurve)(
62 new BRepAdaptor_CompCurve(TopoDS::Wire(shape)));
63 else if (shape.ShapeType() == TopAbs_EDGE)
64 return Handle(BRepAdaptor_Curve)(
65 new BRepAdaptor_Curve(TopoDS::Edge(shape)));
66
67 Assert(false, ExcInternalError());
68 return Handle(BRepAdaptor_Curve)(new BRepAdaptor_Curve());
69 }
70# else
71 Handle_Adaptor3d_HCurve
72 curve_adaptor(const TopoDS_Shape &shape)
73 {
74 Assert((shape.ShapeType() == TopAbs_WIRE) ||
75 (shape.ShapeType() == TopAbs_EDGE),
77 if (shape.ShapeType() == TopAbs_WIRE)
78 return Handle(BRepAdaptor_HCompCurve)(
79 new BRepAdaptor_HCompCurve(TopoDS::Wire(shape)));
80 else if (shape.ShapeType() == TopAbs_EDGE)
81 return Handle(BRepAdaptor_HCurve)(
82 new BRepAdaptor_HCurve(TopoDS::Edge(shape)));
83
84 Assert(false, ExcInternalError());
85 return 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 Handle_Adaptor3d_Curve adapt = curve_adaptor(sh);
97 return GCPnts_AbscissaPoint::Length(*adapt);
98# else
99 Handle_Adaptor3d_HCurve adapt = curve_adaptor(sh);
100 return GCPnts_AbscissaPoint::Length(adapt->GetCurve());
101# endif
102 }
103 } // namespace
104
105 /*======================= NormalProjectionManifold =========================*/
106 template <int dim, int spacedim>
108 const TopoDS_Shape &sh,
109 const double tolerance)
110 : sh(sh)
111 , tolerance(tolerance)
112 {
113 Assert(spacedim == 3, ExcNotImplemented());
114 }
115
116
117
118 template <int dim, int spacedim>
119 std::unique_ptr<Manifold<dim, spacedim>>
121 {
122 return std::unique_ptr<Manifold<dim, spacedim>>(
123 new NormalProjectionManifold(sh, tolerance));
124 }
125
126
127
128 template <int dim, int spacedim>
131 const ArrayView<const Point<spacedim>> &surrounding_points,
132 const Point<spacedim> & candidate) const
133 {
134 (void)surrounding_points;
135# ifdef DEBUG
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(), tolerance),
140 ExcPointNotOnManifold<spacedim>(surrounding_points[i]));
141# endif
142 return closest_point(sh, candidate, tolerance);
143 }
144
145
146 /*===================== DirectionalProjectionManifold ======================*/
147 template <int dim, int spacedim>
149 const TopoDS_Shape & sh,
150 const Tensor<1, spacedim> &direction,
151 const double tolerance)
152 : sh(sh)
153 , direction(direction)
154 , tolerance(tolerance)
155 {
156 Assert(spacedim == 3, ExcNotImplemented());
157 }
158
159
160
161 template <int dim, int spacedim>
162 std::unique_ptr<Manifold<dim, spacedim>>
164 {
165 return std::unique_ptr<Manifold<dim, spacedim>>(
166 new DirectionalProjectionManifold(sh, direction, tolerance));
167 }
168
169
170
171 template <int dim, int spacedim>
174 const ArrayView<const Point<spacedim>> &surrounding_points,
175 const Point<spacedim> & candidate) const
176 {
177 (void)surrounding_points;
178# ifdef DEBUG
179 for (unsigned int i = 0; i < surrounding_points.size(); ++i)
180 Assert(closest_point(sh, surrounding_points[i], tolerance)
181 .distance(surrounding_points[i]) <
182 std::max(tolerance * surrounding_points[i].norm(), tolerance),
183 ExcPointNotOnManifold<spacedim>(surrounding_points[i]));
184# endif
185 return line_intersection(sh, candidate, direction, tolerance);
186 }
187
188
189
190 /*===================== NormalToMeshProjectionManifold =====================*/
191 template <int dim, int spacedim>
193 const TopoDS_Shape &sh,
194 const double tolerance)
195 : sh(sh)
196 , tolerance(tolerance)
197 {
198 Assert(spacedim == 3, ExcNotImplemented());
199 Assert(
200 std::get<0>(count_elements(sh)) > 0,
202 "NormalToMeshProjectionManifold needs a shape containing faces to operate."));
203 }
204
205 template <int dim, int spacedim>
206 std::unique_ptr<Manifold<dim, spacedim>>
208 {
209 return std::unique_ptr<Manifold<dim, spacedim>>(
211 }
212
213
214 namespace
215 {
216 template <int spacedim>
218 internal_project_to_manifold(const TopoDS_Shape &,
219 const double,
220 const ArrayView<const Point<spacedim>> &,
221 const Point<spacedim> &)
222 {
223 Assert(false, ExcNotImplemented());
224 return {};
225 }
226
227 template <>
229 internal_project_to_manifold(
230 const TopoDS_Shape & sh,
231 const double tolerance,
232 const ArrayView<const Point<3>> &surrounding_points,
233 const Point<3> & candidate)
234 {
235 constexpr int spacedim = 3;
236 TopoDS_Shape out_shape;
237 Tensor<1, spacedim> average_normal;
238# ifdef DEBUG
239 for (const auto &point : surrounding_points)
240 {
241 Assert(closest_point(sh, point, tolerance).distance(point) <
242 std::max(tolerance * point.norm(), tolerance),
243 ExcPointNotOnManifold<spacedim>(point));
244 }
245# endif
246
247 switch (surrounding_points.size())
248 {
249 case 2:
250 {
251 for (const auto &point : surrounding_points)
252 {
253 std::tuple<Point<3>, Tensor<1, 3>, double, double>
254 p_and_diff_forms =
256 point,
257 tolerance);
258 average_normal += std::get<1>(p_and_diff_forms);
259 }
260
261 average_normal /= 2.0;
262
263 Assert(
264 average_normal.norm() > 1e-4,
266 "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."));
267
268 Tensor<1, 3> T = surrounding_points[0] - surrounding_points[1];
269 T /= T.norm();
270 average_normal = average_normal - (average_normal * T) * T;
271 average_normal /= average_normal.norm();
272 break;
273 }
274 case 4:
275 {
276 Tensor<1, 3> u = surrounding_points[1] - surrounding_points[0];
277 Tensor<1, 3> v = surrounding_points[2] - surrounding_points[0];
278 const double n1_coords[3] = {u[1] * v[2] - u[2] * v[1],
279 u[2] * v[0] - u[0] * v[2],
280 u[0] * v[1] - u[1] * v[0]};
281 Tensor<1, 3> n1(n1_coords);
282 n1 = n1 / n1.norm();
283 u = surrounding_points[2] - surrounding_points[3];
284 v = surrounding_points[1] - surrounding_points[3];
285 const double n2_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> n2(n2_coords);
289 n2 = n2 / n2.norm();
290
291 average_normal = (n1 + n2) / 2.0;
292
293 Assert(
294 average_normal.norm() > tolerance,
296 "Failed to refine cell: the normal estimated via the surrounding points turns out to be a null vector, making the projection direction undetermined."));
297
298 average_normal /= average_normal.norm();
299 break;
300 }
301 case 8:
302 {
303 Tensor<1, 3> u = surrounding_points[1] - surrounding_points[0];
304 Tensor<1, 3> v = surrounding_points[2] - surrounding_points[0];
305 const double n1_coords[3] = {u[1] * v[2] - u[2] * v[1],
306 u[2] * v[0] - u[0] * v[2],
307 u[0] * v[1] - u[1] * v[0]};
308 Tensor<1, 3> n1(n1_coords);
309 n1 = n1 / n1.norm();
310 u = surrounding_points[2] - surrounding_points[3];
311 v = surrounding_points[1] - surrounding_points[3];
312 const double n2_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> n2(n2_coords);
316 n2 = n2 / n2.norm();
317 u = surrounding_points[4] - surrounding_points[7];
318 v = surrounding_points[6] - surrounding_points[7];
319 const double n3_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> n3(n3_coords);
323 n3 = n3 / n3.norm();
324 u = surrounding_points[6] - surrounding_points[7];
325 v = surrounding_points[5] - surrounding_points[7];
326 const double n4_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> n4(n4_coords);
330 n4 = n4 / n4.norm();
331
332 average_normal = (n1 + n2 + n3 + n4) / 4.0;
333
334 Assert(
335 average_normal.norm() > tolerance,
337 "Failed to refine cell: the normal estimated via the surrounding points turns out to be a null vector, making the projection direction undetermined."));
338
339 average_normal /= average_normal.norm();
340 break;
341 }
342 default:
343 {
344 // Given an arbitrary number of points we compute all the possible
345 // normal vectors
346 for (unsigned int i = 0; i < surrounding_points.size(); ++i)
347 for (unsigned int j = 0; j < surrounding_points.size(); ++j)
348 if (j != i)
349 for (unsigned int k = 0; k < surrounding_points.size(); ++k)
350 if (k != j && k != i)
351 {
352 Tensor<1, 3> u =
353 surrounding_points[i] - surrounding_points[j];
354 Tensor<1, 3> v =
355 surrounding_points[i] - surrounding_points[k];
356 const double n_coords[3] = {u[1] * v[2] - u[2] * v[1],
357 u[2] * v[0] - u[0] * v[2],
358 u[0] * v[1] -
359 u[1] * v[0]};
360 Tensor<1, 3> n1(n_coords);
361 if (n1.norm() > tolerance)
362 {
363 n1 = n1 / n1.norm();
364 if (average_normal.norm() < tolerance)
365 average_normal = n1;
366 else
367 {
368 auto dot_prod = n1 * average_normal;
369 // We check that the direction of the normal
370 // vector w.r.t the current average, and make
371 // sure we flip it if it is opposite
372 if (dot_prod > 0)
373 average_normal += n1;
374 else
375 average_normal -= n1;
376 }
377 }
378 }
379 Assert(
380 average_normal.norm() > tolerance,
382 "Failed to compute a normal: the normal estimated via the surrounding points turns out to be a null vector, making the projection direction undetermined."));
383 average_normal = average_normal / average_normal.norm();
384 break;
385 }
386 }
387
388 return line_intersection(sh, candidate, average_normal, tolerance);
389 }
390 } // namespace
391
392
393 template <int dim, int spacedim>
396 const ArrayView<const Point<spacedim>> &surrounding_points,
397 const Point<spacedim> & candidate) const
398 {
399 return internal_project_to_manifold(sh,
400 tolerance,
401 surrounding_points,
402 candidate);
403 }
404
405
406 /*==================== ArclengthProjectionLineManifold =====================*/
407 template <int dim, int spacedim>
409 ArclengthProjectionLineManifold(const TopoDS_Shape &sh,
410 const double tolerance)
411 :
412
413 ChartManifold<dim, spacedim, 1>(sh.Closed() ? Point<1>(shape_length(sh)) :
414 Point<1>())
415 , sh(sh)
416 , curve(curve_adaptor(sh))
417 , tolerance(tolerance)
418 , length(shape_length(sh))
419 {
420 Assert(spacedim >= 2, ExcImpossibleInDimSpacedim(dim, spacedim));
421 }
422
423
424
425 template <int dim, int spacedim>
426 std::unique_ptr<Manifold<dim, spacedim>>
428 {
429 return std::unique_ptr<Manifold<dim, spacedim>>(
430 new ArclengthProjectionLineManifold(sh, tolerance));
431 }
432
433
434
435 template <int dim, int spacedim>
438 const Point<spacedim> &space_point) const
439 {
440 double t(0.0);
441 ShapeAnalysis_Curve curve_analysis;
442 gp_Pnt proj;
443
444 const double dist = curve_analysis.Project(
446 *curve, point(space_point), tolerance, proj, t, true);
447# else
448 curve->GetCurve(), point(space_point), tolerance, proj, t, true);
449# endif
450
451 (void)dist;
452 Assert(dist < tolerance * length,
453 ExcPointNotOnManifold<spacedim>(space_point));
454
455 return Point<1>(GCPnts_AbscissaPoint::Length(
457 *curve, curve->FirstParameter(), t));
458# else
459 curve->GetCurve(), curve->GetCurve().FirstParameter(), t));
460# endif
461 }
462
463
464
465 template <int dim, int spacedim>
468 const Point<1> &chart_point) const
469 {
470# if DEAL_II_OPENCASCADE_VERSION_GTE(7, 6, 0)
471 GCPnts_AbscissaPoint AP(*curve, chart_point[0], curve->FirstParameter());
472 gp_Pnt P = curve->Value(AP.Parameter());
473# else
474 GCPnts_AbscissaPoint AP(curve->GetCurve(),
475 chart_point[0],
476 curve->GetCurve().FirstParameter());
477 gp_Pnt P = curve->GetCurve().Value(AP.Parameter());
478# endif
479
480 return point<spacedim>(P);
481 }
482
483 template <int dim, int spacedim>
485 const double tolerance)
486 : face(face)
487 , tolerance(tolerance)
488 {}
489
490
491
492 template <int dim, int spacedim>
493 std::unique_ptr<Manifold<dim, spacedim>>
495 {
496 return std::unique_ptr<Manifold<dim, spacedim>>(
497 new NURBSPatchManifold<dim, spacedim>(face, tolerance));
498 }
499
500
501
502 template <int dim, int spacedim>
505 const Point<spacedim> &space_point) const
506 {
507 Handle(Geom_Surface) SurfToProj = BRep_Tool::Surface(face);
508
509 ShapeAnalysis_Surface projector(SurfToProj);
510 gp_Pnt2d proj_params = projector.ValueOfUV(point(space_point), tolerance);
511
512 double u = proj_params.X();
513 double v = proj_params.Y();
514
515 return {u, v};
516 }
517
518 template <int dim, int spacedim>
521 const Point<2> &chart_point) const
522 {
523 return ::OpenCASCADE::push_forward<spacedim>(face,
524 chart_point[0],
525 chart_point[1]);
526 }
527
528 template <int dim, int spacedim>
531 const Point<2> &chart_point) const
532 {
534 Handle(Geom_Surface) surf = BRep_Tool::Surface(face);
535
536 gp_Pnt q;
537 gp_Vec Du, Dv;
538 surf->D1(chart_point[0], chart_point[1], q, Du, Dv);
539
540 DX[0][0] = Du.X();
541 DX[1][0] = Du.Y();
542 if (spacedim > 2)
543 DX[2][0] = Du.Z();
544 else
545 Assert(std::abs(Du.Z()) < tolerance,
547 "Expecting derivative along Z to be zero! Bailing out."));
548 DX[0][1] = Dv.X();
549 DX[1][1] = Dv.Y();
550 if (spacedim > 2)
551 DX[2][1] = Dv.Z();
552 else
553 Assert(std::abs(Dv.Z()) < tolerance,
555 "Expecting derivative along Z to be zero! Bailing out."));
556 return DX;
557 }
558
559 template <int dim, int spacedim>
560 std::tuple<double, double, double, double>
562 {
563 Standard_Real umin, umax, vmin, vmax;
564 BRepTools::UVBounds(face, umin, umax, vmin, vmax);
565 return std::make_tuple(umin, umax, vmin, vmax);
566 }
567
568// Explicit instantiations
569# include "manifold_lib.inst"
570} // end namespace OpenCASCADE
571
573
574#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:112
numbers::NumberTraits< Number >::real_type norm() const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_OPENCASCADE_VERSION_GTE(major, minor, subminor)
Definition config.h:323
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
static ::ExceptionBase & ExcImpossibleInDimSpacedim(int arg1, int arg2)
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcUnsupportedShape()
static ::ExceptionBase & ExcMessage(std::string arg1)
std::tuple< unsigned int, unsigned int, unsigned int > count_elements(const TopoDS_Shape &shape)
Definition utilities.cc:92
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:801
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:189
Point< dim > closest_point(const TopoDS_Shape &in_shape, const Point< dim > &origin, const double tolerance=1e-7)
Definition utilities.cc:791
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:511
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)