[occt.git] / src / MeshTest / MeshTest_CheckTopology.cxx

// File:      MeshTest_CheckTopology.cxx
// Created:   5.10.2004
// Author:    Michael SAZONOV

#include <MeshTest_CheckTopology.hxx>
#include <BRep_Tool.hxx>
#include <TColStd_PackedMapOfInteger.hxx>
#include <TopExp.hxx>
#include <TopExp_Explorer.hxx>
#include <TopoDS_Edge.hxx>
#include <TopoDS.hxx>
#include <TopoDS_Face.hxx>
#include <TopLoc_Location.hxx>
#include <TopTools_IndexedMapOfShape.hxx>
#include <TopTools_IndexedDataMapOfShapeListOfShape.hxx>
#include <TopTools_ListOfShape.hxx>
#include <TopTools_ListIteratorOfListOfShape.hxx>
#include <Poly_Triangulation.hxx>
#include <Poly_PolygonOnTriangulation.hxx>
#include <Poly_Connect.hxx>

//=======================================================================
//function : Perform
//purpose  : Performs checking
//=======================================================================

void MeshTest_CheckTopology::Perform ()
{
  TopTools_IndexedMapOfShape aMapF;
  TopTools_IndexedDataMapOfShapeListOfShape aMapEF;
  TopExp::MapShapes (myShape, TopAbs_FACE, aMapF);
  TopExp::MapShapesAndAncestors (myShape, TopAbs_EDGE, TopAbs_FACE, aMapEF);

  // check polygons
  Standard_Integer ie;
  for (ie=1; ie <= aMapEF.Extent(); ie++) {
    const TopoDS_Edge& aEdge = TopoDS::Edge(aMapEF.FindKey(ie));
    const TopTools_ListOfShape& aFaces = aMapEF(ie);
    if (aFaces.Extent() < 2) continue;

    // get polygon on first face
    const TopoDS_Face& aFace1 = TopoDS::Face(aFaces.First());
    TopLoc_Location aLoc1;
    Handle(Poly_Triangulation) aT1 = BRep_Tool::Triangulation(aFace1, aLoc1);
    Handle(Poly_PolygonOnTriangulation) aPoly1 =
      BRep_Tool::PolygonOnTriangulation(aEdge, aT1, aLoc1);
    if (aPoly1.IsNull() || aT1.IsNull()) {
#ifdef DEB
      cout<<"problem getting PolygonOnTriangulation of edge "<<ie<<endl;
#endif
      continue;
    }
    const TColStd_Array1OfInteger& aNodes1 = aPoly1->Nodes();

    // cycle on other polygons
    TopTools_ListIteratorOfListOfShape it(aFaces);
    it.Next();
    for (; it.More(); it.Next()) {
      const TopoDS_Face& aFace2 = TopoDS::Face(it.Value());
      TopLoc_Location aLoc2;
      Handle(Poly_Triangulation) aT2 = BRep_Tool::Triangulation(aFace2, aLoc2);
      Handle(Poly_PolygonOnTriangulation) aPoly2 =
	BRep_Tool::PolygonOnTriangulation(aEdge, aT2, aLoc2);
      if (aPoly2.IsNull() || aT2.IsNull()) {
#ifdef DEB
	cout<<"problem getting PolygonOnTriangulation of edge "<<ie<<endl;
#endif
	continue;
      }
      const TColStd_Array1OfInteger& aNodes2 = aPoly2->Nodes();

      // check equality of polygons lengths
      if (aNodes2.Length() != aNodes1.Length()) {
	myAsyncEdges.Append(ie);
	break;
      }

      // check distances between corresponding points
      Standard_Real aDefle = Max(aT1->Deflection(), aT2->Deflection());
      const TColgp_Array1OfPnt& aPoints1 = aT1->Nodes();
      const TColgp_Array1OfPnt& aPoints2 = aT2->Nodes();
      Standard_Integer iF1 = aMapF.FindIndex(aFace1);
      Standard_Integer iF2 = aMapF.FindIndex(aFace2);
      Standard_Integer i1 = aNodes1.Lower();
      Standard_Integer i2 = aNodes2.Lower();
      gp_Trsf aTrsf1 = aFace1.Location().Transformation();
      gp_Trsf aTrsf2 = aFace2.Location().Transformation();
      for (; i1 <= aNodes1.Upper(); i1++, i2++) {
	gp_Pnt aP1 = aPoints1(aNodes1(i1)).Transformed(aTrsf1);
	gp_Pnt aP2 = aPoints2(aNodes2(i2)).Transformed(aTrsf2);
	Standard_Real aDist = aP1.Distance(aP2);
	if (aDist > aDefle) {
	  myErrors.Append(iF1);
	  myErrors.Append(i1);
	  myErrors.Append(iF2);
	  myErrors.Append(i2);
	  myErrorsVal.Append(aDist);
	}
      }
    }
  }

  // check triangulations
  Standard_Integer iF;
  for (iF=1; iF <= aMapF.Extent(); iF++) {
    const TopoDS_Face& aFace = TopoDS::Face(aMapF.FindKey(iF));
    TopLoc_Location aLoc;
    Handle(Poly_Triangulation) aT = BRep_Tool::Triangulation(aFace, aLoc);
    if (aT.IsNull()) {
      cout<< "face "<<iF<<" has no triangulation"<<endl;
      continue;
    }

    // remember boundary nodes
    TColStd_PackedMapOfInteger aMapBndNodes;
    TopExp_Explorer ex(aFace, TopAbs_EDGE);
    for (; ex.More(); ex.Next()) {
      const TopoDS_Edge& aEdge = TopoDS::Edge(ex.Current());
      Handle(Poly_PolygonOnTriangulation) aPoly =
	BRep_Tool::PolygonOnTriangulation(aEdge, aT, aLoc);
      if (aPoly.IsNull()) continue;
      const TColStd_Array1OfInteger& aNodes = aPoly->Nodes();
      Standard_Integer i;
      for (i=aNodes.Lower(); i <= aNodes.Upper(); i++)
	aMapBndNodes.Add(aNodes(i));
    }

    TColStd_PackedMapOfInteger aUsedNodes;

    // check of free links and nodes
    Poly_Connect aConn(aT);
    const Poly_Array1OfTriangle& aTriangles = aT->Triangles();
    Standard_Integer nbTri = aT->NbTriangles(), i, j, n[3], t[3];
    for (i = 1; i <= nbTri; i++) {
      aTriangles(i).Get(n[0], n[1], n[2]);
      
      aUsedNodes.Add (n[0]);
      aUsedNodes.Add (n[1]);
      aUsedNodes.Add (n[2]);

      aConn.Triangles(i, t[0], t[1], t[2]);
      for (j = 0; j < 3; j++) {
	if (t[j] == 0) {
	  // free link found
	  Standard_Integer k = (j+1) % 3;  // the following node of the edge
	  Standard_Integer n1 = n[j];
	  Standard_Integer n2 = n[k];
	  // skip if it is on boundary
	  if (aMapBndNodes.Contains(n1) && aMapBndNodes.Contains(n2))
	    continue;
	  if (!myMapFaceLinks.Contains(iF)) {
	    //myMapFaceLinks.Add(iF, TColStd_SequenceOfInteger());
            TColStd_SequenceOfInteger tmpSeq;
	    myMapFaceLinks.Add(iF, tmpSeq);
          }
	  TColStd_SequenceOfInteger& aSeq = myMapFaceLinks.ChangeFromKey(iF);
	  aSeq.Append(n1);
	  aSeq.Append(n2);
	}
      }
    }

    // check of free nodes
    Standard_Integer aNbNodes = aT->NbNodes();
    for (Standard_Integer i = 1; i <= aNbNodes; i++)
      if ( ! aUsedNodes.Contains(i) )
      {
	myFreeNodeFaces.Append (iF);
	myFreeNodeNums.Append (i);
      }
  }
}

