src/SMESH/SMESH_MesherHelper.cxx

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