Salome HOME
Revert "Merge branch 'yan/parallel_mesh2'"
[modules/smesh.git] / src / StdMeshers / StdMeshers_ViscousLayers.cxx
index 4ed66607ada5ded28596be8065e4f7208ce025f0..680708d7ea81b1578e264bf8870c206a848ef67a 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 <unordered_map>
 
 #ifdef _DEBUG_
+#ifndef WIN32
 #define __myDEBUG
+#endif
 //#define __NOT_INVALIDATE_BAD_SMOOTH
 //#define __NODES_AT_POS
 #endif
@@ -2775,13 +2777,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 && !edge->Is( _LayerEdge::BLOCKED ))
+          if ( edge->_cosin > faceMaxCosin && edge->_nodes.size() > 1 )
           {
             faceMaxCosin = edge->_cosin;
             maxCosinEdge = edge;
@@ -3037,6 +3038,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 +3907,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 +3922,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 )
@@ -4739,7 +4744,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 +4935,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 +5449,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);
@@ -6192,9 +6200,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 +6314,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 +6566,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 +6951,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 +6996,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 )
@@ -10841,12 +10853,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 &&
@@ -10857,7 +10877,8 @@ namespace VISCOUS_3D
                  (( !_isShrink ) ||
                   ( *_hyp        == *other._hyp &&
                     vertHyp1()   == other.vertHyp1() &&
-                    vertHyp2()   == other.vertHyp2() ))
+                    vertHyp2()   == other.vertHyp2() &&
+                    IsEqualLengthEWOL( other )))
                  );
       }
       bool CanAppend( const BndPart& other )
@@ -10875,10 +10896,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;
       }
@@ -11026,6 +11049,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 );
@@ -11054,10 +11082,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 );
@@ -11164,6 +11196,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 )
@@ -11193,6 +11228,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 ]);
         }
       }
     }
@@ -11201,6 +11238,7 @@ namespace VISCOUS_3D
               << _shriFace[iSrc]->_subMesh->GetId() << " -> "
               << _shriFace[iTgt]->_subMesh->GetId() << " -- "
               << ( done ? "DONE" : "FAIL"));
+    dumpFunctionEnd();
 
     return done;
   }
@@ -11735,6 +11773,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 );
@@ -11757,8 +11797,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 );
@@ -11773,8 +11824,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 )
@@ -12375,17 +12430,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() )
@@ -12395,16 +12449,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;
 }