bos #18711 [CEA 18704] COTECH action 117.4: NETGEN 6 Integration in SALOME

cleanup + PATCH

diff --git a/src/NETGEN/netgen62ForSalome.patch b/src/NETGEN/netgen62ForSalome.patch
new file mode 100644 (file)
index 0000000..20c187b
--- /dev/null
+++ b/src/NETGEN/netgen62ForSalome.patch
@@ -0,0 +1,293 @@
+diff --git a/libsrc/occ/occgenmesh.cpp b/libsrc/occ/occgenmesh.cpp
+index 314d405a..7c912bcc 100644
+--- a/libsrc/occ/occgenmesh.cpp
++++ b/libsrc/occ/occgenmesh.cpp
+@@ -49,7 +49,7 @@ namespace netgen
+   double ComputeH (double kappa, const MeshingParameters & mparam)\r
+   {\r
+     kappa *= mparam.curvaturesafety;\r
+-    /*\r
++    /**/\r
+     double hret;\r
+ \r
+     if (mparam.maxh * kappa < 1)\r
+@@ -61,7 +61,7 @@ namespace netgen
+       hret = mparam.maxh;\r
+ \r
+     return hret;\r
+-    */\r
++    /**/\r
+     // return min(mparam.maxh, 1/kappa);\r
+     return (mparam.maxh*kappa < 1) ? mparam.maxh : 1/kappa;\r
+   }\r
+@@ -253,10 +253,12 @@ namespace netgen
+       {\r
+         oldpnt = pnt;\r
+         pnt = c->Value(s0+(i/double(DIVIDEEDGESECTIONS))*(s1-s0));\r
++        // !!! no more than 1 segment per <edge length>/DIVIDEEDGESECTIONS\r
+         hvalue[i] = hvalue[i-1] +\r
+-          1.0/mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z()))*\r
+-          pnt.Distance(oldpnt);\r
++          min( 1.0,\r
++               1.0/mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z()))*\r
++               pnt.Distance(oldpnt));\r
+ \r
+         //(*testout) << "mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z())) " << mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z()))\r
+         //       <<  " pnt.Distance(oldpnt) " << pnt.Distance(oldpnt) << endl;\r
+ \r
+@@ -299,7 +301,7 @@ namespace netgen
+         cout << "CORRECTED" << endl;\r
+         ps.SetSize (nsubedges-2);\r
+         params.SetSize (nsubedges);\r
+-        params[nsubedges] = s1;\r
++        params[nsubedges-1] = s1;\r
+       }\r
+   }\r
+ \r
+@@ -323,6 +325,8 @@ namespace netgen
+ \r
+     double eps = 1e-6 * geom.GetBoundingBox().Diam();\r
+ \r
++    int first_vp = mesh.GetNP()+1;\r
++\r
+     tsearch.Start();\r
+     for (int i = 1; i <= nvertices; i++)\r
+       {\r
+@@ -481,22 +485,107 @@ namespace netgen
+ \r
+                 if (!merge_solids)\r
+                   {\r
+-                    pnums[0] = geom.vmap.FindIndex (TopExp::FirstVertex (edge)) + PointIndex::BASE-1;\r
+-                    pnums.Last() = geom.vmap.FindIndex (TopExp::LastVertex (edge)) + PointIndex::BASE-1;\r
++                    //pnums[0] = geom.vmap.FindIndex (TopExp::FirstVertex (edge)) + PointIndex::BASE-1;\r
++                    //pnums.Last() = geom.vmap.FindIndex (TopExp::LastVertex (edge)) + PointIndex::BASE-1;\r
++                    const int dpi = PointIndex::BASE-1;\r
++                    MeshPoint dfltP ( Point<3> ( 0, 0, 0 ) );\r
++                    PointIndex *ipp[] = { &pnums[0], &pnums.Last() };\r
++                    TopoDS_Iterator vIt( edge, false );\r
++                    TopoDS_Vertex v[2];\r
++                    v[0] = TopoDS::Vertex( vIt.Value() ); vIt.Next();\r
++                    v[1] = TopoDS::Vertex( vIt.Value() );\r
++                    if ( v[0].Orientation() == TopAbs_REVERSED )\r
++                      std::swap( v[0], v[1] );\r
++                    for ( int i = 0; i < 2; ++i)\r
++                    {\r
++                      PointIndex &ip = *ipp[i];\r
++                      ip = geom.vmap.FindIndex ( v[i] ) + dpi;\r
++                      if ( ip == dpi || ip > nvertices + dpi )\r
++                      {\r
++                        PointIndex iv = ip;\r
++                        if ( ip == dpi )\r
++                          ip = iv = const_cast<OCCGeometry&>(geom).vmap.Add( v[i] );\r
++                        gp_Pnt pnt = BRep_Tool::Pnt( v[i] );\r
++                        MeshPoint mp( Point<3>(pnt.X(), pnt.Y(), pnt.Z()) );\r
++                        for (PointIndex pi = 1; pi < first_vp; pi++)\r
++                          if ( Dist2 (mesh.Point(pi), Point<3>(mp)) < 1e-100 )\r
++                          {\r
++                            ip = pi;\r
++                            if ( mesh.Point(ip).GetLayer() != dfltP.GetLayer() && mesh.Point(ip).GetLayer() != iv )\r
++                              continue;\r
++                            if ( mesh.Point(ip).GetLayer() == dfltP.GetLayer())\r
++                              mesh.Point(ip) = MeshPoint( mesh.Point(ip), iv );\r
++                            break;\r
++                          }\r
++                      }\r
++                      else\r
++                      {\r
++                        ip += first_vp - 1;\r
++                      }\r
++                    }\r
+                   }\r
+                 else\r
+                   {\r
+-                    Point<3> fp = occ2ng (BRep_Tool::Pnt (TopExp::FirstVertex (edge)));\r
+-                    Point<3> lp = occ2ng (BRep_Tool::Pnt (TopExp::LastVertex (edge)));\r
++                    TopoDS_Iterator vIt( edge, false );\r
++                    TopoDS_Vertex v1 = TopoDS::Vertex( vIt.Value() ); vIt.Next();\r
++                    TopoDS_Vertex v2 = TopoDS::Vertex( vIt.Value() );\r
++                    if ( v1.Orientation() == TopAbs_REVERSED )\r
++                      std::swap( v1, v2 );\r
++                    const bool isClosedEdge = v1.IsSame( v2 );\r
++\r
++                    Point<3> fp = occ2ng (BRep_Tool::Pnt (v1));\r
++                    Point<3> lp = occ2ng (BRep_Tool::Pnt (v2));\r
++                    double tol2 = std::min( eps*eps, 1e-6 * Dist2( fp, lp ));\r
++                    if ( isClosedEdge )\r
++                      tol2 = BRep_Tool::Tolerance( v1 ) * BRep_Tool::Tolerance( v1 );\r
+ \r
+                     pnums[0] = PointIndex::INVALID;\r
+                     pnums.Last() = PointIndex::INVALID;\r
+                     for (PointIndex pi : vertexrange)\r
+                       {\r
+-                        if (Dist2 (mesh[pi], fp) < eps*eps) pnums[0] = pi;\r
+-                        if (Dist2 (mesh[pi], lp) < eps*eps) pnums.Last() = pi;\r
++                        if (Dist2 (mesh[pi], fp) < tol2) pnums[0] = pi;\r
++                        if (Dist2 (mesh[pi], lp) < tol2) pnums.Last() = pi;\r
+                       }\r
+-                  }\r
++                    if (( isClosedEdge && pnums[0] != pnums.Last() ) ||\r
++                        ( !isClosedEdge && pnums[0] == pnums.Last() ))\r
++                      pnums[0] = pnums.Last() = PointIndex::INVALID;\r
++                    if ( pnums[0] < PointIndex::BASE || pnums.Last() < PointIndex::BASE )\r
++                    {\r
++                      // take into account a possible large gap between a vertex and an edge curve\r
++                      // end and a large vertex tolerance covering the whole edge\r
++                      if ( pnums[0] < PointIndex::BASE )\r
++                      {\r
++                        double tol = BRep_Tool::Tolerance( v1 );\r
++                        for (PointIndex pi = 1; pi < first_ep; pi++)\r
++                          if (pi != pnums.Last() && Dist2 (mesh[pi], fp) < 2*tol*tol)\r
++                            pnums[0] = pi;\r
++\r
++                        if ( pnums[0] < PointIndex::BASE )\r
++                          pnums[0] = first_ep-1 - nvertices + geom.vmap.FindIndex ( v1 );\r
++                      }\r
++                      if ( isClosedEdge )\r
++                      {\r
++                        pnums.Last() = pnums[0];\r
++                      }\r
++                      else\r
++                      {\r
++                        if ( pnums.Last() < PointIndex::BASE )\r
++                        {\r
++                          double tol = BRep_Tool::Tolerance( v2 );\r
++                          for (PointIndex pi = 1; pi < first_ep; pi++)\r
++                            if (pi != pnums[0] && Dist2 (mesh[pi], lp) < 2*tol*tol)\r
++                              pnums.Last() = pi;\r
++\r
++                          if ( pnums.Last() < PointIndex::BASE )\r
++                            pnums.Last() = first_ep-1-nvertices + geom.vmap.FindIndex ( v2 );\r
++                        }\r
++\r
++                        if ( Dist2( fp, mesh[PointIndex(pnums[0])]) >\r
++                             Dist2( lp, mesh[PointIndex(pnums.Last())]))\r
++                          std::swap( pnums[0], pnums.Last() );\r
++                      }\r
++                    }\r
++                 }\r
+ \r
+                 for (size_t i = 1; i <= mp.Size(); i++)\r
+                   {\r
+@@ -505,17 +594,19 @@ namespace netgen
+ \r
+                     // for (PointIndex j = first_ep; j < mesh.Points().End(); j++)\r
+                     for (PointIndex j = first_ep; j < *mesh.Points().Range().end(); j++)\r
++                    {\r
++                      if (!merge_solids && mesh.Point(j).GetLayer() != geomedgenr ) continue;\r
+                       if ((mesh.Point(j)-Point<3>(mp[i-1])).Length() < eps)\r
+                         {\r
+                           exists = true;\r
+                           pnums[i] = j;\r
+                           break;\r
+                         }\r
+-                      \r
++                    }\r
+                     tsearch.Stop();\r
+                     \r
+                     if (!exists)\r
+-                      pnums[i] = mesh.AddPoint (mp[i-1]);\r
++                      pnums[i] = mesh.AddPoint (mp[i-1], geomedgenr);\r
+                   }\r
+                 if(geom.enames.Size() && geom.enames[curr-1] != "")\r
+                   mesh.SetCD2Name(geomedgenr, geom.enames[curr-1]);\r
+@@ -599,6 +690,9 @@ namespace netgen
+     //                                << " p1 " << mesh.LineSegment(i)[0] << " p2 " << mesh.LineSegment(i)[1] << endl;\r
+     //        exit(10);\r
+ \r
++    for (int j = 1; j <= mesh.GetNP(); j++) // support SALOME fuse edges: set level to zero\r
++      mesh.Point(j) = MeshPoint( (Point<3>&) mesh.Point(j) );\r
++\r
+     mesh.CalcSurfacesOfNode();\r
+     multithread.task = savetask;\r
+   }\r
+@@ -652,7 +746,7 @@ namespace netgen
+ \r
+         Box<3> bb = geom.GetBoundingBox();\r
+ \r
+-        int projecttype = PLANESPACE;\r
++        //int projecttype = PLANESPACE;\r
+         \r
+         static Timer tinit("init");\r
+         tinit.Start();\r
+diff --git a/libsrc/occ/occmeshsurf.cpp b/libsrc/occ/occmeshsurf.cpp
+index 7fa5d237..9e05fd95 100644
+--- a/libsrc/occ/occmeshsurf.cpp
++++ b/libsrc/occ/occmeshsurf.cpp
+@@ -173,50 +173,6 @@ namespace netgen
+       gp_Vec du, dv;
+       occface->D1 (geominfo1.u, geominfo1.v, pnt, du, dv);
+ 
+-        // static Timer t("occ-defintangplane calculations");
+-        // RegionTimer reg(t);
+-
+-        Mat<3,2> D1_;
+-      D1_(0,0) = du.X(); D1_(1,0) = du.Y(); D1_(2,0) = du.Z();
+-      D1_(0,1) = dv.X(); D1_(1,1) = dv.Y(); D1_(2,1) = dv.Z();
+-        auto D1T_ = Trans(D1_);
+-      auto D1TD1_ = D1T_*D1_;
+-      if (Det (D1TD1_) == 0) throw SingularMatrixException();
+-        Mat<2,2> DDTinv_;
+-        CalcInverse (D1TD1_, DDTinv_);
+-
+-        Mat<3,2> Y_;
+-      Vec<3> y1_ = (ap2-ap1).Normalize();
+-      Vec<3> y2_ = Cross(n, y1_).Normalize();
+-      for (int i = 0; i < 3; i++)
+-        {
+-          Y_(i,0) = y1_(i);
+-          Y_(i,1) = y2_(i);
+-        }
+-
+-        auto A_ = DDTinv_ * D1T_ * Y_;
+-      Mat<2,2> Ainv_;
+-      if (Det(A_) == 0) throw SingularMatrixException();
+-      CalcInverse (A_, Ainv_);
+-
+-      Vec<2> temp_ = Ainv_ * (psp2-psp1);
+-      double r_ = temp_.Length();
+-        Mat<2,2> R_;
+-        R_(0,0) = temp_(0)/r_;
+-        R_(1,0) = temp_(1)/r_;
+-        R_(0,1) = -R_(1,0);
+-        R_(1,1) = R_(0,0);
+-
+-        A_ = A_ * R_;
+-        Ainv_ = Trans(R_) * Ainv_;
+-
+-        Amat = A_;
+-        Amatinv = Ainv_;
+-        
+-      // temp = Amatinv * (psp2-psp1);
+-        
+-
+-#ifdef OLD
+       DenseMatrix D1(3,2), D1T(2,3), DDTinv(2,2);
+       D1(0,0) = du.X(); D1(1,0) = du.Y(); D1(2,0) = du.Z();
+       D1(0,1) = dv.X(); D1(1,1) = dv.Y(); D1(2,1) = dv.Z();
+@@ -289,8 +245,7 @@ namespace netgen
+           }
+         // cout << "=?= Ainv = " << endl << Ainv << endl;
+       temp = Amatinv * (psp2-psp1);
+-        cout << " =?= Amatinv = " << Amatinv << endl;        
+-#endif
++        // cout << " =?= Amatinv = " << Amatinv << endl;        
+       };
+  
+   }
+@@ -380,6 +335,7 @@ namespace netgen
+ 
+     double u = gi.u;
+     double v = gi.v;
++    if ( 0 ) { // exclude Newton's method
+     gp_Pnt x = occface->Value (u,v);
+ 
+     if (p.SquareDistance(x) <= sqr(PROJECTION_TOLERANCE)) return;
+@@ -429,7 +385,7 @@ namespace netgen
+           }
+       }
+     while (count < 20);
+-    
++    }    
+ 
+     // Newton did not converge, use OCC projection
+ 

