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