deal.II version GIT relicensing-2659-g040196caa3 2025-02-18 14:20:01+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
utilities.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2014 - 2023 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#include <deal.II/base/config.h>
16
18
19#ifdef DEAL_II_WITH_OPENCASCADE
20
22# include <deal.II/base/point.h>
24
26
27# include <IGESControl_Controller.hxx>
28# include <IGESControl_Reader.hxx>
29# include <IGESControl_Writer.hxx>
30# include <STEPControl_Controller.hxx>
31# include <STEPControl_Reader.hxx>
32# include <STEPControl_Writer.hxx>
33# include <TopExp_Explorer.hxx>
34# include <TopoDS.hxx>
35# include <TopoDS_Edge.hxx>
36# include <TopoDS_Face.hxx>
37# include <TopoDS_Shape.hxx>
38
39# include <cstdio>
40# include <iostream>
41# if DEAL_II_OPENCASCADE_VERSION_GTE(7, 0, 0)
42# include <Standard_Transient.hxx>
43# else
44# include <Handle_Standard_Transient.hxx>
45# endif
46
47# include <BRepAdaptor_Curve.hxx>
48# if DEAL_II_OPENCASCADE_VERSION_GTE(7, 6, 0)
49# include <BRepAlgoAPI_Section.hxx>
50# else
51# include <BRepAdaptor_HCompCurve.hxx>
52# include <BRepAdaptor_HCurve.hxx>
53# include <BRepAlgo_Section.hxx>
54# endif
55# include <BRepAdaptor_Surface.hxx>
56# include <BRepBndLib.hxx>
57# include <BRepBuilderAPI_MakeEdge.hxx>
58# include <BRepBuilderAPI_Sewing.hxx>
59# include <BRepBuilderAPI_Transform.hxx>
60# include <BRepMesh_IncrementalMesh.hxx>
61# include <BRepTools.hxx>
62# include <BRep_Builder.hxx>
63# include <GCPnts_AbscissaPoint.hxx>
64# include <GeomAPI_Interpolate.hxx>
65# include <GeomAPI_ProjectPointOnCurve.hxx>
66# include <GeomAPI_ProjectPointOnSurf.hxx>
67# include <GeomConvert_CompCurveToBSplineCurve.hxx>
68# include <GeomLProp_SLProps.hxx>
69# include <Geom_BoundedCurve.hxx>
70# include <Geom_Plane.hxx>
71# include <IntCurvesFace_ShapeIntersector.hxx>
72# include <Poly_Triangulation.hxx>
73# include <ShapeAnalysis_Surface.hxx>
74# include <StlAPI_Reader.hxx>
75# include <StlAPI_Writer.hxx>
76# include <TColStd_HSequenceOfTransient.hxx>
77# include <TColStd_SequenceOfTransient.hxx>
78# include <TColgp_HArray1OfPnt.hxx>
79# include <TopLoc_Location.hxx>
80# include <gp_Lin.hxx>
81# include <gp_Pnt.hxx>
82# include <gp_Vec.hxx>
83
84# include <algorithm>
85# include <vector>
86
87
89
90namespace OpenCASCADE
91{
92 std::tuple<unsigned int, unsigned int, unsigned int>
94 {
96 unsigned int n_faces = 0, n_edges = 0, n_vertices = 0;
97 for (exp.Init(shape, TopAbs_FACE); exp.More(); exp.Next(), ++n_faces)
98 {
99 }
100 for (exp.Init(shape, TopAbs_EDGE); exp.More(); exp.Next(), ++n_edges)
101 {
102 }
103 for (exp.Init(shape, TopAbs_VERTEX); exp.More(); exp.Next(), ++n_vertices)
104 {
105 }
106 return std::tuple<unsigned int, unsigned int, unsigned int>(n_faces,
107 n_edges,
108 n_vertices);
109 }
110
111 void
113 std::vector<TopoDS_Face> &faces,
114 std::vector<TopoDS_Edge> &edges,
115 std::vector<TopoDS_Vertex> &vertices)
116 {
117 faces.resize(0);
118 edges.resize(0);
119 vertices.resize(0);
120
121 TopExp_Explorer exp;
122 for (exp.Init(shape, TopAbs_FACE); exp.More(); exp.Next())
123 {
124 faces.push_back(TopoDS::Face(exp.Current()));
125 }
126 for (exp.Init(shape, TopAbs_EDGE); exp.More(); exp.Next())
127 {
128 edges.push_back(TopoDS::Edge(exp.Current()));
129 }
130 for (exp.Init(shape, TopAbs_VERTEX); exp.More(); exp.Next())
131 {
132 vertices.push_back(TopoDS::Vertex(exp.Current()));
133 }
134 }
135
136
137 void
139 std::vector<TopoDS_Compound> &compounds,
140 std::vector<TopoDS_CompSolid> &compsolids,
141 std::vector<TopoDS_Solid> &solids,
142 std::vector<TopoDS_Shell> &shells,
143 std::vector<TopoDS_Wire> &wires)
144 {
145 compounds.resize(0);
146 compsolids.resize(0);
147 solids.resize(0);
148 shells.resize(0);
149 wires.resize(0);
150
151 TopExp_Explorer exp;
152 for (exp.Init(shape, TopAbs_COMPOUND); exp.More(); exp.Next())
153 {
154 compounds.push_back(TopoDS::Compound(exp.Current()));
155 }
156 for (exp.Init(shape, TopAbs_COMPSOLID); exp.More(); exp.Next())
157 {
158 compsolids.push_back(TopoDS::CompSolid(exp.Current()));
159 }
160 for (exp.Init(shape, TopAbs_SOLID); exp.More(); exp.Next())
161 {
162 solids.push_back(TopoDS::Solid(exp.Current()));
163 }
164 for (exp.Init(shape, TopAbs_SHELL); exp.More(); exp.Next())
165 {
166 shells.push_back(TopoDS::Shell(exp.Current()));
167 }
168 for (exp.Init(shape, TopAbs_WIRE); exp.More(); exp.Next())
169 {
170 wires.push_back(TopoDS::Wire(exp.Current()));
171 }
172 }
173
174 template <int spacedim>
175 gp_Pnt
177 {
178 switch (spacedim)
179 {
180 case 1:
181 return gp_Pnt(p[0], 0, 0);
182 case 2:
183 return gp_Pnt(p[0], p[1], 0);
184 case 3:
185 return gp_Pnt(p[0], p[1], p[2]);
186 }
188 return {};
189 }
190
191 template <int spacedim>
193 point(const gp_Pnt &p, const double tolerance)
194 {
195 (void)tolerance;
196 switch (spacedim)
197 {
198 case 1:
199 Assert(std::abs(p.Y()) <= tolerance,
201 "Cannot convert OpenCASCADE point to 1d if p.Y() != 0."));
202 Assert(std::abs(p.Z()) <= tolerance,
204 "Cannot convert OpenCASCADE point to 1d if p.Z() != 0."));
205 return Point<spacedim>(p.X());
206 case 2:
207 Assert(std::abs(p.Z()) <= tolerance,
209 "Cannot convert OpenCASCADE point to 2d if p.Z() != 0."));
210 return Point<spacedim>(p.X(), p.Y());
211 case 3:
212 return Point<spacedim>(p.X(), p.Y(), p.Z());
213 }
215 return {};
216 }
217
218 template <int dim>
219 bool
221 const Point<dim> &p2,
222 const Tensor<1, dim> &direction,
223 const double tolerance)
224 {
225 const double rel_tol =
226 std::max(tolerance, std::max(p1.norm(), p2.norm()) * tolerance);
227 if (direction.norm() > 0.0)
228 return (p1 * direction < p2 * direction - rel_tol);
229 else
230 for (int d = dim; d >= 0; --d)
231 if (p1[d] < p2[d] - rel_tol)
232 return true;
233 else if (p2[d] < p1[d] - rel_tol)
234 return false;
235
236 // If we got here, for all d, none of the conditions above was
237 // satisfied. The two points are equal up to tolerance
238 return false;
239 }
240
241
243 read_IGES(const std::string &filename, const double scale_factor)
244 {
247 stat = reader.ReadFile(filename.c_str());
248 AssertThrow(stat == IFSelect_RetDone, ExcMessage("Error in reading file!"));
249
252 reader.PrintCheckLoad(failsonly, mode);
253
254 Standard_Integer nRoots = reader.TransferRoots();
255 // selects all IGES entities (including non visible ones) in the
256 // file and puts them into a list called MyList,
257
258 AssertThrow(nRoots > 0, ExcMessage("Read nothing from file."));
259
260 // Handle IGES Scale here.
262 gp_Trsf scale;
263 scale.SetScale(Origin, scale_factor);
264
265 TopoDS_Shape sh = reader.OneShape();
267
268 return trans.Shape(); // this is the actual translation
269 }
270
271 void
272 write_IGES(const TopoDS_Shape &shape, const std::string &filename)
273 {
275 IGESControl_Writer ICW("MM", 0);
276 Standard_Boolean ok = ICW.AddShape(shape);
277 AssertThrow(ok, ExcMessage("Failed to add shape to IGES controller."));
278 ICW.ComputeModel();
279 Standard_Boolean OK = ICW.Write(filename.c_str());
280 AssertThrow(OK, ExcMessage("Failed to write IGES file."));
281 }
282
283
285 read_STL(const std::string &filename)
286 {
288 TopoDS_Shape shape;
289 reader.Read(shape, filename.c_str());
290 return shape;
291 }
292
293
294 void
296 const std::string &filename,
297 const double deflection,
298 const bool sew_different_faces,
299 const double sewer_tolerance,
300 const bool is_relative,
301 const double angular_deflection,
302 const bool in_parallel)
303 {
305 std::vector<TopoDS_Vertex> vertices;
306 std::vector<TopoDS_Edge> edges;
307 std::vector<TopoDS_Face> faces;
308 OpenCASCADE::extract_geometrical_shapes(shape, faces, edges, vertices);
309 const bool mesh_is_present =
310 std::none_of(faces.begin(), faces.end(), [&Loc](const TopoDS_Face &face) {
311 Handle(Poly_Triangulation) theTriangulation =
312 BRep_Tool::Triangulation(face, Loc);
313 return theTriangulation.IsNull();
314 });
316 if (!mesh_is_present)
317 {
319 {
322 sewer.Perform();
323 shape_to_be_written = sewer.SewedShape();
324 }
325 else
326 shape_to_be_written = shape;
327 // BRepMesh_IncrementalMesh automatically calls the perform method to
328 // create the triangulation which is stored in the argument
329 // `shape_to_be_written`.
335 }
336
338
339# if DEAL_II_OPENCASCADE_VERSION_GTE(6, 9, 0)
340 // opencascade versions 6.9.0 onwards return an error status
341 const auto error = writer.Write(shape_to_be_written, filename.c_str());
342
343 // which is a custom type between 6.9.0 and 7.1.0
344# if !DEAL_II_OPENCASCADE_VERSION_GTE(7, 2, 0)
346 ExcMessage("Error writing STL from shape."));
347# else
348 // and a boolean from version 7.2.0 onwards
349 AssertThrow(error == true, ExcMessage("Error writing STL from shape."));
350# endif
351
352# else
353
354 // for opencascade versions 6.8.0 and older the return value is void
355 writer.Write(shape_to_be_written, filename.c_str());
356# endif
357 }
358
360 read_STEP(const std::string &filename, const double scale_factor)
361 {
364 stat = reader.ReadFile(filename.c_str());
365 AssertThrow(stat == IFSelect_RetDone, ExcMessage("Error in reading file!"));
366
369 reader.PrintCheckLoad(failsonly, mode);
370
371 Standard_Integer nRoots = reader.TransferRoots();
372 // selects all IGES entities (including non visible ones) in the
373 // file and puts them into a list called MyList,
374
375 AssertThrow(nRoots > 0, ExcMessage("Read nothing from file."));
376
377 // Handle STEP Scale here.
379 gp_Trsf scale;
380 scale.SetScale(Origin, scale_factor);
381
382 TopoDS_Shape sh = reader.OneShape();
384
385 return trans.Shape(); // this is the actual translation
386 }
387
388 void
389 write_STEP(const TopoDS_Shape &shape, const std::string &filename)
390 {
394 status = SCW.Transfer(shape, STEPControl_AsIs);
396 ExcMessage("Failed to add shape to STEP controller."));
397
398 status = SCW.Write(filename.c_str());
399
401 ExcMessage("Failed to write translated shape to STEP file."));
402 }
403
404 double
406 {
407 double tolerance = 0.0;
408
409 std::vector<TopoDS_Face> faces;
410 std::vector<TopoDS_Edge> edges;
411 std::vector<TopoDS_Vertex> vertices;
412
413 extract_geometrical_shapes(shape, faces, edges, vertices);
414
415 for (const auto &vertex : vertices)
416 tolerance = std::fmax(tolerance, BRep_Tool::Tolerance(vertex));
417
418 for (const auto &edge : edges)
419 tolerance = std::fmax(tolerance, BRep_Tool::Tolerance(edge));
420
421 for (const auto &face : faces)
422 tolerance = std::fmax(tolerance, BRep_Tool::Tolerance(face));
423
424
425 return tolerance;
426 }
427
430 const double c_x,
431 const double c_y,
432 const double c_z,
433 const double c,
434 const double /*tolerance*/)
435 {
436 Handle(Geom_Plane) plane = new Geom_Plane(c_x, c_y, c_z, c);
437# if DEAL_II_OPENCASCADE_VERSION_GTE(7, 6, 0)
439# else
441# endif
442 TopoDS_Shape edges = section.Shape();
443 return edges;
444 }
445
447 join_edges(const TopoDS_Shape &in_shape, const double tolerance)
448 {
450 const TopoDS_Shape &edges = in_shape;
451 std::vector<Handle_Geom_BoundedCurve> intersections;
455 gp_Pnt PIn(0.0, 0.0, 0.0);
456 gp_Pnt PFin(0.0, 0.0, 0.0);
457 gp_Pnt PMid(0.0, 0.0, 0.0);
460 while (edgeExplorer.More())
461 {
462 edge = TopoDS::Edge(edgeExplorer.Current());
463 Handle(Geom_Curve) curve = BRep_Tool::Curve(edge, L, First, Last);
464 intersections.push_back(Handle(Geom_BoundedCurve)::DownCast(curve));
465 edgeExplorer.Next();
466 }
467
468 // Now we build a single bspline out of all the geometrical
469 // curves, in Lexycographical order
470 unsigned int numIntersEdges = intersections.size();
471 Assert(numIntersEdges > 0, ExcMessage("No curves to process!"));
472
474
475 bool check = false, one_added = true, one_failed = true;
476 std::vector<bool> added(numIntersEdges, false);
477 added[0] = true;
478 while (one_added == true)
479 {
480 one_added = false;
481 one_failed = false;
482 for (unsigned int i = 1; i < numIntersEdges; ++i)
483 if (added[i] == false)
484 {
485 Handle(Geom_Curve) curve = intersections[i];
486 Handle(Geom_BoundedCurve) bcurve =
487 Handle(Geom_BoundedCurve)::DownCast(curve);
488 check = convert_bspline.Add(bcurve, tolerance, false, true, 0);
489 if (check ==
490 false) // If we failed, try again with the reversed curve
491 {
492 curve->Reverse();
493 Handle(Geom_BoundedCurve) bcurve =
494 Handle(Geom_BoundedCurve)::DownCast(curve);
495 check =
496 convert_bspline.Add(bcurve, tolerance, false, true, 0);
497 }
498 one_failed = one_failed || (check == false);
499 one_added = one_added || (check == true);
500 added[i] = check;
501 }
502 }
503
504 Assert(one_failed == false,
505 ExcMessage("Joining some of the Edges failed."));
506
507 Handle(Geom_Curve) bspline = convert_bspline.BSplineCurve();
508
510 return out_shape;
511 }
512
513 template <int dim>
516 const Point<dim> &origin,
517 const Tensor<1, dim> &direction,
518 const double tolerance)
519 {
520 // translating original Point<dim> to gp point
521
524 gp_Dir(direction[0],
525 dim > 1 ? direction[1] : 0,
526 dim > 2 ? direction[2] : 0));
527 gp_Lin line(gpaxis);
528
529 // destination point
530 gp_Pnt Pproj(0.0, 0.0, 0.0);
531
532 // we prepare now the surface for the projection we get the whole
533 // shape from the iges model
535 Inters.Load(in_shape, tolerance);
536
537 // Keep in mind: PerformNearest sounds pretty but DOESN'T WORK!!!
538 // The closest point must be found by hand
539 Inters.Perform(line, -RealLast(), +RealLast());
540 Assert(Inters.IsDone(), ExcMessage("Could not project point."));
541
542 double minDistance = 1e7;
544 for (int i = 0; i < Inters.NbPnt(); ++i)
545 {
546 const double distance = point(origin).Distance(Inters.Pnt(i + 1));
547 // cout<<"Point "<<i<<": "<<point(Inters.Pnt(i+1))<<" distance:
548 // "<<distance<<endl;
549 if (distance < minDistance)
550 {
551 minDistance = distance;
552 result = point<dim>(Inters.Pnt(i + 1));
553 }
554 }
555
556 return result;
557 }
558
559 template <int dim>
562 const Tensor<1, dim> &direction,
563 const bool closed,
564 const double tolerance)
565 {
566 unsigned int n_vertices = curve_points.size();
567
568 if (direction * direction > 0)
569 {
570 std::sort(curve_points.begin(),
571 curve_points.end(),
572 [&](const Point<dim> &p1, const Point<dim> &p2) {
573 return OpenCASCADE::point_compare(p1,
574 p2,
575 direction,
576 tolerance);
577 });
578 }
579
580 // set up array of vertices
581 Handle(TColgp_HArray1OfPnt) vertices =
582 new TColgp_HArray1OfPnt(1, n_vertices);
583 for (unsigned int vertex = 0; vertex < n_vertices; ++vertex)
584 {
585 vertices->SetValue(vertex + 1, point(curve_points[vertex]));
586 }
587
588
589 GeomAPI_Interpolate bspline_generator(vertices, closed, tolerance);
590 bspline_generator.Perform();
591 Assert((bspline_generator.IsDone()),
592 ExcMessage("Interpolated bspline generation failed"));
593
596 out_shape.Closed(closed);
597 return out_shape;
598 }
599
600
601
602 template <int spacedim>
603 std::vector<TopoDS_Edge>
606 const Mapping<2, spacedim> &mapping)
607
608 {
609 // store maps from global vertex index to pairs of global face indices
610 // and from global face index to pairs of global vertex indices
611 std::map<unsigned int, std::pair<unsigned int, unsigned int>> vert_to_faces;
612 std::map<unsigned int, std::pair<unsigned int, unsigned int>> face_to_verts;
613 std::map<unsigned int, bool> visited_faces;
614 std::map<unsigned int, Point<spacedim>> vert_to_point;
615
616 unsigned int face_index;
617
618 for (const auto &cell : triangulation.active_cell_iterators())
619 for (const unsigned int f : GeometryInfo<2>::face_indices())
620 if (cell->face(f)->at_boundary())
621 {
622 // get global face and vertex indices
623 face_index = cell->face(f)->index();
624 const unsigned int v0 = cell->face(f)->vertex_index(0);
625 const unsigned int v1 = cell->face(f)->vertex_index(1);
626 face_to_verts[face_index].first = v0;
627 face_to_verts[face_index].second = v1;
628 visited_faces[face_index] = false;
629
630 // extract mapped vertex locations
631 const auto verts = mapping.get_vertices(cell);
633 f, 0, true, false, false)];
635 f, 1, true, false, false)];
636
637 // distribute indices into maps
638 if (vert_to_faces.find(v0) == vert_to_faces.end())
639 {
640 vert_to_faces[v0].first = face_index;
641 }
642 else
643 {
644 vert_to_faces[v0].second = face_index;
645 }
646 if (vert_to_faces.find(v1) == vert_to_faces.end())
647 {
648 vert_to_faces[v1].first = face_index;
649 }
650 else
651 {
652 vert_to_faces[v1].second = face_index;
653 }
654 }
655
656 // run through maps in an orderly fashion, i.e., through the
657 // boundary in one cycle and add points to pointlist.
658 std::vector<TopoDS_Edge> interpolation_curves;
659 bool finished = (face_to_verts.empty());
660 face_index = finished ? 0 : face_to_verts.begin()->first;
661
662 while (finished == false)
663 {
664 const unsigned int start_point_index = face_to_verts[face_index].first;
665 unsigned int point_index = start_point_index;
666
667 // point_index and face_index always run together
668 std::vector<Point<spacedim>> pointlist;
669 do
670 {
671 visited_faces[face_index] = true;
673 pointlist.push_back(current_point);
674
675 // Get next point
676 if (face_to_verts[face_index].first != point_index)
677 point_index = face_to_verts[face_index].first;
678 else
679 point_index = face_to_verts[face_index].second;
680
681 // Get next face
682 if (vert_to_faces[point_index].first != face_index)
683 face_index = vert_to_faces[point_index].first;
684 else
685 face_index = vert_to_faces[point_index].second;
686 }
688
689 interpolation_curves.push_back(
691
692 finished = true;
693 for (const auto &f : visited_faces)
694 if (f.second == false)
695 {
696 face_index = f.first;
697 finished = false;
698 break;
699 }
700 }
702 }
703
704
705 template <int dim>
706 std::tuple<Point<dim>, TopoDS_Shape, double, double>
708 const Point<dim> &origin,
709 const double tolerance)
710 {
711 TopExp_Explorer exp;
713
714 double minDistance = 1e7;
715 gp_Pnt tmp_proj(0.0, 0.0, 0.0);
716
717 [[maybe_unused]] unsigned int counter = 0;
718 unsigned int face_counter = 0;
719
721 double u = 0;
722 double v = 0;
723
724 for (exp.Init(in_shape, TopAbs_FACE); exp.More(); exp.Next())
725 {
726 TopoDS_Face face = TopoDS::Face(exp.Current());
727
728 // the projection function needs a surface, so we obtain the
729 // surface upon which the face is defined
731
733 gp_Pnt2d proj_params = projector.ValueOfUV(point(origin), tolerance);
734
736
737 double distance = point<dim>(tmp_proj).distance(origin);
738 if (distance < minDistance)
739 {
740 minDistance = distance;
741 Pproj = tmp_proj;
742 out_shape = face;
743 u = proj_params.X();
744 v = proj_params.Y();
745 ++counter;
746 }
747 ++face_counter;
748 }
749
750 // face counter tells us if the shape contained faces: if it does, there is
751 // no need to loop on edges. Even if the closest point lies on the boundary
752 // of a parametric surface, we need in fact to retain the face and both u
753 // and v, if we want to use this method to retrieve the surface normal
754 if (face_counter == 0)
755 for (exp.Init(in_shape, TopAbs_EDGE); exp.More(); exp.Next())
756 {
757 TopoDS_Edge edge = TopoDS::Edge(exp.Current());
759 {
763
764 // the projection function needs a Curve, so we obtain the
765 // curve upon which the edge is defined
766 Handle(Geom_Curve) CurveToProj =
768
770 unsigned int num_proj_points = Proj.NbPoints();
771 if ((num_proj_points > 0) && (Proj.LowerDistance() < minDistance))
772 {
773 minDistance = Proj.LowerDistance();
774 Pproj = Proj.NearestPoint();
775 out_shape = edge;
776 u = Proj.LowerDistanceParameter();
777 ++counter;
778 }
779 }
780 }
781
782 Assert(counter > 0, ExcMessage("Could not find projection points."));
783 return std::tuple<Point<dim>, TopoDS_Shape, double, double>(
785 }
786
787
788 template <int dim>
791 const Point<dim> &origin,
792 const double tolerance)
793 {
794 std::tuple<Point<dim>, TopoDS_Shape, double, double> ref =
796 return std::get<0>(ref);
797 }
798
799 std::tuple<Point<3>, Tensor<1, 3>, double, double>
801 const Point<3> &origin,
802 const double tolerance)
803
804 {
805 std::tuple<Point<3>, TopoDS_Shape, double, double> shape_and_params =
807
809 double &u = std::get<2>(shape_and_params);
810 double &v = std::get<3>(shape_and_params);
811
812 // just a check here: the number of faces in out_shape must be 1, otherwise
813 // something is wrong
814 std::tuple<unsigned int, unsigned int, unsigned int> numbers =
816 (void)numbers;
817
818 Assert(
819 std::get<0>(numbers) > 0,
821 "Could not find normal: the shape containing the closest point has 0 faces."));
822 Assert(
823 std::get<0>(numbers) < 2,
825 "Could not find normal: the shape containing the closest point has more than 1 face."));
826
827
828 TopExp_Explorer exp;
829 exp.Init(out_shape, TopAbs_FACE);
830 TopoDS_Face face = TopoDS::Face(exp.Current());
831 return push_forward_and_differential_forms(face, u, v, tolerance);
832 }
833
834 template <int dim>
836 push_forward(const TopoDS_Shape &in_shape, const double u, const double v)
837 {
838 switch (in_shape.ShapeType())
839 {
840 case TopAbs_FACE:
841 {
843 return point<dim>(surf.Value(u, v));
844 }
845 case TopAbs_EDGE:
846 {
847 BRepAdaptor_Curve curve(TopoDS::Edge(in_shape));
848 return point<dim>(curve.Value(u));
849 }
850 default:
851 Assert(false, ExcUnsupportedShape());
852 }
853 return {};
854 }
855
856 std::tuple<Point<3>, Tensor<1, 3>, double, double>
857 push_forward_and_differential_forms(const TopoDS_Face &face,
858 const double u,
859 const double v,
860 const double /*tolerance*/)
861 {
863 GeomLProp_SLProps props(SurfToProj, u, v, 1, 1e-7);
864 gp_Pnt Value = props.Value();
865 Assert(props.IsNormalDefined(), ExcMessage("Normal is not well defined!"));
866 gp_Dir Normal = props.Normal();
867 Assert(props.IsCurvatureDefined(),
868 ExcMessage("Curvature is not well defined!"));
869 Standard_Real Min_Curvature = props.MinCurvature();
870 Standard_Real Max_Curvature = props.MaxCurvature();
871 Tensor<1, 3> normal = Point<3>(Normal.X(), Normal.Y(), Normal.Z());
872
873 // In the case your manifold changes from convex to concave or vice-versa
874 // the normal could jump from "inner" to "outer" normal.
875 // However, you should be able to change the normal sense preserving
876 // the manifold orientation:
877 if (face.Orientation() == TopAbs_REVERSED)
878 {
879 normal *= -1;
880 Min_Curvature *= -1;
881 Max_Curvature *= -1;
882 }
883
884 return std::tuple<Point<3>, Tensor<1, 3>, double, double>(point<3>(Value),
885 normal,
888 }
889
890
891
892 template <int spacedim>
893 void
894 create_triangulation(const TopoDS_Face &face,
896 {
898 const double u0 = surf.FirstUParameter();
899 const double u1 = surf.LastUParameter();
900 const double v0 = surf.FirstVParameter();
901 const double v1 = surf.LastVParameter();
902
903 std::vector<CellData<2>> cells;
904 std::vector<Point<spacedim>> vertices;
905 SubCellData t;
906
907 vertices.push_back(point<spacedim>(surf.Value(u0, v0)));
908 vertices.push_back(point<spacedim>(surf.Value(u1, v0)));
909 vertices.push_back(point<spacedim>(surf.Value(u0, v1)));
910 vertices.push_back(point<spacedim>(surf.Value(u1, v1)));
911
912 CellData<2> cell;
913 for (unsigned int i = 0; i < 4; ++i)
914 cell.vertices[i] = i;
915
916 cells.push_back(cell);
917 tria.create_triangulation(vertices, cells, t);
918 }
919
920// We don't build the .inst file if deal.II isn't configured
921// with GMSH, but doxygen doesn't know that and tries to find that
922// file anyway for parsing -- which then of course it fails on. So
923// exclude the following from doxygen consideration.
924# ifndef DOXYGEN
925# include "opencascade/utilities.inst"
926# endif
927
928} // namespace OpenCASCADE
929
931
932#endif
numbers::NumberTraits< Number >::real_type norm() const
friend class Tensor
Definition tensor.h:865
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:518
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:519
Point< 2 > first
Definition grid_out.cc:4629
const unsigned int v0
const unsigned int v1
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcUnsupportedShape()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
TopoDS_Edge interpolation_curve(std::vector< Point< dim > > &curve_points, const Tensor< 1, dim > &direction=Tensor< 1, dim >(), const bool closed=false, const double tolerance=1e-7)
Definition utilities.cc:561
std::tuple< Point< 3 >, Tensor< 1, 3 >, double, double > push_forward_and_differential_forms(const TopoDS_Face &face, const double u, const double v, const double tolerance=1e-7)
Definition utilities.cc:857
void extract_compound_shapes(const TopoDS_Shape &shape, std::vector< TopoDS_Compound > &compounds, std::vector< TopoDS_CompSolid > &compsolids, std::vector< TopoDS_Solid > &solids, std::vector< TopoDS_Shell > &shells, std::vector< TopoDS_Wire > &wires)
Definition utilities.cc:138
Point< dim > push_forward(const TopoDS_Shape &in_shape, const double u, const double v)
Definition utilities.cc:836
std::tuple< unsigned int, unsigned int, unsigned int > count_elements(const TopoDS_Shape &shape)
Definition utilities.cc:93
void write_STEP(const TopoDS_Shape &shape, const std::string &filename)
Definition utilities.cc:389
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)
Definition utilities.cc:800
void extract_geometrical_shapes(const TopoDS_Shape &shape, std::vector< TopoDS_Face > &faces, std::vector< TopoDS_Edge > &edges, std::vector< TopoDS_Vertex > &vertices)
Definition utilities.cc:112
TopoDS_Edge join_edges(const TopoDS_Shape &in_shape, const double tolerance=1e-7)
Definition utilities.cc:447
TopoDS_Shape read_STEP(const std::string &filename, const double scale_factor=1e-3)
Definition utilities.cc:360
bool point_compare(const Point< dim > &p1, const Point< dim > &p2, const Tensor< 1, dim > &direction=Tensor< 1, dim >(), const double tolerance=1e-10)
Definition utilities.cc:220
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:193
Point< dim > closest_point(const TopoDS_Shape &in_shape, const Point< dim > &origin, const double tolerance=1e-7)
Definition utilities.cc:790
std::tuple< Point< dim >, TopoDS_Shape, double, double > project_point_and_pull_back(const TopoDS_Shape &in_shape, const Point< dim > &origin, const double tolerance=1e-7)
Definition utilities.cc:707
Point< dim > line_intersection(const TopoDS_Shape &in_shape, const Point< dim > &origin, const Tensor< 1, dim > &direction, const double tolerance=1e-7)
Definition utilities.cc:515
void write_STL(const TopoDS_Shape &shape, const std::string &filename, const double deflection, const bool sew_different_faces=false, const double sewer_tolerance=1e-6, const bool is_relative=false, const double angular_deflection=0.5, const bool in_parallel=false)
Definition utilities.cc:295
void create_triangulation(const TopoDS_Face &face, Triangulation< 2, spacedim > &tria)
Definition utilities.cc:894
TopoDS_Shape read_STL(const std::string &filename)
Definition utilities.cc:285
TopoDS_Shape intersect_plane(const TopoDS_Shape &in_shape, const double c_x, const double c_y, const double c_z, const double c, const double tolerance=1e-7)
Definition utilities.cc:429
std::vector< TopoDS_Edge > create_curves_from_triangulation_boundary(const Triangulation< 2, spacedim > &triangulation, const Mapping< 2, spacedim > &mapping=StaticMappingQ1< 2, spacedim >::mapping)
Definition utilities.cc:604
double get_shape_tolerance(const TopoDS_Shape &shape)
Definition utilities.cc:405
void write_IGES(const TopoDS_Shape &shape, const std::string &filename)
Definition utilities.cc:272
TopoDS_Shape read_IGES(const std::string &filename, const double scale_factor=1e-3)
Definition utilities.cc:243
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
static unsigned int face_to_cell_vertices(const unsigned int face, const unsigned int vertex, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false)
static std_cxx20::ranges::iota_view< unsigned int, unsigned int > face_indices()