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;
135 for (
unsigned int i = 0; i < surrounding_points.size(); ++i)
137 .distance(surrounding_points[i]) <
138 std::max(tolerance * surrounding_points[i].norm(), tolerance),
146 template <
int dim,
int spacedim>
148 const TopoDS_Shape &sh,
150 const double tolerance)
152 , direction(direction)
153 , tolerance(tolerance)
160 template <
int dim,
int spacedim>
161 std::unique_ptr<Manifold<dim, spacedim>>
164 return std::unique_ptr<Manifold<dim, spacedim>>(
170 template <
int dim,
int spacedim>
176 (void)surrounding_points;
178 for (
unsigned int i = 0; i < surrounding_points.size(); ++i)
180 .distance(surrounding_points[i]) <
181 std::max(tolerance * surrounding_points[i].norm(), tolerance),
190 template <
int dim,
int spacedim>
192 const TopoDS_Shape &sh,
193 const double tolerance)
195 , tolerance(tolerance)
201 "NormalToMeshProjectionManifold needs a shape containing faces to operate."));
204 template <
int dim,
int spacedim>
205 std::unique_ptr<Manifold<dim, spacedim>>
208 return std::unique_ptr<Manifold<dim, spacedim>>(
215 template <
int spacedim>
217 internal_project_to_manifold(
const TopoDS_Shape &,
228 internal_project_to_manifold(
229 const TopoDS_Shape &sh,
230 const double tolerance,
234 constexpr int spacedim = 3;
235 TopoDS_Shape out_shape;
238 for (
const auto &
point : surrounding_points)
246 switch (surrounding_points.size())
250 for (
const auto &
point : surrounding_points)
257 average_normal += std::get<1>(p_and_diff_forms);
260 average_normal /= 2.0;
263 average_normal.
norm() > 1e-4,
265 "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 Tensor<1, 3> T = surrounding_points[0] - surrounding_points[1];
269 average_normal = average_normal - (average_normal * T) * T;
270 average_normal /= average_normal.
norm();
275 Tensor<1, 3> u = surrounding_points[1] - surrounding_points[0];
276 Tensor<1, 3> v = surrounding_points[2] - surrounding_points[0];
277 const double n1_coords[3] = {u[1] * v[2] - u[2] * v[1],
278 u[2] * v[0] - u[0] * v[2],
279 u[0] * v[1] - u[1] * v[0]};
282 u = surrounding_points[2] - surrounding_points[3];
283 v = surrounding_points[1] - surrounding_points[3];
284 const double n2_coords[3] = {u[1] * v[2] - u[2] * v[1],
285 u[2] * v[0] - u[0] * v[2],
286 u[0] * v[1] - u[1] * v[0]};
290 average_normal = (n1 + n2) / 2.0;
293 average_normal.
norm() > tolerance,
295 "Failed to refine cell: the normal estimated via the surrounding points turns out to be a null vector, making the projection direction undetermined."));
297 average_normal /= average_normal.
norm();
302 Tensor<1, 3> u = surrounding_points[1] - surrounding_points[0];
303 Tensor<1, 3> v = surrounding_points[2] - surrounding_points[0];
304 const double n1_coords[3] = {u[1] * v[2] - u[2] * v[1],
305 u[2] * v[0] - u[0] * v[2],
306 u[0] * v[1] - u[1] * v[0]};
309 u = surrounding_points[2] - surrounding_points[3];
310 v = surrounding_points[1] - surrounding_points[3];
311 const double n2_coords[3] = {u[1] * v[2] - u[2] * v[1],
312 u[2] * v[0] - u[0] * v[2],
313 u[0] * v[1] - u[1] * v[0]};
316 u = surrounding_points[4] - surrounding_points[7];
317 v = surrounding_points[6] - surrounding_points[7];
318 const double n3_coords[3] = {u[1] * v[2] - u[2] * v[1],
319 u[2] * v[0] - u[0] * v[2],
320 u[0] * v[1] - u[1] * v[0]};
323 u = surrounding_points[6] - surrounding_points[7];
324 v = surrounding_points[5] - surrounding_points[7];
325 const double n4_coords[3] = {u[1] * v[2] - u[2] * v[1],
326 u[2] * v[0] - u[0] * v[2],
327 u[0] * v[1] - u[1] * v[0]};
331 average_normal = (n1 + n2 + n3 + n4) / 4.0;
334 average_normal.
norm() > tolerance,
336 "Failed to refine cell: the normal estimated via the surrounding points turns out to be a null vector, making the projection direction undetermined."));
338 average_normal /= average_normal.
norm();
345 for (
unsigned int i = 0; i < surrounding_points.size(); ++i)
346 for (
unsigned int j = 0; j < surrounding_points.size(); ++j)
348 for (
unsigned int k = 0; k < surrounding_points.size(); ++k)
349 if (k != j && k != i)
352 surrounding_points[i] - surrounding_points[j];
354 surrounding_points[i] - surrounding_points[k];
355 const double n_coords[3] = {u[1] * v[2] - u[2] * v[1],
356 u[2] * v[0] - u[0] * v[2],
360 if (n1.norm() > tolerance)
363 if (average_normal.
norm() < tolerance)
367 auto dot_prod = n1 * average_normal;
372 average_normal += n1;
374 average_normal -= n1;
379 average_normal.
norm() > tolerance,
381 "Failed to compute a normal: the normal estimated via the surrounding points turns out to be a null vector, making the projection direction undetermined."));
382 average_normal = average_normal / average_normal.
norm();
392 template <
int dim,
int spacedim>
398 return internal_project_to_manifold(sh,
406 template <
int dim,
int spacedim>
409 const double tolerance)
415 , curve(curve_adaptor(sh))
416 , tolerance(tolerance)
417 , length(shape_length(sh))
424 template <
int dim,
int spacedim>
425 std::unique_ptr<Manifold<dim, spacedim>>
428 return std::unique_ptr<Manifold<dim, spacedim>>(
434 template <
int dim,
int spacedim>
440 ShapeAnalysis_Curve curve_analysis;
443 const double dist = curve_analysis.Project(
445 *curve,
point(space_point), tolerance, proj, t,
true);
447 curve->GetCurve(),
point(space_point), tolerance, proj, t,
true);
451 Assert(dist < tolerance * length,
454 return Point<1>(GCPnts_AbscissaPoint::Length(
456 *curve, curve->FirstParameter(), t));
458 curve->GetCurve(), curve->GetCurve().FirstParameter(), t));
464 template <
int dim,
int spacedim>
469# if DEAL_II_OPENCASCADE_VERSION_GTE(7, 6, 0)
470 GCPnts_AbscissaPoint AP(*curve, chart_point[0], curve->FirstParameter());
471 gp_Pnt P = curve->Value(AP.Parameter());
473 GCPnts_AbscissaPoint AP(curve->GetCurve(),
475 curve->GetCurve().FirstParameter());
476 gp_Pnt P = curve->GetCurve().Value(AP.Parameter());
482 template <
int dim,
int spacedim>
484 const double tolerance)
486 , tolerance(tolerance)
491 template <
int dim,
int spacedim>
492 std::unique_ptr<Manifold<dim, spacedim>>
495 return std::unique_ptr<Manifold<dim, spacedim>>(
501 template <
int dim,
int spacedim>
506 Handle(Geom_Surface) SurfToProj = BRep_Tool::Surface(face);
508 ShapeAnalysis_Surface projector(SurfToProj);
509 gp_Pnt2d proj_params = projector.ValueOfUV(
point(space_point), tolerance);
511 double u = proj_params.X();
512 double v = proj_params.Y();
517 template <
int dim,
int spacedim>
522 return ::OpenCASCADE::push_forward<spacedim>(face,
527 template <
int dim,
int spacedim>
533 Handle(Geom_Surface) surf = BRep_Tool::Surface(face);
537 surf->D1(chart_point[0], chart_point[1], q, Du, Dv);
546 "Expecting derivative along Z to be zero! Bailing out."));
554 "Expecting derivative along Z to be zero! Bailing out."));
558 template <
int dim,
int spacedim>
559 std::tuple<double, double, double, double>
562 Standard_Real umin, umax, vmin, vmax;
563 BRepTools::UVBounds(face, umin, umax, vmin, vmax);
564 return std::make_tuple(umin, umax, vmin, vmax);
568# include "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
#define DEAL_II_OPENCASCADE_VERSION_GTE(major, minor, subminor)
#define DEAL_II_NAMESPACE_CLOSE
static ::ExceptionBase & ExcImpossibleInDimSpacedim(int arg1, int arg2)
static ::ExceptionBase & ExcPointNotOnManifold(Point< dim > arg1)
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcUnsupportedShape()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define DEAL_II_ASSERT_UNREACHABLE()
#define DEAL_II_NOT_IMPLEMENTED()
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 > &)