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