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