Salome HOME
IPAL054122: Bad quality prismatic mesh
[modules/smesh.git] / src / SMESH / SMESH_Pattern.cxx
index 0a50cd9dd7b9d1a4a463f6559a5a58d23421ca83..db6e60e8ca629adbb489633dcffe0758783b0583 100644 (file)
@@ -256,7 +256,7 @@ int loadVE( const list< TopoDS_Edge > &          eList,
 //purpose  :
 //=======================================================================
 
-SMESH_Pattern::SMESH_Pattern ()
+SMESH_Pattern::SMESH_Pattern (): myToKeepNodes(false)
 {
 }
 
@@ -267,8 +267,6 @@ SMESH_Pattern::SMESH_Pattern ()
 
 bool SMESH_Pattern::Load (const char* theFileContents)
 {
-  MESSAGE("Load( file ) ");
-
   Kernel_Utils::Localizer loc;
   
   // file structure:
@@ -276,7 +274,7 @@ bool SMESH_Pattern::Load (const char* theFileContents)
   // ! This is a comment
   // NB_POINTS               ! 1 integer - the number of points in the pattern.
   //   X1 Y1 [Z1]            ! 2 or 3 reals - nodes coordinates within 2D or 3D domain:
-  //   X2 Y2 [Z2]            ! the pattern dimention is defined by the number of coordinates
+  //   X2 Y2 [Z2]            ! the pattern dimension is defined by the number of coordinates
   //   ...
   // [ ID1 ID2 ... IDn ]     ! Indices of key-points for a 2D pattern (only).
   // ! elements description goes after all
@@ -299,7 +297,7 @@ bool SMESH_Pattern::Load (const char* theFileContents)
 
   //   X1 Y1 [Z1]            ! 2 or 3 reals - nodes coordinates within 2D or 3D domain:
 
-  // read the first point coordinates to define pattern dimention
+  // read the first point coordinates to define pattern dimension
   int dim = readLine( fields, lineBeg, clearFields );
   if ( dim == 2 )
     myIs2D = true;
@@ -413,8 +411,6 @@ bool SMESH_Pattern::Load (const char* theFileContents)
 
 bool SMESH_Pattern::Save (ostream& theFile)
 {
-  MESSAGE(" ::Save(file) " );
-  
   Kernel_Utils::Localizer loc;
     
   if ( !IsLoaded() ) {
@@ -499,7 +495,7 @@ static gp_XY project (const SMDS_MeshNode* theNode,
     MESSAGE( "SMESH_Pattern: point projection FAILED");
     return gp_XY(0.,0.);
   }
-  double u, v, minVal = DBL_MAX;
+  double u =0, v =0, minVal = DBL_MAX;
   for ( int i = theProjectorPS.NbExt(); i > 0; i-- )
     if ( theProjectorPS.SquareDistance( i ) < minVal ) {
       minVal = theProjectorPS.SquareDistance( i );
@@ -565,11 +561,12 @@ static bool isMeshBoundToShape(SMESHDS_Mesh *     aMeshDS,
 bool SMESH_Pattern::Load (SMESH_Mesh*        theMesh,
                           const TopoDS_Face& theFace,
                           bool               theProject,
-                          TopoDS_Vertex      the1stVertex)
+                          TopoDS_Vertex      the1stVertex,
+                          bool               theKeepNodes)
 {
-  MESSAGE(" ::Load(face) " );
   Clear();
   myIs2D = true;
+  myToKeepNodes = theKeepNodes;
 
   SMESHDS_Mesh * aMeshDS = theMesh->GetMeshDS();
   SMESHDS_SubMesh * fSubMesh = aMeshDS->MeshElements( theFace );
@@ -605,8 +602,7 @@ bool SMESH_Pattern::Load (SMESH_Mesh*        theMesh,
 
   Extrema_GenExtPS projector;
   GeomAdaptor_Surface aSurface( BRep_Tool::Surface( face ));
-  if ( theProject || needProject )
-    projector.Initialize( aSurface, 20,20, 1e-5,1e-5 );
+  projector.Initialize( aSurface, 20,20, 1e-5,1e-5 );
 
   int iPoint = 0;
   TNodePointIDMap nodePointIDMap;
@@ -614,7 +610,6 @@ bool SMESH_Pattern::Load (SMESH_Mesh*        theMesh,
 
   if ( needProject )
   {
-    MESSAGE("Project the submesh");
     // ---------------------------------------------------------------
     // The case where the submesh is projected to theFace
     // ---------------------------------------------------------------
@@ -660,7 +655,7 @@ bool SMESH_Pattern::Load (SMESH_Mesh*        theMesh,
       const TopoDS_Vertex v = TopoDS::Vertex( vExp.Current() );
       gp_Pnt2d uv = BRep_Tool::Parameters( v, face );
       double minDist = DBL_MAX;
-      int index;
+      int index = 0;
       vector< TPoint >::const_iterator pVecIt = myPoints.begin();
       for ( iPoint = 0; pVecIt != myPoints.end(); pVecIt++, iPoint++ ) {
         double dist = uv.SquareDistance( (*pVecIt).myInitUV );
@@ -694,8 +689,28 @@ bool SMESH_Pattern::Load (SMESH_Mesh*        theMesh,
 
     myPoints.resize( nbNodes );
 
+    // care of INTERNAL VERTEXes
+    TopExp_Explorer vExp( face, TopAbs_VERTEX, TopAbs_EDGE );
+    for ( ; vExp.More(); vExp.Next() )
+    {
+      const SMDS_MeshNode* node =
+        SMESH_Algo::VertexNode( TopoDS::Vertex( vExp.Current()), aMeshDS );
+      if ( !node || node->NbInverseElements( SMDSAbs_Face ) == 0 )
+        continue;
+      myPoints.resize( ++nbNodes );
+      list< TPoint* > & fPoints = getShapePoints( face );
+      nodePointIDMap.insert( make_pair( node, iPoint ));
+      TPoint* p = &myPoints[ iPoint++ ];
+      fPoints.push_back( p );
+      gp_XY uv = helper.GetNodeUV( face, node );
+      p->myInitUV.SetCoord( uv.X(), uv.Y() );
+      p->myInitXYZ.SetCoord( p->myInitUV.X(), p->myInitUV.Y(), 0 );
+    }
+
     // Load U of points on edges
 
+    Bnd_Box2d edgesUVBox;
+
     list<int>::iterator nbEinW = myNbKeyPntInBoundary.begin();
     int iE = 0;
     vector< TopoDS_Edge > eVec;
@@ -766,6 +781,7 @@ bool SMESH_Pattern::Load (SMESH_Mesh*        theMesh,
           else
             keyPoint->myInitUV = C2d->Value( isForward ? f : l ).XY();
           keyPoint->myInitXYZ.SetCoord (keyPoint->myInitUV.X(), keyPoint->myInitUV.Y(), 0);
+          edgesUVBox.Add( gp_Pnt2d( keyPoint->myInitUV ));
         }
       }
       if ( !vPoint->empty() )
@@ -792,7 +808,7 @@ bool SMESH_Pattern::Load (SMESH_Mesh*        theMesh,
           double u = epos->GetUParameter();
           paramNodeMap.insert( make_pair( u, node ));
         }
-        if ((int) paramNodeMap.size() != eSubMesh->NbNodes() ) {
+        if ((int) paramNodeMap.size() != eSubMesh->NbNodes() - nbMeduimNodes ) {
           // wrong U on edge, project
           Extrema_ExtPC proj;
           BRepAdaptor_Curve aCurve( edge );
@@ -845,6 +861,7 @@ bool SMESH_Pattern::Load (SMESH_Mesh*        theMesh,
             p->myInitUV = C2d->Value( u ).XY();
           }
           p->myInitXYZ.SetCoord( p->myInitUV.X(), p->myInitUV.Y(), 0 );
+          edgesUVBox.Add( gp_Pnt2d( p->myInitUV ));
           unIt++; unRIt++;
           iPoint++;
         }
@@ -870,6 +887,7 @@ bool SMESH_Pattern::Load (SMESH_Mesh*        theMesh,
           else
             keyPoint->myInitUV = C2d->Value( isForward ? l : f ).XY();
           keyPoint->myInitXYZ.SetCoord( keyPoint->myInitUV.X(), keyPoint->myInitUV.Y(), 0 );
+          edgesUVBox.Add( gp_Pnt2d( keyPoint->myInitUV ));
         }
       }
       if ( !vPoint->empty() )
@@ -910,7 +928,7 @@ bool SMESH_Pattern::Load (SMESH_Mesh*        theMesh,
         nodePointIDMap.insert( make_pair( node, iPoint ));
         TPoint* p = &myPoints[ iPoint++ ];
         fPoints.push_back( p );
-        if ( theProject )
+        if ( theProject || edgesUVBox.IsOut( p->myInitUV ) )
           p->myInitUV = project( node, projector );
         else {
           const SMDS_FacePosition* pos =
@@ -966,6 +984,19 @@ bool SMESH_Pattern::Load (SMESH_Mesh*        theMesh,
     myIsBoundaryPointsFound = true;
   }
 
+  if ( myToKeepNodes )
+  {
+    myInNodes.resize( nodePointIDMap.size() + closeNodePointIDMap.size() );
+
+    TNodePointIDMap::iterator nIdIt = nodePointIDMap.begin();
+    for ( ; nIdIt != nodePointIDMap.end(); nIdIt++ )
+      myInNodes[ nIdIt->second ] = smdsNode( nIdIt->first );
+
+    nIdIt = closeNodePointIDMap.begin();
+    for ( ; nIdIt != closeNodePointIDMap.end(); nIdIt++ )
+      myInNodes[ nIdIt->second ] = smdsNode( nIdIt->first );
+  }
+
   // Assure that U range is proportional to V range
 
   Bnd_Box2d bndBox;
@@ -1089,10 +1120,10 @@ static bool intersectIsolines(const gp_XY& uv11, const gp_XY& uv12, const double
 
 //   resUV /= 2.;
 //     }
-  if ( isDeformed ) {
-    MESSAGE("intersectIsolines(), d1 = " << d1 << ", d2 = " << d2 << ", delta = " << delta <<
-            ", " << (loc1 - loc2).SquareModulus() << " > " << delta * delta);
-  }
+  // if ( isDeformed ) {
+  //   MESSAGE("intersectIsolines(), d1 = " << d1 << ", d2 = " << d2 << ", delta = " << delta <<
+  //           ", " << (loc1 - loc2).SquareModulus() << " > " << delta * delta);
+  // }
   return true;
 }
 
@@ -1185,7 +1216,7 @@ bool SMESH_Pattern::compUVByIsoIntersection (const list< list< TPoint* > >& theB
   }
   if ( !intersectIsolines( uv1[0], uv2[0], ratio[0],
                           uv1[1], uv2[1], ratio[1], theUV, theIsDeformed )) {
-    MESSAGE(" Cant intersect isolines for a point "<<theInitUV.X()<<", "<<theInitUV.Y());
+    MESSAGE(" Can't intersect isolines for a point "<<theInitUV.X()<<", "<<theInitUV.Y());
     return setErrorCode( ERR_APPLF_BAD_TOPOLOGY );
   }
 
@@ -1455,7 +1486,7 @@ static bool checkQuads (const TIsoNode* node,
         return false;
       }
       else {
-        //MESSAGE(" Cant improve UV, uv: "<<uv.X()<<" "<<uv.Y());
+        //MESSAGE(" Can't improve UV, uv: "<<uv.X()<<" "<<uv.Y());
       }
     }
     if ( !oldIsIn && nbOldFix ) {
@@ -1470,7 +1501,7 @@ static bool checkQuads (const TIsoNode* node,
         return false;
       }
       else {
-        //MESSAGE(" Cant fix UV, uv: "<<uv.X()<<" "<<uv.Y());
+        //MESSAGE(" Can't fix UV, uv: "<<uv.X()<<" "<<uv.Y());
       }
     }
     if ( newIsIn && oldIsIn )
@@ -1887,7 +1918,7 @@ bool SMESH_Pattern::
   list < TIsoNode* > internNodes;
   bool needIteration = true;
   if ( startNodes.empty() ) {
-    MESSAGE( " Starting UV by compUVByIsoIntersection()");
+    //MESSAGE( " Starting UV by compUVByIsoIntersection()");
     needIteration = false;
     map < double, TIsoLine >& isos = isoMap[ 0 ];
     map < double, TIsoLine >::iterator isoIt = isos.begin();
@@ -2107,7 +2138,7 @@ bool SMESH_Pattern::
 #endif
   } while ( maxMove > 1e-8 && nbIter++ < maxNbIter );
 
-  MESSAGE( "compUVByElasticIsolines(): Nb iterations " << nbIter << " dist: " << sqrt( maxMove ));
+  //MESSAGE( "compUVByElasticIsolines(): Nb iterations " << nbIter << " dist: " << sqrt( maxMove ));
 
   if ( nbIter >= maxNbIter && sqrt(maxMove) > minUvSize * 0.05 ) {
     MESSAGE( "compUVByElasticIsolines() failed: "<<sqrt(maxMove)<<">"<<minUvSize * 0.05);
@@ -2288,7 +2319,7 @@ bool SMESH_Pattern::sortSameSizeWires (TListOfEdgesList &                theWire
       list< TPoint* > & ePoints = getShapePoints( eID++ );
       TPoint* p = ePoints.front();
       if ( !compUVByIsoIntersection( theEdgesPointsList, p->myInitUV, p->myUV, aBool )) {
-        MESSAGE("cant sortSameSizeWires()");
+        MESSAGE("can't sortSameSizeWires()");
         return false;
       }
       gcVec[iW] += p->myUV;
@@ -2339,7 +2370,7 @@ bool SMESH_Pattern::sortSameSizeWires (TListOfEdgesList &                theWire
 //   " \t vertex: " << vGcVec[iW].X() << " " << vGcVec[iW].Y() << endl;
     double minDist = DBL_MAX;
     gp_XY & wGc = vGcVec[ iW ];
-    int bIndex;
+    int bIndex = 0;
     for ( int iB = 0; iB < nbWires; iB++ ) {
       if ( bndFound[ iB ] ) continue;
       double dist = ( wGc - gcVec[ iB ] ).SquareModulus();
@@ -2392,7 +2423,6 @@ bool SMESH_Pattern::Apply (const TopoDS_Face&   theFace,
                            const TopoDS_Vertex& theVertexOnKeyPoint1,
                            const bool           theReverse)
 {
-  MESSAGE(" ::Apply(face) " );
   TopoDS_Face face  = theReverse ? TopoDS::Face( theFace.Reversed() ) : theFace;
   if ( !setShapeToMesh( face ))
     return false;
@@ -2457,7 +2487,7 @@ bool SMESH_Pattern::Apply (const TopoDS_Face&   theFace,
   if ( nbWires > 1 )
   {
     // compute UV of inner edge-points using the method for in-face points
-    // and devide eList into a list of separate wires
+    // and divide eList into a list of separate wires
     bool aBool;
     list< list< TopoDS_Edge > > wireList;
     list<TopoDS_Edge>::iterator eIt = elIt;
@@ -2474,7 +2504,7 @@ bool SMESH_Pattern::Apply (const TopoDS_Face&   theFace,
         for (  pIt++; pIt != ePoints.end(); pIt++ ) {
           TPoint* p = (*pIt);
           if ( !compUVByIsoIntersection( edgesPointsList, p->myInitUV, p->myUV, aBool )) {
-            MESSAGE("cant Apply(face)");
+            MESSAGE("can't Apply(face)");
             return false;
           }
           // keep the computed UV to compare against by setFirstEdge()
@@ -2593,7 +2623,7 @@ bool SMESH_Pattern::Apply (const TopoDS_Face&   theFace,
   for ( pIt = fPoints.begin(); !isDeformed && pIt != fPoints.end(); pIt++ )
     if ( !compUVByIsoIntersection( edgesPointsList, (*pIt)->myInitUV,
                                   (*pIt)->myUV, isDeformed )) {
-      MESSAGE("cant Apply(face)");
+      MESSAGE("can't Apply(face)");
       return false;
     }
   // try to use a complex algo if it is a difficult case
@@ -2602,7 +2632,7 @@ bool SMESH_Pattern::Apply (const TopoDS_Face&   theFace,
     for ( ; pIt != fPoints.end(); pIt++ ) // continue with the simple algo
       if ( !compUVByIsoIntersection( edgesPointsList, (*pIt)->myInitUV,
                                     (*pIt)->myUV, isDeformed )) {
-        MESSAGE("cant Apply(face)");
+        MESSAGE("can't Apply(face)");
         return false;
       }
   }
@@ -2742,7 +2772,7 @@ bool SMESH_Pattern::Apply (const SMDS_MeshFace* theFace,
   for ( pIt = fPoints.begin(); !isDeformed && pIt != fPoints.end(); pIt++ )
     if ( !compUVByIsoIntersection( edgesPointsList, (*pIt)->myInitUV,
                                   (*pIt)->myUV, isDeformed )) {
-      MESSAGE("cant Apply(face)");
+      MESSAGE("can't Apply(face)");
       return false;
     }
   // try to use a complex algo if it is a difficult case
@@ -2751,7 +2781,7 @@ bool SMESH_Pattern::Apply (const SMDS_MeshFace* theFace,
     for ( ; pIt != fPoints.end(); pIt++ ) // continue with the simple algo
       if ( !compUVByIsoIntersection( edgesPointsList, (*pIt)->myInitUV,
                                     (*pIt)->myUV, isDeformed )) {
-        MESSAGE("cant Apply(face)");
+        MESSAGE("can't Apply(face)");
         return false;
       }
   }
@@ -2895,7 +2925,7 @@ bool SMESH_Pattern::Apply (SMESH_Mesh*          theMesh,
   for ( pIt = fPoints.begin(); !isDeformed && pIt != fPoints.end(); pIt++ )
     if ( !compUVByIsoIntersection( edgesPointsList, (*pIt)->myInitUV,
                                   (*pIt)->myUV, isDeformed )) {
-      MESSAGE("cant Apply(face)");
+      MESSAGE("can't Apply(face)");
       return false;
     }
   // try to use a complex algo if it is a difficult case
@@ -2904,7 +2934,7 @@ bool SMESH_Pattern::Apply (SMESH_Mesh*          theMesh,
     for ( ; pIt != fPoints.end(); pIt++ ) // continue with the simple algo
       if ( !compUVByIsoIntersection( edgesPointsList, (*pIt)->myInitUV,
                                     (*pIt)->myUV, isDeformed )) {
-        MESSAGE("cant Apply(face)");
+        MESSAGE("can't Apply(face)");
         return false;
       }
   }
@@ -3075,8 +3105,6 @@ bool SMESH_Pattern::Apply (std::set<const SMDS_MeshVolume*> & theVolumes,
                            const int                          theNode000Index,
                            const int                          theNode001Index)
 {
-  MESSAGE(" ::Apply(set<MeshVolumes>) " );
-
   if ( !IsLoaded() ) {
     MESSAGE( "Pattern not loaded" );
     return setErrorCode( ERR_APPL_NOT_LOADED );
@@ -3181,11 +3209,12 @@ bool SMESH_Pattern::Apply (std::set<const SMDS_MeshVolume*> & theVolumes,
 //=======================================================================
 
 bool SMESH_Pattern::Load (SMESH_Mesh*         theMesh,
-                          const TopoDS_Shell& theBlock)
+                          const TopoDS_Shell& theBlock,
+                          bool                theKeepNodes)
 {
-  MESSAGE(" ::Load(volume) " );
   Clear();
   myIs2D = false;
+  myToKeepNodes = theKeepNodes;
   SMESHDS_SubMesh * aSubMesh;
 
   const bool isQuadMesh = theMesh->NbVolumes( ORDER_QUADRATIC );
@@ -3219,7 +3248,7 @@ bool SMESH_Pattern::Load (SMESH_Mesh*         theMesh,
     SMDS_NodeIteratorPtr nIt = aSubMesh->GetNodes();
     if ( !nIt->more() ) continue;
 
-      // store a node and a point
+    // store a node and a point
     while ( nIt->more() ) {
       const SMDS_MeshNode* node = smdsNode( nIt->next() );
       if ( isQuadMesh && SMESH_MeshEditor::IsMedium( node, SMDSAbs_Volume ))
@@ -3297,6 +3326,14 @@ bool SMESH_Pattern::Load (SMESH_Mesh*         theMesh,
 
   myIsBoundaryPointsFound = true;
 
+  if ( myToKeepNodes )
+  {
+    myInNodes.resize( nodePointIDMap.size() );
+    TNodePointIDMap::iterator nIdIt = nodePointIDMap.begin();
+    for ( ; nIdIt != nodePointIDMap.end(); nIdIt++ )
+      myInNodes[ nIdIt->second ] = smdsNode( nIdIt->first );
+  }
+
   return setErrorCode( ERR_OK );
 }
 
@@ -3338,8 +3375,6 @@ bool SMESH_Pattern::Apply (const TopoDS_Shell&  theBlock,
                            const TopoDS_Vertex& theVertex000,
                            const TopoDS_Vertex& theVertex001)
 {
-  MESSAGE(" ::Apply(volume) " );
-
   if (!findBoundaryPoints()     || // bind ID to points
       !setShapeToMesh( theBlock )) // check theBlock is a suitable shape
     return false;
@@ -3399,8 +3434,6 @@ bool SMESH_Pattern::Apply (const SMDS_MeshVolume* theVolume,
                            const int              theNode000Index,
                            const int              theNode001Index)
 {
-  //MESSAGE(" ::Apply(MeshVolume) " );
-
   if (!findBoundaryPoints()) // bind ID to points
     return false;
 
@@ -3991,7 +4024,6 @@ bool SMESH_Pattern::MakeMesh(SMESH_Mesh* theMesh,
                              const bool  toCreatePolygons,
                              const bool  toCreatePolyedrs)
 {
-  MESSAGE(" ::MakeMesh() " );
   if ( !myIsComputed )
     return setErrorCode( ERR_MAKEM_NOT_COMPUTED );
 
@@ -4144,6 +4176,9 @@ bool SMESH_Pattern::MakeMesh(SMESH_Mesh* theMesh,
 
   aMeshDS->compactMesh();
 
+  if ( myToKeepNodes )
+    myOutNodes.swap( nodesVector );
+
 //   const map<int,SMESHDS_SubMesh*>& sm = aMeshDS->SubMeshes();
 //   map<int,SMESHDS_SubMesh*>::const_iterator i_sm = sm.begin();
 //   for ( ; i_sm != sm.end(); i_sm++ )
@@ -4524,8 +4559,6 @@ bool SMESH_Pattern::findBoundaryPoints()
 {
   if ( myIsBoundaryPointsFound ) return true;
 
-  MESSAGE(" findBoundaryPoints() ");
-
   myNbKeyPntInBoundary.clear();
 
   if ( myIs2D )
@@ -4767,7 +4800,7 @@ bool SMESH_Pattern::setShapeToMesh(const TopoDS_Shape& theShape)
   TopAbs_ShapeEnum aType = theShape.ShapeType();
   bool dimOk = ( myIs2D ? aType == TopAbs_FACE : aType == TopAbs_SHELL );
   if ( !dimOk ) {
-    MESSAGE( "Pattern dimention mismatch" );
+    MESSAGE( "Pattern dimension mismatch" );
     return setErrorCode( ERR_APPL_BAD_DIMENTION );
   }