src/SMESHUtils/SMESH_ControlPnt.cxx

   1 // Copyright (C) 2007-2024  CEA, EDF, 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   // Iterate on triangles and subdivide them
 167   int nbTriangles = aTri->NbTriangles();
 168   thePoints.reserve( thePoints.size() + nbTriangles );
 169   for ( int i = 1; i <= nbTriangles; i++ )
 170   {
 171     const Poly_Triangle& aTriangle = aTri->Triangle(i);
 172     gp_Pnt p1 = aTri->Node(aTriangle.Value(1));
 173     gp_Pnt p2 = aTri->Node(aTriangle.Value(2));
 174     gp_Pnt p3 = aTri->Node(aTriangle.Value(3));
 175
 176     p1.Transform(aTrsf);
 177     p2.Transform(aTrsf);
 178     p3.Transform(aTrsf);
 179
 180     subdivideTriangle( p1, p2, p3, theSize, thePoints );
 181   }
 182 }
 183
 184 //================================================================================
 185 /*!
 186  * \brief Fills a vector of points with point samples approximately
 187  * \brief spaced with a given size
 188  */
 189 //================================================================================
 190
 191 void SMESHUtils::createPointsSampleFromSolid( const TopoDS_Solid&      theSolid,
 192                                               const double&            theSize,
 193                                               std::vector<ControlPnt>& thePoints )
 194 {
 195   // Compute the bounding box
 196   double Xmin, Ymin, Zmin, Xmax, Ymax, Zmax;
 197   Bnd_Box B;
 198   BRepBndLib::Add(theSolid, B);
 199   B.Get(Xmin, Ymin, Zmin, Xmax, Ymax, Zmax);
 200
 201   // Create the points
 202   double step = theSize;
 203
 204   for ( double x=Xmin; x-Xmax<Precision::Confusion(); x=x+step )
 205   {
 206     for ( double y=Ymin; y-Ymax<Precision::Confusion(); y=y+step )
 207     {
 208       // Step1 : generate the Zmin -> Zmax line
 209       gp_Pnt startPnt(x, y, Zmin);
 210       gp_Pnt endPnt(x, y, Zmax);
 211       gp_Vec aVec(startPnt, endPnt);
 212       gp_Lin aLine(startPnt, aVec);
 213       double endParam = Zmax - Zmin;
 214
 215       // Step2 : for each face of theSolid:
 216       std::set<double> intersections;
 217
 218       for ( TopExp_Explorer Ex( theSolid, TopAbs_FACE ); Ex.More(); Ex.Next() )
 219       {
 220         // check if there is an intersection
 221         IntCurvesFace_Intersector anIntersector(TopoDS::Face(Ex.Current()), Precision::Confusion());
 222         anIntersector.Perform(aLine, 0, endParam);
 223
 224         // get the intersection's parameter and store it
 225         int nbPoints = anIntersector.NbPnt();
 226         for ( int i = 0 ; i < nbPoints; i++ )
 227         {
 228           intersections.insert( anIntersector.WParameter(i+1) );
 229         }
 230       }
 231       // Step3 : go through the line chunk by chunk
 232       if ( intersections.size() > 1 )
 233       {
 234         std::set<double>::iterator intersectionsIterator=intersections.begin();
 235         double first = *intersectionsIterator;
 236         intersectionsIterator++;
 237         bool innerPoints = true;
 238         for ( ; intersectionsIterator!=intersections.end() ; intersectionsIterator++ )
 239         {
 240           double second = *intersectionsIterator;
 241           if ( innerPoints )
 242           {
 243             // If the last chunk was outside of the shape or this is the first chunk
 244             // add the points in the range [first, second] to the points vector
 245             double localStep = (second -first) / ceil( (second - first) / step );
 246             for ( double z = Zmin + first; z < Zmin + second; z = z + localStep )
 247             {
 248               thePoints.emplace_back( x, y, z, theSize );
 249             }
 250             thePoints.emplace_back( x, y, Zmin + second, theSize );
 251           }
 252           first = second;
 253           innerPoints = !innerPoints;
 254         }
 255       }
 256     }
 257   }
 258 }
 259
 260 //================================================================================
 261 /*!
 262  * \brief Subdivides a triangle until it reaches a certain size (recursive function)
 263  */
 264 //================================================================================
 265
 266 void SMESHUtils::subdivideTriangle( const gp_Pnt& p1,
 267                                     const gp_Pnt& p2,
 268                                     const gp_Pnt& p3,
 269                                     const double& theSize,
 270                                     std::vector<ControlPnt>& thePoints)
 271 {
 272   // Size threshold to stop subdividing
 273   // This value ensures that two control points are distant no more than 2*theSize
 274   // as shown below
 275   //
 276   // The greater distance D of the mass center M to each Edge is 1/3 * Median
 277   // and Median < sqrt(3/4) * a  where a is the greater side (by using Apollonius' thorem).
 278   // So D < 1/3 * sqrt(3/4) * a and if a < sqrt(3) * S then D < S/2
 279   // and the distance between two mass centers of two neighbouring triangles
 280   // sharing an edge is < 2 * 1/2 * S = S
 281   // If the traingles share a Vertex and no Edge the distance of the mass centers
 282   // to the Vertices is 2*D < S so the mass centers are distant of less than 2*S
 283
 284   double threshold = sqrt( 3. ) * theSize;
 285
 286   if ( p1.Distance(p2) > threshold ||
 287        p2.Distance(p3) > threshold ||
 288        p3.Distance(p1) > threshold )
 289     try
 290     {
 291       gp_Pnt midPoints[3];
 292       computePointsForSplitting( p1, p2, p3, midPoints );
 293
 294       subdivideTriangle( midPoints[0], midPoints[1], midPoints[2], theSize, thePoints );
 295       subdivideTriangle( midPoints[0], p2, midPoints[1], theSize, thePoints );
 296       subdivideTriangle( midPoints[2], midPoints[1], p3, theSize, thePoints );
 297       subdivideTriangle( p1, midPoints[0], midPoints[2], theSize, thePoints );
 298       return;
 299     }
 300     catch (...)
 301     {
 302     }
 303
 304   gp_Pnt massCenter = ( p1.XYZ() + p2.XYZ() + p3.XYZ() ) / 3.;
 305   thePoints.emplace_back( massCenter, theSize );
 306 }
 307
 308 //================================================================================
 309 /*!
 310  * \brief Returns the appropriate points for splitting a triangle
 311  * the tangency points of the incircle are used in order to have mostly
 312  * well-shaped sub-triangles
 313  */
 314 //================================================================================
 315
 316 void SMESHUtils::computePointsForSplitting( const gp_Pnt& p1,
 317                                             const gp_Pnt& p2,
 318                                             const gp_Pnt& p3,
 319                                             gp_Pnt midPoints[3])
 320 {
 321   //Change coordinates
 322   gp_Trsf Trsf_1;            // Identity transformation
 323   gp_Ax3 reference_system(gp::Origin(), gp::DZ(), gp::DX());   // OXY
 324
 325   gp_Vec Vx(p1, p3);
 326   gp_Vec Vaux(p1, p2);
 327   gp_Dir Dx(Vx);
 328   gp_Dir Daux(Vaux);
 329   gp_Dir Dz = Dx.Crossed(Daux);
 330   gp_Ax3 current_system(p1, Dz, Dx);
 331
 332   Trsf_1.SetTransformation( reference_system, current_system );
 333
 334   gp_Pnt A = p1.Transformed(Trsf_1);
 335   gp_Pnt B = p2.Transformed(Trsf_1);
 336   gp_Pnt C = p3.Transformed(Trsf_1);
 337
 338   double a =  B.Distance(C) ;
 339   double b =  A.Distance(C) ;
 340   double c =  B.Distance(A) ;
 341
 342   // Incenter coordinates
 343   // see http://mathworld.wolfram.com/Incenter.html
 344   double Xi = ( b*B.X() + c*C.X() ) / ( a + b + c );
 345   double Yi = ( b*B.Y() ) / ( a + b + c );
 346   gp_Pnt Center(Xi, Yi, 0);
 347
 348   // Calculate the tangency points of the incircle
 349   gp_Pnt T1 = tangencyPoint( A, B, Center);
 350   gp_Pnt T2 = tangencyPoint( B, C, Center);
 351   gp_Pnt T3 = tangencyPoint( C, A, Center);
 352
 353   midPoints[0] = T1.Transformed(Trsf_1.Inverted());
 354   midPoints[1] = T2.Transformed(Trsf_1.Inverted());
 355   midPoints[2] = T3.Transformed(Trsf_1.Inverted());
 356
 357   return;
 358 }
 359
 360 //================================================================================
 361 /*!
 362  * \brief Computes the tangency points of the circle of center Center with
 363  * \brief the straight line (p1 p2)
 364  */
 365 //================================================================================
 366
 367 gp_Pnt SMESHUtils::tangencyPoint(const gp_Pnt& p1,
 368                                  const gp_Pnt& p2,
 369                                  const gp_Pnt& Center)
 370 {
 371   double Xt = 0;
 372   double Yt = 0;
 373
 374   // The tangency point is the intersection of the straight line (p1 p2)
 375   // and the straight line (Center T) which is orthogonal to (p1 p2)
 376   if ( fabs(p1.X() - p2.X()) <= Precision::Confusion() )
 377   {
 378     Xt=p1.X();     // T is on (p1 p2)
 379     Yt=Center.Y(); // (Center T) is orthogonal to (p1 p2)
 380   }
 381   else if ( fabs(p1.Y() - p2.Y()) <= Precision::Confusion() )
 382   {
 383     Yt=p1.Y();     // T is on (p1 p2)
 384     Xt=Center.X(); // (Center T) is orthogonal to (p1 p2)
 385   }
 386   else
 387   {
 388     // First straight line coefficients (equation y=a*x+b)
 389     double a = (p2.Y() - p1.Y()) / (p2.X() - p1.X())  ;
 390     double b = p1.Y() - a*p1.X();         // p1 is on this straight line
 391
 392     // Second straight line coefficients (equation y=c*x+d)
 393     double c = -1 / a;                    // The 2 lines are orthogonal
 394     double d = Center.Y() - c*Center.X(); // Center is on this straight line
 395
 396     Xt = (d - b) / (a - c);
 397     Yt = a*Xt + b;
 398   }
 399
 400   return gp_Pnt( Xt, Yt, 0 );
 401 }