Salome HOME
Cleanup of parallel meshing + documentation
[modules/smesh.git] / src / StdMeshers / StdMeshers_ViscousLayers.cxx
index 21bcc6c3c7a18aeda8190c0104c2464169112f6e..a8d421abaffe3623b76d2e6171d5168497995deb 100644 (file)
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2021  CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2022  CEA/DEN, EDF R&D, OPEN CASCADE
 //
 // This library is free software; you can redistribute it and/or
 // modify it under the terms of the GNU Lesser General Public
 #include "StdMeshers_Quadrangle_2D.hxx"
 #include "StdMeshers_ViscousLayers2D.hxx"
 
+#include <Basics_OCCTVersion.hxx>
+
+#if OCC_VERSION_LARGE < 0x07070000
 #include <Adaptor3d_HSurface.hxx>
+#else
+#include <Adaptor3d_Surface.hxx>
+#endif
+
 #include <BRepAdaptor_Curve.hxx>
 #include <BRepAdaptor_Curve2d.hxx>
 #include <BRepAdaptor_Surface.hxx>
 #include <unordered_map>
 
 #ifdef _DEBUG_
+#ifndef WIN32
 #define __myDEBUG
+#endif
 //#define __NOT_INVALIDATE_BAD_SMOOTH
 //#define __NODES_AT_POS
 #endif
@@ -1816,8 +1825,13 @@ namespace VISCOUS_3D
     //case GeomAbs_SurfaceOfExtrusion:
     case GeomAbs_OffsetSurface:
     {
+#if OCC_VERSION_LARGE < 0x07070000
       Handle(Adaptor3d_HSurface) base = surface.BasisSurface();
       return getRovolutionAxis( base->Surface(), axis );
+#else
+      Handle(Adaptor3d_Surface) base = surface.BasisSurface();
+      return getRovolutionAxis( *base, axis );
+#endif
     }
     default: return false;
     }
@@ -2429,7 +2443,14 @@ bool _ViscousBuilder::findFacesWithLayers(const bool onlyWith)
         break;
       }
       default:
-        return error("Not yet supported case", _sdVec[i]._index);
+        std::ostringstream msg;
+        msg << "Not yet supported case: vertex bounded by ";
+        msg << facesWOL.size();
+        msg << " faces without layer at coordinates (";
+        TopoDS_Vertex v = TopoDS::Vertex(vertex);
+        gp_Pnt p = BRep_Tool::Pnt(v);
+        msg << p.X() << ", " << p.Y() << ", " << p.Z() << ")";
+        return error(msg.str().c_str(), _sdVec[i]._index);
       }
     }
   }
@@ -2775,13 +2796,12 @@ bool _ViscousBuilder::makeLayer(_SolidData& data)
             if ( !setEdgeData( *edge, edgesByGeom[ shapeID ], helper, data ))
               return false;
 
-            if ( edge->_nodes.size() < 2 )
-              edge->Block( data );
-              //data._noShrinkShapes.insert( shapeID );
+            if ( edge->_nodes.size() < 2 && !noShrink )
+              edge->Block( data ); // a sole node is moved only if noShrink
           }
           dumpMove(edge->_nodes.back());
 
-          if ( edge->_cosin > faceMaxCosin )
+          if ( edge->_cosin > faceMaxCosin && edge->_nodes.size() > 1 )
           {
             faceMaxCosin = edge->_cosin;
             maxCosinEdge = edge;
@@ -3037,6 +3057,7 @@ void _ViscousBuilder::limitStepSizeByCurvature( _SolidData& data )
     cnvFace._face = F;
     cnvFace._normalsFixed = false;
     cnvFace._isTooCurved = false;
+    cnvFace._normalsFixedOnBorders = false;
 
     double maxCurvature = cnvFace.GetMaxCurvature( data, eof, surfProp, helper );
     if ( maxCurvature > 0 )
@@ -3905,7 +3926,7 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge&         edge,
         getMeshDS()->SetNodeOnFace( tgtNode, TopoDS::Face( eos._sWOL ), uv.X(), uv.Y() );
     }
 
