Reference documentation for deal.II version Git 840703ecf3 2019-10-14 17:15:20 -0400
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
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  AssertThrow(false, ExcNotImplemented());
306  break;
307  }
308  }
309 
310  return line_intersection(sh, candidate, average_normal, tolerance);
311  }
312 
313 
314  /*============================== ArclengthProjectionLineManifold
315  * ==============================*/
316  template <int dim, int spacedim>
319  const double tolerance)
320  :
321 
322  ChartManifold<dim, spacedim, 1>(sh.Closed() ? Point<1>(shape_length(sh)) :
323  Point<1>())
324  , sh(sh)
325  , curve(curve_adaptor(sh))
326  , tolerance(tolerance)
327  , length(shape_length(sh))
328  {
329  Assert(spacedim >= 2, ExcImpossibleInDimSpacedim(dim, spacedim));
330  }
331 
332 
333 
334  template <int dim, int spacedim>
335  std::unique_ptr<Manifold<dim, spacedim>>
337  {
338  return std::unique_ptr<Manifold<dim, spacedim>>(
340  }
341 
342 
343 
344  template <int dim, int spacedim>
345  Point<1>
347  const Point<spacedim> &space_point) const
348  {
349  double t(0.0);
350  ShapeAnalysis_Curve curve_analysis;
351  gp_Pnt proj;
352  const double dist = curve_analysis.Project(
353  curve->GetCurve(), point(space_point), tolerance, proj, t, true);
354  Assert(dist < tolerance * length,
355  ExcPointNotOnManifold<spacedim>(space_point));
356  (void)dist; // Silence compiler warning in Release mode.
357  return Point<1>(GCPnts_AbscissaPoint::Length(
358  curve->GetCurve(), curve->GetCurve().FirstParameter(), t));
359  }
360 
361 
362 
363  template <int dim, int spacedim>
366  const Point<1> &chart_point) const
367  {
368  GCPnts_AbscissaPoint AP(curve->GetCurve(),
369  chart_point[0],
370  curve->GetCurve().FirstParameter());
371  gp_Pnt P = curve->GetCurve().Value(AP.Parameter());
372  return point<spacedim>(P);
373  }
374 
375  template <int dim, int spacedim>
377  const double tolerance)
378  : face(face)
379  , tolerance(tolerance)
380  {}
381 
382 
383 
384  template <int dim, int spacedim>
385  std::unique_ptr<Manifold<dim, spacedim>>
387  {
388  return std::unique_ptr<Manifold<dim, spacedim>>(
390  }
391 
392 
393 
394  template <int dim, int spacedim>
395  Point<2>
397  const Point<spacedim> &space_point) const
398  {
399  Handle(Geom_Surface) SurfToProj = BRep_Tool::Surface(face);
400 
401  ShapeAnalysis_Surface projector(SurfToProj);
402  gp_Pnt2d proj_params = projector.ValueOfUV(point(space_point), tolerance);
403 
404  double u = proj_params.X();
405  double v = proj_params.Y();
406 
407  return {u, v};
408  }
409 
410  template <int dim, int spacedim>
413  const Point<2> &chart_point) const
414  {
415  return ::OpenCASCADE::push_forward<spacedim>(face,
416  chart_point[0],
417  chart_point[1]);
418  }
419 
420  template <int dim, int spacedim>
423  const Point<2> &chart_point) const
424  {
426  Handle(Geom_Surface) surf = BRep_Tool::Surface(face);
427 
428  gp_Pnt q;
429  gp_Vec Du, Dv;
430  surf->D1(chart_point[0], chart_point[1], q, Du, Dv);
431 
432  DX[0][0] = Du.X();
433  DX[1][0] = Du.Y();
434  if (spacedim > 2)
435  DX[2][0] = Du.Z();
436  else
437  Assert(std::abs(Du.Z()) < tolerance,
438  ExcMessage(
439  "Expecting derivative along Z to be zero! Bailing out."));
440  DX[0][1] = Dv.X();
441  DX[1][1] = Dv.Y();
442  if (spacedim > 2)
443  DX[2][1] = Dv.Z();
444  else
445  Assert(std::abs(Dv.Z()) < tolerance,
446  ExcMessage(
447  "Expecting derivative along Z to be zero! Bailing out."));
448  return DX;
449  }
450 
451  template <int dim, int spacedim>
452  std::tuple<double, double, double, double>
454  {
455  Standard_Real umin, umax, vmin, vmax;
456  BRepTools::UVBounds(face, umin, umax, vmin, vmax);
457  return std::make_tuple(umin, umax, vmin, vmax);
458  }
459 
460 // Explicit instantiations
461 # include "manifold_lib.inst"
462 } // end namespace OpenCASCADE
463 
464 DEAL_II_NAMESPACE_CLOSE
465 
466 #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:80
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:422
virtual std::unique_ptr< Manifold< dim, spacedim > > clone() const override
#define AssertThrow(cond, exc)
Definition: exceptions.h:1519
virtual Point< spacedim > project_to_manifold(const ArrayView< const Point< spacedim >> &surrounding_points, const Point< spacedim > &candidate) const override
numbers::NumberTraits< Number >::real_type norm() const
Definition: tensor.h:1416
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:1407
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:698
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition: utilities.cc:180
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()
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:708
static ::ExceptionBase & ExcInternalError()