]> SALOME platform Git repositories - plugins/netgenplugin.git/commitdiff
Salome HOME
Using RestrictLocalH netgen function to simulate feed of 1D elements for 2D_Only...
authorcconopoima <cesar.conopoima@gmail.com>
Tue, 30 Apr 2024 11:20:03 +0000 (12:20 +0100)
committercconopoima <cesar.conopoima@gmail.com>
Tue, 30 Apr 2024 11:20:03 +0000 (12:20 +0100)
src/NETGENPlugin/NETGENPlugin_Mesher.cxx
src/NETGENPlugin/NETGENPlugin_Mesher.hxx
src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.cxx
src/NETGENPlugin/NETGENPlugin_Remesher_2D.cxx

index efa5e70c4ed55632742b33c0d389bdf9becd6cd9..7a49e1ff733f338d73e7e6035d2e6be885ea814c 100644 (file)
@@ -52,6 +52,8 @@
 
 #include <utilities.h>
 
+#include <GEOMUtils.hxx>
+
 #include <BRepAdaptor_Surface.hxx>
 #include <BRepBuilderAPI_Copy.hxx>
 #include <BRepLProp_SLProps.hxx>
@@ -957,7 +959,8 @@ void NETGENPlugin_Mesher::PrepareOCCgeometry(netgen::OCCGeometry&     occgeo,
   updateTriangulation( shape );
 
   Bnd_Box bb;
-  BRepBndLib::Add (shape, bb);
+  // BRepBndLib::Add (shape, bb);
+  GEOMUtils::PreciseBoundingBox(shape,bb);
   double x1,y1,z1,x2,y2,z2;
   bb.Get (x1,y1,z1,x2,y2,z2);
   netgen::Point<3> p1 = netgen::Point<3> (x1,y1,z1);
@@ -1026,6 +1029,9 @@ void NETGENPlugin_Mesher::PrepareOCCgeometry(netgen::OCCGeometry&     occgeo,
   occgeo.face_maxh_modified = 0;
   occgeo.face_maxh.SetSize(totNbFaces);
   occgeo.face_maxh = netgen::mparam.maxh;
+  
+  //
+  occgeo.BuildFMapFromPrefilled();
 }
 
 //================================================================================
@@ -1153,10 +1159,11 @@ bool NETGENPlugin_Mesher::FillNgMesh(netgen::OCCGeometry&           occgeom,
     const SMESHDS_SubMesh * smDS = sm->GetSubMeshDS();
     if ( !smDS ) continue;
 
+    int countEdgeId = 0;
     switch ( sm->GetSubShape().ShapeType() )
     {
     case TopAbs_EDGE: { // EDGE
-      // ----------------------
+
       TopoDS_Edge geomEdge  = TopoDS::Edge( sm->GetSubShape() );
       if ( geomEdge.Orientation() >= TopAbs_INTERNAL )
         geomEdge.Orientation( TopAbs_FORWARD ); // issue 0020676
@@ -1180,7 +1187,7 @@ bool NETGENPlugin_Mesher::FillNgMesh(netgen::OCCGeometry&           occgeom,
         // get all meshed EDGEs of the FACE connected to geomEdge (issue 0021140)
         helper.SetSubShape( face );
         list< TopoDS_Edge > edges = getConnectedEdges( geomEdge, face, computedSM, helper,
-                                                       visitedEdgeSM2Faces );
+                                                      visitedEdgeSM2Faces );
         if ( edges.empty() )
           continue; // wrong ancestor?
 
@@ -1203,119 +1210,191 @@ bool NETGENPlugin_Mesher::FillNgMesh(netgen::OCCGeometry&           occgeom,
         if ( points.empty() )
           return false; // invalid node params?
         smIdType i, nbSeg = fSide.NbSegments();
+        for ( i = 0; i < nbSeg; ++i )
+        {
+          const SMDS_MeshNode * n0 = points[ i ].node;
+          const SMDS_MeshNode * n1 = points[ i+1 ].node;
+           double size = SMESH_TNodeXYZ( n0 ).Distance( n1 );
+          netgen::Point3d p0(n0->X(), n0->Y(), n0->Z());
+          netgen::Point3d p1(n1->X(), n1->Y(), n1->Z());
+          ngMesh.RestrictLocalHLine( p0, p1, size, true  );
+        }
 
-        // remember EDGEs of fSide to treat only once
-        for ( int iE = 0; iE < fSide.NbEdges(); ++iE )
-          visitedEdgeSM2Faces[ helper.GetMesh()->GetSubMesh( fSide.Edge(iE )) ].insert(faceSMDSId);
+      }
+///////////////////////////////////////////////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////////////////////////////////////////////
+      {      
+        /*
+          // ----------------------
+          TopoDS_Edge geomEdge  = TopoDS::Edge( sm->GetSubShape() );
+          if ( geomEdge.Orientation() >= TopAbs_INTERNAL )
+            geomEdge.Orientation( TopAbs_FORWARD ); // issue 0020676
+
+          // Add ng segments for each not meshed FACE the EDGE bounds
+          PShapeIteratorPtr fIt = helper.GetAncestors( geomEdge, *sm->GetFather(), TopAbs_FACE );
+          while ( const TopoDS_Shape * anc = fIt->next() )
+          {
+            faceNgID = occgeom.fmap.FindIndex( *anc );
+            if ( faceNgID < 1 )
+              continue; // meshed face
+
+            int faceSMDSId = meshDS->ShapeToIndex( *anc );
+            if ( visitedEdgeSM2Faces[ sm ].count( faceSMDSId ))
+              continue; // already treated EDGE
+
+            TopoDS_Face face = TopoDS::Face( occgeom.fmap( faceNgID ));
+            if ( face.Orientation() >= TopAbs_INTERNAL )
+              face.Orientation( TopAbs_FORWARD ); // issue 0020676
+
+            // get all meshed EDGEs of the FACE connected to geomEdge (issue 0021140)
+            helper.SetSubShape( face );
+            list< TopoDS_Edge > edges = getConnectedEdges( geomEdge, face, computedSM, helper,
+                                                          visitedEdgeSM2Faces );
+            if ( edges.empty() )
+              continue; // wrong ancestor?
+
+            // find out orientation of <edges> within <face>
+            TopoDS_Edge eNotSeam = edges.front();
+            if ( helper.HasSeam() )
+            {
+              list< TopoDS_Edge >::iterator eIt = edges.begin();
+              while ( helper.IsRealSeam( *eIt )) ++eIt;
+              if ( eIt != edges.end() )
+                eNotSeam = *eIt;
+            }
+            TopAbs_Orientation fOri = helper.GetSubShapeOri( face, eNotSeam );
+            bool isForwad = ( fOri == eNotSeam.Orientation() || fOri >= TopAbs_INTERNAL );
 
-        double otherSeamParam = 0;
-        bool isSeam = false;
+            // get all nodes from connected <edges>
+            const bool skipMedium = netgen::mparam.secondorder;//smDS->IsQuadratic();
+            StdMeshers_FaceSide fSide( face, edges, _mesh, isForwad, skipMedium, &helper );
+            const vector<UVPtStruct>& points = fSide.GetUVPtStruct();
+            if ( points.empty() )
+              return false; // invalid node params?
+            smIdType i, nbSeg = fSide.NbSegments();
 
-        // add segments
+            // remember EDGEs of fSide to treat only once
+            for ( int iE = 0; iE < fSide.NbEdges(); ++iE )
+              visitedEdgeSM2Faces[ helper.GetMesh()->GetSubMesh( fSide.Edge(iE )) ].insert(faceSMDSId);
 
-        int prevNgId = ngNodeId( points[0].node, ngMesh, nodeNgIdMap );
+            double otherSeamParam = 0;
+            bool isSeam = false;
 
-        for ( i = 0; i < nbSeg; ++i )
-        {
-          const UVPtStruct& p1 = points[ i ];
-          const UVPtStruct& p2 = points[ i+1 ];
+            // add segments
 
-          if ( p1.node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_VERTEX ) //an EDGE begins
-          {
-            isSeam = false;
-            if ( helper.IsRealSeam( p1.node->GetShapeID() ))
+            int prevNgId = ngNodeId( points[0].node, ngMesh, nodeNgIdMap );
+
+            for ( i = 0; i < nbSeg; ++i )
             {
-              TopoDS_Edge e = fSide.Edge( fSide.EdgeIndex( 0.5 * ( p1.normParam + p2.normParam )));
-              isSeam = helper.IsRealSeam( e );
+              const UVPtStruct& p1 = points[ i ];
+              const UVPtStruct& p2 = points[ i+1 ];
+
+              if ( p1.node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_VERTEX ) //an EDGE begins
+              {
+                isSeam = false;
+                if ( helper.IsRealSeam( p1.node->GetShapeID() ))
+                {
+                  TopoDS_Edge e = fSide.Edge( fSide.EdgeIndex( 0.5 * ( p1.normParam + p2.normParam )));
+                  isSeam = helper.IsRealSeam( e );
+                  if ( isSeam )
+                  {
+                    otherSeamParam = helper.GetOtherParam( helper.GetPeriodicIndex() & 1 ? p2.u : p2.v );
+                  }
+                }
+              }
+              netgen::Segment seg;
+              // ng node ids
+              seg[0] = prevNgId;
+              seg[1] = prevNgId = ngNodeId( p2.node, ngMesh, nodeNgIdMap );
+              seg.edgenr  = countEdgeId+1;
+              seg.si      = countEdgeId+1;
+              seg.epgeominfo[0].edgenr = countEdgeId;
+              seg.epgeominfo[1].edgenr = countEdgeId;
+              // node param on curve
+              seg.epgeominfo[ 0 ].dist = p1.param;
+              seg.epgeominfo[ 1 ].dist = p2.param;
+              // uv on face
+              seg.epgeominfo[ 0 ].u = p1.u;
+              seg.epgeominfo[ 0 ].v = p1.v;
+              seg.epgeominfo[ 1 ].u = p2.u;
+              seg.epgeominfo[ 1 ].v = p2.v;
+
+              //geomEdge = fSide.Edge( fSide.EdgeIndex( 0.5 * ( p1.normParam + p2.normParam )));
+              //seg.epgeominfo[ 0 ].edgenr = seg.epgeominfo[ 1 ].edgenr = occgeom.emap.FindIndex( geomEdge );
+
+              //seg.epgeominfo[ iEnd ].edgenr = edgeID; //  = geom.emap.FindIndex(edge);
+              // seg.si = faceNgID;                   // = geom.fmap.FindIndex (face);
+              // seg.edgenr = ngMesh.GetNSeg() + 1; // segment id
+              ngMesh.AddSegment (seg);
+
+              SMESH_TNodeXYZ np1( p1.node ), np2( p2.node );
+              RestrictLocalSize( ngMesh, 0.5*(np1+np2), (np1-np2).Modulus() );
+
+    #ifdef DUMP_SEGMENTS
+              cout << "Segment: " << seg.edgenr << " on SMESH face " << meshDS->ShapeToIndex( face ) << endl
+                  << "\tface index: " << seg.si << endl
+                  << "\tp1: " << seg[0] << endl
+                  << "\tp2: " << seg[1] << endl
+                  << "\tp0 param: " << seg.epgeominfo[ 0 ].dist << endl
+                  << "\tp0 uv: " << seg.epgeominfo[ 0 ].u <<", "<< seg.epgeominfo[ 0 ].v << endl
+                //<< "\tp0 edge: " << seg.epgeominfo[ 0 ].edgenr << endl
+                  << "\tp1 param: " << seg.epgeominfo[ 1 ].dist << endl
+                  << "\tp1 uv: " << seg.epgeominfo[ 1 ].u <<", "<< seg.epgeominfo[ 1 ].v << endl;
+                //<< "\tp1 edge: " << seg.epgeominfo[ 1 ].edgenr << endl;
+    #endif
               if ( isSeam )
               {
-                otherSeamParam = helper.GetOtherParam( helper.GetPeriodicIndex() & 1 ? p2.u : p2.v );
+                if ( helper.GetPeriodicIndex() && 1 ) {
+                  seg.epgeominfo[ 0 ].u = otherSeamParam;
+                  seg.epgeominfo[ 1 ].u = otherSeamParam;
+                  swap (seg.epgeominfo[0].v, seg.epgeominfo[1].v);
+                } else {
+                  seg.epgeominfo[ 0 ].v = otherSeamParam;
+                  seg.epgeominfo[ 1 ].v = otherSeamParam;
+                  swap (seg.epgeominfo[0].u, seg.epgeominfo[1].u);
+                }
+                swap( seg[0], seg[1] );
+                swap( seg.epgeominfo[0].dist, seg.epgeominfo[1].dist );
+                seg.edgenr = ngMesh.GetNSeg() + 1; // segment id
+                ngMesh.AddSegment( seg );
+    #ifdef DUMP_SEGMENTS
+                cout << "Segment: " << seg.edgenr << endl
+                    << "\t is SEAM (reverse) of the previous. "
+                    << " Other " << (helper.GetPeriodicIndex() && 1 ? "U" : "V")
+                    << " = " << otherSeamParam << endl;
+    #endif
+              }
+              else if ( fOri == TopAbs_INTERNAL )
+              {
+                swap( seg[0], seg[1] );
+                swap( seg.epgeominfo[0], seg.epgeominfo[1] );
+                seg.edgenr = ngMesh.GetNSeg() + 1; // segment id
+                ngMesh.AddSegment( seg );
+    #ifdef DUMP_SEGMENTS
+                cout << "Segment: " << seg.edgenr << endl << "\t is REVERSE of the previous" << endl;
+    #endif
               }
             }
-          }
-          netgen::Segment seg;
-          // ng node ids
-          seg[0] = prevNgId;
-          seg[1] = prevNgId = ngNodeId( p2.node, ngMesh, nodeNgIdMap );
-          // node param on curve
-          seg.epgeominfo[ 0 ].dist = p1.param;
-          seg.epgeominfo[ 1 ].dist = p2.param;
-          // uv on face
-          seg.epgeominfo[ 0 ].u = p1.u;
-          seg.epgeominfo[ 0 ].v = p1.v;
-          seg.epgeominfo[ 1 ].u = p2.u;
-          seg.epgeominfo[ 1 ].v = p2.v;
-
-          //geomEdge = fSide.Edge( fSide.EdgeIndex( 0.5 * ( p1.normParam + p2.normParam )));
-          //seg.epgeominfo[ 0 ].edgenr = seg.epgeominfo[ 1 ].edgenr = occgeom.emap.FindIndex( geomEdge );
-
-          //seg.epgeominfo[ iEnd ].edgenr = edgeID; //  = geom.emap.FindIndex(edge);
-          seg.si = faceNgID;                   // = geom.fmap.FindIndex (face);
-          seg.edgenr = ngMesh.GetNSeg() + 1; // segment id
-          ngMesh.AddSegment (seg);
-
-          SMESH_TNodeXYZ np1( p1.node ), np2( p2.node );
-          RestrictLocalSize( ngMesh, 0.5*(np1+np2), (np1-np2).Modulus() );
+          } // loop on geomEdge ancestors
 
-#ifdef DUMP_SEGMENTS
-          cout << "Segment: " << seg.edgenr << " on SMESH face " << meshDS->ShapeToIndex( face ) << endl
-               << "\tface index: " << seg.si << endl
-               << "\tp1: " << seg[0] << endl
-               << "\tp2: " << seg[1] << endl
-               << "\tp0 param: " << seg.epgeominfo[ 0 ].dist << endl
-               << "\tp0 uv: " << seg.epgeominfo[ 0 ].u <<", "<< seg.epgeominfo[ 0 ].v << endl
-            //<< "\tp0 edge: " << seg.epgeominfo[ 0 ].edgenr << endl
-               << "\tp1 param: " << seg.epgeominfo[ 1 ].dist << endl
-               << "\tp1 uv: " << seg.epgeominfo[ 1 ].u <<", "<< seg.epgeominfo[ 1 ].v << endl;
-            //<< "\tp1 edge: " << seg.epgeominfo[ 1 ].edgenr << endl;
-#endif
-          if ( isSeam )
+          if ( quadHelper ) // remember medium nodes of sub-meshes
           {
-            if ( helper.GetPeriodicIndex() && 1 ) {
-              seg.epgeominfo[ 0 ].u = otherSeamParam;
-              seg.epgeominfo[ 1 ].u = otherSeamParam;
-              swap (seg.epgeominfo[0].v, seg.epgeominfo[1].v);
-            } else {
-              seg.epgeominfo[ 0 ].v = otherSeamParam;
-              seg.epgeominfo[ 1 ].v = otherSeamParam;
-              swap (seg.epgeominfo[0].u, seg.epgeominfo[1].u);
+            SMDS_ElemIteratorPtr edges = smDS->GetElements();
+            while ( edges->more() )
+            {
+              const SMDS_MeshElement* e = edges->next();
+              if ( !quadHelper->AddTLinks( static_cast< const SMDS_MeshEdge*>( e )))
+                break;
             }
-            swap( seg[0], seg[1] );
-            swap( seg.epgeominfo[0].dist, seg.epgeominfo[1].dist );
-            seg.edgenr = ngMesh.GetNSeg() + 1; // segment id
-            ngMesh.AddSegment( seg );
-#ifdef DUMP_SEGMENTS
-            cout << "Segment: " << seg.edgenr << endl
-                 << "\t is SEAM (reverse) of the previous. "
-                 << " Other " << (helper.GetPeriodicIndex() && 1 ? "U" : "V")
-                 << " = " << otherSeamParam << endl;
-#endif
           }
-          else if ( fOri == TopAbs_INTERNAL )
-          {
-            swap( seg[0], seg[1] );
-            swap( seg.epgeominfo[0], seg.epgeominfo[1] );
-            seg.edgenr = ngMesh.GetNSeg() + 1; // segment id
-            ngMesh.AddSegment( seg );
-#ifdef DUMP_SEGMENTS
-            cout << "Segment: " << seg.edgenr << endl << "\t is REVERSE of the previous" << endl;
-#endif
-          }
-        }
-      } // loop on geomEdge ancestors
-
-      if ( quadHelper ) // remember medium nodes of sub-meshes
-      {
-        SMDS_ElemIteratorPtr edges = smDS->GetElements();
-        while ( edges->more() )
-        {
-          const SMDS_MeshElement* e = edges->next();
-          if ( !quadHelper->AddTLinks( static_cast< const SMDS_MeshEdge*>( e )))
-            break;
-        }
-      }
+          
+          countEdgeId++; // Assume edge is iterate incrementally as done by the emap!
+          break;
+        
+        */  
 
-      break;
+      } /* scope of original code*/
     } // case TopAbs_EDGE
 
     case TopAbs_FACE: { // FACE
@@ -1439,6 +1518,7 @@ bool NETGENPlugin_Mesher::FillNgMesh(netgen::OCCGeometry&           occgeom,
       int fID = occgeom.fmap.Add( geomFace );
       if ( occgeom.facemeshstatus.Size() < fID ) occgeom.facemeshstatus.SetSize( fID );
       occgeom.facemeshstatus[ fID-1 ] = netgen::FACE_MESHED_OK;
+
       while ( fID < faceNgID ) // geomFace is already in occgeom.fmap, add a copy
       {
         fID = occgeom.fmap.Add( BRepBuilderAPI_Copy( geomFace, /*copyGeom=*/false ));
@@ -2219,6 +2299,61 @@ void NETGENPlugin_Mesher::AddIntVerticesInSolids(const netgen::OCCGeometry&
   } // loop on solids with internal vertices
 }
 
+
+// Iterate in all edges of the geometry, and use RestrictLocalHLine function to set the value of allow edges in 1D mesh
+SMESH_ComputeErrorPtr
+NETGENPlugin_Mesher::AddSegmentsToMesh2(netgen::Mesh&                    ngMesh,
+                                        netgen::OCCGeometry&             geom,
+                                        const TSideVector&               wires,
+                                        SMESH_MesherHelper&              helper,
+                                        vector< const SMDS_MeshNode* > & nodeVec,
+                                        const bool                       overrideMinH)
+{
+  map<const SMDS_MeshNode*, int > node2ngID;
+
+  for ( size_t iW = 0; iW < wires.size(); ++iW )
+  {
+    StdMeshers_FaceSidePtr       wire = wires[ iW ];
+    const vector<UVPtStruct>& uvPtVec = wire->GetUVPtStruct();
+    const smIdType         nbSegments = wire->NbPoints() - 1;
+
+    for ( int i = 0; i < nbSegments; ++i ) // loop on segments
+    {
+      const SMDS_MeshNode * n0 = uvPtVec[ i ].node;
+      const SMDS_MeshNode * n1 = uvPtVec[ i + 1 ].node;
+      const int posShapeID = n0->GetShapeID();
+      // skip nodes on degenerated edges
+      if ( helper.IsDegenShape( posShapeID ) && helper.IsDegenShape( n1->GetShapeID() ))
+        continue;
+      
+      double size = SMESH_TNodeXYZ( n0 ).Distance( n1 );
+      
+      netgen::Point3d p0(n0->X(), n0->Y(), n0->Z());
+      netgen::Point3d p1(n1->X(), n1->Y(), n1->Z());
+      // ngMesh.RestrictLocalHLine( p0, p1, size, overrideMinH  );
+      ngMesh.RestrictLocalH( p0, size, overrideMinH  );
+
+      // if ( node2ngID.count(n0) == 0 )
+      // {
+      //   netgen::MeshPoint mp0( p0 );
+      //   ngMesh.AddPoint ( mp0, 1, netgen::FIXEDPOINT );
+      //   node2ngID.insert( make_pair( n0, ngMesh.GetNP() + 1 ) );
+      //   nodeVec.push_back( n0 );
+      // }      
+      // if ( node2ngID.count(n1) == 0 )
+      // {
+      //   netgen::MeshPoint mp1( p1 );
+      //   ngMesh.AddPoint ( mp1, 1, netgen::FIXEDPOINT );
+      //   node2ngID.insert( make_pair( n1, ngMesh.GetNP() + 1 ) );
+      //   nodeVec.push_back( n1 );
+      // }
+    
+    }
+  }
+  
+  return TError();
+}
+
 //================================================================================
 /*!
  * \brief Fill netgen mesh with segments of a FACE
@@ -2291,9 +2426,14 @@ NETGENPlugin_Mesher::AddSegmentsToMesh(netgen::Mesh&                    ngMesh,
   }
 
   const int solidID = 0, faceID = geom.fmap.FindIndex( helper.GetSubShape() );
-  if ( ngMesh.GetNFD() < 1 )
-    ngMesh.AddFaceDescriptor( netgen::FaceDescriptor( faceID, solidID, solidID, 0 ));
-
+  // if ( ngMesh.GetNFD() < 1 )
+  //   ngMesh.AddFaceDescriptor( netgen::FaceDescriptor( faceID, solidID, solidID, 0 ));
+  auto numEdges   = geom.GetNEdges();
+  auto numVertex  = geom.GetNVertices();
+  auto numFaces  = geom.GetNFaces();
+  printf("Geom elements: %d, %d, %d\n", numEdges, numVertex, numFaces );
+
+  // Index meshed edges that were already add to the mesh to avoid add then twice as required by netge!
   for ( size_t iW = 0; iW < wires.size(); ++iW )
   {
     StdMeshers_FaceSidePtr       wire = wires[ iW ];
@@ -2303,6 +2443,7 @@ NETGENPlugin_Mesher::AddSegmentsToMesh(netgen::Mesh&                    ngMesh,
     // assure the 1st node to be in node2ngID, which is needed to correctly
     // "close chain of segments" (see below) in case if the 1st node is not
     // onVertex because it is on a Viscous layer
+    // printf("First reference node is: %0.1lf, %0.1lf, %0.1lf\n", uvPtVec[ 0 ].node->X(), uvPtVec[ 0 ].node->Y(), uvPtVec[ 0 ].node->Z() );
     node2ngID.insert( make_pair( uvPtVec[ 0 ].node, ngMesh.GetNP() + 1 ));
 
     // compute length of every segment
@@ -2310,10 +2451,13 @@ NETGENPlugin_Mesher::AddSegmentsToMesh(netgen::Mesh&                    ngMesh,
     for ( int i = 0; i < nbSegments; ++i )
       segLen[i] = SMESH_TNodeXYZ( uvPtVec[ i ].node ).Distance( uvPtVec[ i+1 ].node );
 
-    int edgeID = 1, posID = -2;
+    int edgeID = -1, posID = -2;
     bool isInternalWire = false;
     double vertexNormPar = 0;
     const int prevNbNGSeg = ngMesh.GetNSeg();
+    int previusEdgeId = -1;
+    const SMDS_MeshNode * lastIndexedNode;
+
     for ( int i = 0; i < nbSegments; ++i ) // loop on segments
     {
       // Add the first point of a segment
@@ -2328,11 +2472,29 @@ NETGENPlugin_Mesher::AddSegmentsToMesh(netgen::Mesh&                    ngMesh,
            helper.IsDegenShape( uvPtVec[ i+1 ].node->GetShapeID() ))
         continue;
 
+      double normParam = 0.5 * ( uvPtVec[ i ].normParam + uvPtVec[ i+1 ].normParam );
+      int edgeIndexInWire = wire->EdgeIndex( normParam );
+      const TopoDS_Edge& edge = wire->Edge( edgeIndexInWire );
+      edgeID = geom.GetEdge(edge).nr;    
+
       int ngID1 = ngMesh.GetNP() + 1, ngID2 = ngID1+1;
 
-      if ( onVertex || ( !wasNgMeshEmpty && onEdge ) || helper.IsRealSeam( posShapeID ))
-        ngID1 = node2ngID.insert( make_pair( n, ngID1 )).first->second;
-      
+      // Treat seam edge differently
+      // if ( helper.IsRealSeam( posShapeID ) && node2ngID.count( n ) != 0 )
+      //   continue;
+
+      if ( onVertex || ( !wasNgMeshEmpty && onEdge ) || helper.IsRealSeam( posShapeID ) )
+      {
+        ngID1 = node2ngID.insert( make_pair( n, ngID1 )).first->second;  
+
+        // printf("OnVertex: %d \n", ngID1 );
+        // if ( node2ngID.count( n ) == 0 )
+        //   node2ngID.insert( make_pair( n, ngID1 ));
+        // ngID1 = node2ngID[ n ];        
+
+        // printf("OnVertex: %d, %0.1lf, %0.1lf, %0.1lf\n", ngID1, n->X(), n->Y(), n->Z() );
+      }
+
       if ( ngID1 > ngMesh.GetNP() )
       {
         netgen::MeshPoint mp( netgen::Point<3> (n->X(), n->Y(), n->Z()) );
@@ -2349,70 +2511,64 @@ NETGENPlugin_Mesher::AddSegmentsToMesh(netgen::Mesh&                    ngMesh,
         }
       }
 
-      // Add the segment
-
       netgen::Segment seg;
 
       seg[0]     = ngID1;                // ng node id
       seg[1]     = ngID2;                // ng node id
-      seg.edgenr = ngMesh.GetNSeg() + 1; // ng segment id
-      seg.si     = faceID;               // = geom.fmap.FindIndex (face);
 
-      for ( int iEnd = 0; iEnd < 2; ++iEnd)
-      {
-        const UVPtStruct& pnt = uvPtVec[ i + iEnd ];
-        seg.epgeominfo[ iEnd ].dist = pnt.param; // param on curve
-        seg.epgeominfo[ iEnd ].u    = pnt.u;
-        seg.epgeominfo[ iEnd ].v    = pnt.v;
-        // find out edge id and node parameter on edge
-        onVertex = ( pnt.normParam + 1e-10 > vertexNormPar );
-        if ( onVertex || posShapeID != posID )
+      // netgen6
+      seg.edgenr = edgeID+1;  // ng segment id
+      seg.si     = edgeID+1;  // = geom.fmap.FindIndex (face);
+
+      const UVPtStruct& pnt0 = uvPtVec[ i ];
+      const UVPtStruct& pnt1 = uvPtVec[ i + 1 ];
+      
+      seg.epgeominfo[ 0 ].dist      = pnt0.param; // param on curve
+      seg.epgeominfo[ 0 ].u         = pnt0.u;
+      seg.epgeominfo[ 0 ].v         = pnt0.v;
+      seg.epgeominfo[ 0 ].edgenr    = edgeID;
+    
+      seg.epgeominfo[ 1 ].dist    = pnt1.param; // param on curve
+      seg.epgeominfo[ 1 ].u       = pnt1.u;
+      seg.epgeominfo[ 1 ].v       = pnt1.v;
+      seg.epgeominfo[ 1 ].edgenr  = edgeID; //
+      
+      // seg.singedge_left   = geom.GetEdge(edge).properties.hpref;
+      // seg.singedge_right  = geom.GetEdge(edge).properties.hpref;
+      // seg.domin   = geom.GetEdge(edge).domin+1;
+      // seg.domout  = geom.GetEdge(edge).domout+1;
+        
+      ngMesh.AddSegment (seg);
+
         {
-          // get edge id
-          double normParam = pnt.normParam;
-          if ( onVertex )
-            normParam = 0.5 * ( uvPtVec[ i ].normParam + uvPtVec[ i+1 ].normParam );
-          int edgeIndexInWire = wire->EdgeIndex( normParam );
-          vertexNormPar = wire->LastParameter( edgeIndexInWire );
-          const TopoDS_Edge& edge = wire->Edge( edgeIndexInWire );
-          edgeID = geom.emap.FindIndex( edge );
-          posID  = posShapeID;
-          isInternalWire = ( edge.Orientation() == TopAbs_INTERNAL );
-          // if ( onVertex ) // param on curve is different on each of two edges
-          //   seg.epgeominfo[ iEnd ].dist = helper.GetNodeU( edge, pnt.node );
+          // restrict size of elements near the segment
+          SMESH_TNodeXYZ np1( n ), np2( uvPtVec[ i+1 ].node );
+          // get an average size of adjacent segments to avoid sharp change of
+          // element size (regression on issue 0020452, note 0010898)
+          int   iPrev = SMESH_MesherHelper::WrapIndex( i-1, (int) nbSegments );
+          int   iNext = SMESH_MesherHelper::WrapIndex( i+1, (int) nbSegments );
+          double sumH = segLen[ iPrev ] + segLen[ i ] + segLen[ iNext ];
+          int   nbSeg = ( int( segLen[ iPrev ] > sumH / 100.)  +
+                          int( segLen[ i     ] > sumH / 100.)  +
+                          int( segLen[ iNext ] > sumH / 100.));
+          if ( nbSeg > 0 )
+            RestrictLocalSize( ngMesh, 0.5*(np1+np2), sumH / nbSeg, overrideMinH );
         }
-        seg.epgeominfo[ iEnd ].edgenr = edgeID; //  = geom.emap.FindIndex(edge);
-      }
 
-      ngMesh.AddSegment (seg);
-      {
-        // restrict size of elements near the segment
-        SMESH_TNodeXYZ np1( n ), np2( uvPtVec[ i+1 ].node );
-        // get an average size of adjacent segments to avoid sharp change of
-        // element size (regression on issue 0020452, note 0010898)
-        int   iPrev = SMESH_MesherHelper::WrapIndex( i-1, (int) nbSegments );
-        int   iNext = SMESH_MesherHelper::WrapIndex( i+1, (int) nbSegments );
-        double sumH = segLen[ iPrev ] + segLen[ i ] + segLen[ iNext ];
-        int   nbSeg = ( int( segLen[ iPrev ] > sumH / 100.)  +
-                        int( segLen[ i     ] > sumH / 100.)  +
-                        int( segLen[ iNext ] > sumH / 100.));
-        if ( nbSeg > 0 )
-          RestrictLocalSize( ngMesh, 0.5*(np1+np2), sumH / nbSeg, overrideMinH );
-      }
-      if ( isInternalWire )
-      {
-        swap (seg[0], seg[1]);
-        swap( seg.epgeominfo[0], seg.epgeominfo[1] );
-        seg.edgenr = ngMesh.GetNSeg() + 1; // segment id
-        ngMesh.AddSegment (seg);
-      }
+      // lastIndexedNode = pnt1.node;
+      
+      // printf("Segments init %d %d %0.1lf, %0.1lf, %0.1lf\n", edgeID, ngID1, pnt0.node->X(), pnt0.node->Y(), pnt0.node->Z() );
+      // printf("Segments end  %d %d %0.1lf, %0.1lf, %0.1lf\n", edgeID, ngID2, pnt1.node->X(), pnt1.node->Y(), pnt1.node->Z() );
+
     } // loop on segments on a wire
 
+    // printf("Out of the segment loop \n");
     // close chain of segments
     if ( nbSegments > 0 )
     {
       netgen::Segment& lastSeg = ngMesh.LineSegment( ngMesh.GetNSeg() - int( isInternalWire ));
       const SMDS_MeshNode * lastNode = uvPtVec.back().node;
+
       lastSeg[1] = node2ngID.insert( make_pair( lastNode, lastSeg[1] )).first->second;
       if ( lastSeg[1] > ngMesh.GetNP() )
       {
@@ -2458,16 +2614,34 @@ NETGENPlugin_Mesher::AddSegmentsToMesh(netgen::Mesh&                    ngMesh,
   } // loop on WIREs of a FACE
 
   // add a segment instead of an internal vertex
-  if ( wasNgMeshEmpty )
-  {
-    NETGENPlugin_Internals intShapes( *helper.GetMesh(), helper.GetSubShape(), /*is3D=*/false );
-    AddIntVerticesInFaces( geom, ngMesh, nodeVec, intShapes );
-  }
+    // if ( wasNgMeshEmpty )
+    // {
+    //   NETGENPlugin_Internals intShapes( *helper.GetMesh(), helper.GetSubShape(), /*is3D=*/false );
+    //   AddIntVerticesInFaces( geom, ngMesh, nodeVec, intShapes );
+    // }
+
+
+  // Get the add segment and check 
+
+  printf("Going to call surface\n");
   ngMesh.CalcSurfacesOfNode();
+  printf("End to call surface, %d,%d\n", ngMesh.GetNP(),  ngMesh.GetNSeg() );
 
   return TError();
 }
 
+
+// int NETGENPlugin_Mesher::CheckSegmentChain(const netgen::OCCGeometry&          occgeo,
+//                                             netgen::Mesh&                       ngMesh )
+// {
+
+//   auto segments = occgeo.GetFace(0).GetBoundary();
+//   for (auto s : segments)
+//   {
+
+//   }
+// }
+
 //================================================================================
 /*!
  * \brief Fill SMESH mesh according to contents of netgen mesh
@@ -2494,6 +2668,7 @@ int NETGENPlugin_Mesher::FillSMesh(const netgen::OCCGeometry&          occgeo,
   int nbSeg = ngMesh.GetNSeg();
   int nbFac = ngMesh.GetNSE();
   int nbVol = ngMesh.GetNE();
+  printf( "Numbers: %d,%d,%d,%d\n", nbNod, nbSeg, nbFac, nbVol );
   SMESHDS_Mesh* meshDS = sMesh.GetMeshDS();
 
   // quadHelper is used for either
@@ -3262,6 +3437,7 @@ void NETGENPlugin_Mesher::SetBasicMeshParameters( NETGENPlugin_NetgenLibWrapper&
 
   _ngMesh->ClearFaceDescriptors(); // we make descriptors our-self
 
+  /* This line breaks when using _simpleHyp for 1D2D version*/
   if ( !mparams.uselocalh ) // mparams.grading is not taken into account yet
     _ngMesh->LocalHFunction().SetGrading( mparams.grading );
 
@@ -3273,7 +3449,7 @@ void NETGENPlugin_Mesher::SetBasicMeshParameters( NETGENPlugin_NetgenLibWrapper&
     double segSize = _simpleHyp->GetLocalLength();
     for ( int iE = 1; iE <= occgeo.emap.Extent(); ++iE )
     {
-      const TopoDS_Edge& e = TopoDS::Edge( occgeo.emap(iE));
+      const TopoDS_Edge& e = TopoDS::Edge(occgeo.emap(iE));
       if ( nbSeg )
         segSize = SMESH_Algo::EdgeLength( e ) / ( nbSeg - 0.4 );
       setLocalSize( e, segSize, *_ngMesh );
@@ -3420,7 +3596,8 @@ void NETGENPlugin_Mesher::InitialSetupSA( NETGENPlugin_NetgenLibWrapper& ngLib,
     updateTriangulation( _shape );
 
     Bnd_Box bb;
-    BRepBndLib::Add (_shape, bb);
+    // BRepBndLib::Add (_shape, bb);
+    GEOMUtils::PreciseBoundingBox( _shape, bb );
     double x1,y1,z1,x2,y2,z2;
     bb.Get (x1,y1,z1,x2,y2,z2);
     netgen::Point<3> p1 = netgen::Point<3> (x1,y1,z1);
@@ -3665,6 +3842,7 @@ bool NETGENPlugin_Mesher::Compute()
   //int nbSeg = _ngMesh->GetNSeg();
   int nbFac = _ngMesh->GetNSE();
   int nbVol = _ngMesh->GetNE();
+
   bool isOK = ( !err && (_isVolume ? (nbVol > 0) : (nbFac > 0)) );
 
   // Feed back the SMESHDS with the generated Nodes and Elements
index 173ca086eca2cf89785b4d1d7870f02c5c30c319..b56a9880631622303a44e284919aab140d716e34 100644 (file)
@@ -293,6 +293,13 @@ class NETGENPLUGIN_EXPORT NETGENPlugin_Mesher
                       std::vector< const SMDS_MeshNode* > & nodeVec,
                       const bool                            overrideMinH=true);
 
+  static SMESH_ComputeErrorPtr AddSegmentsToMesh2(netgen::Mesh&                         ngMesh,
+                                                    netgen::OCCGeometry&                  geom,
+                                                    const TSideVector&                    wires,
+                                                    SMESH_MesherHelper&                   helper,
+                                                    std::vector< const SMDS_MeshNode* > & nodeVec,
+                                                    const bool                            overrideMinH=true);
+
   void SetDefaultParameters();
 
   static SMESH_ComputeErrorPtr ReadErrors(const std::vector< const SMDS_MeshNode* >& nodeVec);
index 6d3720590e6fbb03e7573e4862d649f7ae7596d6..1ea1a83b6c58f18fbfdf8279e62a839ff82d7703 100644 (file)
@@ -549,6 +549,9 @@ bool NETGENPlugin_NETGEN_2D_ONLY::MapSegmentsToEdges(SMESH_Mesh& aMesh, const To
     occgeom.face_maxh_modified = 0;
     occgeom.face_maxh.SetSize(1);
     occgeom.face_maxh = netgen::mparam.maxh;
+    occgeom.BuildFMapFromPrefilled();
+    
+    printf("Calling NETGENPlugin_NETGEN_2D_ONLY::MapSegmentsToEdges\n");
 
     // Set the face descriptor
     const int solidID = 0, faceID = 1; /*always 1 because faces are meshed one by one*/
@@ -922,7 +925,7 @@ void NETGENPlugin_NETGEN_2D_ONLY::FillNodesAndElements( SMESH_Mesh& aMesh, SMESH
   for ( int ngID = nbInputNodes + 1; ngID <= nbNodes; ++ngID )
   {
     const MeshPoint& ngPoint = ngMesh->Point( ngID );
-    SMDS_MeshNode * node = meshDS->AddNode(ngPoint(0), ngPoint(1), ngPoint(2));
+    SMDS_MeshNode * node = meshDS->AddNode(ngPoint(0), ngPoint(1), ngPoint(2));  
     nodeVec[ ngID ] = node;
   }
 
@@ -1057,17 +1060,13 @@ bool NETGENPlugin_NETGEN_2D_ONLY::Compute(SMESH_Mesh&         aMesh,
       return setMaxh;
       
     // prepare occgeom
-    netgen::OCCGeometry occgeom;
-    occgeom.shape = F;
-    occgeom.fmap.Add( F );
-    occgeom.CalcBoundingBox();
-    occgeom.facemeshstatus.SetSize(1);
-    occgeom.facemeshstatus = 0;
-    occgeom.face_maxh_modified.SetSize(1);
-    occgeom.face_maxh_modified = 0;
-    occgeom.face_maxh.SetSize(1);
-    occgeom.face_maxh = netgen::mparam.maxh;
+    printf( "going to define occgeom\n");
 
+    netgen::OCCGeometry occgeom;
+    
+    occgeom.shape = F;   
+    occgeom.BuildFMapFace();
+    
     // -------------------------
     // Fill netgen mesh
     // -------------------------
@@ -1102,8 +1101,12 @@ bool NETGENPlugin_NETGEN_2D_ONLY::Compute(SMESH_Mesh&         aMesh,
       }
 
       nodeVec.clear();
-      faceErr = aMesher.AddSegmentsToMesh( *ngMesh, occgeom, wires, helper, nodeVec,
+      // faceErr = aMesher.AddSegmentsToMesh( *ngMesh, occgeom, wires, helper, nodeVec,
+      //                                      /*overrideMinH=*/!_hypParameters);
+
+      faceErr = aMesher.AddSegmentsToMesh2( *ngMesh, occgeom, wires, helper, nodeVec,
                                            /*overrideMinH=*/!_hypParameters);
+
       if ( faceErr && !faceErr->IsOK() )
         break;
 
@@ -1113,15 +1116,19 @@ bool NETGENPlugin_NETGEN_2D_ONLY::Compute(SMESH_Mesh&         aMesh,
       // -------------------------
       // Generate surface mesh
       // -------------------------
-
-      const int startWith = MESHCONST_MESHSURFACE;
+      const int startWith = netgen::MESHCONST_MESHEDGES;
       const int endWith   = toOptimize ? MESHCONST_OPTSURFACE : MESHCONST_MESHSURFACE;
 
+      // const int startWith = MESHCONST_MESHSURFACE;
+      // const int endWith   = toOptimize ? MESHCONST_OPTSURFACE : MESHCONST_MESHSURFACE;
+      // const int endWith   = MESHCONST_MESHSURFACE;
+
       SMESH_Comment str;
       try {
         OCC_CATCH_SIGNALS;
-
+        printf("Going to generate mesh\n");
         err = ngLib.GenerateMesh(occgeom, startWith, endWith, ngMesh);
+        printf("Mesh Generated: %d,%d,%d \n", ngMesh->GetNP(), ngMesh->GetNSeg(), ngMesh->GetNSE() );
 
         if ( netgen::multithread.terminate )
           return false;
index 2162e0460c238b4a15dc3f8a4d6b3d4603f89577..c046ed5dcc7c246368f2ebf38ebf4df8fcea4f28 100644 (file)
@@ -231,7 +231,7 @@ namespace
 
   void HoleFiller::AddHoleBordersAndEdges( Ng_STL_Geometry * ngStlGeo, bool toAddEdges )
   {
-    nglib::readedges.SetSize(0);
+    // nglib::readedges.SetSize(0);
 
     for ( size_t i = 0; i < myHole.size(); ++i )
       for ( size_t iP = 1; iP < myHole[i].size(); ++iP )
@@ -516,9 +516,9 @@ namespace
         n->setIsMarked( true );
 
         SMESH_NodeXYZ p( n );
-        int id = stlGeom->GetPointNum( netgen::Point<3>( p.X(),p.Y(),p.Z() ));
-        if ( id > 0 )
-          stlGeom->SetLineEndPoint( id );
+        // int id = stlGeom->GetPointNum( netgen::Point<3>( p.X(),p.Y(),p.Z() ));
+        // if ( id > 0 )
+        //   stlGeom->SetLineEndPoint( id );
       }
     }
   }