Loading [MathJax]/extensions/TeX/newcommand.js
 Reference documentation for deal.II version 9.3.3
\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\}}
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
manifold_lib.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2014 - 2020 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
18
19#ifdef DEAL_II_WITH_OPENCASCADE
20
21
22# include <BRepAdaptor_CompCurve.hxx>
23# include <BRepAdaptor_Curve.hxx>
24# include <BRepAdaptor_HCompCurve.hxx>
25# include <BRepAdaptor_HCurve.hxx>
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 <Standard_Version.hxx>
32# include <TopoDS.hxx>
33# if (OCC_VERSION_MAJOR < 7)
34# include <Handle_Adaptor3d_HCurve.hxx>
35# endif
36
37
39
40
41namespace OpenCASCADE
42{
43 namespace
44 {
50 Handle_Adaptor3d_HCurve
51 curve_adaptor(const TopoDS_Shape &shape)
52 {
53 Assert((shape.ShapeType() == TopAbs_WIRE) ||
54 (shape.ShapeType() == TopAbs_EDGE),
56 if (shape.ShapeType() == TopAbs_WIRE)
57 return Handle(BRepAdaptor_HCompCurve)(
58 new BRepAdaptor_HCompCurve(TopoDS::Wire(shape)));
59 else if (shape.ShapeType() == TopAbs_EDGE)
60 return Handle(BRepAdaptor_HCurve)(
61 new BRepAdaptor_HCurve(TopoDS::Edge(shape)));
62
63 Assert(false, ExcInternalError());
64 return Handle(BRepAdaptor_HCurve)(new BRepAdaptor_HCurve());
65 }
66
67
68
69 // Helper internal functions.
70 double
71 shape_length(const TopoDS_Shape &sh)
72 {
73 Handle_Adaptor3d_HCurve adapt = curve_adaptor(sh);
74 return GCPnts_AbscissaPoint::Length(adapt->GetCurve());
75 }
76 } // namespace
77
78 /*======================= NormalProjectionManifold =========================*/
79 template <int dim, int spacedim>
81 const TopoDS_Shape &sh,
82 const double tolerance)
83 : sh(sh)
84 , tolerance(tolerance)
85 {
86 Assert(spacedim == 3, ExcNotImplemented());
87 }
88
89
90
91 template <int dim, int spacedim>
92 std::unique_ptr<Manifold<dim, spacedim>>
94 {
95 return std::unique_ptr<Manifold<dim, spacedim>>(
96 new NormalProjectionManifold(sh, tolerance));
97 }
98
99
100
101 template <int dim, int spacedim>
104 const ArrayView<const Point<spacedim>> &surrounding_points,
105 const Point<spacedim> & candidate) const
106 {
107 (void)surrounding_points;
108# ifdef DEBUG
109 for (unsigned int i = 0; i < surrounding_points.size(); ++i)
110 Assert(closest_point(sh, surrounding_points[i], tolerance)
111 .distance(surrounding_points[i]) <
112 std::max(tolerance * surrounding_points[i].norm(), tolerance),
113 ExcPointNotOnManifold<spacedim>(surrounding_points[i]));
114# endif
115 return closest_point(sh, candidate, tolerance);
116 }
117
118
119 /*===================== DirectionalProjectionManifold ======================*/
120 template <int dim, int spacedim>
122 const TopoDS_Shape & sh,
123 const Tensor<1, spacedim> &direction,
124 const double tolerance)
125 : sh(sh)
126 , direction(direction)
127 , tolerance(tolerance)
128 {
129 Assert(spacedim == 3, ExcNotImplemented());
130 }
131
132
133
134 template <int dim, int spacedim>
135 std::unique_ptr<Manifold<dim, spacedim>>
137 {
138 return std::unique_ptr<Manifold<dim, spacedim>>(
139 new DirectionalProjectionManifold(sh, direction, tolerance));
140 }
141
142
143
144 template <int dim, int spacedim>
147 const ArrayView<const Point<spacedim>> &surrounding_points,
148 const Point<spacedim> & candidate) const
149 {
150 (void)surrounding_points;
151# ifdef DEBUG
152 for (unsigned int i = 0; i < surrounding_points.size(); ++i)
153 Assert(closest_point(sh, surrounding_points[i], tolerance)
154 .distance(surrounding_points[i]) <
155 std::max(tolerance * surrounding_points[i].norm(), tolerance),
156 ExcPointNotOnManifold<spacedim>(surrounding_points[i]));
157# endif
158 return line_intersection(sh, candidate, direction, tolerance);
159 }
160
161
162
163 /*===================== NormalToMeshProjectionManifold =====================*/
164 template <int dim, int spacedim>
166 const TopoDS_Shape &sh,
167 const double tolerance)
168 : sh(sh)
169 , tolerance(tolerance)
170 {
171 Assert(spacedim == 3, ExcNotImplemented());
172 Assert(
173 std::get<0>(count_elements(sh)) > 0,
175 "NormalToMeshProjectionManifold needs a shape containing faces to operate."));
176 }
177
178 template <int dim, int spacedim>
179 std::unique_ptr<Manifold<dim, spacedim>>
181 {
182 return std::unique_ptr<Manifold<dim, spacedim>>(
184 }
185
186
187 namespace
188 {
189 template <int spacedim>
191 internal_project_to_manifold(const TopoDS_Shape &,
192 const double,
193 const ArrayView<const Point<spacedim>> &,
194 const Point<spacedim> &)
195 {
196 Assert(false, ExcNotImplemented());
197 return {};
198 }
199
200 template <>
202 internal_project_to_manifold(
203 const TopoDS_Shape & sh,
204 const double tolerance,
205 const ArrayView<const Point<3>> &surrounding_points,
206 const Point<3> & candidate)
207 {
208 constexpr int spacedim = 3;
209 TopoDS_Shape out_shape;
210 Tensor<1, spacedim> average_normal;
211# ifdef DEBUG
212 for (unsigned int i = 0; i < surrounding_points.size(); ++i)
213 {
214 Assert(closest_point(sh, surrounding_points[i], tolerance)
215 .distance(surrounding_points[i]) <
216 std::max(tolerance * surrounding_points[i].norm(),
217 tolerance),
218 ExcPointNotOnManifold<spacedim>(surrounding_points[i]));
219 }
220# endif
221
222 switch (surrounding_points.size())
223 {
224 case 2:
225 {
226 for (unsigned int i = 0; i < surrounding_points.size(); ++i)
227 {
228 std::tuple<Point<3>, Tensor<1, 3>, double, double>
229 p_and_diff_forms = closest_point_and_differential_forms(
230 sh, surrounding_points[i], tolerance);
231 average_normal += std::get<1>(p_and_diff_forms);
232 }
233
234 average_normal /= 2.0;
235
236 Assert(
237 average_normal.norm() > 1e-4,
239 "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."));
240
241 Tensor<1, 3> T = surrounding_points[0] - surrounding_points[1];
242 T /= T.norm();
243 average_normal = average_normal - (average_normal * T) * T;
244 average_normal /= average_normal.norm();
245 break;
246 }
247 case 4:
248 {
249 Tensor<1, 3> u = surrounding_points[1] - surrounding_points[0];
250 Tensor<1, 3> v = surrounding_points[2] - surrounding_points[0];
251 const double n1_coords[3] = {u[1] * v[2] - u[2] * v[1],
252 u[2] * v[0] - u[0] * v[2],
253 u[0] * v[1] - u[1] * v[0]};
254 Tensor<1, 3> n1(n1_coords);
255 n1 = n1 / n1.norm();
256 u = surrounding_points[2] - surrounding_points[3];
257 v = surrounding_points[1] - surrounding_points[3];
258 const double n2_coords[3] = {u[1] * v[2] - u[2] * v[1],
259 u[2] * v[0] - u[0] * v[2],
260 u[0] * v[1] - u[1] * v[0]};
261 Tensor<1, 3> n2(n2_coords);
262 n2 = n2 / n2.norm();
263
264 average_normal = (n1 + n2) / 2.0;
265
266 Assert(
267 average_normal.norm() > tolerance,
269 "Failed to refine cell: the normal estimated via the surrounding points turns out to be a null vector, making the projection direction undetermined."));
270
271 average_normal /= average_normal.norm();
272 break;
273 }
274 case 8:
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 u = surrounding_points[4] - surrounding_points[7];
291 v = surrounding_points[6] - surrounding_points[7];
292 const double n3_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> n3(n3_coords);
296 n3 = n3 / n3.norm();
297 u = surrounding_points[6] - surrounding_points[7];
298 v = surrounding_points[5] - surrounding_points[7];
299 const double n4_coords[3] = {u[1] * v[2] - u[2] * v[1],
300 u[2] * v[0] - u[0] * v[2],
301 u[0] * v[1] - u[1] * v[0]};
302 Tensor<1, 3> n4(n4_coords);
303 n4 = n4 / n4.norm();
304
305 average_normal = (n1 + n2 + n3 + n4) / 4.0;
306
307 Assert(
308 average_normal.norm() > tolerance,
310 "Failed to refine cell: the normal estimated via the surrounding points turns out to be a null vector, making the projection direction undetermined."));
311
312 average_normal /= average_normal.norm();
313 break;
314 }
315 default:
316 {
317 // Given an arbitrary number of points we compute all the possible
318 // normal vectors
319 for (unsigned int i = 0; i < surrounding_points.size(); ++i)
320 for (unsigned int j = 0; j < surrounding_points.size(); ++j)
321 if (j != i)
322 for (unsigned int k = 0; k < surrounding_points.size(); ++k)
323 if (k != j && k != i)
324 {
325 Tensor<1, 3> u =
326 surrounding_points[i] - surrounding_points[j];
327 Tensor<1, 3> v =
328 surrounding_points[i] - surrounding_points[k];
329 const double n_coords[3] = {u[1] * v[2] - u[2] * v[1],
330 u[2] * v[0] - u[0] * v[2],
331 u[0] * v[1] -
332 u[1] * v[0]};
333 Tensor<1, 3> n1(n_coords);
334 if (n1.norm() > tolerance)
335 {
336 n1 = n1 / n1.norm();
337 if (average_normal.norm() < tolerance)
338 average_normal = n1;
339 else
340 {
341 auto dot_prod = n1 * average_normal;
342 // We check that the direction of the normal
343 // vector w.r.t the current average, and make
344 // sure we flip it if it is opposite
345 if (dot_prod > 0)
346 average_normal += n1;
347 else
348 average_normal -= n1;
349 }
350 }
351 }
352 Assert(
353 average_normal.norm() > tolerance,
355 "Failed to compute a normal: the normal estimated via the surrounding points turns out to be a null vector, making the projection direction undetermined."));
356 average_normal = average_normal / average_normal.norm();
357 break;
358 }
359 }
360
361 return line_intersection(sh, candidate, average_normal, tolerance);
362 }
363 } // namespace
364
365
366 template <int dim, int spacedim>
369 const ArrayView<const Point<spacedim>> &surrounding_points,
370 const Point<spacedim> & candidate) const
371 {
372 return internal_project_to_manifold(sh,
373 tolerance,
374 surrounding_points,
375 candidate);
376 }
377
378
379 /*==================== ArclengthProjectionLineManifold =====================*/
380 template <int dim, int spacedim>
382 ArclengthProjectionLineManifold(const TopoDS_Shape &sh,
383 const double tolerance)
384 :
385
386 ChartManifold<dim, spacedim, 1>(sh.Closed() ? Point<1>(shape_length(sh)) :
387 Point<1>())
388 , sh(sh)
389 , curve(curve_adaptor(sh))
390 , tolerance(tolerance)
391 , length(shape_length(sh))
392 {
393 Assert(spacedim >= 2, ExcImpossibleInDimSpacedim(dim, spacedim));
394 }
395
396
397
398 template <int dim, int spacedim>
399 std::unique_ptr<Manifold<dim, spacedim>>
401 {
402 return std::unique_ptr<Manifold<dim, spacedim>>(
403 new ArclengthProjectionLineManifold(sh, tolerance));
404 }
405
406
407
408 template <int dim, int spacedim>
411 const Point<spacedim> &space_point) const
412 {
413 double t(0.0);
414 ShapeAnalysis_Curve curve_analysis;
415 gp_Pnt proj;
416 const double dist = curve_analysis.Project(
417 curve->GetCurve(), point(space_point), tolerance, proj, t, true);
418 Assert(dist < tolerance * length,
419 ExcPointNotOnManifold<spacedim>(space_point));
420 (void)dist; // Silence compiler warning in Release mode.
421 return Point<1>(GCPnts_AbscissaPoint::Length(
422 curve->GetCurve(), curve->GetCurve().FirstParameter(), t));
423 }
424
425
426
427 template <int dim, int spacedim>
430 const Point<1> &chart_point) const
431 {
432 GCPnts_AbscissaPoint AP(curve->GetCurve(),
433 chart_point[0],
434 curve->GetCurve().FirstParameter());
435 gp_Pnt P = curve->GetCurve().Value(AP.Parameter());
436 return point<spacedim>(P);
437 }
438
439 template <int dim, int spacedim>
441 const double tolerance)
442 : face(face)
443 , tolerance(tolerance)
444 {}
445
446
447
448 template <int dim, int spacedim>
449 std::unique_ptr<Manifold<dim, spacedim>>
451 {
452 return std::unique_ptr<Manifold<dim, spacedim>>(
453 new NURBSPatchManifold<dim, spacedim>(face, tolerance));
454 }
455
456
457
458 template <int dim, int spacedim>
461 const Point<spacedim> &space_point) const
462 {
463 Handle(Geom_Surface) SurfToProj = BRep_Tool::Surface(face);
464
465 ShapeAnalysis_Surface projector(SurfToProj);
466 gp_Pnt2d proj_params = projector.ValueOfUV(point(space_point), tolerance);
467
468 double u = proj_params.X();
469 double v = proj_params.Y();
470
471 return {u, v};
472 }
473
474 template <int dim, int spacedim>
477 const Point<2> &chart_point) const
478 {
479 return ::OpenCASCADE::push_forward<spacedim>(face,
480 chart_point[0],
481 chart_point[1]);
482 }
483
484 template <int dim, int spacedim>
487 const Point<2> &chart_point) const
488 {
490 Handle(Geom_Surface) surf = BRep_Tool::Surface(face);
491
492 gp_Pnt q;
493 gp_Vec Du, Dv;
494 surf->D1(chart_point[0], chart_point[1], q, Du, Dv);
495
496 DX[0][0] = Du.X();
497 DX[1][0] = Du.Y();
498 if (spacedim > 2)
499 DX[2][0] = Du.Z();
500 else
501 Assert(std::abs(Du.Z()) < tolerance,
503 "Expecting derivative along Z to be zero! Bailing out."));
504 DX[0][1] = Dv.X();
505 DX[1][1] = Dv.Y();
506 if (spacedim > 2)
507 DX[2][1] = Dv.Z();
508 else
509 Assert(std::abs(Dv.Z()) < tolerance,
511 "Expecting derivative along Z to be zero! Bailing out."));
512 return DX;
513 }
514
515 template <int dim, int spacedim>
516 std::tuple<double, double, double, double>
518 {
519 Standard_Real umin, umax, vmin, vmax;
520 BRepTools::UVBounds(face, umin, umax, vmin, vmax);
521 return std::make_tuple(umin, umax, vmin, vmax);
522 }
523
524// Explicit instantiations
525# include "manifold_lib.inst"
526} // end namespace OpenCASCADE
527
529
530#endif
numbers::NumberTraits< Number >::real_type norm() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:402
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:403
static ::ExceptionBase & ExcImpossibleInDimSpacedim(int arg1, int arg2)
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
Definition: exceptions.h:1465
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcUnsupportedShape()
static ::ExceptionBase & ExcMessage(std::string arg1)
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
std::tuple< double, double, double, double > get_uv_bounds() const
virtual std::unique_ptr< Manifold< dim, spacedim > > clone() const override
virtual Point< spacedim > push_forward(const Point< 1 > &chart_point) const override
DirectionalProjectionManifold(const TopoDS_Shape &sh, const Tensor< 1, spacedim > &direction, 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
Definition: manifold_lib.cc:93
virtual Point< 2 > pull_back(const Point< spacedim > &space_point) const override
virtual Point< spacedim > push_forward(const Point< 2 > &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 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
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
NormalProjectionManifold(const TopoDS_Shape &sh, const double tolerance=1e-7)
Definition: manifold_lib.cc:80
virtual std::unique_ptr< Manifold< dim, spacedim > > clone() const override
NormalToMeshProjectionManifold(const TopoDS_Shape &sh, const double tolerance=1e-7)
static const char T
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
Definition: divergence.h:472
std::tuple< unsigned int, unsigned int, unsigned int > count_elements(const TopoDS_Shape &shape)
Definition: utilities.cc:88
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:791
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition: utilities.cc:188
Point< dim > closest_point(const TopoDS_Shape &in_shape, const Point< dim > &origin, const double tolerance=1e-7)
Definition: utilities.cc:781
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:506
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)