Reference documentation for deal.II version GIT 112f7bbc69 2023-05-28 21:25:02+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\}}\)
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 
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 
44 namespace 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,
201  ExcMessage(
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 <>
228  Point<3>
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 (unsigned int i = 0; i < surrounding_points.size(); ++i)
240  {
241  Assert(closest_point(sh, surrounding_points[i], tolerance)
242  .distance(surrounding_points[i]) <
243  std::max(tolerance * surrounding_points[i].norm(),
244  tolerance),
245  ExcPointNotOnManifold<spacedim>(surrounding_points[i]));
246  }
247 # endif
248 
249  switch (surrounding_points.size())
250  {
251  case 2:
252  {
253  for (unsigned int i = 0; i < surrounding_points.size(); ++i)
254  {
255  std::tuple<Point<3>, Tensor<1, 3>, double, double>
256  p_and_diff_forms = closest_point_and_differential_forms(
257  sh, surrounding_points[i], 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,
265  ExcMessage(
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,
295  ExcMessage(
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,
336  ExcMessage(
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,
381  ExcMessage(
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>
436  Point<1>
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>
503  Point<2>
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,
546  ExcMessage(
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,
554  ExcMessage(
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)
numbers::NumberTraits< Number >::real_type norm() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:474
#define DEAL_II_OPENCASCADE_VERSION_GTE(major, minor, subminor)
Definition: config.h:325
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:475
static ::ExceptionBase & ExcImpossibleInDimSpacedim(int arg1, int arg2)
static ::ExceptionBase & ExcInternalError()
#define Assert(cond, exc)
Definition: exceptions.h:1614
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcUnsupportedShape()
static ::ExceptionBase & ExcMessage(std::string arg1)
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: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:796
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:786
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
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)