0031642: Visualization - crash in Graphic3d_Structure::SetVisual() on redisplaying...
[occt.git] / src / MeshTest / MeshTest_CheckTopology.cxx
CommitLineData
b311480e 1// Created on: 2004-05-10
2// Created by: Michael SAZONOV
973c2be1 3// Copyright (c) 2004-2014 OPEN CASCADE SAS
b311480e 4//
973c2be1 5// This file is part of Open CASCADE Technology software library.
b311480e 6//
d5f74e42 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
973c2be1 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.
b311480e 12//
973c2be1 13// Alternatively, this file may be used under the terms of Open CASCADE
14// commercial license or contractual agreement.
7fd59977 15
16#include <MeshTest_CheckTopology.hxx>
17#include <BRep_Tool.hxx>
446e11f3 18#include <TColStd_PackedMapOfInteger.hxx>
7fd59977 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>
7bd071ed 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//=======================================================================
38static 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//=======================================================================
49static 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}
7fd59977 55
56//=======================================================================
57//function : Perform
58//purpose : Performs checking
59//=======================================================================
60
d2c43192 61void MeshTest_CheckTopology::Perform (Draw_Interpretor& di)
7fd59977 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()) {
0797d9d3 82#ifdef OCCT_DEBUG
04232180 83 std::cout<<"problem getting PolygonOnTriangulation of edge "<<ie<<std::endl;
7fd59977 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()) {
0797d9d3 99#ifdef OCCT_DEBUG
04232180 100 std::cout<<"problem getting PolygonOnTriangulation of edge "<<ie<<std::endl;
7fd59977 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
061cd2d8 113 Standard_Real aSqDefle = BRep_Tool::Tolerance(aEdge);
7bd071ed 114 aSqDefle *= aSqDefle;
7fd59977 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();
7bd071ed 121 const gp_Trsf &aTrsf1 = aFace1.Location().Transformation();
122 const gp_Trsf &aTrsf2 = aFace2.Location().Transformation();
7fd59977 123 for (; i1 <= aNodes1.Upper(); i1++, i2++) {
7bd071ed 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 {
7fd59977 129 myErrors.Append(iF1);
130 myErrors.Append(i1);
131 myErrors.Append(iF2);
132 myErrors.Append(i2);
7bd071ed 133 myErrorsVal.Append(Sqrt(aSqDist));
7fd59977 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()) {
586db386 146 di << "face " <<iF <<" has no triangulation\n";
7fd59977 147 continue;
148 }
446e11f3 149
7bd071ed 150 const gp_Trsf &aTrsf = aLoc.Transformation();
151
7fd59977 152 // remember boundary nodes
446e11f3 153 TColStd_PackedMapOfInteger aMapBndNodes;
7fd59977 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
446e11f3
A
166 TColStd_PackedMapOfInteger aUsedNodes;
167
168 // check of free links and nodes
7fd59977 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]);
446e11f3
A
174
175 aUsedNodes.Add (n[0]);
176 aUsedNodes.Add (n[1]);
177 aUsedNodes.Add (n[2]);
178
7bd071ed 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
7fd59977 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)) {
e3a6386d 213 Handle(TColStd_HSequenceOfInteger) tmpSeq = new TColStd_HSequenceOfInteger;
7fd59977 214 myMapFaceLinks.Add(iF, tmpSeq);
215 }
e3a6386d 216 Handle(TColStd_HSequenceOfInteger)& aSeq = myMapFaceLinks.ChangeFromKey(iF);
217 aSeq->Append(n1);
218 aSeq->Append(n2);
7fd59977 219 }
220 }
221 }
51740958 222
446e11f3
A
223 // check of free nodes
224 Standard_Integer aNbNodes = aT->NbNodes();
51740958 225 for (Standard_Integer k = 1; k <= aNbNodes; k++)
226 if ( ! aUsedNodes.Contains(k) )
446e11f3
A
227 {
228 myFreeNodeFaces.Append (iF);
51740958 229 myFreeNodeNums.Append (k);
446e11f3 230 }
7fd59977 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
240void MeshTest_CheckTopology::GetFreeLink(const Standard_Integer theFaceIndex,
241 const Standard_Integer theLinkIndex,
242 Standard_Integer& theNode1,
243 Standard_Integer& theNode2) const
244{
e3a6386d 245 const Handle(TColStd_HSequenceOfInteger)& aSeq = myMapFaceLinks(theFaceIndex);
7fd59977 246 Standard_Integer aInd = (theLinkIndex-1)*2 + 1;
e3a6386d 247 theNode1 = aSeq->Value(aInd);
248 theNode2 = aSeq->Value(aInd+1);
7fd59977 249}
250
251//=======================================================================
252//function : GetCrossFaceError
253//purpose : gets the attributes of a cross face error with the given index
254//=======================================================================
255
256void 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}