//=======================================================================
//function : GetFreeLink
//purpose  : gets the numbers of nodes of a free link with the given index
//           in the face with the given index
//=======================================================================

void MeshTest_CheckTopology::GetFreeLink(const Standard_Integer theFaceIndex,
	 			 	 const Standard_Integer theLinkIndex,
					 Standard_Integer& theNode1,
					 Standard_Integer& theNode2) const
{
  const TColStd_SequenceOfInteger& aSeq = myMapFaceLinks(theFaceIndex);
  Standard_Integer aInd = (theLinkIndex-1)*2 + 1;
  theNode1 = aSeq(aInd);
  theNode2 = aSeq(aInd+1);
}

//=======================================================================
//function : GetCrossFaceError
//purpose  : gets the attributes of a cross face error with the given index
//=======================================================================

void MeshTest_CheckTopology::GetCrossFaceError(const Standard_Integer theIndex,
					       Standard_Integer& theFace1,
					       Standard_Integer& theNode1,
					       Standard_Integer& theFace2,
					       Standard_Integer& theNode2,
					       Standard_Real&    theValue) const
{
  Standard_Integer aInd = (theIndex-1)*4 + 1;
  theFace1 = myErrors(aInd);
  theNode1 = myErrors(aInd+1);
  theFace2 = myErrors(aInd+2);
  theNode2 = myErrors(aInd+3);
  theValue = myErrorsVal(theIndex);
}
Commit	Line	Data
7fd59977	1	// File: MeshTest_CheckTopology.cxx
	2	// Created: 5.10.2004
	3	// Author: Michael SAZONOV
	4
	5	#include <MeshTest_CheckTopology.hxx>
	6	#include <BRep_Tool.hxx>
