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