src/BRepMesh/BRepMesh_DelaunayDeflectionControlMeshAlgo.hxx

   1 // Created on: 2016-07-07
   2 // Copyright (c) 2016 OPEN CASCADE SAS
   3 // Created by: Oleg AGASHIN
   4 //
   5 // This file is part of Open CASCADE Technology software library.
   6 //
   7 // This library is free software; you can redistribute it and/or modify it under
   8 // the terms of the GNU Lesser General Public License version 2.1 as published
   9 // by the Free Software Foundation, with special exception defined in the file
  10 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
  11 // distribution for complete text of the license and disclaimer of any warranty.
  12 //
  13 // Alternatively, this file may be used under the terms of Open CASCADE
  14 // commercial license or contractual agreement.
  15
  16 #ifndef _BRepMeshTools_DelaunayDeflectionControlMeshAlgo_HeaderFile
  17 #define _BRepMeshTools_DelaunayDeflectionControlMeshAlgo_HeaderFile
  18
  19 #include <BRepMesh_DelaunayNodeInsertionMeshAlgo.hxx>
  20 #include <BRepMesh_GeomTool.hxx>
  21 #include <GeomLib.hxx>
  22
  23 //! Extends node insertion Delaunay meshing algo in order to control
  24 //! deflection of generated trianges. Splits triangles failing the check.
  25 template<class RangeSplitter, class BaseAlgo>
  26 class BRepMesh_DelaunayDeflectionControlMeshAlgo : public BRepMesh_DelaunayNodeInsertionMeshAlgo<RangeSplitter, BaseAlgo>
  27 {
  28 private:
  29   // Typedef for OCCT RTTI
  30   typedef BRepMesh_DelaunayNodeInsertionMeshAlgo<RangeSplitter, BaseAlgo> DelaunayInsertionBaseClass;
  31
  32 public:
  33
  34   //! Constructor.
  35   BRepMesh_DelaunayDeflectionControlMeshAlgo()
  36     : myMaxSqDeflection(-1.),
  37       myIsAllDegenerated(Standard_False)
  38   {
  39   }
  40
  41   //! Destructor.
  42   virtual ~BRepMesh_DelaunayDeflectionControlMeshAlgo()
  43   {
  44   }
  45
  46 protected:
  47
  48   //! Perfroms processing of generated mesh. Generates surface nodes and inserts them into structure.
  49   virtual void postProcessMesh(BRepMesh_Delaun& theMesher) Standard_OVERRIDE
  50   {
  51     // Insert surface nodes.
  52     DelaunayInsertionBaseClass::postProcessMesh(theMesher);
  53
  54     if (this->getParameters().ControlSurfaceDeflection &&
  55         this->getStructure()->ElementsOfDomain().Extent() > 0)
  56     {
  57       optimizeMesh(theMesher);
  58     }
  59   }
  60
  61   //! Checks deviation of a mesh from geometrical surface.
  62   //! Inserts additional nodes in case of huge deviation.
  63   virtual void optimizeMesh(BRepMesh_Delaun& theMesher)
  64   {
  65     Handle(NCollection_IncAllocator) aTmpAlloc =
  66       new NCollection_IncAllocator(IMeshData::MEMORY_BLOCK_SIZE_HUGE);
  67
  68     myCouplesMap   = new IMeshData::MapOfOrientedEdges(3 * this->getStructure()->ElementsOfDomain().Extent(), aTmpAlloc);
  69     myControlNodes = new IMeshData::ListOfPnt2d(aTmpAlloc);
  70     myCircles      = &theMesher.Circles();
  71
  72     const Standard_Integer aIterationsNb = 11;
  73     Standard_Boolean isInserted = Standard_True;
  74     for (Standard_Integer aPass = 1; aPass <= aIterationsNb && isInserted && !myIsAllDegenerated; ++aPass)
  75     {
  76       // Reset stop condition
  77       myMaxSqDeflection = -1.;
  78       myIsAllDegenerated = Standard_True;
  79       myControlNodes->Clear();
  80
  81       if (this->getStructure()->ElementsOfDomain().Extent() < 1)
  82       {
  83         break;
  84       }
  85
  86       // Iterate on current triangles
  87       IMeshData::IteratorOfMapOfInteger aTriangleIt(this->getStructure()->ElementsOfDomain());
  88       for (; aTriangleIt.More(); aTriangleIt.Next())
  89       {
  90         const BRepMesh_Triangle& aTriangle = this->getStructure()->GetElement(aTriangleIt.Key());
  91         splitTriangleGeometry(aTriangle);
  92       }
  93
  94       isInserted = this->insertNodes(myControlNodes, theMesher);
  95     }
  96
  97     myCouplesMap.Nullify();
  98     myControlNodes.Nullify();
  99
 100     if (!(myMaxSqDeflection < 0.))
 101     {
 102       this->getDFace()->SetDeflection(Sqrt(myMaxSqDeflection));
 103     }
 104   }
 105
 106 private:
 107   //! Contains geometrical data related to node of triangle.
 108   struct TriangleNodeInfo
 109   {
 110     gp_XY            Point2d;
 111     gp_XYZ           Point;
 112     Standard_Boolean isFrontierLink;
 113   };
 114
 115   //! Functor computing deflection of a point from surface.
 116   class NormalDeviation
 117   {
 118   public:
 119     NormalDeviation(
 120       const gp_Pnt& theRefPnt,
 121       const gp_Vec& theNormal)
 122       : myRefPnt(theRefPnt),
 123         myNormal(theNormal)
 124     {
 125     }
 126
 127     Standard_Real SquareDeviation(const gp_Pnt& thePoint) const
 128     {
 129       const Standard_Real aDeflection = Abs(myNormal.Dot(gp_Vec(myRefPnt, thePoint)));
 130       return aDeflection * aDeflection;
 131     }
 132
 133   private:
 134
 135     NormalDeviation (const NormalDeviation& theOther);
 136
 137     void operator= (const NormalDeviation& theOther);
 138
 139   private:
 140
 141     const gp_Pnt& myRefPnt;
 142     const gp_Vec& myNormal;
 143   };
 144
 145   //! Functor computing deflection of a point on triangle link from surface.
 146   class LineDeviation
 147   {
 148   public:
 149
 150     LineDeviation(
 151       const gp_Pnt& thePnt1,
 152       const gp_Pnt& thePnt2)
 153       : myPnt1(thePnt1),
 154         myPnt2(thePnt2)
 155     {
 156     }
 157
 158     Standard_Real SquareDeviation(const gp_Pnt& thePoint) const
 159     {
 160       return BRepMesh_GeomTool::SquareDeflectionOfSegment(myPnt1, myPnt2, thePoint);
 161     }
 162
 163   private:
 164
 165     LineDeviation (const LineDeviation& theOther);
 166
 167     void operator= (const LineDeviation& theOther);
 168
 169   private:
 170     const gp_Pnt& myPnt1;
 171     const gp_Pnt& myPnt2;
 172   };
 173
 174   //! Returns nodes info of the given triangle.
 175   inline void getTriangleInfo(
 176     const BRepMesh_Triangle& theTriangle,
 177     const Standard_Integer (&theNodesIndices)[3],
 178     TriangleNodeInfo       (&theInfo)[3]) const
 179   {
 180     const Standard_Integer(&e)[3] = theTriangle.myEdges;
 181     for (Standard_Integer i = 0; i < 3; ++i)
 182     {
 183       const BRepMesh_Vertex& aVertex = this->getStructure()->GetNode(theNodesIndices[i]);
 184       theInfo[i].Point2d        = this->getRangeSplitter().Scale(aVertex.Coord(), Standard_False).XY();
 185       theInfo[i].Point          = this->getNodesMap()->Value(aVertex.Location3d()).XYZ();
 186       theInfo[i].isFrontierLink = (this->getStructure()->GetLink(e[i]).Movability() == BRepMesh_Frontier);
 187     }
 188   }
 189
 190   // Check geometry of the given triangle. If triangle does not suit specified deflection, inserts new point.
 191   void splitTriangleGeometry(const BRepMesh_Triangle& theTriangle)
 192   {
 193     if (theTriangle.Movability() != BRepMesh_Deleted)
 194     {
 195       Standard_Integer aNodexIndices[3];
 196       this->getStructure()->ElementNodes(theTriangle, aNodexIndices);
 197
 198       TriangleNodeInfo aNodesInfo[3];
 199       getTriangleInfo(theTriangle, aNodexIndices, aNodesInfo);
 200
 201       gp_Vec aNormal;
 202       gp_Vec aLinkVec[3];
 203       if (computeTriangleGeometry(aNodesInfo, aLinkVec, aNormal))
 204       {
 205         myIsAllDegenerated = Standard_False;
 206
 207         const gp_XY aCenter2d = (aNodesInfo[0].Point2d +
 208                                  aNodesInfo[1].Point2d +
 209                                  aNodesInfo[2].Point2d) / 3.;
 210
 211         usePoint(aCenter2d, NormalDeviation(aNodesInfo[0].Point, aNormal));
 212         splitLinks(aNodesInfo, aNodexIndices);
 213       }
 214     }
 215   }
 216
 217   //! Updates array of links vectors.
 218   //! @return False on degenerative triangle.
 219   inline Standard_Boolean computeTriangleGeometry(
 220     const TriangleNodeInfo(&theNodesInfo)[3],
 221     gp_Vec                (&theLinks)[3],
 222     gp_Vec                 &theNormal)
 223   {
 224     if (checkTriangleForDegenerativityAndGetLinks(theNodesInfo, theLinks))
 225     {
 226       if (checkTriangleArea2d(theNodesInfo))
 227       {
 228         if (computeNormal(theLinks[0], theLinks[1], theNormal))
 229         {
 230           return Standard_True;
 231         }
 232       }
 233     }
 234
 235     return Standard_False;
 236   }
 237
 238   //! Updates array of links vectors.
 239   //! @return False on degenerative triangle.
 240   inline Standard_Boolean checkTriangleForDegenerativityAndGetLinks(
 241     const TriangleNodeInfo (&theNodesInfo)[3],
 242     gp_Vec                 (&theLinks)[3])
 243   {
 244     const Standard_Real MinimalSqLength3d = 1.e-12;
 245     for (Standard_Integer i = 0; i < 3; ++i)
 246     {
 247       theLinks[i] = theNodesInfo[(i + 1) % 3].Point - theNodesInfo[i].Point;
 248       if (theLinks[i].SquareMagnitude() < MinimalSqLength3d)
 249       {
 250         return Standard_False;
 251       }
 252     }
 253
 254     return Standard_True;
 255   }
 256
 257   //! Checks area of triangle in parametric space for degenerativity.
 258   //! @return False on degenerative triangle.
 259   inline Standard_Boolean checkTriangleArea2d(
 260     const TriangleNodeInfo (&theNodesInfo)[3])
 261   {
 262     const gp_Vec2d aLink2d1(theNodesInfo[0].Point2d, theNodesInfo[1].Point2d);
 263     const gp_Vec2d aLink2d2(theNodesInfo[1].Point2d, theNodesInfo[2].Point2d);
 264
 265     const Standard_Real MinimalArea2d = 1.e-9;
 266     return (Abs(aLink2d1 ^ aLink2d2) > MinimalArea2d);
 267   }
 268
 269   //! Computes normal using two link vectors.
 270   //! @return True on success, False in case of normal of null magnitude.
 271   inline Standard_Boolean computeNormal(const gp_Vec& theLink1,
 272                                         const gp_Vec& theLink2,
 273                                         gp_Vec&       theNormal)
 274   {
 275     const gp_Vec aNormal(theLink1 ^ theLink2);
 276     if (aNormal.SquareMagnitude() > gp::Resolution())
 277     {
 278       theNormal = aNormal.Normalized();
 279       return Standard_True;
 280     }
 281
 282     return Standard_False;
 283   }
 284
 285   //! Computes deflection of midpoints of triangles links.
 286   //! @return True if point fits specified deflection.
 287   inline void splitLinks(
 288     const TriangleNodeInfo (&theNodesInfo)[3],
 289     const Standard_Integer (&theNodesIndices)[3])
 290   {
 291     // Check deflection at triangle links
 292     for (Standard_Integer i = 0; i < 3; ++i)
 293     {
 294       if (theNodesInfo[i].isFrontierLink)
 295       {
 296         continue;
 297       }
 298
 299       const Standard_Integer j = (i + 1) % 3;
 300       // Check if this link was already processed
 301       Standard_Integer aFirstVertex, aLastVertex;
 302       if (theNodesIndices[i] < theNodesIndices[j])
 303       {
 304         aFirstVertex = theNodesIndices[i];
 305         aLastVertex  = theNodesIndices[j];
 306       }
 307       else
 308       {
 309         aFirstVertex = theNodesIndices[j];
 310         aLastVertex  = theNodesIndices[i];
 311       }
 312
 313       if (myCouplesMap->Add(BRepMesh_OrientedEdge(aFirstVertex, aLastVertex)))
 314       {
 315         const gp_XY aMidPnt2d = (theNodesInfo[i].Point2d +
 316                                  theNodesInfo[j].Point2d) / 2.;
 317
 318         if (!usePoint (aMidPnt2d, LineDeviation (theNodesInfo[i].Point,
 319                                                  theNodesInfo[j].Point)))
 320         {
 321           if (!checkLinkEndsForAngularDeviation(theNodesInfo[i],
 322                                                 theNodesInfo[j],
 323                                                 aMidPnt2d))
 324           {
 325             myControlNodes->Append(aMidPnt2d);
 326           }
 327         }
 328       }
 329     }
 330   }
 331
 332   //! Checks the given point (located between the given nodes)
 333   //! for specified angular deviation.
 334   Standard_Boolean checkLinkEndsForAngularDeviation(const TriangleNodeInfo& theNodeInfo1,
 335                                                     const TriangleNodeInfo& theNodeInfo2,
 336                                                     const gp_XY& /*theMidPoint*/)
 337   {
 338     gp_Dir aNorm1, aNorm2;
 339     const Handle(Geom_Surface)& aSurf =
 340       this->getDFace()->GetSurface()->ChangeSurface().Surface().Surface();
 341
 342     if ((GeomLib::NormEstim(aSurf, theNodeInfo1.Point2d, Precision::Confusion(), aNorm1) == 0) &&
 343         (GeomLib::NormEstim(aSurf, theNodeInfo2.Point2d, Precision::Confusion(), aNorm2) == 0))
 344     {
 345       Standard_Real anAngle = aNorm1.Angle(aNorm2);
 346       if (anAngle > this->getParameters().AngleInterior)
 347         return Standard_False;
 348     }
 349 #if 0
 350     else if (GeomLib::NormEstim(aSurf, theMidPoint, Precision::Confusion(), aNorm1) != 0)
 351     {
 352       // It is better to consider the singular point as a node of triangulation.
 353       // However, it leads to hangs up meshing some faces (including faces with
 354       // degenerated edges). E.g. tests "mesh standard_incmesh Q6".
 355       // So, this code fragment is better to implement in the future.
 356       return Standard_False;
 357     }
 358 #endif
 359
 360     return Standard_True;
 361   }
 362
 363   //! Computes deflection of the given point and caches it for
 364   //! insertion in case if it overflows deflection.
 365   //! @return True if point has been cached for insertion.
 366   template<class DeflectionFunctor>
 367   inline Standard_Boolean usePoint(
 368     const gp_XY&             thePnt2d,
 369     const DeflectionFunctor& theDeflectionFunctor)
 370   {
 371     gp_Pnt aPnt;
 372     this->getDFace()->GetSurface()->D0(thePnt2d.X(), thePnt2d.Y(), aPnt);
 373     if (!checkDeflectionOfPointAndUpdateCache(thePnt2d, aPnt, theDeflectionFunctor.SquareDeviation(aPnt)))
 374     {
 375       myControlNodes->Append(thePnt2d);
 376       return Standard_True;
 377     }
 378
 379     return Standard_False;
 380   }
 381
 382   //! Checks the given point for specified linear deflection.
 383   //! Updates value of total mesh defleciton.
 384   Standard_Boolean checkDeflectionOfPointAndUpdateCache(
 385     const gp_XY&        thePnt2d,
 386     const gp_Pnt&       thePnt3d,
 387     const Standard_Real theSqDeflection)
 388   {
 389     if (theSqDeflection > myMaxSqDeflection)
 390     {
 391       myMaxSqDeflection = theSqDeflection;
 392     }
 393
 394     const Standard_Real aSqDeflection =
 395       this->getDFace()->GetDeflection() * this->getDFace()->GetDeflection();
 396     if (theSqDeflection < aSqDeflection)
 397     {
 398       return Standard_True;
 399     }
 400
 401     return rejectByMinSize(thePnt2d, thePnt3d);
 402   }
 403
 404   //! Checks the given node for
 405   Standard_Boolean rejectByMinSize(
 406     const gp_XY&  thePnt2d,
 407     const gp_Pnt& thePnt3d)
 408   {
 409     const Standard_Real aSqMinSize =
 410       this->getParameters().MinSize * this->getParameters().MinSize;
 411
 412     IMeshData::MapOfInteger aUsedNodes;
 413     IMeshData::ListOfInteger& aCirclesList =
 414       const_cast<BRepMesh_CircleTool&>(*myCircles).Select(
 415         this->getRangeSplitter().Scale(thePnt2d, Standard_True).XY());
 416
 417     IMeshData::ListOfInteger::Iterator aCircleIt(aCirclesList);
 418     for (; aCircleIt.More(); aCircleIt.Next())
 419     {
 420       const BRepMesh_Triangle& aTriangle = this->getStructure()->GetElement(aCircleIt.Value());
 421
 422       Standard_Integer aNodes[3];
 423       this->getStructure()->ElementNodes(aTriangle, aNodes);
 424
 425       for (Standard_Integer i = 0; i < 3; ++i)
 426       {
 427         if (!aUsedNodes.Contains(aNodes[i]))
 428         {
 429           aUsedNodes.Add(aNodes[i]);
 430           const BRepMesh_Vertex& aVertex = this->getStructure()->GetNode(aNodes[i]);
 431           const gp_Pnt& aPoint = this->getNodesMap()->Value(aVertex.Location3d());
 432
 433           if (thePnt3d.SquareDistance(aPoint) < aSqMinSize)
 434           {
 435             return Standard_True;
 436           }
 437         }
 438       }
 439     }
 440
 441     return Standard_False;
 442   }
 443
 444 private:
 445   Standard_Real                         myMaxSqDeflection;
 446   Standard_Boolean                      myIsAllDegenerated;
 447   Handle(IMeshData::MapOfOrientedEdges) myCouplesMap;
 448   Handle(IMeshData::ListOfPnt2d)        myControlNodes;
 449   const BRepMesh_CircleTool*            myCircles;
 450 };
 451
 452 #endif