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