X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;ds=sidebyside;f=src%2FStdMeshers%2FStdMeshers_QuadToTriaAdaptor.cxx;h=54c3b7609f79c418c945cc1454a0095b68c07c0a;hb=4cd2499bddcd3da3ec8900fe825bc98669b789b5;hp=1ca7e5b1283995c4f9f1002312b5f00df910c0f8;hpb=9357f5c87098aff2b95b754d69f66c76d2df9c24;p=modules%2Fsmesh.git diff --git a/src/StdMeshers/StdMeshers_QuadToTriaAdaptor.cxx b/src/StdMeshers/StdMeshers_QuadToTriaAdaptor.cxx index 1ca7e5b12..54c3b7609 100644 --- a/src/StdMeshers/StdMeshers_QuadToTriaAdaptor.cxx +++ b/src/StdMeshers/StdMeshers_QuadToTriaAdaptor.cxx @@ -1,32 +1,35 @@ -// Copyright (C) 2007-2010 CEA/DEN, EDF R&D, OPEN CASCADE +// Copyright (C) 2007-2012 CEA/DEN, EDF R&D, OPEN CASCADE // -// 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 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. +// 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 +// 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 +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com // -// SMESH SMESH : implementaion of SMESH idl descriptions // 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 "SMDS_SetIterator.hxx" + +#include "SMESH_Algo.hxx" +#include "SMESH_MesherHelper.hxx" +#include "SMESH_Group.hxx" +#include "SMESHDS_GroupBase.hxx" #include #include @@ -37,40 +40,348 @@ #include #include #include +#include "utilities.h" +#include #include +#include using namespace std; -enum EQuadNature { NOT_QUAD, QUAD, DEGEN_QUAD }; +enum EQuadNature { NOT_QUAD, QUAD, DEGEN_QUAD, PYRAM_APEX = 4, TRIA_APEX = 0 }; + +// std-like iterator used to get coordinates of nodes of mesh element +typedef SMDS_StdIterator< SMESH_TNodeXYZ, SMDS_ElemIteratorPtr > TXyzIterator; + +namespace +{ + //================================================================================ + /*! + * \brief Return true if two nodes of triangles are equal + */ + //================================================================================ + + bool EqualTriangles(const SMDS_MeshElement* F1,const SMDS_MeshElement* F2) + { + return + ( F1->GetNode(1)==F2->GetNode(2) && F1->GetNode(2)==F2->GetNode(1) ) || + ( F1->GetNode(1)==F2->GetNode(1) && F1->GetNode(2)==F2->GetNode(2) ); + } + //================================================================================ + /*! + * \brief Return true if two adjacent pyramids are too close one to another + * so that a tetrahedron to built between them would have too poor quality + */ + //================================================================================ + + bool TooCloseAdjacent( const SMDS_MeshElement* PrmI, + const SMDS_MeshElement* PrmJ, + const bool hasShape) + { + const SMDS_MeshNode* nApexI = PrmI->GetNode(4); + const SMDS_MeshNode* nApexJ = PrmJ->GetNode(4); + if ( nApexI == nApexJ || + nApexI->getshapeId() != nApexJ->getshapeId() ) + return false; + + // Find two common base nodes and their indices within PrmI and PrmJ + const SMDS_MeshNode* baseNodes[2] = { 0,0 }; + int baseNodesIndI[2], baseNodesIndJ[2]; + for ( int i = 0; i < 4 ; ++i ) + { + int j = PrmJ->GetNodeIndex( PrmI->GetNode(i)); + if ( j >= 0 ) + { + int ind = baseNodes[0] ? 1:0; + if ( baseNodes[ ind ]) + return false; // pyramids with a common base face + baseNodes [ ind ] = PrmI->GetNode(i); + baseNodesIndI[ ind ] = i; + baseNodesIndJ[ ind ] = j; + } + } + if ( !baseNodes[1] ) return false; // not adjacent + + // Get normals of triangles sharing baseNodes + gp_XYZ apexI = SMESH_TNodeXYZ( nApexI ); + gp_XYZ apexJ = SMESH_TNodeXYZ( nApexJ ); + gp_XYZ base1 = SMESH_TNodeXYZ( baseNodes[0]); + gp_XYZ base2 = SMESH_TNodeXYZ( baseNodes[1]); + gp_Vec baseVec( base1, base2 ); + gp_Vec baI( base1, apexI ); + gp_Vec baJ( base1, apexJ ); + gp_Vec nI = baseVec.Crossed( baI ); + gp_Vec nJ = baseVec.Crossed( baJ ); + + // Check angle between normals + double angle = nI.Angle( nJ ); + bool tooClose = ( angle < 15. * M_PI / 180. ); + + // Check if pyramids collide + if ( !tooClose && baI * baJ > 0 ) + { + // find out if nI points outside of PrmI or inside + int dInd = baseNodesIndI[1] - baseNodesIndI[0]; + bool isOutI = ( abs(dInd)==1 ) ? dInd < 0 : dInd > 0; + + // find out sign of projection of nJ to baI + double proj = baI * nJ; + + tooClose = isOutI ? proj > 0 : proj < 0; + } + + // Check if PrmI and PrmJ are in same domain + if ( tooClose && !hasShape ) + { + // check order of baseNodes within pyramids, it must be opposite + int dInd; + dInd = baseNodesIndI[1] - baseNodesIndI[0]; + bool isOutI = ( abs(dInd)==1 ) ? dInd < 0 : dInd > 0; + dInd = baseNodesIndJ[1] - baseNodesIndJ[0]; + bool isOutJ = ( abs(dInd)==1 ) ? dInd < 0 : dInd > 0; + if ( isOutJ == isOutI ) + return false; // other domain + + // direct both normals outside pyramid + ( isOutI ? nJ : nI ).Reverse(); + + // check absence of a face separating domains between pyramids + TIDSortedElemSet emptySet, avoidSet; + int i1, i2; + while ( const SMDS_MeshElement* f = + SMESH_MeshEditor::FindFaceInSet( baseNodes[0], baseNodes[1], + emptySet, avoidSet, &i1, &i2 )) + { + avoidSet.insert( f ); + + // face node other than baseNodes + int otherNodeInd = 0; + while ( otherNodeInd == i1 || otherNodeInd == i2 ) otherNodeInd++; + const SMDS_MeshNode* otherFaceNode = f->GetNode( otherNodeInd ); + + if ( otherFaceNode == nApexI || otherFaceNode == nApexJ ) + continue; // f is a temporary triangle + + // check if f is a base face of either of pyramids + if ( f->NbCornerNodes() == 4 && + ( PrmI->GetNodeIndex( otherFaceNode ) >= 0 || + PrmJ->GetNodeIndex( otherFaceNode ) >= 0 )) + continue; // f is a base quadrangle + + // check projections of face direction (baOFN) to triange normals (nI and nJ) + gp_Vec baOFN( base1, SMESH_TNodeXYZ( otherFaceNode )); + if ( nI * baOFN > 0 && nJ * baOFN > 0 ) + { + tooClose = false; // f is between pyramids + break; + } + } + } + + return tooClose; + } + + //================================================================================ + /*! + * \brief Move medium nodes of merged quadratic pyramids + */ + //================================================================================ + + void UpdateQuadraticPyramids(const set& commonApex, + SMESHDS_Mesh* meshDS) + { + typedef SMDS_StdIterator< const SMDS_MeshElement*, SMDS_ElemIteratorPtr > TStdElemIterator; + TStdElemIterator itEnd; + + // shift of node index to get medium nodes between the 4 base nodes and the apex + const int base2MediumShift = 9; + + set::const_iterator nIt = commonApex.begin(); + for ( ; nIt != commonApex.end(); ++nIt ) + { + SMESH_TNodeXYZ apex( *nIt ); + + vector< const SMDS_MeshElement* > pyrams // pyramids sharing the apex node + ( TStdElemIterator( apex._node->GetInverseElementIterator( SMDSAbs_Volume )), itEnd ); + + // Select medium nodes to keep and medium nodes to remove + + typedef map < const SMDS_MeshNode*, const SMDS_MeshNode*, TIDCompare > TN2NMap; + TN2NMap base2medium; // to keep + vector< const SMDS_MeshNode* > nodesToRemove; + + for ( unsigned i = 0; i < pyrams.size(); ++i ) + for ( int baseIndex = 0; baseIndex < PYRAM_APEX; ++baseIndex ) + { + SMESH_TNodeXYZ base = pyrams[i]->GetNode( baseIndex ); + const SMDS_MeshNode* medium = pyrams[i]->GetNode( baseIndex + base2MediumShift ); + TN2NMap::iterator b2m = base2medium.insert( make_pair( base._node, medium )).first; + if ( b2m->second != medium ) + { + nodesToRemove.push_back( medium ); + } + else + { + // move the kept medium node + gp_XYZ newXYZ = 0.5 * ( apex + base ); + meshDS->MoveNode( medium, newXYZ.X(), newXYZ.Y(), newXYZ.Z() ); + } + } + + // Within pyramids, replace nodes to remove by nodes to keep + + for ( unsigned i = 0; i < pyrams.size(); ++i ) + { + vector< const SMDS_MeshNode* > nodes( pyrams[i]->begin_nodes(), + pyrams[i]->end_nodes() ); + for ( int baseIndex = 0; baseIndex < PYRAM_APEX; ++baseIndex ) + { + const SMDS_MeshNode* base = pyrams[i]->GetNode( baseIndex ); + nodes[ baseIndex + base2MediumShift ] = base2medium[ base ]; + } + meshDS->ChangeElementNodes( pyrams[i], &nodes[0], nodes.size()); + } + + // Remove the replaced nodes - // sdt-like iterator used to get coordinates of nodes of mesh element -typedef SMDS_StdIterator< SMESH_MeshEditor::TNodeXYZ, SMDS_ElemIteratorPtr > TXyzIterator; + if ( !nodesToRemove.empty() ) + { + SMESHDS_SubMesh * sm = meshDS->MeshElements( nodesToRemove[0]->getshapeId() ); + for ( unsigned i = 0; i < nodesToRemove.size(); ++i ) + meshDS->RemoveFreeNode( nodesToRemove[i], sm, /*fromGroups=*/false); + } + } + } + +} //================================================================================ /*! - * \brief Destructor + * \brief Merge the two pyramids (i.e. fuse their apex) and others already merged with them */ //================================================================================ -StdMeshers_QuadToTriaAdaptor::~StdMeshers_QuadToTriaAdaptor() +void StdMeshers_QuadToTriaAdaptor::MergePiramids( const SMDS_MeshElement* PrmI, + const SMDS_MeshElement* PrmJ, + set & nodesToMove) { - // delete temporary faces - TQuad2Trias::iterator f_f = myResMap.begin(), ffEnd = myResMap.end(); - for ( ; f_f != ffEnd; ++f_f ) + const SMDS_MeshNode* Nrem = PrmJ->GetNode(4); // node to remove + //int nbJ = Nrem->NbInverseElements( SMDSAbs_Volume ); + SMESH_TNodeXYZ Pj( Nrem ); + + // an apex node to make common to all merged pyramids + SMDS_MeshNode* CommonNode = const_cast(PrmI->GetNode(4)); + if ( CommonNode == Nrem ) return; // already merged + //int nbI = CommonNode->NbInverseElements( SMDSAbs_Volume ); + SMESH_TNodeXYZ Pi( CommonNode ); + gp_XYZ Pnew = /*( nbI*Pi + nbJ*Pj ) / (nbI+nbJ);*/ 0.5 * ( Pi + Pj ); + CommonNode->setXYZ( Pnew.X(), Pnew.Y(), Pnew.Z() ); + + nodesToMove.insert( CommonNode ); + nodesToMove.erase ( Nrem ); + + typedef SMDS_StdIterator< const SMDS_MeshElement*, SMDS_ElemIteratorPtr > TStdElemIterator; + TStdElemIterator itEnd; + + // find and remove coincided faces of merged pyramids + vector< const SMDS_MeshElement* > inverseElems + // copy inverse elements to avoid iteration on changing container + ( TStdElemIterator( CommonNode->GetInverseElementIterator(SMDSAbs_Face)), itEnd); + for ( unsigned i = 0; i < inverseElems.size(); ++i ) { - TTriaList& fList = f_f->second; - TTriaList::iterator f = fList.begin(), fEnd = fList.end(); - for ( ; f != fEnd; ++f ) - delete *f; + const SMDS_MeshElement* FI = inverseElems[i]; + const SMDS_MeshElement* FJEqual = 0; + SMDS_ElemIteratorPtr triItJ = Nrem->GetInverseElementIterator(SMDSAbs_Face); + while ( !FJEqual && triItJ->more() ) + { + const SMDS_MeshElement* FJ = triItJ->next(); + if ( EqualTriangles( FJ, FI )) + FJEqual = FJ; + } + if ( FJEqual ) + { + removeTmpElement( FI ); + removeTmpElement( FJEqual ); + myRemovedTrias.insert( FI ); + myRemovedTrias.insert( FJEqual ); + } } - myResMap.clear(); -// TF2PyramMap::iterator itp = myPyram2Trias.begin(); -// for(; itp!=myPyram2Trias.end(); itp++) -// cout << itp->second << endl; + // set the common apex node to pyramids and triangles merged with J + inverseElems.assign( TStdElemIterator( Nrem->GetInverseElementIterator()), itEnd ); + for ( unsigned i = 0; i < inverseElems.size(); ++i ) + { + const SMDS_MeshElement* elem = inverseElems[i]; + vector< const SMDS_MeshNode* > nodes( elem->begin_nodes(), elem->end_nodes() ); + nodes[ elem->GetType() == SMDSAbs_Volume ? PYRAM_APEX : TRIA_APEX ] = CommonNode; + GetMeshDS()->ChangeElementNodes( elem, &nodes[0], nodes.size()); + } + ASSERT( Nrem->NbInverseElements() == 0 ); + GetMeshDS()->RemoveFreeNode( Nrem, + GetMeshDS()->MeshElements( Nrem->getshapeId()), + /*fromGroups=*/false); } +//================================================================================ +/*! + * \brief Merges adjacent pyramids + */ +//================================================================================ + +void StdMeshers_QuadToTriaAdaptor::MergeAdjacent(const SMDS_MeshElement* PrmI, + set& nodesToMove) +{ + TIDSortedElemSet adjacentPyrams; + bool mergedPyrams = false; + for(int k=0; k<4; k++) // loop on 4 base nodes of PrmI + { + const SMDS_MeshNode* n = PrmI->GetNode(k); + SMDS_ElemIteratorPtr vIt = n->GetInverseElementIterator( SMDSAbs_Volume ); + while ( vIt->more() ) + { + const SMDS_MeshElement* PrmJ = vIt->next(); + if ( PrmJ->NbCornerNodes() != 5 || !adjacentPyrams.insert( PrmJ ).second ) + continue; + if ( PrmI != PrmJ && TooCloseAdjacent( PrmI, PrmJ, GetMesh()->HasShapeToMesh() )) + { + MergePiramids( PrmI, PrmJ, nodesToMove ); + mergedPyrams = true; + // container of inverse elements can change + vIt = n->GetInverseElementIterator( SMDSAbs_Volume ); + } + } + } + if ( mergedPyrams ) + { + TIDSortedElemSet::iterator prm; + for (prm = adjacentPyrams.begin(); prm != adjacentPyrams.end(); ++prm) + MergeAdjacent( *prm, nodesToMove ); + } +} + +//================================================================================ +/*! + * \brief Constructor + */ +//================================================================================ + +StdMeshers_QuadToTriaAdaptor::StdMeshers_QuadToTriaAdaptor(): + myElemSearcher(0) +{ +} + +//================================================================================ +/*! + * \brief Destructor + */ +//================================================================================ + +StdMeshers_QuadToTriaAdaptor::~StdMeshers_QuadToTriaAdaptor() +{ + // temporary faces are deleted by ~SMESH_ProxyMesh() + if ( myElemSearcher ) delete myElemSearcher; + myElemSearcher=0; +} //======================================================================= //function : FindBestPoint @@ -78,30 +389,35 @@ StdMeshers_QuadToTriaAdaptor::~StdMeshers_QuadToTriaAdaptor() // (P, P1, P2) to be equilateral as much as possible // 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); + gp_Pnt Pbest = PC; + const double a = P1.Distance(P2); + const double b = P1.Distance(PC); + const double c = P2.Distance(PC); if( a < (b+c)/2 ) - return PC; + return Pbest; else { // find shift along V in order a to 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.XYZ() + aDir.XYZ() * shift; - return Pbest; + const double Vsize = V.Magnitude(); + if ( fabs( Vsize ) > std::numeric_limits::min() ) + { + const double shift = sqrt( a*a + (b*b-c*c)*(b*b-c*c)/16/a/a - (b*b+c*c)/2 ); + Pbest.ChangeCoord() += shift * V.XYZ() / Vsize; + } } + 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) { @@ -120,7 +436,7 @@ static bool HasIntersection3(const gp_Pnt& P, const gp_Pnt& PC, gp_Pnt& Pint, return false; if( IAICQ.NbPoints() == 1 ) { gp_Pnt PIn = IAICQ.Point(1); - double preci = 1.e-6; + const double preci = 1.e-10 * P.Distance(PC); // check if this point is internal for segment [PC,P] bool IsExternal = ( (PC.X()-PIn.X())*(P.X()-PIn.X()) > preci ) || @@ -133,32 +449,34 @@ static bool HasIntersection3(const gp_Pnt& P, const gp_Pnt& PC, gp_Pnt& Pint, gp_Vec V1(PIn,P1); gp_Vec V2(PIn,P2); gp_Vec V3(PIn,P3); - if( V1.Magnitude()(myElemSearcher); + + //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(); - int nbN = face->NbNodes(); - if( face->IsQuadratic() ) - nbN /= 2; - for ( int i = 0; i < nbN; ++i ) { - const SMDS_MeshNode* node = static_cast( nodeIt->next() ); - aContour->Append(gp_Pnt(node->X(), node->Y(), node->Z())); - } - if( HasIntersection(P, PC, Pres, aContour) ) { - res = true; - double tmp = PC.Distance(Pres); - if(tmp suspectElems; + searcher->GetElementsNearLine( line, SMDSAbs_Face, suspectElems); + + for ( int i = 0; i < suspectElems.size(); ++i ) + { + const SMDS_MeshElement* face = suspectElems[i]; + if ( face == NotCheckedFace ) continue; + Handle(TColgp_HSequenceOfPnt) aContour = new TColgp_HSequenceOfPnt; + for ( int i = 0; i < face->NbCornerNodes(); ++i ) + aContour->Append( SMESH_TNodeXYZ( face->GetNode(i) )); + if( HasIntersection(P, PC, Pres, aContour) ) { + res = true; + double tmp = PC.Distance(Pres); + if(tmpGetNode(1)==F2->GetNode(2) && F1->GetNode(2)==F2->GetNode(1) ) || - ( F1->GetNode(1)==F2->GetNode(1) && F1->GetNode(2)==F2->GetNode(2) ); -} - //================================================================================ /*! * \brief Prepare data for the given face @@ -288,20 +597,20 @@ int StdMeshers_QuadToTriaAdaptor::Preparation(const SMDS_MeshElement* face gp_Vec& VNorm, const SMDS_MeshElement** volumes) { - if( face->NbNodes() != ( face->IsQuadratic() ? 8 : 4 )) - if( face->NbNodes() != 4 ) - return NOT_QUAD; + if( face->NbCornerNodes() != 4 ) + { + return NOT_QUAD; + } int i = 0; gp_XYZ xyzC(0., 0., 0.); for ( i = 0; i < 4; ++i ) { - gp_XYZ p = SMESH_MeshEditor::TNodeXYZ( FNodes[i] = face->GetNode(i) ); + gp_XYZ p = SMESH_TNodeXYZ( FNodes[i] = face->GetNode(i) ); PN->SetValue( i+1, p ); xyzC += p; } PC = xyzC/4; - //cout<<" PC("< myPyramids; SMESHDS_Mesh * meshDS = aMesh.GetMeshDS(); SMESH_MesherHelper helper(aMesh); helper.IsQuadraticSubMesh(aShape); helper.SetElementsOnShape( true ); - for (TopExp_Explorer exp(aShape,TopAbs_FACE);exp.More();exp.Next()) { + if ( myElemSearcher ) delete myElemSearcher; + if ( aProxyMesh ) + myElemSearcher = SMESH_MeshEditor(&aMesh).GetElementSearcher( aProxyMesh->GetFaces(aShape)); + else + myElemSearcher = SMESH_MeshEditor(&aMesh).GetElementSearcher(); + + const SMESHDS_SubMesh * aSubMeshDSFace; + Handle(TColgp_HArray1OfPnt) PN = new TColgp_HArray1OfPnt(1,5); + Handle(TColgp_HArray1OfVec) VN = new TColgp_HArray1OfVec(1,4); + vector FNodes(5); + gp_Pnt PC; + gp_Vec VNorm; + + for (TopExp_Explorer exp(aShape,TopAbs_FACE);exp.More();exp.Next()) + { const TopoDS_Shape& aShapeFace = exp.Current(); - const SMESHDS_SubMesh * aSubMeshDSFace = meshDS->MeshElements( aShapeFace ); - if ( aSubMeshDSFace ) { - bool isRev = SMESH_Algo::IsReversedSubMesh( TopoDS::Face(aShapeFace), meshDS ); + if ( aProxyMesh ) + aSubMeshDSFace = aProxyMesh->GetSubMesh( aShapeFace ); + else + aSubMeshDSFace = meshDS->MeshElements( aShapeFace ); + + vector trias, quads; + bool hasNewTrias = false; + + if ( aSubMeshDSFace ) + { + bool isRev = false; + if ( helper.NbAncestors( aShapeFace, aMesh, aShape.ShapeType() ) > 1 ) + isRev = helper.IsReversedSubMesh( TopoDS::Face(aShapeFace) ); SMDS_ElemIteratorPtr iteratorElem = aSubMeshDSFace->GetElements(); - while ( iteratorElem->more() ) { // loop on elements on a face + while ( iteratorElem->more() ) // loop on elements on a geometrical 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; + // preparation step to get face info + int stat = Preparation(face, PN, VN, FNodes, PC, VNorm); + switch ( stat ) + { + case NOT_QUAD: - if(stat==2) { - // degenerate face - // add triangles to result map - 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] ); - TTriaList aList( 1, NewFace ); - myResMap.insert(make_pair(face,aList)); - continue; - } + trias.push_back( face ); + break; - 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 = PC.XYZ() + VNorm.XYZ() * 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("<AddFace( FNodes[0], FNodes[1], FNodes[2] ); + else + NewFace = meshDS->AddFace( FNodes[0], FNodes[2], FNodes[1] ); + storeTmpElement( NewFace ); + trias.push_back ( NewFace ); + quads.push_back( face ); + hasNewTrias = true; + break; } - else { - gp_Vec VB(PC,PCbest); - gp_Pnt PCbestTmp = PC.XYZ() + VB.XYZ() * 3.0; - bool check = CheckIntersection(PCbestTmp, PC, Pint, aMesh, aShape, aShapeFace); - if(check) { - double dist = PC.Distance(Pint)/3.; - if(distValue(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 = PC.XYZ() + VNorm.XYZ() * 0.001; + } + else { + // check possible intersection with other faces + gp_Pnt Pint; + bool check = CheckIntersection(PCbest, PC, Pint, aMesh, aShape, face); + if(check) { + //cout<<"--PC("<second; - for(i=0; i<4; i++) - triaList.push_back( new SMDS_FaceOfNodes( NewNode, FNodes[i], FNodes[i+1] )); - - // create pyramid - if ( isRev ) swap( FNodes[1], FNodes[3]); - SMDS_MeshVolume* aPyram = - helper.AddVolume( FNodes[0], FNodes[1], FNodes[2], FNodes[3], NewNode ); - myPyram2Trias.insert(make_pair(aPyram, & triaList)); + // create node for PCbest + SMDS_MeshNode* NewNode = helper.AddNode( PCbest.X(), PCbest.Y(), PCbest.Z() ); + + // add triangles to result map + for(i=0; i<4; i++) + { + trias.push_back ( meshDS->AddFace( NewNode, FNodes[i], FNodes[i+1] )); + storeTmpElement( trias.back() ); + } + // create a pyramid + if ( isRev ) swap( FNodes[1], FNodes[3]); + SMDS_MeshVolume* aPyram = + helper.AddVolume( FNodes[0], FNodes[1], FNodes[2], FNodes[3], NewNode ); + myPyramids.push_back(aPyram); + + quads.push_back( face ); + hasNewTrias = true; + break; + + } // case QUAD: + + } // switch ( stat ) } // end loop on elements on a face submesh + + bool sourceSubMeshIsProxy = false; + if ( aProxyMesh ) + { + // move proxy sub-mesh from other proxy mesh to this + sourceSubMeshIsProxy = takeProxySubMesh( aShapeFace, aProxyMesh ); + // move also tmp elements added in mesh + takeTmpElemsInMesh( aProxyMesh ); + } + if ( hasNewTrias ) + { + SMESH_ProxyMesh::SubMesh* prxSubMesh = getProxySubMesh( aShapeFace ); + prxSubMesh->ChangeElements( trias.begin(), trias.end() ); + + // delete tmp quadrangles removed from aProxyMesh + if ( sourceSubMeshIsProxy ) + { + for ( unsigned i = 0; i < quads.size(); ++i ) + removeTmpElement( quads[i] ); + + delete myElemSearcher; + myElemSearcher = + SMESH_MeshEditor(&aMesh).GetElementSearcher( aProxyMesh->GetFaces(aShape)); + } + } } } // end for(TopExp_Explorer exp(aShape,TopAbs_FACE);exp.More();exp.Next()) { - return Compute2ndPart(aMesh); + return Compute2ndPart(aMesh, myPyramids); } - -//======================================================================= -//function : Compute -//purpose : -//======================================================================= +//================================================================================ +/*! + * \brief Computes pyramids in mesh with no shape + */ +//================================================================================ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh) { - myResMap.clear(); - myPyram2Trias.clear(); + SMESH_ProxyMesh::setMesh( aMesh ); + SMESH_ProxyMesh::_allowedTypes.push_back( SMDSEntity_Triangle ); + SMESH_ProxyMesh::_allowedTypes.push_back( SMDSEntity_Quad_Triangle ); + if ( aMesh.NbQuadrangles() < 1 ) + return false; + + // find if there is a group of faces identified as skin faces, with normal going outside the volume + std::string groupName = "skinFaces"; + SMESHDS_GroupBase* groupDS = 0; + SMESH_Mesh::GroupIteratorPtr groupIt = aMesh.GetGroups(); + while ( groupIt->more() ) + { + groupDS = 0; + SMESH_Group * group = groupIt->next(); + if ( !group ) continue; + groupDS = group->GetGroupDS(); + if ( !groupDS || groupDS->IsEmpty() ) + { + groupDS = 0; + continue; + } + if (groupDS->GetType() != SMDSAbs_Face) + { + groupDS = 0; + continue; + } + std::string grpName = group->GetName(); + if (grpName == groupName) + { + MESSAGE("group skinFaces provided"); + break; + } + else + groupDS = 0; + } + + vector myPyramids; SMESH_MesherHelper helper(aMesh); helper.IsQuadraticSubMesh(aMesh.GetShapeToMesh()); helper.SetElementsOnShape( true ); - SMESHDS_Mesh * meshDS = aMesh.GetMeshDS(); + if ( !myElemSearcher ) + myElemSearcher = SMESH_MeshEditor(&aMesh).GetElementSearcher(); + SMESH_ElementSearcher* searcher = const_cast(myElemSearcher); - SMDS_FaceIteratorPtr fIt = meshDS->facesIterator(); - TIDSortedElemSet sortedFaces; // 0020279: control the "random" use when using mesh algorithms - while( fIt->more()) sortedFaces.insert( fIt->next() ); + SMESHDS_Mesh * meshDS = aMesh.GetMeshDS(); + SMESH_ProxyMesh::SubMesh* prxSubMesh = getProxySubMesh(); - TIDSortedElemSet::iterator itFace = sortedFaces.begin(), fEnd = sortedFaces.end(); - for ( ; itFace != fEnd; ++itFace ) + SMDS_FaceIteratorPtr fIt = meshDS->facesIterator(/*idInceasingOrder=*/true); + while( fIt->more()) { - const SMDS_MeshElement* face = *itFace; + const SMDS_MeshElement* face = fIt->next(); if ( !face ) continue; - //cout<GetID() = "<GetID()<Value(1).Distance(PN->Value(2)) + PN->Value(2).Distance(PN->Value(3)); // far points in VNorm direction gp_Pnt Ptmp1 = PC.XYZ() + VNorm.XYZ() * tmp * 1.e6; @@ -566,12 +980,17 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh) double dist1 = RealLast(); double dist2 = RealLast(); gp_Pnt Pres1,Pres2; - for (TIDSortedElemSet::iterator itF = sortedFaces.begin(); itF != fEnd; ++itF ) { - const SMDS_MeshElement* F = *itF; + + gp_Ax1 line( PC, VNorm ); + vector< const SMDS_MeshElement* > suspectElems; + searcher->GetElementsNearLine( line, SMDSAbs_Face, suspectElems); + + for ( int iF = 0; iF < suspectElems.size(); ++iF ) { + const SMDS_MeshElement* F = suspectElems[iF]; if(F==face) continue; Handle(TColgp_HSequenceOfPnt) aContour = new TColgp_HSequenceOfPnt; for ( int i = 0; i < 4; ++i ) - aContour->Append( SMESH_MeshEditor::TNodeXYZ( F->GetNode(i) )); + aContour->Append( SMESH_TNodeXYZ( F->GetNode(i) )); gp_Pnt PPP; if( !volumes[0] && HasIntersection(Ptmp1, PC, PPP, aContour) ) { IsOK1 = true; @@ -610,14 +1029,16 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh) } } if(!IsRev) - NewFace = new SMDS_FaceOfNodes( FNodes[0], FNodes[1], FNodes[2] ); + NewFace = meshDS->AddFace( 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)); + NewFace = meshDS->AddFace( FNodes[0], FNodes[2], FNodes[1] ); + storeTmpElement( NewFace ); + prxSubMesh->AddElement( NewFace ); continue; } + // Case of non-degenerated quadrangle + // Find pyramid peak gp_XYZ PCbest(0., 0., 0.); // pyramid peak @@ -629,12 +1050,13 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh) PCbest /= 4; double height = PC.Distance(PCbest); // pyramid height to precise - if(height<1.e-6) { + if ( height < 1.e-6 ) { // create new PCbest using a bit shift along VNorm PCbest = PC.XYZ() + VNorm.XYZ() * 0.001; height = PC.Distance(PCbest); + if ( height < std::numeric_limits::min() ) + return false; // batterfly element } - //cout<<" PCbest("< suspectElems; + searcher->GetElementsNearLine( line, SMDSAbs_Face, suspectElems); + + for ( int iF = 0; iF < suspectElems.size(); ++iF ) { - const SMDS_MeshElement* F = *itF; + const SMDS_MeshElement* F = suspectElems[iF]; if(F==face) continue; Handle(TColgp_HSequenceOfPnt) aContour = new TColgp_HSequenceOfPnt; int nbN = F->NbNodes() / ( F->IsQuadratic() ? 2 : 1 ); for ( i = 0; i < nbN; ++i ) - aContour->Append( SMESH_MeshEditor::TNodeXYZ( F->GetNode(i) )); + aContour->Append( SMESH_TNodeXYZ( F->GetNode(i) )); gp_Pnt intP; for ( int isRev = 0; isRev < 2; ++isRev ) { @@ -669,6 +1096,10 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh) } } + // if the face belong to the group of skinFaces, do not build a pyramid outside + if (groupDS && groupDS->Contains(face)) + intersected[0] = false; + // Create one or two pyramids for ( int isRev = 0; isRev < 2; ++isRev ) @@ -681,14 +1112,14 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh) SMDS_MeshNode* NewNode = helper.AddNode( PCbest.X(), PCbest.Y(), PCbest.Z() ); // add triangles to result map - TTriaList& aList = myResMap.insert( make_pair( face, TTriaList()))->second; for(i=0; i<4; i++) { - SMDS_FaceOfNodes* NewFace; + SMDS_MeshFace* NewFace; if(isRev) - NewFace = new SMDS_FaceOfNodes( NewNode, FNodes[i], FNodes[i+1] ); + NewFace = meshDS->AddFace( NewNode, FNodes[i], FNodes[i+1] ); else - NewFace = new SMDS_FaceOfNodes( NewNode, FNodes[i+1], FNodes[i] ); - aList.push_back(NewFace); + NewFace = meshDS->AddFace( NewNode, FNodes[i+1], FNodes[i] ); + storeTmpElement( NewFace ); + prxSubMesh->AddElement( NewFace ); } // create a pyramid SMDS_MeshVolume* aPyram; @@ -696,259 +1127,201 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh) aPyram = helper.AddVolume( FNodes[0], FNodes[1], FNodes[2], FNodes[3], NewNode ); else aPyram = helper.AddVolume( FNodes[0], FNodes[3], FNodes[2], FNodes[1], NewNode ); - myPyram2Trias.insert(make_pair(aPyram, & aList)); + myPyramids.push_back(aPyram); } } // end loop on all faces - return Compute2ndPart(aMesh); + return Compute2ndPart(aMesh, myPyramids); } -//======================================================================= -//function : Compute2ndPart -//purpose : Update created pyramids and faces to avoid their intersection -//======================================================================= +//================================================================================ +/*! + * \brief Update created pyramids and faces to avoid their intersection + */ +//================================================================================ -bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh& aMesh) +bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh& aMesh, + const vector& myPyramids) { + if(myPyramids.empty()) + return true; + SMESHDS_Mesh * meshDS = aMesh.GetMeshDS(); + int i, j, k, myShapeID = myPyramids[0]->GetNode(4)->getshapeId(); - // check intersections between created pyramids + if ( myElemSearcher ) delete myElemSearcher; + myElemSearcher = SMESH_MeshEditor(&aMesh).GetElementSearcher(); + SMESH_ElementSearcher* searcher = const_cast(myElemSearcher); - if(myPyram2Trias.empty()) - return true; + set nodesToMove; - int k = 0; + // check adjacent pyramids - // for each pyramid store list of merged pyramids with their faces - typedef map< const SMDS_MeshElement*, list< TPyram2Trias::iterator > > TPyram2Merged; - TPyram2Merged MergesInfo; + for ( i = 0; i < myPyramids.size(); ++i ) + { + const SMDS_MeshElement* PrmI = myPyramids[i]; + MergeAdjacent( PrmI, nodesToMove ); + } // iterate on all pyramids - TPyram2Trias::iterator itPi = myPyram2Trias.begin(), itPEnd = myPyram2Trias.end(); - for ( ; itPi != itPEnd; ++itPi ) + for ( i = 0; i < myPyramids.size(); ++i ) { - const SMDS_MeshElement* PrmI = itPi->first; - TPyram2Merged::iterator pMergesI = MergesInfo.find( PrmI ); - - TXyzIterator xyzIt( PrmI->nodesIterator() ); - vector PsI( xyzIt, TXyzIterator() ); + const SMDS_MeshElement* PrmI = myPyramids[i]; // compare PrmI with all the rest pyramids - bool NeedMove = false; - TPyram2Trias::iterator itPj = itPi; - for ( ++itPj; itPj != itPEnd; ++itPj ) + + // collect adjacent pyramids and nodes coordinates of PrmI + set checkedPyrams; + vector PsI(5); + for(k=0; k<5; k++) // loop on 4 base nodes of PrmI + { + const SMDS_MeshNode* n = PrmI->GetNode(k); + PsI[k] = SMESH_TNodeXYZ( n ); + SMDS_ElemIteratorPtr vIt = n->GetInverseElementIterator( SMDSAbs_Volume ); + while ( vIt->more() ) + { + const SMDS_MeshElement* PrmJ = vIt->next(); + if ( SMESH_Algo::GetCommonNodes( PrmI, PrmJ ).size() > 1 ) + checkedPyrams.insert( PrmJ ); + } + } + + // check intersection with distant pyramids + for(k=0; k<4; k++) // loop on 4 base nodes of PrmI { - const SMDS_MeshElement* PrmJ = itPj->first; - TPyram2Merged::iterator pMergesJ = MergesInfo.find( PrmJ ); - - // check if two pyramids already merged - if ( pMergesJ != MergesInfo.end() && - find(pMergesJ->second.begin(),pMergesJ->second.end(), itPi )!=pMergesJ->second.end()) - continue; // already merged - - xyzIt = TXyzIterator( PrmJ->nodesIterator() ); - vector PsJ( xyzIt, TXyzIterator() ); - - bool hasInt = false; - gp_Pnt Pint; - for(k=0; k<4 && !hasInt; k++) { - gp_Vec Vtmp(PsI[k],PsI[4]); - gp_Pnt Pshift = PsI[k].XYZ() + Vtmp.XYZ() * 0.01; - hasInt = + gp_Vec Vtmp(PsI[k],PsI[4]); + gp_Ax1 line( PsI[k], Vtmp ); + vector< const SMDS_MeshElement* > suspectPyrams; + searcher->GetElementsNearLine( line, SMDSAbs_Volume, suspectPyrams); + + for ( j = 0; j < suspectPyrams.size(); ++j ) + { + const SMDS_MeshElement* PrmJ = suspectPyrams[j]; + if ( PrmJ == PrmI || PrmJ->NbCornerNodes() != 5 ) + continue; + if ( myShapeID != PrmJ->GetNode(4)->getshapeId()) + continue; // pyramid from other SOLID + if ( PrmI->GetNode(4) == PrmJ->GetNode(4) ) + continue; // pyramids PrmI and PrmJ already merged + if ( !checkedPyrams.insert( PrmJ ).second ) + continue; // already checked + + TXyzIterator xyzIt( PrmJ->nodesIterator() ); + vector PsJ( xyzIt, TXyzIterator() ); + + gp_Pnt Pint; + bool hasInt=false; + for(k=0; k<4 && !hasInt; k++) { + gp_Vec Vtmp(PsI[k],PsI[4]); + gp_Pnt Pshift = PsI[k].XYZ() + Vtmp.XYZ() * 0.01; // base node moved a bit to apex + hasInt = ( HasIntersection3( Pshift, PsI[4], Pint, PsJ[0], PsJ[1], PsJ[4]) || HasIntersection3( Pshift, PsI[4], Pint, PsJ[1], PsJ[2], PsJ[4]) || HasIntersection3( Pshift, PsI[4], Pint, PsJ[2], PsJ[3], PsJ[4]) || HasIntersection3( Pshift, PsI[4], Pint, PsJ[3], PsJ[0], PsJ[4]) ); - } - for(k=0; k<4 && !hasInt; k++) { - gp_Vec Vtmp(PsJ[k],PsJ[4]); - gp_Pnt Pshift = PsJ[k].XYZ() + Vtmp.XYZ() * 0.01; - hasInt = - ( HasIntersection3( Pshift, PsJ[4], Pint, PsI[0], PsI[1], PsI[4]) || - HasIntersection3( Pshift, PsJ[4], Pint, PsI[1], PsI[2], PsI[4]) || - HasIntersection3( Pshift, PsJ[4], Pint, PsI[2], PsI[3], PsI[4]) || - HasIntersection3( Pshift, PsJ[4], Pint, PsI[3], PsI[0], PsI[4]) ); - } - if(hasInt) { - // count common nodes of base faces of two pyramids - int nbc = 0; - for(k=0; k<4; k++) - nbc += int ( PrmI->GetNodeIndex( PrmJ->GetNode(k) ) >= 0 ); - //cout<<" nbc = "<0) + if ( hasInt ) { - // Merge the two pyramids and others already merged with them + // count common nodes of base faces of two pyramids + int nbc = 0; + for (k=0; k<4; k++) + nbc += int ( PrmI->GetNodeIndex( PrmJ->GetNode(k) ) >= 0 ); - // initialize merge info of pyramids - if ( pMergesI == MergesInfo.end() ) // first merge of PrmI - { - pMergesI = MergesInfo.insert( make_pair( PrmI, list())).first; - pMergesI->second.push_back( itPi ); - } - if ( pMergesJ == MergesInfo.end() ) // first merge of PrmJ - { - pMergesJ = MergesInfo.insert( make_pair( PrmJ, list())).first; - pMergesJ->second.push_back( itPj ); - } - int nbI = pMergesI->second.size(), nbJ = pMergesJ->second.size(); - - // an apex node to make common to all merged pyramids - SMDS_MeshNode* CommonNode = const_cast(PrmI->GetNode(4)); - CommonNode->setXYZ( ( nbI*PsI[4].X() + nbJ*PsJ[4].X() ) / (nbI+nbJ), - ( nbI*PsI[4].Y() + nbJ*PsJ[4].Y() ) / (nbI+nbJ), - ( nbI*PsI[4].Z() + nbJ*PsJ[4].Z() ) / (nbI+nbJ) ); - NeedMove = true; - const SMDS_MeshNode* Nrem = PrmJ->GetNode(4); // node to remove - - list< TPyram2Trias::iterator >& aMergesI = pMergesI->second; - list< TPyram2Trias::iterator >& aMergesJ = pMergesJ->second; - - // find and remove coincided faces of merged pyramids - list< TPyram2Trias::iterator >::iterator itPttI, itPttJ; - TTriaList::iterator trI, trJ; - for ( itPttI = aMergesI.begin(); itPttI != aMergesI.end(); ++itPttI ) - { - TTriaList* triaListI = (*itPttI)->second; - for ( trI = triaListI->begin(); trI != triaListI->end(); ) - { - const SMDS_FaceOfNodes* FI = *trI; - - for ( itPttJ = aMergesJ.begin(); itPttJ != aMergesJ.end() && FI; ++itPttJ ) - { - TTriaList* triaListJ = (*itPttJ)->second; - for ( trJ = triaListJ->begin(); trJ != triaListJ->end(); ) - { - const SMDS_FaceOfNodes* FJ = *trJ; - - if( EqualTriangles(FI,FJ) ) - { - delete FI; - delete FJ; - FI = FJ = 0; - trI = triaListI->erase( trI ); - trJ = triaListJ->erase( trJ ); - break; // only one triangle of a pyramid can coincide with another pyramid - } - ++trJ; - } - } - if ( FI ) ++trI; // increament if triangle not deleted - } - } + if ( nbc == 4 ) + continue; // pyrams have a common base face - // set the common apex node to pyramids and triangles merged with J - for ( itPttJ = aMergesJ.begin(); itPttJ != aMergesJ.end(); ++itPttJ ) + if(nbc>0) { - const SMDS_MeshElement* Prm = (*itPttJ)->first; - TTriaList* triaList = (*itPttJ)->second; - - vector< const SMDS_MeshNode* > nodes( Prm->begin_nodes(), Prm->end_nodes() ); - nodes[4] = CommonNode; - meshDS->ChangeElementNodes( Prm, &nodes[0], nodes.size()); + // Merge the two pyramids and others already merged with them + MergePiramids( PrmI, PrmJ, nodesToMove ); + } + else { // nbc==0 - for ( TTriaList::iterator trIt = triaList->begin(); trIt != triaList->end(); ++trIt ) - { - SMDS_FaceOfNodes* Ftria = const_cast< SMDS_FaceOfNodes*>( *trIt ); - const SMDS_MeshNode* NF[3] = { CommonNode, Ftria->GetNode(1), Ftria->GetNode(2)}; - Ftria->ChangeNodes(NF, 3); + // decrease height of pyramids + gp_XYZ PCi(0,0,0), PCj(0,0,0); + for(k=0; k<4; k++) { + PCi += PsI[k].XYZ(); + PCj += PsJ[k].XYZ(); } + PCi /= 4; PCj /= 4; + gp_Vec VN1(PCi,PsI[4]); + gp_Vec VN2(PCj,PsJ[4]); + gp_Vec VI1(PCi,Pint); + gp_Vec VI2(PCj,Pint); + double ang1 = fabs(VN1.Angle(VI1)); + double ang2 = fabs(VN2.Angle(VI2)); + double coef1 = 0.5 - (( ang1 < M_PI/3. ) ? cos(ang1)*0.25 : 0 ); + double coef2 = 0.5 - (( ang2 < M_PI/3. ) ? cos(ang2)*0.25 : 0 ); // cos(ang2) ? +// double coef2 = 0.5; +// if(ang2(PrmI->GetNode(4)); + aNode1->setXYZ( PCi.X()+VN1.X(), PCi.Y()+VN1.Y(), PCi.Z()+VN1.Z() ); + SMDS_MeshNode* aNode2 = const_cast(PrmJ->GetNode(4)); + aNode2->setXYZ( PCj.X()+VN2.X(), PCj.Y()+VN2.Y(), PCj.Z()+VN2.Z() ); + nodesToMove.insert( aNode1 ); + nodesToMove.insert( aNode2 ); } + // fix intersections that could appear after apex movement + MergeAdjacent( PrmI, nodesToMove ); + MergeAdjacent( PrmJ, nodesToMove ); - // join MergesInfo of merged pyramids - for ( k = 0, itPttI = aMergesI.begin(); k < nbI; ++itPttI, ++k ) - { - const SMDS_MeshElement* PrmI = (*itPttI)->first; - list< TPyram2Trias::iterator >& merges = MergesInfo[ PrmI ]; - merges.insert( merges.end(), aMergesJ.begin(), aMergesJ.end() ); - } - for ( k = 0, itPttJ = aMergesJ.begin(); k < nbJ; ++itPttJ, ++k ) - { - const SMDS_MeshElement* PrmJ = (*itPttJ)->first; - list< TPyram2Trias::iterator >& merges = MergesInfo[ PrmJ ]; - merges.insert( merges.end(), aMergesI.begin(), aMergesI.end() ); - } + } // end if(hasInt) + } // loop on suspectPyrams + } // loop on 4 base nodes of PrmI - // removing node - meshDS->RemoveNode(Nrem); - } - else { // nbc==0 + } // loop on all pyramids - // decrease height of pyramids - gp_XYZ PC1(0,0,0), PC2(0,0,0); - for(k=0; k<4; k++) { - PC1 += PsI[k].XYZ(); - PC2 += PsJ[k].XYZ(); - } - PC1 /= 4; PC2 /= 4; - gp_Vec VN1(PC1,PsI[4]); - gp_Vec VI1(PC1,Pint); - gp_Vec VN2(PC2,PsJ[4]); - gp_Vec VI2(PC2,Pint); - double ang1 = fabs(VN1.Angle(VI1)); - double ang2 = fabs(VN2.Angle(VI2)); - double h1,h2; - if(ang1>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(PrmI->GetNode(4)); - aNode1->setXYZ( PC1.X()+VN1.X(), PC1.Y()+VN1.Y(), PC1.Z()+VN1.Z() ); - SMDS_MeshNode* aNode2 = const_cast(PrmJ->GetNode(4)); - aNode2->setXYZ( PC2.X()+VN2.X(), PC2.Y()+VN2.Y(), PC2.Z()+VN2.Z() ); - NeedMove = true; - } - } // end if(hasInt) - } - if( NeedMove && !meshDS->IsEmbeddedMode() ) - { - const SMDS_MeshNode* apex = PrmI->GetNode( 4 ); - meshDS->MoveNode( apex, apex->X(), apex->Y(), apex->Z() ); - } + if( !nodesToMove.empty() && !meshDS->IsEmbeddedMode() ) + { + set::iterator n = nodesToMove.begin(); + for ( ; n != nodesToMove.end(); ++n ) + meshDS->MoveNode( *n, (*n)->X(), (*n)->Y(), (*n)->Z() ); } - // rebind triangles of pyramids sharing the same base quadrangle to the first - // entrance of the base quadrangle - TQuad2Trias::iterator q2t = myResMap.begin(), q2tPrev = q2t; - for ( ++q2t; q2t != myResMap.end(); ++q2t, ++q2tPrev ) + // move medium nodes of merged quadratic pyramids + if ( myPyramids[0]->IsQuadratic() ) + UpdateQuadraticPyramids( nodesToMove, GetMeshDS() ); + + // erase removed triangles from the proxy mesh + if ( !myRemovedTrias.empty() ) { - if ( q2t->first == q2tPrev->first ) - q2tPrev->second.splice( q2tPrev->second.end(), q2t->second ); + for ( int i = 0; i <= meshDS->MaxShapeIndex(); ++i ) + if ( SMESH_ProxyMesh::SubMesh* sm = findProxySubMesh(i)) + { + vector faces; + faces.reserve( sm->NbElements() ); + SMDS_ElemIteratorPtr fIt = sm->GetElements(); + while ( fIt->more() ) + { + const SMDS_MeshElement* tria = fIt->next(); + set::iterator rmTria = myRemovedTrias.find( tria ); + if ( rmTria != myRemovedTrias.end() ) + myRemovedTrias.erase( rmTria ); + else + faces.push_back( tria ); + } + sm->ChangeElements( faces.begin(), faces.end() ); + } } - myPyram2Trias.clear(); // no more needed myDegNodes.clear(); - return true; -} + delete myElemSearcher; + myElemSearcher=0; -//================================================================================ -/*! - * \brief Return list of created triangles for given face - */ -//================================================================================ - -const list* StdMeshers_QuadToTriaAdaptor::GetTriangles (const SMDS_MeshElement* aQuad) -{ - TQuad2Trias::iterator it = myResMap.find(aQuad); - if( it != myResMap.end() ) { - return & it->second; - } - return 0; + return true; }