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

// Created on: 2004-05-10
// Created by: Michael SAZONOV
// Copyright (c) 2004-2014 OPEN CASCADE SAS
//
// This file is part of Open CASCADE Technology software library.
//
// This library is free software; you can redistribute it and/or modify it under
// the terms of the GNU Lesser General Public License version 2.1 as published
// by the Free Software Foundation, with special exception defined in the file
// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
// distribution for complete text of the license and disclaimer of any warranty.
//
// Alternatively, this file may be used under the terms of Open CASCADE
// commercial license or contractual agreement.

#include <MeshTest_CheckTopology.hxx>
#include <BRep_Tool.hxx>
#include <TopExp.hxx>
#include <TopExp_Explorer.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 <Poly_Triangulation.hxx>
#include <Poly_PolygonOnTriangulation.hxx>
#include <Poly_Connect.hxx>
#include <Precision.hxx>

//=======================================================================
//function : ComputeArea
//purpose  : Computes area of the triangle given by its three points (either 2D or3D)
//=======================================================================
static Standard_Real ComputeArea(const gp_XYZ& theP1,
                                 const gp_XYZ& theP2,
                                 const gp_XYZ& theP3)
{
  return 0.5*(theP3 - theP1).Crossed(theP2 - theP1).Modulus();
}

//=======================================================================
//function : ComputeArea
//purpose  : Computes area of the triangle given by its three points (either 2D or3D)
//=======================================================================
static Standard_Real ComputeArea(const gp_XY& theP1,
                                 const gp_XY& theP2,
                                 const gp_XY& theP3)
{
  return 0.5*Abs((theP3 - theP1).Crossed(theP2 - theP1));
}

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

void MeshTest_CheckTopology::Perform (Draw_Interpretor& di)
{
  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 OCCT_DEBUG
      std::cout<<"problem getting PolygonOnTriangulation of edge "<<ie<<std::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 OCCT_DEBUG
	std::cout<<"problem getting PolygonOnTriangulation of edge "<<ie<<std::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 aSqDefle = BRep_Tool::Tolerance(aEdge);
      aSqDefle *= aSqDefle;
      Standard_Integer iF1 = aMapF.FindIndex(aFace1);
      Standard_Integer iF2 = aMapF.FindIndex(aFace2);
      Standard_Integer i1 = aNodes1.Lower();
      Standard_Integer i2 = aNodes2.Lower();
      const gp_Trsf &aTrsf1 = aFace1.Location().Transformation();
      const gp_Trsf &aTrsf2 = aFace2.Location().Transformation();
      for (; i1 <= aNodes1.Upper(); i1++, i2++) {
        const gp_Pnt aP1 = aT1->Node (aNodes1[i1]).Transformed (aTrsf1);
        const gp_Pnt aP2 = aT2->Node (aNodes2[i2]).Transformed (aTrsf2);
        const Standard_Real aSqDist = aP1.SquareDistance(aP2);
        if (aSqDist > aSqDefle)
        {
          myErrors.Append(iF1);
          myErrors.Append(i1);
          myErrors.Append(iF2);
          myErrors.Append(i2);
          myErrorsVal.Append(Sqrt(aSqDist));
        }
      }
    }
  }

  // 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()) {
      di << "face " <<iF <<" has no triangulation\n";
      continue;
    }

    const gp_Trsf &aTrsf = aLoc.Transformation();

    // 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);
    Standard_Integer nbTri = aT->NbTriangles(), i, j, n[3], t[3];
    for (i = 1; i <= nbTri; i++) {
      aT->Triangle (i).Get (n[0], n[1], n[2]);
      
      aUsedNodes.Add (n[0]);
      aUsedNodes.Add (n[1]);
      aUsedNodes.Add (n[2]);

      const gp_Pnt aPts[3] = {aT->Node(n[0]).Transformed(aTrsf),
                              aT->Node(n[1]).Transformed(aTrsf),
                              aT->Node(n[2]).Transformed(aTrsf)};

      Standard_Real anArea = ComputeArea(aPts[0].XYZ(), aPts[1].XYZ(), aPts[2].XYZ());
      if (anArea < Precision::SquareConfusion())
      {
        mySmallTrianglesFaces.Append(iF);
        mySmallTrianglesTriangles.Append(i);
      }
      else if (aT->HasUVNodes())
      {
        const gp_XY aPUV[3] = {aT->UVNode(n[0]).XY(),
                               aT->UVNode(n[1]).XY(),
                               aT->UVNode(n[2]).XY()};
        anArea = ComputeArea(aPUV[0], aPUV[1], aPUV[2]);
        if (anArea < Precision::SquarePConfusion())
        {
          mySmallTrianglesFaces.Append(iF);
          mySmallTrianglesTriangles.Append(i);
        }
      }

      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)) {
            Handle(TColStd_HSequenceOfInteger) tmpSeq = new TColStd_HSequenceOfInteger;
	    myMapFaceLinks.Add(iF, tmpSeq);
          }
	  Handle(TColStd_HSequenceOfInteger)& aSeq = myMapFaceLinks.ChangeFromKey(iF);
	  aSeq->Append(n1);
	  aSeq->Append(n2);
	}
      }
    }
    
    // check of free nodes
    Standard_Integer aNbNodes = aT->NbNodes();
    for (Standard_Integer k = 1; k <= aNbNodes; k++)
      if ( ! aUsedNodes.Contains(k) )
      {
	myFreeNodeFaces.Append (iF);
	myFreeNodeNums.Append (k);
      }
  }
}

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