18#ifdef DEAL_II_WITH_OPENCASCADE
20# include <BRepAdaptor_CompCurve.hxx>
21# include <BRepAdaptor_Curve.hxx>
22# if !DEAL_II_OPENCASCADE_VERSION_GTE(7, 6, 0)
23# include <BRepAdaptor_HCompCurve.hxx>
24# 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>
38#ifdef DEAL_II_WITH_OPENCASCADE
50# if DEAL_II_OPENCASCADE_VERSION_GTE(7, 6, 0)
51 opencascade::handle<Adaptor3d_Curve>
52 curve_adaptor(
const TopoDS_Shape &shape)
54 Assert((shape.ShapeType() == TopAbs_WIRE) ||
55 (shape.ShapeType() == TopAbs_EDGE),
57 if (shape.ShapeType() == TopAbs_WIRE)
58 return opencascade::handle<BRepAdaptor_CompCurve>(
59 new BRepAdaptor_CompCurve(TopoDS::Wire(shape)));
60 else if (shape.ShapeType() == TopAbs_EDGE)
61 return opencascade::handle<BRepAdaptor_Curve>(
62 new BRepAdaptor_Curve(TopoDS::Edge(shape)));
65 return opencascade::handle<BRepAdaptor_Curve>(
new BRepAdaptor_Curve());
69 opencascade::handle<Adaptor3d_HCurve>
70 curve_adaptor(
const TopoDS_Shape &shape)
72 Assert((shape.ShapeType() == TopAbs_WIRE) ||
73 (shape.ShapeType() == TopAbs_EDGE),
75 if (shape.ShapeType() == TopAbs_WIRE)
76 return opencascade::handle<BRepAdaptor_HCompCurve>(
77 new BRepAdaptor_HCompCurve(TopoDS::Wire(shape)));
80 else if (shape.ShapeType() == TopAbs_EDGE)
81 return opencascade::handle<BRepAdaptor_HCurve>(
82 new BRepAdaptor_HCurve(TopoDS::Edge(shape)));
85 return opencascade::handle<BRepAdaptor_HCurve>(
new BRepAdaptor_HCurve());
93 shape_length(
const TopoDS_Shape &sh)
95# if DEAL_II_OPENCASCADE_VERSION_GTE(7, 6, 0)
96 opencascade::handle<Adaptor3d_Curve> adapt = curve_adaptor(sh);
97 return GCPnts_AbscissaPoint::Length(*adapt);
100 opencascade::handle<Adaptor3d_HCurve> adapt = curve_adaptor(sh);
102 return GCPnts_AbscissaPoint::Length(adapt->GetCurve());
108 template <
int dim,
int spacedim>
110 const TopoDS_Shape &sh,
111 const double tolerance)
113 , tolerance(tolerance)
120 template <
int dim,
int spacedim>
121 std::unique_ptr<Manifold<dim, spacedim>>
124 return std::unique_ptr<Manifold<dim, spacedim>>(
130 template <
int dim,
int spacedim>
136 (void)surrounding_points;
139 for (
unsigned int i = 0; i < surrounding_points.size(); ++i)
141 .distance(surrounding_points[i]) <
142 std::max(tolerance * surrounding_points[i].norm(),
144 ExcPointNotOnManifold<spacedim>(surrounding_points[i]));
151 template <
int dim,
int spacedim>
153 const TopoDS_Shape &sh,
155 const double tolerance)
157 , direction(direction)
158 , tolerance(tolerance)
165 template <
int dim,
int spacedim>
166 std::unique_ptr<Manifold<dim, spacedim>>
169 return std::unique_ptr<Manifold<dim, spacedim>>(
175 template <
int dim,
int spacedim>
181 (void)surrounding_points;
184 for (
unsigned int i = 0; i < surrounding_points.size(); ++i)
186 .distance(surrounding_points[i]) <
187 std::max(tolerance * surrounding_points[i].norm(),
189 ExcPointNotOnManifold<spacedim>(surrounding_points[i]));
197 template <
int dim,
int spacedim>
199 const TopoDS_Shape &sh,
200 const double tolerance)
202 , tolerance(tolerance)
208 "NormalToMeshProjectionManifold needs a shape containing faces to operate."));
211 template <
int dim,
int spacedim>
212 std::unique_ptr<Manifold<dim, spacedim>>
215 return std::unique_ptr<Manifold<dim, spacedim>>(
222 template <
int spacedim>
224 internal_project_to_manifold(
const TopoDS_Shape &,
235 internal_project_to_manifold(
236 const TopoDS_Shape &sh,
237 const double tolerance,
241 constexpr int spacedim = 3;
242 TopoDS_Shape out_shape;
246 for (
const auto &
point : surrounding_points)
250 ExcPointNotOnManifold<spacedim>(
point));
254 switch (surrounding_points.size())
258 for (
const auto &
point : surrounding_points)
265 average_normal += std::get<1>(p_and_diff_forms);
268 average_normal /= 2.0;
271 average_normal.
norm() > 1e-4,
273 "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."));
275 Tensor<1, 3> T = surrounding_points[0] - surrounding_points[1];
277 average_normal = average_normal - (average_normal *
T) * T;
278 average_normal /= average_normal.
norm();
283 Tensor<1, 3> u = surrounding_points[1] - surrounding_points[0];
284 Tensor<1, 3> v = surrounding_points[2] - surrounding_points[0];
285 const double n1_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]};
290 u = surrounding_points[2] - surrounding_points[3];
291 v = surrounding_points[1] - surrounding_points[3];
292 const double n2_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]};
298 average_normal = (n1 + n2) / 2.0;
301 average_normal.
norm() > tolerance,
303 "Failed to refine cell: the normal estimated via the surrounding points turns out to be a null vector, making the projection direction undetermined."));
305 average_normal /= average_normal.
norm();
310 Tensor<1, 3> u = surrounding_points[1] - surrounding_points[0];
311 Tensor<1, 3> v = surrounding_points[2] - surrounding_points[0];
312 const double n1_coords[3] = {u[1] * v[2] - u[2] * v[1],
313 u[2] * v[0] - u[0] * v[2],
314 u[0] * v[1] - u[1] * v[0]};
317 u = surrounding_points[2] - surrounding_points[3];
318 v = surrounding_points[1] - surrounding_points[3];
319 const double n2_coords[3] = {u[1] * v[2] - u[2] * v[1],
320 u[2] * v[0] - u[0] * v[2],
321 u[0] * v[1] - u[1] * v[0]};
324 u = surrounding_points[4] - surrounding_points[7];
325 v = surrounding_points[6] - surrounding_points[7];
326 const double n3_coords[3] = {u[1] * v[2] - u[2] * v[1],
327 u[2] * v[0] - u[0] * v[2],
328 u[0] * v[1] - u[1] * v[0]};
331 u = surrounding_points[6] - surrounding_points[7];
332 v = surrounding_points[5] - surrounding_points[7];
333 const double n4_coords[3] = {u[1] * v[2] - u[2] * v[1],
334 u[2] * v[0] - u[0] * v[2],
335 u[0] * v[1] - u[1] * v[0]};
339 average_normal = (n1 + n2 + n3 + n4) / 4.0;
342 average_normal.
norm() > tolerance,
344 "Failed to refine cell: the normal estimated via the surrounding points turns out to be a null vector, making the projection direction undetermined."));
346 average_normal /= average_normal.
norm();
353 for (
unsigned int i = 0; i < surrounding_points.size(); ++i)
354 for (
unsigned int j = 0; j < surrounding_points.size(); ++j)
356 for (
unsigned int k = 0; k < surrounding_points.size(); ++k)
357 if (k != j && k != i)
360 surrounding_points[i] - surrounding_points[j];
362 surrounding_points[i] - surrounding_points[k];
363 const double n_coords[3] = {u[1] * v[2] - u[2] * v[1],
364 u[2] * v[0] - u[0] * v[2],
368 if (n1.norm() > tolerance)
371 if (average_normal.
norm() < tolerance)
375 auto dot_prod = n1 * average_normal;
380 average_normal += n1;
382 average_normal -= n1;
387 average_normal.
norm() > tolerance,
389 "Failed to compute a normal: the normal estimated via the surrounding points turns out to be a null vector, making the projection direction undetermined."));
390 average_normal = average_normal / average_normal.
norm();
400 template <
int dim,
int spacedim>
406 return internal_project_to_manifold(sh,
414 template <
int dim,
int spacedim>
417 const double tolerance)
423 , curve(curve_adaptor(sh))
424 , tolerance(tolerance)
425 , length(shape_length(sh))
432 template <
int dim,
int spacedim>
433 std::unique_ptr<Manifold<dim, spacedim>>
436 return std::unique_ptr<Manifold<dim, spacedim>>(
442 template <
int dim,
int spacedim>
448 ShapeAnalysis_Curve curve_analysis;
451 const double dist = curve_analysis.Project(
453 *curve,
point(space_point), tolerance, proj, t,
true);
455 curve->GetCurve(),
point(space_point), tolerance, proj, t,
true);
459 Assert(dist < tolerance * length,
460 ExcPointNotOnManifold<spacedim>(space_point));
462 return Point<1>(GCPnts_AbscissaPoint::Length(
464 *curve, curve->FirstParameter(), t));
466 curve->GetCurve(), curve->GetCurve().FirstParameter(), t));
472 template <
int dim,
int spacedim>
477# if DEAL_II_OPENCASCADE_VERSION_GTE(7, 6, 0)
478 GCPnts_AbscissaPoint AP(*curve, chart_point[0], curve->FirstParameter());
479 gp_Pnt P = curve->Value(AP.Parameter());
481 GCPnts_AbscissaPoint AP(curve->GetCurve(),
483 curve->GetCurve().FirstParameter());
484 gp_Pnt P = curve->GetCurve().Value(AP.Parameter());
487 return point<spacedim>(P);
490 template <
int dim,
int spacedim>
492 const double tolerance)
494 , tolerance(tolerance)
499 template <
int dim,
int spacedim>
500 std::unique_ptr<Manifold<dim, spacedim>>
503 return std::unique_ptr<Manifold<dim, spacedim>>(
509 template <
int dim,
int spacedim>
514 opencascade::handle<Geom_Surface> SurfToProj = BRep_Tool::Surface(face);
516 ShapeAnalysis_Surface projector(SurfToProj);
517 gp_Pnt2d proj_params = projector.ValueOfUV(
point(space_point), tolerance);
519 double u = proj_params.X();
520 double v = proj_params.Y();
525 template <
int dim,
int spacedim>
530 return ::OpenCASCADE::push_forward<spacedim>(face,
535 template <
int dim,
int spacedim>
542 opencascade::handle<Geom_Surface> surf = BRep_Tool::Surface(face);
546 surf->D1(chart_point[0], chart_point[1], q, Du, Dv);
555 "Expecting derivative along Z to be zero! Bailing out."));
563 "Expecting derivative along Z to be zero! Bailing out."));
567 template <
int dim,
int spacedim>
568 std::tuple<double, double, double, double>
571 Standard_Real umin, umax, vmin, vmax;
572 BRepTools::UVBounds(face, umin, umax, vmin, vmax);
573 return std::make_tuple(umin, umax, vmin, vmax);
581# 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 > &)