1 // Copyright (C) 2007-2020 CEA/DEN, EDF R&D, OPEN CASCADE
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.
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.
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
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
20 // Author : Lioka RAZAFINDRAZAKA (CEA)
22 #include "SMESH_ControlPnt.hxx"
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>
40 #include <TopoDS_Edge.hxx>
41 #include <TopoDS_Face.hxx>
42 #include <TopoDS_Iterator.hxx>
43 #include <TopoDS_Solid.hxx>
47 #include <gp_Trsf.hxx>
54 // Some functions for surface sampling
55 void subdivideTriangle( const gp_Pnt& p1,
58 const double& theSize,
59 std::vector<ControlPnt>& thePoints );
61 void computePointsForSplitting( const gp_Pnt& p1,
65 gp_Pnt tangencyPoint(const gp_Pnt& p1,
67 const gp_Pnt& Center);
71 //================================================================================
73 * \brief Fills a vector of points from which a size map input file can be written
75 //================================================================================
77 void SMESHUtils::createControlPoints( const TopoDS_Shape& theShape,
78 const double& theSize,
79 std::vector<ControlPnt>& thePoints )
81 if ( theShape.ShapeType() == TopAbs_VERTEX )
83 gp_Pnt aPnt = BRep_Tool::Pnt( TopoDS::Vertex(theShape) );
84 ControlPnt aControlPnt( aPnt, theSize );
85 thePoints.push_back( aControlPnt );
87 if ( theShape.ShapeType() == TopAbs_EDGE )
89 createPointsSampleFromEdge( TopoDS::Edge( theShape ), theSize, thePoints );
91 else if ( theShape.ShapeType() == TopAbs_WIRE )
94 for (Ex.Init(theShape,TopAbs_EDGE); Ex.More(); Ex.Next())
96 createPointsSampleFromEdge( TopoDS::Edge( Ex.Current() ), theSize, thePoints );
99 else if ( theShape.ShapeType() == TopAbs_FACE )
101 createPointsSampleFromFace( TopoDS::Face( theShape ), theSize, thePoints );
103 else if ( theShape.ShapeType() == TopAbs_SOLID )
105 createPointsSampleFromSolid( TopoDS::Solid( theShape ), theSize, thePoints );
107 else if ( theShape.ShapeType() == TopAbs_COMPOUND )
109 TopoDS_Iterator it( theShape );
110 for(; it.More(); it.Next())
112 createControlPoints( it.Value(), theSize, thePoints );
117 //================================================================================
119 * \brief Fills a vector of points with point samples approximately
120 * \brief spaced with a given size
122 //================================================================================
124 void SMESHUtils::createPointsSampleFromEdge( const TopoDS_Edge& theEdge,
125 const double& theSize,
126 std::vector<ControlPnt>& thePoints )
128 double step = theSize;
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();
136 aPnt.SetSize(theSize);
138 for ( int i = 1; i <= nbPoints; i++ )
140 double param = DiscretisationAlgo.Parameter( i );
141 aCurve->D0( param, aPnt );
142 thePoints.push_back( aPnt );
146 //================================================================================
148 * \brief Fills a vector of points with point samples approximately
149 * \brief spaced with a given size
151 //================================================================================
153 void SMESHUtils::createPointsSampleFromFace( const TopoDS_Face& theFace,
154 const double& theSize,
155 std::vector<ControlPnt>& thePoints )
157 BRepMesh_IncrementalMesh M(theFace, 0.01, Standard_True);
158 TopLoc_Location aLocation;
160 // Triangulate the face
161 Handle(Poly_Triangulation) aTri = BRep_Tool::Triangulation (theFace, aLocation);
163 // Get the transformation associated to the face location
164 gp_Trsf aTrsf = aLocation.Transformation();
167 int nbTriangles = aTri->NbTriangles();
168 const Poly_Array1OfTriangle& triangles = aTri->Triangles();
171 int nbNodes = aTri->NbNodes();
172 TColgp_Array1OfPnt nodes(1,nbNodes);
173 nodes = aTri->Nodes();
175 // Iterate on triangles and subdivide them
176 thePoints.reserve( thePoints.size() + nbTriangles );
177 for ( int i = 1; i <= nbTriangles; i++ )
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));
188 subdivideTriangle( p1, p2, p3, theSize, thePoints );
192 //================================================================================
194 * \brief Fills a vector of points with point samples approximately
195 * \brief spaced with a given size
197 //================================================================================
199 void SMESHUtils::createPointsSampleFromSolid( const TopoDS_Solid& theSolid,
200 const double& theSize,
201 std::vector<ControlPnt>& thePoints )
203 // Compute the bounding box
204 double Xmin, Ymin, Zmin, Xmax, Ymax, Zmax;
206 BRepBndLib::Add(theSolid, B);
207 B.Get(Xmin, Ymin, Zmin, Xmax, Ymax, Zmax);
210 double step = theSize;
212 for ( double x=Xmin; x-Xmax<Precision::Confusion(); x=x+step )
214 for ( double y=Ymin; y-Ymax<Precision::Confusion(); y=y+step )
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;
223 // Step2 : for each face of theSolid:
224 std::set<double> intersections;
226 for ( TopExp_Explorer Ex( theSolid, TopAbs_FACE ); Ex.More(); Ex.Next() )
228 // check if there is an intersection
229 IntCurvesFace_Intersector anIntersector(TopoDS::Face(Ex.Current()), Precision::Confusion());
230 anIntersector.Perform(aLine, 0, endParam);
232 // get the intersection's parameter and store it
233 int nbPoints = anIntersector.NbPnt();
234 for ( int i = 0 ; i < nbPoints; i++ )
236 intersections.insert( anIntersector.WParameter(i+1) );
239 // Step3 : go through the line chunk by chunk
240 if ( intersections.size() > 1 )
242 std::set<double>::iterator intersectionsIterator=intersections.begin();
243 double first = *intersectionsIterator;
244 intersectionsIterator++;
245 bool innerPoints = true;
246 for ( ; intersectionsIterator!=intersections.end() ; intersectionsIterator++ )
248 double second = *intersectionsIterator;
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 )
256 thePoints.emplace_back( x, y, z, theSize );
258 thePoints.emplace_back( x, y, Zmin + second, theSize );
261 innerPoints = !innerPoints;
268 //================================================================================
270 * \brief Subdivides a triangle until it reaches a certain size (recursive function)
272 //================================================================================
274 void SMESHUtils::subdivideTriangle( const gp_Pnt& p1,
277 const double& theSize,
278 std::vector<ControlPnt>& thePoints)
280 // Size threshold to stop subdividing
281 // This value ensures that two control points are distant no more than 2*theSize
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
292 double threshold = sqrt( 3. ) * theSize;
294 if ( p1.Distance(p2) > threshold ||
295 p2.Distance(p3) > threshold ||
296 p3.Distance(p1) > threshold )
300 computePointsForSplitting( p1, p2, p3, midPoints );
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 );
312 gp_Pnt massCenter = ( p1.XYZ() + p2.XYZ() + p3.XYZ() ) / 3.;
313 thePoints.emplace_back( massCenter, theSize );
316 //================================================================================
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
322 //================================================================================
324 void SMESHUtils::computePointsForSplitting( const gp_Pnt& p1,
330 gp_Trsf Trsf_1; // Identity transformation
331 gp_Ax3 reference_system(gp::Origin(), gp::DZ(), gp::DX()); // OXY
337 gp_Dir Dz = Dx.Crossed(Daux);
338 gp_Ax3 current_system(p1, Dz, Dx);
340 Trsf_1.SetTransformation( reference_system, current_system );
342 gp_Pnt A = p1.Transformed(Trsf_1);
343 gp_Pnt B = p2.Transformed(Trsf_1);
344 gp_Pnt C = p3.Transformed(Trsf_1);
346 double a = B.Distance(C) ;
347 double b = A.Distance(C) ;
348 double c = B.Distance(A) ;
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);
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);
361 midPoints[0] = T1.Transformed(Trsf_1.Inverted());
362 midPoints[1] = T2.Transformed(Trsf_1.Inverted());
363 midPoints[2] = T3.Transformed(Trsf_1.Inverted());
368 //================================================================================
370 * \brief Computes the tangency points of the circle of center Center with
371 * \brief the straight line (p1 p2)
373 //================================================================================
375 gp_Pnt SMESHUtils::tangencyPoint(const gp_Pnt& p1,
377 const gp_Pnt& Center)
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() )
386 Xt=p1.X(); // T is on (p1 p2)
387 Yt=Center.Y(); // (Center T) is orthogonal to (p1 p2)
389 else if ( fabs(p1.Y() - p2.Y()) <= Precision::Confusion() )
391 Yt=p1.Y(); // T is on (p1 p2)
392 Xt=Center.X(); // (Center T) is orthogonal to (p1 p2)
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
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
404 Xt = (d - b) / (a - c);
408 return gp_Pnt( Xt, Yt, 0 );