add closeedge meshsize to base geometry (not used)

closedgefac moved to meshingparameters for this

libsrc/meshing/basegeom.cpp

@@ -85,6 +85,23 @@ namespace netgen
       }
   }
 
+  struct Line
+  {
+    Point<3> p0, p1;
+    inline double Length() const { return (p1-p0).Length(); }
+    inline double Dist(const Line& other) const
+    {
+      Vec<3> n = p1-p0;
+      Vec<3> q = other.p1-other.p0;
+      double nq = n*q;
+      Point<3> p = p0 + 0.5*n;
+      double lambda = (p-other.p0)*n / (nq + 1e-10);
+      if (lambda >= 0 && lambda <= 1)
+        return (p-other.p0-lambda*q).Length();
+      return 1e99;
+    }
+  };
+
   void NetgenGeometry :: Analyse(Mesh& mesh,
                                  const MeshingParameters& mparam) const
   {
@@ -100,7 +117,7 @@ namespace netgen
 
     if(mparam.uselocalh)
       {
-        double eps = 1e-12 * bounding_box.Diam();
+        double eps = 1e-10 * bounding_box.Diam();
         const char* savetask = multithread.task;
         multithread.task = "Analyse Edges";
 
@@ -142,8 +159,84 @@ namespace netgen
             const auto& face = faces[i];
             face->RestrictH(mesh, mparam);
           }
+
+        if(mparam.closeedgefac.has_value())
+          {
+            multithread.task = "Analyse close edges";
+            constexpr int sections = 100;
+            Array<Line> lines;
+            lines.SetAllocSize(sections*edges.Size());
+            BoxTree<3> searchtree(bounding_box.PMin(),
+                                  bounding_box.PMax());
+            for(const auto& edge : edges)
+              {
+                if(edge->GetLength() < eps)
+                  continue;
+                double t = 0.;
+                auto p_old = edge->GetPoint(t);
+                auto t_old = edge->GetTangent(t);
+                t_old.Normalize();
+                for(auto i : IntRange(1, sections+1))
+                  {
+                    t = double(i)/sections;
+                    auto p_new = edge->GetPoint(t);
+                    auto t_new = edge->GetTangent(t);
+                    t_new.Normalize();
+                    auto cosalpha = fabs(t_old * t_new);
+                    if((i == sections) || (cosalpha < cos(10./180 * M_PI)))
+                      {
+                        auto index = lines.Append({p_old, p_new});
+                        searchtree.Insert(p_old, p_new, index);
+                        p_old = p_new;
+                        t_old = t_new;
+                      }
+                  }
+              }
+            Array<int> linenums;
+            for(auto i : Range(lines))
+              {
+                const auto& line = lines[i];
+                if(line.Length() < eps) continue;
+                multithread.percent = 100.*i/lines.Size();
+                Box<3> box;
+                box.Set(line.p0);
+                box.Add(line.p1);
+                // box.Increase(max2(mesh.GetH(line.p0), mesh.GetH(line.p1)));
+                box.Increase(line.Length());
+                double mindist = 1e99;
+                linenums.SetSize0();
+                searchtree.GetIntersecting(box.PMin(), box.PMax(),
+                                           linenums);
+                for(auto num : linenums)
+                  {
+                    if(i == num) continue;
+                    const auto & other = lines[num];
+                    if((line.p0 - other.p0).Length2() < eps ||
+                       (line.p0 - other.p1).Length2() < eps ||
+                       (line.p1 - other.p0).Length2() < eps ||
+                       (line.p1 - other.p1).Length2() < eps)
+                      continue;
+                    mindist = min2(mindist, line.Dist(other));
+                  }
+                if(mindist == 1e99) continue;
+                mindist /= mparam.closeedgefac.value() + 1e-10;
+                if(mindist < 1e-3 * bounding_box.Diam())
+                  {
+                    (*testout) << "extremely small local h: " << mindist
+                               << " --> setting to " << 1e-3 * bounding_box.Diam() << endl;
+                    (*testout) << "somewhere near " << line.p0 << " - " << line.p1 << endl
+;
+                    mindist = 1e-3 * bounding_box.Diam();
+                  }
+                mesh.RestrictLocalHLine(line.p0, line.p1, mindist);
+              }
+          }
         multithread.task = savetask;
       }
+
+    for(const auto& mspnt : mparam.meshsize_points)
+      mesh.RestrictLocalH(mspnt.pnt, mspnt.h);
+
     mesh.LoadLocalMeshSize(mparam.meshsizefilename);
   }
 

