0020676: EDF 1212 GEOM: Partition operation creates vertices which causes mesh computation to fail with netgen

* Make NETGEN take into account internal vertices in faces and solids
* Simplify solution for internal faces

src/NETGENPlugin/NETGENPlugin_Mesher.hxx

@@ -31,6 +31,7 @@
 #include "NETGENPlugin_Defs.hxx"
 #include "StdMeshers_FaceSide.hxx"
 #include "SMDS_MeshElement.hxx"
+#include "SMESH_Algo.hxx"
 
 namespace nglib {
 #include <nglib.h>
@@ -45,6 +46,7 @@ class SMESH_Comment;
 class SMESHDS_Mesh;
 class TopoDS_Shape;
 class TopTools_DataMapOfShapeShape;
+class TopTools_IndexedMapOfShape;
 class NETGENPlugin_Hypothesis;
 class NETGENPlugin_SimpleHypothesis_2D;
 class NETGENPlugin_Internals;
@@ -98,11 +100,24 @@ class NETGENPLUGIN_EXPORT NETGENPlugin_Mesher
                        std::vector<const SMDS_MeshNode*>&  nodeVec,
                        SMESH_Comment&                      comment);
 
-  bool fillNgMesh(netgen::OCCGeometry&                occgeom,
+  bool fillNgMesh(const netgen::OCCGeometry&          occgeom,
                   netgen::Mesh&                       ngMesh,
                   std::vector<const SMDS_MeshNode*>&  nodeVec,
-                  const std::list< SMESH_subMesh* > & meshedSM,
+                  const std::list< SMESH_subMesh* > & meshedSM);
-                  NETGENPlugin_Internals*             internalShapes=0);
+
+  static void fixIntFaces(const netgen::OCCGeometry& occgeom,
+                          netgen::Mesh&              ngMesh,
+                          NETGENPlugin_Internals&    internalShapes);
+
+  static void addIntVerticesInFaces(const netgen::OCCGeometry&          occgeom,
+                                    netgen::Mesh&                       ngMesh,
+                                    std::vector<const SMDS_MeshNode*>&  nodeVec,
+                                    NETGENPlugin_Internals&             internalShapes);
+
+  static void addIntVerticesInSolids(const netgen::OCCGeometry&         occgeom,
+                                    netgen::Mesh&                       ngMesh,
+                                    std::vector<const SMDS_MeshNode*>&  nodeVec,
+                                    NETGENPlugin_Internals&             internalShapes);
 
   void defaultParameters();
 
@@ -133,6 +148,10 @@ class NETGENPLUGIN_EXPORT NETGENPlugin_Mesher
  * and their vertices shared by several internal edges. Nodes built on the found
  * shapes and mesh faces built on the found internal faces are to be doubled in
  * netgen mesh to emulate a "crack"
+ *
+ * For internal faces a more simple solution is found, which is just to duplicate
+ * mesh faces on internal geom faces without modeling a "real crack". For this
+ * reason findBorderElements() is no more used anywhere.
  */
 //=============================================================================
 
@@ -141,35 +160,48 @@ class NETGENPLUGIN_EXPORT NETGENPlugin_Internals
   SMESH_Mesh&       _mesh;
   bool              _is3D;
   //2D
-  std::map<int,int> _ev2face; //!< edges and vertices in faces where they are TopAbs_INTERNAL
+  std::map<int,int> _e2face;//!<edges and their vertices in faces where they are TopAbs_INTERNAL
+  std::map<int,std::list<int> > _f2v;//!<faces with internal vertices
   // 3D
   std::set<int>     _intShapes;
   std::set<int>     _borderFaces; //!< non-intrnal faces sharing the internal edge
+  std::map<int,std::list<int> > _s2v;//!<solids with internal vertices
 
 public:
   NETGENPlugin_Internals( SMESH_Mesh& mesh, const TopoDS_Shape& shape, bool is3D );
 
+  SMESH_Mesh& getMesh() const;
+
   bool isShapeToPrecompute(const TopoDS_Shape& s);
 
