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