19 #ifdef DEAL_II_WITH_OPENCASCADE
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>
50 Handle_Adaptor3d_HCurve
51 curve_adaptor(
const TopoDS_Shape &shape)
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)));
64 return Handle(BRepAdaptor_HCurve)(
new BRepAdaptor_HCurve());
71 shape_length(
const TopoDS_Shape &sh)
73 Handle_Adaptor3d_HCurve adapt = curve_adaptor(sh);
74 return GCPnts_AbscissaPoint::Length(adapt->GetCurve());
79 template <
int dim,
int spacedim>
81 const TopoDS_Shape &sh,
82 const double tolerance)
84 , tolerance(tolerance)
91 template <
int dim,
int spacedim>
92 std::unique_ptr<Manifold<dim, spacedim>>
95 return std::unique_ptr<Manifold<dim, spacedim>>(
101 template <
int dim,
int spacedim>
107 (void)surrounding_points;
109 for (
unsigned int i = 0; i < surrounding_points.size(); ++i)
112 std::max(tolerance * surrounding_points[i].
norm(), tolerance),
113 ExcPointNotOnManifold<spacedim>(surrounding_points[i]));
120 template <
int dim,
int spacedim>
122 const TopoDS_Shape & sh,
124 const double tolerance)
126 , direction(direction)
127 , tolerance(tolerance)
134 template <
int dim,
int spacedim>
135 std::unique_ptr<Manifold<dim, spacedim>>
138 return std::unique_ptr<Manifold<dim, spacedim>>(
144 template <
int dim,
int spacedim>
150 (void)surrounding_points;
152 for (
unsigned int i = 0; i < surrounding_points.size(); ++i)
155 std::max(tolerance * surrounding_points[i].
norm(), tolerance),
156 ExcPointNotOnManifold<spacedim>(surrounding_points[i]));
164 template <
int dim,
int spacedim>
166 const TopoDS_Shape &sh,
167 const double tolerance)
169 , tolerance(tolerance)
175 "NormalToMeshProjectionManifold needs a shape containing faces to operate."));
178 template <
int dim,
int spacedim>
179 std::unique_ptr<Manifold<dim, spacedim>>
182 return std::unique_ptr<Manifold<dim, spacedim>>(
189 template <
int spacedim>
191 internal_project_to_manifold(
const TopoDS_Shape &,
202 internal_project_to_manifold(
203 const TopoDS_Shape & sh,
204 const double tolerance,
208 constexpr
int spacedim = 3;
209 TopoDS_Shape out_shape;
212 for (
unsigned int i = 0; i < surrounding_points.size(); ++i)
215 .distance(surrounding_points[i]) <
218 ExcPointNotOnManifold<spacedim>(surrounding_points[i]));
222 switch (surrounding_points.size())
226 for (
unsigned int i = 0; i < surrounding_points.size(); ++i)
230 sh, surrounding_points[i], tolerance);
231 average_normal += std::get<1>(p_and_diff_forms);
234 average_normal /= 2.0;
237 average_normal.
norm() > 1
e-4,
239 "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."));
241 Tensor<1, 3> T = surrounding_points[0] - surrounding_points[1];
243 average_normal = average_normal - (average_normal *
T) *
T;
244 average_normal /= average_normal.
norm();
249 Tensor<1, 3> u = surrounding_points[1] - surrounding_points[0];
250 Tensor<1, 3> v = surrounding_points[2] - surrounding_points[0];
251 const double n1_coords[3] = {u[1] * v[2] - u[2] * v[1],
252 u[2] * v[0] - u[0] * v[2],
253 u[0] * v[1] - u[1] * v[0]};
256 u = surrounding_points[2] - surrounding_points[3];
257 v = surrounding_points[1] - surrounding_points[3];
258 const double n2_coords[3] = {u[1] * v[2] - u[2] * v[1],
259 u[2] * v[0] - u[0] * v[2],
260 u[0] * v[1] - u[1] * v[0]};
264 average_normal = (n1 + n2) / 2.0;
267 average_normal.
norm() > tolerance,
269 "Failed to refine cell: the normal estimated via the surrounding points turns out to be a null vector, making the projection direction undetermined."));
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]};
290 u = surrounding_points[4] - surrounding_points[7];
291 v = surrounding_points[6] - surrounding_points[7];
292 const double n3_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]};
297 u = surrounding_points[6] - surrounding_points[7];
298 v = surrounding_points[5] - surrounding_points[7];
299 const double n4_coords[3] = {u[1] * v[2] - u[2] * v[1],
300 u[2] * v[0] - u[0] * v[2],
301 u[0] * v[1] - u[1] * v[0]};
305 average_normal = (n1 + n2 + n3 + n4) / 4.0;
308 average_normal.
norm() > tolerance,
310 "Failed to refine cell: the normal estimated via the surrounding points turns out to be a null vector, making the projection direction undetermined."));
312 average_normal /= average_normal.
norm();
319 for (
unsigned int i = 0; i < surrounding_points.size(); ++i)
320 for (
unsigned int j = 0; j < surrounding_points.size(); ++j)
322 for (
unsigned int k = 0; k < surrounding_points.size(); ++k)
323 if (k != j && k != i)
326 surrounding_points[i] - surrounding_points[j];
328 surrounding_points[i] - surrounding_points[k];
329 const double n_coords[3] = {u[1] * v[2] - u[2] * v[1],
330 u[2] * v[0] - u[0] * v[2],
334 if (n1.norm() > tolerance)
337 if (average_normal.
norm() < tolerance)
341 auto dot_prod = n1 * average_normal;
346 average_normal += n1;
348 average_normal -= n1;
353 average_normal.
norm() > tolerance,
355 "Failed to compute a normal: the normal estimated via the surrounding points turns out to be a null vector, making the projection direction undetermined."));
356 average_normal = average_normal / average_normal.
norm();
366 template <
int dim,
int spacedim>
372 return internal_project_to_manifold(sh,
380 template <
int dim,
int spacedim>
383 const double tolerance)
389 , curve(curve_adaptor(sh))
390 , tolerance(tolerance)
391 , length(shape_length(sh))
398 template <
int dim,
int spacedim>
399 std::unique_ptr<Manifold<dim, spacedim>>
402 return std::unique_ptr<Manifold<dim, spacedim>>(
408 template <
int dim,
int spacedim>
414 ShapeAnalysis_Curve curve_analysis;
416 const double dist = curve_analysis.Project(
417 curve->GetCurve(),
point(space_point), tolerance, proj, t,
true);
418 Assert(dist < tolerance * length,
419 ExcPointNotOnManifold<spacedim>(space_point));
421 return Point<1>(GCPnts_AbscissaPoint::Length(
422 curve->GetCurve(), curve->GetCurve().FirstParameter(), t));
427 template <
int dim,
int spacedim>
432 GCPnts_AbscissaPoint AP(curve->GetCurve(),
434 curve->GetCurve().FirstParameter());
435 gp_Pnt P = curve->GetCurve().Value(AP.Parameter());
436 return point<spacedim>(P);
439 template <
int dim,
int spacedim>
441 const double tolerance)
443 , tolerance(tolerance)
448 template <
int dim,
int spacedim>
449 std::unique_ptr<Manifold<dim, spacedim>>
452 return std::unique_ptr<Manifold<dim, spacedim>>(
458 template <
int dim,
int spacedim>
463 Handle(Geom_Surface) SurfToProj = BRep_Tool::Surface(face);
465 ShapeAnalysis_Surface projector(SurfToProj);
466 gp_Pnt2d proj_params = projector.ValueOfUV(
point(space_point), tolerance);
468 double u = proj_params.X();
469 double v = proj_params.Y();
474 template <
int dim,
int spacedim>
479 return ::OpenCASCADE::push_forward<spacedim>(face,
484 template <
int dim,
int spacedim>
490 Handle(Geom_Surface) surf = BRep_Tool::Surface(face);
494 surf->D1(chart_point[0], chart_point[1], q, Du, Dv);
501 Assert(std::abs(Du.Z()) < tolerance,
503 "Expecting derivative along Z to be zero! Bailing out."));
509 Assert(std::abs(Dv.Z()) < tolerance,
511 "Expecting derivative along Z to be zero! Bailing out."));
515 template <
int dim,
int spacedim>
516 std::tuple<double, double, double, double>
519 Standard_Real umin, umax, vmin, vmax;
520 BRepTools::UVBounds(face, umin, umax, vmin, vmax);
521 return std::make_tuple(umin, umax, vmin, vmax);
525 # include "manifold_lib.inst"