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