-    if ( edge._nodes.size() > 1 )
+    //if ( edge._nodes.size() > 1 ) -- allow RISKY_SWOL on noShrink shape
     {
       // check if an angle between a FACE with layers and SWOL is sharp,
       // else the edge should not inflate
@@ -3920,8 +3941,11 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge&         edge,
           geomNorm.Reverse(); // inside the SOLID
         if ( geomNorm * edge._normal < -0.001 )
         {
-          getMeshDS()->RemoveFreeNode( tgtNode, 0, /*fromGroups=*/false );
-          edge._nodes.resize( 1 );
+          if ( edge._nodes.size() > 1 )
+          {
+            getMeshDS()->RemoveFreeNode( tgtNode, 0, /*fromGroups=*/false );
+            edge._nodes.resize( 1 );
+          }
         }
         else if ( realLenFactor > 3 ) ///  -- moved to SetCosin()
           //else if ( edge._lenFactor > 3 )
@@ -4639,7 +4663,9 @@ void _Simplex::SortSimplices(vector<_Simplex>& simplices)
 
 void _ViscousBuilder::makeGroupOfLE()
 {
-#ifdef _DEBUG_
+  if (!SALOME::VerbosityActivated())
+    return;
+
   for ( size_t i = 0 ; i < _sdVec.size(); ++i )
   {
     if ( _sdVec[i]._n2eMap.empty() ) continue;
@@ -4695,7 +4721,6 @@ void _ViscousBuilder::makeGroupOfLE()
              << "'%s-%s' % (faceId1+1, faceId2))");
     dumpFunctionEnd();
   }
-#endif
 }
 
 //================================================================================
@@ -4739,7 +4764,7 @@ void _ViscousBuilder::computeGeomSize( _SolidData& data )
     double thinkness = eos._hyp.GetTotalThickness();
     for ( size_t i = 0; i < eos._edges.size(); ++i )
     {
-      if ( eos._edges[i]->Is( _LayerEdge::BLOCKED )) continue;
+      if ( eos._edges[i]->_nodes.size() < 2 ) continue;
       eos._edges[i]->SetMaxLen( thinkness );
       eos._edges[i]->FindIntersection( *searcher, intersecDist, data._epsilon, eos, &face );
       if ( intersecDist > 0 && face )
@@ -4930,7 +4955,7 @@ bool _ViscousBuilder::inflate(_SolidData& data)
         if ( eos._edges[i]->_nodes.size() > 1 )
           avgThick    += Min( 1., eos._edges[i]->_len / shapeTgtThick );
         else
-          avgThick    += shapeTgtThick;
+          avgThick    += 1;
         nbActiveEdges += ( ! eos._edges[i]->Is( _LayerEdge::BLOCKED ));
       }
     }
@@ -5444,8 +5469,11 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data,
 
         // intersection not ignored
 
