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       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 }