X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FStdMeshers%2FStdMeshers_QuadToTriaAdaptor.cxx;h=db32a5f1a5f92274ef99ed98d721e89459bd9711;hb=2c607013a23bd4e7ba07e72e0c04dee2c1209cff;hp=25b472c3eed3d427ee6152666f49627d9f60f2b7;hpb=8fa039a796957b302d86d90b22afc0a998573f83;p=modules%2Fsmesh.git diff --git a/src/StdMeshers/StdMeshers_QuadToTriaAdaptor.cxx b/src/StdMeshers/StdMeshers_QuadToTriaAdaptor.cxx index 25b472c3e..db32a5f1a 100644 --- a/src/StdMeshers/StdMeshers_QuadToTriaAdaptor.cxx +++ b/src/StdMeshers/StdMeshers_QuadToTriaAdaptor.cxx @@ -1,23 +1,22 @@ -// Copyright (C) 2007-2010 CEA/DEN, EDF R&D, OPEN CASCADE +// Copyright (C) 2007-2011 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 @@ -44,14 +43,13 @@ 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_MeshEditor::TNodeXYZ, SMDS_ElemIteratorPtr > TXyzIterator; +typedef SMDS_StdIterator< SMESH_TNodeXYZ, SMDS_ElemIteratorPtr > TXyzIterator; namespace { - //================================================================================ /*! * \brief Return true if two nodes of triangles are equal @@ -64,88 +62,6 @@ namespace ( 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 Merge the two pyramids (i.e. fuse their apex) and others already merged with them - */ - //================================================================================ - - void MergePyramids( const SMDS_MeshElement* PrmI, - const SMDS_MeshElement* PrmJ, - SMESHDS_Mesh* meshDS, - TRemTrias& tempTrias, - set & nodesToMove) - { - const SMDS_MeshNode* Nrem = PrmJ->GetNode(4); // node to remove - int nbJ = Nrem->NbInverseElements( SMDSAbs_Volume ); - SMESH_MeshEditor::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 ); - //cerr << __LINE__ << " " << nbI << " " << nbJ << endl; - SMESH_MeshEditor::TNodeXYZ Pi( CommonNode ); - gp_XYZ Pnew = ( nbI*Pi + nbJ*Pj ) / (nbI+nbJ); - CommonNode->setXYZ( Pnew.X(), Pnew.Y(), Pnew.Z() ); - - nodesToMove.insert( CommonNode ); - nodesToMove.erase ( Nrem ); - - // find and remove coincided faces of merged pyramids - SMDS_ElemIteratorPtr triItI = CommonNode->GetInverseElementIterator(SMDSAbs_Face); - while ( triItI->more() ) - { - const SMDS_MeshElement* FI = triItI->next(); - 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 ) - { - //meshDS->RemoveFreeElement(FI, 0, false); - //meshDS->RemoveFreeElement(FJEqual, 0, false); - tempTrias[FI] = false; - tempTrias[FJEqual] = false; - } - } - - // set the common apex node to pyramids and triangles merged with J - //cerr << __LINE__ << " NbInverseElements " << Nrem->NbInverseElements() << endl; - SMDS_ElemIteratorPtr itJ = Nrem->GetInverseElementIterator(); - while ( itJ->more() ) - { - const SMDS_MeshElement* elem = itJ->next(); - if ( elem->GetType() == SMDSAbs_Volume ) // pyramid - { - vector< const SMDS_MeshNode* > nodes( elem->begin_nodes(), elem->end_nodes() ); - //cerr << __LINE__ << " volId " << elem->GetID() << " nbNodes " << nodes.size() << endl; - nodes[4] = CommonNode; - MESSAGE("ChangeElementNodes"); - meshDS->ChangeElementNodes( elem, &nodes[0], nodes.size()); - } - else if ( tempTrias.count(elem) ) // tmp triangles - { - //cerr << __LINE__ << " triaId " << elem->GetID() << endl; - ((SMDS_VtkFace*) elem )->ChangeApex( CommonNode ); - } -// else -// { -// cerr << __LINE__ << " other " << elem->GetVtkType() << endl; -// } - } - //cerr << __LINE__ << " NbInverseElements " << Nrem->NbInverseElements() << endl; - ASSERT( Nrem->NbInverseElements() == 0 ); - meshDS->RemoveFreeNode( Nrem, - meshDS->MeshElements( Nrem->getshapeId()), - /*fromGroups=*/false); - } - //================================================================================ /*! * \brief Return true if two adjacent pyramids are too close one to another @@ -174,7 +90,7 @@ namespace int ind = baseNodes[0] ? 1:0; if ( baseNodes[ ind ]) return false; // pyramids with a common base face - baseNodes [ ind ] = PrmI->GetNode(i); + baseNodes [ ind ] = PrmI->GetNode(i); baseNodesIndI[ ind ] = i; baseNodesIndJ[ ind ] = j; } @@ -182,10 +98,10 @@ namespace if ( !baseNodes[1] ) return false; // not adjacent // Get normals of triangles sharing baseNodes - gp_XYZ apexI = SMESH_MeshEditor::TNodeXYZ( nApexI ); - gp_XYZ apexJ = SMESH_MeshEditor::TNodeXYZ( nApexJ ); - gp_XYZ base1 = SMESH_MeshEditor::TNodeXYZ( baseNodes[0]); - gp_XYZ base2 = SMESH_MeshEditor::TNodeXYZ( baseNodes[1]); + 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 ); @@ -197,12 +113,11 @@ namespace bool tooClose = ( angle < 15 * PI180 ); // Check if pyramids collide - bool isOutI, isOutJ; if ( !tooClose && baI * baJ > 0 ) { // find out if nI points outside of PrmI or inside int dInd = baseNodesIndI[1] - baseNodesIndI[0]; - isOutI = ( abs(dInd)==1 ) ? dInd < 0 : dInd > 0; + bool isOutI = ( abs(dInd)==1 ) ? dInd < 0 : dInd > 0; // find out sign of projection of nJ to baI double proj = baI * nJ; @@ -214,11 +129,17 @@ namespace if ( tooClose && !hasShape ) { // check order of baseNodes within pyramids, it must be opposite - int dInd = baseNodesIndJ[1] - baseNodesIndJ[0]; - isOutJ = ( abs(dInd)==1 ) ? dInd < 0 : dInd > 0; + 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; @@ -233,6 +154,9 @@ namespace 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 || @@ -240,8 +164,7 @@ namespace continue; // f is a base quadrangle // check projections of face direction (baOFN) to triange normals (nI and nJ) - gp_Vec baOFN( base1, SMESH_MeshEditor::TNodeXYZ( otherFaceNode )); - ( isOutI ? nJ : nI ).Reverse(); + gp_Vec baOFN( base1, SMESH_TNodeXYZ( otherFaceNode )); if ( nI * baOFN > 0 && nJ * baOFN > 0 ) { tooClose = false; // f is between pyramids @@ -255,42 +178,180 @@ namespace //================================================================================ /*! - * \brief Merges adjacent pyramids + * \brief Move medium nodes of merged quadratic pyramids */ //================================================================================ - void MergeAdjacent(const SMDS_MeshElement* PrmI, - SMESH_Mesh& mesh, - TRemTrias & tempTrias, - set& nodesToMove) + void UpdateQuadraticPyramids(const set& commonApex, + SMESHDS_Mesh* meshDS) { - TIDSortedElemSet adjacentPyrams, mergedPyrams; - for(int k=0; k<4; k++) // loop on 4 base nodes of PrmI + 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 ) { - const SMDS_MeshNode* n = PrmI->GetNode(k); - SMDS_ElemIteratorPtr vIt = n->GetInverseElementIterator( SMDSAbs_Volume ); - while ( vIt->more() ) + 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 ) { - const SMDS_MeshElement* PrmJ = vIt->next(); - if ( PrmJ->NbCornerNodes() != 5 || !adjacentPyrams.insert( PrmJ ).second ) - continue; - if ( PrmI != PrmJ && TooCloseAdjacent( PrmI, PrmJ, mesh.HasShapeToMesh() )) + vector< const SMDS_MeshNode* > nodes( pyrams[i]->begin_nodes(), + pyrams[i]->end_nodes() ); + for ( int baseIndex = 0; baseIndex < PYRAM_APEX; ++baseIndex ) { - MergePyramids( PrmI, PrmJ, mesh.GetMeshDS(), tempTrias, nodesToMove ); - mergedPyrams.insert( PrmJ ); + 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 + + 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); } } - if ( !mergedPyrams.empty() ) + } + +} + +//================================================================================ +/*! + * \brief Merge the two pyramids (i.e. fuse their apex) and others already merged with them + */ +//================================================================================ + +void StdMeshers_QuadToTriaAdaptor::MergePiramids( const SMDS_MeshElement* PrmI, + const SMDS_MeshElement* PrmJ, + set & nodesToMove) +{ + 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 ) + { + const SMDS_MeshElement* FI = inverseElems[i]; + const SMDS_MeshElement* FJEqual = 0; + SMDS_ElemIteratorPtr triItJ = Nrem->GetInverseElementIterator(SMDSAbs_Face); + while ( !FJEqual && triItJ->more() ) { - TIDSortedElemSet::iterator prm; -// for (prm = mergedPyrams.begin(); prm != mergedPyrams.end(); ++prm) -// MergeAdjacent( *prm, mesh, tempTrias, nodesToMove ); + const SMDS_MeshElement* FJ = triItJ->next(); + if ( EqualTriangles( FJ, FI )) + FJEqual = FJ; + } + if ( FJEqual ) + { + removeTmpElement( FI ); + removeTmpElement( FJEqual ); + myRemovedTrias.insert( FI ); + myRemovedTrias.insert( FJEqual ); + } + } + + // 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); +} - for (prm = adjacentPyrams.begin(); prm != adjacentPyrams.end(); ++prm) - MergeAdjacent( *prm, mesh, tempTrias, nodesToMove ); +//================================================================================ +/*! + * \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 ); + } } //================================================================================ @@ -300,7 +361,7 @@ namespace //================================================================================ StdMeshers_QuadToTriaAdaptor::StdMeshers_QuadToTriaAdaptor(): - myElemSearcher(0), myNbTriangles(0) + myElemSearcher(0) { } @@ -312,20 +373,7 @@ StdMeshers_QuadToTriaAdaptor::StdMeshers_QuadToTriaAdaptor(): StdMeshers_QuadToTriaAdaptor::~StdMeshers_QuadToTriaAdaptor() { - // delete temporary faces - TQuad2Trias::iterator f_f = myResMap.begin(), ffEnd = myResMap.end(); - for ( ; f_f != ffEnd; ++f_f ) - { - TTriaList& fList = f_f->second; - TTriaList::iterator f = fList.begin(), fEnd = fList.end(); - for ( ; f != fEnd; ++f ) - { - const SMDS_MeshElement *elem = *f; - SMDS_Mesh::_meshList[elem->getMeshId()]->RemoveFreeElement(elem); - } - } - myResMap.clear(); - + // temporary faces are deleted by ~SMESH_ProxyMesh() if ( myElemSearcher ) delete myElemSearcher; myElemSearcher=0; } @@ -499,22 +547,13 @@ bool StdMeshers_QuadToTriaAdaptor::CheckIntersection (const gp_Pnt& P, vector< const SMDS_MeshElement* > suspectElems; searcher->GetElementsNearLine( line, SMDSAbs_Face, suspectElems); -// for (TopExp_Explorer exp(aShape,TopAbs_FACE);exp.More();exp.Next()) { -// const TopoDS_Shape& aShapeFace = exp.Current(); -// if(aShapeFace==NotCheckedFace) -// continue; -// const SMESHDS_SubMesh * aSubMeshDSFace = meshDS->MeshElements(aShapeFace); -// if ( aSubMeshDSFace ) { -// SMDS_ElemIteratorPtr iteratorElem = aSubMeshDSFace->GetElements(); -// while ( iteratorElem->more() ) { // loop on elements on a face -// const SMDS_MeshElement* face = iteratorElem->next(); 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_MeshEditor::TNodeXYZ( face->GetNode(i) )); + aContour->Append( SMESH_TNodeXYZ( face->GetNode(i) )); if( HasIntersection(P, PC, Pres, aContour) ) { res = true; double tmp = PC.Distance(Pres); @@ -552,7 +591,6 @@ int StdMeshers_QuadToTriaAdaptor::Preparation(const SMDS_MeshElement* face { if( face->NbCornerNodes() != 4 ) { - myNbTriangles += int( face->NbCornerNodes() == 3 ); return NOT_QUAD; } @@ -560,12 +598,11 @@ int StdMeshers_QuadToTriaAdaptor::Preparation(const SMDS_MeshElement* face 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 ); + 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 ( aProxyMesh ) + aSubMeshDSFace = aProxyMesh->GetSubMesh( aShapeFace ); + else + aSubMeshDSFace = meshDS->MeshElements( aShapeFace ); + + vector trias, quads; + bool hasNewTrias = false; + if ( aSubMeshDSFace ) { - bool isRev = SMESH_Algo::IsReversedSubMesh( TopoDS::Face(aShapeFace), meshDS ); + bool isRev = false; + if ( helper.NbAncestors( aShapeFace, aMesh, aShape.ShapeType() ) > 1 ) + isRev = SMESH_Algo::IsReversedSubMesh( TopoDS::Face(aShapeFace), meshDS ); SMDS_ElemIteratorPtr iteratorElem = aSubMeshDSFace->GetElements(); 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; - - if(stat==2) + // preparation step to get face info + int stat = Preparation(face, PN, VN, FNodes, PC, VNorm); + switch ( stat ) { - // degenerate face - // add triangles to result map - SMDS_MeshFace* NewFace; - if(!isRev) - NewFace = meshDS->AddFace( FNodes[0], FNodes[1], FNodes[2] ); - else - NewFace = meshDS->AddFace( FNodes[0], FNodes[2], FNodes[1] ); - myTempTriangles[NewFace] =true; - TTriaList aList( 1, NewFace ); - myResMap.insert(make_pair(face,aList)); - continue; - } + case NOT_QUAD: - 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.); + trias.push_back( face ); + break; - // 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("<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; - check = CheckIntersection(PCbestTmp, PC, Pint, aMesh, aShape, face); - 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("<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); - // add triangles to result map - TTriaList& triaList = myResMap.insert( make_pair( face, TTriaList() ))->second; - for(i=0; i<4; i++) - { - SMDS_MeshFace* newFace =meshDS->AddFace( NewNode, FNodes[i], FNodes[i+1] ); - triaList.push_back( newFace ); - myTempTriangles[newFace] = true; - } + quads.push_back( face ); + hasNewTrias = true; + break; - // 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); + } // 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); } //================================================================================ @@ -792,8 +886,13 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh, const TopoDS_Shape bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh) { - myResMap.clear(); - myPyramids.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; + + vector myPyramids; SMESH_MesherHelper helper(aMesh); helper.IsQuadraticSubMesh(aMesh.GetShapeToMesh()); helper.SetElementsOnShape( true ); @@ -803,13 +902,13 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh) SMESH_ElementSearcher* searcher = const_cast(myElemSearcher); SMESHDS_Mesh * meshDS = aMesh.GetMeshDS(); + SMESH_ProxyMesh::SubMesh* prxSubMesh = getProxySubMesh(); SMDS_FaceIteratorPtr fIt = meshDS->facesIterator(/*idInceasingOrder=*/true); while( fIt->more()) { 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; @@ -852,7 +950,7 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh) 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; @@ -894,12 +992,13 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh) NewFace = meshDS->AddFace( FNodes[0], FNodes[1], FNodes[2] ); else NewFace = meshDS->AddFace( FNodes[0], FNodes[2], FNodes[1] ); - myTempTriangles[NewFace] = true; - aList.push_back(NewFace); - myResMap.insert(make_pair(face,aList)); + storeTmpElement( NewFace ); + prxSubMesh->AddElement( NewFace ); continue; } + // Case of non-degenerated quadrangle + // Find pyramid peak gp_XYZ PCbest(0., 0., 0.); // pyramid peak @@ -916,7 +1015,6 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh) PCbest = PC.XYZ() + VNorm.XYZ() * 0.001; height = PC.Distance(PCbest); } - //cout<<" PCbest("<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 ) { @@ -968,15 +1066,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_MeshFace* NewFace; if(isRev) NewFace = meshDS->AddFace( NewNode, FNodes[i], FNodes[i+1] ); else NewFace = meshDS->AddFace( NewNode, FNodes[i+1], FNodes[i] ); - myTempTriangles[NewFace] = true; - aList.push_back(NewFace); + storeTmpElement( NewFace ); + prxSubMesh->AddElement( NewFace ); } // create a pyramid SMDS_MeshVolume* aPyram; @@ -988,7 +1085,7 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh) } } // end loop on all faces - return Compute2ndPart(aMesh); + return Compute2ndPart(aMesh, myPyramids); } //================================================================================ @@ -997,7 +1094,8 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh) */ //================================================================================ -bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh& aMesh) +bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh& aMesh, + const vector& myPyramids) { if(myPyramids.empty()) return true; @@ -1005,8 +1103,8 @@ bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh& aMesh) SMESHDS_Mesh * meshDS = aMesh.GetMeshDS(); int i, j, k, myShapeID = myPyramids[0]->GetNode(4)->getshapeId(); - if ( !myElemSearcher ) - myElemSearcher = SMESH_MeshEditor(&aMesh).GetElementSearcher(); + if ( myElemSearcher ) delete myElemSearcher; + myElemSearcher = SMESH_MeshEditor(&aMesh).GetElementSearcher(); SMESH_ElementSearcher* searcher = const_cast(myElemSearcher); set nodesToMove; @@ -1016,7 +1114,7 @@ bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh& aMesh) for ( i = 0; i < myPyramids.size(); ++i ) { const SMDS_MeshElement* PrmI = myPyramids[i]; - MergeAdjacent( PrmI, aMesh, myTempTriangles, nodesToMove ); + MergeAdjacent( PrmI, nodesToMove ); } // iterate on all pyramids @@ -1032,18 +1130,20 @@ bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh& aMesh) for(k=0; k<5; k++) // loop on 4 base nodes of PrmI { const SMDS_MeshNode* n = PrmI->GetNode(k); - PsI[k] = SMESH_MeshEditor::TNodeXYZ( n ); + PsI[k] = SMESH_TNodeXYZ( n ); SMDS_ElemIteratorPtr vIt = n->GetInverseElementIterator( SMDSAbs_Volume ); while ( vIt->more() ) - checkedPyrams.insert( vIt->next() ); + { + 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 { gp_Vec Vtmp(PsI[k],PsI[4]); - gp_Pnt Pshift = PsI[k].XYZ() + Vtmp.XYZ() * 0.01; // base node moved a bit to apex - gp_Ax1 line( PsI[k], Vtmp ); vector< const SMDS_MeshElement* > suspectPyrams; searcher->GetElementsNearLine( line, SMDSAbs_Volume, suspectPyrams); @@ -1064,12 +1164,16 @@ bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh& aMesh) vector PsJ( xyzIt, TXyzIterator() ); gp_Pnt Pint; - bool hasInt = + 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; @@ -1093,7 +1197,7 @@ bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh& aMesh) if(nbc>0) { // Merge the two pyramids and others already merged with them - MergePyramids( PrmI, PrmJ, meshDS, myTempTriangles, nodesToMove ); + MergePiramids( PrmI, PrmJ, nodesToMove ); } else { // nbc==0 @@ -1126,8 +1230,8 @@ bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh& aMesh) nodesToMove.insert( aNode2 ); } // fix intersections that could appear after apex movement - MergeAdjacent( PrmI, aMesh, myTempTriangles, nodesToMove ); - MergeAdjacent( PrmJ, aMesh, myTempTriangles, nodesToMove ); + MergeAdjacent( PrmI, nodesToMove ); + MergeAdjacent( PrmJ, nodesToMove ); } // end if(hasInt) } // loop on suspectPyrams @@ -1142,42 +1246,32 @@ bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh& aMesh) 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 ) + for ( int i = 0; i <= meshDS->MaxShapeIndex(); ++i ) + if ( SMESH_ProxyMesh::SubMesh* sm = findProxySubMesh(i)) { - //cerr << __LINE__ << " splice" << endl; - q2tPrev->second.splice( q2tPrev->second.end(), q2t->second ); - } - } - // delete removed triangles and count resulting nb of triangles - for (q2t = myResMap.begin(); q2t != myResMap.end(); ++q2t) - { - TTriaList & trias = q2t->second; - vector faceToErase; - faceToErase.clear(); - //cerr << __LINE__ << " " << trias.size() << endl; - for (TTriaList::iterator tri = trias.begin(); tri != trias.end(); ++tri) + vector faces; + faces.reserve( sm->NbElements() ); + SMDS_ElemIteratorPtr fIt = sm->GetElements(); + while ( fIt->more() ) { - //cerr << " " << __LINE__ << endl; - const SMDS_MeshFace* face = *tri; - if (myTempTriangles.count(face) && (myTempTriangles[face] == false)) - faceToErase.push_back(face); + const SMDS_MeshElement* tria = fIt->next(); + set::iterator rmTria = myRemovedTrias.find( tria ); + if ( rmTria != myRemovedTrias.end() ) + myRemovedTrias.erase( rmTria ); else - myNbTriangles++; + faces.push_back( tria ); } - for (vector::iterator it = faceToErase.begin(); it != faceToErase.end(); ++it) - { - const SMDS_MeshFace *face = dynamic_cast(*it); - if (face) trias.remove(face); - meshDS->RemoveFreeElement(face, 0, false); - } - } + sm->ChangeElements( faces.begin(), faces.end() ); + } + } - myPyramids.clear(); // no more needed myDegNodes.clear(); delete myElemSearcher; @@ -1186,14 +1280,3 @@ bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh& aMesh) return true; } -//================================================================================ -/*! - * \brief Return list of created triangles for given face - */ -//================================================================================ - -const list* StdMeshers_QuadToTriaAdaptor::GetTriangles (const SMDS_MeshElement* aQuad) -{ - TQuad2Trias::iterator it = myResMap.find(aQuad); - return ( it != myResMap.end() ? & it->second : 0 ); -}