src/SMESH/SMESH_MesherHelper.cxx

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