0020832: EDF 1359 SMESH : Automatic meshing of boundary layers

for NETGEN_3D to work after StdMeshers_ViscousLayers

+               opt-hypos="ViscousLayers"

resources/NETGENPlugin.xml

@@ -55,6 +55,7 @@
                label-id="Tetrahedron (Netgen)"
                icon-id="mesh_algo_tetra.png"
                hypos="MaxElementVolume"
+               opt-hypos="ViscousLayers"
 	       need-geom="false"
                input="TRIA,QUAD"
                dim="3"/>
@@ -70,6 +71,7 @@
                label-id="Netgen 1D-2D-3D"
                icon-id="mesh_algo_netgen_2d3d.png"
                hypos="NETGEN_Parameters, NETGEN_SimpleParameters_3D"
+               opt-hypos="ViscousLayers"
                dim="3"
                support-submeshes="true"
                />

src/NETGENPlugin/NETGENPlugin_NETGEN_3D.cxx

@@ -39,9 +39,11 @@
 #include "SMESH_ControlsDef.hxx"
 #include "SMESH_Gen.hxx"
 #include "SMESH_Mesh.hxx"
-#include "SMESH_MesherHelper.hxx"
 #include "SMESH_MeshEditor.hxx"
+#include "SMESH_MesherHelper.hxx"
+#include "StdMeshers_MaxElementVolume.hxx"
 #include "StdMeshers_QuadToTriaAdaptor.hxx"
+#include "StdMeshers_ViscousLayers.hxx"
 
 #include <BRepGProp.hxx>
 #include <BRep_Tool.hxx>
@@ -78,14 +80,14 @@ using namespace std;
  */
 //=============================================================================
 
-NETGENPlugin_NETGEN_3D::NETGENPlugin_NETGEN_3D(int hypId, int studyId,
+NETGENPlugin_NETGEN_3D::NETGENPlugin_NETGEN_3D(int hypId, int studyId, SMESH_Gen* gen)
-                             SMESH_Gen* gen)
   : SMESH_3D_Algo(hypId, studyId, gen)
 {
   MESSAGE("NETGENPlugin_NETGEN_3D::NETGENPlugin_NETGEN_3D");
   _name = "NETGEN_3D";
   _shapeType = (1 << TopAbs_SHELL) | (1 << TopAbs_SOLID);// 1 bit /shape type
   _compatibleHypothesis.push_back("MaxElementVolume");
+  _compatibleHypothesis.push_back("ViscousLayers");
 
   _maxElementVolume = 0.;
 
@@ -118,39 +120,35 @@ bool NETGENPlugin_NETGEN_3D::CheckHypothesis (SMESH_Mesh&         aMesh,
   MESSAGE("NETGENPlugin_NETGEN_3D::CheckHypothesis");
 
   _hypMaxElementVolume = NULL;
+  _viscousLayersHyp    = NULL;
   _maxElementVolume = DBL_MAX;
 
-  list<const SMESHDS_Hypothesis*>::const_iterator itl;
+  const list<const SMESHDS_Hypothesis*>& hyps =
-  const SMESHDS_Hypothesis* theHyp;
+    GetUsedHypothesis(aMesh, aShape, /*ignoreAuxiliary=*/false);
-
+  list <const SMESHDS_Hypothesis* >::const_iterator h = hyps.begin();
-  const list<const SMESHDS_Hypothesis*>& hyps = GetUsedHypothesis(aMesh, aShape);
+  if ( h == hyps.end())
-  int nbHyp = hyps.size();
-  if (!nbHyp)
   {
     aStatus = SMESH_Hypothesis::HYP_OK;
-    //aStatus = SMESH_Hypothesis::HYP_MISSING;
     return true;  // can work with no hypothesis
   }
 
-  itl = hyps.begin();
+  aStatus = HYP_OK;
-  theHyp = (*itl); // use only the first hypothesis
+  for ( ; h != hyps.end(); ++h )
-
-  string hypName = theHyp->GetName();
-
-  bool isOk = false;
-
-  if (hypName == "MaxElementVolume")
   {
-    _hypMaxElementVolume = static_cast<const StdMeshers_MaxElementVolume*> (theHyp);
+    if ( !_hypMaxElementVolume )
-    ASSERT(_hypMaxElementVolume);
+      _hypMaxElementVolume = dynamic_cast< const StdMeshers_MaxElementVolume*> ( *h );
-    _maxElementVolume = _hypMaxElementVolume->GetMaxVolume();
+    if ( !_viscousLayersHyp )
-    isOk =true;
+      _viscousLayersHyp = dynamic_cast< const StdMeshers_ViscousLayers*> ( *h );
-    aStatus = SMESH_Hypothesis::HYP_OK;
-  }
-  else
-    aStatus = SMESH_Hypothesis::HYP_INCOMPATIBLE;
 
-  return isOk;
+    if ( *h != _hypMaxElementVolume &&
+         *h != _viscousLayersHyp )
+      aStatus = HYP_INCOMPATIBLE;
+  }
+
+  if ( _hypMaxElementVolume )
+    _maxElementVolume = _hypMaxElementVolume->GetMaxVolume();
+
+  return aStatus == HYP_OK;
 }
 
 //=============================================================================