-        if ( toBlockInfaltion &&
-             dist < ( eos._edges[i]->_len * theThickToIntersection ))
+        double minDist = 0;
+        if ( eos._edges[i]->_maxLen < 0.99 * eos._hyp.GetTotalThickness() ) // limited length
+          minDist = eos._edges[i]->_len * theThickToIntersection;
+
+        if ( toBlockInfaltion && dist < minDist  )
         {
           if ( is1stBlocked ) { is1stBlocked = false; // debug
             dumpFunction(SMESH_Comment("blockIntersected") <<data._index<<"_InfStep"<<infStep);
@@ -5827,27 +5855,26 @@ void _ViscousBuilder::putOnOffsetSurface( _EdgesOnShape&            eos,
     }
   }
 
-
-
-#ifdef _DEBUG_
-  // dumpMove() for debug
-  size_t i = 0;
-  for ( ; i < eos._edges.size(); ++i )
-    if ( eos._edges[i]->Is( _LayerEdge::MARKED ))
-      break;
-  if ( i < eos._edges.size() )
+  if (SALOME::VerbosityActivated())
   {
-    dumpFunction(SMESH_Comment("putOnOffsetSurface_") << eos.ShapeTypeLetter() << eos._shapeID
-                 << "_InfStep" << infStep << "_" << Abs( smooStep ));
+    // dumpMove() for debug
+    size_t i = 0;
     for ( ; i < eos._edges.size(); ++i )
+      if ( eos._edges[i]->Is( _LayerEdge::MARKED ))
+        break;
+    if ( i < eos._edges.size() )
     {
-      if ( eos._edges[i]->Is( _LayerEdge::MARKED )) {
-        dumpMove( eos._edges[i]->_nodes.back() );
+      dumpFunction(SMESH_Comment("putOnOffsetSurface_") << eos.ShapeTypeLetter() << eos._shapeID
+                  << "_InfStep" << infStep << "_" << Abs( smooStep ));
+      for ( ; i < eos._edges.size(); ++i )
+      {
+        if ( eos._edges[i]->Is( _LayerEdge::MARKED )) {
+          dumpMove( eos._edges[i]->_nodes.back() );
+        }
       }
+      dumpFunctionEnd();
     }
-    dumpFunctionEnd();
   }
-#endif
 
   _ConvexFace* cnvFace;
   if ( moveAll != _LayerEdge::UPD_NORMAL_CONV &&
@@ -6192,9 +6219,11 @@ bool _Smoother1D::smoothAnalyticEdge( _SolidData&                    data,
           tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() );
           dumpMove( tgtNode );
 
-          SMDS_FacePositionPtr pos = tgtNode->GetPosition();
-          pos->SetUParameter( newUV.X() );
-          pos->SetVParameter( newUV.Y() );
+          if ( SMDS_FacePositionPtr pos = tgtNode->GetPosition() ) // NULL if F is noShrink
+          {
+            pos->SetUParameter( newUV.X() );
+            pos->SetVParameter( newUV.Y() );
+          }
 
           gp_XYZ newUV0( newUV.X(), newUV.Y(), 0 );
 
@@ -6304,10 +6333,11 @@ bool _Smoother1D::smoothAnalyticEdge( _SolidData&                    data,
         tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() );
         dumpMove( tgtNode );
 
-        SMDS_FacePositionPtr pos = tgtNode->GetPosition();
-        pos->SetUParameter( newUV.X() );
-        pos->SetVParameter( newUV.Y() );
-
+        if ( SMDS_FacePositionPtr pos = tgtNode->GetPosition() ) // NULL if F is noShrink
+        {
+          pos->SetUParameter( newUV.X() );
+          pos->SetVParameter( newUV.Y() );
+        }
         _eos[i]->Set( _LayerEdge::SMOOTHED ); // to check in refine() (IPAL54237)
       }
     }
@@ -6555,6 +6585,7 @@ void _Smoother1D::prepare(_SolidData& data)
   _edgeDir[1] = getEdgeDir( E, leOnV[1]->_nodes[0], data.GetHelper() );
   _leOnV[0]._cosin = Abs( leOnV[0]->_cosin );
   _leOnV[1]._cosin = Abs( leOnV[1]->_cosin );
+  _leOnV[0]._flags = _leOnV[1]._flags = 0;
   if ( _eos._sWOL.IsNull() ) // 3D
     for ( int iEnd = 0; iEnd < 2; ++iEnd )
       _leOnV[iEnd]._cosin = Abs( _edgeDir[iEnd].Normalized() * leOnV[iEnd]->_normal );
@@ -6939,14 +6970,14 @@ void _SolidData::PrepareEdgesToSmoothOnFace( _EdgesOnShape* eos, bool substitute
     eos->_edgeForOffset = 0;
 
     double maxCosin = -1;
-    bool hasNoShrink = false;
+    //bool hasNoShrink = false;
     for ( TopExp_Explorer eExp( eos->_shape, TopAbs_EDGE ); eExp.More(); eExp.Next() )
     {
       _EdgesOnShape* eoe = GetShapeEdges( eExp.Current() );
       if ( !eoe || eoe->_edges.empty() ) continue;
 
-      if ( eos->GetData()._noShrinkShapes.count( eoe->_shapeID ))
-        hasNoShrink = true;
+      // if ( eos->GetData()._noShrinkShapes.count( eoe->_shapeID ))
+      //   hasNoShrink = true;
 
       vector<_LayerEdge*>& eE = eoe->_edges;
       _LayerEdge* e = eE[ eE.size() / 2 ];
@@ -6984,8 +7015,8 @@ void _SolidData::PrepareEdgesToSmoothOnFace( _EdgesOnShape* eos, bool substitute
 
     // Try to initialize _Mapper2D
 
-    if ( hasNoShrink )
-      return;
+    // if ( hasNoShrink )
+    //   return;
 
     SMDS_ElemIteratorPtr fIt = eos->_subMesh->GetSubMeshDS()->GetElements();
     if ( !fIt->more() || fIt->next()->NbCornerNodes() != 4 )
@@ -10497,7 +10528,8 @@ bool _ViscousBuilder::refine(_SolidData& data)
       if ( n2eMap && (( n2e = n2eMap->find( edge._nodes[0] )) != n2eMap->end() ))
       {
         edgeOnSameNode = n2e->second;
-        useExistingPos = ( edgeOnSameNode->_len < edge._len );
+        useExistingPos = ( edgeOnSameNode->_len < edge._len ||
+                           segLen[0] == segLen.back() ); // too short inflation step (bos #20643)
         const gp_XYZ& otherTgtPos = edgeOnSameNode->_pos.back();
         SMDS_PositionPtr  lastPos = tgtNode->GetPosition();
         if ( isOnEdge )
@@ -10840,12 +10872,20 @@ namespace VISCOUS_3D
       std::vector< SMESH_NodeXYZ > _nodes;
       TopAbs_ShapeEnum             _vertSWOLType[2]; // shrink part includes VERTEXes
       AverageHyp*                  _vertHyp[2];
+      double                       _edgeWOLLen[2]; // length of wol EDGE
+      double                       _tol; // to compare _edgeWOLLen's
 
       BndPart():
         _isShrink(0), _isReverse(0), _nbSegments(0), _hyp(0),
-        _vertSWOLType{ TopAbs_WIRE, TopAbs_WIRE }, _vertHyp{ 0, 0 }
+        _vertSWOLType{ TopAbs_WIRE, TopAbs_WIRE }, _vertHyp{ 0, 0 }, _edgeWOLLen{ 0., 0.}
       {}
 
+      bool IsEqualLengthEWOL( const BndPart& other ) const
+      {
+        return ( std::abs( _edgeWOLLen[0] - other._edgeWOLLen[0] ) < _tol &&
+                 std::abs( _edgeWOLLen[1] - other._edgeWOLLen[1] ) < _tol );
+      }
+
       bool operator==( const BndPart& other ) const
       {
         return ( _isShrink       == other._isShrink &&
@@ -10856,7 +10896,8 @@ namespace VISCOUS_3D
                  (( !_isShrink ) ||
                   ( *_hyp        == *other._hyp &&
                     vertHyp1()   == other.vertHyp1() &&
-                    vertHyp2()   == other.vertHyp2() ))
+                    vertHyp2()   == other.vertHyp2() &&
+                    IsEqualLengthEWOL( other )))
                  );
       }
       bool CanAppend( const BndPart& other )
@@ -10874,10 +10915,12 @@ namespace VISCOUS_3D
         bool hasCommonNode = ( _nodes.back()->GetID() == other._nodes.front()->GetID() );
         _nodes.insert( _nodes.end(), other._nodes.begin() + hasCommonNode, other._nodes.end() );
         _vertSWOLType[1] = other._vertSWOLType[1];
-        if ( _isShrink )
-          _vertHyp[1] = other._vertHyp[1];
+        if ( _isShrink ) {
+          _vertHyp[1]    = other._vertHyp[1];
+          _edgeWOLLen[1] = other._edgeWOLLen[1];
+        }
       }
-      const SMDS_MeshNode*    Node(size_t i)  const
+      const SMDS_MeshNode* Node(size_t i)  const
       {
         return _nodes[ _isReverse ? ( _nodes.size() - 1 - i ) : i ]._node;
       }
@@ -11025,6 +11068,11 @@ namespace VISCOUS_3D
       for ( int iE = 0; iE < nbEdgesInWire.front(); ++iE )
       {
         BndPart bndPart;
+
+        std::vector<const SMDS_MeshNode*> nodes = fSide.GetOrderedNodes( iE );
+        bndPart._nodes.assign( nodes.begin(), nodes.end() );
+        bndPart._nbSegments = bndPart._nodes.size() - 1;
+
         _EdgesOnShape*  eos = _data1->GetShapeEdges( fSide.EdgeID( iE ));
 
         bndPart._isShrink = ( eos->SWOLType() == TopAbs_FACE );
@@ -11053,10 +11101,14 @@ namespace VISCOUS_3D
                 bndPart._vertSWOLType[iV] = eov[iV]->SWOLType();
             }
           }
+          bndPart._edgeWOLLen[0] = fSide.EdgeLength( iE - 1 );
+          bndPart._edgeWOLLen[1] = fSide.EdgeLength( iE + 1 );
+
+          bndPart._tol = std::numeric_limits<double>::max(); // tolerance by segment size
+          for ( size_t i = 1; i < bndPart._nodes.size(); ++i )
+            bndPart._tol = Min( bndPart._tol,
+                                ( bndPart._nodes[i-1] - bndPart._nodes[i] ).SquareModulus() );
         }
