src/SMESHUtils/SMESH_ControlPnt.cxx

   1 // Copyright (C) 2007-2019  CEA/DEN, EDF R&D
   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   std::vector<gp_Pnt> computePointsForSplitting( const gp_Pnt& p1,
  62                                                  const gp_Pnt& p2,
  63                                                  const gp_Pnt& p3 );
  64   gp_Pnt tangencyPoint(const gp_Pnt& p1,
  65                        const gp_Pnt& p2,
  66                        const gp_Pnt& Center);
  67
  68 }
  69
  70 //================================================================================
  71 /*!
  72  * \brief Fills a vector of points from which a size map input file can be written
  73  */
  74 //================================================================================
  75
  76 void SMESHUtils::createControlPoints( const TopoDS_Shape&      theShape,
  77                                       const double&            theSize,
  78                                       std::vector<ControlPnt>& thePoints )
  79 {
  80   if ( theShape.ShapeType() == TopAbs_VERTEX )
  81   {
  82     gp_Pnt aPnt = BRep_Tool::Pnt( TopoDS::Vertex(theShape) );
  83     ControlPnt aControlPnt( aPnt, theSize );
  84     thePoints.push_back( aControlPnt );
  85   }
  86   if ( theShape.ShapeType() == TopAbs_EDGE )
  87   {
  88     createPointsSampleFromEdge( TopoDS::Edge( theShape ), theSize, thePoints );
  89   }
  90   else if ( theShape.ShapeType() == TopAbs_WIRE )
  91   {
  92     TopExp_Explorer Ex;
  93     for (Ex.Init(theShape,TopAbs_EDGE); Ex.More(); Ex.Next())
  94     {
  95       createPointsSampleFromEdge( TopoDS::Edge( Ex.Current() ), theSize, thePoints );
  96     }
  97   }
  98   else if ( theShape.ShapeType() ==  TopAbs_FACE )
  99   {
 100     createPointsSampleFromFace( TopoDS::Face( theShape ), theSize, thePoints );
 101   }
 102   else if ( theShape.ShapeType() ==  TopAbs_SOLID )
 103   {
 104     createPointsSampleFromSolid( TopoDS::Solid( theShape ), theSize, thePoints );
 105   }
 106   else if ( theShape.ShapeType() == TopAbs_COMPOUND )
 107   {
 108     TopoDS_Iterator it( theShape );
 109     for(; it.More(); it.Next())
 110     {
 111       createControlPoints( it.Value(), theSize, thePoints );
 112     }
 113   }
 114 }
 115
 116 //================================================================================
 117 /*!
 118  * \brief Fills a vector of points with point samples approximately
 119  * \brief spaced with a given size
 120  */
 121 //================================================================================
 122
 123 void SMESHUtils::createPointsSampleFromEdge( const TopoDS_Edge&       theEdge,
 124                                              const double&            theSize,
 125                                              std::vector<ControlPnt>& thePoints )
 126 {
 127   double step = theSize;
 128   double first, last;
 129   Handle( Geom_Curve ) aCurve = BRep_Tool::Curve( theEdge, first, last );
 130   GeomAdaptor_Curve C ( aCurve );
 131   GCPnts_UniformAbscissa DiscretisationAlgo(C, step , first, last, Precision::Confusion());
 132   int nbPoints = DiscretisationAlgo.NbPoints();
 133
 134   ControlPnt aPnt;
 135   aPnt.SetSize(theSize);
 136
 137   for ( int i = 1; i <= nbPoints; i++ )
 138   {
 139     double param = DiscretisationAlgo.Parameter( i );
 140     aCurve->D0( param, aPnt );
 141     thePoints.push_back( aPnt );
 142   }
 143 }
 144
 145 //================================================================================
 146 /*!
 147  * \brief Fills a vector of points with point samples approximately
 148  * \brief spaced with a given size
 149  */
 150 //================================================================================
 151
 152 void SMESHUtils::createPointsSampleFromFace( const TopoDS_Face&       theFace,
 153                                              const double&            theSize,
 154                                              std::vector<ControlPnt>& thePoints )
 155 {
 156   BRepMesh_IncrementalMesh M(theFace, 0.01, Standard_True);
 157   TopLoc_Location aLocation;
 158
 159   // Triangulate the face
 160   Handle(Poly_Triangulation) aTri = BRep_Tool::Triangulation (theFace, aLocation);
 161
 162   // Get the transformation associated to the face location
 163   gp_Trsf aTrsf = aLocation.Transformation();
 164
 165   // Get triangles
 166   int nbTriangles = aTri->NbTriangles();
 167   Poly_Array1OfTriangle triangles(1,nbTriangles);
 168   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   for(int i=1; i<=nbTriangles; i++)
 177   {
 178     Poly_Triangle aTriangle = triangles.Value(i);
 179     gp_Pnt p1 = nodes.Value(aTriangle.Value(1));
 180     gp_Pnt p2 = nodes.Value(aTriangle.Value(2));
 181     gp_Pnt p3 = nodes.Value(aTriangle.Value(3));
 182
 183     p1.Transform(aTrsf);
 184     p2.Transform(aTrsf);
 185     p3.Transform(aTrsf);
 186
 187     subdivideTriangle(p1, p2, p3, theSize, thePoints);
 188   }
 189 }
 190
 191 //================================================================================
 192 /*!
 193  * \brief Fills a vector of points with point samples approximately
 194  * \brief spaced with a given size
 195  */
 196 //================================================================================
 197
 198 void SMESHUtils::createPointsSampleFromSolid( const TopoDS_Solid&      theSolid,
 199                                               const double&            theSize,
 200                                               std::vector<ControlPnt>& thePoints )
 201 {
 202   // Compute the bounding box
 203   double Xmin, Ymin, Zmin, Xmax, Ymax, Zmax;
 204   Bnd_Box B;
 205   BRepBndLib::Add(theSolid, B);
 206   B.Get(Xmin, Ymin, Zmin, Xmax, Ymax, Zmax);
 207
 208   // Create the points
 209   double step = theSize;
 210
 211   for ( double x=Xmin; x-Xmax<Precision::Confusion(); x=x+step )
 212   {
 213     for ( double y=Ymin; y-Ymax<Precision::Confusion(); y=y+step )
 214     {
 215       // Step1 : generate the Zmin -> Zmax line
 216       gp_Pnt startPnt(x, y, Zmin);
 217       gp_Pnt endPnt(x, y, Zmax);
 218       gp_Vec aVec(startPnt, endPnt);
 219       gp_Lin aLine(startPnt, aVec);
 220       double endParam = Zmax - Zmin;
 221
 222       // Step2 : for each face of theSolid:
 223       std::set<double> intersections;
 224       std::set<double>::iterator it = intersections.begin();
 225
 226       TopExp_Explorer Ex;
 227       for (Ex.Init(theSolid,TopAbs_FACE); Ex.More(); Ex.Next())
 228       {
 229         // check if there is an intersection
 230         IntCurvesFace_Intersector anIntersector(TopoDS::Face(Ex.Current()), Precision::Confusion());
 231         anIntersector.Perform(aLine, 0, endParam);
 232
 233         // get the intersection's parameter and store it
 234         int nbPoints = anIntersector.NbPnt();
 235         for(int i = 0 ; i < nbPoints ; i++ )
 236         {
 237           it = intersections.insert( it, anIntersector.WParameter(i+1) );
 238         }
 239       }
 240       // Step3 : go through the line chunk by chunk
 241       if ( intersections.begin() != intersections.end() )
 242       {
 243         std::set<double>::iterator intersectionsIterator=intersections.begin();
 244         double first = *intersectionsIterator;
 245         intersectionsIterator++;
 246         bool innerPoints = true;
 247         for ( ; intersectionsIterator!=intersections.end() ; intersectionsIterator++ )
 248         {
 249           double second = *intersectionsIterator;
 250           if ( innerPoints )
 251           {
 252             // If the last chunk was outside of the shape or this is the first chunk
 253             // add the points in the range [first, second] to the points vector
 254             double localStep = (second -first) / ceil( (second - first) / step );
 255             for ( double z = Zmin + first; z < Zmin + second; z = z + localStep )
 256             {
 257               thePoints.push_back(ControlPnt( x, y, z, theSize ));
 258             }
 259             thePoints.push_back(ControlPnt( x, y, Zmin + second, theSize ));
 260           }
 261           first = second;
 262           innerPoints = !innerPoints;
 263         }
 264       }
 265     }
 266   }
 267 }
 268
 269 //================================================================================
 270 /*!
 271  * \brief Subdivides a triangle until it reaches a certain size (recursive function)
 272  */
 273 //================================================================================
 274
 275 void SMESHUtils::subdivideTriangle( const gp_Pnt& p1,
 276                                     const gp_Pnt& p2,
 277                                     const gp_Pnt& p3,
 278                                     const double& theSize,
 279                                     std::vector<ControlPnt>& thePoints)
 280 {
 281   // Size threshold to stop subdividing
 282   // This value ensures that two control points are distant no more than 2*theSize
 283   // as shown below
 284   //
 285   // The greater distance D of the mass center M to each Edge is 1/3 * Median
 286   // and Median < sqrt(3/4) * a  where a is the greater side (by using Apollonius' thorem).
 287   // So D < 1/3 * sqrt(3/4) * a and if a < sqrt(3) * S then D < S/2
 288   // and the distance between two mass centers of two neighbouring triangles
 289   // sharing an edge is < 2 * 1/2 * S = S
 290   // If the traingles share a Vertex and no Edge the distance of the mass centers
 291   // to the Vertices is 2*D < S so the mass centers are distant of less than 2*S
 292
 293   double threshold = sqrt( 3. ) * theSize;
 294
 295   if ( (p1.Distance(p2) > threshold ||
 296         p2.Distance(p3) > threshold ||
 297         p3.Distance(p1) > threshold))
 298   {
 299     std::vector<gp_Pnt> midPoints = computePointsForSplitting(p1, p2, p3);
 300
 301     subdivideTriangle( midPoints[0], midPoints[1], midPoints[2], theSize, thePoints );
 302     subdivideTriangle( midPoints[0], p2, midPoints[1], theSize, thePoints );
 303     subdivideTriangle( midPoints[2], midPoints[1], p3, theSize, thePoints );
 304     subdivideTriangle( p1, midPoints[0], midPoints[2], theSize, thePoints );
 305   }
 306   else
 307   {
 308     double x = (p1.X() + p2.X() + p3.X()) / 3 ;
 309     double y = (p1.Y() + p2.Y() + p3.Y()) / 3 ;
 310     double z = (p1.Z() + p2.Z() + p3.Z()) / 3 ;
 311
 312     ControlPnt massCenter( x ,y ,z, theSize );
 313     thePoints.push_back( massCenter );
 314   }
 315 }
 316
 317 //================================================================================
 318 /*!
 319  * \brief Returns the appropriate points for splitting a triangle
 320  * \brief the tangency points of the incircle are used in order to have mostly
 321  * \brief well-shaped sub-triangles
 322  */
 323 //================================================================================
 324
 325 std::vector<gp_Pnt> SMESHUtils::computePointsForSplitting( const gp_Pnt& p1,
 326                                                            const gp_Pnt& p2,
 327                                                            const gp_Pnt& p3 )
 328 {
 329   std::vector<gp_Pnt> midPoints;
 330   //Change coordinates
 331   gp_Trsf Trsf_1;            // Identity transformation
 332   gp_Ax3 reference_system(gp::Origin(), gp::DZ(), gp::DX());   // OXY
 333
 334   gp_Vec Vx(p1, p3);
 335   gp_Vec Vaux(p1, p2);
 336   gp_Dir Dx(Vx);
 337   gp_Dir Daux(Vaux);
 338   gp_Dir Dz = Dx.Crossed(Daux);
 339   gp_Ax3 current_system(p1, Dz, Dx);
 340
 341   Trsf_1.SetTransformation( reference_system, current_system );
 342
 343   gp_Pnt A = p1.Transformed(Trsf_1);
 344   gp_Pnt B = p2.Transformed(Trsf_1);
 345   gp_Pnt C = p3.Transformed(Trsf_1);
 346
 347   double a =  B.Distance(C) ;
 348   double b =  A.Distance(C) ;
 349   double c =  B.Distance(A) ;
 350
 351   // Incenter coordinates
 352   // see http://mathworld.wolfram.com/Incenter.html
 353   double Xi = ( b*B.X() + c*C.X() ) / ( a + b + c );
 354   double Yi = ( b*B.Y() ) / ( a + b + c );
 355   gp_Pnt Center(Xi, Yi, 0);
 356
 357   // Calculate the tangency points of the incircle
 358   gp_Pnt T1 = tangencyPoint( A, B, Center);
 359   gp_Pnt T2 = tangencyPoint( B, C, Center);
 360   gp_Pnt T3 = tangencyPoint( C, A, Center);
 361
 362   gp_Pnt p1_2 = T1.Transformed(Trsf_1.Inverted());
 363   gp_Pnt p2_3 = T2.Transformed(Trsf_1.Inverted());
 364   gp_Pnt p3_1 = T3.Transformed(Trsf_1.Inverted());
 365
 366   midPoints.push_back(p1_2);
 367   midPoints.push_back(p2_3);
 368   midPoints.push_back(p3_1);
 369
 370   return midPoints;
 371 }
 372
 373 //================================================================================
 374 /*!
 375  * \brief Computes the tangency points of the circle of center Center with
 376  * \brief the straight line (p1 p2)
 377  */
 378 //================================================================================
 379
 380 gp_Pnt SMESHUtils::tangencyPoint(const gp_Pnt& p1,
 381                                  const gp_Pnt& p2,
 382                                  const gp_Pnt& Center)
 383 {
 384   double Xt = 0;
 385   double Yt = 0;
 386
 387   // The tangency point is the intersection of the straight line (p1 p2)
 388   // and the straight line (Center T) which is orthogonal to (p1 p2)
 389   if ( fabs(p1.X() - p2.X()) <= Precision::Confusion() )
 390   {
 391     Xt=p1.X();     // T is on (p1 p2)
 392     Yt=Center.Y(); // (Center T) is orthogonal to (p1 p2)
 393   }
 394   else if ( fabs(p1.Y() - p2.Y()) <= Precision::Confusion() )
 395   {
 396     Yt=p1.Y();     // T is on (p1 p2)
 397     Xt=Center.X(); // (Center T) is orthogonal to (p1 p2)
 398   }
 399   else
 400   {
 401     // First straight line coefficients (equation y=a*x+b)
 402     double a = (p2.Y() - p1.Y()) / (p2.X() - p1.X())  ;
 403     double b = p1.Y() - a*p1.X();         // p1 is on this straight line
 404
 405     // Second straight line coefficients (equation y=c*x+d)
 406     double c = -1 / a;                    // The 2 lines are orthogonal
 407     double d = Center.Y() - c*Center.X(); // Center is on this straight line
 408
 409     Xt = (d - b) / (a - c);
 410     Yt = a*Xt + b;
 411   }
 412
 413   return gp_Pnt( Xt, Yt, 0 );
 414 }