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