]> SALOME platform Git repositories - modules/smesh.git/blob - src/SMESHUtils/SMESH_ControlPnt.cxx
Salome HOME
2e644c13566c588da137706d99bdc9fdb1977055
[modules/smesh.git] / src / SMESHUtils / SMESH_ControlPnt.cxx
1 // Copyright (C) 2007-2020  CEA/DEN, EDF R&D, OPEN CASCADE
2 //
3 // This library is free software; you can redistribute it and/or
4 // modify it under the terms of the GNU Lesser General Public
5 // License as published by the Free Software Foundation; either
6 // version 2.1 of the License, or (at your option) any later version.
7 //
8 // This library is distributed in the hope that it will be useful,
9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
11 // Lesser General Public License for more details.
12 //
13 // You should have received a copy of the GNU Lesser General Public
14 // License along with this library; if not, write to the Free Software
15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
16 //
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
18 //
19
20 // Author : Lioka RAZAFINDRAZAKA (CEA)
21
22 #include "SMESH_ControlPnt.hxx"
23
24 #include <BRepBndLib.hxx>
25 #include <BRepMesh_IncrementalMesh.hxx>
26 #include <BRep_Tool.hxx>
27 #include <Bnd_Box.hxx>
28 #include <GCPnts_UniformAbscissa.hxx>
29 #include <GeomAdaptor_Curve.hxx>
30 #include <Geom_Curve.hxx>
31 #include <IntCurvesFace_Intersector.hxx>
32 #include <Poly_Array1OfTriangle.hxx>
33 #include <Poly_Triangle.hxx>
34 #include <Poly_Triangulation.hxx>
35 #include <Precision.hxx>
36 #include <TColgp_Array1OfPnt.hxx>
37 #include <TopExp_Explorer.hxx>
38 #include <TopLoc_Location.hxx>
39 #include <TopoDS.hxx>
40 #include <TopoDS_Edge.hxx>
41 #include <TopoDS_Face.hxx>
42 #include <TopoDS_Iterator.hxx>
43 #include <TopoDS_Solid.hxx>
44 #include <gp_Ax3.hxx>
45 #include <gp_Dir.hxx>
46 #include <gp_Lin.hxx>
47 #include <gp_Trsf.hxx>
48 #include <gp_Vec.hxx>
49
50 #include <set>
51
52 namespace SMESHUtils
53 {
54   // Some functions for surface sampling
55   void subdivideTriangle( const gp_Pnt& p1,
56                           const gp_Pnt& p2,
57                           const gp_Pnt& p3,
58                           const double& theSize,
59                           std::vector<ControlPnt>& thePoints );
60
61   void computePointsForSplitting( const gp_Pnt& p1,
62                                   const gp_Pnt& p2,
63                                   const gp_Pnt& p3,
64                                   gp_Pnt midPoints[3]);
65   gp_Pnt tangencyPoint(const gp_Pnt& p1,
66                        const gp_Pnt& p2,
67                        const gp_Pnt& Center);
68
69 }
70
71 //================================================================================
72 /*!
73  * \brief Fills a vector of points from which a size map input file can be written
74  */
75 //================================================================================
76
77 void SMESHUtils::createControlPoints( const TopoDS_Shape&      theShape,
78                                       const double&            theSize,
79                                       std::vector<ControlPnt>& thePoints )
80 {
81   if ( theShape.ShapeType() == TopAbs_VERTEX )
82   {
83     gp_Pnt aPnt = BRep_Tool::Pnt( TopoDS::Vertex(theShape) );
84     ControlPnt aControlPnt( aPnt, theSize );
85     thePoints.push_back( aControlPnt );
86   }
87   if ( theShape.ShapeType() == TopAbs_EDGE )
88   {
89     createPointsSampleFromEdge( TopoDS::Edge( theShape ), theSize, thePoints );
90   }
91   else if ( theShape.ShapeType() == TopAbs_WIRE )
92   {
93     TopExp_Explorer Ex;
94     for (Ex.Init(theShape,TopAbs_EDGE); Ex.More(); Ex.Next())
95     {
96       createPointsSampleFromEdge( TopoDS::Edge( Ex.Current() ), theSize, thePoints );
97     }
98   }
99   else if ( theShape.ShapeType() ==  TopAbs_FACE )
100   {
101     createPointsSampleFromFace( TopoDS::Face( theShape ), theSize, thePoints );
102   }
103   else if ( theShape.ShapeType() ==  TopAbs_SOLID )
104   {
105     createPointsSampleFromSolid( TopoDS::Solid( theShape ), theSize, thePoints );
106   }
107   else if ( theShape.ShapeType() == TopAbs_COMPOUND )
108   {
109     TopoDS_Iterator it( theShape );
110     for(; it.More(); it.Next())
111     {
112       createControlPoints( it.Value(), theSize, thePoints );
113     }
114   }
115 }
116
117 //================================================================================
118 /*!
119  * \brief Fills a vector of points with point samples approximately
120  * \brief spaced with a given size
121  */
122 //================================================================================
123
124 void SMESHUtils::createPointsSampleFromEdge( const TopoDS_Edge&       theEdge,
125                                              const double&            theSize,
126                                              std::vector<ControlPnt>& thePoints )
127 {
128   double step = theSize;
129   double first, last;
130   Handle( Geom_Curve ) aCurve = BRep_Tool::Curve( theEdge, first, last );
131   GeomAdaptor_Curve C ( aCurve );
132   GCPnts_UniformAbscissa DiscretisationAlgo(C, step , first, last, Precision::Confusion());
133   int nbPoints = DiscretisationAlgo.NbPoints();
134
135   ControlPnt aPnt;
136   aPnt.SetSize(theSize);
137
138   for ( int i = 1; i <= nbPoints; i++ )
139   {
140     double param = DiscretisationAlgo.Parameter( i );
141     aCurve->D0( param, aPnt );
142     thePoints.push_back( aPnt );
143   }
144 }
145
146 //================================================================================
147 /*!
148  * \brief Fills a vector of points with point samples approximately
149  * \brief spaced with a given size
150  */
151 //================================================================================
152
153 void SMESHUtils::createPointsSampleFromFace( const TopoDS_Face&       theFace,
154                                              const double&            theSize,
155                                              std::vector<ControlPnt>& thePoints )
156 {
157   BRepMesh_IncrementalMesh M(theFace, 0.01, Standard_True);
158   TopLoc_Location aLocation;
159
160   // Triangulate the face
161   Handle(Poly_Triangulation) aTri = BRep_Tool::Triangulation (theFace, aLocation);
162
163   // Get the transformation associated to the face location
164   gp_Trsf aTrsf = aLocation.Transformation();
165
166   // Get triangles
167   int nbTriangles = aTri->NbTriangles();
168   const Poly_Array1OfTriangle& triangles = aTri->Triangles();
169
170   // GetNodes
171   int nbNodes = aTri->NbNodes();
172   TColgp_Array1OfPnt nodes(1,nbNodes);
173   nodes = aTri->Nodes();
174
175   // Iterate on triangles and subdivide them
176   thePoints.reserve( thePoints.size() + nbTriangles );
177   for ( int i = 1; i <= nbTriangles; i++ )
178   {
179     const Poly_Triangle& aTriangle = triangles.Value(i);
180     gp_Pnt p1 = nodes.Value(aTriangle.Value(1));
181     gp_Pnt p2 = nodes.Value(aTriangle.Value(2));
182     gp_Pnt p3 = nodes.Value(aTriangle.Value(3));
183
184     p1.Transform(aTrsf);
185     p2.Transform(aTrsf);
186     p3.Transform(aTrsf);
187
188     subdivideTriangle( p1, p2, p3, theSize, thePoints );
189   }
190 }
191
192 //================================================================================
193 /*!
194  * \brief Fills a vector of points with point samples approximately
195  * \brief spaced with a given size
196  */
197 //================================================================================
198
199 void SMESHUtils::createPointsSampleFromSolid( const TopoDS_Solid&      theSolid,
200                                               const double&            theSize,
201                                               std::vector<ControlPnt>& thePoints )
202 {
203   // Compute the bounding box
204   double Xmin, Ymin, Zmin, Xmax, Ymax, Zmax;
205   Bnd_Box B;
206   BRepBndLib::Add(theSolid, B);
207   B.Get(Xmin, Ymin, Zmin, Xmax, Ymax, Zmax);
208
209   // Create the points
210   double step = theSize;
211
212   for ( double x=Xmin; x-Xmax<Precision::Confusion(); x=x+step )
213   {
214     for ( double y=Ymin; y-Ymax<Precision::Confusion(); y=y+step )
215     {
216       // Step1 : generate the Zmin -> Zmax line
217       gp_Pnt startPnt(x, y, Zmin);
218       gp_Pnt endPnt(x, y, Zmax);
219       gp_Vec aVec(startPnt, endPnt);
220       gp_Lin aLine(startPnt, aVec);
221       double endParam = Zmax - Zmin;
222
223       // Step2 : for each face of theSolid:
224       std::set<double> intersections;
225
226       for ( TopExp_Explorer Ex( theSolid, TopAbs_FACE ); Ex.More(); Ex.Next() )
227       {
228         // check if there is an intersection
229         IntCurvesFace_Intersector anIntersector(TopoDS::Face(Ex.Current()), Precision::Confusion());
230         anIntersector.Perform(aLine, 0, endParam);
231
232         // get the intersection's parameter and store it
233         int nbPoints = anIntersector.NbPnt();
234         for ( int i = 0 ; i < nbPoints; i++ )
235         {
236           intersections.insert( anIntersector.WParameter(i+1) );
237         }
238       }
239       // Step3 : go through the line chunk by chunk
240       if ( intersections.size() > 1 )
241       {
242         std::set<double>::iterator intersectionsIterator=intersections.begin();
243         double first = *intersectionsIterator;
244         intersectionsIterator++;
245         bool innerPoints = true;
246         for ( ; intersectionsIterator!=intersections.end() ; intersectionsIterator++ )
247         {
248           double second = *intersectionsIterator;
249           if ( innerPoints )
250           {
251             // If the last chunk was outside of the shape or this is the first chunk
252             // add the points in the range [first, second] to the points vector
253             double localStep = (second -first) / ceil( (second - first) / step );
254             for ( double z = Zmin + first; z < Zmin + second; z = z + localStep )
255             {
256               thePoints.emplace_back( x, y, z, theSize );
257             }
258             thePoints.emplace_back( x, y, Zmin + second, theSize );
259           }
260           first = second;
261           innerPoints = !innerPoints;
262         }
263       }
264     }
265   }
266 }
267
268 //================================================================================
269 /*!
270  * \brief Subdivides a triangle until it reaches a certain size (recursive function)
271  */
272 //================================================================================
273
274 void SMESHUtils::subdivideTriangle( const gp_Pnt& p1,
275                                     const gp_Pnt& p2,
276                                     const gp_Pnt& p3,
277                                     const double& theSize,
278                                     std::vector<ControlPnt>& thePoints)
279 {
280   // Size threshold to stop subdividing
281   // This value ensures that two control points are distant no more than 2*theSize
282   // as shown below
283   //
284   // The greater distance D of the mass center M to each Edge is 1/3 * Median
285   // and Median < sqrt(3/4) * a  where a is the greater side (by using Apollonius' thorem).
286   // So D < 1/3 * sqrt(3/4) * a and if a < sqrt(3) * S then D < S/2
287   // and the distance between two mass centers of two neighbouring triangles
288   // sharing an edge is < 2 * 1/2 * S = S
289   // If the traingles share a Vertex and no Edge the distance of the mass centers
290   // to the Vertices is 2*D < S so the mass centers are distant of less than 2*S
291
292   double threshold = sqrt( 3. ) * theSize;
293
294   if ( p1.Distance(p2) > threshold ||
295        p2.Distance(p3) > threshold ||
296        p3.Distance(p1) > threshold )
297     try
298     {
299       gp_Pnt midPoints[3];
300       computePointsForSplitting( p1, p2, p3, midPoints );
301
302       subdivideTriangle( midPoints[0], midPoints[1], midPoints[2], theSize, thePoints );
303       subdivideTriangle( midPoints[0], p2, midPoints[1], theSize, thePoints );
304       subdivideTriangle( midPoints[2], midPoints[1], p3, theSize, thePoints );
305       subdivideTriangle( p1, midPoints[0], midPoints[2], theSize, thePoints );
306       return;
307     }
308     catch (...)
309     {
310     }
311
312   gp_Pnt massCenter = ( p1.XYZ() + p2.XYZ() + p3.XYZ() ) / 3.;
313   thePoints.emplace_back( massCenter, theSize );
314 }
315
316 //================================================================================
317 /*!
318  * \brief Returns the appropriate points for splitting a triangle
319  * the tangency points of the incircle are used in order to have mostly
320  * well-shaped sub-triangles
321  */
322 //================================================================================
323
324 void SMESHUtils::computePointsForSplitting( const gp_Pnt& p1,
325                                             const gp_Pnt& p2,
326                                             const gp_Pnt& p3,
327                                             gp_Pnt midPoints[3])
328 {
329   //Change coordinates
330   gp_Trsf Trsf_1;            // Identity transformation
331   gp_Ax3 reference_system(gp::Origin(), gp::DZ(), gp::DX());   // OXY
332
333   gp_Vec Vx(p1, p3);
334   gp_Vec Vaux(p1, p2);
335   gp_Dir Dx(Vx);
336   gp_Dir Daux(Vaux);
337   gp_Dir Dz = Dx.Crossed(Daux);
338   gp_Ax3 current_system(p1, Dz, Dx);
339
340   Trsf_1.SetTransformation( reference_system, current_system );
341
342   gp_Pnt A = p1.Transformed(Trsf_1);
343   gp_Pnt B = p2.Transformed(Trsf_1);
344   gp_Pnt C = p3.Transformed(Trsf_1);
345
346   double a =  B.Distance(C) ;
347   double b =  A.Distance(C) ;
348   double c =  B.Distance(A) ;
349
350   // Incenter coordinates
351   // see http://mathworld.wolfram.com/Incenter.html
352   double Xi = ( b*B.X() + c*C.X() ) / ( a + b + c );
353   double Yi = ( b*B.Y() ) / ( a + b + c );
354   gp_Pnt Center(Xi, Yi, 0);
355
356   // Calculate the tangency points of the incircle
357   gp_Pnt T1 = tangencyPoint( A, B, Center);
358   gp_Pnt T2 = tangencyPoint( B, C, Center);
359   gp_Pnt T3 = tangencyPoint( C, A, Center);
360
361   midPoints[0] = T1.Transformed(Trsf_1.Inverted());
362   midPoints[1] = T2.Transformed(Trsf_1.Inverted());
363   midPoints[2] = T3.Transformed(Trsf_1.Inverted());
364
365   return;
366 }
367
368 //================================================================================
369 /*!
370  * \brief Computes the tangency points of the circle of center Center with
371  * \brief the straight line (p1 p2)
372  */
373 //================================================================================
374
375 gp_Pnt SMESHUtils::tangencyPoint(const gp_Pnt& p1,
376                                  const gp_Pnt& p2,
377                                  const gp_Pnt& Center)
378 {
379   double Xt = 0;
380   double Yt = 0;
381
382   // The tangency point is the intersection of the straight line (p1 p2)
383   // and the straight line (Center T) which is orthogonal to (p1 p2)
384   if ( fabs(p1.X() - p2.X()) <= Precision::Confusion() )
385   {
386     Xt=p1.X();     // T is on (p1 p2)
387     Yt=Center.Y(); // (Center T) is orthogonal to (p1 p2)
388   }
389   else if ( fabs(p1.Y() - p2.Y()) <= Precision::Confusion() )
390   {
391     Yt=p1.Y();     // T is on (p1 p2)
392     Xt=Center.X(); // (Center T) is orthogonal to (p1 p2)
393   }
394   else
395   {
396     // First straight line coefficients (equation y=a*x+b)
397     double a = (p2.Y() - p1.Y()) / (p2.X() - p1.X())  ;
398     double b = p1.Y() - a*p1.X();         // p1 is on this straight line
399
400     // Second straight line coefficients (equation y=c*x+d)
401     double c = -1 / a;                    // The 2 lines are orthogonal
402     double d = Center.Y() - c*Center.X(); // Center is on this straight line
403
404     Xt = (d - b) / (a - c);
405     Yt = a*Xt + b;
406   }
407
408   return gp_Pnt( Xt, Yt, 0 );
409 }