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