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