-  // 2D
+  // 2D meshing
-  bool hasInternalEdges() const { return !_ev2face.empty(); }
+  // edges 
-  bool isInternalEdge(int id ) const { return _ev2face.count( id ); }
+  bool hasInternalEdges() const { return !_e2face.empty(); }
-  bool isInternalVertex(int id ) const { return _ev2face.count( id ); }
+  bool isInternalEdge( int id ) const { return _e2face.count( id ); }
-  const std::map<int,int>& getEdgesAndVerticesWithFaces() const { return _ev2face; }
+  const std::map<int,int>& getEdgesAndVerticesWithFaces() const { return _e2face; }
   void getInternalEdges( TopTools_IndexedMapOfShape&  fmap,
                          TopTools_IndexedMapOfShape&  emap,
                          TopTools_IndexedMapOfShape&  vmap,
-                         list< SMESH_subMesh* >& smToPrecompute);
+                         std::list< SMESH_subMesh* >& smToPrecompute);
+  // vertices
+  bool hasInternalVertexInFace() const { return !_f2v.empty(); }
+  const std::map<int,std::list<int> >& getFacesWithVertices() const { return _f2v; }
 
-  // 3D
+  // 3D meshing
+  // faces
   bool hasInternalFaces() const { return !_intShapes.empty(); }
-  bool isInternalShape(int id ) const { return _intShapes.count( id ); }
+  bool isInternalShape( int id ) const { return _intShapes.count( id ); }
   void findBorderElements( std::set< const SMDS_MeshElement*, TIDCompare > & borderElems );
-  bool isBorderFace(int faceID ) const { return _borderFaces.count( faceID ); }
+  bool isBorderFace( int faceID ) const { return _borderFaces.count( faceID ); }
   void getInternalFaces( TopTools_IndexedMapOfShape&  fmap,
                          TopTools_IndexedMapOfShape&  emap,
-                         list< SMESH_subMesh* >&     facesSM,
+                         std::list< SMESH_subMesh* >& facesSM,
-                         list< SMESH_subMesh* >&     boundarySM);
+                         std::list< SMESH_subMesh* >& boundarySM);
+  // vertices
+  bool hasInternalVertexInSolid() const { return !_s2v.empty(); }
+  bool hasInternalVertexInSolid(int soID ) const { return _s2v.count(soID); }
+  const std::map<int,std::list<int> >& getSolidsWithVertices() const { return _s2v; }
+
 
 };
 

src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.cxx

@@ -317,68 +317,8 @@ static TError AddSegmentsToMesh(netgen::Mesh&                    ngMesh,
   } // loop on wires of a face
 
   // add a segment instead of internal vertex