libsrc/meshing/basegeom.hpp

@@ -40,7 +40,7 @@ namespace netgen
                               const EdgePointGeomInfo& gi2,
                               Point<3>& newp,
                               EdgePointGeomInfo& newgi) const = 0;
-    virtual Vec<3> GetTangent(const Point<3>& p) const = 0;
+    virtual Vec<3> GetTangent(double t) const = 0;
   };
 
   class GeometryFace
@@ -177,7 +177,7 @@ namespace netgen
                               int surfi2,
                               const EdgePointGeomInfo & egi) const
     {
-      return edges[egi.edgenr]->GetTangent(p);
+      throw Exception("Base geometry get tangent called");
     }
 
     virtual size_t GetEdgeIndex(const GeometryEdge& edge) const

libsrc/meshing/meshtype.cpp

@@ -2801,7 +2801,9 @@ namespace netgen
         << " elementorder = " <<  elementorder << endl
         << " quad = " <<  quad << endl
         << " inverttets = " <<  inverttets << endl
-        << " inverttrigs = " <<  inverttrigs << endl;
+        << " inverttrigs = " <<  inverttrigs << endl
+        << "closeedge enabled = " << closeedgefac.has_value() << endl
+        << "closeedgefac = " << (closeedgefac.has_value() ? closeedgefac.value() : 0.) << endl;
   }
 
   /*

libsrc/meshing/meshtype.hpp

@@ -1259,6 +1259,8 @@ namespace netgen
     double minh = 0.0;
     /// file for meshsize
     string meshsizefilename = "";
+    /// restrict h based on close edges
+    optional<double> closeedgefac = {};
     /// start surfacemeshing from everywhere in surface
     bool startinsurface = false;
     /// check overlapping surfaces (debug)

libsrc/meshing/python_mesh.hpp

@@ -176,6 +176,8 @@ inline void CreateMPfromKwargs(MeshingParameters& mp, py::kwargs kwargs, bool th
       mp.parallel_meshing = py::cast<bool>(kwargs.attr("pop")("parallel_meshing"));
     if(kwargs.contains("nthreads"))
       mp.nthreads = py::cast<int>(kwargs.attr("pop")("nthreads"));
+    if(kwargs.contains("closeedgefac"))
+      mp.closeedgefac = py::cast<optional<double>>(kwargs.attr("pop")("closeedgefac"));
     if(kwargs.size())
     {
       if(throw_if_not_all_parsed)

libsrc/occ/occgenmesh.cpp

@@ -1174,7 +1174,7 @@ namespace netgen
 
         // setting close edges
 
-        if (occparam.resthcloseedgeenable)
+        if (mparam.closeedgefac.has_value())
           {
             multithread.task = "Setting local mesh size (close edges)";
 
@@ -1259,7 +1259,7 @@ namespace netgen
                     mindist = min (mindist, line.Dist(lines[num]));
                   }
 
-                mindist /= (occparam.resthcloseedgefac + VSMALL);
+                mindist /= (mparam.closeedgefac.value() + VSMALL);
 
                 if (mindist < 1e-3 * bb.Diam())
                   {

libsrc/occ/occgeom.cpp

@@ -1897,8 +1897,6 @@ void STEP_GetEntityName(const TopoDS_Shape & theShape, STEPCAFControl_Reader * a
   void OCCParameters :: Print(ostream & ost) const
    {
       ost << "OCC Parameters:" << endl
-         << "close edges: " << resthcloseedgeenable
-         << ", fac = " << resthcloseedgefac << endl
 		 << "minimum edge length: " << resthminedgelenenable
 		 << ", min len = " << resthminedgelen << endl;
    }

libsrc/occ/occgeom.hpp

@@ -187,11 +187,11 @@ namespace netgen
   {
   public:
 
-    /// Factor for meshing close edges 
+    /// Factor for meshing close edges, moved to meshingparameters
-    double resthcloseedgefac = 2.;
+    // double resthcloseedgefac = 2.;
 
     /// Enable / Disable detection of close edges
-    int resthcloseedgeenable = true;
+    // int resthcloseedgeenable = true;
 
     /// Minimum edge length to be used for dividing edges to mesh points
     double resthminedgelen = 0.001;

libsrc/occ/occpkg.cpp

@@ -45,11 +45,7 @@ namespace netgen
 
     virtual void SetParameters (Tcl_Interp * interp) 
     {
-      occparam.resthcloseedgefac =
+      occparam.resthminedgelen = 
-	atof (Tcl_GetVar (interp, "::stloptions.resthcloseedgefac", 0));
-      occparam.resthcloseedgeenable =
-	atoi (Tcl_GetVar (interp, "::stloptions.resthcloseedgeenable", 0));
-	  occparam.resthminedgelen = 
 	atof (Tcl_GetVar (interp, "::stloptions.resthminedgelen", 0));
 	  occparam.resthminedgelenenable = 
 	atoi (Tcl_GetVar (interp, "::stloptions.resthminedgelenenable", 0));

libsrc/occ/python_occ.cpp

@@ -31,17 +31,6 @@ minedgelen: Optional[float] = 0.001
 
 void CreateOCCParametersFromKwargs(OCCParameters& occparam, py::dict kwargs)
 {
-  if(kwargs.contains("closeedgefac"))
-    {
-      auto val = kwargs.attr("pop")("closeedgefac");
-      if(val.is_none())
-        occparam.resthcloseedgeenable = false;
-      else
-        {
-          occparam.resthcloseedgefac = py::cast<double>(val);
-          occparam.resthcloseedgeenable = true;
-        }
-    }
   if(kwargs.contains("minedgelen"))
     {
       auto val = kwargs.attr("pop")("minedgelen");

libsrc/stlgeom/meshstlsurface.cpp

@@ -22,7 +22,7 @@ static void STLFindEdges (STLGeometry & geom, Mesh & mesh,
   // mark edge points:
   //int ngp = geom.GetNP();
 
-  geom.RestrictLocalH(mesh, h, stlparam);
+  geom.RestrictLocalH(mesh, h, stlparam, mparam);
   
   PushStatusF("Mesh Lines");
 

libsrc/stlgeom/python_stl.cpp

@@ -87,17 +87,6 @@ void CreateSTLParametersFromKwargs(STLParameters& stlparam, py::dict kwargs)
       stlparam.resthchartdistfac = py::cast<double>(val);
     }
   }
-  if(kwargs.contains("closeedgefac"))
-  {
-    auto val = kwargs.attr("pop")("closeedgefac");
-    if(val.is_none())
-      stlparam.resthcloseedgeenable = false;
-    else
-    {
-      stlparam.resthcloseedgeenable = true;
-      stlparam.resthcloseedgefac = py::cast<double>(val);
-    }
-  }
   if(kwargs.contains("edgeanglefac"))
   {
     auto val = kwargs.attr("pop")("edgeanglefac");

libsrc/stlgeom/stlgeom.hpp

@@ -465,7 +465,7 @@ namespace netgen
     int LineEndPointsSet() const {return lineendpoints.Size() == GetNP();}
     void ClearLineEndPoints();
 
-    DLL_HEADER void RestrictLocalH(class Mesh & mesh, double gh, const STLParameters& stlparam);
+    DLL_HEADER void RestrictLocalH(class Mesh & mesh, double gh, const STLParameters& stlparam, const MeshingParameters& mparam);
     void RestrictLocalHCurv(class Mesh & mesh, double gh, const STLParameters& stlparam);
     void RestrictHChartDistOneChart(ChartId chartnum, NgArray<int>& acttrigs, class Mesh & mesh, 
 				    double gh, double fact, double minh, const STLParameters& stlparam);

libsrc/stlgeom/stlgeommesh.cpp

@@ -814,7 +814,7 @@ void STLGeometry :: RestrictLocalHCurv(class Mesh & mesh, double gh, const STLPa
 
 }
   //restrict local h due to near edges and due to outer chart distance
-void STLGeometry :: RestrictLocalH(class Mesh & mesh, double gh, const STLParameters& stlparam)
+void STLGeometry :: RestrictLocalH(class Mesh & mesh, double gh, const STLParameters& stlparam, const MeshingParameters& mparam)
 {
   
   //bei jedem Dreieck alle Nachbardreiecke vergleichen, und, fallskein Kante dazwischen,
@@ -921,12 +921,12 @@ void STLGeometry :: RestrictLocalH(class Mesh & mesh, double gh, const STLParame
       PopStatus();
     }
 
-  if (stlparam.resthcloseedgeenable)
+  if (mparam.closeedgefac.has_value())
     {
       PushStatusF("Restrict H due to close edges");
       //geht nicht für spiralen!!!!!!!!!!!!!!!!!!
       
-      double disttohfact = sqr(10.0 / stlparam.resthcloseedgefac);
+      double disttohfact = sqr(10.0 / mparam.closeedgefac.value());
       int k,l;
       double h1, h2, dist;
       int rc = 0;

libsrc/stlgeom/stlpkg.cpp

@@ -78,11 +78,6 @@ namespace netgen
       stlparam.resthlinelengthenable =
 	atoi (Tcl_GetVar (interp, "::stloptions.resthlinelengthenable", 0));
 
-      stlparam.resthcloseedgefac =
-	atof (Tcl_GetVar (interp, "::stloptions.resthcloseedgefac", 0));
-      stlparam.resthcloseedgeenable =
-	atoi (Tcl_GetVar (interp, "::stloptions.resthcloseedgeenable", 0));
-
       stlparam.resthedgeanglefac =
 	atof (Tcl_GetVar (interp, "::stloptions.resthedgeanglefac", 0));
       stlparam.resthedgeangleenable =
@@ -538,7 +533,7 @@ namespace netgen
 	mesh -> SetLocalH (stlgeometry->GetBoundingBox().PMin() - Vec3d(10, 10, 10),
 			   stlgeometry->GetBoundingBox().PMax() + Vec3d(10, 10, 10),
 			   mparam.grading);
-	stlgeometry -> RestrictLocalH(*mesh, mparam.maxh, stlparam);
+	stlgeometry -> RestrictLocalH(*mesh, mparam.maxh, stlparam, mparam);
 
 	if (stlparam.resthsurfmeshcurvenable)
 	  mesh -> CalcLocalHFromSurfaceCurvature (mparam.grading, 

libsrc/stlgeom/stltool.cpp

@@ -1462,8 +1462,8 @@ STLParameters ::   STLParameters()
   resthchartdistenable = 1;
   resthlinelengthfac = 0.5;
   resthlinelengthenable = 1;
-  resthcloseedgefac = 1;
+  // resthcloseedgefac = 1;
-  resthcloseedgeenable = 1;
+  // resthcloseedgeenable = 1;
   resthedgeanglefac = 1;
   resthedgeangleenable = 0;
   resthsurfmeshcurvfac = 1;
@@ -1488,8 +1488,8 @@ void STLParameters :: Print (ostream & ost) const
       << ", fac = " << resthchartdistfac << endl
       << "line length: " << resthlinelengthenable
       << ", fac = " << resthlinelengthfac << endl
-      << "close edges: " << resthcloseedgeenable
+      // << "close edges: " << resthcloseedgeenable
-      << ", fac = " << resthcloseedgefac << endl
+      // << ", fac = " << resthcloseedgefac << endl
       << "edge angle: " << resthedgeangleenable
       << ", fac = " << resthedgeanglefac << endl;
 }

libsrc/stlgeom/stltool.hpp

@@ -282,8 +282,8 @@ public:
   double resthchartdistfac = 1.2;
   bool resthchartdistenable = true;
 
-  double resthcloseedgefac = 1.;
+  // double resthcloseedgefac = 1.;
-  bool resthcloseedgeenable = true;
+  // bool resthcloseedgeenable = true;
   
   double resthedgeanglefac = 1.;
   bool resthedgeangleenable = false;

ng/ngpkg.cpp

@@ -1283,6 +1283,10 @@ namespace netgen
 
     mparam.parallel_meshing = atoi (Tcl_GetVar (interp, "::options.parallel_meshing", 0));
     mparam.nthreads = atoi (Tcl_GetVar (interp, "::options.nthreads", 0));
+    if(atoi(Tcl_GetVar (interp, "::stloptions.resthcloseedgeenable", 0)))
+      mparam.closeedgefac = atof(Tcl_GetVar (interp, "::stloptions.resthcloseedgefac", 0));
+    else
+      mparam.closeedgefac = {};
 
     //BaseMoveableMem::totalsize = 0;
     // 1048576 * atoi (Tcl_GetVar (interp, "::options.memory", 0));

nglib/nglib.cpp

@@ -856,8 +856,8 @@ namespace nglib
 
       mp->Transfer_Parameters();
       
-      occparam.resthcloseedgeenable = mp->closeedgeenable;
+      if(mp->closeedgeenable)
-      occparam.resthcloseedgefac = mp->closeedgefact;
+        mparam.closeedgefac = mp->closeedgefact;
 
       // Delete the mesh structures in order to start with a clean 
       // slate