| 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 | std::cout<<"problem getting PolygonOnTriangulation of edge "<<ie<<std::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 | std::cout<<"problem getting PolygonOnTriangulation of edge "<<ie<<std::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 = BRep_Tool::Tolerance(aEdge); |
| 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 | } |