Reference documentation for deal.II version 9.2.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\}}\)
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 
41 namespace 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>>(
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>>(
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,
174  ExcMessage(
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 <>
201  Point<3>
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,
238  ExcMessage(
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,
268  ExcMessage(
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,
309  ExcMessage(
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,
354  ExcMessage(
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>
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>>(
404  }
405 
406 
407 
408  template <int dim, int spacedim>
409  Point<1>
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>>(
454  }
455 
456 
457 
458  template <int dim, int spacedim>
459  Point<2>
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,
502  ExcMessage(
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,
510  ExcMessage(
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
NormalToMeshProjectionManifold(const TopoDS_Shape &sh, const double tolerance=1e-7)
virtual Point< 2 > pull_back(const Point< spacedim > &space_point) const override
std::tuple< unsigned int, unsigned int, unsigned int > count_elements(const TopoDS_Shape &shape)
Definition: utilities.cc:88
static ::ExceptionBase & ExcImpossibleInDimSpacedim(int arg1, int arg2)
const Tensor< 1, spacedim > direction
Definition: manifold_lib.h:180
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
virtual Point< spacedim > push_forward(const Point< 1 > &chart_point) const override
NormalProjectionManifold(const TopoDS_Shape &sh, const double tolerance=1e-7)
Definition: manifold_lib.cc:80
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
virtual std::unique_ptr< Manifold< dim, spacedim > > clone() const override
NURBSPatchManifold(const TopoDS_Face &face, const double tolerance=1e-7)
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim >>> &Du)
Definition: divergence.h:548
virtual Point< spacedim > project_to_manifold(const ArrayView< const Point< spacedim >> &surrounding_points, const Point< spacedim > &candidate) const override
virtual Point< spacedim > push_forward(const Point< 2 > &chart_point) const override
static ::ExceptionBase & ExcMessage(std::string arg1)
static const char T
virtual Point< spacedim > project_to_manifold(const ArrayView< const Point< spacedim >> &surrounding_points, const Point< spacedim > &candidate) const override
#define Assert(cond, exc)
Definition: exceptions.h:1419
virtual std::unique_ptr< Manifold< dim, spacedim > > clone() const override
virtual std::unique_ptr< Manifold< dim, spacedim > > clone() const override
virtual std::unique_ptr< Manifold< dim, spacedim > > clone() const override
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
Point< dim > closest_point(const TopoDS_Shape &in_shape, const Point< dim > &origin, const double tolerance=1e-7)
Definition: utilities.cc:782
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition: utilities.cc:188
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const
virtual Point< spacedim > project_to_manifold(const ArrayView< const Point< spacedim >> &surrounding_points, const Point< spacedim > &candidate) const override
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:358
virtual std::unique_ptr< Manifold< dim, spacedim > > clone() const override
Definition: manifold_lib.cc:93
static ::ExceptionBase & ExcNotImplemented()
ArclengthProjectionLineManifold(const TopoDS_Shape &sh, const double tolerance=1e-7)
virtual Point< 1 > pull_back(const Point< spacedim > &space_point) const override
static ::ExceptionBase & ExcUnsupportedShape()
numbers::NumberTraits< Number >::real_type norm() const
std::tuple< double, double, double, double > get_uv_bounds() const
virtual DerivativeForm< 1, 2, spacedim > push_forward_gradient(const Point< 2 > &chart_point) const override
DirectionalProjectionManifold(const TopoDS_Shape &sh, const Tensor< 1, spacedim > &direction, const double tolerance=1e-7)
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:792
T max(const T &t, const MPI_Comm &mpi_communicator)
static ::ExceptionBase & ExcInternalError()