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