From 7d5c34c37c6e0c84799ff91fb82396d0996a4090 Mon Sep 17 00:00:00 2001 From: skl Date: Mon, 26 May 2008 07:24:55 +0000 Subject: [PATCH] Updated for NPAL15716. --- src/SMESH_SWIG/smeshDC.py | 5 +- src/StdMeshers/Makefile.am | 4 +- .../StdMeshers_QuadToTriaAdaptor.cxx | 1161 +++++++++++++++++ .../StdMeshers_QuadToTriaAdaptor.hxx | 75 ++ 4 files changed, 1242 insertions(+), 3 deletions(-) create mode 100644 src/StdMeshers/StdMeshers_QuadToTriaAdaptor.cxx create mode 100644 src/StdMeshers/StdMeshers_QuadToTriaAdaptor.hxx diff --git a/src/SMESH_SWIG/smeshDC.py b/src/SMESH_SWIG/smeshDC.py index 3f77b5cb0..25858311d 100644 --- a/src/SMESH_SWIG/smeshDC.py +++ b/src/SMESH_SWIG/smeshDC.py @@ -823,8 +823,7 @@ class Mesh: def Compute(self, geom=0): if geom == 0 or not isinstance(geom, geompyDC.GEOM._objref_GEOM_Object): if self.geom == 0: - print "Compute impossible: mesh is not constructed on geom shape." - return 0 + geom = self.mesh.GetShapeToMesh() else: geom = self.geom ok = False @@ -931,6 +930,8 @@ class Mesh: pass if not geom: geom = self.geom + if not geom: + geom = self.mesh.GetShapeToMesh() pass status = self.mesh.AddHypothesis(geom, hyp) isAlgo = hyp._narrow( SMESH_Algo ) diff --git a/src/StdMeshers/Makefile.am b/src/StdMeshers/Makefile.am index 74d632a6f..bc886b5df 100644 --- a/src/StdMeshers/Makefile.am +++ b/src/StdMeshers/Makefile.am @@ -63,6 +63,7 @@ salomeinclude_HEADERS = \ StdMeshers_FaceSide.hxx \ StdMeshers_CompositeSegment_1D.hxx \ StdMeshers_UseExisting_1D2D.hxx \ + StdMeshers_QuadToTriaAdaptor.hxx \ SMESH_StdMeshers.hxx # Libraries targets @@ -104,7 +105,8 @@ dist_libStdMeshers_la_SOURCES = \ StdMeshers_SegmentLengthAroundVertex.cxx \ StdMeshers_FaceSide.cxx \ StdMeshers_CompositeSegment_1D.cxx \ - StdMeshers_UseExisting_1D2D.cxx + StdMeshers_UseExisting_1D2D.cxx \ + StdMeshers_QuadToTriaAdaptor.cxx # additionnal information to compil and link file diff --git a/src/StdMeshers/StdMeshers_QuadToTriaAdaptor.cxx b/src/StdMeshers/StdMeshers_QuadToTriaAdaptor.cxx new file mode 100644 index 000000000..52254d698 --- /dev/null +++ b/src/StdMeshers/StdMeshers_QuadToTriaAdaptor.cxx @@ -0,0 +1,1161 @@ +// SMESH SMESH : implementaion of SMESH idl descriptions +// +// Copyright (C) 2003 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN, +// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 of the License. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// +// +// +// File : StdMeshers_QuadToTriaAdaptor.cxx +// Module : SMESH +// Created : Wen May 07 16:37:07 2008 +// Author : Sergey KUUL (skl) + + +#include "StdMeshers_QuadToTriaAdaptor.hxx" + +//#include +//#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +typedef NCollection_Array1 StdMeshers_Array1OfSequenceOfInteger; + + +//======================================================================= +//function : StdMeshers_QuadToTriaAdaptor +//purpose : +//======================================================================= + +StdMeshers_QuadToTriaAdaptor::StdMeshers_QuadToTriaAdaptor() +{ +} + + +//================================================================================ +/*! + * \brief Destructor + */ +//================================================================================ + +StdMeshers_QuadToTriaAdaptor::~StdMeshers_QuadToTriaAdaptor() +{} + + +//======================================================================= +//function : FindBestPoint +//purpose : Auxilare for Compute() +// V - normal to (P1,P2,PC) +//======================================================================= +static gp_Pnt FindBestPoint(const gp_Pnt& P1, const gp_Pnt& P2, + const gp_Pnt& PC, const gp_Vec& V) +{ + double a = P1.Distance(P2); + double b = P1.Distance(PC); + double c = P2.Distance(PC); + if( a < (b+c)/2 ) + return PC; + else { + // find shift along V in order to a became equal to (b+c)/2 + double shift = sqrt( a*a + (b*b-c*c)*(b*b-c*c)/16/a/a - (b*b+c*c)/2 ); + gp_Dir aDir(V); + gp_Pnt Pbest( PC.X() + aDir.X()*shift, PC.Y() + aDir.Y()*shift, + PC.Z() + aDir.Z()*shift ); + return Pbest; + } +} + + +//======================================================================= +//function : HasIntersection3 +//purpose : Auxilare for HasIntersection() +// find intersection point between triangle (P1,P2,P3) +// and segment [PC,P] +//======================================================================= +static bool HasIntersection3(const gp_Pnt& P, const gp_Pnt& PC, gp_Pnt& Pint, + const gp_Pnt& P1, const gp_Pnt& P2, const gp_Pnt& P3) +{ + //cout<<"HasIntersection3"< preci ) || + ( (PC.Y()-PIn.Y())*(P.Y()-PIn.Y()) > preci ) || + ( (PC.Z()-PIn.Z())*(P.Z()-PIn.Z()) > preci ); + if(IsExternal) { + return false; + } + // check if this point is internal for triangle (P1,P2,P3) + gp_Vec V1(PIn,P1); + gp_Vec V2(PIn,P2); + gp_Vec V3(PIn,P3); + if( V1.Magnitude()Length()==3) { + return HasIntersection3( P, PC, Pint, aContour->Value(1), + aContour->Value(2), aContour->Value(3) ); + } + else { + bool check = false; + if( (aContour->Value(1).Distance(aContour->Value(2)) > 1.e-6) && + (aContour->Value(1).Distance(aContour->Value(3)) > 1.e-6) && + (aContour->Value(2).Distance(aContour->Value(3)) > 1.e-6) ) { + check = HasIntersection3( P, PC, Pint, aContour->Value(1), + aContour->Value(2), aContour->Value(3) ); + } + if(check) return true; + if( (aContour->Value(1).Distance(aContour->Value(4)) > 1.e-6) && + (aContour->Value(1).Distance(aContour->Value(3)) > 1.e-6) && + (aContour->Value(4).Distance(aContour->Value(3)) > 1.e-6) ) { + check = HasIntersection3( P, PC, Pint, aContour->Value(1), + aContour->Value(3), aContour->Value(4) ); + } + if(check) return true; + } + + return false; +} + + +//======================================================================= +//function : CheckIntersection +//purpose : Auxilare for Compute() +// NotCheckedFace - for optimization +//======================================================================= +bool StdMeshers_QuadToTriaAdaptor::CheckIntersection + (const gp_Pnt& P, const gp_Pnt& PC, + gp_Pnt& Pint, SMESH_Mesh& aMesh, + const TopoDS_Shape& aShape, + const TopoDS_Shape& NotCheckedFace) +{ + SMESHDS_Mesh * meshDS = aMesh.GetMeshDS(); + //cout<<" CheckIntersection: meshDS->NbFaces() = "<NbFaces()<MeshElements(aShapeFace); + if ( aSubMeshDSFace ) { + SMDS_ElemIteratorPtr iteratorElem = aSubMeshDSFace->GetElements(); + while ( iteratorElem->more() ) { // loop on elements on a face + const SMDS_MeshElement* face = iteratorElem->next(); + Handle(TColgp_HSequenceOfPnt) aContour = new TColgp_HSequenceOfPnt; + SMDS_ElemIteratorPtr nodeIt = face->nodesIterator(); + if( !face->IsQuadratic() ) { + while ( nodeIt->more() ) { + const SMDS_MeshNode* node = static_cast( nodeIt->next() ); + aContour->Append(gp_Pnt(node->X(), node->Y(), node->Z())); + } + } + else { + int nn = 0; + while ( nodeIt->more() ) { + nn++; + const SMDS_MeshNode* node = static_cast( nodeIt->next() ); + aContour->Append(gp_Pnt(node->X(), node->Y(), node->Z())); + if(nn==face->NbNodes()/2) break; + } + } + if( HasIntersection(P, PC, Pres, aContour) ) { + res = true; + double tmp = PC.Distance(Pres); + if(tmpnodesIterator(); + const SMDS_MeshNode* Ns1[3]; + int k = 0; + while( nIt->more() ) { + Ns1[k] = static_cast( nIt->next() ); + k++; + } + nIt = F2->nodesIterator(); + const SMDS_MeshNode* Ns2[3]; + k = 0; + while( nIt->more() ) { + Ns2[k] = static_cast( nIt->next() ); + k++; + } + if( ( Ns1[1]==Ns2[1] && Ns1[2]==Ns2[2] ) || + ( Ns1[1]==Ns2[2] && Ns1[2]==Ns2[1] ) ) + return true; + return false; +} + + +//======================================================================= +//function : IsDegenarate +//purpose : Auxilare for Preparation() +//======================================================================= +static int IsDegenarate(const Handle(TColgp_HArray1OfPnt)& PN) +{ + int i = 1; + for(; i<4; i++) { + int j = i+1; + for(; j<=4; j++) { + if( PN->Value(i).Distance(PN->Value(j)) < 1.e-6 ) + return j; + } + } + return 0; +} + + +//======================================================================= +//function : Preparation +//purpose : Auxilare for Compute() +// : Return 0 if given face is not quad, +// 1 if given face is quad, +// 2 if given face is degenerate quad (two nodes are coincided) +//======================================================================= +int StdMeshers_QuadToTriaAdaptor::Preparation(const SMDS_MeshElement* face, + Handle(TColgp_HArray1OfPnt) PN, + Handle(TColgp_HArray1OfVec) VN, + std::vector& FNodes, + gp_Pnt& PC, gp_Vec& VNorm) +{ + int i = 0; + double xc=0., yc=0., zc=0.; + SMDS_ElemIteratorPtr nodeIt = face->nodesIterator(); + if( !face->IsQuadratic() ) { + if( face->NbNodes() != 4 ) + return 0; + while ( nodeIt->more() ) { + i++; + const SMDS_MeshNode* node = static_cast( nodeIt->next() ); + FNodes[i-1] = node; + PN->SetValue( i, gp_Pnt(node->X(), node->Y(), node->Z()) ); + xc += node->X(); + yc += node->Y(); + zc += node->Z(); + } + } + else { + if( face->NbNodes() != 8) + return 0; + while ( nodeIt->more() ) { + i++; + const SMDS_MeshNode* node = static_cast( nodeIt->next() ); + FNodes[i-1] = node; + PN->SetValue( i, gp_Pnt(node->X(), node->Y(), node->Z()) ); + xc += node->X(); + yc += node->Y(); + zc += node->Z(); + if(i==4) break; + } + } + + int nbp = 4; + + int j = 0; + for(i=1; i<4; i++) { + j = i+1; + for(; j<=4; j++) { + if( PN->Value(i).Distance(PN->Value(j)) < 1.e-6 ) + break; + } + if(j<=4) break; + } + //int deg_num = IsDegenarate(PN); + //if(deg_num>0) { + bool hasdeg = false; + if(i<4) { + //cout<<"find degeneration"<Value(i); + + std::list< const SMDS_MeshNode* >::iterator itdg = myDegNodes.begin(); + const SMDS_MeshNode* DegNode = 0; + for(; itdg!=myDegNodes.end(); itdg++) { + const SMDS_MeshNode* N = (*itdg); + gp_Pnt Ptmp(N->X(),N->Y(),N->Z()); + if(Pdeg.Distance(Ptmp)<1.e-6) { + DegNode = N; + //DegNode = const_cast(N); + break; + } + } + if(!DegNode) { + DegNode = FNodes[i-1]; + myDegNodes.push_back(DegNode); + } + else { + FNodes[i-1] = DegNode; + } + for(i=j; i<4; i++) { + PN->SetValue(i,PN->Value(i+1)); + FNodes[i-1] = FNodes[i]; + } + nbp = 3; + //PC = gp_Pnt( PN->Value(1).X() + PN.Value + } + + PC = gp_Pnt(xc/4., yc/4., zc/4.); + //cout<<" PC("<SetValue(5,PN->Value(1)); + PN->SetValue(nbp+1,PN->Value(1)); + //FNodes[4] = FNodes[0]; + FNodes[nbp] = FNodes[0]; + // find normal direction + //gp_Vec V1(PC,PN->Value(4)); + gp_Vec V1(PC,PN->Value(nbp)); + gp_Vec V2(PC,PN->Value(1)); + VNorm = V1.Crossed(V2); + //VN->SetValue(4,VNorm); + VN->SetValue(nbp,VNorm); + //for(i=1; i<4; i++) { + for(i=1; iValue(i)); + V2 = gp_Vec(PC,PN->Value(i+1)); + gp_Vec Vtmp = V1.Crossed(V2); + VN->SetValue(i,Vtmp); + VNorm += Vtmp; + } + //cout<<" VNorm("<MeshElements( aShapeFace ); + if ( aSubMeshDSFace ) { + bool isRev = SMESH_Algo::IsReversedSubMesh( TopoDS::Face(aShapeFace), meshDS ); + + SMDS_ElemIteratorPtr iteratorElem = aSubMeshDSFace->GetElements(); + while ( iteratorElem->more() ) { // loop on elements on a face + const SMDS_MeshElement* face = iteratorElem->next(); + //cout<GetID() = "<GetID()< FNodes(5); + gp_Pnt PC; + gp_Vec VNorm; + int stat = Preparation(face, PN, VN, FNodes, PC, VNorm); + if(stat==0) + continue; + + if(stat==2) { + // degenerate face + // add triangles to result map + std::list aList; + SMDS_FaceOfNodes* NewFace; + if(!isRev) + NewFace = new SMDS_FaceOfNodes( FNodes[0], FNodes[1], FNodes[2] ); + else + NewFace = new SMDS_FaceOfNodes( FNodes[0], FNodes[2], FNodes[1] ); + aList.push_back(NewFace); + myResMap.insert(make_pair(face,aList)); + continue; + } + + if(!isRev) VNorm.Reverse(); + double xc = 0., yc = 0., zc = 0.; + int i = 1; + for(; i<=4; i++) { + gp_Pnt Pbest; + if(!isRev) + Pbest = FindBestPoint(PN->Value(i), PN->Value(i+1), PC, VN->Value(i).Reversed()); + else + Pbest = FindBestPoint(PN->Value(i), PN->Value(i+1), PC, VN->Value(i)); + xc += Pbest.X(); + yc += Pbest.Y(); + zc += Pbest.Z(); + } + gp_Pnt PCbest(xc/4., yc/4., zc/4.); + + // check PCbest + double height = PCbest.Distance(PC); + if(height<1.e-6) { + // create new PCbest using a bit shift along VNorm + PCbest = gp_Pnt( PC.X() + VNorm.X()*0.001, + PC.Y() + VNorm.Y()*0.001, + PC.Z() + VNorm.Z()*0.001); + } + else { + // check possible intersection with other faces + gp_Pnt Pint; + bool check = CheckIntersection(PCbest, PC, Pint, aMesh, aShape, aShapeFace); + if(check) { + //cout<<"--PC("<AddNode( PCbest.X(), PCbest.Y(), PCbest.Z() ); + // add triangles to result map + std::list aList; + for(i=0; i<4; i++) { + SMDS_FaceOfNodes* NewFace = new SMDS_FaceOfNodes( NewNode, FNodes[i], FNodes[i+1] ); + aList.push_back(NewFace); + } + myResMap.insert(make_pair(face,aList)); + // create pyramid + SMDS_MeshVolume* aPyram = + meshDS->AddVolume( FNodes[0], FNodes[1], FNodes[2], FNodes[3], NewNode ); + myMapFPyram.insert(make_pair(face,aPyram)); + } // end loop on elements on a face + } + } // end for(TopExp_Explorer exp(aShape,TopAbs_FACE);exp.More();exp.Next()) { + + return Compute2ndPart(aMesh); +} + + +//======================================================================= +//function : Compute +//purpose : +//======================================================================= + +bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh) +{ + myResMap.clear(); + myMapFPyram.clear(); + + SMESHDS_Mesh * meshDS = aMesh.GetMeshDS(); + + SMDS_FaceIteratorPtr itFace = meshDS->facesIterator(); + + while(itFace->more()) { + const SMDS_MeshElement* face = itFace->next(); + if ( !face ) continue; + //cout<GetID() = "<GetID()< FNodes(5); + gp_Pnt PC; + gp_Vec VNorm; + + int stat = Preparation(face, PN, VN, FNodes, PC, VNorm); + if(stat==0) + continue; + + if(stat==2) { + // degenerate face + // add triangles to result map + std::list aList; + SMDS_FaceOfNodes* NewFace; + // check orientation + + double tmp = PN->Value(1).Distance(PN->Value(2)) + + PN->Value(2).Distance(PN->Value(3)); + gp_Dir tmpDir(VNorm); + gp_Pnt Ptmp1( PC.X() + tmpDir.X()*tmp*1.e6, + PC.Y() + tmpDir.Y()*tmp*1.e6, + PC.Z() + tmpDir.Z()*tmp*1.e6 ); + gp_Pnt Ptmp2( PC.X() + tmpDir.Reversed().X()*tmp*1.e6, + PC.Y() + tmpDir.Reversed().Y()*tmp*1.e6, + PC.Z() + tmpDir.Reversed().Z()*tmp*1.e6 ); + // check intersection for Ptmp1 and Ptmp2 + bool IsRev = false; + bool IsOK1 = false; + bool IsOK2 = false; + double dist1 = RealLast(); + double dist2 = RealLast(); + gp_Pnt Pres1,Pres2; + SMDS_FaceIteratorPtr itf = meshDS->facesIterator(); + while(itf->more()) { + const SMDS_MeshElement* F = itf->next(); + if(F==face) continue; + Handle(TColgp_HSequenceOfPnt) aContour = new TColgp_HSequenceOfPnt; + SMDS_ElemIteratorPtr nodeIt = F->nodesIterator(); + if( !F->IsQuadratic() ) { + while ( nodeIt->more() ) { + const SMDS_MeshNode* node = static_cast( nodeIt->next() ); + aContour->Append(gp_Pnt(node->X(), node->Y(), node->Z())); + } + } + else { + int nn = 0; + while ( nodeIt->more() ) { + nn++; + const SMDS_MeshNode* node = static_cast( nodeIt->next() ); + aContour->Append(gp_Pnt(node->X(), node->Y(), node->Z())); + if(nn==face->NbNodes()/2) break; + } + } + gp_Pnt PPP; + if( HasIntersection(Ptmp1, PC, PPP, aContour) ) { + IsOK1 = true; + double tmp = PC.Distance(PPP); + if(tmpValue(i), PN->Value(i+1), PC, VN->Value(i)); + xc += Pbest.X(); + yc += Pbest.Y(); + zc += Pbest.Z(); + } + gp_Pnt PCbest(xc/4., yc/4., zc/4.); + double height = PCbest.Distance(PC); + if(height<1.e-6) { + // create new PCbest using a bit shift along VNorm + PCbest = gp_Pnt( PC.X() + VNorm.X()*0.001, + PC.Y() + VNorm.Y()*0.001, + PC.Z() + VNorm.Z()*0.001); + height = PCbest.Distance(PC); + } + //cout<<" PCbest("<Value(1).Distance(PN->Value(3)) + + PN->Value(2).Distance(PN->Value(4)); + gp_Dir tmpDir(V1); + gp_Pnt Ptmp1( PC.X() + tmpDir.X()*tmp*1.e6, + PC.Y() + tmpDir.Y()*tmp*1.e6, + PC.Z() + tmpDir.Z()*tmp*1.e6 ); + gp_Pnt Ptmp2( PC.X() + tmpDir.Reversed().X()*tmp*1.e6, + PC.Y() + tmpDir.Reversed().Y()*tmp*1.e6, + PC.Z() + tmpDir.Reversed().Z()*tmp*1.e6 ); + // check intersection for Ptmp1 and Ptmp2 + bool IsRev = false; + bool IsOK1 = false; + bool IsOK2 = false; + double dist1 = RealLast(); + double dist2 = RealLast(); + gp_Pnt Pres1,Pres2; + SMDS_FaceIteratorPtr itf = meshDS->facesIterator(); + while(itf->more()) { + const SMDS_MeshElement* F = itf->next(); + if(F==face) continue; + Handle(TColgp_HSequenceOfPnt) aContour = new TColgp_HSequenceOfPnt; + SMDS_ElemIteratorPtr nodeIt = F->nodesIterator(); + if( !F->IsQuadratic() ) { + while ( nodeIt->more() ) { + const SMDS_MeshNode* node = static_cast( nodeIt->next() ); + aContour->Append(gp_Pnt(node->X(), node->Y(), node->Z())); + } + } + else { + int nn = 0; + while ( nodeIt->more() ) { + nn++; + const SMDS_MeshNode* node = static_cast( nodeIt->next() ); + aContour->Append(gp_Pnt(node->X(), node->Y(), node->Z())); + if(nn==face->NbNodes()/2) break; + } + } + gp_Pnt PPP; + if( HasIntersection(Ptmp1, PC, PPP, aContour) ) { + IsOK1 = true; + double tmp = PC.Distance(PPP); + if(tmp tmp ) { + height = tmp; + PCbest = gp_Pnt( PC.X() + tmpDir.X()*height, + PC.Y() + tmpDir.Y()*height, + PC.Z() + tmpDir.Z()*height ); + } + } + else if( !IsOK1 && IsOK2 ) { + // using opposite direction + IsRev = true; + double tmp = PC.Distance(Pres2)/3.; + if( height > tmp ) height = tmp; + PCbest = gp_Pnt( PC.X() + tmpDir.Reversed().X()*height, + PC.Y() + tmpDir.Reversed().Y()*height, + PC.Z() + tmpDir.Reversed().Z()*height ); + } + else { // IsOK1 && IsOK2 + double tmp1 = PC.Distance(Pres1)/3.; + double tmp2 = PC.Distance(Pres2)/3.; + if(tmp1 tmp1 ) { + height = tmp1; + PCbest = gp_Pnt( PC.X() + tmpDir.X()*height, + PC.Y() + tmpDir.Y()*height, + PC.Z() + tmpDir.Z()*height ); + } + } + else { + // using opposite direction + IsRev = true; + if( height > tmp2 ) height = tmp2; + PCbest = gp_Pnt( PC.X() + tmpDir.Reversed().X()*height, + PC.Y() + tmpDir.Reversed().Y()*height, + PC.Z() + tmpDir.Reversed().Z()*height ); + } + } + + // create node for PCbest + SMDS_MeshNode* NewNode = meshDS->AddNode( PCbest.X(), PCbest.Y(), PCbest.Z() ); + // add triangles to result map + std::list aList; + for(i=0; i<4; i++) { + SMDS_FaceOfNodes* NewFace; + if(IsRev) + NewFace = new SMDS_FaceOfNodes( NewNode, FNodes[i], FNodes[i+1] ); + else + NewFace = new SMDS_FaceOfNodes( NewNode, FNodes[i+1], FNodes[i] ); + aList.push_back(NewFace); + } + myResMap.insert(make_pair(face,aList)); + // create pyramid + SMDS_MeshVolume* aPyram; + if(IsRev) + aPyram = meshDS->AddVolume( FNodes[0], FNodes[1], FNodes[2], FNodes[3], NewNode ); + else + aPyram = meshDS->AddVolume( FNodes[0], FNodes[3], FNodes[2], FNodes[1], NewNode ); + myMapFPyram.insert(make_pair(face,aPyram)); + } // end loop on elements on a face + + return Compute2ndPart(aMesh); +} + + +//======================================================================= +//function : Compute2ndPart +//purpose : +//======================================================================= + +bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh& aMesh) +{ + SMESHDS_Mesh * meshDS = aMesh.GetMeshDS(); + + // check intersections between created pyramids + int NbPyram = myMapFPyram.size(); + //cout<<"NbPyram = "< Pyrams(NbPyram); + std::vector< const SMDS_MeshElement* > Faces(NbPyram); + std::map< const SMDS_MeshElement*, + const SMDS_MeshElement* >::iterator itp = myMapFPyram.begin(); + int i = 0; + for(; itp!=myMapFPyram.end(); itp++, i++) { + Faces[i] = (*itp).first; + Pyrams[i] = (*itp).second; + } + StdMeshers_Array1OfSequenceOfInteger MergesInfo(0,NbPyram-1); + for(i=0; inodesIterator(); + std::vector Ps1(5); + const SMDS_MeshNode* Ns1[5]; + int k = 0; + while( nIt->more() ) { + const SMDS_MeshNode* node = static_cast( nIt->next() ); + Ns1[k] = node; + Ps1[k] = gp_Pnt(node->X(), node->Y(), node->Z()); + k++; + } + bool NeedMove = false; + for(int j=i+1; jChangeElementNodes(Prm2, Ns2, 5); + // update pyramids for J + for(k=2; k<=nbJ; k++) { + const SMDS_MeshElement* tmpPrm = Pyrams[aMergesJ.Value(k)]; + SMDS_ElemIteratorPtr tmpIt = tmpPrm->nodesIterator(); + const SMDS_MeshNode* Ns[5]; + int m = 0; + while( tmpIt->more() ) { + Ns[m] = static_cast( tmpIt->next() ); + m++; + } + Ns[4] = CommonNode; + meshDS->ChangeElementNodes(tmpPrm, Ns, 5); + } + + // update MergesInfo + for(k=1; k<=nbI; k++) { + int num = aMergesI.Value(k); + const TColStd_SequenceOfInteger& aSeq = MergesInfo.Value(num); + TColStd_SequenceOfInteger tmpSeq; + int m = 1; + for(; m<=aSeq.Length(); m++) { + tmpSeq.Append(aSeq.Value(m)); + } + for(m=1; m<=nbJ; m++) { + tmpSeq.Append(aMergesJ.Value(m)); + } + MergesInfo.SetValue(num,tmpSeq); + } + for(k=1; k<=nbJ; k++) { + int num = aMergesJ.Value(k); + const TColStd_SequenceOfInteger& aSeq = MergesInfo.Value(num); + TColStd_SequenceOfInteger tmpSeq; + int m = 1; + for(; m<=aSeq.Length(); m++) { + tmpSeq.Append(aSeq.Value(m)); + } + for(m=1; m<=nbI; m++) { + tmpSeq.Append(aMergesI.Value(m)); + } + MergesInfo.SetValue(num,tmpSeq); + } + + // update triangles for aMergesJ + for(k=1; k<=nbJ; k++) { + std::list< std::list< const SMDS_MeshNode* > > aFNodes; + std::list< const SMDS_MeshElement* > aFFaces; + int num = aMergesJ.Value(k); + std::map< const SMDS_MeshElement*, + std::list >::iterator itrm = myResMap.find(Faces[num]); + std::list trias = (*itrm).second; + std::list::iterator itt = trias.begin(); + for(; itt!=trias.end(); itt++) { + int nn = -1; + SMDS_ElemIteratorPtr nodeIt = (*itt)->nodesIterator(); + const SMDS_MeshNode* NF[3]; + while ( nodeIt->more() ) { + nn++; + NF[nn] = static_cast( nodeIt->next() ); + } + NF[0] = CommonNode; + SMDS_FaceOfNodes* Ftria = const_cast< SMDS_FaceOfNodes*>( (*itt) ); + Ftria->ChangeNodes(NF, 3); + } + } + + // check and remove coincided faces + TColStd_SequenceOfInteger IdRemovedTrias; + int i1 = 1; + for(; i1<=nbI; i1++) { + int numI = aMergesI.Value(i1); + std::map< const SMDS_MeshElement*, + std::list >::iterator itrmI = myResMap.find(Faces[numI]); + std::list triasI = (*itrmI).second; + std::list::iterator ittI = triasI.begin(); + int nbfI = triasI.size(); + const SMDS_FaceOfNodes* FsI[nbfI]; + k = 0; + for(; ittI!=triasI.end(); ittI++) { + FsI[k] = (*ittI); + k++; + } + int i2 = 0; + for(; i2 >::iterator itrmJ = myResMap.find(Faces[numJ]); + std::list triasJ = (*itrmJ).second; + std::list::iterator ittJ = triasJ.begin(); + int nbfJ = triasJ.size(); + const SMDS_FaceOfNodes* FsJ[nbfJ]; + k = 0; + for(; ittJ!=triasJ.end(); ittJ++) { + FsJ[k] = (*ittJ); + k++; + } + int j2 = 0; + for(; j2GetID() ); + IdRemovedTrias.Append( FJ->GetID() ); + FsI[i2] = 0; + FsJ[j2] = 0; + std::list new_triasI; + for(k=0; k new_triasJ; + for(k=0; kRemoveNode(Nrem); + } + else { // nbc==0 + //cout<<"decrease height of pyramids"<PI/3.) + h1 = VI1.Magnitude()/2; + else + h1 = VI1.Magnitude()*cos(ang1); + if(ang2>PI/3.) + h2 = VI2.Magnitude()/2; + else + h2 = VI2.Magnitude()*cos(ang2); + double coef1 = 0.5; + if(ang1(Ns1[4]); + VN1.Scale(coef1); + aNode1->setXYZ( PC1.X()+VN1.X(), PC1.Y()+VN1.Y(), PC1.Z()+VN1.Z() ); + SMDS_MeshNode* aNode2 = const_cast(Ns2[4]); + VN2.Scale(coef2); + aNode2->setXYZ( PC2.X()+VN2.X(), PC2.Y()+VN2.Y(), PC2.Z()+VN2.Z() ); + NeedMove = true; + } + } // end if(hasInt) + else { + //cout<<" no intersec for i="<