0026106: BRepMesh - revision of data model
[occt.git] / src / MeshTest / MeshTest_CheckTopology.cxx
1 // Created on: 2004-05-10
2 // Created by: Michael SAZONOV
3 // Copyright (c) 2004-2014 OPEN CASCADE SAS
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 #include <MeshTest_CheckTopology.hxx>
17 #include <BRep_Tool.hxx>
18 #include <TColStd_PackedMapOfInteger.hxx>
19 #include <TopExp.hxx>
20 #include <TopExp_Explorer.hxx>
21 #include <TopoDS_Edge.hxx>
22 #include <TopoDS.hxx>
23 #include <TopoDS_Face.hxx>
24 #include <TopLoc_Location.hxx>
25 #include <TopTools_IndexedMapOfShape.hxx>
26 #include <TopTools_IndexedDataMapOfShapeListOfShape.hxx>
27 #include <TopTools_ListOfShape.hxx>
28 #include <TopTools_ListIteratorOfListOfShape.hxx>
29 #include <Poly_Triangulation.hxx>
30 #include <Poly_PolygonOnTriangulation.hxx>
31 #include <Poly_Connect.hxx>
32 #include <Precision.hxx>
33
34 //=======================================================================
35 //function : ComputeArea
36 //purpose  : Computes area of the triangle given by its three points (either 2D or3D)
37 //=======================================================================
38 static Standard_Real ComputeArea(const gp_XYZ& theP1,
39                                  const gp_XYZ& theP2,
40                                  const gp_XYZ& theP3)
41 {
42   return 0.5*(theP3 - theP1).Crossed(theP2 - theP1).Modulus();
43 }
44
45 //=======================================================================
46 //function : ComputeArea
47 //purpose  : Computes area of the triangle given by its three points (either 2D or3D)
48 //=======================================================================
49 static Standard_Real ComputeArea(const gp_XY& theP1,
50                                  const gp_XY& theP2,
51                                  const gp_XY& theP3)
52 {
53   return 0.5*Abs((theP3 - theP1).Crossed(theP2 - theP1));
54 }
55
56 //=======================================================================
57 //function : Perform
58 //purpose  : Performs checking
59 //=======================================================================
60
61 void MeshTest_CheckTopology::Perform (Draw_Interpretor& di)
62 {
63   TopTools_IndexedMapOfShape aMapF;
64   TopTools_IndexedDataMapOfShapeListOfShape aMapEF;
65   TopExp::MapShapes (myShape, TopAbs_FACE, aMapF);
66   TopExp::MapShapesAndAncestors (myShape, TopAbs_EDGE, TopAbs_FACE, aMapEF);
67
68   // check polygons
69   Standard_Integer ie;
70   for (ie=1; ie <= aMapEF.Extent(); ie++) {
71     const TopoDS_Edge& aEdge = TopoDS::Edge(aMapEF.FindKey(ie));
72     const TopTools_ListOfShape& aFaces = aMapEF(ie);
73     if (aFaces.Extent() < 2) continue;
74
75     // get polygon on first face
76     const TopoDS_Face& aFace1 = TopoDS::Face(aFaces.First());
77     TopLoc_Location aLoc1;
78     Handle(Poly_Triangulation) aT1 = BRep_Tool::Triangulation(aFace1, aLoc1);
79     Handle(Poly_PolygonOnTriangulation) aPoly1 =
80       BRep_Tool::PolygonOnTriangulation(aEdge, aT1, aLoc1);
81     if (aPoly1.IsNull() || aT1.IsNull()) {
82 #ifdef OCCT_DEBUG
83       cout<<"problem getting PolygonOnTriangulation of edge "<<ie<<endl;
84 #endif
85       continue;
86     }
87     const TColStd_Array1OfInteger& aNodes1 = aPoly1->Nodes();
88
89     // cycle on other polygons
90     TopTools_ListIteratorOfListOfShape it(aFaces);
91     it.Next();
92     for (; it.More(); it.Next()) {
93       const TopoDS_Face& aFace2 = TopoDS::Face(it.Value());
94       TopLoc_Location aLoc2;
95       Handle(Poly_Triangulation) aT2 = BRep_Tool::Triangulation(aFace2, aLoc2);
96       Handle(Poly_PolygonOnTriangulation) aPoly2 =
97         BRep_Tool::PolygonOnTriangulation(aEdge, aT2, aLoc2);
98       if (aPoly2.IsNull() || aT2.IsNull()) {
99 #ifdef OCCT_DEBUG
100         cout<<"problem getting PolygonOnTriangulation of edge "<<ie<<endl;
101 #endif
102         continue;
103       }
104       const TColStd_Array1OfInteger& aNodes2 = aPoly2->Nodes();
105
106       // check equality of polygons lengths
107       if (aNodes2.Length() != aNodes1.Length()) {
108         myAsyncEdges.Append(ie);
109         break;
110       }
111
112       // check distances between corresponding points
113       Standard_Real aSqDefle = Max(aT1->Deflection(), aT2->Deflection());
114       aSqDefle *= aSqDefle;
115       const TColgp_Array1OfPnt& aPoints1 = aT1->Nodes();
116       const TColgp_Array1OfPnt& aPoints2 = aT2->Nodes();
117       Standard_Integer iF1 = aMapF.FindIndex(aFace1);
118       Standard_Integer iF2 = aMapF.FindIndex(aFace2);
119       Standard_Integer i1 = aNodes1.Lower();
120       Standard_Integer i2 = aNodes2.Lower();
121       const gp_Trsf &aTrsf1 = aFace1.Location().Transformation();
122       const gp_Trsf &aTrsf2 = aFace2.Location().Transformation();
123       for (; i1 <= aNodes1.Upper(); i1++, i2++) {
124         const gp_Pnt aP1 = aPoints1(aNodes1(i1)).Transformed(aTrsf1);
125         const gp_Pnt aP2 = aPoints2(aNodes2(i2)).Transformed(aTrsf2);
126         const Standard_Real aSqDist = aP1.SquareDistance(aP2);
127         if (aSqDist > aSqDefle)
128         {
129           myErrors.Append(iF1);
130           myErrors.Append(i1);
131           myErrors.Append(iF2);
132           myErrors.Append(i2);
133           myErrorsVal.Append(Sqrt(aSqDist));
134         }
135       }
136     }
137   }
138
139   // check triangulations
140   Standard_Integer iF;
141   for (iF=1; iF <= aMapF.Extent(); iF++) {
142     const TopoDS_Face& aFace = TopoDS::Face(aMapF.FindKey(iF));
143     TopLoc_Location aLoc;
144     Handle(Poly_Triangulation) aT = BRep_Tool::Triangulation(aFace, aLoc);
145     if (aT.IsNull()) {
146       di << "face " <<iF <<" has no triangulation\n";
147       continue;
148     }
149
150     const gp_Trsf &aTrsf = aLoc.Transformation();
151
152     // remember boundary nodes
153     TColStd_PackedMapOfInteger aMapBndNodes;
154     TopExp_Explorer ex(aFace, TopAbs_EDGE);
155     for (; ex.More(); ex.Next()) {
156       const TopoDS_Edge& aEdge = TopoDS::Edge(ex.Current());
157       Handle(Poly_PolygonOnTriangulation) aPoly =
158         BRep_Tool::PolygonOnTriangulation(aEdge, aT, aLoc);
159       if (aPoly.IsNull()) continue;
160       const TColStd_Array1OfInteger& aNodes = aPoly->Nodes();
161       Standard_Integer i;
162       for (i=aNodes.Lower(); i <= aNodes.Upper(); i++)
163         aMapBndNodes.Add(aNodes(i));
164     }
165
166     TColStd_PackedMapOfInteger aUsedNodes;
167
168     // check of free links and nodes
169     Poly_Connect aConn(aT);
170     const Poly_Array1OfTriangle& aTriangles = aT->Triangles();
171     Standard_Integer nbTri = aT->NbTriangles(), i, j, n[3], t[3];
172     for (i = 1; i <= nbTri; i++) {
173       aTriangles(i).Get(n[0], n[1], n[2]);
174       
175       aUsedNodes.Add (n[0]);
176       aUsedNodes.Add (n[1]);
177       aUsedNodes.Add (n[2]);
178
179       const gp_Pnt aPts[3] = {aT->Node(n[0]).Transformed(aTrsf),
180                               aT->Node(n[1]).Transformed(aTrsf),
181                               aT->Node(n[2]).Transformed(aTrsf)};
182
183       Standard_Real anArea = ComputeArea(aPts[0].XYZ(), aPts[1].XYZ(), aPts[2].XYZ());
184       if (anArea < Precision::SquareConfusion())
185       {
186         mySmallTrianglesFaces.Append(iF);
187         mySmallTrianglesTriangles.Append(i);
188       }
189       else if (aT->HasUVNodes())
190       {
191         const gp_XY aPUV[3] = {aT->UVNode(n[0]).XY(),
192                                aT->UVNode(n[1]).XY(),
193                                aT->UVNode(n[2]).XY()};
194         anArea = ComputeArea(aPUV[0], aPUV[1], aPUV[2]);
195         if (anArea < Precision::SquarePConfusion())
196         {
197           mySmallTrianglesFaces.Append(iF);
198           mySmallTrianglesTriangles.Append(i);
199         }
200       }
201
202       aConn.Triangles(i, t[0], t[1], t[2]);
203       for (j = 0; j < 3; j++) {
204         if (t[j] == 0) {
205           // free link found
206           Standard_Integer k = (j+1) % 3;  // the following node of the edge
207           Standard_Integer n1 = n[j];
208           Standard_Integer n2 = n[k];
209           // skip if it is on boundary
210           if (aMapBndNodes.Contains(n1) && aMapBndNodes.Contains(n2))
211             continue;
212           if (!myMapFaceLinks.Contains(iF)) {
213             Handle(TColStd_HSequenceOfInteger) tmpSeq = new TColStd_HSequenceOfInteger;
214             myMapFaceLinks.Add(iF, tmpSeq);
215           }
216           Handle(TColStd_HSequenceOfInteger)& aSeq = myMapFaceLinks.ChangeFromKey(iF);
217           aSeq->Append(n1);
218           aSeq->Append(n2);
219         }
220       }
221     }
222     
223     // check of free nodes
224     Standard_Integer aNbNodes = aT->NbNodes();
225     for (Standard_Integer k = 1; k <= aNbNodes; k++)
226       if ( ! aUsedNodes.Contains(k) )
227       {
228         myFreeNodeFaces.Append (iF);
229         myFreeNodeNums.Append (k);
230       }
231   }
232 }
233
234 //=======================================================================
235 //function : GetFreeLink
236 //purpose  : gets the numbers of nodes of a free link with the given index
237 //           in the face with the given index
238 //=======================================================================
239
240 void MeshTest_CheckTopology::GetFreeLink(const Standard_Integer theFaceIndex,
241                                          const Standard_Integer theLinkIndex,
242                                          Standard_Integer& theNode1,
243                                          Standard_Integer& theNode2) const
244 {
245   const Handle(TColStd_HSequenceOfInteger)& aSeq = myMapFaceLinks(theFaceIndex);
246   Standard_Integer aInd = (theLinkIndex-1)*2 + 1;
247   theNode1 = aSeq->Value(aInd);
248   theNode2 = aSeq->Value(aInd+1);
249 }
250
251 //=======================================================================
252 //function : GetCrossFaceError
253 //purpose  : gets the attributes of a cross face error with the given index
254 //=======================================================================
255
256 void MeshTest_CheckTopology::GetCrossFaceError(const Standard_Integer theIndex,
257                                                Standard_Integer& theFace1,
258                                                Standard_Integer& theNode1,
259                                                Standard_Integer& theFace2,
260                                                Standard_Integer& theNode2,
261                                                Standard_Real&    theValue) const
262 {
263   Standard_Integer aInd = (theIndex-1)*4 + 1;
264   theFace1 = myErrors(aInd);
265   theNode1 = myErrors(aInd+1);
266   theFace2 = myErrors(aInd+2);
267   theNode2 = myErrors(aInd+3);
268   theValue = myErrorsVal(theIndex);
269 }