-  // const TopoDS_Face& face = TopoDS::Face( helper.GetSubShape() );
+  NETGENPlugin_Internals intShapes( *helper.GetMesh(), helper.GetSubShape(), /*is3D=*/false );
-//   for ( TopoDS_Iterator sh (face); sh.More(); sh.Next())
+  NETGENPlugin_Mesher::addIntVerticesInFaces( geom, ngMesh, nodeVec, intShapes );
-//   {
-//     if ( sh.Value().ShapeType() != TopAbs_VERTEX ) continue;
-
-//     const TopoDS_Vertex V = TopoDS::Vertex( sh.Value() );
-//     SMESH_subMesh* sm = helper.GetMesh()->GetSubMesh( V );
-//     sm->ComputeStateEngine( SMESH_subMesh::COMPUTE );
-//     const SMDS_MeshNode * nV = SMESH_Algo::VertexNode( V, helper.GetMeshDS() );
-//     if ( !nV ) continue;
-//     double segLen = totalLength / ngMesh.GetNSeg() / 2;
-//     bool uvOK = false;
-//     gp_XY uvV = helper.GetNodeUV( face, nV, 0, &uvOK );
-//     if ( !uvOK ) helper.CheckNodeUV( face, nV, uvV, BRep_Tool::Tolerance( V ),/*force=*/true);
-//     gp_XY uvP( uvV.X() + segLen, uvV.Y() );
-//     TopLoc_Location loc;
-//     Handle(Geom_Surface) surf = BRep_Tool::Surface(face,loc);
-//     gp_Pnt P = surf->Value( uvP.X(), uvP.Y() ).Transformed( loc );
-
-//     MeshPoint mpV( Point<3> (nV->X(), nV->Y(), nV->Z()) );
-//     MeshPoint mpP( Point<3> (P.X(), P.Y(), P.Z()));
-
-//     ngMesh.AddPoint ( mpV, 1, EDGEPOINT );
-//     ngMesh.AddPoint ( mpP, 1, EDGEPOINT );
-
-//     nodeVec.push_back( nV );
-
-//     // Add the segment
-//     Segment seg;
-
-// #ifdef NETGEN_NEW
-//     seg.pnums[0] = ngMesh.GetNP()-1;  // ng node id
-//     seg.pnums[1] = ngMesh.GetNP();  // ng node id
-// #else
-//     seg.p1 = ngMesh.GetNP()-1;  // ng node id
-//     seg.p2 = ngMesh.GetNP(); // ng node id
-// #endif
-//     seg.edgenr = ngMesh.GetNSeg() + 1;// segment id
-//     seg.si = faceID;                  // = geom.fmap.FindIndex (face);
-
-//     seg.epgeominfo[ 0 ].dist = 0; // param on curve
-//     seg.epgeominfo[ 0 ].u    = uvV.X();
-//     seg.epgeominfo[ 0 ].v    = uvV.Y();
-//     seg.epgeominfo[ 1 ].dist = segLen; // param on curve
-//     seg.epgeominfo[ 1 ].u    = uvP.X();
-//     seg.epgeominfo[ 1 ].v    = uvP.Y();
-
-//     seg.epgeominfo[ 0 ].edgenr = 10; //  = geom.emap.FindIndex(edge);
-//     seg.epgeominfo[ 1 ].edgenr = 10; //  = geom.emap.FindIndex(edge);
-
-//     ngMesh.AddSegment (seg);
-
-//     // add reverse segment
-// #ifdef NETGEN_NEW
-//     swap (seg.pnums[0], seg.pnums[1]);
-// #else
-//     swap (seg.p1, seg.p2);
-// #endif
-//     swap( seg.epgeominfo[0], seg.epgeominfo[1] );
-//     seg.edgenr = ngMesh.GetNSeg() + 1; // segment id
-//     ngMesh.AddSegment (seg);
-//   }
 
   ngMesh.CalcSurfacesOfNode();
 

src/NETGENPlugin/NETGENPlugin_NETGEN_3D.cxx

@@ -63,6 +63,8 @@
   Netgen include files
 */
 
+#define OCCGEOMETRY
+#include <occgeom.hpp>
 namespace nglib {
 #include <nglib.h>
 }
@@ -108,10 +110,9 @@ NETGENPlugin_NETGEN_3D::~NETGENPlugin_NETGEN_3D()
  */
 //=============================================================================
 
