X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FSMESH%2FSMESH_Pattern.cxx;h=493799d3551c53b0f4a02aa8338318eb5c75391c;hp=7f20a5cad0f97a57c5ed29fa2597704481299df7;hb=d5f8c39955a2ef0d57238a360a746268994e549a;hpb=e4737e85f0da6d3f90fd08f6be1c2825195fe16f diff --git a/src/SMESH/SMESH_Pattern.cxx b/src/SMESH/SMESH_Pattern.cxx index 7f20a5cad..493799d35 100644 --- a/src/SMESH/SMESH_Pattern.cxx +++ b/src/SMESH/SMESH_Pattern.cxx @@ -1,16 +1,38 @@ -// File : SMESH_Pattern.cxx -// Created : Thu Aug 5 11:09:29 2004 +// 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 + +// File : SMESH_Pattern.hxx +// Created : Mon Aug 2 10:30:00 2004 // Author : Edward AGAPOV (eap) -// Copyright : Open CASCADE - #include "SMESH_Pattern.hxx" -#include #include #include #include +#include +#include +#include +#include +#include #include +#include #include #include #include @@ -24,29 +46,27 @@ #include #include #include -#include #include -#include +#include #include #include #include #include #include -#include -#include -#include -#include -#include -#include -#include #include "SMDS_EdgePosition.hxx" #include "SMDS_FacePosition.hxx" #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_subMesh.hxx" #include "utilities.h" @@ -54,8 +74,6 @@ using namespace std; typedef map< const SMDS_MeshElement*, int > TNodePointIDMap; -#define SQRT_FUNC 1 - //======================================================================= //function : SMESH_Pattern //purpose : @@ -281,8 +299,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 ); @@ -363,11 +381,11 @@ 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; @@ -391,9 +409,6 @@ template struct TSizeCmp { template void sortBySize( list< list < T > > & theListOfList ) { if ( theListOfList.size() > 2 ) { - // keep the car - //list < T > & aFront = theListOfList.front(); - // sort the whole list TSizeCmp< T > SizeCmp; theListOfList.sort( SizeCmp ); } @@ -548,18 +563,20 @@ bool SMESH_Pattern::Load (SMESH_Mesh* theMesh, return setErrorCode( ERR_LOAD_EMPTY_SUBMESH ); } + TopoDS_Face face = TopoDS::Face( theFace.Oriented( TopAbs_FORWARD )); + // check that face is not closed TopoDS_Vertex bidon; list eList; - getOrderedEdges( theFace, bidon, eList, myNbKeyPntInBoundary ); + getOrderedEdges( face, bidon, eList, myNbKeyPntInBoundary ); list::iterator elIt = eList.begin(); for ( ; elIt != eList.end() ; elIt++ ) - if ( BRep_Tool::IsClosed( *elIt , theFace )) + if ( BRep_Tool::IsClosed( *elIt , face )) return setErrorCode( ERR_LOADF_CLOSED_FACE ); Extrema_GenExtPS projector; - GeomAdaptor_Surface aSurface( BRep_Tool::Surface( theFace )); + GeomAdaptor_Surface aSurface( BRep_Tool::Surface( face )); if ( theProject || nbElems == 0 ) projector.Initialize( aSurface, 20,20, 1e-5,1e-5 ); @@ -567,7 +584,7 @@ bool SMESH_Pattern::Load (SMESH_Mesh* theMesh, TNodePointIDMap nodePointIDMap; if ( nbElems == 0 || (theProject && - theMesh->IsMainShape( theFace ) && + theMesh->IsMainShape( face ) && !isMeshBoundToShape( theMesh ))) { MESSAGE("Project the whole mesh"); @@ -579,8 +596,8 @@ bool SMESH_Pattern::Load (SMESH_Mesh* theMesh, SMDS_FaceIteratorPtr fIt = aMeshDS->facesIterator(); while ( fIt->more() ) { - myElemPointIDs.push_back( list< int >() ); - list< int >& elemPoints = myElemPointIDs.back(); + myElemPointIDs.push_back( TElemDef() ); + TElemDef& elemPoints = myElemPointIDs.back(); SMDS_ElemIteratorPtr nIt = fIt->next()->nodesIterator(); while ( nIt->more() ) { @@ -589,7 +606,7 @@ bool SMESH_Pattern::Load (SMESH_Mesh* theMesh, if ( nIdIt == nodePointIDMap.end() ) { elemPoints.push_back( iPoint ); - nodePointIDMap.insert( TNodePointIDMap::value_type( node, iPoint++ )); + nodePointIDMap.insert( make_pair( node, iPoint++ )); } else elemPoints.push_back( (*nIdIt).second ); @@ -608,11 +625,11 @@ bool SMESH_Pattern::Load (SMESH_Mesh* theMesh, 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(); @@ -648,7 +665,7 @@ bool SMESH_Pattern::Load (SMESH_Mesh* theMesh, for ( elIt = eList.begin(); elIt != eList.end(); elIt++ ) myShapeIDMap.Add( *elIt ); // the face - myShapeIDMap.Add( theFace ); + myShapeIDMap.Add( face ); myPoints.resize( nbNodes ); @@ -661,7 +678,7 @@ bool SMESH_Pattern::Load (SMESH_Mesh* theMesh, double f, l; Handle(Geom2d_Curve) C2d; if ( !theProject ) - C2d = BRep_Tool::CurveOnSurface( edge, theFace, f, l ); + C2d = BRep_Tool::CurveOnSurface( edge, face, f, l ); bool isForward = ( edge.Orientation() == TopAbs_FORWARD ); // the forward key-point @@ -674,7 +691,7 @@ bool SMESH_Pattern::Load (SMESH_Mesh* theMesh, myKeyPointIDs.push_back( iPoint ); SMDS_NodeIteratorPtr nIt = vSubMesh->GetNodes(); const SMDS_MeshNode* node = nIt->next(); - nodePointIDMap.insert( TNodePointIDMap::value_type( node, iPoint )); + nodePointIDMap.insert( make_pair( node, iPoint )); TPoint* keyPoint = &myPoints[ iPoint++ ]; vPoint.push_back( keyPoint ); @@ -714,7 +731,7 @@ 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 ( TNodePointIDMap::value_type( node, iPoint )); + nodePointIDMap.insert ( make_pair( node, iPoint )); if ( theProject ) p->myInitUV = project( node, projector ); @@ -738,7 +755,7 @@ bool SMESH_Pattern::Load (SMESH_Mesh* theMesh, myKeyPointIDs.push_back( iPoint ); SMDS_NodeIteratorPtr nIt = vSubMesh->GetNodes(); const SMDS_MeshNode* node = nIt->next(); - nodePointIDMap.insert( TNodePointIDMap::value_type( node, iPoint )); + nodePointIDMap.insert( make_pair( node, iPoint )); TPoint* keyPoint = &myPoints[ iPoint++ ]; vPoint2.push_back( keyPoint ); @@ -777,13 +794,13 @@ 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() ); - nodePointIDMap.insert( TNodePointIDMap::value_type( node, iPoint )); + nodePointIDMap.insert( make_pair( node, iPoint )); TPoint* p = &myPoints[ iPoint++ ]; fPoints.push_back( p ); if ( theProject ) @@ -799,8 +816,8 @@ bool SMESH_Pattern::Load (SMESH_Mesh* theMesh, SMDS_ElemIteratorPtr elemIt = fSubMesh->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() ]); } @@ -2234,7 +2251,6 @@ bool SMESH_Pattern::Apply (const TopoDS_Face& theFace, } int nbVertices = myShapeIDMap.Extent(); - //int nbSeamShapes = 0; // count twice seam edge and its vertices for ( elIt = eList.begin(); elIt != eList.end(); elIt++ ) myShapeIDMap.Add( *elIt ); @@ -2431,163 +2447,387 @@ bool SMESH_Pattern::Apply (const TopoDS_Face& theFace, return setErrorCode( ERR_OK ); } -// ========================================================= -// class calculating coordinates of 3D points by normalized -// parameters inside the block and vice versa -// ========================================================= +//======================================================================= +//function : Apply +//purpose : Compute nodes coordinates applying +// the loaded pattern to . The first key-point +// will be mapped into -th node +//======================================================================= -class TBlock: public math_FunctionSetWithDerivatives +bool SMESH_Pattern::Apply (const SMDS_MeshFace* theFace, + const int theNodeIndexOnKeyPoint1, + const bool theReverse) { - public: - enum TBlockShapeID { // ids of the block sub-shapes - ID_V000 = 1, ID_V100, ID_V010, ID_V110, ID_V001, ID_V101, ID_V011, ID_V111, - - ID_Ex00, ID_Ex10, ID_Ex01, ID_Ex11, - ID_E0y0, ID_E1y0, ID_E0y1, ID_E1y1, - ID_E00z, ID_E10z, ID_E01z, ID_E11z, - - ID_Fxy0, ID_Fxy1, ID_Fx0z, ID_Fx1z, ID_F0yz, ID_F1yz, - - ID_Shell - }; - static inline -bool IsVertexID( int theShapeID ) - { return ( theShapeID >= ID_V000 && theShapeID <= ID_V111 ); } - static inline bool IsEdgeID( int theShapeID ) - { return ( theShapeID >= ID_Ex00 && theShapeID <= ID_E11z ); } - static inline bool IsFaceID( int theShapeID ) - { return ( theShapeID >= ID_Fxy0 && theShapeID <= ID_F1yz ); } - - - TBlock (const TopoDS_Shell& theBlock): - myShell( theBlock ), myNbIterations(0), mySumDist(0.) {} - - bool LoadBlockShapes(const TopoDS_Vertex& theVertex000, - const TopoDS_Vertex& theVertex001, -// TopTools_IndexedMapOfShape& theShapeIDMap - TopTools_IndexedMapOfOrientedShape& theShapeIDMap ); - // add sub-shapes of theBlock to theShapeIDMap so that they get - // IDs acoording to enum TBlockShapeID - - static int GetShapeIDByParams ( const gp_XYZ& theParams ); - // define an id of the block sub-shape by point parameters - - bool VertexPoint( const int theVertexID, gp_XYZ& thePoint ) const { - if ( !IsVertexID( theVertexID )) return false; - thePoint = myPnt[ theVertexID - ID_V000 ]; return true; +// MESSAGE(" ::Apply(MeshFace) " ); + + 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 ); } - // return vertex coordinates - bool EdgePoint( const int theEdgeID, const gp_XYZ& theParams, gp_XYZ& thePoint ) const { - if ( !IsEdgeID( theEdgeID )) return false; - thePoint = myEdge[ theEdgeID - ID_Ex00 ].Point( theParams ); return true; + // 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 ); } - // return coordinates of a point on edge - bool FacePoint( const int theFaceID, const gp_XYZ& theParams, gp_XYZ& thePoint ) const { - if ( !IsFaceID ( theFaceID )) return false; - thePoint = myFace[ theFaceID - ID_Fxy0 ].Point( theParams ); return true; + // Define the nodes order + + 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 = static_cast( noIt->next() ); + nodes.push_back( node ); + if ( iSub++ == theNodeIndexOnKeyPoint1 ) + n = --nodes.end(); } - // return coordinates of a point on face - - bool ShellPoint( const gp_XYZ& theParams, gp_XYZ& thePoint ) const; - // return coordinates of a point in shell - - bool ComputeParameters (const gp_Pnt& thePoint, - gp_XYZ& theParams, - const int theShapeID = ID_Shell); - // compute point parameters in the block - - static void GetFaceEdgesIDs (const int faceID, vector< int >& edgeVec ); - // return edges IDs of a face in the order u0, u1, 0v, 1v - - static int GetCoordIndOnEdge (const int theEdgeID) - { return (theEdgeID < ID_E0y0) ? 1 : (theEdgeID < ID_E00z) ? 2 : 3; } - // return an index of a coordinate which varies along the edge - - static double* GetShapeCoef (const int theShapeID); - // for theShapeID( TBlockShapeID ), returns 3 coefficients used - // to compute an addition of an on-theShape point to coordinates - // of an in-shell point. If an in-shell point has parameters (Px,Py,Pz), - // then the addition of a point P is computed as P*kx*ky*kz and ki is - // defined by the returned coef like this: - // ki = (coef[i] == 0) ? 1 : (coef[i] < 0) ? 1 - Pi : Pi - - static bool IsForwardEdge (const TopoDS_Edge & theEdge, - //TopTools_IndexedMapOfShape& theShapeIDMap - TopTools_IndexedMapOfOrientedShape& theShapeIDMap) { - int v1ID = theShapeIDMap.FindIndex( TopExp::FirstVertex( theEdge ).Oriented( TopAbs_FORWARD )); - int v2ID = theShapeIDMap.FindIndex( TopExp::LastVertex( theEdge ).Oriented( TopAbs_FORWARD )); - return ( v1ID < v2ID ); + 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 ); } - // Return true if an in-block parameter increases along theEdge curve - - static void Swap(double& a, double& b) { double tmp = a; a = b; b = tmp; } - - // methods of math_FunctionSetWithDerivatives - Standard_Integer NbVariables() const; - Standard_Integer NbEquations() const; - Standard_Boolean Value(const math_Vector& X,math_Vector& F) ; - Standard_Boolean Derivatives(const math_Vector& X,math_Matrix& D) ; - Standard_Boolean Values(const math_Vector& X,math_Vector& F,math_Matrix& D) ; - Standard_Integer GetStateNumber (); - - static ostream& DumpShapeID (const int theBlockShapeID, ostream& stream); - // DEBUG: dump an id of a block sub-shape - - private: - - struct TEdge { - int myCoordInd; - double myFirst; - double myLast; - Handle(Geom_Curve) myC3d; - gp_Trsf myTrsf; - double GetU( const gp_XYZ& theParams ) const; - gp_XYZ Point( const gp_XYZ& theParams ) const; - }; - - struct TFace { - // 4 edges in the order u0, u1, 0v, 1v - int myCoordInd[ 4 ]; - double myFirst [ 4 ]; - double myLast [ 4 ]; - Handle(Geom2d_Curve) myC2d [ 4 ]; - // 4 corner points in the order 00, 10, 11, 01 - gp_XY myCorner [ 4 ]; - // surface - Handle(Geom_Surface) myS; - gp_Trsf myTrsf; - gp_XY GetUV( const gp_XYZ& theParams ) const; - gp_XYZ Point( const gp_XYZ& theParams ) const; - int GetUInd() const { return myCoordInd[ 0 ]; } - int GetVInd() const { return myCoordInd[ 2 ]; } - }; - - TopoDS_Shell myShell; - // geometry: - // 8 vertices - gp_XYZ myPnt[ 8 ]; - // 12 edges - TEdge myEdge[ 12 ]; - // 6 faces - TFace myFace[ 6 ]; - - // for param computation - - int myFaceIndex; - double myFaceParam; - int myNbIterations; - double mySumDist; - - gp_XYZ myPoint; // the given point - gp_XYZ myParam; // the best parameters guess - double myValues[ 4 ]; // values computed at myParam - - typedef pair TxyzPair; - TxyzPair my3x3x3GridNodes[ 27 ]; - bool myGridComputed; -}; + list< gp_XYZ > xyzList; + myOrderedNodes.resize( theFace->NbNodes() ); + for ( iSub = 0, n = nodes.begin(); n != nodes.end(); ++n ) { + xyzList.push_back( gp_XYZ( (*n)->X(), (*n)->Y(), (*n)->Z() )); + myOrderedNodes[ iSub++] = *n; + } + + // Define a face plane + + list< gp_XYZ >::iterator xyzIt = xyzList.begin(); + gp_Pnt P ( *xyzIt++ ); + gp_Vec Vx( P, *xyzIt++ ), N; + do { + N = Vx ^ gp_Vec( P, *xyzIt++ ); + } while ( N.SquareMagnitude() <= DBL_MIN && xyzIt != xyzList.end() ); + if ( N.SquareMagnitude() <= DBL_MIN ) + return setErrorCode( ERR_APPLF_BAD_FACE_GEOM ); + gp_Ax2 pos( P, N, Vx ); + + // Compute UV of key-points on a plane + for ( xyzIt = xyzList.begin(), iSub = 1; xyzIt != xyzList.end(); xyzIt++, iSub++ ) + { + gp_Vec vec ( pos.Location(), *xyzIt ); + TPoint* p = getShapePoints( iSub ).front(); + p->myUV.SetX( vec * pos.XDirection() ); + p->myUV.SetY( vec * pos.YDirection() ); + p->myXYZ = *xyzIt; + } + + // 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; + + // compute UV and XYZ of points on edges + + for ( xyzIt = xyzList.begin(); xyzIt != xyzList.end(); iSub++ ) + { + 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(); + while ( *pIt != ePoints.back() ) + { + TPoint* p = *pIt++; + p->myXYZ = xyz1 * ( 1 - p->myInitU ) + xyz2 * p->myInitU; + gp_Vec vec ( pos.Location(), p->myXYZ ); + p->myUV.SetX( vec * pos.XDirection() ); + p->myUV.SetY( vec * pos.YDirection() ); + } + // 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++ ) + { + (*pIt)->myXYZ = ElSLib::PlaneValue( (*pIt)->myUV.X(), (*pIt)->myUV.Y(), pos ); + } + + 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 (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 + + // apply to each face in theFaces set + set::iterator face = theFaces.begin(); + for ( ; face != theFaces.end(); ++face ) + { + if ( !Apply( *face, theNodeIndexOnKeyPoint1, theReverse )) { + 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 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 + 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 ); + } + ind1 += myPoints.size(); + } + + return !myElemXYZIDs.empty(); +} + +//======================================================================= +//function : Apply +//purpose : Compute nodes coordinates applying +// the loaded pattern to . The (0,0,0) key-point +// will be mapped into -th node. The +// (0,0,1) key-point will be mapped into -th +// node. +//======================================================================= + +bool SMESH_Pattern::Apply (std::set & theVolumes, + const int theNode000Index, + const int theNode001Index) +{ + MESSAGE(" ::Apply(set) " ); + + if ( !IsLoaded() ) { + MESSAGE( "Pattern not loaded" ); + return setErrorCode( ERR_APPL_NOT_LOADED ); + } + + // bind ID to points + 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() * theVolumes.size(), undefinedXYZ() ); + myElements.reserve( theVolumes.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 an element + + // apply to each element in theVolumes set + set::iterator vol = theVolumes.begin(); + for ( ; vol != theVolumes.end(); ++vol ) + { + if ( !Apply( *vol, theNode000Index, theNode001Index )) { + MESSAGE( "Failed on " << *vol ); + continue; + } + myElements.push_back( *vol ); + + // 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 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 )) { + subNodes.insert( myOrderedNodes[ Id - 1 ]); + } + else if ( SMESH_Block::IsEdgeID( Id )) { + SMESH_Block::GetEdgeVertexIDs( Id, subIDs ); + subNodes.insert( myOrderedNodes[ subIDs.front() - 1 ]); + subNodes.insert( myOrderedNodes[ subIDs.back() - 1 ]); + } + else { + SMESH_Block::GetFaceEdgesIDs( Id, subIDs ); + int e1 = subIDs[ 0 ], e2 = subIDs[ 1 ]; + SMESH_Block::GetEdgeVertexIDs( e1, subIDs ); + subNodes.insert( myOrderedNodes[ subIDs.front() - 1 ]); + subNodes.insert( myOrderedNodes[ subIDs.back() - 1 ]); + SMESH_Block::GetEdgeVertexIDs( e2, subIDs ); + subNodes.insert( myOrderedNodes[ subIDs.front() - 1 ]); + subNodes.insert( myOrderedNodes[ subIDs.back() - 1 ]); + } + // add points + list< TPoint* > & points = getShapePoints( Id ); + list< TPoint* >::iterator p = points.begin(); + 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(); + } + + return !myElemXYZIDs.empty(); +} //======================================================================= //function : Load @@ -2600,12 +2840,12 @@ 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 - TBlock block( theBlock ); + SMESH_Block block; TopoDS_Vertex v1, v2; - if ( !block.LoadBlockShapes( v1, v2, myShapeIDMap )) + if ( !block.LoadBlockShapes( theBlock, v1, v2, myShapeIDMap )) return setErrorCode( ERR_LOADV_BAD_SHAPE ); // count nodes @@ -2613,7 +2853,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(); } @@ -2626,7 +2866,7 @@ 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; @@ -2634,7 +2874,7 @@ bool SMESH_Pattern::Load (SMESH_Mesh* theMesh, // store a node and a point while ( nIt->more() ) { const SMDS_MeshNode* node = static_cast( nIt->next() ); - nodePointIDMap.insert( TNodePointIDMap::value_type( node, iPoint )); + nodePointIDMap.insert( make_pair( node, iPoint )); if ( block.IsVertexID( shapeID )) myKeyPointIDs.push_back( iPoint ); TPoint* p = & myPoints[ iPoint++ ]; @@ -2662,8 +2902,8 @@ bool SMESH_Pattern::Load (SMESH_Mesh* theMesh, const TopoDS_Edge& edge = TopoDS::Edge( S ); double f,l; BRep_Tool::Range( edge, f, l ); - int iCoord = TBlock::GetCoordIndOnEdge( shapeID ); - bool isForward = TBlock::IsForwardEdge( edge, myShapeIDMap ); + int iCoord = SMESH_Block::GetCoordIndOnEdge( shapeID ); + bool isForward = SMESH_Block::IsForwardEdge( edge, myShapeIDMap ); pIt = shapePoints.begin(); nIt = aSubMesh->GetNodes(); for ( ; nIt->more(); pIt++ ) @@ -2690,14 +2930,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() ]); } @@ -2708,6 +2948,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 @@ -2722,13 +2988,12 @@ bool SMESH_Pattern::Apply (const TopoDS_Shell& theBlock, { MESSAGE(" ::Apply(volume) " ); - TBlock block( theBlock ); - if (!findBoundaryPoints() || // bind ID to points !setShapeToMesh( theBlock )) // check theBlock is a suitable shape return false; - if (!block.LoadBlockShapes( theVertex000, theVertex001, myShapeIDMap )) // bind ID to shape + SMESH_Block block; // bind ID to shape + if (!block.LoadBlockShapes( theBlock, theVertex000, theVertex001, myShapeIDMap )) return setErrorCode( ERR_APPLV_BAD_SHAPE ); // compute XYZ of points on shapes @@ -2770,1541 +3035,1352 @@ bool SMESH_Pattern::Apply (const TopoDS_Shell& theBlock, } //======================================================================= -//function : MakeMesh -//purpose : Create nodes and elements in using nodes -// coordinates computed by either of Apply...() methods +//function : Apply +//purpose : Compute nodes coordinates applying +// the loaded pattern to . The (0,0,0) key-point +// will be mapped into -th node. The +// (0,0,1) key-point will be mapped into -th +// node. //======================================================================= -bool SMESH_Pattern::MakeMesh(SMESH_Mesh* theMesh) +bool SMESH_Pattern::Apply (const SMDS_MeshVolume* theVolume, + const int theNode000Index, + const int theNode001Index) { - MESSAGE(" ::MakeMesh() " ); - if ( !myIsComputed ) - return setErrorCode( ERR_MAKEM_NOT_COMPUTED ); + //MESSAGE(" ::Apply(MeshVolume) " ); - SMESHDS_Mesh* aMeshDS = theMesh->GetMeshDS(); + if (!findBoundaryPoints()) // bind ID to points + return false; - // clear elements and nodes existing on myShape - SMESH_subMesh * aSubMesh = theMesh->GetSubMeshContaining( myShape ); - SMESHDS_SubMesh * aSubMeshDS = aMeshDS->MeshElements( myShape ); - if ( aSubMesh ) - aSubMesh->ComputeStateEngine( SMESH_subMesh::CLEAN ); - else if ( aSubMeshDS ) - { - SMDS_ElemIteratorPtr eIt = aSubMeshDS->GetElements(); - while ( eIt->more() ) - aMeshDS->RemoveElement( eIt->next() ); - SMDS_NodeIteratorPtr nIt = aSubMeshDS->GetNodes(); - while ( nIt->more() ) - aMeshDS->RemoveNode( static_cast( nIt->next() )); - } + SMESH_Block block; // bind ID to shape + if (!block.LoadMeshBlock( theVolume, theNode000Index, theNode001Index, myOrderedNodes )) + return setErrorCode( ERR_APPLV_BAD_SHAPE ); + // compute XYZ of points on shapes - // loop on sub-shapes of myShape: create nodes and build point-node map - typedef map< TPoint*, const SMDS_MeshNode* > TPointNodeMap; - TPointNodeMap pointNodeMap; - map< int, list< TPoint* > >::iterator idPointIt = myShapeIDToPointsMap.begin(); - for ( ; idPointIt != myShapeIDToPointsMap.end(); idPointIt++ ) + for ( int ID = SMESH_Block::ID_V000; ID <= SMESH_Block::ID_Shell; ID++ ) { - const TopoDS_Shape & S = myShapeIDMap( (*idPointIt).first ); - list< TPoint* > & points = (*idPointIt).second; - SMESHDS_SubMesh * subMeshDS = aMeshDS->MeshElements( S ); - SMESH_subMesh * subMesh = theMesh->GetSubMeshContaining( myShape ); - list< TPoint* >::iterator pIt = points.begin(); - for ( ; pIt != points.end(); pIt++ ) - { - TPoint* point = *pIt; - if ( pointNodeMap.find( point ) != pointNodeMap.end() ) - continue; - SMDS_MeshNode* node = aMeshDS->AddNode (point->myXYZ.X(), - point->myXYZ.Y(), - point->myXYZ.Z()); - pointNodeMap.insert( TPointNodeMap::value_type( point, node )); - if ( subMeshDS ) { - switch ( S.ShapeType() ) { - case TopAbs_VERTEX: { - 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; - } - 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; - } - default: - aMeshDS->SetNodeInVolume( node, TopoDS::Shell( S )); - } - } - } - // make that SMESH_subMesh::_computeState = COMPUTE_OK - // so that operations with hypotheses will erase the mesh - // being built - if ( subMesh ) - subMesh->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE ); - } + list< TPoint* > & shapePoints = getShapePoints( ID ); + list< TPoint* >::iterator pIt = shapePoints.begin(); - // create elements - list >::iterator epIt = myElemPointIDs.begin(); - for ( ; epIt != myElemPointIDs.end(); epIt++ ) - { - list< int > & elemPoints = *epIt; - // retrieve nodes - const SMDS_MeshNode* nodes[ 8 ]; - list< int >::iterator iIt = elemPoints.begin(); - int nbNodes; - for ( nbNodes = 0; iIt != elemPoints.end(); iIt++ ) { - nodes[ nbNodes++ ] = pointNodeMap[ & myPoints[ *iIt ]]; - } - // add an element - const SMDS_MeshElement* elem = 0; - if ( myIs2D ) { - 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; - default:; + if ( block.IsVertexID( ID )) + for ( ; pIt != shapePoints.end(); pIt++ ) { + block.VertexPoint( ID, (*pIt)->myXYZ.ChangeCoord() ); } - } - else { - switch ( nbNodes ) { - case 4: - elem = aMeshDS->AddVolume (nodes[0], nodes[1], nodes[2], nodes[3] ); break; - case 5: - elem = aMeshDS->AddVolume (nodes[0], nodes[1], nodes[2], nodes[3], - nodes[4] ); break; - case 6: - elem = aMeshDS->AddVolume (nodes[0], nodes[1], nodes[2], nodes[3], - nodes[4], nodes[5] ); break; - case 8: - elem = aMeshDS->AddVolume (nodes[0], nodes[1], nodes[2], nodes[3], - nodes[4], nodes[5], nodes[6], nodes[7] ); break; - default:; + else if ( block.IsEdgeID( ID )) + for ( ; pIt != shapePoints.end(); pIt++ ) { + block.EdgePoint( ID, (*pIt)->myInitXYZ, (*pIt)->myXYZ.ChangeCoord() ); } - } - if ( elem ) - aMeshDS->SetMeshElementOnShape( elem, myShape ); - } + else if ( block.IsFaceID( ID )) + for ( ; pIt != shapePoints.end(); pIt++ ) { + block.FacePoint( ID, (*pIt)->myInitXYZ, (*pIt)->myXYZ.ChangeCoord() ); + } + else + for ( ; pIt != shapePoints.end(); pIt++ ) + block.ShellPoint( (*pIt)->myInitXYZ, (*pIt)->myXYZ.ChangeCoord() ); + } // loop on block sub-shapes + + myIsComputed = true; return setErrorCode( ERR_OK ); } - //======================================================================= -//function : arrangeBoundaries -//purpose : if there are several wires, arrange boundaryPoints so that -// the outer wire goes first and fix inner wires orientation -// update myKeyPointIDs to correspond to the order of key-points -// in boundaries; sort internal boundaries by the nb of key-points +//function : mergePoints +//purpose : Merge XYZ on edges and/or faces. //======================================================================= -void SMESH_Pattern::arrangeBoundaries (list< list< TPoint* > >& boundaryList) +void SMESH_Pattern::mergePoints (const bool uniteGroups) { - typedef list< list< TPoint* > >::iterator TListOfListIt; - TListOfListIt bndIt; - list< TPoint* >::iterator pIt; - - int nbBoundaries = boundaryList.size(); - if ( nbBoundaries > 1 ) + map< TNodeSet, list< list< int > > >::iterator idListIt = myIdsOnBoundary.begin(); + for ( ; idListIt != myIdsOnBoundary.end(); idListIt++ ) { - // sort boundaries by nb of key-points - if ( nbBoundaries > 2 ) - { - // move boundaries in tmp list - 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 - typedef map< int, TListOfListIt > TNbKpBndPosMap; - TNbKpBndPosMap nbKpBndPosMap; - bndIt = tmpList.begin(); - list< int >::iterator nbKpIt = myNbKeyPntInBoundary.begin(); - for ( ; nbKpIt != myNbKeyPntInBoundary.end(); nbKpIt++, bndIt++ ) { - int nb = *nbKpIt * nbBoundaries; - while ( nbKpBndPosMap.find ( nb ) != nbKpBndPosMap.end() ) - nb++; - nbKpBndPosMap.insert( TNbKpBndPosMap::value_type( nb, bndIt )); - } - // move boundaries back to boundaryList - TNbKpBndPosMap::iterator nbKpBndPosIt = nbKpBndPosMap.begin(); - for ( ; nbKpBndPosIt != nbKpBndPosMap.end(); nbKpBndPosIt++ ) { - TListOfListIt & bndPos2 = (*nbKpBndPosIt).second; - TListOfListIt bndPos1 = bndPos2++; - boundaryList.splice( boundaryList.end(), tmpList, bndPos1, bndPos2 ); - } + list >& groups = idListIt->second; + if ( groups.size() < 2 ) + continue; + + // 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 ); } - // Look for the outer boundary: the one with the point with the least X - double leastX = DBL_MAX; - TListOfListIt outerBndPos; - for ( bndIt = boundaryList.begin(); bndIt != boundaryList.end(); bndIt++ ) + // 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() ); + + // compare points, replace indices + + list< int >::iterator ind1, ind2; + list< list< int > >::iterator grpIt1, grpIt2; + for ( grpIt1 = groups.begin(); grpIt1 != groups.end(); grpIt1++ ) { - list< TPoint* >& boundary = (*bndIt); - for ( pIt = boundary.begin(); pIt != boundary.end(); pIt++) + list< int >& indices1 = *grpIt1; + grpIt2 = grpIt1; + for ( grpIt2++; grpIt2 != groups.end(); grpIt2++ ) { - TPoint* point = *pIt; - if ( point->myInitXYZ.X() < leastX ) { - leastX = point->myInitXYZ.X(); - outerBndPos = bndIt; + 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 +} - if ( outerBndPos != boundaryList.begin() ) - boundaryList.splice( boundaryList.begin(), boundaryList, outerBndPos, ++outerBndPos ); +//======================================================================= +//function : makePolyElements +//purpose : prepare intermediate data to create Polygons and Polyhedrons +//======================================================================= - } // if nbBoundaries > 1 - - // Check boundaries orientation and re-fill myKeyPointIDs +void SMESH_Pattern:: + makePolyElements(const vector< const SMDS_MeshNode* >& theNodes, + const bool toCreatePolygons, + const bool toCreatePolyedrs) +{ + myPolyElemXYZIDs.clear(); + myPolyElems.clear(); + myPolyElems.reserve( myIdsOnBoundary.size() ); - set< TPoint* > keyPointSet; - list< int >::iterator kpIt = myKeyPointIDs.begin(); - for ( ; kpIt != myKeyPointIDs.end(); kpIt++ ) - keyPointSet.insert( & myPoints[ *kpIt ]); - myKeyPointIDs.clear(); + // make a set of refined elements + set< const SMDS_MeshElement* > avoidSet, elemSet; + avoidSet.insert( myElements.begin(), myElements.end() ); - // update myNbKeyPntInBoundary also - list< int >::iterator nbKpIt = myNbKeyPntInBoundary.begin(); + map< TNodeSet, list< list< int > > >::iterator indListIt, nn_IdList; - for ( bndIt = boundaryList.begin(); bndIt != boundaryList.end(); bndIt++, nbKpIt++ ) + if ( toCreatePolygons ) { - // find the point with the least X - double leastX = DBL_MAX; - list< TPoint* >::iterator xpIt; - list< TPoint* >& boundary = (*bndIt); - for ( pIt = boundary.begin(); pIt != boundary.end(); pIt++) - { - TPoint* point = *pIt; - if ( point->myInitXYZ.X() < leastX ) { - leastX = point->myInitXYZ.X(); - xpIt = pIt; - } - } - // find points next to the point with the least X - TPoint* p = *xpIt, *pPrev, *pNext; - if ( p == boundary.front() ) - pPrev = *(++boundary.rbegin()); - else { - xpIt--; - pPrev = *xpIt; - xpIt++; - } - if ( p == boundary.back() ) - pNext = *(++boundary.begin()); - else { - xpIt++; - pNext = *xpIt; - } - // vectors of boundary direction near

