Reference documentation for deal.II version Git 08727cc441 2020-07-02 15:45:42 -0400
\(\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  std::array<Point<spacedim>, GeometryInfo<2>::vertices_per_cell>
623  verts = mapping.get_vertices(cell);
624  vert_to_point[v0] = verts[GeometryInfo<2>::face_to_cell_vertices(
625  f, 0, true, false, false)];
626  vert_to_point[v1] = verts[GeometryInfo<2>::face_to_cell_vertices(
627  f, 1, true, false, false)];
628 
629  // distribute indices into maps
630  if (vert_to_faces.find(v0) == vert_to_faces.end())
631  {
632  vert_to_faces[v0].first = face_index;
633  }
634  else
635  {
636  vert_to_faces[v0].second = face_index;
637  }
638  if (vert_to_faces.find(v1) == vert_to_faces.end())
639  {
640  vert_to_faces[v1].first = face_index;
641  }
642  else
643  {
644  vert_to_faces[v1].second = face_index;
645  }
646  }
647 
648  // run through maps in an orderly fashion, i.e., through the
649  // boundary in one cycle and add points to pointlist.
650  std::vector<TopoDS_Edge> interpolation_curves;
651  bool finished = (face_to_verts.size() == 0);
652  face_index = finished ? 0 : face_to_verts.begin()->first;
653 
654  while (finished == false)
655  {
656  const unsigned int start_point_index = face_to_verts[face_index].first;
657  unsigned int point_index = start_point_index;
658 
659  // point_index and face_index always run together
660  std::vector<Point<spacedim>> pointlist;
661  do
662  {
663  visited_faces[face_index] = true;
664  auto current_point = vert_to_point[point_index];
665  pointlist.push_back(current_point);
666 
667  // Get next point
668  if (face_to_verts[face_index].first != point_index)
669  point_index = face_to_verts[face_index].first;
670  else
671  point_index = face_to_verts[face_index].second;
672 
673  // Get next face
674  if (vert_to_faces[point_index].first != face_index)
675  face_index = vert_to_faces[point_index].first;
676  else
677  face_index = vert_to_faces[point_index].second;
678  }
679  while (point_index != start_point_index);
680 
681  interpolation_curves.push_back(
682  interpolation_curve(pointlist, Tensor<1, spacedim>(), true));
683 
684  finished = true;
685  for (const auto &f : visited_faces)
686  if (f.second == false)
687  {
688  face_index = f.first;
689  finished = false;
690  break;
691  }
692  }
693  return interpolation_curves;
694  }
695 
696 
697  template <int dim>
698  std::tuple<Point<dim>, TopoDS_Shape, double, double>
699  project_point_and_pull_back(const TopoDS_Shape &in_shape,
700  const Point<dim> & origin,
701  const double tolerance)
702  {
703  TopExp_Explorer exp;
704  gp_Pnt Pproj = point(origin);
705 
706  double minDistance = 1e7;
707  gp_Pnt tmp_proj(0.0, 0.0, 0.0);
708 
709  unsigned int counter = 0;
710  unsigned int face_counter = 0;
711 
712  TopoDS_Shape out_shape;
713  double u = 0;
714  double v = 0;
715 
716  for (exp.Init(in_shape, TopAbs_FACE); exp.More(); exp.Next())
717  {
718  TopoDS_Face face = TopoDS::Face(exp.Current());
719 
720  // the projection function needs a surface, so we obtain the
721  // surface upon which the face is defined
722  Handle(Geom_Surface) SurfToProj = BRep_Tool::Surface(face);
723 
724  ShapeAnalysis_Surface projector(SurfToProj);
725  gp_Pnt2d proj_params = projector.ValueOfUV(point(origin), tolerance);
726 
727  SurfToProj->D0(proj_params.X(), proj_params.Y(), tmp_proj);
728 
729  double distance = point<dim>(tmp_proj).distance(origin);
730  if (distance < minDistance)
731  {
732  minDistance = distance;
733  Pproj = tmp_proj;
734  out_shape = face;
735  u = proj_params.X();
736  v = proj_params.Y();
737  ++counter;
738  }
739  ++face_counter;
740  }
741 
742  // face counter tells us if the shape contained faces: if it does, there is
743  // no need to loop on edges. Even if the closest point lies on the boundary
744  // of a parametric surface, we need in fact to retain the face and both u
745  // and v, if we want to use this method to retrieve the surface normal
746  if (face_counter == 0)
747  for (exp.Init(in_shape, TopAbs_EDGE); exp.More(); exp.Next())
748  {
749  TopoDS_Edge edge = TopoDS::Edge(exp.Current());
750  if (!BRep_Tool::Degenerated(edge))
751  {
752  TopLoc_Location L;
753  Standard_Real First;
754  Standard_Real Last;
755 
756  // the projection function needs a Curve, so we obtain the
757  // curve upon which the edge is defined
758  Handle(Geom_Curve) CurveToProj =
759  BRep_Tool::Curve(edge, L, First, Last);
760 
761  GeomAPI_ProjectPointOnCurve Proj(point(origin), CurveToProj);
762  unsigned int num_proj_points = Proj.NbPoints();
763  if ((num_proj_points > 0) && (Proj.LowerDistance() < minDistance))
764  {
765  minDistance = Proj.LowerDistance();
766  Pproj = Proj.NearestPoint();
767  out_shape = edge;
768  u = Proj.LowerDistanceParameter();
769  ++counter;
770  }
771  }
772  }
773 
774  Assert(counter > 0, ExcMessage("Could not find projection points."));
775  return std::tuple<Point<dim>, TopoDS_Shape, double, double>(
776  point<dim>(Pproj), out_shape, u, v);
777  }
778 
779 
780  template <int dim>
781  Point<dim>
782  closest_point(const TopoDS_Shape &in_shape,
783  const Point<dim> & origin,
784  const double tolerance)
785  {
786  std::tuple<Point<dim>, TopoDS_Shape, double, double> ref =
787  project_point_and_pull_back(in_shape, origin, tolerance);
788  return std::get<0>(ref);
789  }
790 
791  std::tuple<Point<3>, Tensor<1, 3>, double, double>
792  closest_point_and_differential_forms(const TopoDS_Shape &in_shape,
793  const Point<3> & origin,
794  const double tolerance)
795 
796  {
797  std::tuple<Point<3>, TopoDS_Shape, double, double> shape_and_params =
798  project_point_and_pull_back(in_shape, origin, tolerance);
799 
800  TopoDS_Shape &out_shape = std::get<1>(shape_and_params);
801  double & u = std::get<2>(shape_and_params);
802  double & v = std::get<3>(shape_and_params);
803 
804  // just a check here: the number of faces in out_shape must be 1, otherwise
805  // something is wrong
806  std::tuple<unsigned int, unsigned int, unsigned int> numbers =
807  count_elements(out_shape);
808  (void)numbers;
809 
810  Assert(
811  std::get<0>(numbers) > 0,
812  ExcMessage(
813  "Could not find normal: the shape containing the closest point has 0 faces."));
814  Assert(
815  std::get<0>(numbers) < 2,
816  ExcMessage(
817  "Could not find normal: the shape containing the closest point has more than 1 face."));
818 
819 
820  TopExp_Explorer exp;
821  exp.Init(out_shape, TopAbs_FACE);
822  TopoDS_Face face = TopoDS::Face(exp.Current());
823  return push_forward_and_differential_forms(face, u, v, tolerance);
824  }
825 
826  template <int dim>
827  Point<dim>
828  push_forward(const TopoDS_Shape &in_shape, const double u, const double v)
829  {
830  switch (in_shape.ShapeType())
831  {
832  case TopAbs_FACE:
833  {
834  BRepAdaptor_Surface surf(TopoDS::Face(in_shape));
835  return point<dim>(surf.Value(u, v));
836  }
837  case TopAbs_EDGE:
838  {
839  BRepAdaptor_Curve curve(TopoDS::Edge(in_shape));
840  return point<dim>(curve.Value(u));
841  }
842  default:
843  Assert(false, ExcUnsupportedShape());
844  }
845  return Point<dim>();
846  }
847 
848  std::tuple<Point<3>, Tensor<1, 3>, double, double>
849  push_forward_and_differential_forms(const TopoDS_Face &face,
850  const double u,
851  const double v,
852  const double /*tolerance*/)
853  {
854  Handle(Geom_Surface) SurfToProj = BRep_Tool::Surface(face);
855  GeomLProp_SLProps props(SurfToProj, u, v, 1, 1e-7);
856  gp_Pnt Value = props.Value();
857  Assert(props.IsNormalDefined(), ExcMessage("Normal is not well defined!"));
858  gp_Dir Normal = props.Normal();
859  Assert(props.IsCurvatureDefined(),
860  ExcMessage("Curvature is not well defined!"));
861  Standard_Real Min_Curvature = props.MinCurvature();
862  Standard_Real Max_Curvature = props.MaxCurvature();
863  Tensor<1, 3> normal = Point<3>(Normal.X(), Normal.Y(), Normal.Z());
864 
865  // In the case your manifold changes from convex to concave or viceversa
866  // the normal could jump from "inner" to "outer" normal.
867  // However, you should be able to change the normal sense preserving
868  // the manifold orientation:
869  if (face.Orientation() == TopAbs_REVERSED)
870  {
871  normal *= -1;
872  Min_Curvature *= -1;
873  Max_Curvature *= -1;
874  }
875 
876  return std::tuple<Point<3>, Tensor<1, 3>, double, double>(point<3>(Value),
877  normal,
878  Min_Curvature,
879  Max_Curvature);
880  }
881 
882 
883 
884  template <int spacedim>
885  void
886  create_triangulation(const TopoDS_Face & face,
888  {
889  BRepAdaptor_Surface surf(face);
890  const double u0 = surf.FirstUParameter();
891  const double u1 = surf.LastUParameter();
892  const double v0 = surf.FirstVParameter();
893  const double v1 = surf.LastVParameter();
894 
895  std::vector<CellData<2>> cells;
896  std::vector<Point<spacedim>> vertices;
897  SubCellData t;
898 
899  vertices.push_back(point<spacedim>(surf.Value(u0, v0)));
900  vertices.push_back(point<spacedim>(surf.Value(u1, v0)));
901  vertices.push_back(point<spacedim>(surf.Value(u0, v1)));
902  vertices.push_back(point<spacedim>(surf.Value(u1, v1)));
903 
904  CellData<2> cell;
905  for (unsigned int i = 0; i < 4; ++i)
906  cell.vertices[i] = i;
907 
908  cells.push_back(cell);
909  tria.create_triangulation(vertices, cells, t);
910  }
911 
912 # include "utilities.inst"
913 
914 } // namespace OpenCASCADE
915 
917 
918 #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:828
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:839
IteratorRange< active_cell_iterator > active_cell_iterators() const
Definition: tria.cc:12400
std::vector< unsigned int > vertices
virtual std::array< Point< spacedim >, GeometryInfo< dim >::vertices_per_cell > get_vertices(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const
Definition: mapping.cc:28
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:849
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:1513
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:886
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:10728
#define Assert(cond, exc)
Definition: exceptions.h:1403
Abstract base class for mapping classes.
Definition: mapping.h:301
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:363
Point< dim > closest_point(const TopoDS_Shape &in_shape, const Point< dim > &origin, const double tolerance=1e-7)
Definition: utilities.cc:782
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:302
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:699
Point< 2 > first
Definition: grid_out.cc:4352
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:362
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
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:792
T max(const T &t, const MPI_Comm &mpi_communicator)