X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FSMESH%2FSMESH_Pattern.cxx;h=d6aee34df540c641a3dd3d648ad5dd98997bf57c;hp=86f1d52ab1688686683fce61988d215bd7a9e1f8;hb=aa67cc96d730566d202d1014d97b7f0b3a4d71f4;hpb=090aff07266d376ae028ae43434bdea7c0a0f9bb diff --git a/src/SMESH/SMESH_Pattern.cxx b/src/SMESH/SMESH_Pattern.cxx index 86f1d52ab..d6aee34df 100644 --- a/src/SMESH/SMESH_Pattern.cxx +++ b/src/SMESH/SMESH_Pattern.cxx @@ -1,44 +1,50 @@ -// 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.opencascade.org/SALOME/ or email : webmaster.salome@opencascade.org +// Copyright (C) 2007-2010 CEA/DEN, EDF R&D, OPEN CASCADE +// +// Copyright (C) 2003-2007 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 : SMESH_Pattern.hxx // Created : Mon Aug 2 10:30:00 2004 // Author : Edward AGAPOV (eap) - +// #include "SMESH_Pattern.hxx" +#include #include #include #include #include #include #include +#include #include #include #include #include #include #include -#include #include #include +#include #include +#include #include #include #include @@ -58,12 +64,13 @@ #include "SMDS_MeshElement.hxx" #include "SMDS_MeshFace.hxx" #include "SMDS_MeshNode.hxx" +#include "SMDS_VolumeTool.hxx" #include "SMESHDS_Group.hxx" #include "SMESHDS_Mesh.hxx" #include "SMESHDS_SubMesh.hxx" #include "SMESH_Block.hxx" #include "SMESH_Mesh.hxx" -#include "SMESH_MeshEditor.hxx" +#include "SMESH_MesherHelper.hxx" #include "SMESH_subMesh.hxx" #include "utilities.h" @@ -72,6 +79,8 @@ using namespace std; typedef map< const SMDS_MeshElement*, int > TNodePointIDMap; +#define smdsNode( elem ) static_cast( elem ) + //======================================================================= //function : SMESH_Pattern //purpose : @@ -243,7 +252,7 @@ bool SMESH_Pattern::Load (const char* theFileContents) MESSAGE(" Too few points "); return setErrorCode( ERR_READ_TOO_FEW_POINTS ); } - + // read the rest points int iPoint; for ( iPoint = 1; iPoint < nbPoints; iPoint++ ) @@ -297,8 +306,8 @@ bool SMESH_Pattern::Load (const char* theFileContents) while ( readLine( fields, lineBeg, clearFields )) { - myElemPointIDs.push_back( list< int >() ); - list< int >& elemPoints = myElemPointIDs.back(); + myElemPointIDs.push_back( TElemDef() ); + TElemDef& elemPoints = myElemPointIDs.back(); for ( fIt = fields.begin(); fIt != fields.end(); fIt++ ) { int pointIndex = getInt( *fIt ); @@ -379,18 +388,18 @@ bool SMESH_Pattern::Save (ostream& theFile) } // elements theFile << "!!! Indices of points of " << myElemPointIDs.size() << " elements:" << endl; - list >::const_iterator epIt = myElemPointIDs.begin(); + list::const_iterator epIt = myElemPointIDs.begin(); for ( ; epIt != myElemPointIDs.end(); epIt++ ) { - const list< int > & elemPoints = *epIt; - list< int >::const_iterator iIt = elemPoints.begin(); + const TElemDef & elemPoints = *epIt; + TElemDef::const_iterator iIt = elemPoints.begin(); for ( ; iIt != elemPoints.end(); iIt++ ) theFile << " " << *iIt; theFile << endl; } theFile << endl; - + return setErrorCode( ERR_OK ); } @@ -412,77 +421,6 @@ template void sortBySize( list< list < T > > & theListOfList ) } } -//======================================================================= -//function : getOrderedEdges -//purpose : return nb wires and a list of oredered edges -//======================================================================= - -static int getOrderedEdges (const TopoDS_Face& theFace, - const TopoDS_Vertex& theFirstVertex, - list< TopoDS_Edge >& theEdges, - list< int > & theNbVertexInWires) -{ - // put wires in a list, so that an outer wire comes first - list aWireList; - TopoDS_Wire anOuterWire = BRepTools::OuterWire( theFace ); - aWireList.push_back( anOuterWire ); - for ( TopoDS_Iterator wIt (theFace); wIt.More(); wIt.Next() ) - if ( !anOuterWire.IsSame( wIt.Value() )) - aWireList.push_back( TopoDS::Wire( wIt.Value() )); - - // loop on edges of wires - theNbVertexInWires.clear(); - list::iterator wlIt = aWireList.begin(); - for ( ; wlIt != aWireList.end(); wlIt++ ) - { - int iE; - BRepTools_WireExplorer wExp( *wlIt, theFace ); - for ( iE = 0; wExp.More(); wExp.Next(), iE++ ) - { - TopoDS_Edge edge = wExp.Current(); - edge = TopoDS::Edge( edge.Oriented( wExp.Orientation() )); - theEdges.push_back( edge ); - } - theNbVertexInWires.push_back( iE ); - iE = 0; - if ( wlIt == aWireList.begin() && theEdges.size() > 1 ) { // the outer wire - // orient closed edges - list< TopoDS_Edge >::iterator eIt, eIt2; - for ( eIt = theEdges.begin(); eIt != theEdges.end(); eIt++ ) - { - TopoDS_Edge& edge = *eIt; - if ( TopExp::FirstVertex( edge ).IsSame( TopExp::LastVertex( edge ) )) - { - eIt2 = eIt; - bool isNext = ( eIt2 == theEdges.begin() ); - TopoDS_Edge edge2 = isNext ? *(++eIt2) : *(--eIt2); - double f1,l1,f2,l2; - Handle(Geom2d_Curve) c1 = BRep_Tool::CurveOnSurface( edge, theFace, f1,l1 ); - Handle(Geom2d_Curve) c2 = BRep_Tool::CurveOnSurface( edge2, theFace, f2,l2 ); - gp_Pnt2d pf = c1->Value( edge.Orientation() == TopAbs_FORWARD ? f1 : l1 ); - gp_Pnt2d pl = c1->Value( edge.Orientation() == TopAbs_FORWARD ? l1 : f1 ); - bool isFirst = ( edge2.Orientation() == TopAbs_FORWARD ? isNext : !isNext ); - gp_Pnt2d p2 = c2->Value( isFirst ? f2 : l2 ); - isFirst = ( p2.SquareDistance( pf ) < p2.SquareDistance( pl )); - if ( isNext ? isFirst : !isFirst ) - edge.Reverse(); - } - } - // rotate theEdges until it begins from theFirstVertex - if ( ! theFirstVertex.IsNull() ) - while ( !theFirstVertex.IsSame( TopExp::FirstVertex( theEdges.front(), true ))) - { - theEdges.splice(theEdges.end(), theEdges, - theEdges.begin(), ++ theEdges.begin()); - if ( iE++ > theNbVertexInWires.back() ) - break; // break infinite loop - } - } - } - - return aWireList.size(); -} - //======================================================================= //function : project //purpose : @@ -507,34 +445,53 @@ static gp_XY project (const SMDS_MeshNode* theNode, } //======================================================================= -//function : isMeshBoundToShape -//purpose : return true if all 2d elements are bound to shape +//function : areNodesBound +//purpose : true if all nodes of faces are bound to shapes //======================================================================= -static bool isMeshBoundToShape(SMESH_Mesh* theMesh) +template bool areNodesBound( TFaceIterator & faceItr ) { - // check faces binding - SMESHDS_Mesh * aMeshDS = theMesh->GetMeshDS(); - SMESHDS_SubMesh * aMainSubMesh = aMeshDS->MeshElements( aMeshDS->ShapeToMesh() ); - if ( aMeshDS->NbFaces() != aMainSubMesh->NbElements() ) - return false; - - // check face nodes binding - SMDS_FaceIteratorPtr fIt = aMeshDS->facesIterator(); - while ( fIt->more() ) + while ( faceItr->more() ) { - SMDS_ElemIteratorPtr nIt = fIt->next()->nodesIterator(); + SMDS_ElemIteratorPtr nIt = faceItr->next()->nodesIterator(); while ( nIt->more() ) { - const SMDS_MeshNode* node = static_cast( nIt->next() ); + const SMDS_MeshNode* node = smdsNode( nIt->next() ); SMDS_PositionPtr pos = node->GetPosition(); - if ( !pos || !pos->GetShapeId() ) + if ( !pos || !pos->GetShapeId() ) { return false; + } } } return true; } +//======================================================================= +//function : isMeshBoundToShape +//purpose : return true if all 2d elements are bound to shape +// if aFaceSubmesh != NULL, then check faces bound to it +// else check all faces in aMeshDS +//======================================================================= + +static bool isMeshBoundToShape(SMESHDS_Mesh * aMeshDS, + SMESHDS_SubMesh * aFaceSubmesh, + const bool isMainShape) +{ + if ( isMainShape ) { + // check that all faces are bound to aFaceSubmesh + if ( aMeshDS->NbFaces() != aFaceSubmesh->NbElements() ) + return false; + } + + // check face nodes binding + if ( aFaceSubmesh ) { + SMDS_ElemIteratorPtr fIt = aFaceSubmesh->GetElements(); + return areNodesBound( fIt ); + } + SMDS_FaceIteratorPtr fIt = aMeshDS->facesIterator(); + return areNodesBound( fIt ); +} + //======================================================================= //function : Load //purpose : Create a pattern from the mesh built on . @@ -552,6 +509,8 @@ bool SMESH_Pattern::Load (SMESH_Mesh* theMesh, SMESHDS_Mesh * aMeshDS = theMesh->GetMeshDS(); SMESHDS_SubMesh * fSubMesh = aMeshDS->MeshElements( theFace ); + SMESH_MesherHelper helper( *theMesh ); + helper.SetSubShape( theFace ); int nbNodes = ( !fSubMesh ? 0 : fSubMesh->NbNodes() ); int nbElems = ( !fSubMesh ? 0 : fSubMesh->NbElements() ); @@ -561,40 +520,64 @@ bool SMESH_Pattern::Load (SMESH_Mesh* theMesh, return setErrorCode( ERR_LOAD_EMPTY_SUBMESH ); } - // check that face is not closed + TopoDS_Face face = TopoDS::Face( theFace.Oriented( TopAbs_FORWARD )); + + // check if face is closed + bool isClosed = helper.HasSeam(); TopoDS_Vertex bidon; list eList; - getOrderedEdges( theFace, bidon, eList, myNbKeyPntInBoundary ); - list::iterator elIt = eList.begin(); - for ( ; elIt != eList.end() ; elIt++ ) - if ( BRep_Tool::IsClosed( *elIt , theFace )) - return setErrorCode( ERR_LOADF_CLOSED_FACE ); - + list::iterator elIt; + SMESH_Block::GetOrderedEdges( face, bidon, eList, myNbKeyPntInBoundary ); + + // check that requested or needed projection is possible + bool isMainShape = theMesh->IsMainShape( face ); + bool needProject = !isMeshBoundToShape( aMeshDS, fSubMesh, isMainShape ); + bool canProject = ( nbElems ? true : isMainShape ); + if ( isClosed ) + canProject = false; // so far + + if ( ( theProject || needProject ) && !canProject ) + return setErrorCode( ERR_LOADF_CANT_PROJECT ); Extrema_GenExtPS projector; - GeomAdaptor_Surface aSurface( BRep_Tool::Surface( theFace )); - if ( theProject || nbElems == 0 ) + GeomAdaptor_Surface aSurface( BRep_Tool::Surface( face )); + if ( theProject || needProject ) projector.Initialize( aSurface, 20,20, 1e-5,1e-5 ); int iPoint = 0; TNodePointIDMap nodePointIDMap; + TNodePointIDMap closeNodePointIDMap; // for nodes on seam edges - if ( nbElems == 0 || (theProject && - theMesh->IsMainShape( theFace ) && - !isMeshBoundToShape( theMesh ))) + if ( needProject ) { - MESSAGE("Project the whole mesh"); + MESSAGE("Project the submesh"); // --------------------------------------------------------------- - // The case where the whole mesh is projected to theFace + // The case where the submesh is projected to theFace // --------------------------------------------------------------- - // put nodes of all faces in the nodePointIDMap and fill myElemPointIDs - SMDS_FaceIteratorPtr fIt = aMeshDS->facesIterator(); - while ( fIt->more() ) + // get all faces + list< const SMDS_MeshElement* > faces; + if ( nbElems > 0 ) { + SMDS_ElemIteratorPtr fIt = fSubMesh->GetElements(); + while ( fIt->more() ) { + const SMDS_MeshElement* f = fIt->next(); + if ( f && f->GetType() == SMDSAbs_Face ) + faces.push_back( f ); + } + } + else { + SMDS_FaceIteratorPtr fIt = aMeshDS->facesIterator(); + while ( fIt->more() ) + faces.push_back( fIt->next() ); + } + + // put nodes of all faces into the nodePointIDMap and fill myElemPointIDs + list< const SMDS_MeshElement* >::iterator fIt = faces.begin(); + for ( ; fIt != faces.end(); ++fIt ) { - myElemPointIDs.push_back( list< int >() ); - list< int >& elemPoints = myElemPointIDs.back(); - SMDS_ElemIteratorPtr nIt = fIt->next()->nodesIterator(); + myElemPointIDs.push_back( TElemDef() ); + TElemDef& elemPoints = myElemPointIDs.back(); + SMDS_ElemIteratorPtr nIt = (*fIt)->nodesIterator(); while ( nIt->more() ) { const SMDS_MeshElement* node = nIt->next(); @@ -614,18 +597,17 @@ bool SMESH_Pattern::Load (SMESH_Mesh* theMesh, TNodePointIDMap::iterator nIdIt = nodePointIDMap.begin(); for ( ; nIdIt != nodePointIDMap.end(); nIdIt++ ) { - const SMDS_MeshNode* node = - static_cast( (*nIdIt).first ); + const SMDS_MeshNode* node = smdsNode( (*nIdIt).first ); TPoint * p = & myPoints[ (*nIdIt).second ]; p->myInitUV = project( node, projector ); p->myInitXYZ.SetCoord( p->myInitUV.X(), p->myInitUV.Y(), 0 ); } // find key-points: the points most close to UV of vertices - TopExp_Explorer vExp( theFace, TopAbs_VERTEX ); + TopExp_Explorer vExp( face, TopAbs_VERTEX ); set foundIndices; for ( ; vExp.More(); vExp.Next() ) { const TopoDS_Vertex v = TopoDS::Vertex( vExp.Current() ); - gp_Pnt2d uv = BRep_Tool::Parameters( v, theFace ); + gp_Pnt2d uv = BRep_Tool::Parameters( v, face ); double minDist = DBL_MAX; int index; vector< TPoint >::const_iterator pVecIt = myPoints.begin(); @@ -652,16 +634,22 @@ bool SMESH_Pattern::Load (SMESH_Mesh* theMesh, // vertices for ( elIt = eList.begin(); elIt != eList.end(); elIt++ ) { + int nbV = myShapeIDMap.Extent(); myShapeIDMap.Add( TopExp::FirstVertex( *elIt, true )); - SMESHDS_SubMesh * eSubMesh = aMeshDS->MeshElements( *elIt ); - if ( eSubMesh ) + bool added = ( nbV < myShapeIDMap.Extent() ); + if ( !added ) { // vertex encountered twice + // a seam vertex have two corresponding key points + myShapeIDMap.Add( TopExp::FirstVertex( *elIt, true ).Reversed()); + ++nbNodes; + } + if ( SMESHDS_SubMesh * eSubMesh = aMeshDS->MeshElements( *elIt )) nbNodes += eSubMesh->NbNodes() + 1; } // edges for ( elIt = eList.begin(); elIt != eList.end(); elIt++ ) myShapeIDMap.Add( *elIt ); // the face - myShapeIDMap.Add( theFace ); + myShapeIDMap.Add( face ); myPoints.resize( nbNodes ); @@ -672,25 +660,58 @@ bool SMESH_Pattern::Load (SMESH_Mesh* theMesh, TopoDS_Edge & edge = *elIt; list< TPoint* > & ePoints = getShapePoints( edge ); double f, l; - Handle(Geom2d_Curve) C2d; - if ( !theProject ) - C2d = BRep_Tool::CurveOnSurface( edge, theFace, f, l ); + Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface( edge, face, f, l ); bool isForward = ( edge.Orientation() == TopAbs_FORWARD ); + TopoDS_Shape v1 = TopExp::FirstVertex( edge, true ); // always FORWARD + TopoDS_Shape v2 = TopExp::LastVertex( edge, true ); // always REVERSED + // to make adjacent edges share key-point, we make v2 FORWARD too + // (as we have different points for same shape with different orienation) + v2.Reverse(); + + // on closed face we must have REVERSED some of seam vertices + if ( isClosed ) { + if ( helper.IsSeamShape( edge ) ) { + if ( helper.IsRealSeam( edge ) && !isForward ) { + // reverse on reversed SEAM edge + v1.Reverse(); + v2.Reverse(); + } + } + else { // on CLOSED edge (i.e. having one vertex with different orienations) + for ( int is2 = 0; is2 < 2; ++is2 ) { + TopoDS_Shape & v = is2 ? v2 : v1; + if ( helper.IsRealSeam( v ) ) { + // reverse or not depending on orientation of adjacent seam + TopoDS_Edge seam; + list::iterator eIt2 = elIt; + if ( is2 ) + seam = ( ++eIt2 == eList.end() ? eList.front() : *eIt2 ); + else + seam = ( eIt2 == eList.begin() ? eList.back() : *(--eIt2) ); + if ( seam.Orientation() == TopAbs_REVERSED ) + v.Reverse(); + } + } + } + } + // the forward key-point - TopoDS_Shape v = TopExp::FirstVertex( edge, true ); - list< TPoint* > & vPoint = getShapePoints( v ); - if ( vPoint.empty() ) + list< TPoint* > * vPoint = & getShapePoints( v1 ); + if ( vPoint->empty() ) { - SMESHDS_SubMesh * vSubMesh = aMeshDS->MeshElements( v ); + SMESHDS_SubMesh * vSubMesh = aMeshDS->MeshElements( v1 ); if ( vSubMesh && vSubMesh->NbNodes() ) { myKeyPointIDs.push_back( iPoint ); SMDS_NodeIteratorPtr nIt = vSubMesh->GetNodes(); const SMDS_MeshNode* node = nIt->next(); - nodePointIDMap.insert( make_pair( node, iPoint )); + if ( v1.Orientation() == TopAbs_REVERSED ) + closeNodePointIDMap.insert( make_pair( node, iPoint )); + else + nodePointIDMap.insert( make_pair( node, iPoint )); TPoint* keyPoint = &myPoints[ iPoint++ ]; - vPoint.push_back( keyPoint ); + vPoint->push_back( keyPoint ); if ( theProject ) keyPoint->myInitUV = project( node, projector ); else @@ -698,8 +719,8 @@ bool SMESH_Pattern::Load (SMESH_Mesh* theMesh, keyPoint->myInitXYZ.SetCoord (keyPoint->myInitUV.X(), keyPoint->myInitUV.Y(), 0); } } - if ( !vPoint.empty() ) - ePoints.push_back( vPoint.front() ); + if ( !vPoint->empty() ) + ePoints.push_back( vPoint->front() ); // on-edge points SMESHDS_SubMesh * eSubMesh = aMeshDS->MeshElements( edge ); @@ -711,14 +732,37 @@ bool SMESH_Pattern::Load (SMESH_Mesh* theMesh, SMDS_NodeIteratorPtr nIt = eSubMesh->GetNodes(); while ( nIt->more() ) { - const SMDS_MeshNode* node = - static_cast( nIt->next() ); + const SMDS_MeshNode* node = smdsNode( nIt->next() ); const SMDS_EdgePosition* epos = - static_cast(node->GetPosition().get()); + static_cast(node->GetPosition()); double u = epos->GetUParameter(); - paramNodeMap.insert( TParamNodeMap::value_type( u, node )); + paramNodeMap.insert( make_pair( u, node )); + } + if ( paramNodeMap.size() != eSubMesh->NbNodes() ) { + // wrong U on edge, project + Extrema_ExtPC proj; + BRepAdaptor_Curve aCurve( edge ); + proj.Initialize( aCurve, f, l ); + paramNodeMap.clear(); + nIt = eSubMesh->GetNodes(); + for ( int iNode = 0; nIt->more(); ++iNode ) { + const SMDS_MeshNode* node = smdsNode( nIt->next() ); + proj.Perform( gp_Pnt( node->X(), node->Y(), node->Z())); + double u = 0; + if ( proj.IsDone() ) { + for ( int i = 1, nb = proj.NbExt(); i <= nb; ++i ) + if ( proj.IsMin( i )) { + u = proj.Point( i ).Parameter(); + break; + } + } else { + u = isForward ? iNode : eSubMesh->NbNodes() - iNode; + } + paramNodeMap.insert( make_pair( u, node )); + } } // put U in [0,1] so that the first key-point has U==0 + bool isSeam = helper.IsRealSeam( edge ); double du = l - f; TParamNodeMap::iterator unIt = paramNodeMap.begin(); TParamNodeMap::reverse_iterator unRIt = paramNodeMap.rbegin(); @@ -727,7 +771,10 @@ bool SMESH_Pattern::Load (SMESH_Mesh* theMesh, TPoint* p = & myPoints[ iPoint ]; ePoints.push_back( p ); const SMDS_MeshNode* node = isForward ? (*unIt).second : (*unRIt).second; - nodePointIDMap.insert ( make_pair( node, iPoint )); + if ( isSeam && !isForward ) + closeNodePointIDMap.insert( make_pair( node, iPoint )); + else + nodePointIDMap.insert ( make_pair( node, iPoint )); if ( theProject ) p->myInitUV = project( node, projector ); @@ -742,19 +789,21 @@ bool SMESH_Pattern::Load (SMESH_Mesh* theMesh, } } // the reverse key-point - v = TopExp::LastVertex( edge, true ).Reversed(); - list< TPoint* > & vPoint2 = getShapePoints( v ); - if ( vPoint2.empty() ) + vPoint = & getShapePoints( v2 ); + if ( vPoint->empty() ) { - SMESHDS_SubMesh * vSubMesh = aMeshDS->MeshElements( v ); + SMESHDS_SubMesh * vSubMesh = aMeshDS->MeshElements( v2 ); if ( vSubMesh && vSubMesh->NbNodes() ) { myKeyPointIDs.push_back( iPoint ); SMDS_NodeIteratorPtr nIt = vSubMesh->GetNodes(); const SMDS_MeshNode* node = nIt->next(); - nodePointIDMap.insert( make_pair( node, iPoint )); + if ( v2.Orientation() == TopAbs_REVERSED ) + closeNodePointIDMap.insert( make_pair( node, iPoint )); + else + nodePointIDMap.insert( make_pair( node, iPoint )); TPoint* keyPoint = &myPoints[ iPoint++ ]; - vPoint2.push_back( keyPoint ); + vPoint->push_back( keyPoint ); if ( theProject ) keyPoint->myInitUV = project( node, projector ); else @@ -762,8 +811,8 @@ bool SMESH_Pattern::Load (SMESH_Mesh* theMesh, keyPoint->myInitXYZ.SetCoord( keyPoint->myInitUV.X(), keyPoint->myInitUV.Y(), 0 ); } } - if ( !vPoint2.empty() ) - ePoints.push_back( vPoint2.front() ); + if ( !vPoint->empty() ) + ePoints.push_back( vPoint->front() ); // compute U of edge-points if ( theProject ) @@ -790,12 +839,11 @@ bool SMESH_Pattern::Load (SMESH_Mesh* theMesh, if ( fSubMesh && fSubMesh->NbElements() ) { - list< TPoint* > & fPoints = getShapePoints( theFace ); + list< TPoint* > & fPoints = getShapePoints( face ); SMDS_NodeIteratorPtr nIt = fSubMesh->GetNodes(); while ( nIt->more() ) { - const SMDS_MeshNode* node = - static_cast( nIt->next() ); + const SMDS_MeshNode* node = smdsNode( nIt->next() ); nodePointIDMap.insert( make_pair( node, iPoint )); TPoint* p = &myPoints[ iPoint++ ]; fPoints.push_back( p ); @@ -803,19 +851,48 @@ bool SMESH_Pattern::Load (SMESH_Mesh* theMesh, p->myInitUV = project( node, projector ); else { const SMDS_FacePosition* pos = - static_cast(node->GetPosition().get()); + static_cast(node->GetPosition()); p->myInitUV.SetCoord( pos->GetUParameter(), pos->GetVParameter() ); } p->myInitXYZ.SetCoord( p->myInitUV.X(), p->myInitUV.Y(), 0 ); } // load elements + TNodePointIDMap::iterator n_id, not_found = closeNodePointIDMap.end(); SMDS_ElemIteratorPtr elemIt = fSubMesh->GetElements(); - while ( elemIt->more() ) { - SMDS_ElemIteratorPtr nIt = elemIt->next()->nodesIterator(); - myElemPointIDs.push_back( list< int >() ); - list< int >& elemPoints = myElemPointIDs.back(); + while ( elemIt->more() ) + { + const SMDS_MeshElement* elem = elemIt->next(); + SMDS_ElemIteratorPtr nIt = elem->nodesIterator(); + myElemPointIDs.push_back( TElemDef() ); + TElemDef& elemPoints = myElemPointIDs.back(); + // find point indices corresponding to element nodes while ( nIt->more() ) - elemPoints.push_back( nodePointIDMap[ nIt->next() ]); + { + const SMDS_MeshNode* node = smdsNode( nIt->next() ); + iPoint = nodePointIDMap[ node ]; // point index of interest + // for a node on a seam edge there are two points + if ( helper.IsRealSeam( node->GetPosition()->GetShapeId() ) && + ( n_id = closeNodePointIDMap.find( node )) != not_found ) + { + TPoint & p1 = myPoints[ iPoint ]; + TPoint & p2 = myPoints[ n_id->second ]; + // Select point closest to the rest nodes of element in UV space + SMDS_ElemIteratorPtr nIt2 = elem->nodesIterator(); + const SMDS_MeshNode* notSeamNode = 0; + // find node not on a seam edge + while ( nIt2->more() && !notSeamNode ) { + const SMDS_MeshNode* n = smdsNode( nIt2->next() ); + if ( !helper.IsSeamShape( n->GetPosition()->GetShapeId() )) + notSeamNode = n; + } + gp_Pnt2d uv = helper.GetNodeUV( theFace, node, notSeamNode ); + double dist1 = uv.SquareDistance( p1.myInitUV ); + double dist2 = uv.SquareDistance( p2.myInitUV ); + if ( dist2 < dist1 ) + iPoint = n_id->second; + } + elemPoints.push_back( iPoint ); + } } } @@ -833,7 +910,21 @@ bool SMESH_Pattern::Load (SMESH_Mesh* theMesh, double dU = maxU - minU, dV = maxV - minV; if ( dU <= DBL_MIN || dV <= DBL_MIN ) { Clear(); - return setErrorCode( ERR_LOADF_NARROW_FACE ); + bndBox.SetVoid(); + // define where is the problem, in the face or in the mesh + TopExp_Explorer vExp( face, TopAbs_VERTEX ); + for ( ; vExp.More(); vExp.Next() ) { + gp_Pnt2d uv = BRep_Tool::Parameters( TopoDS::Vertex( vExp.Current() ), face ); + bndBox.Add( uv ); + } + bndBox.Get( minU, minV, maxU, maxV ); + dU = maxU - minU, dV = maxV - minV; + if ( dU <= DBL_MIN || dV <= DBL_MIN ) + // face problem + return setErrorCode( ERR_LOADF_NARROW_FACE ); + else + // mesh is projected onto a line, e.g. + return setErrorCode( ERR_LOADF_CANT_PROJECT ); } double ratio = dU / dV, maxratio = 3, scale; int iCoord = 0; @@ -900,13 +991,20 @@ static bool intersectIsolines(const gp_XY& uv11, const gp_XY& uv12, const double gp_XY loc1 = uv11 * ( 1 - r1 ) + uv12 * r1; gp_XY loc2 = uv21 * ( 1 - r2 ) + uv22 * r2; resUV = 0.5 * ( loc1 + loc2 ); - isDeformed = ( loc1 - loc2 ).SquareModulus() > 1e-8; + //isDeformed = ( loc1 - loc2 ).SquareModulus() > 1e-8; + // SKL 26.07.2007 for NPAL16567 + double d1 = (uv11-uv12).Modulus(); + double d2 = (uv21-uv22).Modulus(); + // double delta = d1*d2*1e-6; PAL17233 + double delta = min( d1, d2 ) / 10.; + isDeformed = ( loc1 - loc2 ).SquareModulus() > delta * delta; + // double len1 = ( uv11 - uv12 ).Modulus(); // double len2 = ( uv21 - uv22 ).Modulus(); // resUV = loc1 * len2 / ( len1 + len2 ) + loc2 * len1 / ( len1 + len2 ); // return true; - + // gp_Lin2d line1( uv11, uv12 - uv11 ); // gp_Lin2d line2( uv21, uv22 - uv21 ); // double angle = Abs( line1.Angle( line2 ) ); @@ -920,9 +1018,13 @@ static bool intersectIsolines(const gp_XY& uv11, const gp_XY& uv12, const double // inter.Perform( line1, line2 ); // interUV = inter.Point(1).Value(); // resUV += interUV.XY(); - + // resUV /= 2.; // } + if ( isDeformed ) { + MESSAGE("intersectIsolines(), d1 = " << d1 << ", d2 = " << d2 << ", delta = " << delta << + ", " << (loc1 - loc2).SquareModulus() << " > " << delta * delta); + } return true; } @@ -956,7 +1058,7 @@ bool SMESH_Pattern::compUVByIsoIntersection (const list< list< TPoint* > >& theB const list< TPoint* > & bndPoints = * bndIt; TPoint* prevP = bndPoints.back(); // this is the first point list< TPoint* >::const_iterator pIt = bndPoints.begin(); - bool coincPrev = false; + bool coincPrev = false; // loop on the edge-points for ( ; pIt != bndPoints.end(); pIt++ ) { @@ -1221,7 +1323,7 @@ static bool checkQuads (const TIsoNode* node, gp_XY uv1, uv2 = node->myUV; for ( i = isTriangle ? 2 : 0; i < 3; i++ ) // mark not computed vectors if ( wasOk[i] ) - moveVec[ i ].SetCoord( 1, 2e100); // not use this vector + moveVec[ i ].SetCoord( 1, 2e100); // not use this vector while ( !isOldOk ) { // find the least moveVec int i, iMin = 4; @@ -1322,6 +1424,7 @@ bool SMESH_Pattern:: compUVByElasticIsolines(const list< list< TPoint* > >& theBndPoints, const list< TPoint* >& thePntToCompute) { + return false; // PAL17233 //cout << "============================== KEY POINTS =============================="<::iterator kpIt = myKeyPointIDs.begin(); // for ( ; kpIt != myKeyPointIDs.end(); kpIt++ ) { @@ -1549,6 +1652,11 @@ bool SMESH_Pattern:: // " dir0: "<myDir[0].X()<<" "<myDir[0].Y() << // " dir1: "<myDir[1].X()<<" "<myDir[1].Y() << endl; } + else { + /// WHAT IN THIS CASE ????????????? MAY BE THIS, I AM NOT SURE :( + node->SetBoundaryNode( 0, iDir, 0 ); + node->SetBoundaryNode( 0, iDir, 1 ); + } } nIt++; nPrevIt++; if ( nNextIt != isoLine.end() ) nNextIt++; @@ -1646,7 +1754,7 @@ bool SMESH_Pattern:: aNorm[1-iDir].Normalize(); double r = Abs ( ratio[iDir] - 0.5 ) * 2.0; // [0,1] - distance from the middle r *= r; - + node->myDir[iDir] = //aTgt[iDir]; aNorm[1-iDir] * r + aTgt[iDir] * ( 1. - r ); } @@ -1798,6 +1906,10 @@ bool SMESH_Pattern:: // dir = node->myDir[ 1 - iDir ].XY() * ( isEnd ? -1. : 1. ); //cout << "__________"<myInitUV.X()<<" "<myInitUV.Y()<GetBoundaryNode( iDir, isEnd ); + if ( !bndNode ) { + MESSAGE("Why we are here?"); + continue; + } gp_XY tgt( bndNode->myDir[0].XY() + bndNode->myDir[1].XY() ); dir.SetCoord( 1, tgt.Y() * ( reversed ? 1 : -1 )); dir.SetCoord( 2, tgt.X() * ( reversed ? -1 : 1 )); @@ -1827,6 +1939,7 @@ bool SMESH_Pattern:: nbComp++; } } + if ( !nbComp ) continue; newUV /= nbComp; node->myUV = newUV; //cout << "NODE: "<myInitUV.X()<<" "<myInitUV.Y()<myDir[ iDir ] ); } // define ratio + bool ok = true; // <- stupid fix TO AVOID PB OF NODES WITH NULL BND NODES double locR[2] = { 0, 0 }; for ( iDir = 0; iDir < 2; iDir++ ) { const int iCoord = 2 - iDir; // coord changing along an isoline TIsoNode* bndNode1 = node->GetBoundaryNode( iDir, 0 ); TIsoNode* bndNode2 = node->GetBoundaryNode( iDir, 1 ); + if ( !bndNode1 || !bndNode2 ) { + ok = false; break; + } double par1 = bndNode1->myInitUV.Coord( iCoord ); double par2 = node->myInitUV.Coord( iCoord ); double par3 = bndNode2->myInitUV.Coord( iCoord ); @@ -1899,7 +2016,7 @@ bool SMESH_Pattern:: //locR[0] = locR[1] = 0.25; // intersect the 2 lines and move a node //IntAna2d_AnaIntersection inter( line[0], line[1] ); - if ( /*inter.IsDone() && inter.NbPoints() ==*/ 1 ) + if ( ok /*inter.IsDone() && inter.NbPoints() ==*/ ) { // double intR = 1 - locR[0] - locR[1]; // gp_XY newUV = inter.Point(1).Value().XY(); @@ -1947,8 +2064,7 @@ bool SMESH_Pattern:: } } } - - + return true; } @@ -1974,7 +2090,7 @@ double SMESH_Pattern::setFirstEdge (list< TopoDS_Edge > & theWire, int theFirstE int eID = theFirstEdgeID; for ( iE = 0; iE < nbEdges; iE++ ) maxNbPnt = Max ( maxNbPnt, getShapePoints( eID++ ).size() ); - + // compute bnd boxes TopoDS_Face face = TopoDS::Face( myShape ); Bnd_Box2d bndBox, eBndBox; @@ -2002,8 +2118,8 @@ double SMESH_Pattern::setFirstEdge (list< TopoDS_Edge > & theWire, int theFirstE bndBox.Get( minPar[0], minPar[1], maxPar[0], maxPar[1] ); eBndBox.Get( eMinPar[0], eMinPar[1], eMaxPar[0], eMaxPar[1] ); #ifdef DBG_SETFIRSTEDGE - cout << "EDGES: X: " << eMinPar[0] << " - " << eMaxPar[0] << " Y: " - << eMinPar[1] << " - " << eMaxPar[1] << endl; + MESSAGE ( "EDGES: X: " << eMinPar[0] << " - " << eMaxPar[0] << " Y: " + << eMinPar[1] << " - " << eMaxPar[1] ); #endif for ( int iC = 1, i = 0; i < 2; iC++, i++ ) // loop on 2 coordinates { @@ -2029,7 +2145,7 @@ double SMESH_Pattern::setFirstEdge (list< TopoDS_Edge > & theWire, int theFirstE for ( iE = 0 ; iE < nbEdges; iE++ ) { #ifdef DBG_SETFIRSTEDGE - cout << " VARIANT " << iE << endl; + MESSAGE ( " VARIANT " << iE ); #endif // evaluate the distance between UV computed by the 2 methods: // by isos intersection ( myXYZ ) and by edge p-curves ( myUV ) @@ -2043,13 +2159,13 @@ double SMESH_Pattern::setFirstEdge (list< TopoDS_Edge > & theWire, int theFirstE TPoint* p = (*pIt); dist += ( p->myUV - gp_XY( p->myXYZ.X(), p->myXYZ.Y() )).SquareModulus(); #ifdef DBG_SETFIRSTEDGE - cout << " ISO : ( " << p->myXYZ.X() << ", "<< p->myXYZ.Y() << " ) PCURVE : ( " << - p->myUV.X() << ", " << p->myUV.Y() << ") " << endl; + MESSAGE ( " ISO : ( " << p->myXYZ.X() << ", "<< p->myXYZ.Y() << " ) PCURVE : ( " << + p->myUV.X() << ", " << p->myUV.Y() << ") " ); #endif } } #ifdef DBG_SETFIRSTEDGE - cout << "dist -- " << dist << endl; + MESSAGE ( "dist -- " << dist ); #endif if ( dist < minDist ) { minDist = dist; @@ -2169,7 +2285,7 @@ bool SMESH_Pattern::sortSameSizeWires (TListOfEdgesList & theWire bndIndWirePosMap.insert( TIntWirePosMap::value_type( bIndex, wlIt )); } - // Treat each wire + // Treat each wire TIntWirePosMap::iterator bIndWPosIt = bndIndWirePosMap.begin(); eID = theFirstEdgeID; @@ -2180,7 +2296,7 @@ bool SMESH_Pattern::sortSameSizeWires (TListOfEdgesList & theWire // choose the best first edge of a wire setFirstEdge( wire, eID ); - + // compute eventual UV and fill theEdgesPointsList theEdgesPointsList.push_back( list< TPoint* >() ); list< TPoint* > & edgesPoints = theEdgesPointsList.back(); @@ -2223,7 +2339,7 @@ bool SMESH_Pattern::Apply (const TopoDS_Face& theFace, list< TopoDS_Edge > eList; list< int > nbVertexInWires; - int nbWires = getOrderedEdges( face, theVertexOnKeyPoint1, eList, nbVertexInWires); + int nbWires = SMESH_Block::GetOrderedEdges( face, theVertexOnKeyPoint1, eList, nbVertexInWires); if ( !theVertexOnKeyPoint1.IsSame( TopExp::FirstVertex( eList.front(), true ))) { MESSAGE( " theVertexOnKeyPoint1 not found in the outer wire "); @@ -2242,8 +2358,24 @@ bool SMESH_Pattern::Apply (const TopoDS_Face& theFace, list::iterator elIt = eList.begin(); for ( ; elIt != eList.end(); elIt++ ) { myShapeIDMap.Add( TopExp::FirstVertex( *elIt, true )); - if ( BRep_Tool::IsClosed( *elIt, theFace ) ) - myShapeIDMap.Add( TopExp::LastVertex( *elIt, true )); + bool isClosed1 = BRep_Tool::IsClosed( *elIt, theFace ); + // BEGIN: jfa for bug 0019943 + if (isClosed1) { + isClosed1 = false; + for (TopExp_Explorer expw (theFace, TopAbs_WIRE); expw.More() && !isClosed1; expw.Next()) { + const TopoDS_Wire& wire = TopoDS::Wire(expw.Current()); + int nbe = 0; + for (BRepTools_WireExplorer we (wire, theFace); we.More() && !isClosed1; we.Next()) { + if (we.Current().IsSame(*elIt)) { + nbe++; + if (nbe == 2) isClosed1 = true; + } + } + } + } + // END: jfa for bug 0019943 + if (isClosed1) + myShapeIDMap.Add( TopExp::LastVertex( *elIt, true ));// vertex orienation is REVERSED } int nbVertices = myShapeIDMap.Extent(); @@ -2252,7 +2384,7 @@ bool SMESH_Pattern::Apply (const TopoDS_Face& theFace, myShapeIDMap.Add( face ); - if ( myShapeIDToPointsMap.size() != myShapeIDMap.Extent()/* + nbSeamShapes*/ ) { + if ( myShapeIDToPointsMap.size() != myShapeIDMap.Extent() ) { MESSAGE( myShapeIDToPointsMap.size() <<" != " << myShapeIDMap.Extent()); return setErrorCode( ERR_APPLF_INTERNAL_EEROR ); } @@ -2349,7 +2481,7 @@ bool SMESH_Pattern::Apply (const TopoDS_Face& theFace, } // find boundary - wire correspondence for several wires of same size - + id1 = nbVertices + nbEdgesInOuterWire + 1; wlIt = wireList.begin(); while ( wlIt != wireList.end() ) @@ -2369,7 +2501,7 @@ bool SMESH_Pattern::Apply (const TopoDS_Face& theFace, } // add well-ordered edges to eList - + for ( wlIt = wireList.begin(); wlIt != wireList.end(); wlIt++ ) { list< TopoDS_Edge >& wire = (*wlIt); @@ -2384,7 +2516,7 @@ bool SMESH_Pattern::Apply (const TopoDS_Face& theFace, for ( elIt = eList.begin(); elIt != eList.end(); elIt++ ) myShapeIDMap.Add( *elIt ); myShapeIDMap.Add( face ); - + } // there are inner wires // Compute XYZ of on-edge points @@ -2392,17 +2524,13 @@ bool SMESH_Pattern::Apply (const TopoDS_Face& theFace, TopLoc_Location loc; for ( iE = nbVertices + 1, elIt = eList.begin(); elIt != eList.end(); elIt++ ) { - double f,l; - Handle(Geom_Curve) C3d = BRep_Tool::Curve( *elIt, loc, f, l ); - const gp_Trsf & aTrsf = loc.Transformation(); + BRepAdaptor_Curve C3d( *elIt ); list< TPoint* > & ePoints = getShapePoints( iE++ ); pIt = ePoints.begin(); for ( pIt++; pIt != ePoints.end(); pIt++ ) { TPoint* point = *pIt; - point->myXYZ = C3d->Value( point->myU ); - if ( !loc.IsIdentity() ) - aTrsf.Transforms( point->myXYZ.ChangeCoord() ); + point->myXYZ = C3d.Value( point->myU ); } } @@ -2454,7 +2582,7 @@ bool SMESH_Pattern::Apply (const SMDS_MeshFace* theFace, const int theNodeIndexOnKeyPoint1, const bool theReverse) { - MESSAGE(" ::Apply(MeshFace) " ); +// MESSAGE(" ::Apply(MeshFace) " ); if ( !IsLoaded() ) { MESSAGE( "Pattern not loaded" ); @@ -2483,7 +2611,7 @@ bool SMESH_Pattern::Apply (const SMDS_MeshFace* theFace, SMDS_ElemIteratorPtr noIt = theFace->nodesIterator(); int iSub = 0; while ( noIt->more() ) { - const SMDS_MeshNode* node = static_cast( noIt->next() ); + const SMDS_MeshNode* node = smdsNode( noIt->next() ); nodes.push_back( node ); if ( iSub++ == theNodeIndexOnKeyPoint1 ) n = --nodes.end(); @@ -2538,7 +2666,7 @@ bool SMESH_Pattern::Apply (const SMDS_MeshFace* theFace, { gp_XYZ& xyz1 = *xyzIt++; gp_XYZ& xyz2 = ( xyzIt != xyzList.end() ) ? *xyzIt : xyzList.front(); - + list< TPoint* > & ePoints = getShapePoints( iSub ); ePoints.back()->myInitU = 1.0; list< TPoint* >::const_iterator pIt = ++ePoints.begin(); @@ -2586,119 +2714,39 @@ bool SMESH_Pattern::Apply (const SMDS_MeshFace* theFace, return setErrorCode( ERR_OK ); } -//======================================================================= -//function : undefinedXYZ -//purpose : -//======================================================================= - -static const gp_XYZ& undefinedXYZ() -{ - static gp_XYZ xyz( 1.e100, 0., 0. ); - return xyz; -} - -//======================================================================= -//function : isDefined -//purpose : -//======================================================================= - -inline static bool isDefined(const gp_XYZ& theXYZ) -{ - return theXYZ.X() < 1.e100; -} - -//======================================================================= -//function : mergePoints -//purpose : Look for coincident points between myXYZs indexed with -// list of each element of xyzIndGroups. Coincident indices -// are merged in myElemXYZIDs. -//======================================================================= - -void SMESH_Pattern::mergePoints (map > >& indGroups, - map< int, list< list< int >* > > & reverseConnectivity) -{ - map< TNodeSet, list< list< int > > >::iterator indListIt; - for ( indListIt = indGroups.begin(); indListIt != indGroups.end(); indListIt++ ) - { - list > groups = indListIt->second; - if ( groups.size() < 2 ) - continue; - -// const TNodeSet & nodes = indListIt->first; -// TNodeSet::const_iterator n = nodes.begin(); -// for ( ; n != nodes.end(); n++ ) -// cout << *n ; - - // find tolerance - Bnd_Box box; - list< int >& indices = groups.front(); - list< int >::iterator ind, ind1, ind2; - for ( ind = indices.begin(); ind != indices.end(); ind++ ) - box.Add( gp_Pnt( myXYZ[ *ind ])); - double x, y, z, X, Y, Z; - box.Get( x, y, z, X, Y, Z ); - gp_Pnt p( x, y, z ), P( X, Y, Z ); - double tol2 = 1.e-4 * p.SquareDistance( P ); - - // compare points, replace indices - - list< list< int > >::iterator grpIt1, grpIt2; - for ( grpIt1 = groups.begin(); grpIt1 != groups.end(); grpIt1++ ) - { - list< int >& indices1 = *grpIt1; - grpIt2 = grpIt1; - for ( grpIt2++; grpIt2 != groups.end(); grpIt2++ ) - { - list< int >& indices2 = *grpIt2; - for ( ind1 = indices1.begin(); ind1 != indices1.end(); ind1++ ) - { - gp_XYZ& p1 = myXYZ[ *ind1 ]; - ind2 = indices2.begin(); - while ( ind2 != indices2.end() ) - { - gp_XYZ& p2 = myXYZ[ *ind2 ]; - //MESSAGE("COMP: " << *ind1 << " " << *ind2 << " X: " << p2.X() << " tol2: " << tol2); - if ( ( p1 - p2 ).SquareModulus() <= tol2 ) - { - ASSERT( reverseConnectivity.find( *ind2 ) != reverseConnectivity.end() ); - list< list< int >* > & elemXYZIDsList = reverseConnectivity[ *ind2 ]; - list< list< int >* >::iterator elemXYZIDs = elemXYZIDsList.begin(); - for ( ; elemXYZIDs != elemXYZIDsList.end(); elemXYZIDs++ ) - { - ind = find( (*elemXYZIDs)->begin(), (*elemXYZIDs)->end(), *ind2 ); - //MESSAGE( " Replace " << *ind << " with " << *ind1 ); - myXYZ[ *ind ] = undefinedXYZ(); - *ind = *ind1; - } - ind2 = indices2.erase( ind2 ); - } - else - ind2++; - } - } - } - } - } -} - //======================================================================= //function : Apply //purpose : Compute nodes coordinates applying -// the loaded pattern to . The first key-point +// the loaded pattern to . The first key-point // will be mapped into -th node //======================================================================= -bool SMESH_Pattern::Apply (std::set theFaces, - const int theNodeIndexOnKeyPoint1, - const bool theReverse) +bool SMESH_Pattern::Apply (SMESH_Mesh* theMesh, + const SMDS_MeshFace* theFace, + const TopoDS_Shape& theSurface, + const int theNodeIndexOnKeyPoint1, + const bool theReverse) { - MESSAGE(" ::Apply(set) " ); +// MESSAGE(" ::Apply(MeshFace) " ); + if ( theSurface.IsNull() || theSurface.ShapeType() != TopAbs_FACE ) { + return Apply( theFace, theNodeIndexOnKeyPoint1, theReverse); + } + const TopoDS_Face& face = TopoDS::Face( theSurface ); + TopLoc_Location loc; + Handle(Geom_Surface) surface = BRep_Tool::Surface( face, loc ); + const gp_Trsf & aTrsf = loc.Transformation(); if ( !IsLoaded() ) { MESSAGE( "Pattern not loaded" ); return setErrorCode( ERR_APPL_NOT_LOADED ); } + // check nb of nodes + if (theFace->NbNodes() != myNbKeyPntInBoundary.front() ) { + MESSAGE( myKeyPointIDs.size() << " != " << theFace->NbNodes() ); + return setErrorCode( ERR_APPL_BAD_NB_VERTICES ); + } + // find points on edges, it fills myNbKeyPntInBoundary if ( !findBoundaryPoints() ) return false; @@ -2708,68 +2756,255 @@ bool SMESH_Pattern::Apply (std::set theFaces, return setErrorCode( ERR_APPL_BAD_NB_VERTICES ); } - myXYZ.clear(); - myElemXYZIDs.clear(); - myXYZIdToNodeMap.clear(); - myElements.clear(); + // Define the nodes order - myXYZ.resize( myPoints.size() * theFaces.size(), undefinedXYZ() ); - myElements.reserve( theFaces.size() ); + list< const SMDS_MeshNode* > nodes; + list< const SMDS_MeshNode* >::iterator n = nodes.end(); + SMDS_ElemIteratorPtr noIt = theFace->nodesIterator(); + int iSub = 0; + while ( noIt->more() ) { + const SMDS_MeshNode* node = smdsNode( noIt->next() ); + nodes.push_back( node ); + if ( iSub++ == theNodeIndexOnKeyPoint1 ) + n = --nodes.end(); + } + if ( n != nodes.end() ) { + if ( theReverse ) { + if ( n != --nodes.end() ) + nodes.splice( nodes.begin(), nodes, ++n, nodes.end() ); + nodes.reverse(); + } + else if ( n != nodes.begin() ) + nodes.splice( nodes.end(), nodes, nodes.begin(), n ); + } - // to find point index - map< TPoint*, int > pointIndex; - for ( int i = 0; i < myPoints.size(); i++ ) - pointIndex.insert( make_pair( & myPoints[ i ], i )); + // find a node not on a seam edge, if necessary + SMESH_MesherHelper helper( *theMesh ); + helper.SetSubShape( theSurface ); + const SMDS_MeshNode* inFaceNode = 0; + if ( helper.GetNodeUVneedInFaceNode() ) + { + SMESH_MeshEditor editor( theMesh ); + for ( n = nodes.begin(); ( !inFaceNode && n != nodes.end()); ++n ) { + int shapeID = editor.FindShape( *n ); + if ( !shapeID ) + return Apply( theFace, theNodeIndexOnKeyPoint1, theReverse); + if ( !helper.IsSeamShape( shapeID )) + inFaceNode = *n; + } + } + + // Set UV of key-points (i.e. of nodes of theFace ) + vector< gp_XY > keyUV( theFace->NbNodes() ); + myOrderedNodes.resize( theFace->NbNodes() ); + for ( iSub = 1, n = nodes.begin(); n != nodes.end(); ++n, ++iSub ) + { + TPoint* p = getShapePoints( iSub ).front(); + p->myUV = helper.GetNodeUV( face, *n, inFaceNode ); + p->myXYZ = gp_XYZ( (*n)->X(), (*n)->Y(), (*n)->Z() ); - // to merge nodes on edges of the elements being refined - typedef set TLink; - map< TLink, list< list< int > > > linkPointIndListMap; - map< int, list< list< int >* > > reverseConnectivity; + keyUV[ iSub-1 ] = p->myUV; + myOrderedNodes[ iSub-1 ] = *n; + } - int ind1 = 0; // lowest point index for a face + // points on edges to be used for UV computation of in-face points + list< list< TPoint* > > edgesPointsList; + edgesPointsList.push_back( list< TPoint* >() ); + list< TPoint* > * edgesPoints = & edgesPointsList.back(); + list< TPoint* >::iterator pIt; - // apply to each face in theFaces set - set::iterator face = theFaces.begin(); - for ( ; face != theFaces.end(); ++face ) + // compute UV and XYZ of points on edges + + for ( int i = 0; i < myOrderedNodes.size(); ++i, ++iSub ) { - if ( !Apply( *face, theNodeIndexOnKeyPoint1, theReverse )) { - MESSAGE( "Failed on " << *face ); - continue; - } - myElements.push_back( *face ); + gp_XY& uv1 = keyUV[ i ]; + gp_XY& uv2 = ( i+1 < keyUV.size() ) ? keyUV[ i+1 ] : keyUV[ 0 ]; - // store computed points belonging to elements - list< list< int > >::iterator ll = myElemPointIDs.begin(); - for ( ; ll != myElemPointIDs.end(); ++ll ) + list< TPoint* > & ePoints = getShapePoints( iSub ); + ePoints.back()->myInitU = 1.0; + list< TPoint* >::const_iterator pIt = ++ePoints.begin(); + while ( *pIt != ePoints.back() ) { - myElemXYZIDs.push_back(); - list< int >& xyzIds = myElemXYZIDs.back(); - list< int >& pIds = *ll; - for ( list::iterator id = pIds.begin(); id != pIds.end(); id++ ) { - int pIndex = *id + ind1; - xyzIds.push_back( pIndex ); - myXYZ[ pIndex ] = myPoints[ *id ].myXYZ.XYZ(); - reverseConnectivity[ pIndex ].push_back( & xyzIds ); - } + TPoint* p = *pIt++; + p->myUV = uv1 * ( 1 - p->myInitU ) + uv2 * p->myInitU; + p->myXYZ = surface->Value( p->myUV.X(), p->myUV.Y() ); + if ( !loc.IsIdentity() ) + aTrsf.Transforms( p->myXYZ.ChangeCoord() ); } - // put points on links to linkPointIndListMap - int nbNodes = (*face)->NbNodes(), eID = nbNodes + 1; - for ( int i = 0; i < nbNodes; i++ ) - { + // collect on-edge points (excluding the last one) + edgesPoints->insert( edgesPoints->end(), ePoints.begin(), --ePoints.end()); + } + + // Compute UV and XYZ of in-face points + + // try to use a simple algo to compute UV + list< TPoint* > & fPoints = getShapePoints( iSub ); + bool isDeformed = false; + for ( pIt = fPoints.begin(); !isDeformed && pIt != fPoints.end(); pIt++ ) + if ( !compUVByIsoIntersection( edgesPointsList, (*pIt)->myInitUV, + (*pIt)->myUV, isDeformed )) { + MESSAGE("cant Apply(face)"); + return false; + } + // try to use a complex algo if it is a difficult case + if ( isDeformed && !compUVByElasticIsolines( edgesPointsList, fPoints )) + { + for ( ; pIt != fPoints.end(); pIt++ ) // continue with the simple algo + if ( !compUVByIsoIntersection( edgesPointsList, (*pIt)->myInitUV, + (*pIt)->myUV, isDeformed )) { + MESSAGE("cant Apply(face)"); + return false; + } + } + + for ( pIt = fPoints.begin(); pIt != fPoints.end(); pIt++ ) + { + TPoint * point = *pIt; + point->myXYZ = surface->Value( point->myUV.X(), point->myUV.Y() ); + if ( !loc.IsIdentity() ) + aTrsf.Transforms( point->myXYZ.ChangeCoord() ); + } + + myIsComputed = true; + + return setErrorCode( ERR_OK ); +} + +//======================================================================= +//function : undefinedXYZ +//purpose : +//======================================================================= + +static const gp_XYZ& undefinedXYZ() +{ + static gp_XYZ xyz( 1.e100, 0., 0. ); + return xyz; +} + +//======================================================================= +//function : isDefined +//purpose : +//======================================================================= + +inline static bool isDefined(const gp_XYZ& theXYZ) +{ + return theXYZ.X() < 1.e100; +} + +//======================================================================= +//function : Apply +//purpose : Compute nodes coordinates applying +// the loaded pattern to . The first key-point +// will be mapped into -th node +//======================================================================= + +bool SMESH_Pattern::Apply (SMESH_Mesh* theMesh, + std::set& theFaces, + const int theNodeIndexOnKeyPoint1, + const bool theReverse) +{ + MESSAGE(" ::Apply(set) " ); + + if ( !IsLoaded() ) { + MESSAGE( "Pattern not loaded" ); + return setErrorCode( ERR_APPL_NOT_LOADED ); + } + + // find points on edges, it fills myNbKeyPntInBoundary + if ( !findBoundaryPoints() ) + return false; + + // check that there are no holes in a pattern + if (myNbKeyPntInBoundary.size() > 1 ) { + return setErrorCode( ERR_APPL_BAD_NB_VERTICES ); + } + + myShape.Nullify(); + myXYZ.clear(); + myElemXYZIDs.clear(); + myXYZIdToNodeMap.clear(); + myElements.clear(); + myIdsOnBoundary.clear(); + myReverseConnectivity.clear(); + + myXYZ.resize( myPoints.size() * theFaces.size(), undefinedXYZ() ); + myElements.reserve( theFaces.size() ); + + // to find point index + map< TPoint*, int > pointIndex; + for ( int i = 0; i < myPoints.size(); i++ ) + pointIndex.insert( make_pair( & myPoints[ i ], i )); + + int ind1 = 0; // lowest point index for a face + + // meshed geometry + TopoDS_Shape shape; +// int shapeID = 0; +// SMESH_MeshEditor editor( theMesh ); + + // apply to each face in theFaces set + set::iterator face = theFaces.begin(); + for ( ; face != theFaces.end(); ++face ) + { +// int curShapeId = editor.FindShape( *face ); +// if ( curShapeId != shapeID ) { +// if ( curShapeId ) +// shape = theMesh->GetMeshDS()->IndexToShape( curShapeId ); +// else +// shape.Nullify(); +// shapeID = curShapeId; +// } + bool ok; + if ( shape.IsNull() ) + ok = Apply( *face, theNodeIndexOnKeyPoint1, theReverse ); + else + ok = Apply( theMesh, *face, shape, theNodeIndexOnKeyPoint1, theReverse ); + if ( !ok ) { + MESSAGE( "Failed on " << *face ); + continue; + } + myElements.push_back( *face ); + + // store computed points belonging to elements + list< TElemDef >::iterator ll = myElemPointIDs.begin(); + for ( ; ll != myElemPointIDs.end(); ++ll ) + { + myElemXYZIDs.push_back(TElemDef()); + TElemDef& xyzIds = myElemXYZIDs.back(); + TElemDef& pIds = *ll; + for ( TElemDef::iterator id = pIds.begin(); id != pIds.end(); id++ ) { + int pIndex = *id + ind1; + xyzIds.push_back( pIndex ); + myXYZ[ pIndex ] = myPoints[ *id ].myXYZ.XYZ(); + myReverseConnectivity[ pIndex ].push_back( & xyzIds ); + } + } + // put points on links to myIdsOnBoundary, + // they will be used to sew new elements on adjacent refined elements + int nbNodes = (*face)->NbNodes(), eID = nbNodes + 1; + for ( int i = 0; i < nbNodes; i++ ) + { + list< TPoint* > & linkPoints = getShapePoints( eID++ ); const SMDS_MeshNode* n1 = myOrderedNodes[ i ]; const SMDS_MeshNode* n2 = myOrderedNodes[ i + 1 == nbNodes ? 0 : i + 1 ]; - // make a link of node pointers - TLink link; - link.insert( n1 ); - link.insert( n2 ); - // add the link to the map - list< list< int > >& groups = linkPointIndListMap[ link ]; - groups.push_back(); - list< int >& indList = groups.back(); - list< TPoint* > & linkPoints = getShapePoints( eID++ ); + // make a link and a node set + TNodeSet linkSet, node1Set; + linkSet.insert( n1 ); + linkSet.insert( n2 ); + node1Set.insert( n1 ); list< TPoint* >::iterator p = linkPoints.begin(); - // map the first link point to n1 - myXYZIdToNodeMap[ pointIndex[ *p ] + ind1 ] = n1; + { + // map the first link point to n1 + int nId = pointIndex[ *p ] + ind1; + myXYZIdToNodeMap[ nId ] = n1; + list< list< int > >& groups = myIdsOnBoundary[ node1Set ]; + groups.push_back(list< int > ()); + groups.back().push_back( nId ); + } + // add the linkSet to the map + list< list< int > >& groups = myIdsOnBoundary[ linkSet ]; + groups.push_back(list< int > ()); + list< int >& indList = groups.back(); // add points to the map excluding the end points for ( p++; *p != linkPoints.back(); p++ ) indList.push_back( pointIndex[ *p ] + ind1 ); @@ -2777,8 +3012,6 @@ bool SMESH_Pattern::Apply (std::set theFaces, ind1 += myPoints.size(); } - mergePoints( linkPointIndListMap, reverseConnectivity ); - return !myElemXYZIDs.empty(); } @@ -2791,9 +3024,9 @@ bool SMESH_Pattern::Apply (std::set theFaces, // node. //======================================================================= -bool SMESH_Pattern::Apply (std::set theVolumes, - const int theNode000Index, - const int theNode001Index) +bool SMESH_Pattern::Apply (std::set & theVolumes, + const int theNode000Index, + const int theNode001Index) { MESSAGE(" ::Apply(set) " ); @@ -2811,10 +3044,13 @@ bool SMESH_Pattern::Apply (std::set theVolumes, return setErrorCode( ERR_APPL_BAD_NB_VERTICES ); } + myShape.Nullify(); myXYZ.clear(); myElemXYZIDs.clear(); myXYZIdToNodeMap.clear(); myElements.clear(); + myIdsOnBoundary.clear(); + myReverseConnectivity.clear(); myXYZ.resize( myPoints.size() * theVolumes.size(), undefinedXYZ() ); myElements.reserve( theVolumes.size() ); @@ -2824,10 +3060,6 @@ bool SMESH_Pattern::Apply (std::set theVolumes, for ( int i = 0; i < myPoints.size(); i++ ) pointIndex.insert( make_pair( & myPoints[ i ], i )); - // to merge nodes on edges and faces of the elements being refined - map< TNodeSet, list< list< int > > > subPointIndListMap; - map< int, list< list< int >* > > reverseConnectivity; - int ind1 = 0; // lowest point index for an element // apply to each element in theVolumes set @@ -2841,27 +3073,28 @@ bool SMESH_Pattern::Apply (std::set theVolumes, myElements.push_back( *vol ); // store computed points belonging to elements - list< list< int > >::iterator ll = myElemPointIDs.begin(); + list< TElemDef >::iterator ll = myElemPointIDs.begin(); for ( ; ll != myElemPointIDs.end(); ++ll ) { - myElemXYZIDs.push_back(); - list< int >& xyzIds = myElemXYZIDs.back(); - list< int >& pIds = *ll; - for ( list::iterator id = pIds.begin(); id != pIds.end(); id++ ) { + myElemXYZIDs.push_back(TElemDef()); + TElemDef& xyzIds = myElemXYZIDs.back(); + TElemDef& pIds = *ll; + for ( TElemDef::iterator id = pIds.begin(); id != pIds.end(); id++ ) { int pIndex = *id + ind1; xyzIds.push_back( pIndex ); myXYZ[ pIndex ] = myPoints[ *id ].myXYZ.XYZ(); - reverseConnectivity[ pIndex ].push_back( & xyzIds ); + myReverseConnectivity[ pIndex ].push_back( & xyzIds ); } } - // put points on edges and faces to subPointIndListMap + // put points on edges and faces to myIdsOnBoundary, + // they will be used to sew new elements on adjacent refined elements for ( int Id = SMESH_Block::ID_V000; Id <= SMESH_Block::ID_F1yz; Id++ ) { // make a set of sub-points TNodeSet subNodes; vector< int > subIDs; if ( SMESH_Block::IsVertexID( Id )) { - // use nodes of refined volumes for merge + subNodes.insert( myOrderedNodes[ Id - 1 ]); } else if ( SMESH_Block::IsEdgeID( Id )) { SMESH_Block::GetEdgeVertexIDs( Id, subIDs ); @@ -2878,23 +3111,20 @@ bool SMESH_Pattern::Apply (std::set theVolumes, subNodes.insert( myOrderedNodes[ subIDs.front() - 1 ]); subNodes.insert( myOrderedNodes[ subIDs.back() - 1 ]); } - list< list< int > >& groups = subPointIndListMap[ subNodes ]; - groups.push_back(); - list< int >& indList = groups.back(); // add points list< TPoint* > & points = getShapePoints( Id ); list< TPoint* >::iterator p = points.begin(); - if ( subNodes.empty() ) // vertex case - myXYZIdToNodeMap[ pointIndex[ *p ] + ind1 ] = myOrderedNodes[ Id - 1 ]; - else - for ( ; p != points.end(); p++ ) - indList.push_back( pointIndex[ *p ] + ind1 ); + list< list< int > >& groups = myIdsOnBoundary[ subNodes ]; + groups.push_back(list< int > ()); + list< int >& indList = groups.back(); + for ( ; p != points.end(); p++ ) + indList.push_back( pointIndex[ *p ] + ind1 ); + if ( subNodes.size() == 1 ) // vertex case + myXYZIdToNodeMap[ indList.back() ] = myOrderedNodes[ Id - 1 ]; } ind1 += myPoints.size(); } - mergePoints( subPointIndListMap, reverseConnectivity ); - return !myElemXYZIDs.empty(); } @@ -2909,7 +3139,7 @@ bool SMESH_Pattern::Load (SMESH_Mesh* theMesh, MESSAGE(" ::Load(volume) " ); Clear(); myIs2D = false; - SMESHDS_Mesh * aMeshDS = theMesh->GetMeshDS(); + SMESHDS_SubMesh * aSubMesh; // load shapes in myShapeIDMap SMESH_Block block; @@ -2922,7 +3152,7 @@ bool SMESH_Pattern::Load (SMESH_Mesh* theMesh, for ( shapeID = 1; shapeID <= myShapeIDMap.Extent(); shapeID++ ) { const TopoDS_Shape& S = myShapeIDMap( shapeID ); - SMESHDS_SubMesh * aSubMesh = aMeshDS->MeshElements( S ); + aSubMesh = getSubmeshWithElements( theMesh, S ); if ( aSubMesh ) nbNodes += aSubMesh->NbNodes(); } @@ -2935,14 +3165,14 @@ bool SMESH_Pattern::Load (SMESH_Mesh* theMesh, { const TopoDS_Shape& S = myShapeIDMap( shapeID ); list< TPoint* > & shapePoints = getShapePoints( shapeID ); - SMESHDS_SubMesh * aSubMesh = aMeshDS->MeshElements( S ); + aSubMesh = getSubmeshWithElements( theMesh, S ); if ( ! aSubMesh ) continue; SMDS_NodeIteratorPtr nIt = aSubMesh->GetNodes(); if ( !nIt->more() ) continue; // store a node and a point while ( nIt->more() ) { - const SMDS_MeshNode* node = static_cast( nIt->next() ); + const SMDS_MeshNode* node = smdsNode( nIt->next() ); nodePointIDMap.insert( make_pair( node, iPoint )); if ( block.IsVertexID( shapeID )) myKeyPointIDs.push_back( iPoint ); @@ -2977,10 +3207,9 @@ bool SMESH_Pattern::Load (SMESH_Mesh* theMesh, nIt = aSubMesh->GetNodes(); for ( ; nIt->more(); pIt++ ) { - const SMDS_MeshNode* node = - static_cast( nIt->next() ); + const SMDS_MeshNode* node = smdsNode( nIt->next() ); const SMDS_EdgePosition* epos = - static_cast(node->GetPosition().get()); + static_cast(node->GetPosition()); double u = ( epos->GetUParameter() - f ) / ( l - f ); (*pIt)->myInitXYZ.SetCoord( iCoord, isForward ? u : 1 - u ); } @@ -2999,14 +3228,14 @@ bool SMESH_Pattern::Load (SMESH_Mesh* theMesh, // load elements - SMESHDS_SubMesh * aSubMesh = aMeshDS->MeshElements( theBlock ); + aSubMesh = getSubmeshWithElements( theMesh, theBlock ); if ( aSubMesh ) { SMDS_ElemIteratorPtr elemIt = aSubMesh->GetElements(); while ( elemIt->more() ) { SMDS_ElemIteratorPtr nIt = elemIt->next()->nodesIterator(); - myElemPointIDs.push_back( list< int >() ); - list< int >& elemPoints = myElemPointIDs.back(); + myElemPointIDs.push_back( TElemDef() ); + TElemDef& elemPoints = myElemPointIDs.back(); while ( nIt->more() ) elemPoints.push_back( nodePointIDMap[ nIt->next() ]); } @@ -3017,6 +3246,32 @@ bool SMESH_Pattern::Load (SMESH_Mesh* theMesh, return setErrorCode( ERR_OK ); } +//======================================================================= +//function : getSubmeshWithElements +//purpose : return submesh containing elements bound to theBlock in theMesh +//======================================================================= + +SMESHDS_SubMesh * SMESH_Pattern::getSubmeshWithElements(SMESH_Mesh* theMesh, + const TopoDS_Shape& theShape) +{ + SMESHDS_SubMesh * aSubMesh = theMesh->GetMeshDS()->MeshElements( theShape ); + if ( aSubMesh && ( aSubMesh->GetElements()->more() || aSubMesh->GetNodes()->more() )) + return aSubMesh; + + if ( theShape.ShapeType() == TopAbs_SHELL ) + { + // look for submesh of VOLUME + TopTools_ListIteratorOfListOfShape it( theMesh->GetAncestors( theShape )); + for (; it.More(); it.Next()) { + aSubMesh = theMesh->GetMeshDS()->MeshElements( it.Value() ); + if ( aSubMesh && ( aSubMesh->GetElements()->more() || aSubMesh->GetNodes()->more() )) + return aSubMesh; + } + } + return 0; +} + + //======================================================================= //function : Apply //purpose : Compute nodes coordinates applying @@ -3090,7 +3345,7 @@ bool SMESH_Pattern::Apply (const SMDS_MeshVolume* theVolume, const int theNode000Index, const int theNode001Index) { - MESSAGE(" ::Apply(MeshVolume) " ); + //MESSAGE(" ::Apply(MeshVolume) " ); if (!findBoundaryPoints()) // bind ID to points return false; @@ -3128,53 +3383,509 @@ bool SMESH_Pattern::Apply (const SMDS_MeshVolume* theVolume, } //======================================================================= -//function : MakeMesh -//purpose : Create nodes and elements in using nodes -// coordinates computed by either of Apply...() methods +//function : mergePoints +//purpose : Merge XYZ on edges and/or faces. //======================================================================= -bool SMESH_Pattern::MakeMesh(SMESH_Mesh* theMesh) +void SMESH_Pattern::mergePoints (const bool uniteGroups) { - MESSAGE(" ::MakeMesh() " ); - if ( !myIsComputed ) - return setErrorCode( ERR_MAKEM_NOT_COMPUTED ); + map< TNodeSet, list< list< int > > >::iterator idListIt = myIdsOnBoundary.begin(); + for ( ; idListIt != myIdsOnBoundary.end(); idListIt++ ) + { + list >& groups = idListIt->second; + if ( groups.size() < 2 ) + continue; - SMESHDS_Mesh* aMeshDS = theMesh->GetMeshDS(); - SMESH_MeshEditor editor( theMesh ); + // find tolerance + const TNodeSet& nodes = idListIt->first; + double tol2 = 1.e-10; + if ( nodes.size() > 1 ) { + Bnd_Box box; + TNodeSet::const_iterator n = nodes.begin(); + for ( ; n != nodes.end(); ++n ) + box.Add( gp_Pnt( (*n)->X(), (*n)->Y(), (*n)->Z() )); + double x, y, z, X, Y, Z; + box.Get( x, y, z, X, Y, Z ); + gp_Pnt p( x, y, z ), P( X, Y, Z ); + tol2 = 1.e-4 * p.SquareDistance( P ); + } - // clear elements and nodes existing on myShape + // to unite groups on link + bool unite = ( uniteGroups && nodes.size() == 2 ); + map< double, int > distIndMap; + const SMDS_MeshNode* node = *nodes.begin(); + gp_Pnt P( node->X(), node->Y(), node->Z() ); - if ( !myShape.IsNull() ) + // compare points, replace indices + + list< int >::iterator ind1, ind2; + list< list< int > >::iterator grpIt1, grpIt2; + for ( grpIt1 = groups.begin(); grpIt1 != groups.end(); grpIt1++ ) + { + list< int >& indices1 = *grpIt1; + grpIt2 = grpIt1; + for ( grpIt2++; grpIt2 != groups.end(); grpIt2++ ) + { + list< int >& indices2 = *grpIt2; + for ( ind1 = indices1.begin(); ind1 != indices1.end(); ind1++ ) + { + gp_XYZ& p1 = myXYZ[ *ind1 ]; + ind2 = indices2.begin(); + while ( ind2 != indices2.end() ) + { + gp_XYZ& p2 = myXYZ[ *ind2 ]; + //MESSAGE("COMP: " << *ind1 << " " << *ind2 << " X: " << p2.X() << " tol2: " << tol2); + if ( ( p1 - p2 ).SquareModulus() <= tol2 ) + { + ASSERT( myReverseConnectivity.find( *ind2 ) != myReverseConnectivity.end() ); + list< TElemDef* > & elemXYZIDsList = myReverseConnectivity[ *ind2 ]; + list< TElemDef* >::iterator elemXYZIDs = elemXYZIDsList.begin(); + for ( ; elemXYZIDs != elemXYZIDsList.end(); elemXYZIDs++ ) + { + //MESSAGE( " Replace " << *ind2 << " with " << *ind1 ); + myXYZ[ *ind2 ] = undefinedXYZ(); + replace( (*elemXYZIDs)->begin(), (*elemXYZIDs)->end(), *ind2, *ind1 ); + } + ind2 = indices2.erase( ind2 ); + } + else + ind2++; + } + } + } + if ( unite ) { // sort indices using distIndMap + for ( ind1 = indices1.begin(); ind1 != indices1.end(); ind1++ ) + { + ASSERT( isDefined( myXYZ[ *ind1 ] )); + double dist = P.SquareDistance( myXYZ[ *ind1 ]); + distIndMap.insert( make_pair( dist, *ind1 )); + } + } + } + if ( unite ) { // put all sorted indices into the first group + list< int >& g = groups.front(); + g.clear(); + map< double, int >::iterator dist_ind = distIndMap.begin(); + for ( ; dist_ind != distIndMap.end(); dist_ind++ ) + g.push_back( dist_ind->second ); + } + } // loop on myIdsOnBoundary +} + +//======================================================================= +//function : makePolyElements +//purpose : prepare intermediate data to create Polygons and Polyhedrons +//======================================================================= + +void SMESH_Pattern:: + makePolyElements(const vector< const SMDS_MeshNode* >& theNodes, + const bool toCreatePolygons, + const bool toCreatePolyedrs) +{ + myPolyElemXYZIDs.clear(); + myPolyElems.clear(); + myPolyElems.reserve( myIdsOnBoundary.size() ); + + // make a set of refined elements + TIDSortedElemSet avoidSet, elemSet; + std::vector::iterator itv = myElements.begin(); + for(; itv!=myElements.end(); itv++) { + const SMDS_MeshElement* el = (*itv); + avoidSet.insert( el ); + } + //avoidSet.insert( myElements.begin(), myElements.end() ); + + map< TNodeSet, list< list< int > > >::iterator indListIt, nn_IdList; + + if ( toCreatePolygons ) { - SMESH_subMesh * aSubMesh = theMesh->GetSubMeshContaining( myShape ); - SMESHDS_SubMesh * aSubMeshDS = aMeshDS->MeshElements( myShape ); - if ( aSubMesh ) + int lastFreeId = myXYZ.size(); + + // loop on links of refined elements + indListIt = myIdsOnBoundary.begin(); + for ( ; indListIt != myIdsOnBoundary.end(); indListIt++ ) + { + const TNodeSet & linkNodes = indListIt->first; + if ( linkNodes.size() != 2 ) + continue; // skip face + const SMDS_MeshNode* n1 = * linkNodes.begin(); + const SMDS_MeshNode* n2 = * linkNodes.rbegin(); + + list >& idGroups = indListIt->second; // ids of nodes to build + if ( idGroups.empty() || idGroups.front().empty() ) + continue; + + // find not refined face having n1-n2 link + + while (true) + { + const SMDS_MeshElement* face = + SMESH_MeshEditor::FindFaceInSet( n1, n2, elemSet, avoidSet ); + if ( face ) + { + avoidSet.insert ( face ); + myPolyElems.push_back( face ); + + // some links of are split; + // make list of xyz for + myPolyElemXYZIDs.push_back(TElemDef()); + TElemDef & faceNodeIds = myPolyElemXYZIDs.back(); + // loop on links of a + SMDS_ElemIteratorPtr nIt = face->nodesIterator(); + int i = 0, nbNodes = face->NbNodes(); + vector nodes( nbNodes + 1 ); + while ( nIt->more() ) + nodes[ i++ ] = smdsNode( nIt->next() ); + nodes[ i ] = nodes[ 0 ]; + for ( i = 0; i < nbNodes; ++i ) + { + // look for point mapped on a link + TNodeSet faceLinkNodes; + faceLinkNodes.insert( nodes[ i ] ); + faceLinkNodes.insert( nodes[ i + 1 ] ); + if ( faceLinkNodes == linkNodes ) + nn_IdList = indListIt; + else + nn_IdList = myIdsOnBoundary.find( faceLinkNodes ); + // add face point ids + faceNodeIds.push_back( ++lastFreeId ); + myXYZIdToNodeMap.insert( make_pair( lastFreeId, nodes[ i ])); + if ( nn_IdList != myIdsOnBoundary.end() ) + { + // there are points mapped on a link + list< int >& mappedIds = nn_IdList->second.front(); + if ( isReversed( nodes[ i ], mappedIds )) + faceNodeIds.insert (faceNodeIds.end(),mappedIds.rbegin(), mappedIds.rend() ); + else + faceNodeIds.insert (faceNodeIds.end(),mappedIds.begin(), mappedIds.end() ); + } + } // loop on links of a + } // if ( face ) + else + break; + } // while (true) + + if ( myIs2D && idGroups.size() > 1 ) { + + // sew new elements on 2 refined elements sharing n1-n2 link + + list< int >& idsOnLink = idGroups.front(); + // temporarily add ids of link nodes to idsOnLink + bool rev = isReversed( n1, idsOnLink ); + for ( int i = 0; i < 2; ++i ) + { + TNodeSet nodeSet; + nodeSet.insert( i ? n2 : n1 ); + ASSERT( myIdsOnBoundary.find( nodeSet ) != myIdsOnBoundary.end() ); + list >& groups = myIdsOnBoundary[ nodeSet ]; + int nodeId = groups.front().front(); + bool append = i; + if ( rev ) append = !append; + if ( append ) + idsOnLink.push_back( nodeId ); + else + idsOnLink.push_front( nodeId ); + } + list< int >::iterator id = idsOnLink.begin(); + for ( ; id != idsOnLink.end(); ++id ) // loop on XYZ ids on a link + { + list< TElemDef* >& elemDefs = myReverseConnectivity[ *id ]; // elems sharing id + list< TElemDef* >::iterator pElemDef = elemDefs.begin(); + for ( ; pElemDef != elemDefs.end(); pElemDef++ ) // loop on elements sharing id + { + TElemDef* pIdList = *pElemDef; // ptr on list of ids making element up + // look for in element definition + TElemDef::iterator idDef = find( pIdList->begin(), pIdList->end(), *id ); + ASSERT ( idDef != pIdList->end() ); + // look for 2 neighbour ids of in element definition + for ( int prev = 0; prev < 2; ++prev ) { + TElemDef::iterator idDef2 = idDef; + if ( prev ) + idDef2 = ( idDef2 == pIdList->begin() ) ? --pIdList->end() : --idDef2; + else + idDef2 = ( ++idDef2 == pIdList->end() ) ? pIdList->begin() : idDef2; + // look for idDef2 on a link starting from id + list< int >::iterator id2 = find( id, idsOnLink.end(), *idDef2 ); + if ( id2 != idsOnLink.end() && id != --id2 ) { // found not next to id + // insert ids located on link between and + // into the element definition between idDef and idDef2 + if ( prev ) + for ( ; id2 != id; --id2 ) + pIdList->insert( idDef, *id2 ); + else { + list< int >::iterator id1 = id; + for ( ++id1, ++id2; id1 != id2; ++id1 ) + pIdList->insert( idDef2, *id1 ); + } + } + } + } + } + // remove ids of link nodes + idsOnLink.pop_front(); + idsOnLink.pop_back(); + } + } // loop on myIdsOnBoundary + } // if ( toCreatePolygons ) + + if ( toCreatePolyedrs ) + { + // check volumes adjacent to the refined elements + SMDS_VolumeTool volTool; + vector::iterator refinedElem = myElements.begin(); + for ( ; refinedElem != myElements.end(); ++refinedElem ) + { + // loop on nodes of refinedElem + SMDS_ElemIteratorPtr nIt = (*refinedElem)->nodesIterator(); + while ( nIt->more() ) { + const SMDS_MeshNode* node = smdsNode( nIt->next() ); + // loop on inverse elements of node + SMDS_ElemIteratorPtr eIt = node->GetInverseElementIterator(); + while ( eIt->more() ) + { + const SMDS_MeshElement* elem = eIt->next(); + if ( !volTool.Set( elem ) || !avoidSet.insert( elem ).second ) + continue; // skip faces or refined elements + // add polyhedron definition + myPolyhedronQuantities.push_back(vector ()); + myPolyElemXYZIDs.push_back(TElemDef()); + vector& quantity = myPolyhedronQuantities.back(); + TElemDef & elemDef = myPolyElemXYZIDs.back(); + // get definitions of new elements on volume faces + bool makePoly = false; + for ( int iF = 0; iF < volTool.NbFaces(); ++iF ) + { + if ( getFacesDefinition(volTool.GetFaceNodes( iF ), + volTool.NbFaceNodes( iF ), + theNodes, elemDef, quantity)) + makePoly = true; + } + if ( makePoly ) + myPolyElems.push_back( elem ); + else { + myPolyhedronQuantities.pop_back(); + myPolyElemXYZIDs.pop_back(); + } + } + } + } + } +} + +//======================================================================= +//function : getFacesDefinition +//purpose : return faces definition for a volume face defined by theBndNodes +//======================================================================= + +bool SMESH_Pattern:: + getFacesDefinition(const SMDS_MeshNode** theBndNodes, + const int theNbBndNodes, + const vector< const SMDS_MeshNode* >& theNodes, + list< int >& theFaceDefs, + vector& theQuantity) +{ + bool makePoly = false; +// cout << "FROM FACE NODES: " < bndNodeSet; + for ( int i = 0; i < theNbBndNodes; ++i ) + bndNodeSet.insert( theBndNodes[ i ]); + + map< TNodeSet, list< list< int > > >::iterator nn_IdList; + + // make a set of all nodes on a face + set< int > ids; + if ( !myIs2D ) { // for 2D, merge only edges + nn_IdList = myIdsOnBoundary.find( bndNodeSet ); + if ( nn_IdList != myIdsOnBoundary.end() ) { + makePoly = true; + list< int > & faceIds = nn_IdList->second.front(); + ids.insert( faceIds.begin(), faceIds.end() ); + } + } + //bool hasIdsInFace = !ids.empty(); + + // add ids on links and bnd nodes + int lastFreeId = Max( myXYZIdToNodeMap.rbegin()->first, theNodes.size() ); + TElemDef faceDef; // definition for the case if there is no new adjacent volumes + for ( int iN = 0; iN < theNbBndNodes; ++iN ) + { + // add id of iN-th bnd node + TNodeSet nSet; + nSet.insert( theBndNodes[ iN ] ); + nn_IdList = myIdsOnBoundary.find( nSet ); + int bndId = ++lastFreeId; + if ( nn_IdList != myIdsOnBoundary.end() ) { + bndId = nn_IdList->second.front().front(); + ids.insert( bndId ); + } + else + myXYZIdToNodeMap.insert( make_pair( bndId, theBndNodes[ iN ] )); + faceDef.push_back( bndId ); + // add ids on a link + TNodeSet linkNodes; + linkNodes.insert( theBndNodes[ iN ]); + linkNodes.insert( theBndNodes[ iN + 1 == theNbBndNodes ? 0 : iN + 1 ]); + nn_IdList = myIdsOnBoundary.find( linkNodes ); + if ( nn_IdList != myIdsOnBoundary.end() ) { + makePoly = true; + list< int > & linkIds = nn_IdList->second.front(); + ids.insert( linkIds.begin(), linkIds.end() ); + if ( isReversed( theBndNodes[ iN ], linkIds )) + faceDef.insert( faceDef.end(), linkIds.begin(), linkIds.end() ); + else + faceDef.insert( faceDef.end(), linkIds.rbegin(), linkIds.rend() ); + } + } + + // find faces definition of new volumes + + bool defsAdded = false; + if ( !myIs2D ) { // for 2D, merge only edges + SMDS_VolumeTool vol; + set< TElemDef* > checkedVolDefs; + set< int >::iterator id = ids.begin(); + for ( ; id != ids.end(); ++id ) + { + // definitions of volumes sharing id + list< TElemDef* >& defList = myReverseConnectivity[ *id ]; + ASSERT( !defList.empty() ); + // loop on volume definitions + list< TElemDef* >::iterator pIdList = defList.begin(); + for ( ; pIdList != defList.end(); ++pIdList) + { + if ( !checkedVolDefs.insert( *pIdList ).second ) + continue; // skip already checked volume definition + vector< int > idVec; + idVec.reserve( (*pIdList)->size() ); + idVec.insert( idVec.begin(), (*pIdList)->begin(), (*pIdList)->end() ); + // loop on face defs of a volume + SMDS_VolumeTool::VolumeType volType = vol.GetType( idVec.size() ); + if ( volType == SMDS_VolumeTool::UNKNOWN ) + continue; + int nbFaces = vol.NbFaces( volType ); + for ( int iF = 0; iF < nbFaces; ++iF ) + { + const int* nodeInds = vol.GetFaceNodesIndices( volType, iF, true ); + int iN, nbN = vol.NbFaceNodes( volType, iF ); + // check if all nodes of a faces are in + bool all = true; + for ( iN = 0; iN < nbN && all; ++iN ) { + int nodeId = idVec[ nodeInds[ iN ]]; + all = ( ids.find( nodeId ) != ids.end() ); + } + if ( all ) { + // store a face definition + for ( iN = 0; iN < nbN; ++iN ) { + theFaceDefs.push_back( idVec[ nodeInds[ iN ]]); + } + theQuantity.push_back( nbN ); + defsAdded = true; + } + } + } + } + } + if ( !defsAdded ) { + theQuantity.push_back( faceDef.size() ); + theFaceDefs.splice( theFaceDefs.end(), faceDef, faceDef.begin(), faceDef.end() ); + } + + return makePoly; +} + +//======================================================================= +//function : clearSubMesh +//purpose : +//======================================================================= + +static bool clearSubMesh( SMESH_Mesh* theMesh, + const TopoDS_Shape& theShape) +{ + bool removed = false; + if ( SMESH_subMesh * aSubMesh = theMesh->GetSubMeshContaining( theShape )) + { + removed = !aSubMesh->IsEmpty(); + if ( removed ) aSubMesh->ComputeStateEngine( SMESH_subMesh::CLEAN ); - else if ( aSubMeshDS ) + } + else { + SMESHDS_Mesh* aMeshDS = theMesh->GetMeshDS(); + if ( SMESHDS_SubMesh* aSubMeshDS = aMeshDS->MeshElements( theShape )) { SMDS_ElemIteratorPtr eIt = aSubMeshDS->GetElements(); + removed = eIt->more(); while ( eIt->more() ) aMeshDS->RemoveElement( eIt->next() ); SMDS_NodeIteratorPtr nIt = aSubMeshDS->GetNodes(); + removed = removed || nIt->more(); while ( nIt->more() ) - aMeshDS->RemoveNode( static_cast( nIt->next() )); + aMeshDS->RemoveNode( smdsNode( nIt->next() )); + } + } + return removed; +} + +//======================================================================= +//function : clearMesh +//purpose : clear mesh elements existing on myShape in theMesh +//======================================================================= + +void SMESH_Pattern::clearMesh(SMESH_Mesh* theMesh) const +{ + + if ( !myShape.IsNull() ) + { + if ( !clearSubMesh( theMesh, myShape ) && !myIs2D ) { // myShape is SHELL but volumes may be bound to SOLID + TopTools_ListIteratorOfListOfShape it( theMesh->GetAncestors( myShape )); + for (; it.More() && it.Value().ShapeType() == TopAbs_SOLID; it.Next()) + { + clearSubMesh( theMesh, it.Value() ); + } } } +} + +//======================================================================= +//function : MakeMesh +//purpose : Create nodes and elements in using nodes +// coordinates computed by either of Apply...() methods +// WARNING : StdMeshers_Projection_... relies on MakeMesh() behavior: that +// it does not care of nodes and elements already existing on +// subshapes. DO NOT MERGE them or modify also StdMeshers_Projection_.. +//======================================================================= + +bool SMESH_Pattern::MakeMesh(SMESH_Mesh* theMesh, + const bool toCreatePolygons, + const bool toCreatePolyedrs) +{ + MESSAGE(" ::MakeMesh() " ); + if ( !myIsComputed ) + return setErrorCode( ERR_MAKEM_NOT_COMPUTED ); + + mergePoints( toCreatePolygons ); + + SMESHDS_Mesh* aMeshDS = theMesh->GetMeshDS(); + + // clear elements and nodes existing on myShape + clearMesh(theMesh); bool onMeshElements = ( !myElements.empty() ); - // loop on sub-shapes of myShape: create nodes and build point-node map + // Create missing nodes - vector< const SMDS_MeshNode* > nodesVector; - map< TPoint*, const SMDS_MeshNode* > pointNodeMap; + vector< const SMDS_MeshNode* > nodesVector; // i-th point/xyz -> node if ( onMeshElements ) { - nodesVector.resize( myXYZ.size() ); + nodesVector.resize( Max( myXYZ.size(), myXYZIdToNodeMap.rbegin()->first ), 0 ); + map< int, const SMDS_MeshNode*>::iterator i_node = myXYZIdToNodeMap.begin(); + for ( ; i_node != myXYZIdToNodeMap.end(); i_node++ ) { + nodesVector[ i_node->first ] = i_node->second; + } for ( int i = 0; i < myXYZ.size(); ++i ) { - map< int, const SMDS_MeshNode*>::iterator idNode = myXYZIdToNodeMap.find( i ); - if ( idNode != myXYZIdToNodeMap.end() ) - nodesVector[ i ] = idNode->second; - else if ( isDefined( myXYZ[ i ] )) + if ( !nodesVector[ i ] && isDefined( myXYZ[ i ] ) ) nodesVector[ i ] = aMeshDS->AddNode (myXYZ[ i ].X(), myXYZ[ i ].Y(), myXYZ[ i ].Z()); @@ -3182,46 +3893,48 @@ bool SMESH_Pattern::MakeMesh(SMESH_Mesh* theMesh) } else { + nodesVector.resize( myPoints.size(), 0 ); + + // to find point index + map< TPoint*, int > pointIndex; + for ( int i = 0; i < myPoints.size(); i++ ) + pointIndex.insert( make_pair( & myPoints[ i ], i )); + + // loop on sub-shapes of myShape: create nodes map< int, list< TPoint* > >::iterator idPointIt = myShapeIDToPointsMap.begin(); for ( ; idPointIt != myShapeIDToPointsMap.end(); idPointIt++ ) { TopoDS_Shape S; - SMESHDS_SubMesh * subMeshDS = 0; + //SMESHDS_SubMesh * subMeshDS = 0; if ( !myShapeIDMap.IsEmpty() ) { S = myShapeIDMap( idPointIt->first ); - subMeshDS = aMeshDS->MeshElements( S ); + //subMeshDS = aMeshDS->MeshElements( S ); } list< TPoint* > & points = idPointIt->second; list< TPoint* >::iterator pIt = points.begin(); for ( ; pIt != points.end(); pIt++ ) { TPoint* point = *pIt; - if ( pointNodeMap.find( point ) != pointNodeMap.end() ) + int pIndex = pointIndex[ point ]; + if ( nodesVector [ pIndex ] ) continue; SMDS_MeshNode* node = aMeshDS->AddNode (point->myXYZ.X(), point->myXYZ.Y(), point->myXYZ.Z()); - pointNodeMap.insert( make_pair( point, node )); - if ( subMeshDS ) { + nodesVector [ pIndex ] = node; + + if ( !S.IsNull() /*subMeshDS*/ ) { + // !!!!! do not merge new nodes with ones existing on submeshes (see method comment) switch ( S.ShapeType() ) { case TopAbs_VERTEX: { - aMeshDS->SetNodeOnVertex( node, TopoDS::Vertex( S )); - break; + aMeshDS->SetNodeOnVertex( node, TopoDS::Vertex( S )); break; } case TopAbs_EDGE: { - aMeshDS->SetNodeOnEdge( node, TopoDS::Edge( S )); - SMDS_EdgePosition* epos = - dynamic_cast(node->GetPosition().get()); - epos->SetUParameter( point->myU ); - break; + aMeshDS->SetNodeOnEdge( node, TopoDS::Edge( S ), point->myU ); break; } case TopAbs_FACE: { - aMeshDS->SetNodeOnFace( node, TopoDS::Face( S )); - SMDS_FacePosition* pos = - dynamic_cast(node->GetPosition().get()); - pos->SetUParameter( point->myUV.X() ); - pos->SetVParameter( point->myUV.Y() ); - break; + aMeshDS->SetNodeOnFace( node, TopoDS::Face( S ), + point->myUV.X(), point->myUV.Y() ); break; } default: aMeshDS->SetNodeInVolume( node, TopoDS::Shell( S )); @@ -3230,62 +3943,136 @@ bool SMESH_Pattern::MakeMesh(SMESH_Mesh* theMesh) } } } - + // create elements - // shapes and groups myElements are on + if ( onMeshElements ) + { + // prepare data to create poly elements + makePolyElements( nodesVector, toCreatePolygons, toCreatePolyedrs ); + + // refine elements + createElements( theMesh, nodesVector, myElemXYZIDs, myElements ); + // sew old and new elements + createElements( theMesh, nodesVector, myPolyElemXYZIDs, myPolyElems ); + } + else + { + createElements( theMesh, nodesVector, myElemPointIDs, myElements ); + } + +// const map& sm = aMeshDS->SubMeshes(); +// map::const_iterator i_sm = sm.begin(); +// for ( ; i_sm != sm.end(); i_sm++ ) +// { +// cout << " SM " << i_sm->first << " "; +// TopAbs::Print( aMeshDS->IndexToShape( i_sm->first ).ShapeType(), cout)<< " "; +// //SMDS_ElemIteratorPtr GetElements(); +// SMDS_NodeIteratorPtr nit = i_sm->second->GetNodes(); +// while ( nit->more() ) +// cout << nit->next()->GetID() << " "; +// cout << endl; +// } + return setErrorCode( ERR_OK ); +} + +//======================================================================= +//function : createElements +//purpose : add elements to the mesh +//======================================================================= + +void SMESH_Pattern::createElements(SMESH_Mesh* theMesh, + const vector& theNodesVector, + const list< TElemDef > & theElemNodeIDs, + const vector& theElements) +{ + SMESHDS_Mesh* aMeshDS = theMesh->GetMeshDS(); + SMESH_MeshEditor editor( theMesh ); + + bool onMeshElements = !theElements.empty(); + + // shapes and groups theElements are on vector< int > shapeIDs; vector< list< SMESHDS_Group* > > groups; + set< const SMDS_MeshNode* > shellNodes; if ( onMeshElements ) { - shapeIDs.resize( myElements.size() ); - groups.resize( myElements.size() ); + shapeIDs.resize( theElements.size() ); + groups.resize( theElements.size() ); const set& allGroups = aMeshDS->GetGroups(); set::const_iterator grIt; - for ( int i = 0; i < myElements.size(); i++ ) + for ( int i = 0; i < theElements.size(); i++ ) { - shapeIDs[ i ] = editor.FindShape( myElements[ i ] ); + shapeIDs[ i ] = editor.FindShape( theElements[ i ] ); for ( grIt = allGroups.begin(); grIt != allGroups.end(); grIt++ ) { SMESHDS_Group* group = dynamic_cast( *grIt ); - if ( group && group->SMDSGroup().Contains( myElements[ i ] )) + if ( group && group->SMDSGroup().Contains( theElements[ i ] )) groups[ i ].push_back( group ); } } + // get all nodes bound to shells because their SpacePosition is not set + // by SMESHDS_Mesh::SetNodeInVolume() + TopoDS_Shape aMainShape = aMeshDS->ShapeToMesh(); + if ( !aMainShape.IsNull() ) { + TopExp_Explorer shellExp( aMainShape, TopAbs_SHELL ); + for ( ; shellExp.More(); shellExp.Next() ) + { + SMESHDS_SubMesh * sm = aMeshDS->MeshElements( shellExp.Current() ); + if ( sm ) { + SMDS_NodeIteratorPtr nIt = sm->GetNodes(); + while ( nIt->more() ) + shellNodes.insert( nIt->next() ); + } + } + } } - int nbElems = myElemPointIDs.size(); // nb elements in a pattern + // nb new elements per a refined element + int nbNewElemsPerOld = 1; + if ( onMeshElements ) + nbNewElemsPerOld = theElemNodeIDs.size() / theElements.size(); - list >::iterator epIt, epEnd; - if ( onMeshElements ) { - epIt = myElemXYZIDs.begin(); - epEnd = myElemXYZIDs.end(); - } - else { - epIt = myElemPointIDs.begin(); - epEnd = myElemPointIDs.end(); - } - for ( int iElem = 0; epIt != epEnd; epIt++, iElem++ ) + bool is2d = myIs2D; + + list< TElemDef >::const_iterator enIt = theElemNodeIDs.begin(); + list< vector >::iterator quantity = myPolyhedronQuantities.begin(); + for ( int iElem = 0; enIt != theElemNodeIDs.end(); enIt++, iElem++ ) { - list< int > & elemPoints = *epIt; + const TElemDef & elemNodeInd = *enIt; // retrieve nodes - const SMDS_MeshNode* nodes[ 8 ]; - list< int >::iterator iIt = elemPoints.begin(); + vector< const SMDS_MeshNode* > nodes( elemNodeInd.size() ); + TElemDef::const_iterator id = elemNodeInd.begin(); int nbNodes; - for ( nbNodes = 0; iIt != elemPoints.end(); iIt++ ) { - if ( onMeshElements ) - nodes[ nbNodes++ ] = nodesVector[ *iIt ]; + for ( nbNodes = 0; id != elemNodeInd.end(); id++ ) { + if ( *id < theNodesVector.size() ) + nodes[ nbNodes++ ] = theNodesVector[ *id ]; else - nodes[ nbNodes++ ] = pointNodeMap[ & myPoints[ *iIt ]]; + nodes[ nbNodes++ ] = myXYZIdToNodeMap[ *id ]; + } + // dim of refined elem + int elemIndex = iElem / nbNewElemsPerOld; // refined element index + if ( onMeshElements ) { + is2d = ( theElements[ elemIndex ]->GetType() == SMDSAbs_Face ); } // add an element const SMDS_MeshElement* elem = 0; - if ( myIs2D ) { + if ( is2d ) { switch ( nbNodes ) { case 3: elem = aMeshDS->AddFace( nodes[0], nodes[1], nodes[2] ); break; case 4: elem = aMeshDS->AddFace( nodes[0], nodes[1], nodes[2], nodes[3] ); break; + case 6: + if ( !onMeshElements ) {// create a quadratic face + elem = aMeshDS->AddFace (nodes[0], nodes[1], nodes[2], nodes[3], + nodes[4], nodes[5] ); break; + } // else do not break but create a polygon + case 8: + if ( !onMeshElements ) {// create a quadratic face + elem = aMeshDS->AddFace (nodes[0], nodes[1], nodes[2], nodes[3], + nodes[4], nodes[5], nodes[6], nodes[7] ); break; + } // else do not break but create a polygon default: - ASSERT( nbNodes < 8 ); + elem = aMeshDS->AddPolygonalFace( nodes ); } } else { @@ -3302,51 +4089,96 @@ bool SMESH_Pattern::MakeMesh(SMESH_Mesh* theMesh) elem = aMeshDS->AddVolume (nodes[0], nodes[1], nodes[2], nodes[3], nodes[4], nodes[5], nodes[6], nodes[7] ); break; default: - ASSERT( nbNodes < 8 ); + elem = aMeshDS->AddPolyhedralVolume( nodes, *quantity++ ); } } // set element on a shape - if ( elem && onMeshElements ) + if ( elem && onMeshElements ) // applied to mesh elements { - int elemIndex = iElem / nbElems; - if ( shapeIDs[ elemIndex ] > 0 ) - aMeshDS->SetMeshElementOnShape( elem, shapeIDs[ elemIndex ] ); + int shapeID = shapeIDs[ elemIndex ]; + if ( shapeID > 0 ) { + aMeshDS->SetMeshElementOnShape( elem, shapeID ); + // set nodes on a shape + TopoDS_Shape S = aMeshDS->IndexToShape( shapeID ); + if ( S.ShapeType() == TopAbs_SOLID ) { + TopoDS_Iterator shellIt( S ); + if ( shellIt.More() ) + shapeID = aMeshDS->ShapeToIndex( shellIt.Value() ); + } + SMDS_ElemIteratorPtr noIt = elem->nodesIterator(); + while ( noIt->more() ) { + SMDS_MeshNode* node = const_cast(smdsNode( noIt->next() )); + if (!node->GetPosition()->GetShapeId() && + shellNodes.find( node ) == shellNodes.end() ) { + if ( S.ShapeType() == TopAbs_FACE ) + aMeshDS->SetNodeOnFace( node, shapeID ); + else { + aMeshDS->SetNodeInVolume( node, shapeID ); + shellNodes.insert( node ); + } + } + } + } // add elem in groups list< SMESHDS_Group* >::iterator g = groups[ elemIndex ].begin(); for ( ; g != groups[ elemIndex ].end(); ++g ) (*g)->SMDSGroup().Add( elem ); } - if ( elem && !myShape.IsNull() ) + if ( elem && !myShape.IsNull() ) // applied to shape aMeshDS->SetMeshElementOnShape( elem, myShape ); } - // make that SMESH_subMesh::_computeState = COMPUTE_OK + // make that SMESH_subMesh::_computeState == COMPUTE_OK // so that operations with hypotheses will erase the mesh being built SMESH_subMesh * subMesh; if ( !myShape.IsNull() ) { - subMesh = theMesh->GetSubMeshContaining( myShape ); + subMesh = theMesh->GetSubMesh( myShape ); if ( subMesh ) subMesh->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE ); } if ( onMeshElements ) { list< int > elemIDs; - for ( int i = 0; i < myElements.size(); i++ ) + for ( int i = 0; i < theElements.size(); i++ ) { - int shapeID = shapeIDs[ i ]; - if ( shapeID > 0 ) { - TopoDS_Shape S = aMeshDS->IndexToShape( shapeID ); - subMesh = theMesh->GetSubMeshContaining( S ); - if ( subMesh ) - subMesh->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE ); - } - elemIDs.push_back( myElements[ i ]->GetID() ); + subMesh = theMesh->GetSubMeshContaining( shapeIDs[ i ] ); + if ( subMesh ) + subMesh->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE ); + + elemIDs.push_back( theElements[ i ]->GetID() ); } - // remove refined elements and their nodes + // remove refined elements editor.Remove( elemIDs, false ); } +} - return setErrorCode( ERR_OK ); +//======================================================================= +//function : isReversed +//purpose : check xyz ids order in theIdsList taking into account +// theFirstNode on a link +//======================================================================= + +bool SMESH_Pattern::isReversed(const SMDS_MeshNode* theFirstNode, + const list< int >& theIdsList) const +{ + if ( theIdsList.size() < 2 ) + return false; + + gp_Pnt Pf ( theFirstNode->X(), theFirstNode->Y(), theFirstNode->Z() ); + gp_Pnt P[2]; + list::const_iterator id = theIdsList.begin(); + for ( int i = 0; i < 2; ++i, ++id ) { + if ( *id < myXYZ.size() ) + P[ i ] = myXYZ[ *id ]; + else { + map< int, const SMDS_MeshNode*>::const_iterator i_n; + i_n = myXYZIdToNodeMap.find( *id ); + ASSERT( i_n != myXYZIdToNodeMap.end() ); + const SMDS_MeshNode* n = i_n->second; + P[ i ].SetCoord( n->X(), n->Y(), n->Z() ); + } + } + return Pf.SquareDistance( P[ 1 ] ) < Pf.SquareDistance( P[ 0 ] ); } @@ -3371,7 +4203,7 @@ void SMESH_Pattern::arrangeBoundaries (list< list< TPoint* > >& boundaryList) if ( nbBoundaries > 2 ) { // move boundaries in tmp list - list< list< TPoint* > > tmpList; + list< list< TPoint* > > tmpList; tmpList.splice( tmpList.begin(), boundaryList, boundaryList.begin(), boundaryList.end()); // make a map nb-key-points to boundary-position-in-tmpList, // boundary-positions get ordered in it @@ -3414,7 +4246,7 @@ void SMESH_Pattern::arrangeBoundaries (list< list< TPoint* > >& boundaryList) boundaryList.splice( boundaryList.begin(), boundaryList, outerBndPos, ++outerBndPos ); } // if nbBoundaries > 1 - + // Check boundaries orientation and re-fill myKeyPointIDs set< TPoint* > keyPointSet; @@ -3508,6 +4340,8 @@ bool SMESH_Pattern::findBoundaryPoints() MESSAGE(" findBoundaryPoints() "); + myNbKeyPntInBoundary.clear(); + if ( myIs2D ) { set< TPoint* > pointsInElems; @@ -3517,11 +4351,11 @@ bool SMESH_Pattern::findBoundaryPoints() typedef pair< TPoint*, TPoint*> TLink; set< TLink > linkSet; - list >::iterator epIt = myElemPointIDs.begin(); + list::iterator epIt = myElemPointIDs.begin(); for ( ; epIt != myElemPointIDs.end(); epIt++ ) { - list< int > & elemPoints = *epIt; - list< int >::iterator pIt = elemPoints.begin(); + TElemDef & elemPoints = *epIt; + TElemDef::iterator pIt = elemPoints.begin(); int prevP = elemPoints.back(); for ( ; pIt != elemPoints.end(); pIt++ ) { TPoint* p1 = & myPoints[ prevP ]; @@ -3683,7 +4517,7 @@ bool SMESH_Pattern::findBoundaryPoints() getShapePoints( shapeID ).push_back( point ); // detect key-points if ( SMESH_Block::IsVertexID( shapeID )) - myKeyPointIDs.push_back( i ); + myKeyPointIDs.push_back( i ); } } @@ -3731,21 +4565,30 @@ bool SMESH_Pattern::setShapeToMesh(const TopoDS_Shape& theShape) // check if a face is closed int nbNodeOnSeamEdge = 0; if ( myIs2D ) { + TopTools_MapOfShape seamVertices; TopoDS_Face face = TopoDS::Face( theShape ); TopExp_Explorer eExp( theShape, TopAbs_EDGE ); - for ( ; eExp.More() && nbNodeOnSeamEdge == 0; eExp.Next() ) - if ( BRep_Tool::IsClosed( TopoDS::Edge( eExp.Current() ), face )) - nbNodeOnSeamEdge = 2; + for ( ; eExp.More() && nbNodeOnSeamEdge == 0; eExp.Next() ) { + const TopoDS_Edge& ee = TopoDS::Edge(eExp.Current()); + if ( BRep_Tool::IsClosed(ee, face) ) { + // seam edge and vertices encounter twice in theFace + if ( !seamVertices.Add( TopExp::FirstVertex( ee ))) nbNodeOnSeamEdge++; + if ( !seamVertices.Add( TopExp::LastVertex( ee ))) nbNodeOnSeamEdge++; + } + } } - + // check nb of vertices TopTools_IndexedMapOfShape vMap; TopExp::MapShapes( theShape, TopAbs_VERTEX, vMap ); if ( vMap.Extent() + nbNodeOnSeamEdge != myKeyPointIDs.size() ) { - MESSAGE( myKeyPointIDs.size() << " != " << vMap.Extent() ); + MESSAGE( myKeyPointIDs.size() + nbNodeOnSeamEdge << " != " << vMap.Extent() ); return setErrorCode( ERR_APPL_BAD_NB_VERTICES ); } + myElements.clear(); // not refine elements + myElemXYZIDs.clear(); + myShapeIDMap.Clear(); myShape = theShape; return true; @@ -3768,9 +4611,13 @@ bool SMESH_Pattern::GetMappedPoints ( list< const gp_XYZ * > & thePoints ) const thePoints.push_back( & (*pVecIt).myXYZ.XYZ() ); } else { // applied to mesh elements + const gp_XYZ * definedXYZ = & myPoints[ myKeyPointIDs.front() ].myXYZ.XYZ(); vector::const_iterator xyz = myXYZ.begin(); for ( ; xyz != myXYZ.end(); ++xyz ) - thePoints.push_back( & (*xyz) ); + if ( !isDefined( *xyz )) + thePoints.push_back( definedXYZ ); + else + thePoints.push_back( & (*xyz) ); } return !thePoints.empty(); } @@ -3832,7 +4679,7 @@ void SMESH_Pattern::DumpPoints() const #ifdef _DEBUG_ vector< TPoint >::const_iterator pVecIt = myPoints.begin(); for ( int i = 0; pVecIt != myPoints.end(); pVecIt++, i++ ) - cout << i << ": " << *pVecIt; + MESSAGE_ADD ( std::endl << i << ": " << *pVecIt ); #endif } @@ -3872,6 +4719,6 @@ ostream & operator <<(ostream & OS, const SMESH_Pattern::TPoint& p) OS << " uv( " << xy.X() << " " << xy.Y() << " )"; u = p.myU; OS << " u( " << u << " ))" << endl; - + return OS; }