21#ifdef DEAL_II_WITH_OPENCASCADE
24# include <BRepAdaptor_CompCurve.hxx>
25# include <BRepAdaptor_Curve.hxx>
26# if !DEAL_II_OPENCASCADE_VERSION_GTE(7, 6, 0)
27# include <BRepAdaptor_HCompCurve.hxx>
28# include <BRepAdaptor_HCurve.hxx>
30# include <BRepTools.hxx>
31# include <BRep_Tool.hxx>
32# include <GCPnts_AbscissaPoint.hxx>
33# include <ShapeAnalysis_Curve.hxx>
34# include <ShapeAnalysis_Surface.hxx>
36# if !DEAL_II_OPENCASCADE_VERSION_GTE(7, 0, 0)
37# include <Handle_Adaptor3d_HCurve.hxx>
53# if DEAL_II_OPENCASCADE_VERSION_GTE(7, 6, 0)
54 Handle_Adaptor3d_Curve
55 curve_adaptor(
const TopoDS_Shape &shape)
57 Assert((shape.ShapeType() == TopAbs_WIRE) ||
58 (shape.ShapeType() == TopAbs_EDGE),
60 if (shape.ShapeType() == TopAbs_WIRE)
61 return Handle(BRepAdaptor_CompCurve)(
62 new BRepAdaptor_CompCurve(TopoDS::Wire(shape)));
63 else if (shape.ShapeType() == TopAbs_EDGE)
64 return Handle(BRepAdaptor_Curve)(
65 new BRepAdaptor_Curve(TopoDS::Edge(shape)));
68 return Handle(BRepAdaptor_Curve)(
new BRepAdaptor_Curve());
71 Handle_Adaptor3d_HCurve
72 curve_adaptor(
const TopoDS_Shape &shape)
74 Assert((shape.ShapeType() == TopAbs_WIRE) ||
75 (shape.ShapeType() == TopAbs_EDGE),
77 if (shape.ShapeType() == TopAbs_WIRE)
78 return Handle(BRepAdaptor_HCompCurve)(
79 new BRepAdaptor_HCompCurve(TopoDS::Wire(shape)));
80 else if (shape.ShapeType() == TopAbs_EDGE)
81 return Handle(BRepAdaptor_HCurve)(
82 new BRepAdaptor_HCurve(TopoDS::Edge(shape)));
85 return Handle(BRepAdaptor_HCurve)(
new BRepAdaptor_HCurve());
93 shape_length(
const TopoDS_Shape &sh)
95# if DEAL_II_OPENCASCADE_VERSION_GTE(7, 6, 0)
96 Handle_Adaptor3d_Curve adapt = curve_adaptor(sh);
97 return GCPnts_AbscissaPoint::Length(*adapt);
99 Handle_Adaptor3d_HCurve adapt = curve_adaptor(sh);
100 return GCPnts_AbscissaPoint::Length(adapt->GetCurve());
106 template <
int dim,
int spacedim>
108 const TopoDS_Shape &sh,
109 const double tolerance)
111 , tolerance(tolerance)
118 template <
int dim,
int spacedim>
119 std::unique_ptr<Manifold<dim, spacedim>>
122 return std::unique_ptr<Manifold<dim, spacedim>>(
128 template <
int dim,
int spacedim>
134 (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(), tolerance),
140 ExcPointNotOnManifold<spacedim>(surrounding_points[i]));
147 template <
int dim,
int spacedim>
149 const TopoDS_Shape & sh,
151 const double tolerance)
153 , direction(direction)
154 , tolerance(tolerance)
161 template <
int dim,
int spacedim>
162 std::unique_ptr<Manifold<dim, spacedim>>
165 return std::unique_ptr<Manifold<dim, spacedim>>(
171 template <
int dim,
int spacedim>
177 (void)surrounding_points;
179 for (
unsigned int i = 0; i < surrounding_points.size(); ++i)
181 .distance(surrounding_points[i]) <
182 std::max(tolerance * surrounding_points[i].norm(), tolerance),
183 ExcPointNotOnManifold<spacedim>(surrounding_points[i]));
191 template <
int dim,
int spacedim>
193 const TopoDS_Shape &sh,
194 const double tolerance)
196 , tolerance(tolerance)
202 "NormalToMeshProjectionManifold needs a shape containing faces to operate."));
205 template <
int dim,
int spacedim>
206 std::unique_ptr<Manifold<dim, spacedim>>
209 return std::unique_ptr<Manifold<dim, spacedim>>(
216 template <
int spacedim>
218 internal_project_to_manifold(
const TopoDS_Shape &,
229 internal_project_to_manifold(
230 const TopoDS_Shape & sh,
231 const double tolerance,
235 constexpr int spacedim = 3;
236 TopoDS_Shape out_shape;
239 for (
unsigned int i = 0; i < surrounding_points.size(); ++i)
242 .distance(surrounding_points[i]) <
245 ExcPointNotOnManifold<spacedim>(surrounding_points[i]));
249 switch (surrounding_points.size())
253 for (
unsigned int i = 0; i < surrounding_points.size(); ++i)
257 sh, surrounding_points[i], tolerance);
258 average_normal += std::get<1>(p_and_diff_forms);
261 average_normal /= 2.0;
264 average_normal.
norm() > 1e-4,
266 "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."));
268 Tensor<1, 3> T = surrounding_points[0] - surrounding_points[1];
270 average_normal = average_normal - (average_normal * T) * T;
271 average_normal /= average_normal.
norm();
276 Tensor<1, 3> u = surrounding_points[1] - surrounding_points[0];
277 Tensor<1, 3> v = surrounding_points[2] - surrounding_points[0];
278 const double n1_coords[3] = {u[1] * v[2] - u[2] * v[1],
279 u[2] * v[0] - u[0] * v[2],
280 u[0] * v[1] - u[1] * v[0]};
283 u = surrounding_points[2] - surrounding_points[3];
284 v = surrounding_points[1] - surrounding_points[3];
285 const double n2_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]};
291 average_normal = (n1 + n2) / 2.0;
294 average_normal.
norm() > tolerance,
296 "Failed to refine cell: the normal estimated via the surrounding points turns out to be a null vector, making the projection direction undetermined."));
298 average_normal /= average_normal.
norm();
303 Tensor<1, 3> u = surrounding_points[1] - surrounding_points[0];
304 Tensor<1, 3> v = surrounding_points[2] - surrounding_points[0];
305 const double n1_coords[3] = {u[1] * v[2] - u[2] * v[1],
306 u[2] * v[0] - u[0] * v[2],
307 u[0] * v[1] - u[1] * v[0]};
310 u = surrounding_points[2] - surrounding_points[3];
311 v = surrounding_points[1] - surrounding_points[3];
312 const double n2_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[4] - surrounding_points[7];
318 v = surrounding_points[6] - surrounding_points[7];
319 const double n3_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[6] - surrounding_points[7];
325 v = surrounding_points[5] - surrounding_points[7];
326 const double n4_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]};
332 average_normal = (n1 + n2 + n3 + n4) / 4.0;
335 average_normal.
norm() > tolerance,
337 "Failed to refine cell: the normal estimated via the surrounding points turns out to be a null vector, making the projection direction undetermined."));
339 average_normal /= average_normal.
norm();
346 for (
unsigned int i = 0; i < surrounding_points.size(); ++i)
347 for (
unsigned int j = 0; j < surrounding_points.size(); ++j)
349 for (
unsigned int k = 0; k < surrounding_points.size(); ++k)
350 if (k != j && k != i)
353 surrounding_points[i] - surrounding_points[j];
355 surrounding_points[i] - surrounding_points[k];
356 const double n_coords[3] = {u[1] * v[2] - u[2] * v[1],
357 u[2] * v[0] - u[0] * v[2],
361 if (n1.norm() > tolerance)
364 if (average_normal.
norm() < tolerance)
368 auto dot_prod = n1 * average_normal;
373 average_normal += n1;
375 average_normal -= n1;
380 average_normal.
norm() > tolerance,
382 "Failed to compute a normal: the normal estimated via the surrounding points turns out to be a null vector, making the projection direction undetermined."));
383 average_normal = average_normal / average_normal.
norm();
393 template <
int dim,
int spacedim>
399 return internal_project_to_manifold(sh,
407 template <
int dim,
int spacedim>
410 const double tolerance)
416 , curve(curve_adaptor(sh))
417 , tolerance(tolerance)
418 , length(shape_length(sh))
425 template <
int dim,
int spacedim>
426 std::unique_ptr<Manifold<dim, spacedim>>
429 return std::unique_ptr<Manifold<dim, spacedim>>(
435 template <
int dim,
int spacedim>
441 ShapeAnalysis_Curve curve_analysis;
444 const double dist = curve_analysis.Project(
446 *curve,
point(space_point), tolerance, proj, t,
true);
448 curve->GetCurve(),
point(space_point), tolerance, proj, t,
true);
452 Assert(dist < tolerance * length,
453 ExcPointNotOnManifold<spacedim>(space_point));
455 return Point<1>(GCPnts_AbscissaPoint::Length(
457 *curve, curve->FirstParameter(), t));
459 curve->GetCurve(), curve->GetCurve().FirstParameter(), t));
465 template <
int dim,
int spacedim>
470# if DEAL_II_OPENCASCADE_VERSION_GTE(7, 6, 0)
471 GCPnts_AbscissaPoint AP(*curve, chart_point[0], curve->FirstParameter());
472 gp_Pnt P = curve->Value(AP.Parameter());
474 GCPnts_AbscissaPoint AP(curve->GetCurve(),
476 curve->GetCurve().FirstParameter());
477 gp_Pnt P = curve->GetCurve().Value(AP.Parameter());
480 return point<spacedim>(P);
483 template <
int dim,
int spacedim>
485 const double tolerance)
487 , tolerance(tolerance)
492 template <
int dim,
int spacedim>
493 std::unique_ptr<Manifold<dim, spacedim>>
496 return std::unique_ptr<Manifold<dim, spacedim>>(
502 template <
int dim,
int spacedim>
507 Handle(Geom_Surface) SurfToProj = BRep_Tool::Surface(face);
509 ShapeAnalysis_Surface projector(SurfToProj);
510 gp_Pnt2d proj_params = projector.ValueOfUV(
point(space_point), tolerance);
512 double u = proj_params.X();
513 double v = proj_params.Y();
518 template <
int dim,
int spacedim>
523 return ::OpenCASCADE::push_forward<spacedim>(face,
528 template <
int dim,
int spacedim>
534 Handle(Geom_Surface) surf = BRep_Tool::Surface(face);
538 surf->D1(chart_point[0], chart_point[1], q, Du, Dv);
547 "Expecting derivative along Z to be zero! Bailing out."));
555 "Expecting derivative along Z to be zero! Bailing out."));
559 template <
int dim,
int spacedim>
560 std::tuple<double, double, double, double>
563 Standard_Real umin, umax, vmin, vmax;
564 BRepTools::UVBounds(face, umin, umax, vmin, vmax);
565 return std::make_tuple(umin, umax, vmin, vmax);
569# include "manifold_lib.inst"
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 & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcUnsupportedShape()
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
virtual std::unique_ptr< Manifold< dim, spacedim > > clone() const override
std::tuple< double, double, double, double > get_uv_bounds() const
virtual std::unique_ptr< Manifold< dim, spacedim > > clone() const override
virtual Point< spacedim > push_forward(const Point< 1 > &chart_point) const override
DirectionalProjectionManifold(const TopoDS_Shape &sh, const Tensor< 1, spacedim > &direction, 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
virtual Point< 2 > pull_back(const Point< spacedim > &space_point) const override
virtual Point< spacedim > push_forward(const Point< 2 > &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 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
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
NormalProjectionManifold(const TopoDS_Shape &sh, const double tolerance=1e-7)
virtual std::unique_ptr< Manifold< dim, spacedim > > clone() const override
NormalToMeshProjectionManifold(const TopoDS_Shape &sh, const double tolerance=1e-7)
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
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 > &)