src/SMESH/SMESH_MesherHelper.cxx

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