20#ifdef DEAL_II_WITH_OPENCASCADE
23# include <BRepAdaptor_CompCurve.hxx>
24# include <BRepAdaptor_Curve.hxx>
25# if !DEAL_II_OPENCASCADE_VERSION_GTE(7, 6, 0)
26# include <BRepAdaptor_HCompCurve.hxx>
27# include <BRepAdaptor_HCurve.hxx>
29# include <BRepTools.hxx>
30# include <BRep_Tool.hxx>
31# include <GCPnts_AbscissaPoint.hxx>
32# include <ShapeAnalysis_Curve.hxx>
33# include <ShapeAnalysis_Surface.hxx>
35# if !DEAL_II_OPENCASCADE_VERSION_GTE(7, 0, 0)
36# include <Handle_Adaptor3d_HCurve.hxx>
52# if DEAL_II_OPENCASCADE_VERSION_GTE(7, 6, 0)
53 Handle_Adaptor3d_Curve
54 curve_adaptor(
const TopoDS_Shape &shape)
56 Assert((shape.ShapeType() == TopAbs_WIRE) ||
57 (shape.ShapeType() == TopAbs_EDGE),
59 if (shape.ShapeType() == TopAbs_WIRE)
60 return Handle(BRepAdaptor_CompCurve)(
61 new BRepAdaptor_CompCurve(TopoDS::Wire(shape)));
62 else if (shape.ShapeType() == TopAbs_EDGE)
63 return Handle(BRepAdaptor_Curve)(
64 new BRepAdaptor_Curve(TopoDS::Edge(shape)));
67 return Handle(BRepAdaptor_Curve)(
new BRepAdaptor_Curve());
70 Handle_Adaptor3d_HCurve
71 curve_adaptor(
const TopoDS_Shape &shape)
73 Assert((shape.ShapeType() == TopAbs_WIRE) ||
74 (shape.ShapeType() == TopAbs_EDGE),
76 if (shape.ShapeType() == TopAbs_WIRE)
77 return Handle(BRepAdaptor_HCompCurve)(
78 new BRepAdaptor_HCompCurve(TopoDS::Wire(shape)));
79 else if (shape.ShapeType() == TopAbs_EDGE)
80 return Handle(BRepAdaptor_HCurve)(
81 new BRepAdaptor_HCurve(TopoDS::Edge(shape)));
84 return Handle(BRepAdaptor_HCurve)(
new BRepAdaptor_HCurve());
92 shape_length(
const TopoDS_Shape &sh)
94# if DEAL_II_OPENCASCADE_VERSION_GTE(7, 6, 0)
95 Handle_Adaptor3d_Curve adapt = curve_adaptor(sh);
96 return GCPnts_AbscissaPoint::Length(*adapt);
98 Handle_Adaptor3d_HCurve adapt = curve_adaptor(sh);
99 return GCPnts_AbscissaPoint::Length(adapt->GetCurve());
105 template <
int dim,
int spacedim>
107 const TopoDS_Shape &sh,
108 const double tolerance)
110 , tolerance(tolerance)
117 template <
int dim,
int spacedim>
118 std::unique_ptr<Manifold<dim, spacedim>>
121 return std::unique_ptr<Manifold<dim, spacedim>>(
127 template <
int dim,
int spacedim>
133 (void)surrounding_points;
136 for (
unsigned int i = 0; i < surrounding_points.size(); ++i)
138 .distance(surrounding_points[i]) <
139 std::max(tolerance * surrounding_points[i].norm(),
141 ExcPointNotOnManifold<spacedim>(surrounding_points[i]));
148 template <
int dim,
int spacedim>
150 const TopoDS_Shape &sh,
152 const double tolerance)
154 , direction(direction)
155 , tolerance(tolerance)
162 template <
int dim,
int spacedim>
163 std::unique_ptr<Manifold<dim, spacedim>>
166 return std::unique_ptr<Manifold<dim, spacedim>>(
172 template <
int dim,
int spacedim>
178 (void)surrounding_points;
181 for (
unsigned int i = 0; i < surrounding_points.size(); ++i)
183 .distance(surrounding_points[i]) <
184 std::max(tolerance * surrounding_points[i].norm(),
186 ExcPointNotOnManifold<spacedim>(surrounding_points[i]));
194 template <
int dim,
int spacedim>
196 const TopoDS_Shape &sh,
197 const double tolerance)
199 , tolerance(tolerance)
205 "NormalToMeshProjectionManifold needs a shape containing faces to operate."));
208 template <
int dim,
int spacedim>
209 std::unique_ptr<Manifold<dim, spacedim>>
212 return std::unique_ptr<Manifold<dim, spacedim>>(
219 template <
int spacedim>
221 internal_project_to_manifold(
const TopoDS_Shape &,
232 internal_project_to_manifold(
233 const TopoDS_Shape &sh,
234 const double tolerance,
238 constexpr int spacedim = 3;
239 TopoDS_Shape out_shape;
243 for (
const auto &
point : surrounding_points)
247 ExcPointNotOnManifold<spacedim>(
point));
251 switch (surrounding_points.size())
255 for (
const auto &
point : surrounding_points)
262 average_normal += std::get<1>(p_and_diff_forms);
265 average_normal /= 2.0;
268 average_normal.
norm() > 1e-4,
270 "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."));
272 Tensor<1, 3> T = surrounding_points[0] - surrounding_points[1];
274 average_normal = average_normal - (average_normal *
T) * T;
275 average_normal /= average_normal.
norm();
280 Tensor<1, 3> u = surrounding_points[1] - surrounding_points[0];
281 Tensor<1, 3> v = surrounding_points[2] - surrounding_points[0];
282 const double n1_coords[3] = {u[1] * v[2] - u[2] * v[1],
283 u[2] * v[0] - u[0] * v[2],
284 u[0] * v[1] - u[1] * v[0]};
287 u = surrounding_points[2] - surrounding_points[3];
288 v = surrounding_points[1] - surrounding_points[3];
289 const double n2_coords[3] = {u[1] * v[2] - u[2] * v[1],
290 u[2] * v[0] - u[0] * v[2],
291 u[0] * v[1] - u[1] * v[0]};
295 average_normal = (n1 + n2) / 2.0;
298 average_normal.
norm() > tolerance,
300 "Failed to refine cell: the normal estimated via the surrounding points turns out to be a null vector, making the projection direction undetermined."));
302 average_normal /= average_normal.
norm();
307 Tensor<1, 3> u = surrounding_points[1] - surrounding_points[0];
308 Tensor<1, 3> v = surrounding_points[2] - surrounding_points[0];
309 const double n1_coords[3] = {u[1] * v[2] - u[2] * v[1],
310 u[2] * v[0] - u[0] * v[2],
311 u[0] * v[1] - u[1] * v[0]};
314 u = surrounding_points[2] - surrounding_points[3];
315 v = surrounding_points[1] - surrounding_points[3];
316 const double n2_coords[3] = {u[1] * v[2] - u[2] * v[1],
317 u[2] * v[0] - u[0] * v[2],
318 u[0] * v[1] - u[1] * v[0]};
321 u = surrounding_points[4] - surrounding_points[7];
322 v = surrounding_points[6] - surrounding_points[7];
323 const double n3_coords[3] = {u[1] * v[2] - u[2] * v[1],
324 u[2] * v[0] - u[0] * v[2],
325 u[0] * v[1] - u[1] * v[0]};
328 u = surrounding_points[6] - surrounding_points[7];
329 v = surrounding_points[5] - surrounding_points[7];
330 const double n4_coords[3] = {u[1] * v[2] - u[2] * v[1],
331 u[2] * v[0] - u[0] * v[2],
332 u[0] * v[1] - u[1] * v[0]};
336 average_normal = (n1 + n2 + n3 + n4) / 4.0;
339 average_normal.
norm() > tolerance,
341 "Failed to refine cell: the normal estimated via the surrounding points turns out to be a null vector, making the projection direction undetermined."));
343 average_normal /= average_normal.
norm();
350 for (
unsigned int i = 0; i < surrounding_points.size(); ++i)
351 for (
unsigned int j = 0; j < surrounding_points.size(); ++j)
353 for (
unsigned int k = 0; k < surrounding_points.size(); ++k)
354 if (k != j && k != i)
357 surrounding_points[i] - surrounding_points[j];
359 surrounding_points[i] - surrounding_points[k];
360 const double n_coords[3] = {u[1] * v[2] - u[2] * v[1],
361 u[2] * v[0] - u[0] * v[2],
365 if (n1.norm() > tolerance)
368 if (average_normal.
norm() < tolerance)
372 auto dot_prod = n1 * average_normal;
377 average_normal += n1;
379 average_normal -= n1;
384 average_normal.
norm() > tolerance,
386 "Failed to compute a normal: the normal estimated via the surrounding points turns out to be a null vector, making the projection direction undetermined."));
387 average_normal = average_normal / average_normal.
norm();
397 template <
int dim,
int spacedim>
403 return internal_project_to_manifold(sh,
411 template <
int dim,
int spacedim>
414 const double tolerance)
420 , curve(curve_adaptor(sh))
421 , tolerance(tolerance)
422 , length(shape_length(sh))
429 template <
int dim,
int spacedim>
430 std::unique_ptr<Manifold<dim, spacedim>>
433 return std::unique_ptr<Manifold<dim, spacedim>>(
439 template <
int dim,
int spacedim>
445 ShapeAnalysis_Curve curve_analysis;
448 const double dist = curve_analysis.Project(
450 *curve,
point(space_point), tolerance, proj, t,
true);
452 curve->GetCurve(),
point(space_point), tolerance, proj, t,
true);
456 Assert(dist < tolerance * length,
457 ExcPointNotOnManifold<spacedim>(space_point));
459 return Point<1>(GCPnts_AbscissaPoint::Length(
461 *curve, curve->FirstParameter(), t));
463 curve->GetCurve(), curve->GetCurve().FirstParameter(), t));
469 template <
int dim,
int spacedim>
474# if DEAL_II_OPENCASCADE_VERSION_GTE(7, 6, 0)
475 GCPnts_AbscissaPoint AP(*curve, chart_point[0], curve->FirstParameter());
476 gp_Pnt P = curve->Value(AP.Parameter());
478 GCPnts_AbscissaPoint AP(curve->GetCurve(),
480 curve->GetCurve().FirstParameter());
481 gp_Pnt P = curve->GetCurve().Value(AP.Parameter());
484 return point<spacedim>(P);
487 template <
int dim,
int spacedim>
489 const double tolerance)
491 , tolerance(tolerance)
496 template <
int dim,
int spacedim>
497 std::unique_ptr<Manifold<dim, spacedim>>
500 return std::unique_ptr<Manifold<dim, spacedim>>(
506 template <
int dim,
int spacedim>
511 Handle(Geom_Surface) SurfToProj = BRep_Tool::Surface(face);
513 ShapeAnalysis_Surface projector(SurfToProj);
514 gp_Pnt2d proj_params = projector.ValueOfUV(
point(space_point), tolerance);
516 double u = proj_params.X();
517 double v = proj_params.Y();
522 template <
int dim,
int spacedim>
527 return ::OpenCASCADE::push_forward<spacedim>(face,
532 template <
int dim,
int spacedim>
538 Handle(Geom_Surface) surf = BRep_Tool::Surface(face);
542 surf->D1(chart_point[0], chart_point[1], q, Du, Dv);
551 "Expecting derivative along Z to be zero! Bailing out."));
559 "Expecting derivative along Z to be zero! Bailing out."));
563 template <
int dim,
int spacedim>
564 std::tuple<double, double, double, double>
567 Standard_Real umin, umax, vmin, vmax;
568 BRepTools::UVBounds(face, umin, umax, vmin, vmax);
569 return std::make_tuple(umin, umax, vmin, vmax);
577# include "opencascade/manifold_lib.inst"
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
constexpr bool running_in_debug_mode()
#define DEAL_II_OPENCASCADE_VERSION_GTE(major, minor, subminor)
#define DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_ASSERT_UNREACHABLE()
#define DEAL_II_NOT_IMPLEMENTED()
static ::ExceptionBase & ExcImpossibleInDimSpacedim(int arg1, int arg2)
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcUnsupportedShape()
static ::ExceptionBase & ExcMessage(std::string arg1)
std::tuple< unsigned int, unsigned int, unsigned int > count_elements(const TopoDS_Shape &shape)
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)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Point< dim > closest_point(const TopoDS_Shape &in_shape, const Point< dim > &origin, const double tolerance=1e-7)
Point< dim > line_intersection(const TopoDS_Shape &in_shape, const Point< dim > &origin, const Tensor< 1, dim > &direction, const double tolerance=1e-7)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)