@@ -187,9 +185,6 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh&         aMesh,
   {
     const int invalid_ID = -1;
 
-    SMESH::Controls::Area areaControl;
-    SMESH::Controls::TSequenceOfXYZ nodesCoords;
-
     // maps nodes to ng ID
     typedef map< const SMDS_MeshNode*, int, TIDCompare > TNodeToIDMap;
     typedef TNodeToIDMap::value_type                     TN2ID;
@@ -205,9 +200,19 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh&         aMesh,
     TopAbs_ShapeEnum mainType = aMesh.GetShapeToMesh().ShapeType();
     bool checkReverse = ( mainType == TopAbs_COMPOUND || mainType == TopAbs_COMPSOLID );
 
-    StdMeshers_QuadToTriaAdaptor Adaptor;
+    StdMeshers_ProxyMesh::Ptr proxyMesh( new StdMeshers_ProxyMesh( aMesh ));
+    if ( _viscousLayersHyp )
+    {
+      proxyMesh = _viscousLayersHyp->Compute( aMesh, aShape );
+      if ( !proxyMesh )
+        return false;
+    }
     if ( aMesh.NbQuadrangles() > 0 )
-      Adaptor.Compute(aMesh,aShape);
+    {
+      StdMeshers_QuadToTriaAdaptor* Adaptor = new StdMeshers_QuadToTriaAdaptor;
+      Adaptor->Compute(aMesh,aShape,proxyMesh.get());
+      proxyMesh.reset( Adaptor );
+    }
 
     for ( TopExp_Explorer exFa( aShape, TopAbs_FACE ); exFa.More(); exFa.Next())
     {
@@ -221,7 +226,7 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh&         aMesh,
         // so we use it as less as possible
         isRev = SMESH_Algo::IsReversedSubMesh( TopoDS::Face(aShapeFace), meshDS );
 
-      const SMESHDS_SubMesh * aSubMeshDSFace = meshDS->MeshElements( aShapeFace );
+      const SMESHDS_SubMesh * aSubMeshDSFace = proxyMesh->GetSubMesh( aShapeFace );
       if ( !aSubMeshDSFace ) continue;
       SMDS_ElemIteratorPtr iteratorElem = aSubMeshDSFace->GetElements();
       while ( iteratorElem->more() ) // loop on elements on a geom face
@@ -230,72 +235,49 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh&         aMesh,
         const SMDS_MeshElement* elem = iteratorElem->next();
         if ( !elem )
           return error( COMPERR_BAD_INPUT_MESH, "Null element encounters");
-        vector< const SMDS_MeshElement* > trias;
+        if ( elem->NbCornerNodes() != 3 )
-        bool isTraingle = ( elem->NbNodes() == ( elem->IsQuadratic() ? 6 : 3 ));
+          return error( COMPERR_BAD_INPUT_MESH, "Not triangle element encounters");
-        if ( !isTraingle )
+
-        {
-          // use adaptor to convert quadrangle face into triangles
-          const list<const SMDS_MeshFace*>* faces = Adaptor.GetTriangles(elem);
-          if(faces==0)
-            return error( COMPERR_BAD_INPUT_MESH,
-                          SMESH_Comment("No triangles in adaptor for element ")<<elem->GetID());
-          trias.assign( faces->begin(), faces->end() );
-        }
-        else
-        {
-          trias.push_back( elem );
-        }
         // Add nodes of triangles and triangles them-selves to netgen mesh
 
-        for ( int i = 0; i < trias.size(); ++i )
+        // add three nodes of triangle
+        bool hasDegen = false;
+        for ( int iN = 0; iN < 3; ++iN )
         {
-          // add three nodes of triangle
+          const SMDS_MeshNode* node = elem->GetNode( iN );
-          bool hasDegen = false;
+          const int shapeID = node->GetPosition()->GetShapeId();
-          for ( int iN = 0; iN < 3; ++iN )
+          if ( node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_EDGE &&
+               helper.IsDegenShape( shapeID ))
           {
-            const SMDS_MeshNode* node = trias[i]->GetNode( iN );
+            // ignore all nodes on degeneraged edge and use node on its vertex instead
-            int shapeID = node->GetPosition()->GetShapeId();
+            TopoDS_Shape vertex = TopoDS_Iterator( meshDS->IndexToShape( shapeID )).Value();
-            if ( node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_EDGE &&
+            node = SMESH_Algo::VertexNode( TopoDS::Vertex( vertex ), meshDS );
-                 helper.IsDegenShape( shapeID ))
+            hasDegen = true;
-            {
-              // ignore all nodes on degeneraged edge and use node on its vertex instead
-              TopoDS_Shape vertex = TopoDS_Iterator( meshDS->IndexToShape( shapeID )).Value();
-              node = SMESH_Algo::VertexNode( TopoDS::Vertex( vertex ), meshDS );
-              hasDegen = true;
-            }
-            int& ngID = nodeToNetgenID.insert(TN2ID( node, invalid_ID )).first->second;
-            if ( ngID == invalid_ID )
-            {
-              ngID = ++Netgen_NbOfNodes;
-              Netgen_point [ 0 ] = node->X();
-              Netgen_point [ 1 ] = node->Y();
-              Netgen_point [ 2 ] = node->Z();
-              Ng_AddPoint(Netgen_mesh, Netgen_point);
-            }
-            Netgen_triangle[ isRev ? 2-iN : iN ] = ngID;
           }
-          // add triangle
+          int& ngID = nodeToNetgenID.insert(TN2ID( node, invalid_ID )).first->second;
-          if ( hasDegen && (Netgen_triangle[0] == Netgen_triangle[1] ||
+          if ( ngID == invalid_ID )
-                            Netgen_triangle[0] == Netgen_triangle[2] ||
+          {
-                            Netgen_triangle[2] == Netgen_triangle[1] ))
+            ngID = ++Netgen_NbOfNodes;
-            continue;
+            Netgen_point [ 0 ] = node->X();
+            Netgen_point [ 1 ] = node->Y();
+            Netgen_point [ 2 ] = node->Z();
+            Ng_AddPoint(Netgen_mesh, Netgen_point);
+          }
+          Netgen_triangle[ isRev ? 2-iN : iN ] = ngID;
+        }
+        // add triangle
+        if ( hasDegen && (Netgen_triangle[0] == Netgen_triangle[1] ||
+                          Netgen_triangle[0] == Netgen_triangle[2] ||
+                          Netgen_triangle[2] == Netgen_triangle[1] ))
+          continue;
 
+        Ng_AddSurfaceElement(Netgen_mesh, NG_TRIG, Netgen_triangle);
+
+        if ( isInternalFace && !proxyMesh->IsTemporary( elem ))
+        {
+          swap( Netgen_triangle[1], Netgen_triangle[2] );
           Ng_AddSurfaceElement(Netgen_mesh, NG_TRIG, Netgen_triangle);
-
-          if ( isInternalFace && isTraingle )
-          {
-            swap( Netgen_triangle[1], Netgen_triangle[2] );
-            Ng_AddSurfaceElement(Netgen_mesh, NG_TRIG, Netgen_triangle);
-          }
         }
-#ifdef _DEBUG_
-        // check if a trainge is degenerated
-        areaControl.GetPoints( elem, nodesCoords );
-        double area = areaControl.GetValue( nodesCoords );
-        if ( area <= DBL_MIN ) {
-          MESSAGE( "Warning: Degenerated " << elem );
-        }
-#endif
       } // loop on elements on a face
     } // loop on faces of a SOLID or SHELL
 