- gp_Vec2d v1( pPrev->myInitUV, p->myInitUV ), v2( p->myInitUV, pNext->myInitUV ); - // rotate vectors counterclockwise - v1.SetCoord( -v1.Y(), v1.X() ); - v2.SetCoord( -v2.Y(), v2.X() ); - // in-face direction - gp_Vec2d dirInFace = v1 + v2; - // for the outer boundary dirInFace should go to the right - bool reverse; - if ( bndIt == boundaryList.begin() ) // outer boundary - reverse = dirInFace.X() < 0; - else - reverse = dirInFace.X() > 0; - if ( reverse ) - boundary.reverse(); + int lastFreeId = myXYZ.size(); - // Put key-point IDs of a well-oriented boundary in myKeyPointIDs - (*nbKpIt) = 0; // count nb of key-points again - pIt = boundary.begin(); - for ( ; pIt != boundary.end(); pIt++) + // loop on links of refined elements + indListIt = myIdsOnBoundary.begin(); + for ( ; indListIt != myIdsOnBoundary.end(); indListIt++ ) { - TPoint* point = *pIt; - if ( keyPointSet.find( point ) == keyPointSet.end() ) + 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 an index of a keypoint - int index = 0; - vector< TPoint >::const_iterator pVecIt = myPoints.begin(); - for ( ; pVecIt != myPoints.end(); pVecIt++, index++ ) - if ( &(*pVecIt) == point ) + + // 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++ ] = static_cast( 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; - myKeyPointIDs.push_back( index ); - (*nbKpIt)++; - } - myKeyPointIDs.pop_back(); // remove the first key-point from the back - (*nbKpIt)--; + } // while (true) - } // loop on a list of boundaries + if ( myIs2D && idGroups.size() > 1 ) { - ASSERT( myKeyPointIDs.size() == keyPointSet.size() ); + // 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 = static_cast( 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 : findBoundaryPoints -//purpose : if loaded from file, find points to map on edges and faces and -// compute their parameters +//function : getFacesDefinition +//purpose : return faces definition for a volume face defined by theBndNodes //======================================================================= -bool SMESH_Pattern::findBoundaryPoints() +bool SMESH_Pattern:: + getFacesDefinition(const SMDS_MeshNode** theBndNodes, + const int theNbBndNodes, + const vector< const SMDS_MeshNode* >& theNodes, + list< int >& theFaceDefs, + vector& theQuantity) { - if ( myIsBoundaryPointsFound ) return true; - - MESSAGE(" findBoundaryPoints() "); + 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(); - if ( myIs2D ) + // 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 ) { - set< TPoint* > pointsInElems; + // 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 free links of elements: - // put links of all elements in a set and remove links encountered twice + // find faces definition of new volumes - typedef pair< TPoint*, TPoint*> TLink; - set< TLink > linkSet; - list >::iterator epIt = myElemPointIDs.begin(); - for ( ; epIt != myElemPointIDs.end(); epIt++ ) + 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 ) { - list< int > & elemPoints = *epIt; - list< int >::iterator pIt = elemPoints.begin(); - int prevP = elemPoints.back(); - for ( ; pIt != elemPoints.end(); pIt++ ) { - TPoint* p1 = & myPoints[ prevP ]; - TPoint* p2 = & myPoints[ *pIt ]; - TLink link(( p1 < p2 ? p1 : p2 ), ( p1 < p2 ? p2 : p1 )); - ASSERT( link.first != link.second ); - pair::iterator,bool> itUniq = linkSet.insert( link ); - if ( !itUniq.second ) - linkSet.erase( itUniq.first ); - prevP = *pIt; - - pointsInElems.insert( p1 ); + // 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; + } + } } } - // Now linkSet contains only free links, - // find the points order that they have in boundaries + } + if ( !defsAdded ) { + theQuantity.push_back( faceDef.size() ); + theFaceDefs.splice( theFaceDefs.end(), faceDef, faceDef.begin(), faceDef.end() ); + } - // 1. make a map of key-points - set< TPoint* > keyPointSet; - list< int >::iterator kpIt = myKeyPointIDs.begin(); - for ( ; kpIt != myKeyPointIDs.end(); kpIt++ ) - keyPointSet.insert( & myPoints[ *kpIt ]); + return makePoly; +} - // 2. chain up boundary points - list< list< TPoint* > > boundaryList; - boundaryList.push_back( list< TPoint* >() ); - list< TPoint* > * boundary = & boundaryList.back(); +//======================================================================= +//function : clearMesh +//purpose : clear mesh elements existing on myShape in theMesh +//======================================================================= - TPoint *point1, *point2, *keypoint1; - kpIt = myKeyPointIDs.begin(); - point1 = keypoint1 = & myPoints[ *kpIt++ ]; - // loop on free links: look for the next point - int iKeyPoint = 0; - set< TLink >::iterator lIt = linkSet.begin(); - while ( lIt != linkSet.end() ) - { - if ( (*lIt).first == point1 ) - point2 = (*lIt).second; - else if ( (*lIt).second == point1 ) - point2 = (*lIt).first; - else { - lIt++; - continue; - } - linkSet.erase( lIt ); - lIt = linkSet.begin(); +void SMESH_Pattern::clearMesh(SMESH_Mesh* theMesh) const +{ - if ( keyPointSet.find( point2 ) == keyPointSet.end() ) // not a key-point + if ( !myShape.IsNull() ) + { + if ( SMESH_subMesh * aSubMesh = theMesh->GetSubMesh/*Containing*/( myShape )) + { + aSubMesh->ComputeStateEngine( SMESH_subMesh::CLEANDEP ); + } + else { + SMESHDS_Mesh* aMeshDS = theMesh->GetMeshDS(); + if ( SMESHDS_SubMesh* aSubMeshDS = aMeshDS->MeshElements( myShape )) { - boundary->push_back( point2 ); + SMDS_ElemIteratorPtr eIt = aSubMeshDS->GetElements(); + while ( eIt->more() ) + aMeshDS->RemoveElement( eIt->next() ); + SMDS_NodeIteratorPtr nIt = aSubMeshDS->GetNodes(); + while ( nIt->more() ) + aMeshDS->RemoveNode( static_cast( nIt->next() )); } - else // a key-point found - { - keyPointSet.erase( point2 ); // keyPointSet contains not found key-points only - iKeyPoint++; - if ( point2 != keypoint1 ) // its not the boundary end - { - boundary->push_back( point2 ); - } - else // the boundary end reached - { - boundary->push_front( keypoint1 ); - boundary->push_back( keypoint1 ); - myNbKeyPntInBoundary.push_back( iKeyPoint ); - if ( keyPointSet.empty() ) - break; // all boundaries containing key-points are found + } + } +} - // prepare to search for the next boundary - boundaryList.push_back( list< TPoint* >() ); - boundary = & boundaryList.back(); - point2 = keypoint1 = (*keyPointSet.begin()); - } - } - point1 = point2; - } // loop on the free links set +//======================================================================= +//function : MakeMesh +//purpose : Create nodes and elements in using nodes +// coordinates computed by either of Apply...() methods +//======================================================================= - if ( boundary->empty() ) { - MESSAGE(" a separate key-point"); - return setErrorCode( ERR_READ_BAD_KEY_POINT ); - } +bool SMESH_Pattern::MakeMesh(SMESH_Mesh* theMesh, + const bool toCreatePolygons, + const bool toCreatePolyedrs) +{ + MESSAGE(" ::MakeMesh() " ); + if ( !myIsComputed ) + return setErrorCode( ERR_MAKEM_NOT_COMPUTED ); - // if there are several wires, arrange boundaryPoints so that - // the outer wire goes first and fix inner wires orientation; - // sort myKeyPointIDs to correspond to the order of key-points - // in boundaries - arrangeBoundaries( boundaryList ); + mergePoints( toCreatePolygons ); - // Find correspondence shape ID - points, - // compute points parameter on edge + SMESHDS_Mesh* aMeshDS = theMesh->GetMeshDS(); - keyPointSet.clear(); - for ( kpIt = myKeyPointIDs.begin(); kpIt != myKeyPointIDs.end(); kpIt++ ) - keyPointSet.insert( & myPoints[ *kpIt ]); + // clear elements and nodes existing on myShape + clearMesh(theMesh); - set< TPoint* > edgePointSet; // to find in-face points - int vertexID = 1; // the first index in TopTools_IndexedMapOfShape - int edgeID = myKeyPointIDs.size() + 1; + bool onMeshElements = ( !myElements.empty() ); - list< list< TPoint* > >::iterator bndIt = boundaryList.begin(); - for ( ; bndIt != boundaryList.end(); bndIt++ ) + // Create missing nodes + + vector< const SMDS_MeshNode* > nodesVector; // i-th point/xyz -> node + if ( onMeshElements ) + { + 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 ) { + if ( !nodesVector[ i ] && isDefined( myXYZ[ i ] ) ) + nodesVector[ i ] = aMeshDS->AddNode (myXYZ[ i ].X(), + myXYZ[ i ].Y(), + myXYZ[ i ].Z()); + } + } + 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++ ) { - boundary = & (*bndIt); - double edgeLength = 0; - list< TPoint* >::iterator pIt = boundary->begin(); - getShapePoints( edgeID ).push_back( *pIt ); - for ( pIt++; pIt != boundary->end(); pIt++) + TopoDS_Shape S; + SMESHDS_SubMesh * subMeshDS = 0; + if ( !myShapeIDMap.IsEmpty() ) { + S = myShapeIDMap( idPointIt->first ); + subMeshDS = aMeshDS->MeshElements( S ); + } + list< TPoint* > & points = idPointIt->second; + list< TPoint* >::iterator pIt = points.begin(); + for ( ; pIt != points.end(); pIt++ ) { - list< TPoint* > & edgePoints = getShapePoints( edgeID ); - TPoint* prevP = edgePoints.empty() ? 0 : edgePoints.back(); TPoint* point = *pIt; - edgePointSet.insert( point ); - if ( keyPointSet.find( point ) == keyPointSet.end() ) // inside-edge point - { - edgePoints.push_back( point ); - edgeLength += ( point->myInitUV - prevP->myInitUV ).Modulus(); - point->myInitU = edgeLength; - } - else // a key-point - { - // treat points on the edge which ends up: compute U [0,1] - edgePoints.push_back( point ); - if ( edgePoints.size() > 2 ) { - edgeLength += ( point->myInitUV - prevP->myInitUV ).Modulus(); - list< TPoint* >::iterator epIt = edgePoints.begin(); - for ( ; epIt != edgePoints.end(); epIt++ ) - (*epIt)->myInitU /= edgeLength; + int pIndex = pointIndex[ point ]; + if ( nodesVector [ pIndex ] ) + continue; + SMDS_MeshNode* node = aMeshDS->AddNode (point->myXYZ.X(), + point->myXYZ.Y(), + point->myXYZ.Z()); + nodesVector [ pIndex ] = node; + + if ( subMeshDS ) { + switch ( S.ShapeType() ) { + case TopAbs_VERTEX: { + 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; + } + 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; + } + default: + aMeshDS->SetNodeInVolume( node, TopoDS::Shell( S )); } - // begin the next edge treatment - edgeLength = 0; - getShapePoints( vertexID++ ).push_back( point ); - edgeID++; - if ( point != boundary->front() ) - getShapePoints( edgeID ).push_back( point ); } } } + } - // find in-face points - list< TPoint* > & facePoints = getShapePoints( edgeID ); - vector< TPoint >::iterator pVecIt = myPoints.begin(); - for ( ; pVecIt != myPoints.end(); pVecIt++ ) { - TPoint* point = &(*pVecIt); - if ( edgePointSet.find( point ) == edgePointSet.end() && - pointsInElems.find( point ) != pointsInElems.end()) - facePoints.push_back( point ); - } + // create elements - } // 2D case + if ( onMeshElements ) + { + // prepare data to create poly elements + makePolyElements( nodesVector, toCreatePolygons, toCreatePolyedrs ); - else // 3D case + // refine elements + createElements( theMesh, nodesVector, myElemXYZIDs, myElements ); + // sew old and new elements + createElements( theMesh, nodesVector, myPolyElemXYZIDs, myPolyElems ); + } + else { - // bind points to shapes according to point parameters - vector< TPoint >::iterator pVecIt = myPoints.begin(); - for ( int i = 0; pVecIt != myPoints.end(); pVecIt++, i++ ) { - TPoint* point = &(*pVecIt); - int shapeID = TBlock::GetShapeIDByParams( point->myInitXYZ ); - getShapePoints( shapeID ).push_back( point ); - // detect key-points - if ( TBlock::IsVertexID( shapeID )) - myKeyPointIDs.push_back( i ); - } + createElements( theMesh, nodesVector, myElemPointIDs, myElements ); } - myIsBoundaryPointsFound = true; - return myIsBoundaryPointsFound; +// 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 : Clear -//purpose : clear fields +//function : createElements +//purpose : add elements to the mesh //======================================================================= -void SMESH_Pattern::Clear() +void SMESH_Pattern::createElements(SMESH_Mesh* theMesh, + const vector& theNodesVector, + const list< TElemDef > & theElemNodeIDs, + const vector& theElements) { - myIsComputed = myIsBoundaryPointsFound = false; - - myPoints.clear(); - myKeyPointIDs.clear(); - myElemPointIDs.clear(); - myShapeIDToPointsMap.clear(); - myShapeIDMap.Clear(); - myShape.Nullify(); - myNbKeyPntInBoundary.clear(); -} + SMESHDS_Mesh* aMeshDS = theMesh->GetMeshDS(); + SMESH_MeshEditor editor( theMesh ); -//======================================================================= -//function : setShapeToMesh -//purpose : set a shape to be meshed. Return True if meshing is possible -//======================================================================= + bool onMeshElements = !theElements.empty(); -bool SMESH_Pattern::setShapeToMesh(const TopoDS_Shape& theShape) -{ - if ( !IsLoaded() ) { - MESSAGE( "Pattern not loaded" ); - return setErrorCode( ERR_APPL_NOT_LOADED ); + // shapes and groups theElements are on + vector< int > shapeIDs; + vector< list< SMESHDS_Group* > > groups; + set< const SMDS_MeshNode* > shellNodes; + if ( onMeshElements ) + { + shapeIDs.resize( theElements.size() ); + groups.resize( theElements.size() ); + const set& allGroups = aMeshDS->GetGroups(); + set::const_iterator grIt; + for ( int i = 0; i < theElements.size(); 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( 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() ); + } + } + } } + // nb new elements per a refined element + int nbNewElemsPerOld = 1; + if ( onMeshElements ) + nbNewElemsPerOld = theElemNodeIDs.size() / theElements.size(); - TopAbs_ShapeEnum aType = theShape.ShapeType(); - bool dimOk = ( myIs2D ? aType == TopAbs_FACE : aType == TopAbs_SHELL ); - if ( !dimOk ) { - MESSAGE( "Pattern dimention mismatch" ); - return setErrorCode( ERR_APPL_BAD_DIMENTION ); - } + bool is2d = myIs2D; - // check if a face is closed - int nbNodeOnSeamEdge = 0; - if ( myIs2D ) { - 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; - } - - // check nb of vertices - TopTools_IndexedMapOfShape vMap; - TopExp::MapShapes( theShape, TopAbs_VERTEX, vMap ); - if ( vMap.Extent() + nbNodeOnSeamEdge != myKeyPointIDs.size() ) { - MESSAGE( myKeyPointIDs.size() << " != " << vMap.Extent() ); - return setErrorCode( ERR_APPL_BAD_NB_VERTICES ); + list< TElemDef >::const_iterator enIt = theElemNodeIDs.begin(); + list< vector >::iterator quantity = myPolyhedronQuantities.begin(); + for ( int iElem = 0; enIt != theElemNodeIDs.end(); enIt++, iElem++ ) + { + const TElemDef & elemNodeInd = *enIt; + // retrieve nodes + vector< const SMDS_MeshNode* > nodes( elemNodeInd.size() ); + TElemDef::const_iterator id = elemNodeInd.begin(); + int nbNodes; + for ( nbNodes = 0; id != elemNodeInd.end(); id++ ) { + if ( *id < theNodesVector.size() ) + nodes[ nbNodes++ ] = theNodesVector[ *id ]; + else + 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 ( 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; + default: + elem = aMeshDS->AddPolygonalFace( nodes ); + } + } + else { + switch ( nbNodes ) { + case 4: + elem = aMeshDS->AddVolume (nodes[0], nodes[1], nodes[2], nodes[3] ); break; + case 5: + elem = aMeshDS->AddVolume (nodes[0], nodes[1], nodes[2], nodes[3], + nodes[4] ); break; + case 6: + elem = aMeshDS->AddVolume (nodes[0], nodes[1], nodes[2], nodes[3], + nodes[4], nodes[5] ); break; + case 8: + elem = aMeshDS->AddVolume (nodes[0], nodes[1], nodes[2], nodes[3], + nodes[4], nodes[5], nodes[6], nodes[7] ); break; + default: + elem = aMeshDS->AddPolyhedralVolume( nodes, *quantity++ ); + } + } + // set element on a shape + if ( elem && onMeshElements ) // applied to mesh elements + { + 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 + ( static_cast( 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() ) // applied to shape + aMeshDS->SetMeshElementOnShape( elem, myShape ); } - myShapeIDMap.Clear(); - myShape = theShape; - return true; -} - -//======================================================================= -//function : GetMappedPoints -//purpose : Return nodes coordinates computed by Apply() method -//======================================================================= - -bool SMESH_Pattern::GetMappedPoints ( list< const gp_XYZ * > & thePoints ) -{ - thePoints.clear(); - if ( !myIsComputed ) - return false; + // make that SMESH_subMesh::_computeState == COMPUTE_OK + // so that operations with hypotheses will erase the mesh being built - vector< TPoint >::iterator pVecIt = myPoints.begin(); - for ( ; pVecIt != myPoints.end(); pVecIt++ ) - thePoints.push_back( & (*pVecIt).myXYZ.XYZ() ); + SMESH_subMesh * subMesh; + if ( !myShape.IsNull() ) { + subMesh = theMesh->GetSubMeshContaining( myShape ); + if ( subMesh ) + subMesh->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE ); + } + if ( onMeshElements ) { + list< int > elemIDs; + for ( int i = 0; i < theElements.size(); i++ ) + { + subMesh = theMesh->GetSubMeshContaining( shapeIDs[ i ] ); + if ( subMesh ) + subMesh->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE ); - return ( thePoints.size() > 0 ); + elemIDs.push_back( theElements[ i ]->GetID() ); + } + // remove refined elements + editor.Remove( elemIDs, false ); + } } - //======================================================================= -//function : GetPoints -//purpose : Return nodes coordinates of the pattern +//function : isReversed +//purpose : check xyz ids order in theIdsList taking into account +// theFirstNode on a link //======================================================================= -bool SMESH_Pattern::GetPoints ( list< const gp_XYZ * > & thePoints ) const +bool SMESH_Pattern::isReversed(const SMDS_MeshNode* theFirstNode, + const list< int >& theIdsList) const { - thePoints.clear(); - - if ( !IsLoaded() ) + if ( theIdsList.size() < 2 ) return false; - vector< TPoint >::const_iterator pVecIt = myPoints.begin(); - for ( ; pVecIt != myPoints.end(); pVecIt++ ) - thePoints.push_back( & (*pVecIt).myInitXYZ ); - - return ( thePoints.size() > 0 ); + 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 ] ); } + //======================================================================= -//function : getShapePoints -//purpose : return list of points located on theShape +//function : arrangeBoundaries +//purpose : if there are several wires, arrange boundaryPoints so that +// the outer wire goes first and fix inner wires orientation +// update myKeyPointIDs to correspond to the order of key-points +// in boundaries; sort internal boundaries by the nb of key-points //======================================================================= -list< SMESH_Pattern::TPoint* > & - SMESH_Pattern::getShapePoints(const TopoDS_Shape& theShape) +void SMESH_Pattern::arrangeBoundaries (list< list< TPoint* > >& boundaryList) { - int aShapeID; - if ( !myShapeIDMap.Contains( theShape )) - aShapeID = myShapeIDMap.Add( theShape ); - else - aShapeID = myShapeIDMap.FindIndex( theShape ); + typedef list< list< TPoint* > >::iterator TListOfListIt; + TListOfListIt bndIt; + list< TPoint* >::iterator pIt; - return myShapeIDToPointsMap[ aShapeID ]; -} - -//======================================================================= -//function : getShapePoints -//purpose : return list of points located on the shape -//======================================================================= - -list< SMESH_Pattern::TPoint* > & SMESH_Pattern::getShapePoints(const int theShapeID) -{ - return myShapeIDToPointsMap[ theShapeID ]; -} - -//======================================================================= -//function : DumpPoints -//purpose : Debug -//======================================================================= - -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; -#endif -} - -//======================================================================= -//function : TPoint() -//purpose : -//======================================================================= - -SMESH_Pattern::TPoint::TPoint() -{ -#ifdef _DEBUG_ - myInitXYZ.SetCoord(0,0,0); - myInitUV.SetCoord(0.,0.); - myInitU = 0; - myXYZ.SetCoord(0,0,0); - myUV.SetCoord(0.,0.); - myU = 0; -#endif -} - -//======================================================================= -//function : operator << -//purpose : -//======================================================================= - -ostream & operator <<(ostream & OS, const SMESH_Pattern::TPoint& p) -{ - gp_XYZ xyz = p.myInitXYZ; - OS << "\tinit( xyz( " << xyz.X() << " " << xyz.Y() << " " << xyz.Z() << " )"; - gp_XY xy = p.myInitUV; - OS << " uv( " << xy.X() << " " << xy.Y() << " )"; - double u = p.myInitU; - OS << " u( " << u << " )) " << &p << endl; - xyz = p.myXYZ.XYZ(); - OS << "\t ( xyz( " << xyz.X() << " " << xyz.Y() << " " << xyz.Z() << " )"; - xy = p.myUV; - OS << " uv( " << xy.X() << " " << xy.Y() << " )"; - u = p.myU; - OS << " u( " << u << " ))" << endl; - - return OS; -} - -//======================================================================= -//function : TBlock::TEdge::GetU -//purpose : -//======================================================================= - -double TBlock::TEdge::GetU( const gp_XYZ& theParams ) const -{ - double u = theParams.Coord( myCoordInd ); - return ( 1 - u ) * myFirst + u * myLast; -} - -//======================================================================= -//function : TBlock::TEdge::Point -//purpose : -//======================================================================= - -gp_XYZ TBlock::TEdge::Point( const gp_XYZ& theParams ) const -{ - gp_XYZ p = myC3d->Value( GetU( theParams )).XYZ(); - if ( myTrsf.Form() != gp_Identity ) - myTrsf.Transforms( p ); - return p; -} - -//======================================================================= -//function : TBlock::TFace::GetUV -//purpose : -//======================================================================= - -gp_XY TBlock::TFace::GetUV( const gp_XYZ& theParams ) const -{ - gp_XY uv(0.,0.); - double dU = theParams.Coord( GetUInd() ); - double dV = theParams.Coord( GetVInd() ); - for ( int iE = 0; iE < 4; iE++ ) // loop on 4 edges + int nbBoundaries = boundaryList.size(); + if ( nbBoundaries > 1 ) { - double Ecoef = 0, Vcoef = 0; - switch ( iE ) { - case 0: - Ecoef = ( 1 - dV ); // u0 - Vcoef = ( 1 - dU ) * ( 1 - dV ); break; // 00 - case 1: - Ecoef = dV; // u1 - Vcoef = dU * ( 1 - dV ); break; // 10 - case 2: - Ecoef = ( 1 - dU ); // 0v - Vcoef = dU * dV ; break; // 11 - case 3: - Ecoef = dU ; // 1v - Vcoef = ( 1 - dU ) * dV ; break; // 01 - default:; + // sort boundaries by nb of key-points + if ( nbBoundaries > 2 ) + { + // move boundaries in tmp list + 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 + typedef map< int, TListOfListIt > TNbKpBndPosMap; + TNbKpBndPosMap nbKpBndPosMap; + bndIt = tmpList.begin(); + list< int >::iterator nbKpIt = myNbKeyPntInBoundary.begin(); + for ( ; nbKpIt != myNbKeyPntInBoundary.end(); nbKpIt++, bndIt++ ) { + int nb = *nbKpIt * nbBoundaries; + while ( nbKpBndPosMap.find ( nb ) != nbKpBndPosMap.end() ) + nb++; + nbKpBndPosMap.insert( TNbKpBndPosMap::value_type( nb, bndIt )); + } + // move boundaries back to boundaryList + TNbKpBndPosMap::iterator nbKpBndPosIt = nbKpBndPosMap.begin(); + for ( ; nbKpBndPosIt != nbKpBndPosMap.end(); nbKpBndPosIt++ ) { + TListOfListIt & bndPos2 = (*nbKpBndPosIt).second; + TListOfListIt bndPos1 = bndPos2++; + boundaryList.splice( boundaryList.end(), tmpList, bndPos1, bndPos2 ); + } } - // edge addition - double u = theParams.Coord( myCoordInd[ iE ] ); - u = ( 1 - u ) * myFirst[ iE ] + u * myLast[ iE ]; - uv += Ecoef * myC2d[ iE ]->Value( u ).XY(); - // corner addition - uv -= Vcoef * myCorner[ iE ]; - } - return uv; -} -//======================================================================= -//function : TBlock::TFace::Point -//purpose : -//======================================================================= + // Look for the outer boundary: the one with the point with the least X + double leastX = DBL_MAX; + TListOfListIt outerBndPos; + for ( bndIt = boundaryList.begin(); bndIt != boundaryList.end(); bndIt++ ) + { + list< TPoint* >& boundary = (*bndIt); + for ( pIt = boundary.begin(); pIt != boundary.end(); pIt++) + { + TPoint* point = *pIt; + if ( point->myInitXYZ.X() < leastX ) { + leastX = point->myInitXYZ.X(); + outerBndPos = bndIt; + } + } + } -gp_XYZ TBlock::TFace::Point( const gp_XYZ& theParams ) const -{ - gp_XY uv = GetUV( theParams ); - gp_XYZ p = myS->Value( uv.X(), uv.Y() ).XYZ(); - if ( myTrsf.Form() != gp_Identity ) - myTrsf.Transforms( p ); - return p; -} + if ( outerBndPos != boundaryList.begin() ) + boundaryList.splice( boundaryList.begin(), boundaryList, outerBndPos, ++outerBndPos ); -//======================================================================= -//function : GetShapeCoef -//purpose : -//======================================================================= + } // if nbBoundaries > 1 + + // Check boundaries orientation and re-fill myKeyPointIDs -double* TBlock::GetShapeCoef (const int theShapeID) -{ - static double shapeCoef[][3] = { - // V000, V100, V010, V110 - { -1,-1,-1 }, { 1,-1,-1 }, { -1, 1,-1 }, { 1, 1,-1 }, - // V001, V101, V011, V111, - { -1,-1, 1 }, { 1,-1, 1 }, { -1, 1, 1 }, { 1, 1, 1 }, - // Ex00, Ex10, Ex01, Ex11, - { 0,-1,-1 }, { 0, 1,-1 }, { 0,-1, 1 }, { 0, 1, 1 }, - // E0y0, E1y0, E0y1, E1y1, - { -1, 0,-1 }, { 1, 0,-1 }, { -1, 0, 1 }, { 1, 0, 1 }, - // E00z, E10z, E01z, E11z, - { -1,-1, 0 }, { 1,-1, 0 }, { -1, 1, 0 }, { 1, 1, 0 }, - // Fxy0, Fxy1, Fx0z, Fx1z, F0yz, F1yz, - { 0, 0,-1 }, { 0, 0, 1 }, { 0,-1, 0 }, { 0, 1, 0 }, { -1, 0, 0 }, { 1, 0, 0 }, - // ID_Shell - { 0, 0, 0 } - }; - if ( theShapeID < ID_V000 || theShapeID > ID_F1yz ) - return shapeCoef[ ID_Shell - 1 ]; - - return shapeCoef[ theShapeID - 1 ]; -} + set< TPoint* > keyPointSet; + list< int >::iterator kpIt = myKeyPointIDs.begin(); + for ( ; kpIt != myKeyPointIDs.end(); kpIt++ ) + keyPointSet.insert( & myPoints[ *kpIt ]); + myKeyPointIDs.clear(); -//======================================================================= -//function : ShellPoint -//purpose : return coordinates of a point in shell -//======================================================================= + // update myNbKeyPntInBoundary also + list< int >::iterator nbKpIt = myNbKeyPntInBoundary.begin(); -bool TBlock::ShellPoint( const gp_XYZ& theParams, gp_XYZ& thePoint ) const -{ - thePoint.SetCoord( 0., 0., 0. ); - for ( int shapeID = ID_V000; shapeID < ID_Shell; shapeID++ ) + for ( bndIt = boundaryList.begin(); bndIt != boundaryList.end(); bndIt++, nbKpIt++ ) { - // coef - double* coefs = GetShapeCoef( shapeID ); - double k = 1; - for ( int iCoef = 0; iCoef < 3; iCoef++ ) { - if ( coefs[ iCoef ] != 0 ) { - if ( coefs[ iCoef ] < 0 ) - k *= ( 1. - theParams.Coord( iCoef + 1 )); - else - k *= theParams.Coord( iCoef + 1 ); + // find the point with the least X + double leastX = DBL_MAX; + list< TPoint* >::iterator xpIt; + list< TPoint* >& boundary = (*bndIt); + for ( pIt = boundary.begin(); pIt != boundary.end(); pIt++) + { + TPoint* point = *pIt; + if ( point->myInitXYZ.X() < leastX ) { + leastX = point->myInitXYZ.X(); + xpIt = pIt; } } - // point on a shape - gp_XYZ Ps; - if ( shapeID < ID_Ex00 ) // vertex - VertexPoint( shapeID, Ps ); - else if ( shapeID < ID_Fxy0 ) { // edge - EdgePoint( shapeID, theParams, Ps ); - k = -k; - } else // face - FacePoint( shapeID, theParams, Ps ); - - thePoint += k * Ps; - } - return true; -} + // find points next to the point with the least X + TPoint* p = *xpIt, *pPrev, *pNext; + if ( p == boundary.front() ) + pPrev = *(++boundary.rbegin()); + else { + xpIt--; + pPrev = *xpIt; + xpIt++; + } + if ( p == boundary.back() ) + pNext = *(++boundary.begin()); + else { + xpIt++; + pNext = *xpIt; + } + // vectors of boundary direction near