446e11f3	7	#include <TColStd_PackedMapOfInteger.hxx>
7fd59977	8	#include <TopExp.hxx>
	9	#include <TopExp_Explorer.hxx>
	10	#include <TopoDS_Edge.hxx>
	11	#include <TopoDS.hxx>
	12	#include <TopoDS_Face.hxx>
	13	#include <TopLoc_Location.hxx>
	14	#include <TopTools_IndexedMapOfShape.hxx>
	15	#include <TopTools_IndexedDataMapOfShapeListOfShape.hxx>
	16	#include <TopTools_ListOfShape.hxx>
	17	#include <TopTools_ListIteratorOfListOfShape.hxx>
	18	#include <Poly_Triangulation.hxx>
	19	#include <Poly_PolygonOnTriangulation.hxx>
	20	#include <Poly_Connect.hxx>
	21
	22	//=======================================================================
	23	//function : Perform
	24	//purpose : Performs checking
	25	//=======================================================================
	26
	27	void MeshTest_CheckTopology::Perform ()
	28	{
	29	TopTools_IndexedMapOfShape aMapF;
	30	TopTools_IndexedDataMapOfShapeListOfShape aMapEF;
	31	TopExp::MapShapes (myShape, TopAbs_FACE, aMapF);
	32	TopExp::MapShapesAndAncestors (myShape, TopAbs_EDGE, TopAbs_FACE, aMapEF);
	33
	34	// check polygons
	35	Standard_Integer ie;
	36	for (ie=1; ie <= aMapEF.Extent(); ie++) {
	37	const TopoDS_Edge& aEdge = TopoDS::Edge(aMapEF.FindKey(ie));
	38	const TopTools_ListOfShape& aFaces = aMapEF(ie);
	39	if (aFaces.Extent() < 2) continue;
	40
	41	// get polygon on first face
	42	const TopoDS_Face& aFace1 = TopoDS::Face(aFaces.First());
	43	TopLoc_Location aLoc1;
	44	Handle(Poly_Triangulation) aT1 = BRep_Tool::Triangulation(aFace1, aLoc1);
	45	Handle(Poly_PolygonOnTriangulation) aPoly1 =
	46	BRep_Tool::PolygonOnTriangulation(aEdge, aT1, aLoc1);
	47	if (aPoly1.IsNull() \|\| aT1.IsNull()) {
	48	#ifdef DEB
	49	cout<<"problem getting PolygonOnTriangulation of edge "<<ie<<endl;
	50	#endif
	51	continue;
	52	}
	53	const TColStd_Array1OfInteger& aNodes1 = aPoly1->Nodes();
	54
	55	// cycle on other polygons
	56	TopTools_ListIteratorOfListOfShape it(aFaces);
	57	it.Next();
	58	for (; it.More(); it.Next()) {
	59	const TopoDS_Face& aFace2 = TopoDS::Face(it.Value());
	60	TopLoc_Location aLoc2;
	61	Handle(Poly_Triangulation) aT2 = BRep_Tool::Triangulation(aFace2, aLoc2);
	62	Handle(Poly_PolygonOnTriangulation) aPoly2 =
	63	BRep_Tool::PolygonOnTriangulation(aEdge, aT2, aLoc2);
	64	if (aPoly2.IsNull() \|\| aT2.IsNull()) {
	65	#ifdef DEB
	66	cout<<"problem getting PolygonOnTriangulation of edge "<<ie<<endl;
	67	#endif
	68	continue;
	69	}
	70	const TColStd_Array1OfInteger& aNodes2 = aPoly2->Nodes();
	71
72	// check equality of polygons lengths
73	if (aNodes2.Length() != aNodes1.Length()) {
74	myAsyncEdges.Append(ie);
75	break;
76	}
77
78	// check distances between corresponding points
79	Standard_Real aDefle = Max(aT1->Deflection(), aT2->Deflection());
80	const TColgp_Array1OfPnt& aPoints1 = aT1->Nodes();
81	const TColgp_Array1OfPnt& aPoints2 = aT2->Nodes();
82	Standard_Integer iF1 = aMapF.FindIndex(aFace1);
83	Standard_Integer iF2 = aMapF.FindIndex(aFace2);
84	Standard_Integer i1 = aNodes1.Lower();
85	Standard_Integer i2 = aNodes2.Lower();
86	gp_Trsf aTrsf1 = aFace1.Location().Transformation();
87	gp_Trsf aTrsf2 = aFace2.Location().Transformation();
88	for (; i1 <= aNodes1.Upper(); i1++, i2++) {
89	gp_Pnt aP1 = aPoints1(aNodes1(i1)).Transformed(aTrsf1);
90	gp_Pnt aP2 = aPoints2(aNodes2(i2)).Transformed(aTrsf2);
91	Standard_Real aDist = aP1.Distance(aP2);
92	if (aDist > aDefle) {
93	myErrors.Append(iF1);
94	myErrors.Append(i1);
95	myErrors.Append(iF2);
96	myErrors.Append(i2);
97	myErrorsVal.Append(aDist);
98	}
99	}
100	}
101	}
102
103	// check triangulations
104	Standard_Integer iF;
105	for (iF=1; iF <= aMapF.Extent(); iF++) {
106	const TopoDS_Face& aFace = TopoDS::Face(aMapF.FindKey(iF));
107	TopLoc_Location aLoc;
108	Handle(Poly_Triangulation) aT = BRep_Tool::Triangulation(aFace, aLoc);
109	if (aT.IsNull()) {
e1bcbb52	110	cout<< "face "<<iF<<" has no triangulation"<<endl;
7fd59977	111	continue;
7fd59977	112	}
446e11f3	113
7fd59977	114	// remember boundary nodes
446e11f3	115	TColStd_PackedMapOfInteger aMapBndNodes;
7fd59977	116	TopExp_Explorer ex(aFace, TopAbs_EDGE);
	117	for (; ex.More(); ex.Next()) {
	118	const TopoDS_Edge& aEdge = TopoDS::Edge(ex.Current());
	119	Handle(Poly_PolygonOnTriangulation) aPoly =
	120	BRep_Tool::PolygonOnTriangulation(aEdge, aT, aLoc);
	121	if (aPoly.IsNull()) continue;
	122	const TColStd_Array1OfInteger& aNodes = aPoly->Nodes();
	123	Standard_Integer i;
	124	for (i=aNodes.Lower(); i <= aNodes.Upper(); i++)
	125	aMapBndNodes.Add(aNodes(i));
	126	}
	127
446e11f3 A	128	TColStd_PackedMapOfInteger aUsedNodes;
	129
	130	// check of free links and nodes
7fd59977	131	Poly_Connect aConn(aT);
	132	const Poly_Array1OfTriangle& aTriangles = aT->Triangles();
	133	Standard_Integer nbTri = aT->NbTriangles(), i, j, n[3], t[3];
	134	for (i = 1; i <= nbTri; i++) {
	135	aTriangles(i).Get(n[0], n[1], n[2]);
446e11f3 A	136
	137	aUsedNodes.Add (n[0]);
	138	aUsedNodes.Add (n[1]);
	139	aUsedNodes.Add (n[2]);
	140
7fd59977	141	aConn.Triangles(i, t[0], t[1], t[2]);
	142	for (j = 0; j < 3; j++) {
	143	if (t[j] == 0) {
	144	// free link found
	145	Standard_Integer k = (j+1) % 3; // the following node of the edge
	146	Standard_Integer n1 = n[j];
	147	Standard_Integer n2 = n[k];
	148	// skip if it is on boundary
	149	if (aMapBndNodes.Contains(n1) && aMapBndNodes.Contains(n2))
	150	continue;
	151	if (!myMapFaceLinks.Contains(iF)) {
	152	//myMapFaceLinks.Add(iF, TColStd_SequenceOfInteger());
	153	TColStd_SequenceOfInteger tmpSeq;
	154	myMapFaceLinks.Add(iF, tmpSeq);
	155	}
	156	TColStd_SequenceOfInteger& aSeq = myMapFaceLinks.ChangeFromKey(iF);
	157	aSeq.Append(n1);
	158	aSeq.Append(n2);
	159	}
	160	}
	161	}
446e11f3 A	162
	163	// check of free nodes
	164	Standard_Integer aNbNodes = aT->NbNodes();
	165	for (Standard_Integer i = 1; i <= aNbNodes; i++)
	166	if ( ! aUsedNodes.Contains(i) )
	167	{
	168	myFreeNodeFaces.Append (iF);
	169	myFreeNodeNums.Append (i);
	170	}
7fd59977	171	}
	172	}
	173
	174	//=======================================================================
	175	//function : GetFreeLink
	176	//purpose : gets the numbers of nodes of a free link with the given index
	177	// in the face with the given index
	178	//=======================================================================
	179
	180	void MeshTest_CheckTopology::GetFreeLink(const Standard_Integer theFaceIndex,
	181	const Standard_Integer theLinkIndex,
	182	Standard_Integer& theNode1,
	183	Standard_Integer& theNode2) const
	184	{
	185	const TColStd_SequenceOfInteger& aSeq = myMapFaceLinks(theFaceIndex);
	186	Standard_Integer aInd = (theLinkIndex-1)*2 + 1;
	187	theNode1 = aSeq(aInd);
	188	theNode2 = aSeq(aInd+1);
	189	}
	190
	191	//=======================================================================
	192	//function : GetCrossFaceError
	193	//purpose : gets the attributes of a cross face error with the given index
	194	//=======================================================================
	195
	196	void MeshTest_CheckTopology::GetCrossFaceError(const Standard_Integer theIndex,
	197	Standard_Integer& theFace1,
	198	Standard_Integer& theNode1,
	199	Standard_Integer& theFace2,
	200	Standard_Integer& theNode2,
	201	Standard_Real& theValue) const
	202	{
	203	Standard_Integer aInd = (theIndex-1)*4 + 1;
	204	theFace1 = myErrors(aInd);
	205	theNode1 = myErrors(aInd+1);
	206	theFace2 = myErrors(aInd+2);
	207	theNode2 = myErrors(aInd+3);
	208	theValue = myErrorsVal(theIndex);
	209	}