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