-        std::vector<const SMDS_MeshNode*> nodes = fSide.GetOrderedNodes( iE );
-        bndPart._nodes.assign( nodes.begin(), nodes.end() );
-        bndPart._nbSegments = bndPart._nodes.size() - 1;
 
         if ( _boundary.empty() || ! _boundary.back().CanAppend( bndPart ))
           _boundary.push_back( bndPart );
@@ -11163,6 +11215,9 @@ namespace VISCOUS_3D
     }
     SMESHDS_Mesh* meshDS = dataSrc->GetHelper().GetMeshDS();
 
+    dumpFunction(SMESH_Comment("periodicMoveNodes_F")
+                               << _shriFace[iSrc]->_subMesh->GetId() << "_F"
+                               << _shriFace[iTgt]->_subMesh->GetId() );
     TNode2Edge::iterator n2e;
     TNodeNodeMap::iterator n2n = _nnMap.begin();
     for ( ; n2n != _nnMap.end(); ++n2n )
@@ -11192,6 +11247,8 @@ namespace VISCOUS_3D
           SMESH_NodeXYZ pSrc = leSrc->_nodes[ iN ];
           gp_XYZ pTgt = trsf->Transform( pSrc );
           meshDS->MoveNode( leTgt->_nodes[ iN ], pTgt.X(), pTgt.Y(), pTgt.Z() );
