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