src/StdMeshers/StdMeshers_Cartesian_VL.cxx

   1 // Copyright (C) 2007-2021  CEA/DEN, EDF R&D, OPEN CASCADE
   2 //
   3 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   5 //
   6 // This library is free software; you can redistribute it and/or
   7 // modify it under the terms of the GNU Lesser General Public
   8 // License as published by the Free Software Foundation; either
   9 // version 2.1 of the License, or (at your option) any later version.
  10 //
  11 // This library is distributed in the hope that it will be useful,
  12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 // Lesser General Public License for more details.
  15 //
  16 // You should have received a copy of the GNU Lesser General Public
  17 // License along with this library; if not, write to the Free Software
  18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  19 //
  20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  21 //
  22 // File      : StdMeshers_Cartesian_VL.cxx
  23 // Created   : Tue May 24 13:03:09 2022
  24 // Author    : Edward AGAPOV (eap)
  25
  26 #include "StdMeshers_Cartesian_VL.hxx"
  27 #include "StdMeshers_ViscousLayers.hxx"
  28
  29 #include <SMDS_MeshGroup.hxx>
  30 #include <SMESHDS_Mesh.hxx>
  31 #include <SMESHDS_SubMesh.hxx>
  32 #include <SMESH_Algo.hxx>
  33 #include <SMESH_Mesh.hxx>
  34 #include <SMESH_MeshEditor.hxx>
  35 #include <SMESH_MesherHelper.hxx>
  36 #include <SMESH_TypeDefs.hxx>
  37 #include <SMESH_subMesh.hxx>
  38
  39 #include <BRepAdaptor_Curve.hxx>
  40 #include <BRepTopAdaptor_FClass2d.hxx>
  41 #include <BRep_Builder.hxx>
  42 #include <BRep_Tool.hxx>
  43 #include <ShapeAnalysis_Curve.hxx>
  44 #include <ShapeAnalysis_Surface.hxx>
  45 #include <TopExp.hxx>
  46 #include <TopExp_Explorer.hxx>
  47 #include <TopoDS.hxx>
  48
  49 using namespace StdMeshers_Cartesian_VL;
  50
  51 namespace
  52 {
  53   typedef int                     TGeomID; // IDs of sub-shapes
  54
  55   /*!
  56    * \brief Temporary mesh
  57    */
  58   struct TmpMesh: public SMESH_Mesh
  59   {
  60     TmpMesh() {
  61       _isShapeToMesh = (_id = 0);
  62       _meshDS  = new SMESHDS_Mesh( _id, true );
  63     }
  64   };
  65   // --------------------------------------------------------------------------
  66   /*!
  67    * \brief Edge of viscous layers; goes from a grid node by normal to boundary
  68    */
  69   struct VLEdge
  70   {
  71     std::vector< SMESH_NodeXYZ > _nodes;
  72     double                       _uv[2]; // position of TgtNode() on geometry
  73     double                       _length;
  74
  75     const SMDS_MeshNode* TgtNode() const { return _nodes.back().Node(); }
  76   };
  77   typedef NCollection_DataMap< const SMDS_MeshNode*, VLEdge* > TNode2VLEdge;
  78   // --------------------------------------------------------------------------
  79   /*!
  80    * \brief VLEdge's of one shape
  81    */
  82   struct VLEdgesOnShape
  83   {
  84     std::vector< VLEdge > _edges;
  85     TopoDS_Shape          _initShape;
  86     TopoDS_Shape          _offsetShape;
  87     int                   _initShapeID;
  88     bool                  _hasVL;
  89     bool                  _toCheckCoinc; // to check if nodes on shape boundary coincide
  90   };
  91
  92   //================================================================================
  93   /*!
  94    * \brief Project a point on offset FACE to EDGEs of an initial FACE
  95    *  \param [in] offP - point to project
  96    *  \param [in] initFace - FACE to project on
  97    *  \return gp_Pnt - projection point
  98    */
  99   //================================================================================
 100
 101   gp_Pnt projectToWire( const gp_Pnt&      offP,
 102                         const TopoDS_Face& initFace,
 103                         gp_Pnt2d &         uv )
 104   {
 105     double   minDist = Precision::Infinite();
 106     const double tol = Precision::Confusion();
 107     gp_Pnt              proj, projction;
 108     Standard_Real       param;
 109     ShapeAnalysis_Curve projector;
 110
 111     for ( TopExp_Explorer eExp( initFace, TopAbs_EDGE ); eExp.More(); eExp.Next() )
 112     {
 113       BRepAdaptor_Curve initCurve( TopoDS::Edge( eExp.Current() ));
 114       //const double f = initCurve.FirstParameter(), l = initCurve.LastParameter();
 115       double dist = projector.Project( initCurve, offP, tol, proj, param, /*adjustToEnds=*/false );
 116       if ( dist < minDist )
 117       {
 118         projction = proj;
 119         minDist = dist;
 120       }
 121     }
 122     uv.SetCoord(0.,0.); // !!!!!!!
 123     return projction;
 124   }
 125
 126   //================================================================================
 127   /*!
 128    * \brief Project nodes from offset FACE to initial FACE
 129    *  \param [inout] theEOS - VL edges on a geom FACE
 130    *  \param [inout] theOffsetMDS - offset mesh to fill in
 131    */
 132   //================================================================================
 133
 134   void projectToFace( VLEdgesOnShape & theEOS,
 135                       SMESHDS_Mesh*    theOffsetMDS,
 136                       TNode2VLEdge &   theN2E )
 137   {
 138     SMESHDS_SubMesh* sm = theOffsetMDS->MeshElements( theEOS._offsetShape );
 139     if ( !sm || sm->NbElements() == 0 || sm->NbNodes() == 0 )
 140       return;
 141     theEOS._edges.resize( sm->NbNodes() );
 142
 143     const TopoDS_Face& initFace = TopoDS::Face( theEOS._initShape );
 144     ShapeAnalysis_Surface projector( BRep_Tool::Surface( initFace ));
 145     const double clsfTol = 1e2 * BRep_Tool::MaxTolerance( initFace, TopAbs_VERTEX );
 146     BRepTopAdaptor_FClass2d classifier( initFace, clsfTol );
 147
 148     const double        tol = Precision::Confusion();
 149     //const double f = initCurve.FirstParameter(), l = initCurve.LastParameter();
 150     gp_Pnt              proj;
 151
 152     int iN = 0;
 153     for ( SMDS_NodeIteratorPtr nIt = sm->GetNodes(); nIt->more(); ++iN )
 154     {
 155       SMESH_NodeXYZ offP = nIt->next();
 156       gp_Pnt2d uv = projector.ValueOfUV( offP, tol );
 157       TopAbs_State st = classifier.Perform( uv );
 158       if ( st == TopAbs_IN || st == TopAbs_ON )
 159       {
 160         proj = projector.Value( uv );
 161       }
 162       else
 163       {
 164         proj = projectToWire( offP, initFace, uv );
 165       }
 166       if ( st == TopAbs_ON || st == TopAbs_OUT )
 167         theEOS._toCheckCoinc = true;
 168
 169       VLEdge & vlEdge = theEOS._edges[ iN ];
 170       vlEdge._nodes.push_back( offP.Node() );
 171       vlEdge._nodes.push_back( theOffsetMDS->AddNode( proj.X(), proj.Y(), proj.Z() ));
 172       vlEdge._uv[0] = uv.X();
 173       vlEdge._uv[1] = uv.Y();
 174       //vlEdge._length = proj.Distance( offP );
 175
 176       theN2E.Bind( offP.Node(), &vlEdge );
 177     }
 178     return;
 179
 180   }
 181
 182   //================================================================================
 183   /*!
 184    * \brief Project nodes from offset EDGE to initial EDGE
 185    *  \param [inout] theEOS - VL edges on a geom EDGE
 186    *  \param [inout] theOffsetMDS - offset mesh to add new nodes to
 187    */
 188   //================================================================================
 189
 190   void projectToEdge( VLEdgesOnShape & theEOS,
 191                       SMESHDS_Mesh*    theOffsetMDS,
 192                       TNode2VLEdge &   theN2E )
 193   {
 194     SMESHDS_SubMesh* sm = theOffsetMDS->MeshElements( theEOS._offsetShape );
 195     if ( !sm || sm->NbElements() == 0 )
 196       return;
 197     theEOS._edges.resize( sm->NbNodes() );
 198
 199     ShapeAnalysis_Curve projector;
 200     BRepAdaptor_Curve   initCurve( TopoDS::Edge( theEOS._initShape ));
 201     const double        tol = Precision::Confusion();
 202     const double f = initCurve.FirstParameter(), l = initCurve.LastParameter();
 203     gp_Pnt              proj;
 204     Standard_Real       param;
 205
 206     int iN = 0;
 207     for ( SMDS_NodeIteratorPtr nIt = sm->GetNodes(); nIt->more(); ++iN )
 208     {
 209       SMESH_NodeXYZ offP = nIt->next();
 210       double dist = projector.Project( initCurve, offP, tol, proj, param, /*adjustToEnds=*/false );
 211       bool paramOK = ( f < param && param < l );
 212       if ( !paramOK )
 213       {
 214         if ( param < f )
 215           param = f;
 216         else if ( param > l )
 217           param = l;
 218         theEOS._toCheckCoinc = true;
 219       }
 220       proj = initCurve.Value( param );
 221
 222       VLEdge & vlEdge = theEOS._edges[ iN ];
 223       vlEdge._nodes.push_back( offP.Node() );
 224       vlEdge._nodes.push_back( theOffsetMDS->AddNode( proj.X(), proj.Y(), proj.Z() ));
 225       vlEdge._uv[0] = param;
 226       vlEdge._length = paramOK ? dist : proj.Distance( offP );
 227
 228       theN2E.Bind( offP.Node(), &vlEdge );
 229     }
 230     return;
 231   }
 232
 233   //================================================================================
 234   /*!
 235    * \brief Compute heights of viscous layers and finds FACEs with VL
 236    *  \param [in] hyp - viscous layers hypothesis
 237    *  \param [in] mesh - the main mesh
 238    *  \param [in] shape - the main shape
 239    *  \param [out] vlH - heights of viscous layers to compute
 240    *  \param [out] shapesWVL - IDs of shapes with VL
 241    *  \return bool - isOK
 242    */
 243   //================================================================================
 244
 245   void computeVLHeight( const StdMeshers_ViscousLayers* hyp,
 246                         std::vector< double > &         vlH )
 247   {
 248     const double  T = hyp->GetTotalThickness();
 249     const double  f = hyp->GetStretchFactor();
 250     const int     N = hyp->GetNumberLayers();
 251     const double h0 = hyp->Get1stLayerThickness( T, f, N );
 252
 253     vlH.reserve( hyp->GetNumberLayers() );
 254     vlH.push_back( h0 );
 255     for ( double h = h0 * f; (int)vlH.size() < N-1; h *= f )
 256       vlH.push_back( h + vlH.back() );
 257     vlH.push_back( T );
 258   }
 259
 260   //================================================================================
 261   /*!
 262    * \brief Create intermediate nodes on VLEdge
 263    *  \param [inout] vlEdge - VLEdge to divide
 264    *  \param [in] vlH - thicknesses of layers
 265    *  \param [inout] mesh - the mesh no fill in
 266    */
 267   //================================================================================
 268
 269   void divideVLEdge( VLEdge*                      vlEdge,
 270                      const std::vector< double >& vlH,
 271                      SMESHDS_Mesh*                mesh )
 272   {
 273     SMESH_NodeXYZ lastNode = vlEdge->_nodes.back();
 274     SMESH_NodeXYZ frstNode = vlEdge->_nodes[0];
 275     gp_XYZ dir = frstNode - lastNode;
 276
 277     vlEdge->_nodes.resize( vlH.size() + 1 );
 278
 279     for ( size_t i = 1; i < vlH.size(); ++i )
 280     {
 281       double r = vlH[ i-1 ] / vlH.back();
 282       gp_XYZ p = lastNode + dir * r;
 283       vlEdge->_nodes[ vlH.size() - i ] = mesh->AddNode( p.X(), p.Y(), p.Z() );
 284     }
 285     vlEdge->_nodes.back() = lastNode;
 286   }
 287
 288   //================================================================================
 289   /*!
 290    * \brief Create a polyhedron from nodes of VLEdge's
 291    *  \param [in] edgesVec - the VLEdge's
 292    *  \param [in] vNodes - node buffer
 293    *  \param [in] editor - editor of offset mesh
 294    */
 295   //================================================================================
 296
 297   bool makePolyhedron( const std::vector< VLEdge*> &         edgesVec,
 298                        std::vector< const SMDS_MeshNode* > & vNodes,
 299                        SMESH_MeshEditor&                     editor,
 300                        SMESH_MeshEditor::ElemFeatures        elemType)
 301   {
 302     elemType.SetPoly( true );
 303     elemType.myPolyhedQuantities.clear();
 304     elemType.myNodes.clear();
 305
 306     // add facets of walls
 307     size_t nbBaseNodes = edgesVec.size();
 308     for ( size_t iN = 0; iN < nbBaseNodes; ++iN )
 309     {
 310       VLEdge* e0 = edgesVec[ iN ];
 311       VLEdge* e1 = edgesVec[( iN + 1 ) % nbBaseNodes ];
 312       size_t nbN0 = e0->_nodes.size();
 313       size_t nbN1 = e1->_nodes.size();
 314       if ( nbN0 == nbN1 )
 315       {
 316         for ( size_t i = 1; i < nbN0; ++i )
 317         {
 318           elemType.myNodes.push_back( e1->_nodes[ i - 1 ].Node());
 319           elemType.myNodes.push_back( e1->_nodes[ i     ].Node());
 320           elemType.myNodes.push_back( e0->_nodes[ i     ].Node());
 321           elemType.myNodes.push_back( e0->_nodes[ i - 1 ].Node());
 322           elemType.myPolyhedQuantities.push_back( 4 );
 323         }
 324       }
 325       else
 326       {
 327         for ( size_t i = 0; i < nbN1; ++i )
 328           elemType.myNodes.push_back( e1->_nodes[ i ].Node() );
 329         for ( size_t i = 0; i < nbN0; ++i )
 330           elemType.myNodes.push_back( e0->_nodes[ nbN0 - 1 - i ].Node() );
 331         elemType.myPolyhedQuantities.push_back( nbN0 + nbN1 );
 332       }
 333     }
 334
 335     // add facets of top
 336     vNodes.clear();
 337     for ( size_t iN = 0; iN < nbBaseNodes; ++iN )
 338     {
 339       VLEdge* e0 = edgesVec[ iN ];
 340       elemType.myNodes.push_back( e0->_nodes.back().Node() );
 341       vNodes.push_back( e0->_nodes[ 0 ].Node());
 342     }
 343     elemType.myPolyhedQuantities.push_back( nbBaseNodes );
 344
 345     // add facets of bottom
 346     elemType.myNodes.insert( elemType.myNodes.end(), vNodes.rbegin(), vNodes.rend() );
 347     elemType.myPolyhedQuantities.push_back( nbBaseNodes );
 348
 349     const SMDS_MeshElement* vol = editor.AddElement( elemType.myNodes, elemType );
 350     vol->setIsMarked( true ); // to add to group
 351
 352     return vol;
 353   }
 354
 355   //================================================================================
 356   /*!
 357    * \brief Transform prism nodal connectivity to that of polyhedron
 358    *  \param [inout] nodes - nodes of prism, return nodes of polyhedron
 359    *  \param [inout] elemType - return quantities of polyhedron,
 360    */
 361   //================================================================================
 362
 363   void prism2polyhedron( std::vector< const SMDS_MeshNode* > & nodes,
 364                          SMESH_MeshEditor::ElemFeatures &      elemType )
 365   {
 366     elemType.myPolyhedQuantities.clear();
 367     elemType.myNodes.clear();
 368
 369     // walls
 370     size_t nbBaseNodes = nodes.size() / 2;
 371     for ( size_t i = 0; i < nbBaseNodes; ++i )
 372     {
 373       int iNext = ( i + 1 ) % nbBaseNodes;
 374       elemType.myNodes.push_back( nodes[ iNext ]);
 375       elemType.myNodes.push_back( nodes[ i ]);
 376       elemType.myNodes.push_back( nodes[ i + nbBaseNodes ]);
 377       elemType.myNodes.push_back( nodes[ iNext + nbBaseNodes ]);
 378       elemType.myPolyhedQuantities.push_back( 4 );
 379     }
 380
 381     // base
 382     elemType.myNodes.insert( elemType.myNodes.end(), nodes.begin(), nodes.begin() + nbBaseNodes );
 383     elemType.myPolyhedQuantities.push_back( nbBaseNodes );
 384
 385     // top
 386     elemType.myNodes.insert( elemType.myNodes.end(), nodes.rbegin(), nodes.rbegin() + nbBaseNodes );
 387     elemType.myPolyhedQuantities.push_back( nbBaseNodes );
 388
 389     nodes.swap( elemType.myNodes );
 390   }
 391
 392   //================================================================================
 393   /*!
 394    * \brief Create prisms from faces
 395    *  \param [in] theEOS - shape to treat
 396    *  \param [inout] theMesh - offset mesh to fill in
 397    */
 398   //================================================================================
 399
 400   bool makePrisms( VLEdgesOnShape & theEOS,
 401                    SMESH_Mesh*      theMesh,
 402                    TNode2VLEdge   & theN2E )
 403   {
 404     SMESHDS_SubMesh* sm = theMesh->GetMeshDS()->MeshElements( theEOS._offsetShape );
 405     if ( !sm || sm->NbElements() == 0 )
 406       return true;
 407
 408     SMESH_MeshEditor editor( theMesh );
 409     SMESH_MeshEditor::ElemFeatures volumElem( SMDSAbs_Volume );
 410
 411     std::vector< const SMDS_MeshNode* > vNodes;
 412     std::vector< VLEdge*> edgesVec;
 413     for ( SMDS_ElemIteratorPtr eIt = sm->GetElements(); eIt->more(); )
 414     {
 415       const SMDS_MeshElement* face = eIt->next();
 416       if ( face->GetType() != SMDSAbs_Face )
 417         continue;
 418
 419       const int nbNodes = face->NbCornerNodes();
 420       edgesVec.resize( nbNodes );
 421       vNodes.resize( nbNodes * 2 );
 422       int maxNbLayer = 0, minNbLayer = IntegerLast();
 423       for ( int i = 0; i < nbNodes; ++i )
 424       {
 425         const SMDS_MeshNode* n = face->GetNode( i );
 426         if ( !theN2E.IsBound( n ))
 427           return false;
 428         edgesVec[ i ] = theN2E( n );
 429         maxNbLayer = Max( maxNbLayer, edgesVec[i]->_nodes.size() - 1 );
 430         minNbLayer = Min( minNbLayer, edgesVec[i]->_nodes.size() - 1 );
 431       }
 432       if ( maxNbLayer == minNbLayer )
 433       {
 434         size_t nbPrism = edgesVec[0]->_nodes.size() - 1;
 435         for ( size_t iP = 0; iP < nbPrism; ++iP )
 436         {
 437           vNodes.resize( nbNodes * 2 );
 438           for ( int i = 0; i < nbNodes; ++i )
 439           {
 440             vNodes[ i           ] = edgesVec[ i ]->_nodes[ iP + 1 ].Node();
 441             vNodes[ i + nbNodes ] = edgesVec[ i ]->_nodes[ iP ].Node();
 442           }
 443           volumElem.SetPoly( nbNodes > 4 );
 444           if ( volumElem.myIsPoly )
 445             prism2polyhedron( vNodes, volumElem );
 446
 447           if ( const SMDS_MeshElement* vol = editor.AddElement( vNodes, volumElem ))
 448             vol->setIsMarked( true ); // to add to group
 449         }
 450       }
 451       else // at inlet/outlet
 452       {
 453         makePolyhedron( edgesVec, vNodes, editor, volumElem );
 454       }
 455       editor.ClearLastCreated();
 456
 457       // move the face to the top of prisms, on mesh boundary
 458       //theMesh->GetMeshDS()->ChangeElementNodes( face, fNodes.data(), nbNodes );
 459     }
 460     return true;
 461   }
 462
 463   //================================================================================
 464   /*!
 465    * \brief Create faces from edges
 466    *  \param [inout] theEOS - shape to treat
 467    *  \param [inout] theMesh - offset mesh to fill in
 468    *  \param [inout] theN2E - map of node to VLEdge
 469    *  \param [inout] theFaceID - ID of WOVL FACE for new faces to set on
 470    *  \return bool - ok
 471    */
 472   //================================================================================
 473
 474   bool makeFaces( VLEdgesOnShape & theEOS,
 475                   SMESH_Mesh*      theMesh,
 476                   TNode2VLEdge   & theN2E,
 477                   const TGeomID    theFaceID)
 478   {
 479     SMESHDS_SubMesh* sm = theMesh->GetMeshDS()->MeshElements( theEOS._offsetShape );
 480     if ( !sm || sm->NbElements() == 0 )
 481       return true;
 482
 483     SMESH_MeshEditor editor( theMesh );
 484     SMESH_MeshEditor::ElemFeatures elem2D( SMDSAbs_Face );
 485
 486     std::vector< const SMDS_MeshNode* > fNodes( 4 );
 487     std::vector<const SMDS_MeshElement *> foundVolum;
 488     std::vector< VLEdge*> edgesVec;
 489     for ( SMDS_ElemIteratorPtr eIt = sm->GetElements(); eIt->more(); )
 490     {
 491       const SMDS_MeshElement* edge = eIt->next();
 492       if ( edge->GetType() != SMDSAbs_Edge )
 493         continue;
 494
 495       const int nbNodes = 2; //edge->NbCornerNodes();
 496       edgesVec.resize( nbNodes );
 497       fNodes.resize( nbNodes * 2 );
 498       for ( int i = 0; i < nbNodes; ++i )
 499       {
 500         const SMDS_MeshNode* n = edge->GetNode( i );
 501         if ( !theN2E.IsBound( n ))
 502           return false;
 503         edgesVec[ i ] = theN2E( n );
 504       }
 505       size_t nbFaces = edgesVec[0]->_nodes.size() - 1;
 506       for ( size_t iF = 0; iF < nbFaces; ++iF )
 507       {
 508         fNodes[ 0 ] = edgesVec[ 0 ]->_nodes[ iF ].Node();
 509         fNodes[ 1 ] = edgesVec[ 1 ]->_nodes[ iF ].Node();
 510         fNodes[ 2 ] = edgesVec[ 1 ]->_nodes[ iF + 1 ].Node();
 511         fNodes[ 3 ] = edgesVec[ 0 ]->_nodes[ iF + 1 ].Node();
 512
 513         if ( const SMDS_MeshElement* face = editor.AddElement( fNodes, elem2D ))
 514         {
 515           theMesh->GetMeshDS()->SetMeshElementOnShape( face, theFaceID );
 516
 517           if ( SMDS_Mesh::GetElementsByNodes( fNodes, foundVolum, SMDSAbs_Volume ))
 518           {
 519             TIDSortedElemSet faces = { face }, volumes = { foundVolum[0] };
 520             editor.Reorient2DBy3D( faces, volumes, /*outside=*/true );
 521           }
 522         }
 523       }
 524       editor.ClearLastCreated();
 525
 526       // move the face to the top of prisms, on mesh boundary
 527       //theMesh->GetMeshDS()->ChangeElementNodes( face, fNodes.data(), nbNodes );
 528     }
 529     return true;
 530   }
 531
 532   //================================================================================
 533   /*!
 534    * \brief Append the offset mesh to the initial one
 535    */
 536   //================================================================================
 537
 538   void copyMesh( SMESH_Mesh* theFromMesh,
 539                  SMESH_Mesh* theToMesh,
 540                  int         theSolidID)
 541   {
 542     SMESHDS_Mesh* offsetMDS = theFromMesh->GetMeshDS();
 543     SMESHDS_Mesh*   initMDS = theToMesh->GetMeshDS();
 544
 545     const smIdType nShift = initMDS->NbNodes();
 546
 547     for ( SMDS_NodeIteratorPtr nIt = offsetMDS->nodesIterator(); nIt->more(); )
 548     {
 549       SMESH_NodeXYZ pOff = nIt->next();
 550       const SMDS_MeshNode* n = initMDS->AddNodeWithID( pOff.X(), pOff.Y(), pOff.Z(),
 551                                                        nShift + pOff.Node()->GetID() );
 552       initMDS->SetNodeInVolume( n, theSolidID );
 553     }
 554
 555     SMESH_MeshEditor editor( theToMesh );
 556     SMESH_MeshEditor::ElemFeatures elemType;
 557     std::vector<smIdType> nIniIDs;
 558
 559     for ( SMDS_ElemIteratorPtr eIt = offsetMDS->elementsIterator(); eIt->more(); )
 560     {
 561       const SMDS_MeshElement* offElem = eIt->next();
 562       const int nbNodes = offElem->NbNodes();
 563       nIniIDs.resize( nbNodes );
 564       for ( int i = 0; i < nbNodes; ++i )
 565       {
 566         const SMDS_MeshNode* offNode = offElem->GetNode( i );
 567         nIniIDs[ i ] = offNode->GetID() + nShift;
 568       }
 569       elemType.Init( offElem, /*basicOnly=*/false );
 570       const SMDS_MeshElement* iniElem = editor.AddElement( nIniIDs, elemType );
 571       initMDS->SetMeshElementOnShape( iniElem, theSolidID );
 572       iniElem->setIsMarked( offElem->isMarked() );
 573       editor.ClearLastCreated();
 574     }
 575     return;
 576   }
 577
 578   //================================================================================
 579   /*!
 580    * \brief set elements of layers boundary to sub-meshes
 581    *  \param [in] theEOS - vlEdges of a sub-shape
 582    *  \param [inout] theMesh - the mesh
 583    *  \param [in] theN2E - map of node to vlEdge
 584    *  \param [in] theNShift - nb of nodes in the mesh before adding the offset mesh
 585    */
 586   //================================================================================
 587
 588   void setBnd2Sub( VLEdgesOnShape &     theEOS,
 589                    SMESH_Mesh*          theInitMesh,
 590                    SMESH_Mesh*          theOffsMesh,
 591                    const TNode2VLEdge & theN2E,
 592                    const smIdType       theNShift,
 593                    TIDSortedNodeSet&    theNodesToCheckCoinc )
 594   {
 595     SMESHDS_Mesh* offsetMDS = theOffsMesh->GetMeshDS();
 596     SMESHDS_Mesh*   initMDS = theInitMesh->GetMeshDS();
 597
 598     SMESHDS_SubMesh* sm = offsetMDS->MeshElements( theEOS._offsetShape );
 599     if ( !sm || ( sm->NbElements() + sm->NbNodes() == 0 ))
 600       return;
 601
 602     for ( SMDS_NodeIteratorPtr nIt = sm->GetNodes(); nIt->more(); )
 603     {
 604       const SMDS_MeshNode* offNode = nIt->next();
 605       VLEdge* vlEd = theN2E( offNode );
 606       const SMDS_MeshNode* nOffBnd = vlEd->_nodes.back().Node();
 607       const SMDS_MeshNode* nIniBnd = initMDS->FindNode( nOffBnd->GetID() + theNShift );
 608       initMDS->UnSetNodeOnShape( nIniBnd );
 609
 610       switch( theEOS._initShape.ShapeType() ) {
 611       case TopAbs_FACE:   initMDS->SetNodeOnFace( nIniBnd, TopoDS::Face( theEOS._initShape ),
 612                                                   vlEd->_uv[0], vlEd->_uv[1] );
 613         break;
 614       case TopAbs_EDGE:   initMDS->SetNodeOnEdge( nIniBnd, TopoDS::Edge( theEOS._initShape ),
 615                                                   vlEd->_uv[0]);
 616         break;
 617       case TopAbs_VERTEX: initMDS->SetNodeOnVertex( nIniBnd, TopoDS::Vertex( theEOS._initShape ));
 618         break;
 619       default:;
 620       }
 621     }
 622
 623     std::vector< const SMDS_MeshNode* > iniNodes, iniNodesBnd;
 624     std::vector< VLEdge*> edgesVec;
 625     for ( SMDS_ElemIteratorPtr eIt = sm->GetElements(); eIt->more(); )
 626     {
 627       const SMDS_MeshElement* offElem = eIt->next();
 628       if ( offElem->GetType() != SMDSAbs_Face &&
 629            offElem->GetType() != SMDSAbs_Edge )
 630         continue;
 631
 632       const int nbNodes = offElem->NbCornerNodes();
 633       iniNodes.resize( nbNodes );
 634       iniNodesBnd.resize( nbNodes );
 635       for ( int i = 0; i < nbNodes; ++i )
 636       {
 637         const SMDS_MeshNode* nOff = offElem->GetNode( i );
 638         const SMDS_MeshNode* nIni = initMDS->FindNode( nOff->GetID() + theNShift );
 639         iniNodes[ i ] = nIni;
 640         if ( !theN2E.IsBound( nOff ))
 641           break;
 642         VLEdge* vlEd = theN2E( nOff );
 643         const SMDS_MeshNode* nOffBnd = vlEd->_nodes.back().Node();
 644         const SMDS_MeshNode* nIniBnd = initMDS->FindNode( nOffBnd->GetID() + theNShift );
 645         iniNodesBnd[ i ] = nIniBnd;
 646       }
 647       if ( const SMDS_MeshElement* iniElem = initMDS->FindElement( iniNodes ))
 648       {
 649         initMDS->UnSetElementOnShape( iniElem );
 650         initMDS->SetMeshElementOnShape( iniElem, theEOS._initShape );
 651
 652         // move the face to the top of prisms, on init mesh boundary
 653         initMDS->ChangeElementNodes( iniElem, iniNodesBnd.data(), nbNodes );
 654
 655         if ( theEOS._toCheckCoinc )
 656           theNodesToCheckCoinc.insert( iniNodesBnd.begin(), iniNodesBnd.end() );
 657       }
 658     }
 659     return;
 660   }
 661
 662   //================================================================================
 663   /*!
 664    * \brief set elements of WOVL face to sub-meshes
 665    *  \param [in] theEOS - vlEdges of a sub-shape
 666    *  \param [inout] theMesh - the mesh
 667    *  \param [in] theN2E - map of node to vlEdge
 668    *  \param [in] theNShift - nb of nodes in the mesh before adding the offset mesh
 669    */
 670   //================================================================================
 671
 672   void setBnd2FVWL( VLEdgesOnShape &     theEOS,
 673                     SMESH_Mesh*          theInitMesh,
 674                     SMESH_Mesh*          theOffsMesh,
 675                     //const TNode2VLEdge & theN2E,
 676                     const smIdType       theNShift )
 677   {
 678     SMESHDS_Mesh* offsetMDS = theOffsMesh->GetMeshDS();
 679     SMESHDS_Mesh*   initMDS = theInitMesh->GetMeshDS();
 680
 681     SMESHDS_SubMesh* sm = offsetMDS->MeshElements( theEOS._offsetShape );
 682     if ( !sm || ( sm->NbElements() + sm->NbNodes() == 0 ))
 683       return;
 684
 685     for ( SMDS_NodeIteratorPtr nIt = sm->GetNodes(); nIt->more(); )
 686     {
 687       const SMDS_MeshNode* offNode = nIt->next();
 688       const SMDS_MeshNode* iniNode = initMDS->FindNode( offNode->GetID() + theNShift );
 689       initMDS->UnSetNodeOnShape( iniNode );
 690
 691       switch( theEOS._initShape.ShapeType() ) {
 692       case TopAbs_FACE:   initMDS->SetNodeOnFace( iniNode, TopoDS::Face( theEOS._initShape ));
 693         break;
 694       case TopAbs_EDGE:   initMDS->SetNodeOnEdge( iniNode, TopoDS::Edge( theEOS._initShape ));
 695         break;
 696       case TopAbs_VERTEX: initMDS->SetNodeOnVertex( iniNode, TopoDS::Vertex( theEOS._initShape ));
 697         break;
 698       default:;
 699       }
 700     }
 701
 702     std::vector< const SMDS_MeshNode* > iniNodes;
 703     for ( SMDS_ElemIteratorPtr eIt = sm->GetElements(); eIt->more(); )
 704     {
 705       const SMDS_MeshElement* offElem = eIt->next();
 706       // if ( offElem->GetType() != SMDSAbs_Face )
 707       //   continue;
 708
 709       int nbNodes = offElem->NbCornerNodes();
 710       iniNodes.resize( nbNodes );
 711       for ( int i = 0; i < nbNodes; ++i )
 712       {
 713         const SMDS_MeshNode* nOff = offElem->GetNode( i );
 714         const SMDS_MeshNode* nIni = initMDS->FindNode( nOff->GetID() + theNShift );
 715         iniNodes[ i ] = nIni;
 716       }
 717       if ( const SMDS_MeshElement* iniElem = initMDS->FindElement( iniNodes ))
 718       {
 719         initMDS->UnSetElementOnShape( iniElem );
 720         initMDS->SetMeshElementOnShape( iniElem, theEOS._initShapeID );
 721
 722         for ( const SMDS_MeshNode* n : iniNodes )
 723           if ( n->GetPosition()->GetDim() == 3 )
 724           {
 725             initMDS->UnSetNodeOnShape( n );
 726             if ( theEOS._initShape.ShapeType() == TopAbs_FACE )
 727               initMDS->SetNodeOnFace( n, theEOS._initShapeID );
 728             else
 729               initMDS->SetNodeOnEdge( n, theEOS._initShapeID );
 730           }
 731       }
 732     }
 733     return;
 734   }
 735 } // namespace
 736
 737
 738 //================================================================================
 739 /*!
 740  * \brief Create a temporary mesh used to hold elements of the offset shape
 741  */
 742 //================================================================================
 743
 744 SMESH_Mesh* StdMeshers_Cartesian_VL::ViscousBuilder::MakeOffsetMesh()
 745 {
 746   return _offsetMesh;
 747 }
 748
 749 //================================================================================
 750 /*!
 751  * \brief Constructor
 752  */
 753 //================================================================================
 754
 755 StdMeshers_Cartesian_VL::ViscousBuilder::ViscousBuilder( const StdMeshers_ViscousLayers* hyp,
 756                                                          const SMESH_Mesh &              mainMesh,
 757                                                          const TopoDS_Shape &            mainShape)
 758   : _hyp( hyp ), _offsetMesh( new TmpMesh )
 759 {
 760   // find shapes with viscous layers
 761   TopTools_IndexedMapOfShape faces;
 762   TopExp::MapShapes( mainShape, TopAbs_FACE, faces );
 763   const SMESHDS_Mesh* iniMDS = mainMesh.GetMeshDS();
 764
 765   if ( hyp->IsToIgnoreShapes() )
 766   {
 767     for ( int i = 1; i <= faces.Size(); ++i )
 768       _shapesWVL.insert( iniMDS->ShapeToIndex( faces( i )));
 769
 770     for ( TGeomID& sID : hyp->GetBndShapes() )
 771       _shapesWVL.erase( sID );
 772   }
 773   else
 774   {
 775     for ( TGeomID& sID : hyp->GetBndShapes() )
 776       _shapesWVL.insert( sID );
 777   }
 778
 779   for ( TGeomID fID : _shapesWVL )
 780   {
 781     const TopoDS_Shape & face = iniMDS->IndexToShape( fID );
 782     for ( TopExp_Explorer exp( face, TopAbs_EDGE ); exp.More(); exp.Next() )
 783       _shapesWVL.insert( iniMDS->ShapeToIndex( exp.Current() ));
 784     for ( TopExp_Explorer exp( face, TopAbs_VERTEX ); exp.More(); exp.Next() )
 785       _shapesWVL.insert( iniMDS->ShapeToIndex( exp.Current() ));
 786   }
 787
 788   // find EDGE and adjacent FACEs w/o VL
 789   for ( int i = 1; i <= faces.Size(); ++i )
 790   {
 791     const TopoDS_Shape& face = faces( i );
 792     TGeomID fID = iniMDS->ShapeToIndex( face );
 793     if ( _shapesWVL.count( fID ))
 794       continue;
 795     for ( TopExp_Explorer exp( face, TopAbs_EDGE ); exp.More(); exp.Next() )
 796       _edge2facesWOVL[ iniMDS->ShapeToIndex( exp.Current() )].push_back( fID );
 797   }
 798   return;
 799 }
 800
 801 //================================================================================
 802 /*!
 803  * \brief Destructor
 804  */
 805 //================================================================================
 806
 807 StdMeshers_Cartesian_VL::ViscousBuilder::~ViscousBuilder()
 808 {
 809   delete _offsetMesh; _offsetMesh = 0;
 810 }
 811
 812 //================================================================================
 813 /*!
 814  * \brief Create an offset shape from a given one
 815  *  \param [in] theShape - input shape
 816  *  \param [in] theMesh - main mesh
 817  *  \param [out] theError - error description
 818  *  \return TopoDS_Shape - result offset shape
 819  */
 820 //================================================================================
 821
 822 TopoDS_Shape
 823 StdMeshers_Cartesian_VL::ViscousBuilder::MakeOffsetShape(const TopoDS_Shape & theShape,
 824                                                          SMESH_Mesh &         theMesh,
 825                                                          std::string &        theError )
 826 {
 827   double offset = -_hyp->GetTotalThickness();
 828   double    tol = Precision::Confusion();
 829   _makeOffset.Initialize( theShape, offset, tol, BRepOffset_Skin, /*Intersection=*/false,
 830                           /*selfInter=*/false, GeomAbs_Intersection );
 831   // exclude inlet FACEs
 832   SMESHDS_Mesh* meshDS = theMesh.GetMeshDS();
 833   for ( TopExp_Explorer fEx( theShape, TopAbs_FACE ); fEx.More(); fEx.Next() )
 834   {
 835     TGeomID fID = meshDS->ShapeToIndex( fEx.Current() );
 836     if ( !_shapesWVL.count( fID ))
 837       _makeOffset.SetOffsetOnFace( TopoDS::Face( fEx.Current()), 0 );
 838   }
 839
 840   _makeOffset.MakeOffsetShape();
 841   if ( _makeOffset.IsDone() )
 842   {
 843     _offsetShape = _makeOffset.Shape();
 844     SMESH_MesherHelper::WriteShape( _offsetShape );////
 845
 846     _offsetMesh->ShapeToMesh( _offsetShape );
 847     _offsetMesh->GetSubMesh( _offsetShape )->DependsOn();
 848     return _offsetShape;
 849   }
 850
 851   switch ( _makeOffset.Error() )
 852   {
 853   case BRepOffset_NoError:
 854     theError = "OK. Offset performed successfully.";break;
 855   case BRepOffset_BadNormalsOnGeometry:
 856     theError = "Degenerated normal on input data.";break;
 857   case BRepOffset_C0Geometry:
 858     theError = "C0 continuity of input data.";break;
 859   case BRepOffset_NullOffset:
 860     theError = "Null offset of all faces.";break;
 861   case BRepOffset_NotConnectedShell:
 862     theError = "Incorrect set of faces to remove, the remaining shell is not connected.";break;
 863   case BRepOffset_CannotTrimEdges:
 864     theError = "Can not trim edges.";break;
 865   case BRepOffset_CannotFuseVertices:
 866     theError = "Can not fuse vertices.";break;
 867   case BRepOffset_CannotExtentEdge:
 868     theError = "Can not extent edge.";break;
 869   default:
 870     theError = "operation not done.";
 871   }
 872   theError = "BRepOffset_MakeOffset error: " + theError;
 873
 874   return TopoDS_Shape();
 875 }
 876
 877 //================================================================================
 878 /*!
 879  * \brief Return a sub-shape of the offset shape generated from a given initial sub-shape
 880  */
 881 //================================================================================
 882
 883 TopoDS_Shape StdMeshers_Cartesian_VL::ViscousBuilder::getOffsetSubShape( const TopoDS_Shape& S )
 884 {
 885   const TopTools_ListOfShape& newShapes = _makeOffset.Generated( S );
 886   for ( const TopoDS_Shape& ns : newShapes )
 887     if ( ns.ShapeType() == S.ShapeType() )
 888       return ns;
 889   return TopoDS_Shape();
 890 }
 891
 892 //================================================================================
 893 /*!
 894  * \brief Create prismatic mesh between _offsetShape and theShape
 895  *  \param [out] theMesh - mesh to fill in
 896  *  \param [in] theShape - initial shape
 897  *  \return bool - is Ok
 898  */
 899 //================================================================================
 900
 901 bool StdMeshers_Cartesian_VL::ViscousBuilder::MakeViscousLayers( SMESH_Mesh &         theMesh,
 902                                                                  const TopoDS_Shape & theShape )
 903 {
 904   SMESHDS_Mesh* offsetMDS = _offsetMesh->GetMeshDS();
 905   SMESHDS_Mesh*   initMDS = theMesh.GetMeshDS();
 906   offsetMDS->SetAllCellsNotMarked();
 907
 908   // Compute heights of viscous layers
 909   std::vector< double > vlH;
 910   computeVLHeight( _hyp, vlH );
 911
 912   std::vector< VLEdgesOnShape > edgesOnShape;
 913   edgesOnShape.reserve( offsetMDS->MaxShapeIndex() + 1 );
 914   TNode2VLEdge n2e;
 915
 916   // loop on sub-shapes to project nodes from offset boundary to initial boundary
 917   TopAbs_ShapeEnum types[3] = { TopAbs_VERTEX, TopAbs_EDGE, TopAbs_FACE };
 918   for ( TopAbs_ShapeEnum shType : types )
 919   {
 920     TopTools_IndexedMapOfShape shapes;
 921     TopExp::MapShapes( theShape, shType, shapes );
 922     for ( int i = 1; i <= shapes.Size(); ++i )
 923     {
 924       edgesOnShape.resize( edgesOnShape.size() + 1 );
 925       VLEdgesOnShape& EOS = edgesOnShape.back();
 926
 927       EOS._initShape    = shapes( i );
 928       EOS._offsetShape  = getOffsetSubShape( EOS._initShape );
 929       EOS._initShapeID  = initMDS->ShapeToIndex( EOS._initShape );
 930       EOS._hasVL        = _shapesWVL.count( EOS._initShapeID );
 931       EOS._toCheckCoinc = false;
 932       if ( !EOS._hasVL )
 933         continue;
 934
 935       // project boundary nodes of offset mesh to boundary of init mesh
 936       // (new nodes are created in the offset mesh)
 937       switch( EOS._offsetShape.ShapeType() ) {
 938       case TopAbs_VERTEX:
 939       {
 940         EOS._edges.resize( 1 );
 941         EOS._edges[0]._nodes.resize( 2 );
 942         EOS._edges[0]._nodes[0] = SMESH_Algo::VertexNode( TopoDS::Vertex( EOS._offsetShape ),
 943                                                           offsetMDS );
 944         gp_Pnt offP = BRep_Tool::Pnt( TopoDS::Vertex( EOS._initShape ));
 945         EOS._edges[0]._nodes[1] = offsetMDS->AddNode( offP.X(), offP.Y(), offP.Z() );
 946         //EOS._edges[0]._length   = offP.Distance( EOS._edges[0]._nodes[0] );
 947         n2e.Bind( EOS._edges[0]._nodes[0].Node(), & EOS._edges[0] );
 948         break;
 949       }
 950       case TopAbs_EDGE:
 951       {
 952         projectToEdge( EOS, offsetMDS, n2e );
 953         break;
 954       }
 955       case TopAbs_FACE:
 956       {
 957         projectToFace( EOS, offsetMDS, n2e );
 958         break;
 959       }
 960       default:;
 961       }
 962
 963       // create nodes of layers
 964       if ( _hyp->GetNumberLayers() > 1 )
 965       {
 966         //if ( _shapesWVL.count( EOS._initShapeID ))
 967         for ( size_t i = 0; i < EOS._edges.size(); ++i )
 968         {
 969           divideVLEdge( &EOS._edges[ i ], vlH, offsetMDS );
 970         }
 971       }
 972     } // loop on shapes
 973   } // loop on shape types
 974
 975   // create prisms
 976   bool prismsOk = true;
 977   for ( size_t i = 0; i < edgesOnShape.size(); ++i )
 978   {
 979     VLEdgesOnShape& EOS = edgesOnShape[ i ];
 980     if ( EOS._initShape.ShapeType() == TopAbs_FACE && EOS._hasVL )
 981     {
 982       if ( !makePrisms( EOS, _offsetMesh, n2e ))
 983         prismsOk = false;
 984     }
 985   }
 986   if ( prismsOk )
 987   {
 988     // create faces on FACEs WOVL
 989     for ( size_t i = 0; i < edgesOnShape.size(); ++i )
 990     {
 991       VLEdgesOnShape& EOS = edgesOnShape[ i ];
 992       if ( EOS._initShape.ShapeType() == TopAbs_EDGE && EOS._hasVL )
 993       {
 994         auto e2f = _edge2facesWOVL.find( EOS._initShapeID );
 995         if ( e2f != _edge2facesWOVL.end() && !e2f->second.empty() )
 996         {
 997           TopoDS_Shape f = initMDS->IndexToShape( e2f->second[0] );
 998           TopoDS_Shape f2 = getOffsetSubShape( f );
 999           //cout << e2f->second[0] << " OFF " << offsetMDS->ShapeToIndex( f2 ) << endl;
1000           makeFaces( EOS, _offsetMesh, n2e, offsetMDS->ShapeToIndex( f2 ) );
1001         }
1002       }
1003     }
1004   }
1005
1006   // copy offset mesh to the main one
1007   initMDS->Modified();
1008   initMDS->CompactMesh();
1009   smIdType nShift = initMDS->NbNodes();
1010   TGeomID solidID = initMDS->ShapeToIndex( theShape );
1011   copyMesh( _offsetMesh, & theMesh, solidID );
1012
1013
1014   if ( !prismsOk )
1015   {
1016     if ( SMESH_subMesh * sm = theMesh.GetSubMesh( theShape ))
1017     {
1018       sm->GetComputeError() =
1019         SMESH_ComputeError::New( COMPERR_ALGO_FAILED,
1020                                  "Viscous layers construction error: bad mesh on offset geometry" );
1021     }
1022     return prismsOk;
1023   }
1024
1025   // set elements of layers boundary to sub-meshes
1026   TIDSortedNodeSet nodesToCheckCoinc;
1027   for ( size_t i = 0; i < edgesOnShape.size(); ++i )
1028   {
1029     VLEdgesOnShape& EOS = edgesOnShape[ i ];
1030     if ( EOS._hasVL )
1031       setBnd2Sub( EOS, &theMesh, _offsetMesh, n2e, nShift, nodesToCheckCoinc );
1032     else
1033       setBnd2FVWL( EOS, &theMesh, _offsetMesh, nShift );
1034   }
1035
1036   // merge coincident nodes
1037   SMESH_MeshEditor editor( &theMesh );
1038   SMESH_MeshEditor::TListOfListOfNodes nodesToMerge;
1039   editor.FindCoincidentNodes( nodesToCheckCoinc, vlH[0]/50., nodesToMerge, false );
1040   editor.MergeNodes( nodesToMerge );
1041
1042   // create a group
1043   if ( !_hyp->GetGroupName().empty() )
1044   {
1045     SMDS_MeshGroup* group = _hyp->CreateGroup( _hyp->GetGroupName(), theMesh, SMDSAbs_Volume );
1046
1047     for ( SMDS_ElemIteratorPtr it = initMDS->elementsIterator(); it->more(); )
1048     {
1049       const SMDS_MeshElement* el = it->next();
1050       if ( el->isMarked() )
1051         group->Add( el );
1052     }
1053   }
1054
1055   return prismsOk;
1056 }