@@ -409,54 +391,14 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh&         aMesh,
   MESSAGE("NETGENPlugin_NETGEN_3D::Compute with maxElmentsize = " << _maxElementVolume);  
   const int invalid_ID = -1;
   bool _quadraticMesh = false;
-  typedef map< const SMDS_MeshNode*, int, TIDCompare > TNodeToIDMap;
-  TNodeToIDMap nodeToNetgenID;
-  list< const SMDS_MeshElement* > triangles;
-  SMESHDS_Mesh* MeshDS = aHelper->GetMeshDS();
 
   SMESH_MesherHelper::MType MeshType = aHelper->IsQuadraticMesh();
-  
+
   if(MeshType == SMESH_MesherHelper::COMP)
     return error( COMPERR_BAD_INPUT_MESH,
                   SMESH_Comment("Mesh with linear and quadratic elements given."));
   else if (MeshType == SMESH_MesherHelper::QUADRATIC)
     _quadraticMesh = true;
-    
-  StdMeshers_QuadToTriaAdaptor Adaptor;
-  if ( aMesh.NbQuadrangles() > 0 )
-    Adaptor.Compute(aMesh);
-
-  SMDS_FaceIteratorPtr fIt = MeshDS->facesIterator(/*idInceasingOrder=*/true);
-  while( fIt->more())
-  {
-    // check element
-    const SMDS_MeshElement* elem = fIt->next();
-    if ( !elem )
-      return error( COMPERR_BAD_INPUT_MESH, "Null element encounters");
-
-    vector< const SMDS_MeshElement* > trias;
-    bool isTraingle = ( elem->NbCornerNodes() == 3 );
-    if ( !isTraingle ) {
-      // using adaptor
-      const list<const SMDS_MeshFace*>* faces = Adaptor.GetTriangles(elem);
-      if(faces==0)
-        continue; // Issue 0020682. There already can be 3d mesh
-      trias.assign( faces->begin(), faces->end() );
-    }
-    else {
-      trias.push_back( elem );
-    }
-    for ( int i = 0; i < trias.size(); ++i )
-    {
-      triangles.push_back( trias[i] );
-      for ( int iN = 0; iN < 3; ++iN )
-      {
-        const SMDS_MeshNode* node = trias[i]->GetNode( iN );
-        // put elem nodes to nodeToNetgenID map
-        nodeToNetgenID.insert( make_pair( node, invalid_ID ));
-      }
-    }
-  }
 
   // ---------------------------------
   // Feed the Netgen with surface mesh