+
+          dumpMove( leTgt->_nodes[ iN ]);
         }
       }
     }
@@ -11200,6 +11257,7 @@ namespace VISCOUS_3D
               << _shriFace[iSrc]->_subMesh->GetId() << " -> "
               << _shriFace[iTgt]->_subMesh->GetId() << " -- "
               << ( done ? "DONE" : "FAIL"));
+    dumpFunctionEnd();
 
     return done;
   }
@@ -11734,6 +11792,8 @@ bool _ViscousBuilder::shrink(_SolidData& theData)
       {
         _EdgesOnShape& eos = * subEOS[ iS ];
         if ( eos.ShapeType() != TopAbs_EDGE ) continue;
+        if ( eos.size() == 0 )
+          continue;
 
         const TopoDS_Edge& E = TopoDS::Edge( eos._shape );
         data.SortOnEdge( E, eos._edges );
@@ -11756,8 +11816,19 @@ bool _ViscousBuilder::shrink(_SolidData& theData)
           uvPtVec[ i ].param = helper.GetNodeU( E, edges[i]->_nodes[0] );
           uvPtVec[ i ].SetUV( helper.GetNodeUV( F, edges[i]->_nodes.back() ));
         }
-        if ( edges.empty() )
-          continue;
+        if ( uvPtVec[ 0 ].node == uvPtVec.back().node &&            // closed
+             helper.IsSeamShape( uvPtVec[ 0 ].node->GetShapeID() ))
+        {
+          uvPtVec[ 0 ].SetUV( helper.GetNodeUV( F,
+                                                edges[0]->_nodes.back(),
+                                                edges[1]->_nodes.back() ));
+          size_t i = edges.size() - 1;
+          uvPtVec[ i ].SetUV( helper.GetNodeUV( F,
+                                                edges[i  ]->_nodes.back(),
+                                                edges[i-1]->_nodes.back() ));
+        }
+        // if ( edges.empty() )
+        //   continue;
         BRep_Tool::Range( E, uvPtVec[0].param, uvPtVec.back().param );
         StdMeshers_FaceSide fSide( uvPtVec, F, E, _mesh );
         StdMeshers_ViscousLayers2D::SetProxyMeshOfEdge( fSide );
