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