Reference documentation for deal.II version Git f40be01994 2020-04-09 07:13:12 +0200
\(\newcommand{\vcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\vcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
manifold_lib.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2014 - 2019 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/opencascade/manifold_lib.h>
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 
38 DEAL_II_NAMESPACE_OPEN
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  * ==============================*/
80  template <int dim, int spacedim>
82  const TopoDS_Shape &sh,
83  const double tolerance)
84  : sh(sh)
85  , tolerance(tolerance)
86  {
87  Assert(spacedim == 3, ExcNotImplemented());
88  }
89 
90 
91 
92  template <int dim, int spacedim>
93  std::unique_ptr<Manifold<dim, spacedim>>
95  {
96  return std::unique_ptr<Manifold<dim, spacedim>>(
98  }
99 
100 
101 
102  template <int dim, int spacedim>
105  const ArrayView<const Point<spacedim>> &surrounding_points,
106  const Point<spacedim> & candidate) const
107  {
108  (void)surrounding_points;
109 # ifdef DEBUG
110  for (unsigned int i = 0; i < surrounding_points.size(); ++i)
111  Assert(closest_point(sh, surrounding_points[i], tolerance)
112  .distance(surrounding_points[i]) <
113  std::max(tolerance * surrounding_points[i].norm(), tolerance),
114  ExcPointNotOnManifold<spacedim>(surrounding_points[i]));
115 # endif
116  return closest_point(sh, candidate, tolerance);
117  }
118 
119 
120  /*============================== DirectionalProjectionManifold
121  * ==============================*/
122  template <int dim, int spacedim>
124  const TopoDS_Shape & sh,
125  const Tensor<1, spacedim> &direction,
126  const double tolerance)
127  : sh(sh)
128  , direction(direction)
129  , tolerance(tolerance)
130  {
131  Assert(spacedim == 3, ExcNotImplemented());
132  }
133 
134 
135 
136  template <int dim, int spacedim>
137  std::unique_ptr<Manifold<dim, spacedim>>
139  {
140  return std::unique_ptr<Manifold<dim, spacedim>>(
142  }
143 
144 
145 
146  template <int dim, int spacedim>
149  const ArrayView<const Point<spacedim>> &surrounding_points,
150  const Point<spacedim> & candidate) const
151  {
152  (void)surrounding_points;
153 # ifdef DEBUG
154  for (unsigned int i = 0; i < surrounding_points.size(); ++i)
155  Assert(closest_point(sh, surrounding_points[i], tolerance)
156  .distance(surrounding_points[i]) <
157  std::max(tolerance * surrounding_points[i].norm(), tolerance),
158  ExcPointNotOnManifold<spacedim>(surrounding_points[i]));
159 # endif
160  return line_intersection(sh, candidate, direction, tolerance);
161  }
162 
163 
164 
165  /*============================== NormalToMeshProjectionManifold
166  * ==============================*/
167  template <int dim, int spacedim>
169  const TopoDS_Shape &sh,
170  const double tolerance)
171  : sh(sh)
172  , tolerance(tolerance)
173  {
174  Assert(spacedim == 3, ExcNotImplemented());
175  Assert(
176  std::get<0>(count_elements(sh)) > 0,
177  ExcMessage(
178  "NormalToMeshProjectionManifold needs a shape containing faces to operate."));
179  }
180 
181  template <int dim, int spacedim>
182  std::unique_ptr<Manifold<dim, spacedim>>
184  {
185  return std::unique_ptr<Manifold<dim, spacedim>>(
187  }
188 
189 
190  template <int dim, int spacedim>
193  const ArrayView<const Point<spacedim>> &surrounding_points,
194  const Point<spacedim> & candidate) const
195  {
196  TopoDS_Shape out_shape;
197  Tensor<1, 3> average_normal;
198 # ifdef DEBUG
199  for (unsigned int i = 0; i < surrounding_points.size(); ++i)
200  {
201  Assert(closest_point(sh, surrounding_points[i], tolerance)
202  .distance(surrounding_points[i]) <
203  std::max(tolerance * surrounding_points[i].norm(), tolerance),
204  ExcPointNotOnManifold<spacedim>(surrounding_points[i]));
205  }
206 # endif
207 
208  switch (surrounding_points.size())
209  {
210  case 2:
211  {
212  for (unsigned int i = 0; i < surrounding_points.size(); ++i)
213  {
214  std::tuple<Point<3>, Tensor<1, 3>, double, double>
215  p_and_diff_forms =
217  surrounding_points[i],
218  tolerance);
219  average_normal += std::get<1>(p_and_diff_forms);
220  }
221 
222  average_normal /= 2.0;
223 
224  Assert(
225  average_normal.norm() > 1e-4,
226  ExcMessage(
227  "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."));
228 
229  Tensor<1, 3> T = surrounding_points[0] - surrounding_points[1];
230  T /= T.norm();
231  average_normal = average_normal - (average_normal * T) * T;
232  average_normal /= average_normal.norm();
233  break;
234  }
235  case 4:
236  {
237  Tensor<1, 3> u = surrounding_points[1] - surrounding_points[0];
238  Tensor<1, 3> v = surrounding_points[2] - surrounding_points[0];
239  const double n1_coords[3] = {u[1] * v[2] - u[2] * v[1],
240  u[2] * v[0] - u[0] * v[2],
241  u[0] * v[1] - u[1] * v[0]};
242  Tensor<1, 3> n1(n1_coords);
243  n1 = n1 / n1.norm();
244  u = surrounding_points[2] - surrounding_points[3];
245  v = surrounding_points[1] - surrounding_points[3];
246  const double n2_coords[3] = {u[1] * v[2] - u[2] * v[1],
247  u[2] * v[0] - u[0] * v[2],
248  u[0] * v[1] - u[1] * v[0]};
249  Tensor<1, 3> n2(n2_coords);
250  n2 = n2 / n2.norm();
251 
252  average_normal = (n1 + n2) / 2.0;
253 
254  Assert(
255  average_normal.norm() > tolerance,
256  ExcMessage(
257  "Failed to refine cell: the normal estimated via the surrounding points turns out to be a null vector, making the projection direction undetermined."));
258 
259  average_normal /= average_normal.norm();
260  break;
261  }
262  case 8:
263  {
264  Tensor<1, 3> u = surrounding_points[1] - surrounding_points[0];
265  Tensor<1, 3> v = surrounding_points[2] - surrounding_points[0];
266  const double n1_coords[3] = {u[1] * v[2] - u[2] * v[1],
267  u[2] * v[0] - u[0] * v[2],
268  u[0] * v[1] - u[1] * v[0]};
269  Tensor<1, 3> n1(n1_coords);
270  n1 = n1 / n1.norm();
271  u = surrounding_points[2] - surrounding_points[3];
272  v = surrounding_points[1] - surrounding_points[3];
273  const double n2_coords[3] = {u[1] * v[2] - u[2] * v[1],
274  u[2] * v[0] - u[0] * v[2],
275  u[0] * v[1] - u[1] * v[0]};
276  Tensor<1, 3> n2(n2_coords);
277  n2 = n2 / n2.norm();
278  u = surrounding_points[4] - surrounding_points[7];
279  v = surrounding_points[6] - surrounding_points[7];
280  const double n3_coords[3] = {u[1] * v[2] - u[2] * v[1],
281  u[2] * v[0] - u[0] * v[2],
282  u[0] * v[1] - u[1] * v[0]};
283  Tensor<1, 3> n3(n3_coords);
284  n3 = n3 / n3.norm();
285  u = surrounding_points[6] - surrounding_points[7];
286  v = surrounding_points[5] - surrounding_points[7];
287  const double n4_coords[3] = {u[1] * v[2] - u[2] * v[1],
288  u[2] * v[0] - u[0] * v[2],
289  u[0] * v[1] - u[1] * v[0]};
290  Tensor<1, 3> n4(n4_coords);
291  n4 = n4 / n4.norm();
292 
293  average_normal = (n1 + n2 + n3 + n4) / 4.0;
294 
295  Assert(
296  average_normal.norm() > tolerance,
297  ExcMessage(
298  "Failed to refine cell: the normal estimated via the surrounding points turns out to be a null vector, making the projection direction undetermined."));
299 
300  average_normal /= average_normal.norm();
301  break;
302  }
303  default:
304  {
305  // Given an arbitrary number of points we compute all the possible
306  // normal vectors
307  for (unsigned int i = 0; i < surrounding_points.size(); ++i)
308  for (unsigned int j = 0; j < surrounding_points.size(); ++j)
309  if (j != i)
310  for (unsigned int k = 0; k < surrounding_points.size(); ++k)
311  if (k != j && k != i)
312  {
313  Tensor<1, 3> u =
314  surrounding_points[i] - surrounding_points[j];
315  Tensor<1, 3> v =
316  surrounding_points[i] - surrounding_points[k];
317  const double n_coords[3] = {u[1] * v[2] - u[2] * v[1],
318  u[2] * v[0] - u[0] * v[2],
319  u[0] * v[1] - u[1] * v[0]};
320  Tensor<1, 3> n1(n_coords);
321  if (n1.norm() > tolerance)
322  {
323  n1 = n1 / n1.norm();
324  if (average_normal.norm() < tolerance)
325  average_normal = n1;
326  else
327  {
328  auto dot_prod = n1 * average_normal;
329  // We check that the direction of the normal
330  // vector w.r.t the current average, and make
331  // sure we flip it if it is opposite
332  if (dot_prod > 0)
333  average_normal += n1;
334  else
335  average_normal -= n1;
336  }
337  }
338  }
339  Assert(
340  average_normal.norm() > tolerance,
341  ExcMessage(
342  "Failed to compute a normal: the normal estimated via the surrounding points turns out to be a null vector, making the projection direction undetermined."));
343  average_normal = average_normal / average_normal.norm();
344  break;
345  }
346  }
347 
348  return line_intersection(sh, candidate, average_normal, tolerance);
349  }
350 
351 
352  /*============================== ArclengthProjectionLineManifold
353  * ==============================*/
354  template <int dim, int spacedim>
357  const double tolerance)
358  :
359 
360  ChartManifold<dim, spacedim, 1>(sh.Closed() ? Point<1>(shape_length(sh)) :
361  Point<1>())
362  , sh(sh)
363  , curve(curve_adaptor(sh))
364  , tolerance(tolerance)
365  , length(shape_length(sh))
366  {
367  Assert(spacedim >= 2, ExcImpossibleInDimSpacedim(dim, spacedim));
368  }
369 
370 
371 
372  template <int dim, int spacedim>
373  std::unique_ptr<Manifold<dim, spacedim>>
375  {
376  return std::unique_ptr<Manifold<dim, spacedim>>(
378  }
379 
380 
381 
382  template <int dim, int spacedim>
383  Point<1>
385  const Point<spacedim> &space_point) const
386  {
387  double t(0.0);
388  ShapeAnalysis_Curve curve_analysis;
389  gp_Pnt proj;
390  const double dist = curve_analysis.Project(
391  curve->GetCurve(), point(space_point), tolerance, proj, t, true);
392  Assert(dist < tolerance * length,
393  ExcPointNotOnManifold<spacedim>(space_point));
394  (void)dist; // Silence compiler warning in Release mode.
395  return Point<1>(GCPnts_AbscissaPoint::Length(
396  curve->GetCurve(), curve->GetCurve().FirstParameter(), t));
397  }
398 
399 
400 
401  template <int dim, int spacedim>
404  const Point<1> &chart_point) const
405  {
406  GCPnts_AbscissaPoint AP(curve->GetCurve(),
407  chart_point[0],
408  curve->GetCurve().FirstParameter());
409  gp_Pnt P = curve->GetCurve().Value(AP.Parameter());
410  return point<spacedim>(P);
411  }
412 
413  template <int dim, int spacedim>
415  const double tolerance)
416  : face(face)
417  , tolerance(tolerance)
418  {}
419 
420 
421 
422  template <int dim, int spacedim>
423  std::unique_ptr<Manifold<dim, spacedim>>
425  {
426  return std::unique_ptr<Manifold<dim, spacedim>>(
428  }
429 
430 
431 
432  template <int dim, int spacedim>
433  Point<2>
435  const Point<spacedim> &space_point) const
436  {
437  Handle(Geom_Surface) SurfToProj = BRep_Tool::Surface(face);
438 
439  ShapeAnalysis_Surface projector(SurfToProj);
440  gp_Pnt2d proj_params = projector.ValueOfUV(point(space_point), tolerance);
441 
442  double u = proj_params.X();
443  double v = proj_params.Y();
444 
445  return {u, v};
446  }
447 
448  template <int dim, int spacedim>
451  const Point<2> &chart_point) const
452  {
453  return ::OpenCASCADE::push_forward<spacedim>(face,
454  chart_point[0],
455  chart_point[1]);
456  }
457 
458  template <int dim, int spacedim>
461  const Point<2> &chart_point) const
462  {
464  Handle(Geom_Surface) surf = BRep_Tool::Surface(face);
465 
466  gp_Pnt q;
467  gp_Vec Du, Dv;
468  surf->D1(chart_point[0], chart_point[1], q, Du, Dv);
469 
470  DX[0][0] = Du.X();
471  DX[1][0] = Du.Y();
472  if (spacedim > 2)
473  DX[2][0] = Du.Z();
474  else
475  Assert(std::abs(Du.Z()) < tolerance,
476  ExcMessage(
477  "Expecting derivative along Z to be zero! Bailing out."));
478  DX[0][1] = Dv.X();
479  DX[1][1] = Dv.Y();
480  if (spacedim > 2)
481  DX[2][1] = Dv.Z();
482  else
483  Assert(std::abs(Dv.Z()) < tolerance,
484  ExcMessage(
485  "Expecting derivative along Z to be zero! Bailing out."));
486  return DX;
487  }
488 
489  template <int dim, int spacedim>
490  std::tuple<double, double, double, double>
492  {
493  Standard_Real umin, umax, vmin, vmax;
494  BRepTools::UVBounds(face, umin, umax, vmin, vmax);
495  return std::make_tuple(umin, umax, vmin, vmax);
496  }
497 
498 // Explicit instantiations
499 # include "manifold_lib.inst"
500 } // end namespace OpenCASCADE
501 
502 DEAL_II_NAMESPACE_CLOSE
503 
504 #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
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:81
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
virtual std::unique_ptr< Manifold< dim, spacedim > > clone() const override
NURBSPatchManifold(const TopoDS_Face &face, 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 Point< spacedim > push_forward(const Point< 2 > &chart_point) const override
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
#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
Point< dim > closest_point(const TopoDS_Shape &in_shape, const Point< dim > &origin, const double tolerance=1e-7)
Definition: utilities.cc:787
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
virtual std::unique_ptr< Manifold< dim, spacedim > > clone() const override
Definition: manifold_lib.cc:94
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:797
static ::ExceptionBase & ExcInternalError()