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