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