-bool NETGENPlugin_NETGEN_3D::CheckHypothesis
+bool NETGENPlugin_NETGEN_3D::CheckHypothesis (SMESH_Mesh&         aMesh,
-                         (SMESH_Mesh& aMesh,
                                               const TopoDS_Shape& aShape,
-                          SMESH_Hypothesis::Hypothesis_Status& aStatus)
+                                              Hypothesis_Status&  aStatus)
 {
   MESSAGE("NETGENPlugin_NETGEN_3D::CheckHypothesis");
 
@@ -166,6 +167,7 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh&         aMesh,
 
   SMESH_MesherHelper helper(aMesh);
   bool _quadraticMesh = helper.IsQuadraticSubMesh(aShape);
+  helper.SetElementsOnShape( true );
 
   int Netgen_NbOfNodes     = 0;
   int Netgen_param2ndOrder = 0;
@@ -179,29 +181,22 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh&         aMesh,
   NETGENPlugin_NetgenLibWrapper ngLib;
   Ng_Mesh * Netgen_mesh = ngLib._ngMesh;
 
-  // maps of 1) ordinary nodes and 2) doubled nodes on internal shapes
+  // vector of nodes in which node index == netgen ID
-  typedef map< const SMDS_MeshNode*, int, TIDCompare > TNodeToIDMap;
+  vector< const SMDS_MeshNode* > nodeVec;
-  typedef TNodeToIDMap::value_type                     TN2ID;
-  TNodeToIDMap nodeToNetgenID[2];
-
   {
     const int invalid_ID = -1;
 
     SMESH::Controls::Area areaControl;
     SMESH::Controls::TSequenceOfXYZ nodesCoords;
 
-    // Issue 0020676 (StudyFiss_bugNetgen3D.hdf). Pb with internal face.
+    // maps nodes to ng ID
-    // Find internal geom faces, edges and vertices. 
+    typedef map< const SMDS_MeshNode*, int, TIDCompare > TNodeToIDMap;
-    // Nodes and faces built on the found internal shapes
+    typedef TNodeToIDMap::value_type                     TN2ID;
-    // will be doubled in Netgen input to make two borders of the "crack".
+    TNodeToIDMap nodeToNetgenID;
 
+    // find internal shapes
     NETGENPlugin_Internals internals( aMesh, aShape, /*is3D=*/true );
 
-    // mesh faces on non-internal geom faces sharing internal edge, whose some nodes
-    // are on internal edge and are to be replaced by doubled nodes
-    TIDSortedElemSet borderElems;
-    internals.findBorderElements( borderElems );
-
     // ---------------------------------
     // Feed the Netgen with surface mesh
     // ---------------------------------
@@ -218,7 +213,6 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh&         aMesh,
       const TopoDS_Shape& aShapeFace = exFa.Current();
       int faceID = meshDS->ShapeToIndex( aShapeFace );
       bool isInternalFace = internals.isInternalShape( faceID );
-      bool isBorderFace   = internals.isBorderFace( faceID );
       bool isRev = false;
       if ( checkReverse && !isInternalFace &&
            helper.NbAncestors(aShapeFace, aMesh, aShape.ShapeType()) > 1 )
@@ -252,15 +246,7 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh&         aMesh,
         }
         // Add nodes of triangles and triangles them-selves to netgen mesh
 
-        // a triangle on internal face is added twice,
-        // on border face, once but with doubled nodes
-        bool isBorder = ( isBorderFace && borderElems.count( elem ));
-        int nbDblLoops = ( isInternalFace && isTraingle || isBorder ) ? 2 : 1;
-
         for ( int i = 0; i < trias.size(); ++i )
-        {
-          bool reverse = isRev;
-          for ( int isDblF = isBorder; isDblF < nbDblLoops; ++isDblF, reverse = !reverse )
         {
           // add three nodes of triangle
           bool hasDegen = false;
@@ -276,8 +262,7 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh&         aMesh,
               node = SMESH_Algo::VertexNode( TopoDS::Vertex( vertex ), meshDS );
               hasDegen = true;
             }
-              bool isDblN = isDblF && internals.isInternalShape( shapeID );
+            int& ngID = nodeToNetgenID.insert(TN2ID( node, invalid_ID )).first->second;
-              int& ngID = nodeToNetgenID[isDblN].insert(TN2ID( node, invalid_ID )).first->second;
             if ( ngID == invalid_ID )
             {
               ngID = ++Netgen_NbOfNodes;
@@ -286,17 +271,20 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh&         aMesh,
               Netgen_point [ 2 ] = node->Z();
               Ng_AddPoint(Netgen_mesh, Netgen_point);
             }
-              Netgen_triangle[ iN ] = ngID;
+            Netgen_triangle[ isRev ? 3-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] ))
-              break;
+            continue;
-            if ( reverse )
-              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_
@@ -310,6 +298,22 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh&         aMesh,
       } // loop on elements on a face
     } // loop on faces of a SOLID or SHELL
 