diff --git a/src/NETGENPlugin/NETGENPlugin_Mesher.cxx b/src/NETGENPlugin/NETGENPlugin_Mesher.cxx
index ed6021089369bdb10565e735407cb2bb8d4492e6..5853a7e4a870a2cb051e71ad3c4cfd4196aada55 100644 (file)
--- a/src/NETGENPlugin/NETGENPlugin_Mesher.cxx
+++ b/src/NETGENPlugin/NETGENPlugin_Mesher.cxx
@@ -669,6 +669,12 @@ void NETGENPlugin_Mesher::SetParameters(const NETGENPlugin_Hypothesis* hyp)
       }
     }
   }
+
+#ifdef NETGEN_V6
+
+  netgen::mparam.closeedgefac = 2;
+
+#endif
 }
 
 //=============================================================================
@@ -732,7 +738,7 @@ void NETGENPlugin_Mesher::SetLocalSize( netgen::OCCGeometry& occgeo,
     if ( faceNgID >= 1 )
     {
 #ifdef NETGEN_V6
-      occgeo.SetFaceMaxH(faceNgID-1, val);
+      occgeo.SetFaceMaxH(faceNgID-1, val, netgen::mparam);
 #else
       occgeo.SetFaceMaxH(faceNgID, val);
 #endif
@@ -801,7 +807,7 @@ void NETGENPlugin_Mesher::SetLocalSizeForChordalError( netgen::OCCGeometry& occg
       double maxCurv = Max( Abs( surfProp.MaxCurvature()), Abs( surfProp.MinCurvature() ));
       double    size = elemSizeForChordalError( _chordalError, 1 / maxCurv );
 #ifdef NETGEN_V6
-      occgeo.SetFaceMaxH( i-1, size * sizeCoef );
+      occgeo.SetFaceMaxH( i-1, size * sizeCoef, netgen::mparam );
 #else
       occgeo.SetFaceMaxH( i, size * sizeCoef );
 #endif
@@ -3814,7 +3820,7 @@ NETGENPlugin_Mesher::ReadErrors(const vector<const SMDS_MeshNode* >& nodeVec)
       ok = ok && file.getInts( three2 );
       for ( int i = 0; ok && i < 3; ++i )
         ok = ( three1[i] < nbNodes && nodeVec[ three1[i]]);
-      for ( int i = 0; ok && i < 3; ++i ) 
+      for ( int i = 0; ok && i < 3; ++i )
         ok = ( three2[i] < nbNodes && nodeVec[ three2[i]]);
       if ( ok )
       {
@@ -3865,6 +3871,59 @@ void NETGENPlugin_Mesher::toPython( const netgen::Mesh* ngMesh )
           << "mesh = smesh.Mesh()" << std::endl << std::endl;
 
   using namespace netgen;
+
+#ifdef NETGEN_V6
+
+  for ( int i = 1; i <= ngMesh->GetNP(); i++)
+  {
+    const Point3d & p = ngMesh->Point(i);
+    outfile << "mesh.AddNode( ";
+    outfile << p.X() << ", ";
+    outfile << p.Y() << ", ";
+    outfile << p.Z() << ") ## "<< i << std::endl;
+  }
+
+  int nbDom = ngMesh->GetNDomains();
+  for ( int i = 0; i < nbDom; ++i )
+    outfile<< "grp" << i+1 << " = mesh.CreateEmptyGroup( SMESH.FACE, 'domain"<< i+1 << "')"<< std::endl;
+
+  for (int i = 1; i <= ngMesh->GetNSE(); i++)
+  {
+    outfile << "mesh.AddFace([ ";
+    Element2d sel = ngMesh->SurfaceElement(i);
+    for (int j = 1; j <= sel.GetNP(); j++)
+      outfile << sel.PNum(j) << ( j < sel.GetNP() ? ", " : " ])");
+    if ( sel.IsDeleted() ) outfile << " ## IsDeleted ";
+    outfile << std::endl;
+
+    if (sel.GetIndex())
+    {
+      if ( int dom1 = ngMesh->GetFaceDescriptor(sel.GetIndex ()).DomainIn())
+        outfile << "grp"<< dom1 <<".Add([ " << i << " ])" << std::endl;
+      if ( int dom2 = ngMesh->GetFaceDescriptor(sel.GetIndex ()).DomainOut())
+        outfile << "grp"<< dom2 <<".Add([ " << i << " ])" << std::endl;
+    }
+  }
+
+  for (int i = 1; i <= ngMesh->GetNE(); i++)
+  {
+    Element el = ngMesh->VolumeElement(i);
+    outfile << "mesh.AddVolume([ ";
+    for (int j = 1; j <= el.GetNP(); j++)
+      outfile << el.PNum(j) << ( j < el.GetNP() ? ", " : " ])");
+    outfile << std::endl;
+  }
+
+  for (int i = 1; i <= ngMesh->GetNSeg(); i++)
+  {
+    const Segment & seg = ngMesh->LineSegment (i);
+    outfile << "mesh.AddEdge([ "
+            << seg[0]+1 << ", "
+            << seg[1]+1 << " ])" << std::endl;
+  }
+
+#else  //////// V 5
+
   PointIndex pi;
   for (pi = PointIndex::BASE; 
        pi < ngMesh->GetNP()+PointIndex::BASE; pi++)
@@ -3914,6 +3973,9 @@ void NETGENPlugin_Mesher::toPython( const netgen::Mesh* ngMesh )
             << seg[0] << ", "
             << seg[1] << " ])" << std::endl;
   }
+
+#endif
+
   std::cout << "Write " << pyFile << std::endl;
 }
 

diff --git a/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.cxx b/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.cxx
index 84d4363fbf46589663fb73ddbfe79fbafc860752..cef4819b00f05b122ff916d5e952e0ad03d03a48 100644 (file)
--- a/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.cxx
+++ b/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.cxx
@@ -65,7 +65,9 @@ namespace nglib {
 namespace netgen {
   NETGENPLUGIN_DLL_HEADER
   extern MeshingParameters mparam;
+#ifdef NETGEN_V5
   extern void OCCSetLocalMeshSize(OCCGeometry & geom, Mesh & mesh);
+#endif
 }
 
 using namespace std;
@@ -280,15 +282,19 @@ bool NETGENPlugin_NETGEN_2D_ONLY::Compute(SMESH_Mesh&         aMesh,
     }
     // set local size depending on curvature and NOT closeness of EDGEs
 #ifdef NETGEN_V6
-    const double factor = *netgen::mparam.closeedgefac;
-    netgen::mparam.closeedgefac = std::nullopt;
+    const double factor = 2; //netgen::occparam.resthcloseedgefac;
 #else
     const double factor = netgen::occparam.resthcloseedgefac;
     netgen::occparam.resthcloseedgeenable = false;
-    //netgen::occparam.resthcloseedgefac = 1.0 + netgen::mparam.grading;
+    netgen::occparam.resthcloseedgefac = 1.0 + netgen::mparam.grading;
 #endif
     occgeoComm.face_maxh = netgen::mparam.maxh;
+#ifdef NETGEN_V6
+    netgen::OCCParameters occparam;
+    netgen::OCCSetLocalMeshSize( occgeoComm, *ngMeshes[0], netgen::mparam, occparam );
+#else
     netgen::OCCSetLocalMeshSize( occgeoComm, *ngMeshes[0] );
+#endif
     occgeoComm.emap.Clear();
     occgeoComm.vmap.Clear();
 

diff --git a/src/NETGENPlugin/NETGENPlugin_NETGEN_3D.cxx b/src/NETGENPlugin/NETGENPlugin_NETGEN_3D.cxx
index de45d54d0bd11cc0ece20cbe0d069de69c7812dd..88abe103ead5411e68e3a1ace13582dd166a94dc 100644 (file)
--- a/src/NETGENPlugin/NETGENPlugin_NETGEN_3D.cxx
+++ b/src/NETGENPlugin/NETGENPlugin_NETGEN_3D.cxx
@@ -76,18 +76,13 @@
 #include <ngexception.hpp>
 #endif
 #ifdef NETGEN_V6
-#include <exception.hpp>
+#include <core/exception.hpp>
 #endif
 
 namespace nglib {
 #include <nglib.h>
 }
 namespace netgen {
-#ifdef NETGEN_V5
-  extern int OCCGenerateMesh (OCCGeometry&, Mesh*&, MeshingParameters&, int, int);
-#else
-  extern int OCCGenerateMesh (OCCGeometry&, Mesh*&, int, int, char*);
-#endif
 
   NETGENPLUGIN_DLL_HEADER
   extern MeshingParameters mparam;