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