WaveBEM: Unsteady Nonlinear Potential Flow Solver for Ship-Wave Interaction.
occ_normal_projection.cc
Go to the documentation of this file.
1 #include "occ_normal_projection.h"
2 #include "occ_utilities.h"
3 
4 #include <iostream>
5 #include <TopoDS.hxx>
6 #include <TopoDS_Shape.hxx>
7 #include <TopoDS_Shell.hxx>
8 #include <TopoDS_Face.hxx>
9 #include <BRepTools.hxx>
10 #include <TopTools_SequenceOfShape.hxx>
11 #include <Handle_Standard_Transient.hxx>
12 #include <TColStd_SequenceOfTransient.hxx>
13 #include <TColStd_HSequenceOfTransient.hxx>
14 #include <TopExp_Explorer.hxx>
15 #include <gp_Pnt.hxx>
16 #include <gp_Vec.hxx>
17 #include <GeomAPI_ProjectPointOnSurf.hxx>
18 #include <GeomAPI_ProjectPointOnCurve.hxx>
19 #include <Standard_Real.hxx>
20 #include <Standard_Integer.hxx>
21 #include <BRep_Tool.hxx>
22 #include <Geom_Surface.hxx>
23 #include <Geom_Plane.hxx>
24 #include <Prs3d_ShapeTool.hxx>
25 #include <Bnd_Box.hxx>
26 #include <gp_Ax3.hxx>
27 #include <gp_Pln.hxx>
28 #include <BRepBuilderAPI_MakeFace.hxx>
29 #include <BRepBuilderAPI_MakeEdge.hxx>
30 #include <GeomAPI_IntSS.hxx>
31 #include <TopoDS_Wire.hxx>
32 #include <BRepBuilderAPI_MakePolygon.hxx>
33 #include <Geom_Curve.hxx>
34 #include <Geom_BoundedCurve.hxx>
35 #include <Geom_TrimmedCurve.hxx>
36 #include <TColGeom_Array1OfCurve.hxx>
37 #include <GeomAPI_IntCS.hxx>
38 #include <BRepBuilderAPI_MakeEdge.hxx>
39 #include <GeomLib_Tool.hxx>
40 #include <GCPnts_AbscissaPoint.hxx>
41 #include <GeomConvert_CompCurveToBSplineCurve.hxx>
42 #include <ShapeAnalysis_Surface.hxx>
43 #include <GeomLProp_SLProps.hxx>
44 #include <TColgp_Array1OfPnt.hxx>
45 #include <TColgp_Array1OfPnt2d.hxx>
46 #include <BRepAdaptor_Curve.hxx>
47 #include <BRepAdaptor_Surface.hxx>
48 #include <vector>
49 #define FALSE 0
50 #define TRUE 1
51 
52 
53 #include <deal.II/grid/tria_iterator.h>
54 #include <deal.II/grid/tria_accessor.h>
55 
56 namespace OpenCascade
57 {
58 
59  template <int dim>
60  NormalProjection<dim>::NormalProjection(const TopoDS_Shape &sh) :
61  sh(sh)
62  {
63  if (dim == 1)
64  {
65  // Check that we have at least
66  // one edge.
67  TopExp_Explorer edgeExplorer(sh , TopAbs_EDGE);
68  unsigned int n_edges = 0;
69  while (edgeExplorer.More())
70  {
71  n_edges++;
72  edgeExplorer.Next();
73  }
74  AssertThrow(n_edges > 0, ExcMessage("We have no edges to process"));
75  }
76  else if (dim == 2)
77  {
78  // Check that we have at least
79  // one face.
80  TopExp_Explorer faceExplorer(sh , TopAbs_FACE);
81  unsigned int n_faces = 0;
82  while (faceExplorer.More())
83  {
84  n_faces++;
85  faceExplorer.Next();
86  }
87  AssertThrow(n_faces > 0, ExcMessage("We have no faces to process"));
88  }
89  else
90  AssertThrow(false, ExcNotImplemented());
91  }
92 
93 
94  template <>
95  void NormalProjection<1>::normal_projection
96  (Point<3> &projection, const Point<3> &origin) const
97  {
98  TopLoc_Location L = sh.Location();
99  TopLoc_Location L_inv = L.Inverted();
100  // translating original
101  // Point<dim> to gp point
102  gp_Pnt P0(origin(0),origin(1),origin(2));
103  P0.Transform(L_inv);
104  // destination point
105  gp_Pnt Pproj(0.0,0.0,0.0);
106  // we prepare now the surface
107  // for the projection we get
108  // the whole shape from the
109  // iges model
110 
111 // and here we loop on the faces of the shape
112  TopExp_Explorer edgeExplorer(sh , TopAbs_EDGE);
113  unsigned int tot_proj_points = 0;
114  double minDistance = 1e7;
115  TopoDS_Edge edge;
116  while (edgeExplorer.More())
117  {
118  edge = TopoDS::Edge(edgeExplorer.Current());
119  edge.Location(L);
120  BRepAdaptor_Curve AC(edge);
121  Standard_Real First;
122  Standard_Real Last;
123  // the projection
124  // function needs a
125  // surface, so we obtain
126  // the surface upon which
127  // the face is defined
128  Handle(Geom_Curve) CurveToProj = BRep_Tool::Curve(edge,L,First,Last);
129 
130  TopoDS_Shape shnew = BRepBuilderAPI_MakeEdge(CurveToProj);
131 
132  GeomAPI_ProjectPointOnCurve Proj(P0,CurveToProj);
133  unsigned int num_proj_points = Proj.NbPoints();
134  tot_proj_points+=num_proj_points;
135  if ((num_proj_points > 0) && (Proj.LowerDistance() < minDistance))
136  {
137  Pproj = Proj.NearestPoint();
138  Pproj.Transform(L);
139  minDistance = Proj.LowerDistance();
140  }
141  edgeExplorer.Next();
142  }
143 
144  Assert(tot_proj_points > 0,
145  ExcMessage("Point projection on curve in normal direction does not exist"));
146 
147 // translating destination point
148  projection(0) = Pproj.X();
149  projection(1) = Pproj.Y();
150  projection(2) = Pproj.Z();
151 
152  }
153 
154  template <>
155  void NormalProjection<2>::normal_projection
156  (Point<3> &projection, const Point<3> &origin) const
157  {
158  TopLoc_Location L = sh.Location();
159  TopLoc_Location L_inv = L.Inverted();
160  // translating original
161  // Point<dim> to gp point
162  gp_Pnt P0(origin(0),origin(1),origin(2));
163  P0.Transform(L_inv);
164  // destination point
165  gp_Pnt Pproj(0.0,0.0,0.0);
166  // we prepare now the surface
167  // for the projection we get
168  // the whole shape from the
169  // iges model
170 
171 
172  // and here we loop on the faces of the shape
173  unsigned int face_count=0;
174  TopExp_Explorer faceExplorer(sh, TopAbs_FACE);
175  double minDistance = 1e7;
176  TopoDS_Face face;
177  gp_Pnt face_proj(0.0,0.0,0.0);
178 
179  while (faceExplorer.More())
180  {
181  face = TopoDS::Face(faceExplorer.Current());
182  // the projection
183  // function needs a
184  // surface, so we obtain
185  // the surface upon which
186  // the face is defined
187  Handle(Geom_Surface) SurfToProj = BRep_Tool::Surface(face);
188 
189 
190  ShapeAnalysis_Surface projector(SurfToProj);
191  gp_Pnt2d proj_params = projector.ValueOfUV(P0, 1e-7);
192  SurfToProj->D0(proj_params.X(),proj_params.Y(),face_proj);
193  if (Pnt(face_proj).distance(origin) < minDistance)
194  {
195  minDistance = Pnt(face_proj).distance(origin);
196  Pproj = face_proj;
197  }
198  faceExplorer.Next();
199  ++face_count;
200  }
201  Pproj.Transform(L);
202 
203 
204 // translating destination point
205  projection = Pnt(Pproj);
206 
207 
208  }
209 
210 
211 
212  template <>
213  void NormalProjection<1>::normal_projection_and_diff_forms(Point<3> &/*projection*/,
214  Point<3> &/*normal*/,
215  double &/*mean_curvature*/,
216  const Point<3> &/*origin*/) const
217  {
218  Assert(0 > 1,
219  ExcMessage("Method is not implemented for dim = 1"));
220  }
221 
222 
223  template <>
224  void NormalProjection<2>::normal_projection_and_diff_forms(Point<3> &projection,
225  Point<3> &normal,
226  double &mean_curvature,
227  const Point<3> &origin) const
228  {
229  TopLoc_Location L = sh.Location();
230  TopLoc_Location L_inv = L.Inverted();
231  // translating original
232  // Point<dim> to gp point
233  gp_Pnt P0 = Pnt(origin);
234 //P0.Transform(L_inv);
235  // destination point
236  gp_Pnt Pproj(0.0,0.0,0.0);
237 
238  // we prepare now the surface
239  // for the projection we get
240  // the whole shape from the
241  // iges model
242 
243 
244 // and here we loop on the faces of the shape
245  unsigned int face_count=0;
246  TopExp_Explorer faceExplorer(sh, TopAbs_FACE);
247  double minDistance = 1e7;
248  TopoDS_Face face;
249  gp_Pnt face_proj(0.0,0.0,0.0);
250  gp_Dir Normal;
251  Standard_Real Mean_Curvature = 0;
252 
253  while (faceExplorer.More())
254  {
255  face = TopoDS::Face(faceExplorer.Current());
256  face.Location(L);
257  BRepAdaptor_Surface AF(face);
258  //cout<<"** "<<Pnt(AF.Value(AF.FirstUParameter(),AF.FirstVParameter()))<<endl;
259  //cout<<"** "<<Pnt(AF.Value(AF.LastUParameter(),AF.FirstVParameter()))<<endl;
260  //cout<<"** "<<Pnt(AF.Value(AF.LastUParameter(),AF.LastVParameter()))<<endl;
261  //cout<<"** "<<Pnt(AF.Value(AF.FirstUParameter(),AF.LastVParameter()))<<endl;
262  // the projection
263  // function needs a
264  // surface, so we obtain
265  // the surface upon which
266  // the face is defined
267 
268 
269  Handle(Geom_Surface) SurfToProj = BRep_Tool::Surface(face);
270  Standard_Real uufirst,uulast,vvfirst,vvlast;
271  SurfToProj->Bounds (uufirst, uulast, vvfirst, vvlast);
272  //cout<<"*#* "<<Pnt(SurfToProj->Value(uufirst,vvfirst))<<endl;
273  //cout<<"*#* "<<Pnt(SurfToProj->Value(uulast,vvfirst))<<endl;
274  //cout<<"*#* "<<Pnt(SurfToProj->Value(uulast,vvlast))<<endl;
275  //cout<<"*#* "<<Pnt(SurfToProj->Value(uufirst,vvlast))<<endl;
276  //cout<<Pnt(P0)<<endl;
277  ShapeAnalysis_Surface projector(SurfToProj);
278  if (projector.IsDegenerated(P0,1e-7) )
279  {
280  //cout<<"Gotcha!"<<endl;
281  Standard_Real ufirst,ulast,vfirst,vlast;
282  projector.Bounds (ufirst, ulast, vfirst, vlast);
283  //cout<<"Ubound "<<ufirst<<" "<<ulast<<endl;
284  //cout<<"Vbound "<<vfirst<<" "<<vlast<<endl;
285  gp_Pnt Pmid = SurfToProj->Value((ufirst+ulast)/2.0, (vfirst+vlast)/2.0);
286 
287  gp_Pnt P0Mod(P0.X()+(Pmid.X()-P0.X())/1000.0,
288  P0.Y()+(Pmid.Y()-P0.Y())/1000.0,
289  P0.Z()+(Pmid.Z()-P0.Z())/1000.0);
290  P0=P0Mod;
291  //cout<<"Point "<<Pnt(P0)<<endl;
292  }
293  gp_Pnt2d proj_params = projector.ValueOfUV(P0, 1e-7);
294  SurfToProj->D0(proj_params.X(),proj_params.Y(),face_proj);
295  if (Pnt(face_proj).distance(origin) < minDistance)
296  {
297  minDistance = Pnt(face_proj).distance(origin);
298  Pproj = face_proj;
299 
300  GeomLProp_SLProps props(SurfToProj, proj_params.X(), proj_params.Y(),1, 1e-7);
301  Normal = props.Normal();
302  Mean_Curvature = props.MeanCurvature();
303 
304  // adjusting normal orientation
305  if (face.Orientation()==TopAbs_REVERSED)
306  {
307  Normal.Reverse();
308  Mean_Curvature*=-1;
309  }
310  }
311  faceExplorer.Next();
312  ++face_count;
313  }
314  //Pproj.Transform(L);
315 
316 // translating destination point
317  projection = Pnt(Pproj);
318 
319 // translating normal vector
320  //Normal.Transform(L);
321  normal(0) = Normal.X();
322  normal(1) = Normal.Y();
323  normal(2) = Normal.Z();
324  // translating mean curvature
325  //cout<<"????? "<<normal<<endl;
326 
327  mean_curvature = double(Mean_Curvature);
328 
329  }
330 
331  template <int dim>
332  Point<3> NormalProjection<dim>::get_new_point_on_line
333  (const Triangulation< 2,3 >::line_iterator &line) const
334  {
335  Point<3> projected_point;
336  Point<3> source_point = StraightBoundary<2,3>::get_new_point_on_line(line);
337  normal_projection(projected_point, source_point);
338  return projected_point;
339  }
340 
341 
342  template <int dim>
343  Point<3> NormalProjection<dim>::get_new_point_on_quad
344  (const Triangulation< 2,3 >::quad_iterator &quad) const
345  {
346  AssertThrow(dim == 2, ExcMessage("Cannot project on quad on this"
347  " geometry type."));
348  Point<3> proj_point;
349  Point<3> init_point =
350  StraightBoundary<2,3>::get_new_point_on_quad(quad);
351  normal_projection(proj_point, init_point);
352  return proj_point;
353  }
354 
355 
356 
357  template class NormalProjection<1>;
358  template class NormalProjection<2>;
359 
360 }
OpenCascade::NormalProjection::NormalProjection
NormalProjection(const TopoDS_Shape &sh)
Definition: occ_normal_projection.cc:60
edge
Definition: geometry_flags.h:41
AssertThrow
#define AssertThrow(cond, exc)
occ_normal_projection.h
ExcMessage
static::ExceptionBase & ExcMessage(std::string arg1)
tria_accessor.h
Assert
#define Assert(cond, exc)
StraightBoundary::get_new_point_on_line
virtual Point< spacedim > get_new_point_on_line(const typename Triangulation< dim, spacedim >::line_iterator &line) const
OpenCascade::Pnt
dealii::Point< 3 > Pnt(const gp_Pnt &p)
Convert OpenCascade point into.
Definition: occ_utilities.h:76
OpenCascade::NormalProjection::normal_projection
void normal_projection(Point< 3 > &projection, const Point< 3 > &origin) const
tria_iterator.h
OpenCascade
We collect in this namespace all utilities which operate on OpenCascade entities which don't need cla...
Definition: occ_arclength_projection.cc:53
ExcNotImplemented
static::ExceptionBase & ExcNotImplemented()
OpenCascade::NormalProjection::normal_projection_and_diff_forms
void normal_projection_and_diff_forms(Point< 3 > &projection, Point< 3 > &normal, double &mean_curvature, const Point< 3 > &origin) const
Handle
Handle(Geom_Curve) NumericalTowingTank
Definition: numerical_towing_tank.cc:4028
StraightBoundary::get_new_point_on_quad
virtual Point< spacedim > get_new_point_on_quad(const typename Triangulation< dim, spacedim >::quad_iterator &quad) const
occ_utilities.h