src/StdMeshers/StdMeshers_Cartesian_VL.cxx

   1 // Copyright (C) 2007-2023  CEA, EDF, 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_MeshAlgos.hxx>
  34 #include <SMESH_Mesh.hxx>
  35 #include <SMESH_MeshEditor.hxx>
  36 #include <SMESH_MesherHelper.hxx>
  37 #include <SMESH_TypeDefs.hxx>
  38 #include <SMESH_subMesh.hxx>
  39
  40 #include <BRepAdaptor_Curve.hxx>
  41 #include <BRep_Builder.hxx>
  42 #include <BRepBuilderAPI_MakeFace.hxx>
  43 #include <BRepGProp.hxx>
  44 #include <BRepTopAdaptor_FClass2d.hxx>
  45 #include <BRep_Tool.hxx>
  46 #include <GProp_GProps.hxx>
  47 #include <ShapeAnalysis_Curve.hxx>
  48 #include <ShapeAnalysis_Surface.hxx>
  49 #include <TopExp.hxx>
  50 #include <TopExp_Explorer.hxx>
  51 #include <TopoDS.hxx>
  52
  53 using namespace StdMeshers_Cartesian_VL;
  54
  55 namespace
  56 {
  57   typedef int                     TGeomID; // IDs of sub-shapes
  58
  59   /*!
  60    * \brief Temporary mesh
  61    */
  62   struct TmpMesh: public SMESH_Mesh
  63   {
  64     TmpMesh() {
  65       _isShapeToMesh = (_id = 0);
  66       _meshDS  = new SMESHDS_Mesh( _id, true );
  67     }
  68   };
  69   // --------------------------------------------------------------------------
  70   /*!
  71    * \brief Edge of viscous layers; goes from a grid node by normal to boundary
  72    */
  73   struct VLEdge
  74   {
  75     std::vector< SMESH_NodeXYZ > _nodes;
  76     double                       _uv[2]; // position of TgtNode() on geometry
  77     double                       _length;
  78
  79     const SMDS_MeshNode* TgtNode() const { return _nodes.back().Node(); }
  80   };
  81   typedef NCollection_DataMap< const SMDS_MeshNode*, VLEdge* > TNode2VLEdge;
  82   // --------------------------------------------------------------------------
  83   /*!
  84    * \brief VLEdge's of one shape
  85    */
  86   struct VLEdgesOnShape
  87   {
  88     std::vector< VLEdge > _edges;
  89     TopoDS_Shape          _initShape;
  90     TopoDS_Shape          _offsetShape;
  91     int                   _initShapeID;
  92     bool                  _hasVL;
  93     bool                  _toCheckCoinc; // to check if nodes on shape boundary coincide
  94   };
  95
  96   //================================================================================
  97   /*!
  98    * \brief Project a point on offset FACE to EDGEs of an initial FACE
  99    *  \param [in] offP - point to project
 100    *  \param [in] initFace - FACE to project on
 101    *  \return gp_Pnt - projection point
 102    */
 103   //================================================================================
 104
 105   gp_Pnt projectToWire( const gp_Pnt&      offP,
 106                         const TopoDS_Face& initFace,
 107                         gp_Pnt2d &         uv )
 108   {
 109     double   minDist = Precision::Infinite();
 110     const double tol = Precision::Confusion();
 111     gp_Pnt              proj, projction;
 112     Standard_Real       param;
 113     ShapeAnalysis_Curve projector;
 114
 115     for ( TopExp_Explorer eExp( initFace, TopAbs_EDGE ); eExp.More(); eExp.Next() )
 116     {
 117       BRepAdaptor_Curve initCurve( TopoDS::Edge( eExp.Current() ));
 118       //const double f = initCurve.FirstParameter(), l = initCurve.LastParameter();
 119       double dist = projector.Project( initCurve, offP, tol, proj, param, /*adjustToEnds=*/false );
 120       if ( dist < minDist )
 121       {
 122         projction = proj;
 123         minDist = dist;
 124       }
 125     }
 126     uv.SetCoord(0.,0.); // !!!!!!!
 127     return projction;
 128   }
 129
 130   //================================================================================
 131   /*!
 132    * \brief Project nodes from offset FACE to initial FACE
 133    *  \param [inout] theEOS - VL edges on a geom FACE
 134    *  \param [inout] theOffsetMDS - offset mesh to fill in
 135    */
 136   //================================================================================
 137
 138   void projectToFace( VLEdgesOnShape & theEOS,
 139                       SMESHDS_Mesh*    theOffsetMDS,
 140                       TNode2VLEdge &   theN2E )
 141   {
 142     SMESHDS_SubMesh* sm = theOffsetMDS->MeshElements( theEOS._offsetShape );
 143     if ( !sm || sm->NbElements() == 0 || sm->NbNodes() == 0 )
 144       return;
 145
 146     theEOS._edges.resize( sm->NbNodes() );
 147
 148     const TopoDS_Face& initFace = TopoDS::Face( theEOS._initShape );
 149     ShapeAnalysis_Surface projector( BRep_Tool::Surface( initFace ));
 150     const double clsfTol = 1e2 * BRep_Tool::MaxTolerance( initFace, TopAbs_VERTEX );
 151     BRepTopAdaptor_FClass2d classifier( initFace, clsfTol );
 152
 153     const double        tol = Precision::Confusion();
 154     //const double f = initCurve.FirstParameter(), l = initCurve.LastParameter();
 155     gp_Pnt              proj;
 156
 157     int iN = 0;
 158     for ( SMDS_NodeIteratorPtr nIt = sm->GetNodes(); nIt->more(); ++iN )
 159     {
 160       SMESH_NodeXYZ offP = nIt->next();
 161       gp_Pnt2d uv = projector.ValueOfUV( offP, tol );
 162       TopAbs_State st = classifier.Perform( uv );
 163       if ( st == TopAbs_IN || st == TopAbs_ON )
 164       {
 165         proj = projector.Value( uv );
 166       }
 167       else
 168       {
 169         proj = projectToWire( offP, initFace, uv );
 170       }
 171       if ( st == TopAbs_ON || st == TopAbs_OUT )
 172         theEOS._toCheckCoinc = true;
 173
 174       VLEdge & vlEdge = theEOS._edges[ iN ];
 175       vlEdge._nodes.push_back( offP.Node() );
 176       vlEdge._nodes.push_back( theOffsetMDS->AddNode( proj.X(), proj.Y(), proj.Z() ));
 177       vlEdge._uv[0] = uv.X();
 178       vlEdge._uv[1] = uv.Y();
 179       //vlEdge._length = proj.Distance( offP );
 180
 181       theN2E.Bind( offP.Node(), &vlEdge );
 182     }
 183     return;
 184
 185   }
 186
 187   //================================================================================
 188   /*!
 189    * \brief Project nodes from offset EDGE to initial EDGE
 190    *  \param [inout] theEOS - VL edges on a geom EDGE
 191    *  \param [inout] theOffsetMDS - offset mesh to add new nodes to
 192    */
 193   //================================================================================
 194
 195   void projectToEdge( VLEdgesOnShape & theEOS,
 196                       SMESHDS_Mesh*    theOffsetMDS,
 197                       TNode2VLEdge &   theN2E,
 198                       bool createVertex )
 199   {
 200     SMESHDS_SubMesh* sm = theOffsetMDS->MeshElements( theEOS._offsetShape );
 201
 202     if ( !sm || sm->NbElements() == 0 )
 203       return;
 204
 205     int addVertexNode = createVertex ? 1 : 0;
 206     theEOS._edges.resize( sm->NbNodes() + addVertexNode ); // +1 to set the vertex
 207
 208     ShapeAnalysis_Curve projector;
 209     BRepAdaptor_Curve   initCurve( TopoDS::Edge( theEOS._initShape ));
 210     const double        tol = Precision::Confusion();
 211     const double f = initCurve.FirstParameter(), l = initCurve.LastParameter();
 212     gp_Pnt              proj;
 213     Standard_Real       param;
 214
 215     int iN = 0;
 216     for ( SMDS_NodeIteratorPtr nIt = sm->GetNodes(); nIt->more(); ++iN )
 217     {
 218       SMESH_NodeXYZ offP = nIt->next();
 219       double dist = projector.Project( initCurve, offP, tol, proj, param, /*adjustToEnds=*/false );
 220       bool paramOK = ( f < param && param < l );
 221       if ( !paramOK )
 222       {
 223         if ( param < f )
 224           param = f;
 225         else if ( param > l )
 226           param = l;
 227         theEOS._toCheckCoinc = true;
 228       }
 229       proj = initCurve.Value( param );
 230
 231       VLEdge & vlEdge = theEOS._edges[ iN ];
 232       vlEdge._nodes.push_back( offP.Node() );
 233       vlEdge._nodes.push_back( theOffsetMDS->AddNode( proj.X(), proj.Y(), proj.Z() ));
 234       vlEdge._uv[0] = param;
 235       vlEdge._length = paramOK ? dist : proj.Distance( offP );
 236
 237       theN2E.Bind( offP.Node(), &vlEdge );
 238     }
 239
 240     if ( createVertex )
 241     {
 242       // It is possible to define the vertex projections from the existing edges
 243       // EOS._offsetShape the edge generated from the original edge
 244       // Get the first vertex of both edges to define the connecting edges
 245       auto offsetEdge = TopoDS::Edge( theEOS._offsetShape );
 246       auto initEdge   = TopoDS::Edge( theEOS._initShape );
 247       TopoDS_Vertex offsetVertex;
 248       TopoDS_Vertex initVertex;
 249
 250       if ( offsetEdge.Orientation() == TopAbs_FORWARD )
 251       {
 252         offsetVertex  = TopExp::FirstVertex ( offsetEdge );
 253         initVertex    = TopExp::FirstVertex ( initEdge );
 254       }
 255       else
 256       {
 257         offsetVertex  = TopExp::LastVertex ( offsetEdge );
 258         initVertex    = TopExp::LastVertex ( initEdge );
 259       }
 260
 261       VLEdge & vlEdge = theEOS._edges[ iN ];
 262       vlEdge._nodes.resize( 2 );
 263       vlEdge._nodes[0] = SMESH_Algo::VertexNode( offsetVertex, theOffsetMDS );
 264
 265       gp_Pnt offP = BRep_Tool::Pnt( initVertex );
 266
 267       vlEdge._nodes[1] = theOffsetMDS->AddNode( offP.X(), offP.Y(), offP.Z() );
 268       theN2E.Bind( vlEdge._nodes[0].Node(), &vlEdge );
 269     }
 270     return;
 271   }
 272
 273   //================================================================================
 274   /*!
 275    * \brief Compute heights of viscous layers and finds FACEs with VL
 276    *  \param [in] hyp - viscous layers hypothesis
 277    *  \param [in] mesh - the main mesh
 278    *  \param [in] shape - the main shape
 279    *  \param [out] vlH - heights of viscous layers to compute
 280    *  \param [out] shapesWVL - IDs of shapes with VL
 281    *  \return bool - isOK
 282    */
 283   //================================================================================
 284
 285   void computeVLHeight( const StdMeshers_ViscousLayers* hyp,
 286                         std::vector< double > &         vlH )
 287   {
 288     const double  T = hyp->GetTotalThickness();
 289     const double  f = hyp->GetStretchFactor();
 290     const int     N = hyp->GetNumberLayers();
 291     const double h0 = hyp->Get1stLayerThickness( T, f, N );
 292
 293     vlH.reserve( hyp->GetNumberLayers() );
 294     vlH.push_back( h0 );
 295     for ( double h = h0 * f; (int)vlH.size() < N-1; h *= f )
 296       vlH.push_back( h + vlH.back() );
 297     vlH.push_back( T );
 298   }
 299
 300   //================================================================================
 301   /*!
 302    * \brief Create intermediate nodes on VLEdge
 303    *  \param [inout] vlEdge - VLEdge to divide
 304    *  \param [in] vlH - thicknesses of layers
 305    *  \param [inout] mesh - the mesh no fill in
 306    */
 307   //================================================================================
 308
 309   void divideVLEdge( VLEdge*                      vlEdge,
 310                      const std::vector< double >& vlH,
 311                      SMESHDS_Mesh*                mesh )
 312   {
 313     SMESH_NodeXYZ lastNode = vlEdge->_nodes.back();
 314     SMESH_NodeXYZ frstNode = vlEdge->_nodes[0];
 315     gp_XYZ dir = frstNode - lastNode;
 316
 317     vlEdge->_nodes.resize( vlH.size() + 1 );
 318
 319     for ( size_t i = 1; i < vlH.size(); ++i )
 320     {
 321       double r = vlH[ i-1 ] / vlH.back();
 322       gp_XYZ p = lastNode + dir * r;
 323       vlEdge->_nodes[ vlH.size() - i ] = mesh->AddNode( p.X(), p.Y(), p.Z() );
 324     }
 325     vlEdge->_nodes.back() = lastNode;
 326   }
 327
 328   //================================================================================
 329   /*!
 330    * \brief Create a polyhedron from nodes of VLEdge's
 331    *  \param [in] edgesVec - the VLEdge's
 332    *  \param [in] vNodes - node buffer
 333    *  \param [in] editor - editor of offset mesh
 334    */
 335   //================================================================================
 336
 337   bool makePolyhedron( const std::vector< VLEdge*> &         edgesVec,
 338                        std::vector< const SMDS_MeshNode* > & vNodes,
 339                        SMESH_MeshEditor&                     editor,
 340                        SMESH_MeshEditor::ElemFeatures        elemType)
 341   {
 342     elemType.SetPoly( true );
 343     elemType.myPolyhedQuantities.clear();
 344     elemType.myNodes.clear();
 345
 346     // add facets of walls
 347     size_t nbBaseNodes = edgesVec.size();
 348     for ( size_t iN = 0; iN < nbBaseNodes; ++iN )
 349     {
 350       VLEdge* e0 = edgesVec[ iN ];
 351       VLEdge* e1 = edgesVec[( iN + 1 ) % nbBaseNodes ];
 352       size_t nbN0 = e0->_nodes.size();
 353       size_t nbN1 = e1->_nodes.size();
 354       if ( nbN0 == nbN1 )
 355       {
 356         for ( size_t i = 1; i < nbN0; ++i )
 357         {
 358           elemType.myNodes.push_back( e1->_nodes[ i - 1 ].Node());
 359           elemType.myNodes.push_back( e1->_nodes[ i     ].Node());
 360           elemType.myNodes.push_back( e0->_nodes[ i     ].Node());
 361           elemType.myNodes.push_back( e0->_nodes[ i - 1 ].Node());
 362           elemType.myPolyhedQuantities.push_back( 4 );
 363         }
 364       }
 365       else
 366       {
 367         for ( size_t i = 0; i < nbN1; ++i )
 368           elemType.myNodes.push_back( e1->_nodes[ i ].Node() );
 369         for ( size_t i = 0; i < nbN0; ++i )
 370           elemType.myNodes.push_back( e0->_nodes[ nbN0 - 1 - i ].Node() );
 371         elemType.myPolyhedQuantities.push_back( nbN0 + nbN1 );
 372       }
 373     }
 374
 375     // add facets of top
 376     vNodes.clear();
 377     for ( size_t iN = 0; iN < nbBaseNodes; ++iN )
 378     {
 379       VLEdge* e0 = edgesVec[ iN ];
 380       elemType.myNodes.push_back( e0->_nodes.back().Node() );
 381       vNodes.push_back( e0->_nodes[ 0 ].Node());
 382     }
 383     elemType.myPolyhedQuantities.push_back( nbBaseNodes );
 384
 385     // add facets of bottom
 386     elemType.myNodes.insert( elemType.myNodes.end(), vNodes.rbegin(), vNodes.rend() );
 387     elemType.myPolyhedQuantities.push_back( nbBaseNodes );
 388
 389     const SMDS_MeshElement* vol = editor.AddElement( elemType.myNodes, elemType );
 390     vol->setIsMarked( true ); // to add to group
 391
 392     return vol;
 393   }
 394
 395   //================================================================================
 396   /*!
 397    * \brief Transform prism nodal connectivity to that of polyhedron
 398    *  \param [inout] nodes - nodes of prism, return nodes of polyhedron
 399    *  \param [inout] elemType - return quantities of polyhedron,
 400    */
 401   //================================================================================
 402
 403   void prism2polyhedron( std::vector< const SMDS_MeshNode* > & nodes,
 404                          SMESH_MeshEditor::ElemFeatures &      elemType )
 405   {
 406     elemType.myPolyhedQuantities.clear();
 407     elemType.myNodes.clear();
 408
 409     // walls
 410     size_t nbBaseNodes = nodes.size() / 2;
 411     for ( size_t i = 0; i < nbBaseNodes; ++i )
 412     {
 413       int iNext = ( i + 1 ) % nbBaseNodes;
 414       elemType.myNodes.push_back( nodes[ iNext ]);
 415       elemType.myNodes.push_back( nodes[ i ]);
 416       elemType.myNodes.push_back( nodes[ i + nbBaseNodes ]);
 417       elemType.myNodes.push_back( nodes[ iNext + nbBaseNodes ]);
 418       elemType.myPolyhedQuantities.push_back( 4 );
 419     }
 420
 421     // base
 422     elemType.myNodes.insert( elemType.myNodes.end(), nodes.begin(), nodes.begin() + nbBaseNodes );
 423     elemType.myPolyhedQuantities.push_back( nbBaseNodes );
 424
 425     // top
 426     elemType.myNodes.insert( elemType.myNodes.end(), nodes.rbegin(), nodes.rbegin() + nbBaseNodes );
 427     elemType.myPolyhedQuantities.push_back( nbBaseNodes );
 428
 429     nodes.swap( elemType.myNodes );
 430   }
 431
 432   //================================================================================
 433   /*!
 434    * \brief Create prisms from faces
 435    *  \param [in] theEOS - shape to treat
 436    *  \param [inout] theMesh - offset mesh to fill in
 437    */
 438   //================================================================================
 439
 440   bool makePrisms( VLEdgesOnShape & theEOS,
 441                    SMESH_Mesh*      theMesh,
 442                    TNode2VLEdge   & theN2E )
 443   {
 444     SMESHDS_SubMesh* sm = theMesh->GetMeshDS()->MeshElements( theEOS._offsetShape );
 445     if ( !sm || sm->NbElements() == 0 )
 446       return true;
 447
 448     SMESH_MeshEditor editor( theMesh );
 449     SMESH_MeshEditor::ElemFeatures volumElem( SMDSAbs_Volume );
 450
 451     std::vector< const SMDS_MeshNode* > vNodes;
 452     std::vector< VLEdge*> edgesVec;
 453     for ( SMDS_ElemIteratorPtr eIt = sm->GetElements(); eIt->more(); )
 454     {
 455       const SMDS_MeshElement* face = eIt->next();
 456       if ( face->GetType() != SMDSAbs_Face )
 457         continue;
 458
 459       const int nbNodes = face->NbCornerNodes();
 460       edgesVec.resize( nbNodes );
 461       vNodes.resize( nbNodes * 2 );
 462       int maxNbLayer = 0, minNbLayer = IntegerLast();
 463       for ( int i = 0; i < nbNodes; ++i )
 464       {
 465         const SMDS_MeshNode* n = face->GetNode( i );
 466         if ( !theN2E.IsBound( n ))
 467           return false;
 468         edgesVec[ i ] = theN2E( n );
 469         maxNbLayer = Max( maxNbLayer, edgesVec[i]->_nodes.size() - 1 );
 470         minNbLayer = Min( minNbLayer, edgesVec[i]->_nodes.size() - 1 );
 471       }
 472       if ( maxNbLayer == minNbLayer )
 473       {
 474         size_t nbPrism = edgesVec[0]->_nodes.size() - 1;
 475         for ( size_t iP = 0; iP < nbPrism; ++iP )
 476         {
 477           vNodes.resize( nbNodes * 2 );
 478           for ( int i = 0; i < nbNodes; ++i )
 479           {
 480             vNodes[ i           ] = edgesVec[ i ]->_nodes[ iP + 1 ].Node();
 481             vNodes[ i + nbNodes ] = edgesVec[ i ]->_nodes[ iP ].Node();
 482           }
 483           volumElem.SetPoly( nbNodes > 4 );
 484           if ( volumElem.myIsPoly )
 485             prism2polyhedron( vNodes, volumElem );
 486
 487           if ( const SMDS_MeshElement* vol = editor.AddElement( vNodes, volumElem ))
 488             vol->setIsMarked( true ); // to add to group
 489         }
 490       }
 491       else // at inlet/outlet
 492         makePolyhedron( edgesVec, vNodes, editor, volumElem );
 493
 494       editor.ClearLastCreated();
 495
 496       // move the face to the top of prisms, on mesh boundary
 497       //theMesh->GetMeshDS()->ChangeElementNodes( face, fNodes.data(), nbNodes );
 498     }
 499     return true;
 500   }
 501
 502   //================================================================================
 503   /*!
 504    * \brief Create faces from edges
 505    *  \param [inout] theEOS - shape to treat
 506    *  \param [inout] theMesh - offset mesh to fill in
 507    *  \param [inout] theN2E - map of node to VLEdge
 508    *  \param [inout] theFaceID - ID of WOVL FACE for new faces to set on
 509    *  \param [in]    isMainShape2D - used to identify the geometry where the new elements are included
 510    *  \return bool - ok
 511    */
 512   //================================================================================
 513
 514   bool makeFaces( VLEdgesOnShape & theEOS,
 515                   SMESH_Mesh*      theMesh,
 516                   TNode2VLEdge   & theN2E,
 517                   const TGeomID    theFaceID,
 518                   bool isMainShape2D = false )
 519   {
 520     SMESHDS_SubMesh* sm = theMesh->GetMeshDS()->MeshElements( theEOS._offsetShape );
 521     if ( !sm || sm->NbElements() == 0 )
 522       return true;
 523
 524     SMESH_MeshEditor editor( theMesh );
 525     SMESH_MeshEditor::ElemFeatures elem2D( SMDSAbs_Face );
 526
 527     std::vector< const SMDS_MeshNode* > fNodes( 4 );
 528     std::vector<const SMDS_MeshElement *> foundVolum;
 529     std::vector< VLEdge*> edgesVec;
 530     TIDSortedElemSet  refSetFace;
 531
 532     // Check orientation of face and re
 533     gp_XYZ refNormalVector(0.0,0.0,0.0);
 534     if ( isMainShape2D )
 535     {
 536       SMESHDS_Mesh* offsetMDS = theMesh->GetMeshDS();
 537       for ( SMDS_ElemIteratorPtr eIt = offsetMDS->elementsIterator(); eIt->more(); )
 538       {
 539         const SMDS_MeshElement* refFace = eIt->next(); // define the ref
 540         if ( refFace->GetType() == SMDSAbs_Face )
 541         {
 542           SMESH_MeshAlgos::FaceNormal( refFace, refNormalVector, /*normalized=*/true );
 543           break;
 544         }
 545       }
 546       if ( refNormalVector.X() == 0.0 && refNormalVector.Y() == 0.0 && refNormalVector.Z() == 0.0 )
 547         throw SALOME_Exception("No 2D element found in the mesh!\n");
 548
 549     }
 550
 551     for ( SMDS_ElemIteratorPtr eIt = sm->GetElements(); eIt->more(); )
 552     {
 553       const SMDS_MeshElement* edge = eIt->next();
 554       if ( edge->GetType() != SMDSAbs_Edge )
 555         continue;
 556
 557       const int nbNodes = 2; //edge->NbCornerNodes();
 558       edgesVec.resize( nbNodes );
 559       fNodes.resize( nbNodes * 2 );
 560       for ( int i = 0; i < nbNodes; ++i )
 561       {
 562         const SMDS_MeshNode* n = edge->GetNode( i );
 563         if ( !theN2E.IsBound( n ))
 564           return false;
 565         edgesVec[ i ] = theN2E( n );
 566       }
 567       size_t nbFaces = edgesVec[0]->_nodes.size() - 1;
 568
 569       for ( size_t iF = 0; iF < nbFaces; ++iF )
 570       {
 571         fNodes[ 0 ] = edgesVec[ 0 ]->_nodes[ iF ].Node();
 572         fNodes[ 1 ] = edgesVec[ 1 ]->_nodes[ iF ].Node();
 573         fNodes[ 2 ] = edgesVec[ 1 ]->_nodes[ iF + 1 ].Node();
 574         fNodes[ 3 ] = edgesVec[ 0 ]->_nodes[ iF + 1 ].Node();
 575
 576         if ( const SMDS_MeshElement* face = editor.AddElement( fNodes, elem2D ))
 577         {
 578           theMesh->GetMeshDS()->SetMeshElementOnShape( face, theFaceID );
 579
 580           if ( SMDS_Mesh::GetElementsByNodes( fNodes, foundVolum, SMDSAbs_Volume ))
 581           {
 582             TIDSortedElemSet faces = { face }, volumes = { foundVolum[0] };
 583             editor.Reorient2DBy3D( faces, volumes, /*outside=*/true );
 584           }
 585           else if ( isMainShape2D )
 586           {
 587             gp_XYZ elementNormal;
 588             SMESH_MeshAlgos::FaceNormal( face, elementNormal, /*normalized=*/true );
 589             if ( elementNormal * refNormalVector < 0.0 /* diff orientation from the ref element */)
 590               editor.Reorient( face );
 591           }
 592         }
 593       }
 594       editor.ClearLastCreated();
 595
 596       // move the face to the top of prisms, on mesh boundary
 597       //theMesh->GetMeshDS()->ChangeElementNodes( face, fNodes.data(), nbNodes );
 598     }
 599     return true;
 600   }
 601
 602   //================================================================================
 603   /*!
 604    * \brief Append the offset mesh to the initial one
 605    */
 606   //================================================================================
 607
 608   void copyMesh( SMESH_Mesh* theFromMesh,
 609                  SMESH_Mesh* theToMesh,
 610                  int         theSolidID)
 611   {
 612     SMESHDS_Mesh* offsetMDS = theFromMesh->GetMeshDS();
 613     SMESHDS_Mesh*   initMDS = theToMesh->GetMeshDS();
 614
 615     const smIdType nShift = initMDS->NbNodes();
 616
 617     for ( SMDS_NodeIteratorPtr nIt = offsetMDS->nodesIterator(); nIt->more(); )
 618     {
 619       SMESH_NodeXYZ pOff = nIt->next();
 620       const SMDS_MeshNode* n = initMDS->AddNodeWithID( pOff.X(), pOff.Y(), pOff.Z(),
 621                                                        nShift + pOff.Node()->GetID() );
 622       initMDS->SetNodeInVolume( n, theSolidID );
 623     }
 624
 625     SMESH_MeshEditor editor( theToMesh );
 626     SMESH_MeshEditor::ElemFeatures elemType;
 627     std::vector<smIdType> nIniIDs;
 628
 629     for ( SMDS_ElemIteratorPtr eIt = offsetMDS->elementsIterator(); eIt->more(); )
 630     {
 631       const SMDS_MeshElement* offElem = eIt->next();
 632       const int nbNodes = offElem->NbNodes();
 633       nIniIDs.resize( nbNodes );
 634       for ( int i = 0; i < nbNodes; ++i )
 635       {
 636         const SMDS_MeshNode* offNode = offElem->GetNode( i );
 637         nIniIDs[ i ] = offNode->GetID() + nShift;
 638       }
 639       elemType.Init( offElem, /*basicOnly=*/false );
 640       const SMDS_MeshElement* iniElem = editor.AddElement( nIniIDs, elemType );
 641       initMDS->SetMeshElementOnShape( iniElem, theSolidID );
 642       iniElem->setIsMarked( offElem->isMarked() );
 643       editor.ClearLastCreated();
 644     }
 645     return;
 646   }
 647
 648   //================================================================================
 649   /*!
 650    * \brief set elements of layers boundary to sub-meshes
 651    *  \param [in] theEOS - vlEdges of a sub-shape
 652    *  \param [inout] theMesh - the mesh
 653    *  \param [in] theN2E - map of node to vlEdge
 654    *  \param [in] theNShift - nb of nodes in the mesh before adding the offset mesh
 655    */
 656   //================================================================================
 657
 658   void setBnd2Sub( VLEdgesOnShape &     theEOS,
 659                    SMESH_Mesh*          theInitMesh,
 660                    SMESH_Mesh*          theOffsMesh,
 661                    const TNode2VLEdge & theN2E,
 662                    const smIdType       theNShift,
 663                    TIDSortedNodeSet&    theNodesToCheckCoinc )
 664   {
 665     SMESHDS_Mesh* offsetMDS = theOffsMesh->GetMeshDS();
 666     SMESHDS_Mesh*   initMDS = theInitMesh->GetMeshDS();
 667
 668     SMESHDS_SubMesh* sm = offsetMDS->MeshElements( theEOS._offsetShape );
 669     if ( !sm || ( sm->NbElements() + sm->NbNodes() == 0 ))
 670       return;
 671
 672     for ( SMDS_NodeIteratorPtr nIt = sm->GetNodes(); nIt->more(); )
 673     {
 674       const SMDS_MeshNode* offNode = nIt->next();
 675       VLEdge* vlEd = theN2E( offNode );
 676       const SMDS_MeshNode* nOffBnd = vlEd->_nodes.back().Node();
 677       const SMDS_MeshNode* nIniBnd = initMDS->FindNode( nOffBnd->GetID() + theNShift );
 678       initMDS->UnSetNodeOnShape( nIniBnd );
 679
 680       switch( theEOS._initShape.ShapeType() ) {
 681       case TopAbs_FACE:   initMDS->SetNodeOnFace( nIniBnd, TopoDS::Face( theEOS._initShape ),
 682                                                   vlEd->_uv[0], vlEd->_uv[1] );
 683         break;
 684       case TopAbs_EDGE:   initMDS->SetNodeOnEdge( nIniBnd, TopoDS::Edge( theEOS._initShape ),
 685                                                   vlEd->_uv[0]);
 686         break;
 687       case TopAbs_VERTEX: initMDS->SetNodeOnVertex( nIniBnd, TopoDS::Vertex( theEOS._initShape ));
 688         break;
 689       default:;
 690       }
 691     }
 692
 693     std::vector< const SMDS_MeshNode* > iniNodes, iniNodesBnd;
 694     std::vector< VLEdge*> edgesVec;
 695     for ( SMDS_ElemIteratorPtr eIt = sm->GetElements(); eIt->more(); )
 696     {
 697       const SMDS_MeshElement* offElem = eIt->next();
 698       if ( offElem->GetType() != SMDSAbs_Face &&
 699            offElem->GetType() != SMDSAbs_Edge )
 700         continue;
 701
 702       const int nbNodes = offElem->NbCornerNodes();
 703       iniNodes.resize( nbNodes );
 704       iniNodesBnd.resize( nbNodes );
 705       for ( int i = 0; i < nbNodes; ++i )
 706       {
 707         const SMDS_MeshNode* nOff = offElem->GetNode( i );
 708         const SMDS_MeshNode* nIni = initMDS->FindNode( nOff->GetID() + theNShift );
 709         iniNodes[ i ] = nIni;
 710         if ( !theN2E.IsBound( nOff ))
 711           break;
 712         VLEdge* vlEd = theN2E( nOff );
 713         const SMDS_MeshNode* nOffBnd = vlEd->_nodes.back().Node();
 714         const SMDS_MeshNode* nIniBnd = initMDS->FindNode( nOffBnd->GetID() + theNShift );
 715         iniNodesBnd[ i ] = nIniBnd;
 716       }
 717       if ( const SMDS_MeshElement* iniElem = initMDS->FindElement( iniNodes ))
 718       {
 719         initMDS->UnSetElementOnShape( iniElem );
 720         initMDS->SetMeshElementOnShape( iniElem, theEOS._initShape );
 721
 722         // move the face to the top of prisms, on init mesh boundary
 723         initMDS->ChangeElementNodes( iniElem, iniNodesBnd.data(), nbNodes );
 724
 725         if ( theEOS._toCheckCoinc )
 726           theNodesToCheckCoinc.insert( iniNodesBnd.begin(), iniNodesBnd.end() );
 727       }
 728     }
 729     return;
 730   }
 731
 732   //================================================================================
 733   /*!
 734    * \brief set elements of WOVL face to sub-meshes
 735    *  \param [in] theEOS - vlEdges of a sub-shape
 736    *  \param [inout] theMesh - the mesh
 737    *  \param [in] theN2E - map of node to vlEdge
 738    *  \param [in] theNShift - nb of nodes in the mesh before adding the offset mesh
 739    */
 740   //================================================================================
 741
 742   void setBnd2FVWL( VLEdgesOnShape &     theEOS,
 743                     SMESH_Mesh*          theInitMesh,
 744                     SMESH_Mesh*          theOffsMesh,
 745                     //const TNode2VLEdge & theN2E,
 746                     const smIdType       theNShift )
 747   {
 748     SMESHDS_Mesh* offsetMDS = theOffsMesh->GetMeshDS();
 749     SMESHDS_Mesh*   initMDS = theInitMesh->GetMeshDS();
 750
 751     SMESHDS_SubMesh* sm = offsetMDS->MeshElements( theEOS._offsetShape );
 752     if ( !sm || ( sm->NbElements() + sm->NbNodes() == 0 ))
 753       return;
 754
 755     for ( SMDS_NodeIteratorPtr nIt = sm->GetNodes(); nIt->more(); )
 756     {
 757       const SMDS_MeshNode* offNode = nIt->next();
 758       const SMDS_MeshNode* iniNode = initMDS->FindNode( offNode->GetID() + theNShift );
 759       initMDS->UnSetNodeOnShape( iniNode );
 760
 761       switch( theEOS._initShape.ShapeType() ) {
 762       case TopAbs_FACE:   initMDS->SetNodeOnFace( iniNode, TopoDS::Face( theEOS._initShape ));
 763         break;
 764       case TopAbs_EDGE:   initMDS->SetNodeOnEdge( iniNode, TopoDS::Edge( theEOS._initShape ));
 765         break;
 766       case TopAbs_VERTEX: initMDS->SetNodeOnVertex( iniNode, TopoDS::Vertex( theEOS._initShape ));
 767         break;
 768       default:;
 769       }
 770     }
 771
 772     std::vector< const SMDS_MeshNode* > iniNodes;
 773     for ( SMDS_ElemIteratorPtr eIt = sm->GetElements(); eIt->more(); )
 774     {
 775       const SMDS_MeshElement* offElem = eIt->next();
 776       // if ( offElem->GetType() != SMDSAbs_Face )
 777       //   continue;
 778
 779       int nbNodes = offElem->NbCornerNodes();
 780       iniNodes.resize( nbNodes );
 781       for ( int i = 0; i < nbNodes; ++i )
 782       {
 783         const SMDS_MeshNode* nOff = offElem->GetNode( i );
 784         const SMDS_MeshNode* nIni = initMDS->FindNode( nOff->GetID() + theNShift );
 785         iniNodes[ i ] = nIni;
 786       }
 787       if ( const SMDS_MeshElement* iniElem = initMDS->FindElement( iniNodes ))
 788       {
 789         initMDS->UnSetElementOnShape( iniElem );
 790         initMDS->SetMeshElementOnShape( iniElem, theEOS._initShapeID );
 791
 792         for ( const SMDS_MeshNode* n : iniNodes )
 793           if ( n->GetPosition()->GetDim() == 3 )
 794           {
 795             initMDS->UnSetNodeOnShape( n );
 796             if ( theEOS._initShape.ShapeType() == TopAbs_FACE )
 797               initMDS->SetNodeOnFace( n, theEOS._initShapeID );
 798             else
 799               initMDS->SetNodeOnEdge( n, theEOS._initShapeID );
 800           }
 801       }
 802     }
 803     return;
 804   }
 805 } // namespace
 806
 807
 808 //================================================================================
 809 /*!
 810  * \brief Create a temporary mesh used to hold elements of the offset shape
 811  */
 812 //================================================================================
 813
 814 SMESH_Mesh* StdMeshers_Cartesian_VL::ViscousBuilder::MakeOffsetMesh()
 815 {
 816   return _offsetMesh;
 817 }
 818
 819 //================================================================================
 820 /*!
 821  * \brief Constructor
 822  */
 823 //================================================================================
 824
 825 StdMeshers_Cartesian_VL::ViscousBuilder::ViscousBuilder( const StdMeshers_ViscousLayers* hyp,
 826                                                          const SMESH_Mesh &              mainMesh,
 827                                                          const TopoDS_Shape &            mainShape)
 828   : _hyp( hyp ), _offsetMesh( new TmpMesh )
 829 {
 830   // find shapes with viscous layers
 831   TopTools_IndexedMapOfShape faces;
 832   TopExp::MapShapes( mainShape, TopAbs_FACE, faces );
 833   const SMESHDS_Mesh* iniMDS = mainMesh.GetMeshDS();
 834
 835   if ( hyp->IsToIgnoreShapes() )
 836   {
 837     for ( int i = 1; i <= faces.Size(); ++i )
 838       _shapesWVL.insert( iniMDS->ShapeToIndex( faces( i )));
 839
 840     for ( TGeomID& sID : hyp->GetBndShapes() )
 841       _shapesWVL.erase( sID );
 842   }
 843   else
 844   {
 845     for ( TGeomID& sID : hyp->GetBndShapes() )
 846       _shapesWVL.insert( sID );
 847   }
 848
 849   for ( TGeomID fID : _shapesWVL )
 850   {
 851     const TopoDS_Shape & face = iniMDS->IndexToShape( fID );
 852     for ( TopExp_Explorer exp( face, TopAbs_EDGE ); exp.More(); exp.Next() )
 853       _shapesWVL.insert( iniMDS->ShapeToIndex( exp.Current() ));
 854     for ( TopExp_Explorer exp( face, TopAbs_VERTEX ); exp.More(); exp.Next() )
 855       _shapesWVL.insert( iniMDS->ShapeToIndex( exp.Current() ));
 856   }
 857
 858   // find EDGE and adjacent FACEs w/o VL
 859   for ( int i = 1; i <= faces.Size(); ++i )
 860   {
 861     const TopoDS_Shape& face = faces( i );
 862     TGeomID fID = iniMDS->ShapeToIndex( face );
 863     bool isMainShape2D = (mainShape.ShapeType() == TopAbs_FACE) ? true : false;
 864
 865     if ( _shapesWVL.count( fID ) && !isMainShape2D )
 866       continue;
 867
 868     for ( TopExp_Explorer exp( face, TopAbs_EDGE ); exp.More(); exp.Next() )
 869       _edge2facesWOVL[ iniMDS->ShapeToIndex( exp.Current() )].push_back( fID );
 870   }
 871
 872   // When 2D meshing Need to add edges where segments need to be added due to geometry shrink
 873   return;
 874 }
 875
 876 //================================================================================
 877 /*!
 878  * \brief Destructor
 879  */
 880 //================================================================================
 881
 882 StdMeshers_Cartesian_VL::ViscousBuilder::~ViscousBuilder()
 883 {
 884   delete _offsetMesh; //_offsetMesh = 0;
 885 }
 886
 887 //================================================================================
 888 /*!
 889  * \brief Create an offset solid from a given one
 890  *  \param [in] theShape - input shape can be a solid, solidcompound or a compound with solids
 891  *  \param [in] theMesh - main mesh
 892  *  \param [out] theError - error description
 893  *  \return TopoDS_Shape - result offset shape of the same type as the received shape
 894  */
 895 //================================================================================
 896
 897 TopoDS_Shape StdMeshers_Cartesian_VL::ViscousBuilder::MakeOffsetSolid(const TopoDS_Shape & theShape,
 898                                                                         SMESH_Mesh &         theMesh,
 899                                                                         std::string &        theError )
 900 {
 901   double offset = -_hyp->GetTotalThickness();
 902   double    tol = Precision::Confusion();
 903   TopAbs_ShapeEnum typeOfShape = theShape.ShapeType();
 904
 905   TopTools_IndexedMapOfShape shapeList;
 906   TopExp::MapShapes( theShape, TopAbs_SOLID, shapeList );
 907   std::vector<TopoDS_Shape> shrinkBodies;
 908
 909   for ( int i = 1; i <= shapeList.Size(); ++i )
 910   {
 911     auto solid = shapeList( i );
 912       // If Shape is solid call direct
 913     BRepOffset_MakeOffset * makeOffset = new BRepOffset_MakeOffset();
 914     makeOffset->Initialize( solid, offset, tol, BRepOffset_Skin, /*Intersection=*/false,
 915                           /*selfInter=*/false, GeomAbs_Intersection);
 916
 917     // exclude inlet FACEs
 918     SMESHDS_Mesh* meshDS = theMesh.GetMeshDS();
 919     for ( TopExp_Explorer fEx( theShape, TopAbs_FACE ); fEx.More(); fEx.Next() )
 920     {
 921       TGeomID fID = meshDS->ShapeToIndex( fEx.Current() );
 922       if ( !_shapesWVL.count( fID ))
 923         makeOffset->SetOffsetOnFace( TopoDS::Face( fEx.Current()), 0 );
 924     }
 925
 926     makeOffset->MakeOffsetShape();
 927     if ( makeOffset->IsDone() )
 928     {
 929       shrinkBodies.push_back( makeOffset->Shape() );
 930       _makeOffsetCollection.push_back( makeOffset );
 931     }
 932     else
 933     {
 934       switch ( makeOffset->Error() )
 935       {
 936       case BRepOffset_NoError:
 937         theError = "OK. Offset performed successfully.";break;
 938       case BRepOffset_BadNormalsOnGeometry:
 939         theError = "Degenerated normal on input data.";break;
 940       case BRepOffset_C0Geometry:
 941         theError = "C0 continuity of input data.";break;
 942       case BRepOffset_NullOffset:
 943         theError = "Null offset of all faces.";break;
 944       case BRepOffset_NotConnectedShell:
 945         theError = "Incorrect set of faces to remove, the remaining shell is not connected.";break;
 946       case BRepOffset_CannotTrimEdges:
 947         theError = "Can not trim edges.";break;
 948       case BRepOffset_CannotFuseVertices:
 949         theError = "Can not fuse vertices.";break;
 950       case BRepOffset_CannotExtentEdge:
 951         theError = "Can not extent edge.";break;
 952       default:
 953         theError = "operation not done.";
 954       }
 955       theError = "BRepOffset_MakeOffset error: " + theError;
 956
 957       return TopoDS_Shape();
 958     }
 959   }
 960
 961   if ( typeOfShape == TopAbs_COMPOUND || typeOfShape == TopAbs_COMPSOLID )
 962   {
 963     _solidCompound.SetGlue( BOPAlgo_GlueFull );
 964     _solidCompound.SetToFillHistory( true );
 965     for ( auto solid : shrinkBodies )
 966       _solidCompound.AddArgument( solid );
 967
 968     _solidCompound.Perform();
 969     return _solidCompound.Shape();
 970   }
 971   else
 972     return shrinkBodies[ 0 ]; // return one solid
 973
 974 }
 975
 976 //================================================================================
 977 /*!
 978  * \brief Create an offset shape from a given one
 979  *  \param [in] theShape - input shape
 980  *  \param [in] theMesh - main mesh
 981  *  \param [out] theError - error description
 982  *  \return TopoDS_Shape - result offset shape
 983  */
 984 //================================================================================
 985
 986 TopoDS_Shape
 987 StdMeshers_Cartesian_VL::ViscousBuilder::MakeOffsetShape(const TopoDS_Shape & theShape,
 988                                                          SMESH_Mesh &         theMesh,
 989                                                          std::string &        theError )
 990 {
 991   double offset                 = -_hyp->GetTotalThickness();
 992   double    tol                 = Precision::Confusion();
 993   TopAbs_ShapeEnum typeOfShape  = theShape.ShapeType();
 994
 995   // Switch here for the treatment of faces
 996   if ( typeOfShape == TopAbs_FACE )
 997   {
 998     TopoDS_Face face = TopoDS::Face( theShape );
 999     GProp_GProps gprops;
1000     BRepGProp::SurfaceProperties(face, gprops); // Stores results in gprops
1001     double faceArea = gprops.Mass();
1002
1003     _makeFaceOffset = BRepOffsetAPI_MakeOffset( face, GeomAbs_Intersection );
1004     _makeFaceOffset.Perform( offset );
1005     TopoDS_Wire wireFrame   = TopoDS::Wire( _makeFaceOffset.Shape() );
1006     TopoDS_Face shrinkFace  = TopoDS::Face( BRepBuilderAPI_MakeFace( wireFrame, false ) );
1007     BRepGProp::SurfaceProperties(shrinkFace, gprops); // Stores results in gprops
1008     double sArea = gprops.Mass();
1009
1010     if ( sArea > faceArea /*recompute the shrink face because offset was done in the contrary direction as expected*/)
1011     {
1012       _makeFaceOffset.Perform( -offset );
1013       wireFrame  = TopoDS::Wire( _makeFaceOffset.Shape() );
1014       shrinkFace = TopoDS::Face( BRepBuilderAPI_MakeFace( wireFrame, false ) );
1015     }
1016     _offsetShape = shrinkFace;
1017     return _offsetShape;
1018   }
1019   else
1020   {
1021     _offsetShape = MakeOffsetSolid( theShape, theMesh, theError );
1022     return _offsetShape;
1023   }
1024 }
1025
1026 //================================================================================
1027 /*!
1028  * \brief Return the list of sub-shape of the same type of the offset shape generated from a given initial sub-shape
1029  */
1030 //================================================================================
1031
1032 void StdMeshers_Cartesian_VL::ViscousBuilder::getOffsetSubShape( const TopoDS_Shape& S, std::vector<TopoDS_Shape>& subShapeList )
1033 {
1034   for( auto offset : _makeOffsetCollection )
1035   {
1036     const TopTools_ListOfShape& newShapes = offset->Generated( S );
1037     if ( newShapes.Size() == 0 )
1038       continue; // keep searching
1039
1040     for ( const TopoDS_Shape& ns : newShapes )
1041     {
1042       if ( ns.ShapeType() == S.ShapeType() )
1043       {
1044         if  ( _solidCompound.Arguments().Size() == 0 /* only one solid shrank*/ )
1045         {
1046           subShapeList.push_back( ns );
1047         }
1048         else
1049         {
1050           // In boolean operations the shapes are modified or deleted
1051           const TopTools_ListOfShape& newGlueShapes = _solidCompound.Modified( ns );
1052           for ( TopoDS_Shape& ngs : newGlueShapes )
1053             if ( ngs.ShapeType() == ns.ShapeType() /*&& !ngs.Checked()*/  )
1054               subShapeList.push_back( ngs );
1055
1056           if ( newGlueShapes.Size() == 0 && !_solidCompound.IsDeleted( ns ) )
1057             subShapeList.push_back( ns );
1058         }
1059       }
1060     }
1061   }
1062
1063   // check for _makeFaceOffset in face shrink
1064   if ( _makeOffsetCollection.size()  == 0 )
1065   {
1066     const TopTools_ListOfShape& newShapes = _makeFaceOffset.Generated( S );
1067     for ( const TopoDS_Shape& ns : newShapes )
1068     {
1069       if ( ns.ShapeType() == S.ShapeType() )
1070         return subShapeList.push_back( ns );
1071     }
1072   }
1073 }
1074
1075 bool StdMeshers_Cartesian_VL::ViscousBuilder::CheckGeometryMaps( SMESH_Mesh &        offsetMesh,
1076                                                                   const TopoDS_Shape & theShape  )
1077 {
1078   SMESHDS_Mesh* offsetMDS       = offsetMesh.GetMeshDS();
1079   TopoDS_Shape shrinkGeomToMesh = offsetMDS->ShapeToMesh();
1080
1081   TopTools_IndexedMapOfShape shrinkGeomMap;
1082   TopExp::MapShapes( shrinkGeomToMesh, shrinkGeomMap );
1083   TopTools_IndexedMapOfShape offsetGeomMap;
1084   TopExp::MapShapes( _offsetShape, offsetGeomMap );
1085
1086   // loop on sub-shapes to project nodes from offset boundary to initial boundary
1087   TopAbs_ShapeEnum types[3] = { TopAbs_VERTEX, TopAbs_EDGE, TopAbs_FACE };
1088   for ( TopAbs_ShapeEnum shType : types )
1089   {
1090     TopTools_IndexedMapOfShape shapes;
1091     TopExp::MapShapes( theShape, shType, shapes );
1092     for ( int i = 1; i <= shapes.Size(); ++i )
1093     {
1094       // For each type of geometry check the existence of one or more equivalents
1095       std::vector<TopoDS_Shape> listOfShapes;
1096       getOffsetSubShape( shapes(i), listOfShapes );
1097       if ( listOfShapes.size() == 0 ) return false;
1098     }
1099   }
1100   return true;
1101 }
1102
1103 //================================================================================
1104 /*!
1105  * \brief Create prismatic mesh between _offsetShape and theShape
1106  *  \remark Build the viscous layer from the iteration of shrink geometry
1107  *  \param [out] theMesh - mesh to fill in
1108  *  \param [in] theShape - initial shape
1109  *  \return bool - is Ok
1110  */
1111 //================================================================================
1112
1113 bool StdMeshers_Cartesian_VL::ViscousBuilder::MakeViscousLayers( SMESH_Mesh &        offsetMesh,
1114                                                                  SMESH_Mesh &         theMesh,
1115                                                                  const TopoDS_Shape & theShape )
1116 {
1117   SMESHDS_Mesh* offsetMDS       = offsetMesh.GetMeshDS();
1118   SMESHDS_Mesh*   initMDS       = theMesh.GetMeshDS();
1119   TopoDS_Shape shrinkGeomToMesh = offsetMDS->ShapeToMesh();
1120   bool isMainShape2D = (theShape.ShapeType() == TopAbs_FACE) ? true : false;
1121
1122   // Validate map of shrink+joint geometry elements
1123   if ( !CheckGeometryMaps(offsetMesh, theShape ) && !isMainShape2D )
1124     throw SALOME_Exception("The shrink geometry does not match or respect the original topology.The viscous layer can't be build");
1125
1126
1127   initMDS->ClearMesh(); // avoid mesh superposition on multiple calls of addLayers
1128   offsetMDS->SetAllCellsNotMarked();
1129
1130   TopTools_IndexedMapOfShape shrinkGeomMap;
1131   TopExp::MapShapes( shrinkGeomToMesh, shrinkGeomMap );
1132
1133   TopTools_IndexedMapOfShape offsetGeomMap;
1134   TopExp::MapShapes( _offsetShape, offsetGeomMap );
1135
1136   // Compute heights of viscous layers
1137   std::vector< double > vlH;
1138   computeVLHeight( _hyp, vlH );
1139
1140   std::vector< VLEdgesOnShape > edgesOnShape;
1141   edgesOnShape.reserve( offsetMDS->MaxShapeIndex() + 1 );
1142   TNode2VLEdge n2e;
1143
1144   // loop on sub-shapes to project nodes from offset boundary to initial boundary
1145   TopAbs_ShapeEnum types[3] = { TopAbs_VERTEX, TopAbs_EDGE, TopAbs_FACE };
1146   for ( TopAbs_ShapeEnum shType : types )
1147   {
1148     TopTools_IndexedMapOfShape shapes;
1149     TopExp::MapShapes( theShape, shType, shapes );
1150     for ( int i = 1; i <= shapes.Size(); ++i )
1151     {
1152       edgesOnShape.resize( edgesOnShape.size() + 1 );
1153       VLEdgesOnShape& EOS = edgesOnShape.back();
1154       std::vector<TopoDS_Shape> listOfShapes;
1155       EOS._initShape    = shapes( i );
1156
1157       // Get a list of subShapes of the same type generated from the same face
1158       // It is the case with split objects.
1159       getOffsetSubShape( EOS._initShape, listOfShapes );
1160       for ( TopoDS_Shape& shrinkShape : listOfShapes )
1161       {
1162         int shapeId       = offsetGeomMap.FindIndex( shrinkShape );
1163         EOS._offsetShape  = shrinkGeomMap.FindKey( shapeId );
1164         EOS._initShapeID  = initMDS->ShapeToIndex( EOS._initShape );
1165         EOS._hasVL        = _shapesWVL.count( EOS._initShapeID );
1166
1167         EOS._toCheckCoinc = false;
1168         if ( !EOS._hasVL )
1169           continue;
1170
1171         // project boundary nodes of offset mesh to boundary of init mesh
1172         // (new nodes are created in the offset mesh)
1173         switch( EOS._offsetShape.ShapeType() ) {
1174         case TopAbs_VERTEX:
1175         {
1176           EOS._edges.resize( 1 );
1177           EOS._edges[0]._nodes.resize( 2 );
1178           EOS._edges[0]._nodes[0] = SMESH_Algo::VertexNode( TopoDS::Vertex( EOS._offsetShape ),
1179                                                             offsetMDS );
1180           gp_Pnt offP = BRep_Tool::Pnt( TopoDS::Vertex( EOS._initShape ));
1181           EOS._edges[0]._nodes[1] = offsetMDS->AddNode( offP.X(), offP.Y(), offP.Z() );
1182           //EOS._edges[0]._length   = offP.Distance( EOS._edges[0]._nodes[0] );
1183           n2e.Bind( EOS._edges[0]._nodes[0].Node(), & EOS._edges[0] );
1184           break;
1185         }
1186         case TopAbs_EDGE:
1187         {
1188           projectToEdge( EOS, offsetMDS, n2e, isMainShape2D /* add vertex from edges*/ );
1189           break;
1190         }
1191         case TopAbs_FACE:
1192         {
1193           projectToFace( EOS, offsetMDS, n2e );
1194           break;
1195         }
1196         default:;
1197         }
1198
1199         // create nodes of layers
1200         if ( _hyp->GetNumberLayers() > 1 )
1201         {
1202           //if ( _shapesWVL.count( EOS._initShapeID ))
1203           for ( size_t i = 0; i < EOS._edges.size(); ++i )
1204           {
1205             divideVLEdge( &EOS._edges[ i ], vlH, offsetMDS );
1206           }
1207         }
1208       } // loop on generated shrink shape
1209     }//loop on original shape
1210   } // loop on shape types
1211
1212   // create prisms
1213   bool prismsOk = true;
1214   for ( size_t i = 0; i < edgesOnShape.size(); ++i )
1215   {
1216     VLEdgesOnShape& EOS = edgesOnShape[ i ];
1217     if ( EOS._initShape.ShapeType() == TopAbs_FACE && EOS._hasVL )
1218     {
1219       if ( !makePrisms( EOS, &offsetMesh, n2e ))
1220         prismsOk = false;
1221     }
1222   }
1223
1224   if ( prismsOk )
1225   {
1226     // create faces on FACEs WOVL
1227     for ( size_t i = 0; i < edgesOnShape.size(); ++i )
1228     {
1229       VLEdgesOnShape& EOS = edgesOnShape[ i ];
1230       if ( EOS._initShape.ShapeType() == TopAbs_EDGE && EOS._hasVL )
1231       {
1232         auto e2f = _edge2facesWOVL.find( EOS._initShapeID );
1233         if ( e2f != _edge2facesWOVL.end() && !e2f->second.empty() )
1234         {
1235           TopoDS_Shape  f = initMDS->IndexToShape( e2f->second[0] );
1236           std::vector<TopoDS_Shape> listOfShapes;
1237           getOffsetSubShape( f, listOfShapes );
1238           for( TopoDS_Shape& subShape : listOfShapes )
1239           {
1240             int shapeId   = offsetGeomMap.FindIndex( subShape );
1241             TopoDS_Shape f2 = shrinkGeomMap.FindKey( shapeId );
1242             makeFaces( EOS, & offsetMesh, n2e, offsetMDS->ShapeToIndex( f2 ) );
1243           }
1244         }
1245       }
1246     }
1247   }
1248
1249   if ( isMainShape2D )
1250   {
1251      // create faces on FACEs of the inflate viscous layer in 2D faces
1252     for ( size_t i = 0; i < edgesOnShape.size(); ++i )
1253     {
1254       VLEdgesOnShape& EOS = edgesOnShape[ i ];
1255       if ( EOS._initShape.ShapeType() == TopAbs_EDGE && EOS._hasVL /* iterate in market edges with viscous layer*/)
1256       {
1257         int shapeId = offsetMDS->ShapeToIndex( shrinkGeomToMesh );
1258         makeFaces( EOS, & offsetMesh, n2e, shapeId, isMainShape2D ); // pass face Id of shrink geometry
1259       }
1260     }
1261   }
1262
1263    // copy offset mesh to the main one
1264   initMDS->Modified();
1265   initMDS->CompactMesh();
1266   smIdType nShift = initMDS->NbNodes();
1267   TGeomID solidID = initMDS->ShapeToIndex( theShape );
1268   copyMesh( & offsetMesh, & theMesh, solidID );
1269
1270   if ( !prismsOk )
1271   {
1272     if ( SMESH_subMesh * sm = theMesh.GetSubMesh( theShape ))
1273     {
1274       sm->GetComputeError() =
1275         SMESH_ComputeError::New( COMPERR_ALGO_FAILED,
1276                                  "Viscous layers construction error: bad mesh on offset geometry" );
1277     }
1278     return prismsOk;
1279   }
1280
1281   // set elements of layers boundary to sub-meshes
1282   TIDSortedNodeSet nodesToCheckCoinc;
1283   for ( size_t i = 0; i < edgesOnShape.size(); ++i )
1284   {
1285     VLEdgesOnShape& EOS = edgesOnShape[ i ];
1286     if ( EOS._hasVL )
1287       setBnd2Sub( EOS, &theMesh, &offsetMesh, n2e, nShift, nodesToCheckCoinc );
1288     else
1289       setBnd2FVWL( EOS, &theMesh, &offsetMesh, nShift );
1290   }
1291
1292   // merge coincident nodes
1293   SMESH_MeshEditor editor( &theMesh );
1294   SMESH_MeshEditor::TListOfListOfNodes nodesToMerge;
1295   editor.FindCoincidentNodes( nodesToCheckCoinc, vlH[0]/50., nodesToMerge, false );
1296   editor.MergeNodes( nodesToMerge );
1297
1298   // create a group
1299   if ( !_hyp->GetGroupName().empty() )
1300   {
1301     SMDS_MeshGroup* group = _hyp->CreateGroup( _hyp->GetGroupName(), theMesh, SMDSAbs_Volume );
1302
1303     for ( SMDS_ElemIteratorPtr it = initMDS->elementsIterator(); it->more(); )
1304     {
1305       const SMDS_MeshElement* el = it->next();
1306       if ( el->isMarked() )
1307         group->Add( el );
1308     }
1309   }
1310
1311
1312   return prismsOk;
1313 }