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. BRepMesh_DelaunayDeflectionControlMeshAlgo()35 BRepMesh_DelaunayDeflectionControlMeshAlgo() 36 : myMaxSqDeflection(-1.), 37 mySqMinSize(-1.), 38 myIsAllDegenerated(Standard_False), 39 myCircles(NULL) 40 { 41 } 42 43 //! Destructor. ~BRepMesh_DelaunayDeflectionControlMeshAlgo()44 virtual ~BRepMesh_DelaunayDeflectionControlMeshAlgo() 45 { 46 } 47 48 protected: 49 50 //! Performs processing of generated mesh. Generates surface nodes and inserts them into structure. postProcessMesh(BRepMesh_Delaun & theMesher,const Message_ProgressRange & theRange)51 virtual void postProcessMesh (BRepMesh_Delaun& theMesher, 52 const Message_ProgressRange& theRange) Standard_OVERRIDE 53 { 54 Message_ProgressScope aPS(theRange, "Post process mesh", 2); 55 // Insert surface nodes. 56 DelaunayInsertionBaseClass::postProcessMesh (theMesher, aPS.Next()); 57 if (!aPS.More()) 58 { 59 return; 60 } 61 62 if (this->getParameters().ControlSurfaceDeflection && 63 this->getStructure()->ElementsOfDomain().Extent() > 0) 64 { 65 optimizeMesh(theMesher, aPS.Next()); 66 } 67 else 68 { 69 aPS.Next(); 70 } 71 } 72 73 //! Checks deviation of a mesh from geometrical surface. 74 //! Inserts additional nodes in case of huge deviation. optimizeMesh(BRepMesh_Delaun & theMesher,const Message_ProgressRange & theRange)75 virtual void optimizeMesh (BRepMesh_Delaun& theMesher, 76 const Message_ProgressRange& theRange) 77 { 78 Handle(NCollection_IncAllocator) aTmpAlloc = 79 new NCollection_IncAllocator(IMeshData::MEMORY_BLOCK_SIZE_HUGE); 80 81 mySqMinSize = this->getParameters().MinSize * this->getParameters().MinSize; 82 myCouplesMap = new IMeshData::MapOfOrientedEdges(3 * this->getStructure()->ElementsOfDomain().Extent(), aTmpAlloc); 83 myControlNodes = new IMeshData::ListOfPnt2d(aTmpAlloc); 84 myCircles = &theMesher.Circles(); 85 86 const Standard_Integer aIterationsNb = 11; 87 Standard_Boolean isInserted = Standard_True; 88 Message_ProgressScope aPS(theRange, "Iteration", aIterationsNb); 89 for (Standard_Integer aPass = 1; aPass <= aIterationsNb && isInserted && !myIsAllDegenerated; ++aPass) 90 { 91 if (!aPS.More()) 92 { 93 return; 94 } 95 // Reset stop condition 96 myMaxSqDeflection = -1.; 97 myIsAllDegenerated = Standard_True; 98 myControlNodes->Clear(); 99 100 if (this->getStructure()->ElementsOfDomain().Extent() < 1) 101 { 102 break; 103 } 104 // Iterate on current triangles 105 IMeshData::IteratorOfMapOfInteger aTriangleIt(this->getStructure()->ElementsOfDomain()); 106 for (; aTriangleIt.More(); aTriangleIt.Next()) 107 { 108 const BRepMesh_Triangle& aTriangle = this->getStructure()->GetElement(aTriangleIt.Key()); 109 splitTriangleGeometry(aTriangle); 110 } 111 112 isInserted = this->insertNodes(myControlNodes, theMesher, aPS.Next()); 113 } 114 115 myCouplesMap.Nullify(); 116 myControlNodes.Nullify(); 117 118 if (!(myMaxSqDeflection < 0.)) 119 { 120 this->getDFace()->SetDeflection(Sqrt(myMaxSqDeflection)); 121 } 122 } 123 124 private: 125 //! Contains geometrical data related to node of triangle. 126 struct TriangleNodeInfo 127 { TriangleNodeInfoBRepMesh_DelaunayDeflectionControlMeshAlgo::TriangleNodeInfo128 TriangleNodeInfo() 129 : isFrontierLink(Standard_False) 130 { 131 } 132 133 gp_XY Point2d; 134 gp_XYZ Point; 135 Standard_Boolean isFrontierLink; 136 }; 137 138 //! Functor computing deflection of a point from surface. 139 class NormalDeviation 140 { 141 public: NormalDeviation(const gp_Pnt & theRefPnt,const gp_Vec & theNormal)142 NormalDeviation( 143 const gp_Pnt& theRefPnt, 144 const gp_Vec& theNormal) 145 : myRefPnt(theRefPnt), 146 myNormal(theNormal) 147 { 148 } 149 SquareDeviation(const gp_Pnt & thePoint) const150 Standard_Real SquareDeviation(const gp_Pnt& thePoint) const 151 { 152 const Standard_Real aDeflection = Abs(myNormal.Dot(gp_Vec(myRefPnt, thePoint))); 153 return aDeflection * aDeflection; 154 } 155 156 private: 157 158 NormalDeviation (const NormalDeviation& theOther); 159 160 void operator= (const NormalDeviation& theOther); 161 162 private: 163 164 const gp_Pnt& myRefPnt; 165 const gp_Vec& myNormal; 166 }; 167 168 //! Functor computing deflection of a point on triangle link from surface. 169 class LineDeviation 170 { 171 public: 172 LineDeviation(const gp_Pnt & thePnt1,const gp_Pnt & thePnt2)173 LineDeviation( 174 const gp_Pnt& thePnt1, 175 const gp_Pnt& thePnt2) 176 : myPnt1(thePnt1), 177 myPnt2(thePnt2) 178 { 179 } 180 SquareDeviation(const gp_Pnt & thePoint) const181 Standard_Real SquareDeviation(const gp_Pnt& thePoint) const 182 { 183 return BRepMesh_GeomTool::SquareDeflectionOfSegment(myPnt1, myPnt2, thePoint); 184 } 185 186 private: 187 188 LineDeviation (const LineDeviation& theOther); 189 190 void operator= (const LineDeviation& theOther); 191 192 private: 193 const gp_Pnt& myPnt1; 194 const gp_Pnt& myPnt2; 195 }; 196 197 //! Returns nodes info of the given triangle. getTriangleInfo(const BRepMesh_Triangle & theTriangle,const Standard_Integer (& theNodesIndices)[3],TriangleNodeInfo (& theInfo)[3]) const198 void getTriangleInfo( 199 const BRepMesh_Triangle& theTriangle, 200 const Standard_Integer (&theNodesIndices)[3], 201 TriangleNodeInfo (&theInfo)[3]) const 202 { 203 const Standard_Integer(&e)[3] = theTriangle.myEdges; 204 for (Standard_Integer i = 0; i < 3; ++i) 205 { 206 const BRepMesh_Vertex& aVertex = this->getStructure()->GetNode(theNodesIndices[i]); 207 theInfo[i].Point2d = this->getRangeSplitter().Scale(aVertex.Coord(), Standard_False).XY(); 208 theInfo[i].Point = this->getNodesMap()->Value(aVertex.Location3d()).XYZ(); 209 theInfo[i].isFrontierLink = (this->getStructure()->GetLink(e[i]).Movability() == BRepMesh_Frontier); 210 } 211 } 212 213 // Check geometry of the given triangle. If triangle does not suit specified deflection, inserts new point. splitTriangleGeometry(const BRepMesh_Triangle & theTriangle)214 void splitTriangleGeometry(const BRepMesh_Triangle& theTriangle) 215 { 216 if (theTriangle.Movability() != BRepMesh_Deleted) 217 { 218 Standard_Integer aNodexIndices[3]; 219 this->getStructure()->ElementNodes(theTriangle, aNodexIndices); 220 221 TriangleNodeInfo aNodesInfo[3]; 222 getTriangleInfo(theTriangle, aNodexIndices, aNodesInfo); 223 224 gp_Vec aNormal; 225 gp_Vec aLinkVec[3]; 226 if (computeTriangleGeometry(aNodesInfo, aLinkVec, aNormal)) 227 { 228 myIsAllDegenerated = Standard_False; 229 230 const gp_XY aCenter2d = (aNodesInfo[0].Point2d + 231 aNodesInfo[1].Point2d + 232 aNodesInfo[2].Point2d) / 3.; 233 234 usePoint(aCenter2d, NormalDeviation(aNodesInfo[0].Point, aNormal)); 235 splitLinks(aNodesInfo, aNodexIndices); 236 } 237 } 238 } 239 240 //! Updates array of links vectors. 241 //! @return False on degenerative triangle. computeTriangleGeometry(const TriangleNodeInfo (& theNodesInfo)[3],gp_Vec (& theLinks)[3],gp_Vec & theNormal)242 Standard_Boolean computeTriangleGeometry( 243 const TriangleNodeInfo(&theNodesInfo)[3], 244 gp_Vec (&theLinks)[3], 245 gp_Vec &theNormal) 246 { 247 if (checkTriangleForDegenerativityAndGetLinks(theNodesInfo, theLinks)) 248 { 249 if (checkTriangleArea2d(theNodesInfo)) 250 { 251 if (computeNormal(theLinks[0], theLinks[1], theNormal)) 252 { 253 return Standard_True; 254 } 255 } 256 } 257 258 return Standard_False; 259 } 260 261 //! Updates array of links vectors. 262 //! @return False on degenerative triangle. checkTriangleForDegenerativityAndGetLinks(const TriangleNodeInfo (& theNodesInfo)[3],gp_Vec (& theLinks)[3])263 Standard_Boolean checkTriangleForDegenerativityAndGetLinks( 264 const TriangleNodeInfo (&theNodesInfo)[3], 265 gp_Vec (&theLinks)[3]) 266 { 267 const Standard_Real MinimalSqLength3d = 1.e-12; 268 for (Standard_Integer i = 0; i < 3; ++i) 269 { 270 theLinks[i] = theNodesInfo[(i + 1) % 3].Point - theNodesInfo[i].Point; 271 if (theLinks[i].SquareMagnitude() < MinimalSqLength3d) 272 { 273 return Standard_False; 274 } 275 } 276 277 return Standard_True; 278 } 279 280 //! Checks area of triangle in parametric space for degenerativity. 281 //! @return False on degenerative triangle. checkTriangleArea2d(const TriangleNodeInfo (& theNodesInfo)[3])282 Standard_Boolean checkTriangleArea2d( 283 const TriangleNodeInfo (&theNodesInfo)[3]) 284 { 285 const gp_Vec2d aLink2d1(theNodesInfo[0].Point2d, theNodesInfo[1].Point2d); 286 const gp_Vec2d aLink2d2(theNodesInfo[1].Point2d, theNodesInfo[2].Point2d); 287 288 const Standard_Real MinimalArea2d = 1.e-9; 289 return (Abs(aLink2d1 ^ aLink2d2) > MinimalArea2d); 290 } 291 292 //! Computes normal using two link vectors. 293 //! @return True on success, False in case of normal of null magnitude. computeNormal(const gp_Vec & theLink1,const gp_Vec & theLink2,gp_Vec & theNormal)294 Standard_Boolean computeNormal(const gp_Vec& theLink1, 295 const gp_Vec& theLink2, 296 gp_Vec& theNormal) 297 { 298 const gp_Vec aNormal(theLink1 ^ theLink2); 299 if (aNormal.SquareMagnitude() > gp::Resolution()) 300 { 301 theNormal = aNormal.Normalized(); 302 return Standard_True; 303 } 304 305 return Standard_False; 306 } 307 308 //! Computes deflection of midpoints of triangles links. 309 //! @return True if point fits specified deflection. splitLinks(const TriangleNodeInfo (& theNodesInfo)[3],const Standard_Integer (& theNodesIndices)[3])310 void splitLinks( 311 const TriangleNodeInfo (&theNodesInfo)[3], 312 const Standard_Integer (&theNodesIndices)[3]) 313 { 314 // Check deflection at triangle links 315 for (Standard_Integer i = 0; i < 3; ++i) 316 { 317 if (theNodesInfo[i].isFrontierLink) 318 { 319 continue; 320 } 321 322 const Standard_Integer j = (i + 1) % 3; 323 // Check if this link was already processed 324 Standard_Integer aFirstVertex, aLastVertex; 325 if (theNodesIndices[i] < theNodesIndices[j]) 326 { 327 aFirstVertex = theNodesIndices[i]; 328 aLastVertex = theNodesIndices[j]; 329 } 330 else 331 { 332 aFirstVertex = theNodesIndices[j]; 333 aLastVertex = theNodesIndices[i]; 334 } 335 336 if (myCouplesMap->Add(BRepMesh_OrientedEdge(aFirstVertex, aLastVertex))) 337 { 338 const gp_XY aMidPnt2d = (theNodesInfo[i].Point2d + 339 theNodesInfo[j].Point2d) / 2.; 340 341 if (!usePoint (aMidPnt2d, LineDeviation (theNodesInfo[i].Point, 342 theNodesInfo[j].Point))) 343 { 344 if (!rejectSplitLinksForMinSize (theNodesInfo[i], 345 theNodesInfo[j], 346 aMidPnt2d)) 347 { 348 if (!checkLinkEndsForAngularDeviation (theNodesInfo[i], 349 theNodesInfo[j], 350 aMidPnt2d)) 351 { 352 myControlNodes->Append(aMidPnt2d); 353 } 354 } 355 } 356 } 357 } 358 } 359 360 //! Checks that two links produced as the result of a split of 361 //! the given link by the middle point fit MinSize requirement. rejectSplitLinksForMinSize(const TriangleNodeInfo & theNodeInfo1,const TriangleNodeInfo & theNodeInfo2,const gp_XY & theMidPoint)362 Standard_Boolean rejectSplitLinksForMinSize (const TriangleNodeInfo& theNodeInfo1, 363 const TriangleNodeInfo& theNodeInfo2, 364 const gp_XY& theMidPoint) 365 { 366 const gp_Pnt aPnt = getPoint3d (theMidPoint); 367 return ((theNodeInfo1.Point - aPnt.XYZ()).SquareModulus() < mySqMinSize || 368 (theNodeInfo2.Point - aPnt.XYZ()).SquareModulus() < mySqMinSize); 369 } 370 371 //! Checks the given point (located between the given nodes) 372 //! for specified angular deviation. checkLinkEndsForAngularDeviation(const TriangleNodeInfo & theNodeInfo1,const TriangleNodeInfo & theNodeInfo2,const gp_XY &)373 Standard_Boolean checkLinkEndsForAngularDeviation(const TriangleNodeInfo& theNodeInfo1, 374 const TriangleNodeInfo& theNodeInfo2, 375 const gp_XY& /*theMidPoint*/) 376 { 377 gp_Dir aNorm1, aNorm2; 378 const Handle(Geom_Surface)& aSurf = this->getDFace()->GetSurface()->Surface().Surface(); 379 380 if ((GeomLib::NormEstim(aSurf, theNodeInfo1.Point2d, Precision::Confusion(), aNorm1) == 0) && 381 (GeomLib::NormEstim(aSurf, theNodeInfo2.Point2d, Precision::Confusion(), aNorm2) == 0)) 382 { 383 Standard_Real anAngle = aNorm1.Angle(aNorm2); 384 if (anAngle > this->getParameters().AngleInterior) 385 { 386 return Standard_False; 387 } 388 } 389 #if 0 390 else if (GeomLib::NormEstim(aSurf, theMidPoint, Precision::Confusion(), aNorm1) != 0) 391 { 392 // It is better to consider the singular point as a node of triangulation. 393 // However, it leads to hangs up meshing some faces (including faces with 394 // degenerated edges). E.g. tests "mesh standard_incmesh Q6". 395 // So, this code fragment is better to implement in the future. 396 return Standard_False; 397 } 398 #endif 399 400 return Standard_True; 401 } 402 403 //! Returns 3d point corresponding to the given one in 2d space. getPoint3d(const gp_XY & thePnt2d)404 gp_Pnt getPoint3d (const gp_XY& thePnt2d) 405 { 406 gp_Pnt aPnt; 407 this->getDFace()->GetSurface()->D0(thePnt2d.X(), thePnt2d.Y(), aPnt); 408 return aPnt; 409 } 410 411 //! Computes deflection of the given point and caches it for 412 //! insertion in case if it overflows deflection. 413 //! @return True if point has been cached for insertion. 414 template<class DeflectionFunctor> usePoint(const gp_XY & thePnt2d,const DeflectionFunctor & theDeflectionFunctor)415 Standard_Boolean usePoint( 416 const gp_XY& thePnt2d, 417 const DeflectionFunctor& theDeflectionFunctor) 418 { 419 const gp_Pnt aPnt = getPoint3d (thePnt2d); 420 if (!checkDeflectionOfPointAndUpdateCache(thePnt2d, aPnt, theDeflectionFunctor.SquareDeviation(aPnt))) 421 { 422 myControlNodes->Append(thePnt2d); 423 return Standard_True; 424 } 425 426 return Standard_False; 427 } 428 429 //! Checks the given point for specified linear deflection. 430 //! Updates value of total mesh defleciton. checkDeflectionOfPointAndUpdateCache(const gp_XY & thePnt2d,const gp_Pnt & thePnt3d,const Standard_Real theSqDeflection)431 Standard_Boolean checkDeflectionOfPointAndUpdateCache( 432 const gp_XY& thePnt2d, 433 const gp_Pnt& thePnt3d, 434 const Standard_Real theSqDeflection) 435 { 436 if (theSqDeflection > myMaxSqDeflection) 437 { 438 myMaxSqDeflection = theSqDeflection; 439 } 440 441 const Standard_Real aSqDeflection = 442 this->getDFace()->GetDeflection() * this->getDFace()->GetDeflection(); 443 if (theSqDeflection < aSqDeflection) 444 { 445 return Standard_True; 446 } 447 448 return rejectByMinSize(thePnt2d, thePnt3d); 449 } 450 451 //! Checks distance between the given node and nodes of triangles 452 //! shot by it for MinSize criteria. 453 //! This check is expected to roughly estimate and prevent 454 //! generation of triangles with sides smaller than MinSize. rejectByMinSize(const gp_XY & thePnt2d,const gp_Pnt & thePnt3d)455 Standard_Boolean rejectByMinSize( 456 const gp_XY& thePnt2d, 457 const gp_Pnt& thePnt3d) 458 { 459 IMeshData::MapOfInteger aUsedNodes; 460 IMeshData::ListOfInteger& aCirclesList = 461 const_cast<BRepMesh_CircleTool&>(*myCircles).Select( 462 this->getRangeSplitter().Scale(thePnt2d, Standard_True).XY()); 463 464 IMeshData::ListOfInteger::Iterator aCircleIt(aCirclesList); 465 for (; aCircleIt.More(); aCircleIt.Next()) 466 { 467 const BRepMesh_Triangle& aTriangle = this->getStructure()->GetElement(aCircleIt.Value()); 468 469 Standard_Integer aNodes[3]; 470 this->getStructure()->ElementNodes(aTriangle, aNodes); 471 472 for (Standard_Integer i = 0; i < 3; ++i) 473 { 474 if (!aUsedNodes.Contains(aNodes[i])) 475 { 476 aUsedNodes.Add(aNodes[i]); 477 const BRepMesh_Vertex& aVertex = this->getStructure()->GetNode(aNodes[i]); 478 const gp_Pnt& aPoint = this->getNodesMap()->Value(aVertex.Location3d()); 479 480 if (thePnt3d.SquareDistance(aPoint) < mySqMinSize) 481 { 482 return Standard_True; 483 } 484 } 485 } 486 } 487 488 return Standard_False; 489 } 490 491 private: 492 Standard_Real myMaxSqDeflection; 493 Standard_Real mySqMinSize; 494 Standard_Boolean myIsAllDegenerated; 495 Handle(IMeshData::MapOfOrientedEdges) myCouplesMap; 496 Handle(IMeshData::ListOfPnt2d) myControlNodes; 497 const BRepMesh_CircleTool* myCircles; 498 }; 499 500 #endif 501