Salome HOME
Regression of 21397: EDF SMESH: a quadrangle face mesh can't be projected to a cylinder
[modules/smesh.git] / src / StdMeshers / StdMeshers_ViscousLayers.cxx
index 61acb13c5c9178699857ce8e53186b6aa84f9245..8167f864a2b4098759020098080ea1fccfe34517 100644 (file)
@@ -1590,7 +1590,7 @@ SMESH_ComputeErrorPtr _ViscousBuilder::CheckHypotheses( SMESH_Mesh&         mesh
 
 
   findSolidsWithLayers();
-  bool ok = findFacesWithLayers();
+  bool ok = findFacesWithLayers( true );
 
   // remove _MeshOfSolid's of _SolidData's
   for ( size_t i = 0; i < _sdVec.size(); ++i )
@@ -2371,14 +2371,14 @@ void _ViscousBuilder::limitStepSizeByCurvature( _SolidData& data )
 
   for ( size_t iS = 0; iS < data._edgesOnShape.size(); ++iS )
   {
-    _EdgesOnShape& eos = data._edgesOnShape[iS];
-    if ( eos.ShapeType() != TopAbs_FACE ||
-         data._ignoreFaceIds.count( eos._shapeID ))
+    _EdgesOnShape& eof = data._edgesOnShape[iS];
+    if ( eof.ShapeType() != TopAbs_FACE ||
+         data._ignoreFaceIds.count( eof._shapeID ))
       continue;
 
-    TopoDS_Face        F = TopoDS::Face( eos._shape );
-    SMESH_subMesh *   sm = eos._subMesh;
-    const TGeomID faceID = eos._shapeID;
+    TopoDS_Face        F = TopoDS::Face( eof._shape );
+    SMESH_subMesh *   sm = eof._subMesh;
+    const TGeomID faceID = eof._shapeID;
 
     BRepAdaptor_Surface surface( F, false );
     surfProp.SetSurface( surface );
@@ -2393,18 +2393,18 @@ void _ViscousBuilder::limitStepSizeByCurvature( _SolidData& data )
       sm = smIt->next();
       const TGeomID subID = sm->GetId();
       // find _LayerEdge's of a sub-shape
-      size_t edgesEnd;
-      if ( _EdgesOnShape* eos = data.GetShapeEdges( subID ))
+      _EdgesOnShape* eos;
+      if (( eos = data.GetShapeEdges( subID )))
         cnvFace._subIdToEOS.insert( make_pair( subID, eos ));
       else
         continue;
       // check concavity and curvature and limit data._stepSize
       const double minCurvature =
-        1. / ( eos._hyp.GetTotalThickness() * ( 1+theThickToIntersection ));
-      size_t iStep = Max( 1, eos._edges.size() / nbTestPnt );
-      for ( size_t i = 0; i < eos._edges.size(); i += iStep )
+        1. / ( eos->_hyp.GetTotalThickness() * ( 1+theThickToIntersection ));
+      size_t iStep = Max( 1, eos->_edges.size() / nbTestPnt );
+      for ( size_t i = 0; i < eos->_edges.size(); i += iStep )
       {
-        gp_XY uv = helper.GetNodeUV( F, eos._edges[ i ]->_nodes[0] );
+        gp_XY uv = helper.GetNodeUV( F, eos->_edges[ i ]->_nodes[0] );
         surfProp.SetParameters( uv.X(), uv.Y() );
         if ( !surfProp.IsCurvatureDefined() )
           continue;
@@ -2458,8 +2458,8 @@ void _ViscousBuilder::limitStepSizeByCurvature( _SolidData& data )
       set< const SMDS_MeshNode* > usedNodes;
 
       // look for _LayerEdge's with null _sWOL
-      map< TGeomID, _EdgesOnShape* >::iterator id2oes = convFace._subIdToEOS.begin();
-      for ( ; id2oes != convFace._subIdToEOS.end(); ++id2oes )
+      id2eos = convFace._subIdToEOS.begin();
+      for ( ; id2eos != convFace._subIdToEOS.end(); ++id2eos )
       {
         _EdgesOnShape& eos = * id2eos->second;
         if ( !eos._sWOL.IsNull() )
@@ -2920,6 +2920,7 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge&         edge,
   }
   else // !useGeometry - get _normal using surrounding mesh faces
   {
+    set<TGeomID> faceIds;
 
     SMDS_ElemIteratorPtr fIt = node->GetInverseElementIterator(SMDSAbs_Face);
     while ( fIt->more() )
@@ -2927,6 +2928,8 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge&         edge,
       const SMDS_MeshElement* face = fIt->next();
       if ( eos.GetNormal( face, geomNorm ))
       {
+        if ( onShrinkShape && !faceIds.insert( face->getshapeId() ).second )
+          continue; // use only one mesh face on FACE
         edge._normal += geomNorm.XYZ();
         totalNbFaces++;
       }
@@ -3236,74 +3239,87 @@ gp_XYZ _ViscousBuilder::getWeigthedNormal( const SMDS_MeshNode*             n,
   }
 
   // exclude equal normals
-  //int nbUniqNorms = nbFaces;
-  for ( int i = 0; i < nbFaces; ++i )
+  int nbUniqNorms = nbFaces;
+  for ( int i = 0; i < nbFaces; ++i ) {
     for ( int j = i+1; j < nbFaces; ++j )
       if ( fId2Normal[i].second.IsEqual( fId2Normal[j].second, 0.1 ))
       {
         fId2Normal[i].second.SetCoord( 0,0,0 );
-        //--nbUniqNorms;
+        --nbUniqNorms;
         break;
       }
-  //if ( nbUniqNorms < 3 )
-  {
-    for ( int i = 0; i < nbFaces; ++i )
-      resNorm += fId2Normal[i].second;
-    return resNorm;
   }
-
-  double angles[30];
   for ( int i = 0; i < nbFaces; ++i )
-  {
-    const TopoDS_Face& F = fId2Normal[i].first;
+    resNorm += fId2Normal[i].second;
 
-    // look for two EDGEs shared by F and other FACEs within fId2Normal
-    TopoDS_Edge ee[2];
-    int nbE = 0;
-    PShapeIteratorPtr eIt = SMESH_MesherHelper::GetAncestors( V, *_mesh, TopAbs_EDGE );
-    while ( const TopoDS_Shape* E = eIt->next() )
-    {
-      if ( !SMESH_MesherHelper::IsSubShape( *E, F ))
-        continue;
-      bool isSharedEdge = false;
-      for ( int j = 0; j < nbFaces && !isSharedEdge; ++j )
-      {
-        if ( i == j ) continue;
-        const TopoDS_Shape& otherF = fId2Normal[j].first;
-        isSharedEdge = SMESH_MesherHelper::IsSubShape( *E, otherF );
-      }
-      if ( !isSharedEdge )
-        continue;
-      ee[ nbE ] = TopoDS::Edge( *E );
-      ee[ nbE ].Orientation( SMESH_MesherHelper::GetSubShapeOri( F, *E ));
-      if ( ++nbE == 2 )
-        break;
-    }
-
-    // get an angle between the two EDGEs
-    angles[i] = 0;
-    if ( nbE < 1 ) continue;
-    if ( nbE == 1 )
-    {
-      ee[ 1 ] == ee[ 0 ];
-    }
-    else
+  // assure that resNorm is visible by every FACE (IPAL0052675)
+  if ( nbUniqNorms > 3 )
+  {
+    bool change = false;
+    for ( int nbAttempts = 0; nbAttempts < nbFaces; ++nbAttempts)
     {
-      if ( !V.IsSame( SMESH_MesherHelper::IthVertex( 0, ee[ 1 ] )))
-        std::swap( ee[0], ee[1] );
+      for ( int i = 0; i < nbFaces; ++i )
+        if ( resNorm * fId2Normal[i].second < 0.5 )
+        {
+          resNorm += fId2Normal[i].second;
+          change = true;
+        }
+      if ( !change ) break;
     }
-    angles[i] = SMESH_MesherHelper::GetAngle( ee[0], ee[1], F, TopoDS::Vertex( V ));
   }
 
-  // compute a weighted normal
-  double sumAngle = 0;
-  for ( int i = 0; i < nbFaces; ++i )
-  {
-    angles[i] = ( angles[i] > 2*M_PI )  ?  0  :  M_PI - angles[i];
-    sumAngle += angles[i];
-  }
-  for ( int i = 0; i < nbFaces; ++i )
-    resNorm += angles[i] / sumAngle * fId2Normal[i].second;
+  // double angles[30];
+  // for ( int i = 0; i < nbFaces; ++i )
+  // {
+  //   const TopoDS_Face& F = fId2Normal[i].first;
+
+  //   // look for two EDGEs shared by F and other FACEs within fId2Normal
+  //   TopoDS_Edge ee[2];
+  //   int nbE = 0;
+  //   PShapeIteratorPtr eIt = SMESH_MesherHelper::GetAncestors( V, *_mesh, TopAbs_EDGE );
+  //   while ( const TopoDS_Shape* E = eIt->next() )
+  //   {
+  //     if ( !SMESH_MesherHelper::IsSubShape( *E, F ))
+  //       continue;
+  //     bool isSharedEdge = false;
+  //     for ( int j = 0; j < nbFaces && !isSharedEdge; ++j )
+  //     {
+  //       if ( i == j ) continue;
+  //       const TopoDS_Shape& otherF = fId2Normal[j].first;
+  //       isSharedEdge = SMESH_MesherHelper::IsSubShape( *E, otherF );
+  //     }
+  //     if ( !isSharedEdge )
+  //       continue;
+  //     ee[ nbE ] = TopoDS::Edge( *E );
+  //     ee[ nbE ].Orientation( SMESH_MesherHelper::GetSubShapeOri( F, *E ));
+  //     if ( ++nbE == 2 )
+  //       break;
+  //   }
+
+  //   // get an angle between the two EDGEs
+  //   angles[i] = 0;
+  //   if ( nbE < 1 ) continue;
+  //   if ( nbE == 1 )
+  //   {
+  //     ee[ 1 ] == ee[ 0 ];
+  //   }
+  //   else
+  //   {
+  //     if ( !V.IsSame( SMESH_MesherHelper::IthVertex( 0, ee[ 1 ] )))
+  //       std::swap( ee[0], ee[1] );
+  //   }
+  //   angles[i] = SMESH_MesherHelper::GetAngle( ee[0], ee[1], F, TopoDS::Vertex( V ));
+  // }
+
+  // // compute a weighted normal
+  // double sumAngle = 0;
+  // for ( int i = 0; i < nbFaces; ++i )
+  // {
+  //   angles[i] = ( angles[i] > 2*M_PI )  ?  0  :  M_PI - angles[i];
+  //   sumAngle += angles[i];
+  // }
+  // for ( int i = 0; i < nbFaces; ++i )
+  //   resNorm += angles[i] / sumAngle * fId2Normal[i].second;
 
   return resNorm;
 }
@@ -4731,10 +4747,10 @@ bool _ViscousBuilder::updateNormalsOfConvexFaces( _SolidData&         data,
     Bnd_B3d nodesBox;
     gp_Pnt  center;
 
-    map< TGeomID, _EdgesOnShape* >::iterator id2oes = convFace._subIdToEOS.begin();
-    for ( ; id2oes != convFace._subIdToEOS.end(); ++id2oes )
+    map< TGeomID, _EdgesOnShape* >::iterator id2eos = convFace._subIdToEOS.begin();
+    for ( ; id2eos != convFace._subIdToEOS.end(); ++id2eos )
     {
-      _EdgesOnShape& eos = *(id2oes->second);
+      _EdgesOnShape& eos = *(id2eos->second);
       if ( eos.ShapeType() == TopAbs_VERTEX )
       {
         _LayerEdge* ledge = eos._edges[ 0 ];
@@ -4764,10 +4780,10 @@ bool _ViscousBuilder::updateNormalsOfConvexFaces( _SolidData&         data,
 
       gp_XYZ avgNormal( 0,0,0 );
       nbEdges = 0;
-      id2oes = convFace._subIdToEOS.begin();
-      for ( ; id2oes != convFace._subIdToEOS.end(); ++id2oes )
+      id2eos = convFace._subIdToEOS.begin();
+      for ( ; id2eos != convFace._subIdToEOS.end(); ++id2eos )
       {
-        _EdgesOnShape& eos = *(id2oes->second);
+        _EdgesOnShape& eos = *(id2eos->second);
         // set data of _CentralCurveOnEdge
         if ( eos.ShapeType() == TopAbs_EDGE )
         {
@@ -4817,10 +4833,10 @@ bool _ViscousBuilder::updateNormalsOfConvexFaces( _SolidData&         data,
         avgCosin /= nbCosin;
 
       // set _LayerEdge::_normal = avgNormal
-      id2oes = convFace._subIdToEOS.begin();
-      for ( ; id2oes != convFace._subIdToEOS.end(); ++id2oes )
+      id2eos = convFace._subIdToEOS.begin();
+      for ( ; id2eos != convFace._subIdToEOS.end(); ++id2eos )
       {
-        _EdgesOnShape& eos = *(id2oes->second);
+        _EdgesOnShape& eos = *(id2eos->second);
         if ( eos.ShapeType() != TopAbs_EDGE )
           for ( size_t i = 0; i < eos._edges.size(); ++i )
             eos._edges[ i ]->_cosin = avgCosin;
@@ -4978,12 +4994,12 @@ bool _ViscousBuilder::updateNormalsOfConvexFaces( _SolidData&         data,
       if ( nbCosin > 0 )
         avgCosin /= nbCosin;
       const TGeomID faceID = meshDS->ShapeToIndex( convFace._face );
-      map< TGeomID, _EdgesOnShape* >::iterator id2oes = convFace._subIdToEOS.find( faceID );
-      if ( id2oes != convFace._subIdToEOS.end() )
+      map< TGeomID, _EdgesOnShape* >::iterator id2eos = convFace._subIdToEOS.find( faceID );
+      if ( id2eos != convFace._subIdToEOS.end() )
       {
         int iE = 0;
         gp_XYZ newNorm;
-        _EdgesOnShape& eos = * ( id2oes->second );
+        _EdgesOnShape& eos = * ( id2eos->second );
         for ( size_t i = 0; i < eos._edges.size(); ++i )
         {
           _LayerEdge* ledge = eos._edges[ i ];
@@ -5012,10 +5028,10 @@ bool _ViscousBuilder::updateNormalsOfConvexFaces( _SolidData&         data,
     dumpFunction(SMESH_Comment("updateNormalsOfConvexFaces")<<data._index
                  <<"_F"<<meshDS->ShapeToIndex( convFace._face ));
 
-    id2oes = convFace._subIdToEOS.begin();
-    for ( ; id2oes != convFace._subIdToEOS.end(); ++id2oes )
+    id2eos = convFace._subIdToEOS.begin();
+    for ( ; id2eos != convFace._subIdToEOS.end(); ++id2eos )
     {
-      _EdgesOnShape& eos = * ( id2oes->second );
+      _EdgesOnShape& eos = * ( id2eos->second );
       for ( size_t i = 0; i < eos._edges.size(); ++i )
       {
         _LayerEdge* & ledge = eos._edges[ i ];
@@ -6600,7 +6616,8 @@ bool _ViscousBuilder::shrink()
     subEOS.clear();
     lEdges.clear();
     {
-      SMESH_subMeshIteratorPtr subIt = sm->getDependsOnIterator(/*includeSelf=*/false);
+      SMESH_subMeshIteratorPtr subIt = sm->getDependsOnIterator(/*includeSelf=*/false,
+                                                                /*complexFirst=*/true); //!!!
       while ( subIt->more() )
       {
         const TGeomID subID = subIt->next()->GetId();
@@ -6762,10 +6779,12 @@ bool _ViscousBuilder::shrink()
         int oldBadNb = badNb;
         badNb = 0;
         moved = false;
+        // '% 5' minimizes NB FUNCTIONS on viscous_layers_00/B2 case
+        _SmoothNode::SmoothType smooTy = ( smooStep % 5 ) ? smoothType : _SmoothNode::LAPLACIAN;
         for ( size_t i = 0; i < nodesToSmooth.size(); ++i )
         {
           moved |= nodesToSmooth[i].Smooth( badNb, surface, helper, refSign,
-                                            smoothType, /*set3D=*/isConcaveFace);
+                                            smooTy, /*set3D=*/isConcaveFace);
         }
         if ( badNb < oldBadNb )
           nbNoImpSteps = 0;