@@ -11772,8 +11843,12 @@ bool _ViscousBuilder::shrink(_SolidData& theData)
       smDS->Clear();
 
       // compute the mesh on the FACE
+      TopTools_IndexedMapOfShape allowed(1);
+      allowed.Add( sm->GetSubShape() );
+      sm->SetAllowedSubShapes( &allowed );
       sm->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE );
       sm->ComputeStateEngine( SMESH_subMesh::COMPUTE_SUBMESH );
+      sm->SetAllowedSubShapes( nullptr );
 
       // re-fill proxy sub-meshes of the FACE
       for ( size_t i = 0 ; i < _sdVec.size(); ++i )
@@ -12374,17 +12449,16 @@ gp_XY _SmoothNode::computeAngularPos(vector<gp_XY>& uv,
   edgeSize.back() = edgeSize.front();
 
   gp_XY  newPos(0,0);
-  //int    nbEdges = 0;
-  double sumSize = 0;
+  double sumWgt = 0;
   for ( size_t i = 1; i < edgeDir.size(); ++i )
   {
-    if ( edgeDir[i-1].X() > 1. ) continue;
-    int i1 = i-1;
+    const int i1 = i-1;
+    if ( edgeDir[i1].X() > 1. ) continue;
     while ( edgeDir[i].X() > 1. && ++i < edgeDir.size() );
     if ( i == edgeDir.size() ) break;
     gp_XY p = uv[i];
     gp_XY norm1( -edgeDir[i1].Y(), edgeDir[i1].X() );
-    gp_XY norm2( -edgeDir[i].Y(),  edgeDir[i].X() );
+    gp_XY norm2( -edgeDir[i ].Y(), edgeDir[i ].X() );
     gp_XY bisec = norm1 + norm2;
     double bisecSize = bisec.Modulus();
     if ( bisecSize < numeric_limits<double>::min() )
@@ -12394,16 +12468,16 @@ gp_XY _SmoothNode::computeAngularPos(vector<gp_XY>& uv,
     }
     bisec /= bisecSize;
 
-    gp_XY  dirToN  = uvToFix - p;
-    double distToN = dirToN.Modulus();
+    gp_XY   dirToN = uvToFix - p;
+    double distToN = bisec * dirToN;
     if ( bisec * dirToN < 0 )
       distToN = -distToN;
 
-    newPos += ( p + bisec * distToN ) * ( edgeSize[i1] + edgeSize[i] );
-    //++nbEdges;
-    sumSize += edgeSize[i1] + edgeSize[i];
+    double wgt = edgeSize[i1] + edgeSize[i];
+    newPos += ( p + bisec * distToN ) * wgt;
+    sumWgt += wgt;
   }
-  newPos /= /*nbEdges * */sumSize;
+  newPos /= sumWgt;
   return newPos;
 }