+    // insert old nodes into nodeVec
+    nodeVec.resize( nodeToNetgenID.size() + 1, 0 );
+    TNodeToIDMap::iterator n_id = nodeToNetgenID.begin();
+    for ( ; n_id != nodeToNetgenID.end(); ++n_id )
+      nodeVec[ n_id->second ] = n_id->first;
+    nodeToNetgenID.clear();
+
+    if ( internals.hasInternalVertexInSolid() )
+    {
+      netgen::OCCGeometry occgeo;
+      NETGENPlugin_Mesher::addIntVerticesInSolids( occgeo,
+                                                   (netgen::Mesh&) *Netgen_mesh,
+                                                   nodeVec,
+                                                   internals);
+      Netgen_NbOfNodes = Ng_GetNP(Netgen_mesh);
+    }
   }
 
   // -------------------------
@@ -349,7 +353,6 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh&         aMesh,
   }
 
   int Netgen_NbOfNodesNew = Ng_GetNP(Netgen_mesh);
-
   int Netgen_NbOfTetra    = Ng_GetNE(Netgen_mesh);
 
   MESSAGE("End of Volume Mesh Generation. status=" << status <<
@@ -360,16 +363,6 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh&         aMesh,
   // Feed back the SMESHDS with the generated Nodes and Volume Elements
   // -------------------------------------------------------------------
 
-  // vector of nodes in which node index == netgen ID
-  vector< const SMDS_MeshNode* > nodeVec ( Netgen_NbOfNodesNew + 1 );
-  // insert old nodes into nodeVec
-  for ( int isDbl = 0; isDbl < 2; ++isDbl )
-  {
-    TNodeToIDMap::iterator n_id = nodeToNetgenID[isDbl].begin();
-    for ( ; n_id != nodeToNetgenID[isDbl].end(); ++n_id )
-      nodeVec[ n_id->second ] = n_id->first;
-    nodeToNetgenID[isDbl].clear();
-  }
   if ( status == NG_VOLUME_FAILURE )
   {
     SMESH_ComputeErrorPtr err = NETGENPlugin_Mesher::readErrors(nodeVec);
@@ -381,27 +374,22 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh&         aMesh,
   if ( isOK )
   {
     // create and insert new nodes into nodeVec
+    nodeVec.resize( Netgen_NbOfNodesNew + 1, 0 );
     int nodeIndex = Netgen_NbOfNodes + 1;
-    int shapeID = meshDS->ShapeToIndex( aShape );
     for ( ; nodeIndex <= Netgen_NbOfNodesNew; ++nodeIndex )
     {
       Ng_GetPoint( Netgen_mesh, nodeIndex, Netgen_point );
-      SMDS_MeshNode * node = meshDS->AddNode(Netgen_point[0],
+      nodeVec.at(nodeIndex) = helper.AddNode(Netgen_point[0], Netgen_point[1], Netgen_point[2]);
-                                             Netgen_point[1],
-                                             Netgen_point[2]);
-      meshDS->SetNodeInVolume(node, shapeID);
-      nodeVec.at(nodeIndex) = node;
     }
 
     // create tetrahedrons
     for ( int elemIndex = 1; elemIndex <= Netgen_NbOfTetra; ++elemIndex )
     {
       Ng_GetVolumeElement(Netgen_mesh, elemIndex, Netgen_tetrahedron);
-      SMDS_MeshVolume * elt = helper.AddVolume (nodeVec.at( Netgen_tetrahedron[0] ),
+      helper.AddVolume (nodeVec.at( Netgen_tetrahedron[0] ),
                         nodeVec.at( Netgen_tetrahedron[1] ),
                         nodeVec.at( Netgen_tetrahedron[2] ),
                         nodeVec.at( Netgen_tetrahedron[3] ));
-      meshDS->SetMeshElementOnShape(elt, shapeID );
     }
   }