+ gp_Vec2d v1( pPrev->myInitUV, p->myInitUV ), v2( p->myInitUV, pNext->myInitUV ); + double sqMag1 = v1.SquareMagnitude(), sqMag2 = v2.SquareMagnitude(); + if ( sqMag1 > DBL_MIN && sqMag2 > DBL_MIN ) { + double yPrev = v1.Y() / sqrt( sqMag1 ); + double yNext = v2.Y() / sqrt( sqMag2 ); + double sumY = yPrev + yNext; + bool reverse; + if ( bndIt == boundaryList.begin() ) // outer boundary + reverse = sumY > 0; + else + reverse = sumY < 0; + if ( reverse ) + boundary.reverse(); + } -//======================================================================= -//function : NbVariables -//purpose : -//======================================================================= + // Put key-point IDs of a well-oriented boundary in myKeyPointIDs + (*nbKpIt) = 0; // count nb of key-points again + pIt = boundary.begin(); + for ( ; pIt != boundary.end(); pIt++) + { + TPoint* point = *pIt; + if ( keyPointSet.find( point ) == keyPointSet.end() ) + continue; + // find an index of a keypoint + int index = 0; + vector< TPoint >::const_iterator pVecIt = myPoints.begin(); + for ( ; pVecIt != myPoints.end(); pVecIt++, index++ ) + if ( &(*pVecIt) == point ) + break; + myKeyPointIDs.push_back( index ); + (*nbKpIt)++; + } + myKeyPointIDs.pop_back(); // remove the first key-point from the back + (*nbKpIt)--; -Standard_Integer TBlock::NbVariables() const -{ - return 3; + } // loop on a list of boundaries + + ASSERT( myKeyPointIDs.size() == keyPointSet.size() ); } //======================================================================= -//function : NbEquations -//purpose : +//function : findBoundaryPoints +//purpose : if loaded from file, find points to map on edges and faces and +// compute their parameters //======================================================================= -Standard_Integer TBlock::NbEquations() const +bool SMESH_Pattern::findBoundaryPoints() { - return 1; -} + if ( myIsBoundaryPointsFound ) return true; -//======================================================================= -//function : Value -//purpose : -//======================================================================= + MESSAGE(" findBoundaryPoints() "); -Standard_Boolean TBlock::Value(const math_Vector& theXYZ, math_Vector& theFxyz) -{ - gp_XYZ P, params( theXYZ(1), theXYZ(2), theXYZ(3) ); - if ( params.IsEqual( myParam, DBL_MIN )) { // same param - theFxyz( 1 ) = myValues[ 0 ]; - } - else { - ShellPoint( params, P ); - gp_Vec dP( P - myPoint ); - theFxyz(1) = SQRT_FUNC ? dP.SquareMagnitude() : dP.Magnitude(); - } - return true; -} + if ( myIs2D ) + { + set< TPoint* > pointsInElems; -//======================================================================= -//function : Derivatives -//purpose : -//======================================================================= + // Find free links of elements: + // put links of all elements in a set and remove links encountered twice -Standard_Boolean TBlock::Derivatives(const math_Vector& XYZ,math_Matrix& Df) -{ - MESSAGE( "TBlock::Derivatives()"); - math_Vector F(1,3); - return Values(XYZ,F,Df); -} + typedef pair< TPoint*, TPoint*> TLink; + set< TLink > linkSet; + list::iterator epIt = myElemPointIDs.begin(); + for ( ; epIt != myElemPointIDs.end(); epIt++ ) + { + TElemDef & elemPoints = *epIt; + TElemDef::iterator pIt = elemPoints.begin(); + int prevP = elemPoints.back(); + for ( ; pIt != elemPoints.end(); pIt++ ) { + TPoint* p1 = & myPoints[ prevP ]; + TPoint* p2 = & myPoints[ *pIt ]; + TLink link(( p1 < p2 ? p1 : p2 ), ( p1 < p2 ? p2 : p1 )); + ASSERT( link.first != link.second ); + pair::iterator,bool> itUniq = linkSet.insert( link ); + if ( !itUniq.second ) + linkSet.erase( itUniq.first ); + prevP = *pIt; -//======================================================================= -//function : Values -//purpose : -//======================================================================= + pointsInElems.insert( p1 ); + } + } + // Now linkSet contains only free links, + // find the points order that they have in boundaries -Standard_Boolean TBlock::Values(const math_Vector& theXYZ, - math_Vector& theFxyz, - math_Matrix& theDf) -{ -// MESSAGE( endl<<"TBlock::Values( "< keyPointSet; + list< int >::iterator kpIt = myKeyPointIDs.begin(); + for ( ; kpIt != myKeyPointIDs.end(); kpIt++ ) + keyPointSet.insert( & myPoints[ *kpIt ]); - ShellPoint( params, P ); - //myNbIterations++; // how many time call ShellPoint() - - gp_Vec dP( P - myPoint ); - theFxyz(1) = SQRT_FUNC ? dP.SquareMagnitude() : dP.Magnitude(); - if ( theFxyz(1) < 1e-6 ) { - myParam = params; - myValues[ 0 ]= 0; - theDf( 1,1 ) = 0; - theDf( 1,2 ) = 0; - theDf( 1,3 ) = 0; - return true; - } + // 2. chain up boundary points + list< list< TPoint* > > boundaryList; + boundaryList.push_back( list< TPoint* >() ); + list< TPoint* > * boundary = & boundaryList.back(); - if ( theFxyz(1) < myValues[0] ) - { - // 3 partial derivatives - gp_Vec drv[ 3 ]; - for ( int iP = 1; iP <= 3; iP++ ) { - gp_XYZ Pi; - params.SetCoord( iP, theXYZ( iP ) + 0.001 ); - ShellPoint( params, Pi ); - params.SetCoord( iP, theXYZ( iP ) ); // restore params - gp_Vec dPi ( P, Pi ); - double mag = dPi.Magnitude(); - if ( mag > DBL_MIN ) - dPi /= mag; - drv[ iP - 1 ] = dPi; - } - for ( int iP = 0; iP < 3; iP++ ) { - if ( iP == myFaceIndex ) - theDf( 1, iP + 1 ) = myFaceParam; + TPoint *point1, *point2, *keypoint1; + kpIt = myKeyPointIDs.begin(); + point1 = keypoint1 = & myPoints[ *kpIt++ ]; + // loop on free links: look for the next point + int iKeyPoint = 0; + set< TLink >::iterator lIt = linkSet.begin(); + while ( lIt != linkSet.end() ) + { + if ( (*lIt).first == point1 ) + point2 = (*lIt).second; + else if ( (*lIt).second == point1 ) + point2 = (*lIt).first; else { - // like IntAna_IntConicQuad::Perform (const gp_Lin& L, const gp_Pln& P) - // where L is (P -> myPoint), P is defined by the 2 other derivative direction - int iPrev = ( iP ? iP - 1 : 2 ); - int iNext = ( iP == 2 ? 0 : iP + 1 ); - gp_Vec plnNorm = drv[ iPrev ].Crossed( drv [ iNext ] ); - double Direc = plnNorm * drv[ iP ]; - if ( Abs(Direc) <= DBL_MIN ) - theDf( 1, iP + 1 ) = dP * drv[ iP ]; - else { - double Dis = plnNorm * P - plnNorm * myPoint; - theDf( 1, iP + 1 ) = Dis/Direc; - } + lIt++; + continue; } - } - //myNbIterations +=3; // how many time call ShellPoint() - - // store better values - myParam = params; - myValues[0]= theFxyz(1); - myValues[1]= theDf(1,1); - myValues[2]= theDf(1,2); - myValues[3]= theDf(1,3); - -// SCRUTE( theFxyz(1) ); -// SCRUTE( theDf( 1,1 )); -// SCRUTE( theDf( 1,2 )); -// SCRUTE( theDf( 1,3 )); - } + linkSet.erase( lIt ); + lIt = linkSet.begin(); - return true; -} + if ( keyPointSet.find( point2 ) == keyPointSet.end() ) // not a key-point + { + boundary->push_back( point2 ); + } + else // a key-point found + { + keyPointSet.erase( point2 ); // keyPointSet contains not found key-points only + iKeyPoint++; + if ( point2 != keypoint1 ) // its not the boundary end + { + boundary->push_back( point2 ); + } + else // the boundary end reached + { + boundary->push_front( keypoint1 ); + boundary->push_back( keypoint1 ); + myNbKeyPntInBoundary.push_back( iKeyPoint ); + if ( keyPointSet.empty() ) + break; // all boundaries containing key-points are found -//======================================================================= -//function : ComputeParameters -//purpose : compute point parameters in the block -//======================================================================= + // prepare to search for the next boundary + boundaryList.push_back( list< TPoint* >() ); + boundary = & boundaryList.back(); + point2 = keypoint1 = (*keyPointSet.begin()); + } + } + point1 = point2; + } // loop on the free links set -bool TBlock::ComputeParameters(const gp_Pnt& thePoint, - gp_XYZ& theParams, - const int theShapeID) -{ -// MESSAGE( endl<<"TBlock::ComputeParameters( " -// <empty() ) { + MESSAGE(" a separate key-point"); + return setErrorCode( ERR_READ_BAD_KEY_POINT ); + } - myPoint = thePoint.XYZ(); + // if there are several wires, arrange boundaryPoints so that + // the outer wire goes first and fix inner wires orientation; + // sort myKeyPointIDs to correspond to the order of key-points + // in boundaries + arrangeBoundaries( boundaryList ); - myParam.SetCoord( -1,-1,-1 ); - myValues[0] = 1e100; + // Find correspondence shape ID - points, + // compute points parameter on edge - const bool isOnFace = IsFaceID( theShapeID ); - double * coef = GetShapeCoef( theShapeID ); + keyPointSet.clear(); + for ( kpIt = myKeyPointIDs.begin(); kpIt != myKeyPointIDs.end(); kpIt++ ) + keyPointSet.insert( & myPoints[ *kpIt ]); - // the first guess - math_Vector start( 1, 3, 0.0 ); - if ( !myGridComputed ) - { - // define the first guess by thePoint projection on lines - // connecting vertices - bool needGrid = false; - gp_XYZ par000( 0, 0, 0 ), par111( 1, 1, 1 ); - double zero = DBL_MIN * DBL_MIN; - for ( int iEdge = 0, iParam = 1; iParam <= 3 && !needGrid; iParam++ ) + set< TPoint* > edgePointSet; // to find in-face points + int vertexID = 1; // the first index in TopTools_IndexedMapOfShape + int edgeID = myKeyPointIDs.size() + 1; + + list< list< TPoint* > >::iterator bndIt = boundaryList.begin(); + for ( ; bndIt != boundaryList.end(); bndIt++ ) { - if ( isOnFace && coef[ iParam - 1 ] != 0 ) { - iEdge += 4; - continue; - } - for ( int iE = 0; iE < 4; iE++, iEdge++ ) { // loop on 4 parallel edges - gp_Pnt p0 = myEdge[ iEdge ].Point( par000 ); - gp_Pnt p1 = myEdge[ iEdge ].Point( par111 ); - gp_Vec v01( p0, p1 ), v0P( p0, thePoint ); - double len2 = v01.SquareMagnitude(); - double par = 0; - if ( len2 > zero ) { - par = v0P.Dot( v01 ) / len2; - if ( par < 0 || par > 1 ) { - needGrid = true; - break; + boundary = & (*bndIt); + double edgeLength = 0; + list< TPoint* >::iterator pIt = boundary->begin(); + getShapePoints( edgeID ).push_back( *pIt ); + getShapePoints( vertexID++ ).push_back( *pIt ); + for ( pIt++; pIt != boundary->end(); pIt++) + { + list< TPoint* > & edgePoints = getShapePoints( edgeID ); + TPoint* prevP = edgePoints.empty() ? 0 : edgePoints.back(); + TPoint* point = *pIt; + edgePointSet.insert( point ); + if ( keyPointSet.find( point ) == keyPointSet.end() ) // inside-edge point + { + edgePoints.push_back( point ); + edgeLength += ( point->myInitUV - prevP->myInitUV ).Modulus(); + point->myInitU = edgeLength; + } + else // a key-point + { + // treat points on the edge which ends up: compute U [0,1] + edgePoints.push_back( point ); + if ( edgePoints.size() > 2 ) { + edgeLength += ( point->myInitUV - prevP->myInitUV ).Modulus(); + list< TPoint* >::iterator epIt = edgePoints.begin(); + for ( ; epIt != edgePoints.end(); epIt++ ) + (*epIt)->myInitU /= edgeLength; + } + // begin the next edge treatment + edgeLength = 0; + edgeID++; + if ( point != boundary->front() ) { // not the first key-point again + getShapePoints( edgeID ).push_back( point ); + getShapePoints( vertexID++ ).push_back( point ); } } - start( iParam ) += par; - } - start( iParam ) /= 4.; - } - if ( needGrid ) { - // compute nodes of 3 x 3 x 3 grid - int iNode = 0; - for ( double x = 0.25; x < 0.9; x += 0.25 ) - for ( double y = 0.25; y < 0.9; y += 0.25 ) - for ( double z = 0.25; z < 0.9; z += 0.25 ) { - TxyzPair & prmPtn = my3x3x3GridNodes[ iNode++ ]; - prmPtn.first.SetCoord( x, y, z ); - ShellPoint( prmPtn.first, prmPtn.second ); - } - myGridComputed = true; - } - } - if ( myGridComputed ) { - double minDist = DBL_MAX; - gp_XYZ* bestParam = 0; - for ( int iNode = 0; iNode < 27; iNode++ ) { - TxyzPair & prmPtn = my3x3x3GridNodes[ iNode ]; - double dist = ( thePoint.XYZ() - prmPtn.second ).SquareModulus(); - if ( dist < minDist ) { - minDist = dist; - bestParam = & prmPtn.first; } } - start( 1 ) = bestParam->X(); - start( 2 ) = bestParam->Y(); - start( 3 ) = bestParam->Z(); - } - int myFaceIndex = -1; - if ( isOnFace ) { - // put a point on the face - for ( int iCoord = 0; iCoord < 3; iCoord++ ) - if ( coef[ iCoord ] ) { - myFaceIndex = iCoord; - myFaceParam = ( coef[ myFaceIndex ] < 0.5 ) ? 0.0 : 1.0; - start( iCoord + 1 ) = myFaceParam; - } - } - math_Vector low ( 1, 3, 0.0 ); - math_Vector up ( 1, 3, 1.0 ); - math_Vector tol ( 1, 3, 1e-4 ); - math_FunctionSetRoot paramSearch( *this, tol ); - - int nbLoops = 0; - while ( myValues[0] > 1e-1 && nbLoops++ < 10 ) { - paramSearch.Perform ( *this, start, low, up ); - if ( !paramSearch.IsDone() ) { - //MESSAGE( " !paramSearch.IsDone() " ); - } - else { - //MESSAGE( " NB ITERATIONS: " << paramSearch.NbIterations() ); + // find in-face points + list< TPoint* > & facePoints = getShapePoints( edgeID ); + vector< TPoint >::iterator pVecIt = myPoints.begin(); + for ( ; pVecIt != myPoints.end(); pVecIt++ ) { + TPoint* point = &(*pVecIt); + if ( edgePointSet.find( point ) == edgePointSet.end() && + pointsInElems.find( point ) != pointsInElems.end()) + facePoints.push_back( point ); } - start( 1 ) = myParam.X(); - start( 2 ) = myParam.Y(); - start( 3 ) = myParam.Z(); - //MESSAGE( "Distance: " << ( SQRT_FUNC ? sqrt(myValues[0]) : myValues[0] )); - } -// MESSAGE( endl << myParam.X() << " " << myParam.Y() << " " << myParam.Z() << endl); -// mySumDist += myValues[0]; -// MESSAGE( " TOTAL NB ITERATIONS: " << myNbIterations << -// " DIST: " << ( SQRT_FUNC ? sqrt(mySumDist) : mySumDist )); + } // 2D case - theParams = myParam; + else // 3D case + { + // bind points to shapes according to point parameters + vector< TPoint >::iterator pVecIt = myPoints.begin(); + for ( int i = 0; pVecIt != myPoints.end(); pVecIt++, i++ ) { + TPoint* point = &(*pVecIt); + int shapeID = SMESH_Block::GetShapeIDByParams( point->myInitXYZ ); + getShapePoints( shapeID ).push_back( point ); + // detect key-points + if ( SMESH_Block::IsVertexID( shapeID )) + myKeyPointIDs.push_back( i ); + } + } - return true; + myIsBoundaryPointsFound = true; + return myIsBoundaryPointsFound; } //======================================================================= -//function : GetStateNumber -//purpose : +//function : Clear +//purpose : clear fields //======================================================================= -Standard_Integer TBlock::GetStateNumber () +void SMESH_Pattern::Clear() { -// MESSAGE( endl<<"TBlock::GetStateNumber( "< 26 || id < 0 ) { - MESSAGE( "GetShapeIDByParams() = " << id - <<" "<< theCoord.X() <<" "<< theCoord.Y() <<" "<< theCoord.Z() ); + + // check nb of vertices + TopTools_IndexedMapOfShape vMap; + TopExp::MapShapes( theShape, TopAbs_VERTEX, vMap ); + if ( vMap.Extent() + nbNodeOnSeamEdge != myKeyPointIDs.size() ) { + MESSAGE( myKeyPointIDs.size() << " != " << vMap.Extent() ); + return setErrorCode( ERR_APPL_BAD_NB_VERTICES ); } - return id + 1; // shape ids start at 1 + myElements.clear(); // not refine elements + myElemXYZIDs.clear(); + + myShapeIDMap.Clear(); + myShape = theShape; + return true; } //======================================================================= -//function : LoadBlockShapes -//purpose : add sub-shapes of theBlock to theShapeIDMap so that they get -// IDs acoording to enum TBlockShapeID +//function : GetMappedPoints +//purpose : Return nodes coordinates computed by Apply() method //======================================================================= -bool TBlock::LoadBlockShapes(const TopoDS_Vertex& theVertex000, - const TopoDS_Vertex& theVertex001, -// TopTools_IndexedMapOfShape& theShapeIDMap - TopTools_IndexedMapOfOrientedShape& theShapeIDMap ) +bool SMESH_Pattern::GetMappedPoints ( list< const gp_XYZ * > & thePoints ) const { - MESSAGE(" ::LoadBlockShapes()"); - myGridComputed = false; - - // 6 vertices - TopoDS_Shape V000, V100, V010, V110, V001, V101, V011, V111; - // 12 edges - TopoDS_Shape Ex00, Ex10, Ex01, Ex11; - TopoDS_Shape E0y0, E1y0, E0y1, E1y1; - TopoDS_Shape E00z, E10z, E01z, E11z; - // 6 faces - TopoDS_Shape Fxy0, Fx0z, F0yz, Fxy1, Fx1z, F1yz; - - // nb of faces bound to a vertex in TopTools_IndexedDataMapOfShapeListOfShape - // filled by TopExp::MapShapesAndAncestors() - const int NB_FACES_BY_VERTEX = 6; - - TopTools_IndexedDataMapOfShapeListOfShape vfMap; - TopExp::MapShapesAndAncestors( myShell, TopAbs_VERTEX, TopAbs_FACE, vfMap ); - if ( vfMap.Extent() != 8 ) { - MESSAGE(" Wrong nb of vertices in the block: " << vfMap.Extent() ); + thePoints.clear(); + if ( !myIsComputed ) return false; - } - - V000 = theVertex000; - V001 = theVertex001; - - if ( V000.IsNull() ) { - // find vertex 000 - the one with smallest coordinates - double minVal = DBL_MAX, minX, val; - for ( int i = 1; i <= 8; i++ ) { - const TopoDS_Vertex& v = TopoDS::Vertex( vfMap.FindKey( i )); - gp_Pnt P = BRep_Tool::Pnt( v ); - val = P.X() + P.Y() + P.Z(); - if ( val < minVal || ( val == minVal && P.X() < minX )) { - V000 = v; - minVal = val; - minX = P.X(); - } - } - // find vertex 001 - the one on the most vertical edge passing through V000 - TopTools_IndexedDataMapOfShapeListOfShape veMap; - TopExp::MapShapesAndAncestors( myShell, TopAbs_VERTEX, TopAbs_EDGE, veMap ); - gp_Vec dir001 = gp::DZ(); - gp_Pnt p000 = BRep_Tool::Pnt( TopoDS::Vertex( V000 )); - double maxVal = -DBL_MAX; - TopTools_ListIteratorOfListOfShape eIt ( veMap.FindFromKey( V000 )); - for ( ; eIt.More(); eIt.Next() ) { - const TopoDS_Edge& e = TopoDS::Edge( eIt.Value() ); - TopoDS_Vertex v = TopExp::FirstVertex( e ); - if ( v.IsSame( V000 )) - v = TopExp::LastVertex( e ); - val = dir001 * gp_Vec( p000, BRep_Tool::Pnt( v )).Normalized(); - if ( val > maxVal ) { - V001 = v; - maxVal = val; - } - } - } - - // find the bottom (Fxy0), Fx0z and F0yz faces - const TopTools_ListOfShape& f000List = vfMap.FindFromKey( V000 ); - const TopTools_ListOfShape& f001List = vfMap.FindFromKey( V001 ); - if (f000List.Extent() != NB_FACES_BY_VERTEX || - f001List.Extent() != NB_FACES_BY_VERTEX ) { - MESSAGE(" LoadBlockShapes() " << f000List.Extent() << " " << f001List.Extent()); - return false; - } - TopTools_ListIteratorOfListOfShape f001It, f000It ( f000List ); - int i, j, iFound1, iFound2; - for ( j = 0; f000It.More(); f000It.Next(), j++ ) - { - if ( NB_FACES_BY_VERTEX == 6 && j % 2 ) continue; // each face encounters twice - const TopoDS_Shape& F = f000It.Value(); - for ( i = 0, f001It.Initialize( f001List ); f001It.More(); f001It.Next(), i++ ) { - if ( NB_FACES_BY_VERTEX == 6 && i % 2 ) continue; // each face encounters twice - if ( F.IsSame( f001It.Value() )) - break; - } - if ( f001It.More() ) // Fx0z or F0yz found - if ( Fx0z.IsNull() ) { - Fx0z = F; - iFound1 = i; - } else { - F0yz = F; - iFound2 = i; - } - else // F is the bottom face - Fxy0 = F; + if ( myElements.empty() ) { // applied to shape + vector< TPoint >::const_iterator pVecIt = myPoints.begin(); + for ( ; pVecIt != myPoints.end(); pVecIt++ ) + thePoints.push_back( & (*pVecIt).myXYZ.XYZ() ); } - if ( Fxy0.IsNull() || Fx0z.IsNull() || F0yz.IsNull() ) { - MESSAGE( Fxy0.IsNull() <<" "<< Fx0z.IsNull() <<" "<< F0yz.IsNull() ); - return false; + 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 ) + if ( !isDefined( *xyz )) + thePoints.push_back( definedXYZ ); + else + thePoints.push_back( & (*xyz) ); } + return !thePoints.empty(); +} - // choose the top face (Fxy1) - for ( i = 0, f001It.Initialize( f001List ); f001It.More(); f001It.Next(), i++ ) { - if ( NB_FACES_BY_VERTEX == 6 && i % 2 ) continue; // each face encounters twice - if ( i != iFound1 && i != iFound2 ) - break; - } - Fxy1 = f001It.Value(); - if ( Fxy1.IsNull() ) { - MESSAGE(" LoadBlockShapes() error "); - return false; - } - // find bottom edges and veritices - list< TopoDS_Edge > eList; - list< int > nbVertexInWires; - getOrderedEdges( TopoDS::Face( Fxy0 ), TopoDS::Vertex( V000 ), eList, nbVertexInWires ); - if ( nbVertexInWires.size() != 1 || nbVertexInWires.front() != 4 ) { - MESSAGE(" LoadBlockShapes() error "); - return false; - } - list< TopoDS_Edge >::iterator elIt = eList.begin(); - for ( i = 0; elIt != eList.end(); elIt++, i++ ) - switch ( i ) { - case 0: E0y0 = *elIt; V010 = TopExp::LastVertex( *elIt, true ); break; - case 1: Ex10 = *elIt; V110 = TopExp::LastVertex( *elIt, true ); break; - case 2: E1y0 = *elIt; V100 = TopExp::LastVertex( *elIt, true ); break; - case 3: Ex00 = *elIt; break; - default:; - } - if ( i != 4 || E0y0.IsNull() || Ex10.IsNull() || E1y0.IsNull() || Ex00.IsNull() ) { - MESSAGE(" LoadBlockShapes() error, eList.size()=" << eList.size()); - return false; - } +//======================================================================= +//function : GetPoints +//purpose : Return nodes coordinates of the pattern +//======================================================================= +bool SMESH_Pattern::GetPoints ( list< const gp_XYZ * > & thePoints ) const +{ + thePoints.clear(); - // find top edges and veritices - eList.clear(); - getOrderedEdges( TopoDS::Face( Fxy1 ), TopoDS::Vertex( V001 ), eList, nbVertexInWires ); - if ( nbVertexInWires.size() != 1 || nbVertexInWires.front() != 4 ) { - MESSAGE(" LoadBlockShapes() error "); - return false; - } - for ( i = 0, elIt = eList.begin(); elIt != eList.end(); elIt++, i++ ) - switch ( i ) { - case 0: Ex01 = *elIt; V101 = TopExp::LastVertex( *elIt, true ); break; - case 1: E1y1 = *elIt; V111 = TopExp::LastVertex( *elIt, true ); break; - case 2: Ex11 = *elIt; V011 = TopExp::LastVertex( *elIt, true ); break; - case 3: E0y1 = *elIt; break; - default:; - } - if ( i != 4 || Ex01.IsNull() || E1y1.IsNull() || Ex11.IsNull() || E0y1.IsNull() ) { - MESSAGE(" LoadBlockShapes() error, eList.size()=" << eList.size()); + if ( !IsLoaded() ) return false; - } - // swap Fx0z and F0yz if necessary - TopExp_Explorer exp( Fx0z, TopAbs_VERTEX ); - for ( ; exp.More(); exp.Next() ) // Fx0z shares V101 and V100 - if ( V101.IsSame( exp.Current() ) || V100.IsSame( exp.Current() )) - break; // V101 or V100 found - if ( !exp.More() ) { // not found - TopoDS_Shape f = Fx0z; Fx0z = F0yz; F0yz = f; - } - - // find Fx1z and F1yz faces - const TopTools_ListOfShape& f111List = vfMap.FindFromKey( V111 ); - const TopTools_ListOfShape& f110List = vfMap.FindFromKey( V110 ); - if (f111List.Extent() != NB_FACES_BY_VERTEX || - f110List.Extent() != NB_FACES_BY_VERTEX ) { - MESSAGE(" LoadBlockShapes() " << f111List.Extent() << " " << f110List.Extent()); - return false; - } - TopTools_ListIteratorOfListOfShape f111It, f110It ( f110List); - for ( j = 0 ; f110It.More(); f110It.Next(), j++ ) { - if ( NB_FACES_BY_VERTEX == 6 && j % 2 ) continue; // each face encounters twice - const TopoDS_Shape& F = f110It.Value(); - for ( i = 0, f111It.Initialize( f111List ); f111It.More(); f111It.Next(), i++ ) { - if ( NB_FACES_BY_VERTEX == 6 && i % 2 ) continue; // each face encounters twice - if ( F.IsSame( f111It.Value() )) { // Fx1z or F1yz found - if ( Fx1z.IsNull() ) - Fx1z = F; - else - F1yz = F; - } - } - } - if ( Fx1z.IsNull() || F1yz.IsNull() ) { - MESSAGE(" LoadBlockShapes() error "); - return false; - } + vector< TPoint >::const_iterator pVecIt = myPoints.begin(); + for ( ; pVecIt != myPoints.end(); pVecIt++ ) + thePoints.push_back( & (*pVecIt).myInitXYZ ); - // swap Fx1z and F1yz if necessary - for ( exp.Init( Fx1z, TopAbs_VERTEX ); exp.More(); exp.Next() ) - if ( V010.IsSame( exp.Current() ) || V011.IsSame( exp.Current() )) - break; - if ( !exp.More() ) { - TopoDS_Shape f = Fx1z; Fx1z = F1yz; F1yz = f; - } + return ( thePoints.size() > 0 ); +} - // find vertical edges - for ( exp.Init( Fx0z, TopAbs_EDGE ); exp.More(); exp.Next() ) { - const TopoDS_Edge& edge = TopoDS::Edge( exp.Current() ); - const TopoDS_Shape& vFirst = TopExp::FirstVertex( edge, true ); - if ( vFirst.IsSame( V001 )) - E00z = edge; - else if ( vFirst.IsSame( V100 )) - E10z = edge; - } - if ( E00z.IsNull() || E10z.IsNull() ) { - MESSAGE(" LoadBlockShapes() error "); - return false; - } - for ( exp.Init( Fx1z, TopAbs_EDGE ); exp.More(); exp.Next() ) { - const TopoDS_Edge& edge = TopoDS::Edge( exp.Current() ); - const TopoDS_Shape& vFirst = TopExp::FirstVertex( edge, true ); - if ( vFirst.IsSame( V111 )) - E11z = edge; - else if ( vFirst.IsSame( V010 )) - E01z = edge; - } - if ( E01z.IsNull() || E11z.IsNull() ) { - MESSAGE(" LoadBlockShapes() error "); - return false; - } +//======================================================================= +//function : getShapePoints +//purpose : return list of points located on theShape +//======================================================================= - // load shapes in theShapeIDMap +list< SMESH_Pattern::TPoint* > & + SMESH_Pattern::getShapePoints(const TopoDS_Shape& theShape) +{ + int aShapeID; + if ( !myShapeIDMap.Contains( theShape )) + aShapeID = myShapeIDMap.Add( theShape ); + else + aShapeID = myShapeIDMap.FindIndex( theShape ); - theShapeIDMap.Clear(); - - theShapeIDMap.Add(V000.Oriented( TopAbs_FORWARD )); - theShapeIDMap.Add(V100.Oriented( TopAbs_FORWARD )); - theShapeIDMap.Add(V010.Oriented( TopAbs_FORWARD )); - theShapeIDMap.Add(V110.Oriented( TopAbs_FORWARD )); - theShapeIDMap.Add(V001.Oriented( TopAbs_FORWARD )); - theShapeIDMap.Add(V101.Oriented( TopAbs_FORWARD )); - theShapeIDMap.Add(V011.Oriented( TopAbs_FORWARD )); - theShapeIDMap.Add(V111.Oriented( TopAbs_FORWARD )); - - theShapeIDMap.Add(Ex00); - theShapeIDMap.Add(Ex10); - theShapeIDMap.Add(Ex01); - theShapeIDMap.Add(Ex11); - - theShapeIDMap.Add(E0y0); - theShapeIDMap.Add(E1y0); - theShapeIDMap.Add(E0y1); - theShapeIDMap.Add(E1y1); - - theShapeIDMap.Add(E00z); - theShapeIDMap.Add(E10z); - theShapeIDMap.Add(E01z); - theShapeIDMap.Add(E11z); - - theShapeIDMap.Add(Fxy0); - theShapeIDMap.Add(Fxy1); - theShapeIDMap.Add(Fx0z); - theShapeIDMap.Add(Fx1z); - theShapeIDMap.Add(F0yz); - theShapeIDMap.Add(F1yz); - - theShapeIDMap.Add(myShell); + return myShapeIDToPointsMap[ aShapeID ]; +} - if ( theShapeIDMap.Extent() != 27 ) { - MESSAGE("LoadBlockShapes() " << theShapeIDMap.Extent() ); - return false; - } +//======================================================================= +//function : getShapePoints +//purpose : return list of points located on the shape +//======================================================================= - // store shapes geometry - for ( int shapeID = 1; shapeID < theShapeIDMap.Extent(); shapeID++ ) - { - const TopoDS_Shape& S = theShapeIDMap( shapeID ); - switch ( S.ShapeType() ) - { - case TopAbs_VERTEX: { +list< SMESH_Pattern::TPoint* > & SMESH_Pattern::getShapePoints(const int theShapeID) +{ + return myShapeIDToPointsMap[ theShapeID ]; +} - if ( shapeID > ID_V111 ) { - MESSAGE(" shapeID =" << shapeID ); - return false; - } - myPnt[ shapeID - ID_V000 ] = - BRep_Tool::Pnt( TopoDS::Vertex( S )).XYZ(); - break; - } - case TopAbs_EDGE: { +//======================================================================= +//function : DumpPoints +//purpose : Debug +//======================================================================= - const TopoDS_Edge& edge = TopoDS::Edge( S ); - if ( shapeID < ID_Ex00 || shapeID > ID_E11z || edge.IsNull() ) { - MESSAGE(" shapeID =" << shapeID ); - return false; - } - TEdge& tEdge = myEdge[ shapeID - ID_Ex00 ]; - tEdge.myCoordInd = GetCoordIndOnEdge( shapeID ); - TopLoc_Location loc; - tEdge.myC3d = BRep_Tool::Curve( edge, loc, tEdge.myFirst, tEdge.myLast ); - if ( !IsForwardEdge( edge, theShapeIDMap )) - Swap( tEdge.myFirst, tEdge.myLast ); - tEdge.myTrsf = loc; - break; - } - case TopAbs_FACE: { +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; +#endif +} - const TopoDS_Face& face = TopoDS::Face( S ); - if ( shapeID < ID_Fxy0 || shapeID > ID_F1yz || face.IsNull() ) { - MESSAGE(" shapeID =" << shapeID ); - return false; - } - TFace& tFace = myFace[ shapeID - ID_Fxy0 ]; - // pcurves - vector< int > edgeIdVec(4, -1); - GetFaceEdgesIDs( shapeID, edgeIdVec ); - for ( int iE = 0; iE < 4; iE++ ) // loop on 4 edges - { - const TopoDS_Edge& edge = TopoDS::Edge( theShapeIDMap( edgeIdVec[ iE ])); - tFace.myCoordInd[ iE ] = GetCoordIndOnEdge( edgeIdVec[ iE ] ); - tFace.myC2d[ iE ] = - BRep_Tool::CurveOnSurface( edge, face, tFace.myFirst[iE], tFace.myLast[iE] ); - if ( !IsForwardEdge( edge, theShapeIDMap )) - Swap( tFace.myFirst[ iE ], tFace.myLast[ iE ] ); - } - // 2d corners - tFace.myCorner[ 0 ] = tFace.myC2d[ 0 ]->Value( tFace.myFirst[0] ).XY(); - tFace.myCorner[ 1 ] = tFace.myC2d[ 0 ]->Value( tFace.myLast[0] ).XY(); - tFace.myCorner[ 2 ] = tFace.myC2d[ 1 ]->Value( tFace.myLast[1] ).XY(); - tFace.myCorner[ 3 ] = tFace.myC2d[ 1 ]->Value( tFace.myFirst[1] ).XY(); - // sufrace - TopLoc_Location loc; - tFace.myS = BRep_Tool::Surface( face, loc ); - tFace.myTrsf = loc; - break; - } - default: break; - } - } // loop on shapes in theShapeIDMap +//======================================================================= +//function : TPoint() +//purpose : +//======================================================================= - return true; +SMESH_Pattern::TPoint::TPoint() +{ +#ifdef _DEBUG_ + myInitXYZ.SetCoord(0,0,0); + myInitUV.SetCoord(0.,0.); + myInitU = 0; + myXYZ.SetCoord(0,0,0); + myUV.SetCoord(0.,0.); + myU = 0; +#endif } //======================================================================= -//function : GetFaceEdgesIDs -//purpose : return edges IDs in the order u0, u1, 0v, 1v -// u0 means "|| u, v == 0" +//function : operator << +//purpose : //======================================================================= -void TBlock::GetFaceEdgesIDs (const int faceID, vector< int >& edgeVec ) +ostream & operator <<(ostream & OS, const SMESH_Pattern::TPoint& p) { - switch ( faceID ) { - case ID_Fxy0: - edgeVec[ 0 ] = ID_Ex00; - edgeVec[ 1 ] = ID_Ex10; - edgeVec[ 2 ] = ID_E0y0; - edgeVec[ 3 ] = ID_E1y0; - break; - case ID_Fxy1: - edgeVec[ 0 ] = ID_Ex01; - edgeVec[ 1 ] = ID_Ex11; - edgeVec[ 2 ] = ID_E0y1; - edgeVec[ 3 ] = ID_E1y1; - break; - case ID_Fx0z: - edgeVec[ 0 ] = ID_Ex00; - edgeVec[ 1 ] = ID_Ex01; - edgeVec[ 2 ] = ID_E00z; - edgeVec[ 3 ] = ID_E10z; - break; - case ID_Fx1z: - edgeVec[ 0 ] = ID_Ex10; - edgeVec[ 1 ] = ID_Ex11; - edgeVec[ 2 ] = ID_E01z; - edgeVec[ 3 ] = ID_E11z; - break; - case ID_F0yz: - edgeVec[ 0 ] = ID_E0y0; - edgeVec[ 1 ] = ID_E0y1; - edgeVec[ 2 ] = ID_E00z; - edgeVec[ 3 ] = ID_E01z; - break; - case ID_F1yz: - edgeVec[ 0 ] = ID_E1y0; - edgeVec[ 1 ] = ID_E1y1; - edgeVec[ 2 ] = ID_E10z; - edgeVec[ 3 ] = ID_E11z; - break; - default: - MESSAGE(" GetFaceEdgesIDs(), wrong face ID: " << faceID ); - } + gp_XYZ xyz = p.myInitXYZ; + OS << "\tinit( xyz( " << xyz.X() << " " << xyz.Y() << " " << xyz.Z() << " )"; + gp_XY xy = p.myInitUV; + OS << " uv( " << xy.X() << " " << xy.Y() << " )"; + double u = p.myInitU; + OS << " u( " << u << " )) " << &p << endl; + xyz = p.myXYZ.XYZ(); + OS << "\t ( xyz( " << xyz.X() << " " << xyz.Y() << " " << xyz.Z() << " )"; + xy = p.myUV; + OS << " uv( " << xy.X() << " " << xy.Y() << " )"; + u = p.myU; + OS << " u( " << u << " ))" << endl; + + return OS; }