@@ -474,41 +416,53 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh&         aMesh,
   NETGENPlugin_NetgenLibWrapper ngLib;
   Ng_Mesh * Netgen_mesh = ngLib._ngMesh;
 
-  // set nodes and remember thier netgen IDs
+  StdMeshers_ProxyMesh::Ptr proxyMesh( new StdMeshers_ProxyMesh( aMesh ));
-  
+  if ( aMesh.NbQuadrangles() > 0 )
-  TNodeToIDMap::iterator n_id = nodeToNetgenID.begin();
-  for ( ; n_id != nodeToNetgenID.end(); ++n_id )
   {
-    const SMDS_MeshNode* node = n_id->first;
+    StdMeshers_QuadToTriaAdaptor* Adaptor = new StdMeshers_QuadToTriaAdaptor;
-
+    Adaptor->Compute(aMesh);
-    Netgen_point [ 0 ] = node->X();
+    proxyMesh.reset( Adaptor );
-    Netgen_point [ 1 ] = node->Y();
-    Netgen_point [ 2 ] = node->Z();
-    Ng_AddPoint(Netgen_mesh, Netgen_point);
-    n_id->second = ++Netgen_NbOfNodes; // set netgen ID
   }
 
-  // set triangles
+  // maps nodes to ng ID
-  list< const SMDS_MeshElement* >::iterator tria = triangles.begin();
+  typedef map< const SMDS_MeshNode*, int, TIDCompare > TNodeToIDMap;
-  for ( ; tria != triangles.end(); ++tria)
+  typedef TNodeToIDMap::value_type                     TN2ID;
+  TNodeToIDMap nodeToNetgenID;
+
+  SMDS_ElemIteratorPtr fIt = proxyMesh->GetFaces();
+  while( fIt->more())
   {
-    int i = 0;
+    // check element
-    SMDS_ElemIteratorPtr triangleNodesIt = (*tria)->nodesIterator();
+    const SMDS_MeshElement* elem = fIt->next();
-    while ( triangleNodesIt->more() ) {
+    if ( !elem )
-      const SMDS_MeshNode * node =
+      return error( COMPERR_BAD_INPUT_MESH, "Null element encounters");
-        static_cast<const SMDS_MeshNode *>(triangleNodesIt->next());
+    if ( elem->NbCornerNodes() != 3 )
-      if(aHelper->IsMedium(node))
+      return error( COMPERR_BAD_INPUT_MESH, "Not triangle element encounters");
-        continue;
+      
-      Netgen_triangle[ i ] = nodeToNetgenID[ node ];
+    // add three nodes of triangle
-      ++i;
+    for ( int iN = 0; iN < 3; ++iN )
+    {
+      const SMDS_MeshNode* node = elem->GetNode( iN );
+      int& ngID = nodeToNetgenID.insert(TN2ID( node, invalid_ID )).first->second;
+      if ( ngID == invalid_ID )
+      {
+        ngID = ++Netgen_NbOfNodes;
+        Netgen_point [ 0 ] = node->X();
+        Netgen_point [ 1 ] = node->Y();
+        Netgen_point [ 2 ] = node->Z();
+        Ng_AddPoint(Netgen_mesh, Netgen_point);
+      }
+      Netgen_triangle[ iN ] = ngID;
     }
     Ng_AddSurfaceElement(Netgen_mesh, NG_TRIG, Netgen_triangle);
   }
+  proxyMesh.reset(); // delete tmp faces
 
   // vector of nodes in which node index == netgen ID
   vector< const SMDS_MeshNode* > nodeVec ( nodeToNetgenID.size() + 1 );
   // insert old nodes into nodeVec
-  for ( n_id = nodeToNetgenID.begin(); n_id != nodeToNetgenID.end(); ++n_id )
+  TNodeToIDMap::iterator n_id = nodeToNetgenID.begin();
+  for ( ; n_id != nodeToNetgenID.end(); ++n_id )
     nodeVec.at( n_id->second ) = n_id->first;
   nodeToNetgenID.clear();
 

src/NETGENPlugin/NETGENPlugin_NETGEN_3D.hxx

@@ -26,7 +26,6 @@
 // Created   : lundi 27 Janvier 2003
 // Author    : Nadir BOUHAMOU (CEA)
 // Project   : SALOME
-// $Header$
 //=============================================================================
 //
 #ifndef _NETGENPlugin_NETGEN_3D_HXX_
@@ -35,10 +34,11 @@
 #include "NETGENPlugin_Defs.hxx"
 
 #include "SMESH_3D_Algo.hxx"
-#include "SMESH_Mesh.hxx"
-#include "StdMeshers_MaxElementVolume.hxx"
 #include "Utils_SALOME_Exception.hxx"
 
+class StdMeshers_ViscousLayers;
+class StdMeshers_MaxElementVolume;
+
 class NETGENPLUGIN_EXPORT NETGENPlugin_NETGEN_3D: public SMESH_3D_Algo
 {
 public:
@@ -63,6 +63,7 @@ protected:
   double _maxElementVolume;
 
   const StdMeshers_MaxElementVolume* _hypMaxElementVolume;
+  const StdMeshers_ViscousLayers*    _viscousLayersHyp;
 };
 
 #endif