src/SMESH/SMESH_MesherHelper.cxx

   1 // Copyright (C) 2005  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   2 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   3 //
   4 // This library is free software; you can redistribute it and/or
   5 // modify it under the terms of the GNU Lesser General Public
   6 // License as published by the Free Software Foundation; either
   7 // version 2.1 of the License.
   8 //
   9 // This library is distributed in the hope that it will be useful
  10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  12 // Lesser General Public License for more details.
  13 //
  14 // You should have received a copy of the GNU Lesser General Public
  15 // License along with this library; if not, write to the Free Software
  16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  17 //
  18 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  19 //
  20 // File:      SMESH_MesherHelper.cxx
  21 // Created:   15.02.06 15:22:41
  22 // Author:    Sergey KUUL
  23 // Copyright: Open CASCADE 2006
  24
  25
  26 #include "SMESH_MesherHelper.hxx"
  27
  28 #include "SMDS_FacePosition.hxx"
  29 #include "SMDS_EdgePosition.hxx"
  30 #include "SMESH_MeshEditor.hxx"
  31
  32 #include <BRepAdaptor_Surface.hxx>
  33 #include <BRepTools.hxx>
  34 #include <BRep_Tool.hxx>
  35 #include <Geom2d_Curve.hxx>
  36 #include <Geom_Curve.hxx>
  37 #include <Geom_Surface.hxx>
  38 #include <TopExp_Explorer.hxx>
  39 #include <TopTools_MapOfShape.hxx>
  40 #include <gp_Pnt2d.hxx>
  41 #include <ShapeAnalysis.hxx>
  42 #include <TopTools_ListIteratorOfListOfShape.hxx>
  43
  44 #include <Standard_Failure.hxx>
  45 #include <Standard_ErrorHandler.hxx>
  46
  47 #include <utilities.h>
  48
  49 #define RETURN_BAD_RESULT(msg) { MESSAGE(msg); return false; }
  50
  51 //================================================================================
  52 /*!
  53  * \brief Constructor
  54  */
  55 //================================================================================
  56
  57 SMESH_MesherHelper::SMESH_MesherHelper(SMESH_Mesh& theMesh)
  58   : myMesh(&theMesh), myShapeID(-1), myCreateQuadratic(false)
  59 {
  60   mySetElemOnShape = ( ! myMesh->HasShapeToMesh() );
  61 }
  62
  63 //=======================================================================
  64 //function : CheckShape
  65 //purpose  :
  66 //=======================================================================
  67
  68 bool SMESH_MesherHelper::IsQuadraticSubMesh(const TopoDS_Shape& aSh)
  69 {
  70   SMESHDS_Mesh* meshDS = GetMeshDS();
  71   // we can create quadratic elements only if all elements
  72   // created on subshapes of given shape are quadratic
  73   // also we have to fill myNLinkNodeMap
  74   myCreateQuadratic = true;
  75   mySeamShapeIds.clear();
  76   TopAbs_ShapeEnum subType( aSh.ShapeType()==TopAbs_FACE ? TopAbs_EDGE : TopAbs_FACE );
  77   SMDSAbs_ElementType elemType( subType==TopAbs_FACE ? SMDSAbs_Face : SMDSAbs_Edge );
  78
  79   TopExp_Explorer exp( aSh, subType );
  80   for (; exp.More() && myCreateQuadratic; exp.Next()) {
  81     if ( SMESHDS_SubMesh * subMesh = meshDS->MeshElements( exp.Current() )) {
  82       if ( SMDS_ElemIteratorPtr it = subMesh->GetElements() ) {
  83         while(it->more()) {
  84           const SMDS_MeshElement* e = it->next();
  85           if ( e->GetType() != elemType || !e->IsQuadratic() ) {
  86             myCreateQuadratic = false;
  87             break;
  88           }
  89           else {
  90             // fill NLinkNodeMap
  91             switch ( e->NbNodes() ) {
  92             case 3:
  93               AddNLinkNode(e->GetNode(0),e->GetNode(1),e->GetNode(2)); break;
  94             case 6:
  95               AddNLinkNode(e->GetNode(0),e->GetNode(1),e->GetNode(3));
  96               AddNLinkNode(e->GetNode(1),e->GetNode(2),e->GetNode(4));
  97               AddNLinkNode(e->GetNode(2),e->GetNode(0),e->GetNode(5)); break;
  98             case 8:
  99               AddNLinkNode(e->GetNode(0),e->GetNode(1),e->GetNode(4));
 100               AddNLinkNode(e->GetNode(1),e->GetNode(2),e->GetNode(5));
 101               AddNLinkNode(e->GetNode(2),e->GetNode(3),e->GetNode(6));
 102               AddNLinkNode(e->GetNode(3),e->GetNode(0),e->GetNode(7));
 103               break;
 104             default:
 105               myCreateQuadratic = false;
 106               break;
 107             }
 108           }
 109         }
 110       }
 111     }
 112   }
 113
 114   if(!myCreateQuadratic) {
 115     myNLinkNodeMap.clear();
 116   }
 117   SetSubShape( aSh );
 118
 119   return myCreateQuadratic;
 120 }
 121
 122 //================================================================================
 123 /*!
 124  * \brief Set geomerty to make elements on
 125   * \param aSh - geomertic shape
 126  */
 127 //================================================================================
 128
 129 void SMESH_MesherHelper::SetSubShape(const int aShID)
 130 {
 131   if ( aShID == myShapeID )
 132     return;
 133   if ( aShID > 1 )
 134     SetSubShape( GetMeshDS()->IndexToShape( aShID ));
 135   else
 136     SetSubShape( TopoDS_Shape() );
 137 }
 138
 139 //================================================================================
 140 /*!
 141  * \brief Set geomerty to make elements on
 142   * \param aSh - geomertic shape
 143  */
 144 //================================================================================
 145
 146 void SMESH_MesherHelper::SetSubShape(const TopoDS_Shape& aSh)
 147 {
 148   if ( myShape.IsSame( aSh ))
 149     return;
 150
 151   myShape = aSh;
 152   mySeamShapeIds.clear();
 153
 154   if ( myShape.IsNull() ) {
 155     myShapeID  = -1;
 156     return;
 157   }
 158   SMESHDS_Mesh* meshDS = GetMeshDS();
 159   myShapeID = meshDS->ShapeToIndex(aSh);
 160
 161   // treatment of periodic faces
 162   for ( TopExp_Explorer eF( aSh, TopAbs_FACE ); eF.More(); eF.Next() )
 163   {
 164     const TopoDS_Face& face = TopoDS::Face( eF.Current() );
 165     BRepAdaptor_Surface surface( face );
 166     if ( surface.IsUPeriodic() || surface.IsVPeriodic() )
 167     {
 168       // look for a seam edge
 169       for ( TopExp_Explorer exp( face, TopAbs_EDGE ); exp.More(); exp.Next()) {
 170         const TopoDS_Edge& edge = TopoDS::Edge( exp.Current() );
 171         if ( BRep_Tool::IsClosed( edge, face )) {
 172           // initialize myPar1, myPar2 and myParIndex
 173           if ( mySeamShapeIds.empty() ) {
 174             gp_Pnt2d uv1, uv2;
 175             BRep_Tool::UVPoints( edge, face, uv1, uv2 );
 176             if ( Abs( uv1.Coord(1) - uv2.Coord(1) ) < Abs( uv1.Coord(2) - uv2.Coord(2) ))
 177             {
 178               myParIndex = 1; // U periodic
 179               myPar1 = surface.FirstUParameter();
 180               myPar2 = surface.LastUParameter();
 181             }
 182             else {
 183               myParIndex = 2;  // V periodic
 184               myPar1 = surface.FirstVParameter();
 185               myPar2 = surface.LastVParameter();
 186             }
 187           }
 188           // store shapes indices
 189           mySeamShapeIds.insert( meshDS->ShapeToIndex( exp.Current() ));
 190           for ( TopExp_Explorer v( exp.Current(), TopAbs_VERTEX ); v.More(); v.Next() )
 191             mySeamShapeIds.insert( meshDS->ShapeToIndex( v.Current() ));
 192         }
 193       }
 194     }
 195   }
 196 }
 197
 198 //================================================================================
 199   /*!
 200    * \brief Check if inFaceNode argument is necessary for call GetNodeUV(F,..)
 201     * \param F - the face
 202     * \retval bool - return true if the face is periodic
 203    */
 204 //================================================================================
 205
 206 bool SMESH_MesherHelper::GetNodeUVneedInFaceNode(const TopoDS_Face& F) const
 207 {
 208   if ( F.IsNull() ) return !mySeamShapeIds.empty();
 209
 210   if ( !F.IsNull() && !myShape.IsNull() && myShape.IsSame( F ))
 211     return !mySeamShapeIds.empty();
 212
 213   Handle(Geom_Surface) aSurface = BRep_Tool::Surface( F );
 214   if ( !aSurface.IsNull() )
 215     return ( aSurface->IsUPeriodic() || aSurface->IsVPeriodic() );
 216
 217   return false;
 218 }
 219
 220 //=======================================================================
 221 //function : IsMedium
 222 //purpose  :
 223 //=======================================================================
 224
 225 bool SMESH_MesherHelper::IsMedium(const SMDS_MeshNode*      node,
 226                                  const SMDSAbs_ElementType typeToCheck)
 227 {
 228   return SMESH_MeshEditor::IsMedium( node, typeToCheck );
 229 }
 230
 231 //=======================================================================
 232 //function : AddNLinkNode
 233 //purpose  :
 234 //=======================================================================
 235 /*!
 236  * Auxilary function for filling myNLinkNodeMap
 237  */
 238 void SMESH_MesherHelper::AddNLinkNode(const SMDS_MeshNode* n1,
 239                                      const SMDS_MeshNode* n2,
 240                                      const SMDS_MeshNode* n12)
 241 {
 242   NLink link( n1, n2 );
 243   if ( n1 > n2 ) link = NLink( n2, n1 );
 244   // add new record to map
 245   myNLinkNodeMap.insert( make_pair(link,n12));
 246 }
 247
 248 //=======================================================================
 249 /*!
 250  * \brief Select UV on either of 2 pcurves of a seam edge, closest to the given UV
 251  * \param uv1 - UV on the seam
 252  * \param uv2 - UV within a face
 253  * \retval gp_Pnt2d - selected UV
 254  */
 255 //=======================================================================
 256
 257 gp_Pnt2d SMESH_MesherHelper::GetUVOnSeam( const gp_Pnt2d& uv1, const gp_Pnt2d& uv2 ) const
 258 {
 259   double p1 = uv1.Coord( myParIndex );
 260   double p2 = uv2.Coord( myParIndex );
 261   double p3 = ( Abs( p1 - myPar1 ) < Abs( p1 - myPar2 )) ? myPar2 : myPar1;
 262   if ( Abs( p2 - p1 ) > Abs( p2 - p3 ))
 263     p1 = p3;
 264   gp_Pnt2d result = uv1;
 265   result.SetCoord( myParIndex, p1 );
 266   return result;
 267 }
 268
 269 //=======================================================================
 270 /*!
 271  * \brief Return node UV on face
 272  * \param F - the face
 273  * \param n - the node
 274  * \param n2 - a node of element being created located inside a face
 275  * \retval gp_XY - resulting UV
 276  *
 277  * Auxilary function called form GetMediumNode()
 278  */
 279 //=======================================================================
 280
 281 gp_XY SMESH_MesherHelper::GetNodeUV(const TopoDS_Face&   F,
 282                                     const SMDS_MeshNode* n,
 283                                     const SMDS_MeshNode* n2) const
 284 {
 285   gp_Pnt2d uv( 1e100, 1e100 );
 286   const SMDS_PositionPtr Pos = n->GetPosition();
 287   if(Pos->GetTypeOfPosition()==SMDS_TOP_FACE)
 288   {
 289     // node has position on face
 290     const SMDS_FacePosition* fpos =
 291       static_cast<const SMDS_FacePosition*>(n->GetPosition().get());
 292     uv = gp_Pnt2d(fpos->GetUParameter(),fpos->GetVParameter());
 293   }
 294   else if(Pos->GetTypeOfPosition()==SMDS_TOP_EDGE)
 295   {
 296     // node has position on edge => it is needed to find
 297     // corresponding edge from face, get pcurve for this
 298     // edge and recieve value from this pcurve
 299     const SMDS_EdgePosition* epos =
 300       static_cast<const SMDS_EdgePosition*>(n->GetPosition().get());
 301     SMESHDS_Mesh* meshDS = GetMeshDS();
 302     int edgeID = Pos->GetShapeId();
 303     TopoDS_Edge E = TopoDS::Edge(meshDS->IndexToShape(edgeID));
 304     double f, l;
 305     TopLoc_Location loc;
 306     Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(E, F, f, l);
 307     uv = C2d->Value( epos->GetUParameter() );
 308     // for a node on a seam edge select one of UVs on 2 pcurves
 309     if ( n2 && mySeamShapeIds.find( edgeID ) != mySeamShapeIds.end() )
 310       uv = GetUVOnSeam( uv, GetNodeUV( F, n2, 0 ));
 311   }
 312   else if(Pos->GetTypeOfPosition()==SMDS_TOP_VERTEX)
 313   {
 314     if ( int vertexID = n->GetPosition()->GetShapeId() ) {
 315       bool ok = true;
 316       const TopoDS_Vertex& V = TopoDS::Vertex(GetMeshDS()->IndexToShape(vertexID));
 317       try {
 318         uv = BRep_Tool::Parameters( V, F );
 319       }
 320       catch (Standard_Failure& exc) {
 321         ok = false;
 322       }
 323       if ( !ok ) {
 324         for ( TopExp_Explorer vert(F,TopAbs_VERTEX); !ok && vert.More(); vert.Next() )
 325           ok = ( V == vert.Current() );
 326         if ( !ok ) {
 327 #ifdef _DEBUG_
 328           cout << "SMESH_MesherHelper::GetNodeUV(); Vertex " << vertexID
 329                << " not in face " << GetMeshDS()->ShapeToIndex( F ) << endl;
 330 #endif
 331           // get UV of a vertex closest to the node
 332           double dist = 1e100;
 333           gp_Pnt pn ( n->X(),n->Y(),n->Z() );
 334           for ( TopExp_Explorer vert(F,TopAbs_VERTEX); !ok && vert.More(); vert.Next() ) {
 335             TopoDS_Vertex curV = TopoDS::Vertex( vert.Current() );
 336             gp_Pnt p = BRep_Tool::Pnt( curV );
 337             double curDist = p.SquareDistance( pn );
 338             if ( curDist < dist ) {
 339               dist = curDist;
 340               uv = BRep_Tool::Parameters( curV, F );
 341               if ( dist < DBL_MIN ) break;
 342             }
 343           }
 344         }
 345         else {
 346           TopTools_ListIteratorOfListOfShape it( myMesh->GetAncestors( V ));
 347           for ( ; it.More(); it.Next() ) {
 348             if ( it.Value().ShapeType() == TopAbs_EDGE ) {
 349               const TopoDS_Edge & edge = TopoDS::Edge( it.Value() );
 350               double f,l;
 351               Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(edge, F, f, l);
 352               if ( !C2d.IsNull() ) {
 353                 double u = ( V == TopExp::FirstVertex( edge ) ) ?  f : l;
 354                 uv = C2d->Value( u );
 355                 break;
 356               }
 357             }
 358           }
 359         }
 360       }
 361       if ( n2 && mySeamShapeIds.find( vertexID ) != mySeamShapeIds.end() )
 362         uv = GetUVOnSeam( uv, GetNodeUV( F, n2, 0 ));
 363     }
 364   }
 365   return uv.XY();
 366 }
 367
 368 //=======================================================================
 369 /*!
 370  * \brief Return node U on edge
 371  * \param E - the Edge
 372  * \param n - the node
 373  * \retval double - resulting U
 374  *
 375  * Auxilary function called form GetMediumNode()
 376  */
 377 //=======================================================================
 378
 379 double SMESH_MesherHelper::GetNodeU(const TopoDS_Edge&   E,
 380                                     const SMDS_MeshNode* n)
 381 {
 382   double param = 0;
 383   const SMDS_PositionPtr Pos = n->GetPosition();
 384   if(Pos->GetTypeOfPosition()==SMDS_TOP_EDGE) {
 385     const SMDS_EdgePosition* epos =
 386       static_cast<const SMDS_EdgePosition*>(n->GetPosition().get());
 387     param =  epos->GetUParameter();
 388   }
 389   else if(Pos->GetTypeOfPosition()==SMDS_TOP_VERTEX) {
 390     SMESHDS_Mesh * meshDS = GetMeshDS();
 391     int vertexID = n->GetPosition()->GetShapeId();
 392     const TopoDS_Vertex& V = TopoDS::Vertex(meshDS->IndexToShape(vertexID));
 393     param =  BRep_Tool::Parameter( V, E );
 394   }
 395   return param;
 396 }
 397
 398 //=======================================================================
 399 //function : GetMediumNode
 400 //purpose  :
 401 //=======================================================================
 402 /*!
 403  * Special function for search or creation medium node
 404  */
 405 const SMDS_MeshNode* SMESH_MesherHelper::GetMediumNode(const SMDS_MeshNode* n1,
 406                                                        const SMDS_MeshNode* n2,
 407                                                        bool force3d)
 408 {
 409   TopAbs_ShapeEnum shapeType = myShape.IsNull() ? TopAbs_SHAPE : myShape.ShapeType();
 410
 411   NLink link(( n1 < n2 ? n1 : n2 ), ( n1 < n2 ? n2 : n1 ));
 412   ItNLinkNode itLN = myNLinkNodeMap.find( link );
 413   if ( itLN != myNLinkNodeMap.end() ) {
 414     return (*itLN).second;
 415   }
 416   else {
 417     // create medium node
 418     SMDS_MeshNode* n12;
 419     SMESHDS_Mesh* meshDS = GetMeshDS();
 420     int faceID = -1, edgeID = -1;
 421     const SMDS_PositionPtr Pos1 = n1->GetPosition();
 422     const SMDS_PositionPtr Pos2 = n2->GetPosition();
 423
 424     if( myShape.IsNull() )
 425     {
 426       if( Pos1->GetTypeOfPosition()==SMDS_TOP_FACE ) {
 427         faceID = Pos1->GetShapeId();
 428       }
 429       else if( Pos2->GetTypeOfPosition()==SMDS_TOP_FACE ) {
 430         faceID = Pos2->GetShapeId();
 431       }
 432
 433       if( Pos1->GetTypeOfPosition()==SMDS_TOP_EDGE ) {
 434         edgeID = Pos1->GetShapeId();
 435       }
 436       if( Pos2->GetTypeOfPosition()==SMDS_TOP_EDGE ) {
 437         edgeID = Pos2->GetShapeId();
 438       }
 439     }
 440
 441     if(!force3d) {
 442       // we try to create medium node using UV parameters of
 443       // nodes, else - medium between corresponding 3d points
 444       if(faceID>-1 || shapeType == TopAbs_FACE) {
 445         // obtaining a face and 2d points for nodes
 446         TopoDS_Face F;
 447         if( myShape.IsNull() )
 448           F = TopoDS::Face(meshDS->IndexToShape(faceID));
 449         else {
 450           F = TopoDS::Face(myShape);
 451           faceID = myShapeID;
 452         }
 453
 454         gp_XY p1 = GetNodeUV(F,n1,n2);
 455         gp_XY p2 = GetNodeUV(F,n2,n1);
 456
 457         //checking if surface is periodic
 458         Handle(Geom_Surface) S = BRep_Tool::Surface(F);
 459         Standard_Real UF,UL,VF,VL;
 460         S->Bounds(UF,UL,VF,VL);
 461
 462         Standard_Real u,v;
 463         Standard_Boolean isUPeriodic = S->IsUPeriodic();
 464         if(isUPeriodic) {
 465           Standard_Real UPeriod = S->UPeriod();
 466           Standard_Real p2x = p2.X()+ShapeAnalysis::AdjustByPeriod(p2.X(),p1.X(),UPeriod);
 467           Standard_Real pmid = (p1.X()+p2x)/2.;
 468           u = pmid+ShapeAnalysis::AdjustToPeriod(pmid,UF,UL);
 469         }
 470         else
 471           u= (p1.X()+p2.X())/2.;
 472
 473         Standard_Boolean isVPeriodic = S->IsVPeriodic();
 474         if(isVPeriodic) {
 475           Standard_Real VPeriod = S->VPeriod();
 476           Standard_Real p2y = p2.Y()+ShapeAnalysis::AdjustByPeriod(p2.Y(),p1.Y(),VPeriod);
 477           Standard_Real pmid = (p1.Y()+p2y)/2.;
 478           v = pmid+ShapeAnalysis::AdjustToPeriod(pmid,VF,VL);
 479         }
 480         else
 481           v = (p1.Y()+p2.Y())/2.;
 482
 483         gp_Pnt P = S->Value(u, v);
 484         n12 = meshDS->AddNode(P.X(), P.Y(), P.Z());
 485         meshDS->SetNodeOnFace(n12, faceID, u, v);
 486         myNLinkNodeMap.insert(NLinkNodeMap::value_type(link,n12));
 487         return n12;
 488       }
 489       if (edgeID>-1 || shapeType == TopAbs_EDGE) {
 490
 491         TopoDS_Edge E;
 492         if( myShape.IsNull() )
 493           E = TopoDS::Edge(meshDS->IndexToShape(edgeID));
 494         else {
 495           E = TopoDS::Edge(myShape);
 496           edgeID = myShapeID;
 497         }
 498
 499         double p1 = GetNodeU(E,n1);
 500         double p2 = GetNodeU(E,n2);
 501
 502         double f,l;
 503         Handle(Geom_Curve) C = BRep_Tool::Curve(E, f, l);
 504         if(!C.IsNull()) {
 505
 506           Standard_Boolean isPeriodic = C->IsPeriodic();
 507           double u;
 508           if(isPeriodic) {
 509             Standard_Real Period = C->Period();
 510             Standard_Real p = p2+ShapeAnalysis::AdjustByPeriod(p2,p1,Period);
 511             Standard_Real pmid = (p1+p)/2.;
 512             u = pmid+ShapeAnalysis::AdjustToPeriod(pmid,C->FirstParameter(),C->LastParameter());
 513           }
 514           else
 515             u = (p1+p2)/2.;
 516
 517           gp_Pnt P = C->Value( u );
 518           n12 = meshDS->AddNode(P.X(), P.Y(), P.Z());
 519           meshDS->SetNodeOnEdge(n12, edgeID, u);
 520           myNLinkNodeMap.insert(NLinkNodeMap::value_type(link,n12));
 521           return n12;
 522         }
 523       }
 524     }
 525     // 3d variant
 526     double x = ( n1->X() + n2->X() )/2.;
 527     double y = ( n1->Y() + n2->Y() )/2.;
 528     double z = ( n1->Z() + n2->Z() )/2.;
 529     n12 = meshDS->AddNode(x,y,z);
 530     if(edgeID>-1)
 531         meshDS->SetNodeOnEdge(n12, edgeID);
 532     else if(faceID>-1)
 533         meshDS->SetNodeOnFace(n12, faceID);
 534     else
 535       meshDS->SetNodeInVolume(n12, myShapeID);
 536     myNLinkNodeMap.insert(NLinkNodeMap::value_type(link,n12));
 537     return n12;
 538   }
 539 }
 540
 541 //=======================================================================
 542 /*!
 543  * Creates a node
 544  */
 545 //=======================================================================
 546
 547 SMDS_MeshNode* SMESH_MesherHelper::AddNode(double x, double y, double z, int ID)
 548 {
 549   SMESHDS_Mesh * meshDS = GetMeshDS();
 550   SMDS_MeshNode* node = 0;
 551   if ( ID )
 552     node = meshDS->AddNodeWithID( x, y, z, ID );
 553   else
 554     node = meshDS->AddNode( x, y, z );
 555   if ( mySetElemOnShape && myShapeID > 0 ) {
 556     switch ( myShape.ShapeType() ) {
 557     case TopAbs_SOLID:  meshDS->SetNodeInVolume( node, myShapeID); break;
 558     case TopAbs_SHELL:  meshDS->SetNodeInVolume( node, myShapeID); break;
 559     case TopAbs_FACE:   meshDS->SetNodeOnFace(   node, myShapeID); break;
 560     case TopAbs_EDGE:   meshDS->SetNodeOnEdge(   node, myShapeID); break;
 561     case TopAbs_VERTEX: meshDS->SetNodeOnVertex( node, myShapeID); break;
 562     default: ;
 563     }
 564   }
 565   return node;
 566 }
 567
 568 //=======================================================================
 569 /*!
 570  * Creates quadratic or linear edge
 571  */
 572 //=======================================================================
 573
 574 SMDS_MeshEdge* SMESH_MesherHelper::AddEdge(const SMDS_MeshNode* n1,
 575                                                 const SMDS_MeshNode* n2,
 576                                                 const int id,
 577                                                 const bool force3d)
 578 {
 579   SMESHDS_Mesh * meshDS = GetMeshDS();
 580
 581   SMDS_MeshEdge* edge = 0;
 582   if (myCreateQuadratic) {
 583     const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
 584     if(id)
 585       edge = meshDS->AddEdgeWithID(n1, n2, n12, id);
 586     else
 587       edge = meshDS->AddEdge(n1, n2, n12);
 588   }
 589   else {
 590     if(id)
 591       edge = meshDS->AddEdgeWithID(n1, n2, id);
 592     else
 593       edge = meshDS->AddEdge(n1, n2);
 594   }
 595
 596   if ( mySetElemOnShape && myShapeID > 0 )
 597     meshDS->SetMeshElementOnShape( edge, myShapeID );
 598
 599   return edge;
 600 }
 601
 602 //=======================================================================
 603 /*!
 604  * Creates quadratic or linear triangle
 605  */
 606 //=======================================================================
 607
 608 SMDS_MeshFace* SMESH_MesherHelper::AddFace(const SMDS_MeshNode* n1,
 609                                            const SMDS_MeshNode* n2,
 610                                            const SMDS_MeshNode* n3,
 611                                            const int id,
 612                                            const bool force3d)
 613 {
 614   SMESHDS_Mesh * meshDS = GetMeshDS();
 615   SMDS_MeshFace* elem = 0;
 616   if(!myCreateQuadratic) {
 617     if(id)
 618       elem = meshDS->AddFaceWithID(n1, n2, n3, id);
 619     else
 620       elem = meshDS->AddFace(n1, n2, n3);
 621   }
 622   else {
 623     const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
 624     const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,force3d);
 625     const SMDS_MeshNode* n31 = GetMediumNode(n3,n1,force3d);
 626
 627     if(id)
 628       elem = meshDS->AddFaceWithID(n1, n2, n3, n12, n23, n31, id);
 629     else
 630       elem = meshDS->AddFace(n1, n2, n3, n12, n23, n31);
 631   }
 632   if ( mySetElemOnShape && myShapeID > 0 )
 633     meshDS->SetMeshElementOnShape( elem, myShapeID );
 634
 635   return elem;
 636 }
 637
 638 //=======================================================================
 639 /*!
 640  * Creates quadratic or linear quadrangle
 641  */
 642 //=======================================================================
 643
 644 SMDS_MeshFace* SMESH_MesherHelper::AddFace(const SMDS_MeshNode* n1,
 645                                            const SMDS_MeshNode* n2,
 646                                            const SMDS_MeshNode* n3,
 647                                            const SMDS_MeshNode* n4,
 648                                            const int id,
 649                                            const bool force3d)
 650 {
 651   SMESHDS_Mesh * meshDS = GetMeshDS();
 652   SMDS_MeshFace* elem = 0;
 653   if(!myCreateQuadratic) {
 654     if(id)
 655       elem = meshDS->AddFaceWithID(n1, n2, n3, n4, id);
 656     else
 657       elem = meshDS->AddFace(n1, n2, n3, n4);
 658   }
 659   else {
 660     const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
 661     const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,force3d);
 662     const SMDS_MeshNode* n34 = GetMediumNode(n3,n4,force3d);
 663     const SMDS_MeshNode* n41 = GetMediumNode(n4,n1,force3d);
 664
 665     if(id)
 666       elem = meshDS->AddFaceWithID(n1, n2, n3, n4, n12, n23, n34, n41, id);
 667     else
 668       elem = meshDS->AddFace(n1, n2, n3, n4, n12, n23, n34, n41);
 669   }
 670   if ( mySetElemOnShape && myShapeID > 0 )
 671     meshDS->SetMeshElementOnShape( elem, myShapeID );
 672
 673   return elem;
 674 }
 675
 676 //=======================================================================
 677 /*!
 678  * Creates quadratic or linear volume
 679  */
 680 //=======================================================================
 681
 682 SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1,
 683                                                const SMDS_MeshNode* n2,
 684                                                const SMDS_MeshNode* n3,
 685                                                const SMDS_MeshNode* n4,
 686                                                const SMDS_MeshNode* n5,
 687                                                const SMDS_MeshNode* n6,
 688                                                const int id,
 689                                                const bool force3d)
 690 {
 691   SMESHDS_Mesh * meshDS = GetMeshDS();
 692   SMDS_MeshVolume* elem = 0;
 693   if(!myCreateQuadratic) {
 694     if(id)
 695       elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, n5, n6, id);
 696     else
 697       elem = meshDS->AddVolume(n1, n2, n3, n4, n5, n6);
 698   }
 699   else {
 700     const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
 701     const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,force3d);
 702     const SMDS_MeshNode* n31 = GetMediumNode(n3,n1,force3d);
 703
 704     const SMDS_MeshNode* n45 = GetMediumNode(n4,n5,force3d);
 705     const SMDS_MeshNode* n56 = GetMediumNode(n5,n6,force3d);
 706     const SMDS_MeshNode* n64 = GetMediumNode(n6,n4,force3d);
 707
 708     const SMDS_MeshNode* n14 = GetMediumNode(n1,n4,force3d);
 709     const SMDS_MeshNode* n25 = GetMediumNode(n2,n5,force3d);
 710     const SMDS_MeshNode* n36 = GetMediumNode(n3,n6,force3d);
 711
 712     if(id)
 713       elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, n5, n6,
 714                                      n12, n23, n31, n45, n56, n64, n14, n25, n36, id);
 715     else
 716       elem = meshDS->AddVolume(n1, n2, n3, n4, n5, n6,
 717                                n12, n23, n31, n45, n56, n64, n14, n25, n36);
 718   }
 719   if ( mySetElemOnShape && myShapeID > 0 )
 720     meshDS->SetMeshElementOnShape( elem, myShapeID );
 721
 722   return elem;
 723 }
 724
 725 //=======================================================================
 726 /*!
 727  * Creates quadratic or linear volume
 728  */
 729 //=======================================================================
 730
 731 SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1,
 732                                                const SMDS_MeshNode* n2,
 733                                                const SMDS_MeshNode* n3,
 734                                                const SMDS_MeshNode* n4,
 735                                                const int id,
 736                                                const bool force3d)
 737 {
 738   SMESHDS_Mesh * meshDS = GetMeshDS();
 739   SMDS_MeshVolume* elem = 0;
 740   if(!myCreateQuadratic) {
 741     if(id)
 742       elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, id);
 743     else
 744       elem = meshDS->AddVolume(n1, n2, n3, n4);
 745   }
 746   else {
 747     const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
 748     const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,force3d);
 749     const SMDS_MeshNode* n31 = GetMediumNode(n3,n1,force3d);
 750
 751     const SMDS_MeshNode* n14 = GetMediumNode(n1,n4,force3d);
 752     const SMDS_MeshNode* n24 = GetMediumNode(n2,n4,force3d);
 753     const SMDS_MeshNode* n34 = GetMediumNode(n3,n4,force3d);
 754
 755     if(id)
 756       elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, n12, n23, n31, n14, n24, n34, id);
 757     else
 758       elem = meshDS->AddVolume(n1, n2, n3, n4, n12, n23, n31, n14, n24, n34);
 759   }
 760   if ( mySetElemOnShape && myShapeID > 0 )
 761     meshDS->SetMeshElementOnShape( elem, myShapeID );
 762
 763   return elem;
 764 }
 765
 766 //=======================================================================
 767 /*!
 768  * Creates quadratic or linear pyramid
 769  */
 770 //=======================================================================
 771
 772 SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1,
 773                                                const SMDS_MeshNode* n2,
 774                                                const SMDS_MeshNode* n3,
 775                                                const SMDS_MeshNode* n4,
 776                                                const SMDS_MeshNode* n5,
 777                                                const int id,
 778                                                const bool force3d)
 779 {
 780   SMDS_MeshVolume* elem = 0;
 781   if(!myCreateQuadratic) {
 782     if(id)
 783       elem = GetMeshDS()->AddVolumeWithID(n1, n2, n3, n4, n5, id);
 784     else
 785       elem = GetMeshDS()->AddVolume(n1, n2, n3, n4, n5);
 786   }
 787   else {
 788     const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
 789     const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,force3d);
 790     const SMDS_MeshNode* n34 = GetMediumNode(n3,n4,force3d);
 791     const SMDS_MeshNode* n41 = GetMediumNode(n4,n1,force3d);
 792
 793     const SMDS_MeshNode* n15 = GetMediumNode(n1,n5,force3d);
 794     const SMDS_MeshNode* n25 = GetMediumNode(n2,n5,force3d);
 795     const SMDS_MeshNode* n35 = GetMediumNode(n3,n5,force3d);
 796     const SMDS_MeshNode* n45 = GetMediumNode(n4,n5,force3d);
 797
 798     if(id)
 799       elem = GetMeshDS()->AddVolumeWithID ( n1,  n2,  n3,  n4,  n5,
 800                                             n12, n23, n34, n41,
 801                                             n15, n25, n35, n45,
 802                                             id);
 803     else
 804       elem = GetMeshDS()->AddVolume( n1,  n2,  n3,  n4,  n5,
 805                                      n12, n23, n34, n41,
 806                                      n15, n25, n35, n45);
 807   }
 808   if ( mySetElemOnShape && myShapeID > 0 )
 809     GetMeshDS()->SetMeshElementOnShape( elem, myShapeID );
 810
 811   return elem;
 812 }
 813
 814 //=======================================================================
 815 /*!
 816  * Creates quadratic or linear hexahedron
 817  */
 818 //=======================================================================
 819
 820 SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1,
 821                                                const SMDS_MeshNode* n2,
 822                                                const SMDS_MeshNode* n3,
 823                                                const SMDS_MeshNode* n4,
 824                                                const SMDS_MeshNode* n5,
 825                                                const SMDS_MeshNode* n6,
 826                                                const SMDS_MeshNode* n7,
 827                                                const SMDS_MeshNode* n8,
 828                                                const int id,
 829                                                const bool force3d)
 830 {
 831   SMESHDS_Mesh * meshDS = GetMeshDS();
 832   SMDS_MeshVolume* elem = 0;
 833   if(!myCreateQuadratic) {
 834     if(id)
 835       elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, n5, n6, n7, n8, id);
 836     else
 837       elem = meshDS->AddVolume(n1, n2, n3, n4, n5, n6, n7, n8);
 838   }
 839   else {
 840     const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
 841     const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,force3d);
 842     const SMDS_MeshNode* n34 = GetMediumNode(n3,n4,force3d);
 843     const SMDS_MeshNode* n41 = GetMediumNode(n4,n1,force3d);
 844
 845     const SMDS_MeshNode* n56 = GetMediumNode(n5,n6,force3d);
 846     const SMDS_MeshNode* n67 = GetMediumNode(n6,n7,force3d);
 847     const SMDS_MeshNode* n78 = GetMediumNode(n7,n8,force3d);
 848     const SMDS_MeshNode* n85 = GetMediumNode(n8,n5,force3d);
 849
 850     const SMDS_MeshNode* n15 = GetMediumNode(n1,n5,force3d);
 851     const SMDS_MeshNode* n26 = GetMediumNode(n2,n6,force3d);
 852     const SMDS_MeshNode* n37 = GetMediumNode(n3,n7,force3d);
 853     const SMDS_MeshNode* n48 = GetMediumNode(n4,n8,force3d);
 854
 855     if(id)
 856       elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, n5, n6, n7, n8,
 857                                      n12, n23, n34, n41, n56, n67,
 858                                      n78, n85, n15, n26, n37, n48, id);
 859     else
 860       elem = meshDS->AddVolume(n1, n2, n3, n4, n5, n6, n7, n8,
 861                                n12, n23, n34, n41, n56, n67,
 862                                n78, n85, n15, n26, n37, n48);
 863   }
 864   if ( mySetElemOnShape && myShapeID > 0 )
 865     meshDS->SetMeshElementOnShape( elem, myShapeID );
 866
 867   return elem;
 868 }
 869
 870 //=======================================================================
 871   /*!
 872    * \brief Load nodes bound to face into a map of node columns
 873     * \param theParam2ColumnMap - map of node columns to fill
 874     * \param theFace - the face on which nodes are searched for
 875     * \param theBaseEdge - the edge nodes of which are columns' bases
 876     * \param theMesh - the mesh containing nodes
 877     * \retval bool - false if something is wrong
 878    *
 879    * The key of the map is a normalized parameter of each
 880    * base node on theBaseEdge.
 881    * This method works in supposition that nodes on the face
 882    * forms a rectangular grid and elements can be quardrangles or triangles
 883    */
 884 //=======================================================================
 885
 886 bool SMESH_MesherHelper::LoadNodeColumns(TParam2ColumnMap & theParam2ColumnMap,
 887                                          const TopoDS_Face& theFace,
 888                                          const TopoDS_Edge& theBaseEdge,
 889                                          SMESHDS_Mesh*      theMesh)
 890 {
 891   // get vertices of theBaseEdge
 892   TopoDS_Vertex vfb, vlb, vft; // first and last, bottom and top vertices
 893   TopoDS_Edge eFrw = TopoDS::Edge( theBaseEdge.Oriented( TopAbs_FORWARD ));
 894   TopExp::Vertices( eFrw, vfb, vlb );
 895
 896   // find the other edges of theFace and orientation of e1
 897   TopoDS_Edge e1, e2, eTop;
 898   bool rev1, CumOri = false;
 899   TopExp_Explorer exp( theFace, TopAbs_EDGE );
 900   int nbEdges = 0;
 901   for ( ; exp.More(); exp.Next() ) {
 902     if ( ++nbEdges > 4 ) {
 903       return false; // more than 4 edges in theFace
 904     }
 905     TopoDS_Edge e = TopoDS::Edge( exp.Current() );
 906     if ( theBaseEdge.IsSame( e ))
 907       continue;
 908     TopoDS_Vertex vCommon;
 909     if ( !TopExp::CommonVertex( theBaseEdge, e, vCommon ))
 910       eTop = e;
 911     else if ( vCommon.IsSame( vfb )) {
 912       e1 = e;
 913       vft = TopExp::LastVertex( e1, CumOri );
 914       rev1 = vfb.IsSame( vft );
 915       if ( rev1 )
 916         vft = TopExp::FirstVertex( e1, CumOri );
 917     }
 918     else
 919       e2 = e;
 920   }
 921   if ( nbEdges < 4 ) {
 922     return false; // less than 4 edges in theFace
 923   }
 924   if ( e2.IsNull() && vfb.IsSame( vlb ))
 925     e2 = e1;
 926
 927   // submeshes corresponding to shapes
 928   SMESHDS_SubMesh* smFace = theMesh->MeshElements( theFace );
 929   SMESHDS_SubMesh* smb = theMesh->MeshElements( theBaseEdge );
 930   SMESHDS_SubMesh* smt = theMesh->MeshElements( eTop );
 931   SMESHDS_SubMesh* sm1 = theMesh->MeshElements( e1 );
 932   SMESHDS_SubMesh* sm2 = theMesh->MeshElements( e2 );
 933   SMESHDS_SubMesh* smVfb = theMesh->MeshElements( vfb );
 934   SMESHDS_SubMesh* smVlb = theMesh->MeshElements( vlb );
 935   SMESHDS_SubMesh* smVft = theMesh->MeshElements( vft );
 936   if (!smFace || !smb || !smt || !sm1 || !sm2 || !smVfb || !smVlb || !smVft ) {
 937     RETURN_BAD_RESULT( "NULL submesh " <<smFace<<" "<<smb<<" "<<smt<<" "<<
 938                        sm1<<" "<<sm2<<" "<<smVfb<<" "<<smVlb<<" "<<smVft);
 939   }
 940   if ( smb->NbNodes() != smt->NbNodes() || sm1->NbNodes() != sm2->NbNodes() ) {
 941     RETURN_BAD_RESULT(" Diff nb of nodes on opposite edges" );
 942   }
 943   if (smVfb->NbNodes() != 1 || smVlb->NbNodes() != 1 || smVft->NbNodes() != 1) {
 944     RETURN_BAD_RESULT("Empty submesh of vertex");
 945   }
 946   // define whether mesh is quadratic
 947   bool isQuadraticMesh = false;
 948   SMDS_ElemIteratorPtr eIt = smFace->GetElements();
 949   if ( !eIt->more() ) {
 950     RETURN_BAD_RESULT("No elements on the face");
 951   }
 952   const SMDS_MeshElement* e = eIt->next();
 953   isQuadraticMesh = e->IsQuadratic();
 954
 955   if ( sm1->NbNodes() * smb->NbNodes() != smFace->NbNodes() ) {
 956     // check quadratic case
 957     if ( isQuadraticMesh ) {
 958       // what if there are quadrangles and triangles mixed?
 959 //       int n1 = sm1->NbNodes()/2;
 960 //       int n2 = smb->NbNodes()/2;
 961 //       int n3 = sm1->NbNodes() - n1;
 962 //       int n4 = smb->NbNodes() - n2;
 963 //       int nf = sm1->NbNodes()*smb->NbNodes() - n3*n4;
 964 //       if( nf != smFace->NbNodes() ) {
 965 //         MESSAGE( "Wrong nb face nodes: " <<
 966 //                 sm1->NbNodes()<<" "<<smb->NbNodes()<<" "<<smFace->NbNodes());
 967 //         return false;
 968 //       }
 969     }
 970     else {
 971       RETURN_BAD_RESULT( "Wrong nb face nodes: " <<
 972                          sm1->NbNodes()<<" "<<smb->NbNodes()<<" "<<smFace->NbNodes());
 973     }
 974   }
 975   // IJ size
 976   int vsize = sm1->NbNodes() + 2;
 977   int hsize = smb->NbNodes() + 2;
 978   if(isQuadraticMesh) {
 979     vsize = vsize - sm1->NbNodes()/2 -1;
 980     hsize = hsize - smb->NbNodes()/2 -1;
 981   }
 982
 983   // load nodes from theBaseEdge
 984
 985   set<const SMDS_MeshNode*> loadedNodes;
 986   const SMDS_MeshNode* nullNode = 0;
 987
 988   vector<const SMDS_MeshNode*> & nVecf = theParam2ColumnMap[ 0.];
 989   nVecf.resize( vsize, nullNode );
 990   loadedNodes.insert( nVecf[ 0 ] = smVfb->GetNodes()->next() );
 991
 992   vector<const SMDS_MeshNode*> & nVecl = theParam2ColumnMap[ 1.];
 993   nVecl.resize( vsize, nullNode );
 994   loadedNodes.insert( nVecl[ 0 ] = smVlb->GetNodes()->next() );
 995
 996   double f, l;
 997   BRep_Tool::Range( eFrw, f, l );
 998   double range = l - f;
 999   SMDS_NodeIteratorPtr nIt = smb->GetNodes();
1000   const SMDS_MeshNode* node;
1001   while ( nIt->more() ) {
1002     node = nIt->next();
1003     if(IsMedium(node, SMDSAbs_Edge))
1004       continue;
1005     const SMDS_EdgePosition* pos =
1006       dynamic_cast<const SMDS_EdgePosition*>( node->GetPosition().get() );
1007     if ( !pos ) {
1008       return false;
1009     }
1010     double u = ( pos->GetUParameter() - f ) / range;
1011     vector<const SMDS_MeshNode*> & nVec = theParam2ColumnMap[ u ];
1012     nVec.resize( vsize, nullNode );
1013     loadedNodes.insert( nVec[ 0 ] = node );
1014   }
1015   if ( theParam2ColumnMap.size() != hsize ) {
1016     RETURN_BAD_RESULT( "Wrong node positions on theBaseEdge" );
1017   }
1018
1019   // load nodes from e1
1020
1021   map< double, const SMDS_MeshNode*> sortedNodes; // sort by param on edge
1022   nIt = sm1->GetNodes();
1023   while ( nIt->more() ) {
1024     node = nIt->next();
1025     if(IsMedium(node))
1026       continue;
1027     const SMDS_EdgePosition* pos =
1028       dynamic_cast<const SMDS_EdgePosition*>( node->GetPosition().get() );
1029     if ( !pos ) {
1030       return false;
1031     }
1032     sortedNodes.insert( make_pair( pos->GetUParameter(), node ));
1033   }
1034   loadedNodes.insert( nVecf[ vsize - 1 ] = smVft->GetNodes()->next() );
1035   map< double, const SMDS_MeshNode*>::iterator u_n = sortedNodes.begin();
1036   int row  = rev1 ? vsize - 1 : 0;
1037   int dRow = rev1 ? -1 : +1;
1038   for ( ; u_n != sortedNodes.end(); u_n++ ) {
1039     row += dRow;
1040     loadedNodes.insert( nVecf[ row ] = u_n->second );
1041   }
1042
1043   // try to load the rest nodes
1044
1045   // get all faces from theFace
1046   TIDSortedElemSet allFaces, foundFaces;
1047   eIt = smFace->GetElements();
1048   while ( eIt->more() ) {
1049     const SMDS_MeshElement* e = eIt->next();
1050     if ( e->GetType() == SMDSAbs_Face )
1051       allFaces.insert( e );
1052   }
1053   // Starting from 2 neighbour nodes on theBaseEdge, look for a face
1054   // the nodes belong to, and between the nodes of the found face,
1055   // look for a not loaded node considering this node to be the next
1056   // in a column of the starting second node. Repeat, starting
1057   // from nodes next to the previous starting nodes in their columns,
1058   // and so on while a face can be found. Then go the the next pair
1059   // of nodes on theBaseEdge.
1060   TParam2ColumnMap::iterator par_nVec_1 = theParam2ColumnMap.begin();
1061   TParam2ColumnMap::iterator par_nVec_2 = par_nVec_1;
1062   // loop on columns
1063   int col = 0;
1064   for ( par_nVec_2++; par_nVec_2 != theParam2ColumnMap.end(); par_nVec_1++, par_nVec_2++ ) {
1065     col++;
1066     row = 0;
1067     const SMDS_MeshNode* n1 = par_nVec_1->second[ row ];
1068     const SMDS_MeshNode* n2 = par_nVec_2->second[ row ];
1069     const SMDS_MeshElement* face = 0;
1070     bool lastColOnClosedFace = ( nVecf[ row ] == n2 );
1071     do {
1072       // look for a face by 2 nodes
1073       face = SMESH_MeshEditor::FindFaceInSet( n1, n2, allFaces, foundFaces );
1074       if ( face ) {
1075         int nbFaceNodes = face->NbNodes();
1076         if ( face->IsQuadratic() )
1077           nbFaceNodes /= 2;
1078         if ( nbFaceNodes>4 ) {
1079           RETURN_BAD_RESULT(" Too many nodes in a face: " << nbFaceNodes );
1080         }
1081         // look for a not loaded node of the <face>
1082         bool found = false;
1083         const SMDS_MeshNode* n3 = 0; // a node defferent from n1 and n2
1084         for ( int i = 0; i < nbFaceNodes && !found; ++i ) {
1085           node = face->GetNode( i );
1086           found = loadedNodes.insert( node ).second;
1087           if ( !found && node != n1 && node != n2 )
1088             n3 = node;
1089         }
1090         if ( lastColOnClosedFace && row + 1 < vsize ) {
1091           node = nVecf[ row + 1 ];
1092           found = ( face->GetNodeIndex( node ) >= 0 );
1093         }
1094         if ( found ) {
1095           if ( ++row > vsize - 1 ) {
1096             RETURN_BAD_RESULT( "Too many nodes in column "<< col <<": "<< row+1);
1097           }
1098           par_nVec_2->second[ row ] = node;
1099           foundFaces.insert( face );
1100           n2 = node;
1101           if ( nbFaceNodes==4 ) {
1102             n1 = par_nVec_1->second[ row ];
1103           }
1104         }
1105         else if ( nbFaceNodes==3 && n3 == par_nVec_1->second[ row + 1 ] ) {
1106           n1 = n3;
1107         }
1108         else  {
1109           RETURN_BAD_RESULT( "Not quad mesh, column "<< col );
1110         }
1111       }
1112     }
1113     while ( face && n1 && n2 );
1114
1115     if ( row < vsize - 1 ) {
1116       MESSAGE( "Too few nodes in column "<< col <<": "<< row+1);
1117       MESSAGE( "Base node 1: "<< par_nVec_1->second[0]);
1118       MESSAGE( "Base node 2: "<< par_nVec_2->second[0]);
1119       if ( n1 ) { MESSAGE( "Current node 1: "<< n1); }
1120       else      { MESSAGE( "Current node 1: NULL");  }
1121       if ( n2 ) { MESSAGE( "Current node 2: "<< n2); }
1122       else      { MESSAGE( "Current node 2: NULL");  }
1123       MESSAGE( "first base node: "<< theParam2ColumnMap.begin()->second[0]);
1124       MESSAGE( "last base node: "<< theParam2ColumnMap.rbegin()->second[0]);
1125       return false;
1126     }
1127   } // loop on columns
1128
1129   return true;
1130 }