Merge branch 'cleanup_geometry' into 'master'

Cleanup geometry

See merge request jschoeberl/netgen!292

libsrc/csg/csgeom.cpp

@@ -72,11 +72,12 @@ namespace netgen
     Clean();
   }
 
-  void CSGeometry :: ProjectPoint(int surfind, Point<3> & p) const
+  PointGeomInfo CSGeometry :: ProjectPoint(int surfind, Point<3> & p) const
   {
     Point<3> hp = p;
     GetSurface(surfind)->Project (hp);
     p = hp;
+    return PointGeomInfo();
   }
 
   bool CSGeometry :: ProjectPointGI(int surfind, Point<3> & p, PointGeomInfo & gi) const
@@ -86,7 +87,7 @@ namespace netgen
   }
 
   void CSGeometry :: ProjectPointEdge(int surfind, INDEX surfind2,
-                                                Point<3> & p) const
+                                      Point<3> & p, EdgePointGeomInfo* /*unused*/) const
   {
     Point<3> hp = p;
     ProjectToEdge (GetSurface(surfind),
@@ -95,7 +96,8 @@ namespace netgen
   }
 
 
-  Vec<3> CSGeometry :: GetNormal(int surfind, const Point<3> & p) const
+  Vec<3> CSGeometry :: GetNormal(int surfind, const Point<3> & p,
+                                 const PointGeomInfo* /*unused*/) const
   {
     Vec<3> hn;
     GetSurface(surfind)->CalcGradient(p, hn);

libsrc/csg/csgeom.hpp

@@ -188,10 +188,12 @@ namespace netgen
 
     virtual void SaveToMeshFile (ostream & ost) const override;
 
-    void ProjectPoint(INDEX surfind, Point<3> & p) const override;
+    PointGeomInfo ProjectPoint(INDEX surfind, Point<3> & p) const override;
     bool ProjectPointGI (int surfind, Point<3> & p, PointGeomInfo & gi) const override;
-    void ProjectPointEdge(INDEX surfind, INDEX surfind2, Point<3> & p) const override;
+    void ProjectPointEdge(INDEX surfind, INDEX surfind2, Point<3> & p,
-    Vec<3> GetNormal(int surfind, const Point<3> & p) const override;
+                          EdgePointGeomInfo* gi = nullptr) const override;
+    Vec<3> GetNormal(int surfind, const Point<3> & p, const PointGeomInfo* gi = nullptr) const override;
+
     void PointBetween(const Point<3> & p1, const Point<3> & p2,
                       double secpoint, int surfi,
                       const PointGeomInfo & gi1,

libsrc/geom2d/geometry2d.cpp

@@ -86,7 +86,7 @@ namespace netgen
   }
 
   Vec<3> SplineGeometry2d :: GetNormal(int surfi1, const Point<3> & p,
-                                       const PointGeomInfo & gi) const
+                                       const PointGeomInfo* gi) const
   {
     return Vec<3> (0,0,1);
   }

libsrc/geom2d/geometry2d.hpp

@@ -182,7 +182,7 @@ namespace netgen
     Vec<3> GetTangent (const Point<3> & p, int surfi1, int surfi2,
                        const EdgePointGeomInfo & ap1) const override;
     Vec<3> GetNormal(int surfi1, const Point<3> & p,
-                     const PointGeomInfo & gi) const override;
+                     const PointGeomInfo* gi) const override;
 
     const SplineSegExt & GetSpline (const int i) const 
     { 

libsrc/meshing/basegeom.cpp

@@ -10,14 +10,6 @@ namespace netgen
   GeometryRegister :: ~GeometryRegister()
   { ; }
 
-  Array<Point<3>> GeometryEdge :: GetEquidistantPointArray(size_t npoints) const
-  {
-    Array<Point<3>> pts(npoints);
-    for(auto i : Range(npoints))
-      pts[i] = GetPoint(double(i)/(npoints-1));
-    return pts;
-  }
-
   void GeometryFace :: RestrictHTrig(Mesh& mesh,
                                      const PointGeomInfo& gi0,
                                      const PointGeomInfo& gi1,
@@ -103,16 +95,23 @@ namespace netgen
     mesh.SetLocalH(bounding_box.PMin(), bounding_box.PMax(),
                    mparam.grading);
 
+    // only set meshsize for edges longer than this
+    double mincurvelength = 1e-3 * bounding_box.Diam();
+
     if(mparam.uselocalh)
       {
         double eps = 1e-12 * bounding_box.Diam();
+        const char* savetask = multithread.task;
+        multithread.task = "Analyse Edges";
 
         // restrict meshsize on edges
-        for(const auto & edge : edges)
+        for(auto i : Range(edges))
           {
+            multithread.percent = 100. * i/edges.Size();
+            const auto & edge = edges[i];
             auto length = edge->GetLength();
             // skip very short edges
-            if(length < eps)
+            if(length < mincurvelength)
               continue;
             static constexpr int npts = 20;
             // restrict mesh size based on edge length
@@ -135,9 +134,15 @@ namespace netgen
               }
           }
 
+        multithread.task = "Analyse Faces";
         // restrict meshsize on faces
-        for(const auto& face : faces)
+        for(auto i : Range(faces))
-          face->RestrictH(mesh, mparam);
+          {
+            multithread.percent = 100. * i/faces.Size();
+            const auto& face = faces[i];
+            face->RestrictH(mesh, mparam);
+          }
+        multithread.task = savetask;
       }
     mesh.LoadLocalMeshSize(mparam.meshsizefilename);
   }
@@ -148,6 +153,8 @@ namespace netgen
     static Timer t1("MeshEdges"); RegionTimer regt(t1);
     static Timer tdivide("Divide Edges");
     static Timer tdivedgesections("Divide edge sections");
+    const char* savetask = multithread.task;
+    multithread.task = "Mesh Edges";
 
     // create face descriptors and set bc names
     mesh.SetNBCNames(faces.Size());
@@ -177,6 +184,7 @@ namespace netgen
             auto boundary = face.GetBoundary(facebndnr);
             for(auto enr : Range(boundary))
               {
+                multithread.percent = 100. * ((double(enr)/boundary.Size() + facebndnr)/face.GetNBoundaries() + facenr)/faces.Size();
                 const auto& oriented_edge = *boundary[enr];
                 auto edgenr = GetEdgeIndex(oriented_edge);
                 const auto& edge = edges[edgenr];
@@ -196,11 +204,15 @@ namespace netgen
                 double hvalue[divide_edge_sections+1];
                 hvalue[0] = 0;
 
+                Point<3> old_pt = edge->GetPoint(0.);
                 // calc local h for edge
                 tdivedgesections.Start();
-                auto edgepts = edge->GetEquidistantPointArray(divide_edge_sections+1);
+                for(auto i : Range(divide_edge_sections))
-                for(auto i : Range(edgepts.Size()-1))
+                  {
-                  hvalue[i+1] = hvalue[i] + 1./mesh.GetH(edgepts[i+1]) * (edgepts[i+1]-edgepts[i]).Length();
+                    auto pt = edge->GetPoint(double(i+1)/divide_edge_sections);
+                    hvalue[i+1] = hvalue[i] + 1./mesh.GetH(pt) * (pt-old_pt).Length();
+                    old_pt = pt;
+                  }
                 int nsubedges = max2(1, int(floor(hvalue[divide_edge_sections]+0.5)));
                 tdivedgesections.Stop();
                 mps.SetSize(nsubedges-1);
@@ -235,7 +247,7 @@ namespace netgen
                     cout << "CORRECTED" << endl;
                     mps.SetSize (nsubedges-2);
                     params.SetSize (nsubedges);
-                    params[nsubedges] = 1.;
+                    params[nsubedges-1] = 1.;
                   }
                 tdivide.Stop();
                 // ----------- Add Points to mesh and create segments -----
@@ -292,16 +304,21 @@ namespace netgen
               }
           }
       }
+    mesh.CalcSurfacesOfNode();
+    multithread.task = savetask;
   }
 
   void NetgenGeometry :: MeshSurface(Mesh& mesh,
                                      const MeshingParameters& mparam) const
   {
     static Timer t1("Surface Meshing"); RegionTimer regt(t1);
+    const char* savetask = multithread.task;
+    multithread.task = "Mesh Surface";
 
     Array<int, PointIndex> glob2loc(mesh.GetNP());
     for(auto k : Range(faces))
       {
+        multithread.percent = 100. * k/faces.Size();
         const auto& face = *faces[k];
         auto bb = face.GetBoundingBox();
         bb.Increase(bb.Diam()/10);
@@ -354,6 +371,7 @@ namespace netgen
             mesh.SurfaceElements()[i].SetIndex(k+1);
           }
       }
+    multithread.task = savetask;
   }
 
   void NetgenGeometry :: OptimizeSurface(Mesh& mesh, const MeshingParameters& mparam) const
@@ -366,9 +384,11 @@ namespace netgen
     auto meshopt = MeshOptimize2d(mesh);
     for(auto i : Range(mparam.optsteps2d))
       {
-        PrintMessage(2, "Optimization step ", i);
+        PrintMessage(3, "Optimization step ", i);
+        int innerstep = 0;
         for(auto optstep : mparam.optimize2d)
           {
+            multithread.percent = 100. * (double(innerstep++)/mparam.optimize2d.size() + i)/mparam.optsteps2d;
             switch(optstep)
               {
               case 's':

libsrc/meshing/basegeom.hpp

@@ -32,7 +32,15 @@ namespace netgen
     virtual double CalcStep(double t, double sag) const = 0;
     virtual bool OrientedLikeGlobal() const = 0;
     virtual size_t GetHash() const = 0;
-    virtual Array<Point<3>> GetEquidistantPointArray(size_t npoints) const;
+    virtual void ProjectPoint(Point<3>& p, EdgePointGeomInfo* gi) const = 0;
+    virtual void PointBetween(const Point<3>& p1,
+                              const Point<3>& p2,
+                              double secpoint,
+                              const EdgePointGeomInfo& gi1,
+                              const EdgePointGeomInfo& gi2,
+                              Point<3>& newp,
+                              EdgePointGeomInfo& newgi) const = 0;
+    virtual Vec<3> GetTangent(const Point<3>& p) const = 0;
   };
 
   class GeometryFace
@@ -42,9 +50,11 @@ namespace netgen
     virtual size_t GetNBoundaries() const = 0;
     virtual Array<unique_ptr<GeometryEdge>> GetBoundary(size_t index) const = 0;
     virtual string GetName() const { return "default"; }
+    virtual PointGeomInfo Project(Point<3>& p) const = 0;
     // Project point using geo info. Fast if point is close to
     // parametrization in geo info.
     virtual bool ProjectPointGI(Point<3>& p, PointGeomInfo& gi) const =0;
+    virtual bool CalcPointGeomInfo(const Point<3>& p, PointGeomInfo& gi) const = 0;
     virtual Point<3> GetPoint(const PointGeomInfo& gi) const = 0;
     virtual void CalcEdgePointGI(const GeometryEdge& edge,
                                  double t,
@@ -55,6 +65,15 @@ namespace netgen
     virtual double GetCurvature(const PointGeomInfo& gi) const = 0;
 
     virtual void RestrictH(Mesh& mesh, const MeshingParameters& mparam) const = 0;
+    virtual Vec<3> GetNormal(const Point<3>& p, const PointGeomInfo* gi = nullptr) const = 0;
+
+    virtual void PointBetween(const Point<3>& p1,
+                              const Point<3>& p2,
+                              double secpoint,
+                              const PointGeomInfo& gi1,
+                              const PointGeomInfo& gi2,
+                              Point<3>& newp,
+                              PointGeomInfo& newgi) const = 0;
 
   protected:
     void RestrictHTrig(Mesh& mesh,
@@ -99,27 +118,27 @@ namespace netgen
 
     virtual void FinalizeMesh(Mesh& mesh) const {}
 
-    virtual void ProjectPoint (int surfind, Point<3> & p) const
+    virtual PointGeomInfo ProjectPoint (int surfind, Point<3> & p) const
-    { }
+    {
-    virtual void ProjectPointEdge (int surfind, int surfind2, Point<3> & p) const { }
+      return faces[surfind-1]->Project(p);
-  virtual void ProjectPointEdge (int surfind, int surfind2, Point<3> & p, EdgePointGeomInfo& gi) const
+    }
-  { ProjectPointEdge(surfind, surfind2, p); }
 
-    virtual bool CalcPointGeomInfo(int surfind, PointGeomInfo& gi, const Point<3> & p3) const {return false;}
+  virtual void ProjectPointEdge (int surfind, int surfind2, Point<3> & p, EdgePointGeomInfo* gi = nullptr) const
+  {
+    edges[gi->edgenr]->ProjectPoint(p, gi);
+  }
+
+    virtual bool CalcPointGeomInfo(int surfind, PointGeomInfo& gi, const Point<3> & p3) const
+    {
+      return faces[surfind-1]->CalcPointGeomInfo(p3, gi);
+    }
     virtual bool ProjectPointGI (int surfind, Point<3> & p, PointGeomInfo & gi) const
     {
-      throw Exception("ProjectPointGI not overloaded in class" + Demangle(typeid(*this).name()));
+      return faces[surfind-1]->ProjectPointGI(p, gi);
     }
 
-    virtual Vec<3> GetNormal(int surfind, const Point<3> & p) const
+    virtual Vec<3> GetNormal(int surfind, const Point<3> & p, const PointGeomInfo* gi = nullptr) const
-    { return {0.,0.,1.}; }
+    { return faces[surfind-1]->GetNormal(p, gi); }
-    virtual Vec<3> GetNormal(int surfind, const Point<3> & p, const PointGeomInfo & gi) const
-    { return GetNormal(surfind, p); }
-    [[deprecated]]
-    void GetNormal(int surfind, const Point<3> & p, Vec<3> & n) const
-    {
-      n = GetNormal(surfind, p);
-    }
 
     virtual void PointBetween (const Point<3> & p1,
                                const Point<3> & p2, double secpoint,
@@ -129,6 +148,11 @@ namespace netgen
                                Point<3> & newp,
                                PointGeomInfo & newgi) const
     {
+      if(faces.Size())
+        {
+          faces[surfi-1]->PointBetween(p1, p2, secpoint, gi1, gi2, newp, newgi);
+          return;
+        }
       newp = p1 + secpoint * (p2-p1);
     }
 
@@ -140,13 +164,21 @@ namespace netgen
                                   Point<3> & newp,
                                   EdgePointGeomInfo & newgi) const
     {
+      if(edges.Size())
+        {
+          edges[ap1.edgenr]->PointBetween(p1, p2, secpoint,
+                                          ap1, ap2, newp, newgi);
+          return;
+        }
       newp = p1+secpoint*(p2-p1);
     }
 
     virtual Vec<3> GetTangent(const Point<3> & p, int surfi1,
                               int surfi2,
                               const EdgePointGeomInfo & egi) const
-    { throw Exception("Call GetTangent of " + Demangle(typeid(*this).name())); }
+    {
+      return edges[egi.edgenr]->GetTangent(p);
+    }
 
     virtual size_t GetEdgeIndex(const GeometryEdge& edge) const
     {

libsrc/meshing/curvedelems.cpp

@@ -839,8 +839,8 @@ namespace netgen
 		{
 		  Point<3> pm = Center (p1, p2);
 
-		  Vec<3> n1 = geo.GetNormal (surfnr[e], p1, gi0[e]);
+		  Vec<3> n1 = geo.GetNormal (surfnr[e], p1, &gi0[e]);
-		  Vec<3> n2 = geo.GetNormal (surfnr[e], p2, gi1[e]);
+		  Vec<3> n2 = geo.GetNormal (surfnr[e], p2, &gi1[e]);
 
 		  // p3 = pm + alpha1 n1 + alpha2 n2
 		
@@ -1084,7 +1084,7 @@ namespace netgen
 		  v05 /= 1 + (w-1) * 0.5;
 		  Point<3> p05 (v05), pp05(v05);
 		  geo.ProjectPointEdge(edge_surfnr1[edgenr], edge_surfnr2[edgenr], pp05,
-                                       edge_gi0[edgenr]);
+                                       &edge_gi0[edgenr]);
 		  double d = Dist (pp05, p05);
 
 		  if (d < dold)

libsrc/meshing/improve2.cpp

@@ -85,11 +85,11 @@ namespace netgen
     nv1.Normalize();
     nv2.Normalize();
 
-    auto nvp3 = geo.GetNormal (surfnr, mesh.Point(pi3), gi3);
+    auto nvp3 = geo.GetNormal (surfnr, mesh.Point(pi3), &gi3);
 
     nvp3.Normalize();
 
-    auto nvp4 = geo.GetNormal (surfnr, mesh.Point(pi4), gi4);
+    auto nvp4 = geo.GetNormal (surfnr, mesh.Point(pi4), &gi4);
 
     nvp4.Normalize();
 
@@ -648,7 +648,7 @@ namespace netgen
                     {
                       const int faceindex = hel.GetIndex();
                       const int surfnr = mesh.GetFaceDescriptor (faceindex).SurfNr();
-                      normals[pi] = geo.GetNormal (surfnr, mesh[pi], hel.GeomInfoPi(k+1));
+                      normals[pi] = geo.GetNormal (surfnr, mesh[pi], &hel.GeomInfoPi(k+1));
                       break;
                     }
               }

libsrc/meshing/improve2gen.cpp

@@ -397,7 +397,7 @@ namespace netgen
 		// calc metric badness
 		double bad1 = 0, bad2 = 0;
 		// SelectSurfaceOfPoint (mesh.Point(pmap.Get(1)), pgi.Get(1));
-		auto n = geo.GetNormal(surfnr, mesh.Point(pmap.Get(1)), pgi.Elem(1));
+		auto n = geo.GetNormal(surfnr, mesh.Point(pmap.Get(1)), &pgi.Elem(1));
 		  
 		for (int j = 0; j < rule.oldels.Size(); j++)
 		  bad1 += mesh[elmap[j]].CalcJacobianBadness (mesh.Points(), n);

libsrc/meshing/improve3.cpp

@@ -183,7 +183,7 @@ void MeshOptimize3d :: CombineImproveSequential (Mesh & mesh,
 
   //  mesh.CalcSurfacesOfNode ();
   const char * savetask = multithread.task;
-  multithread.task = "Combine Improve";
+  multithread.task = "Optimize Volume: Combine Improve";
 
 
   double totalbad = 0;
@@ -435,7 +435,7 @@ void MeshOptimize3d :: CombineImprove (Mesh & mesh,
 
   //  mesh.CalcSurfacesOfNode ();
   const char * savetask = multithread.task;
-  multithread.task = "Combine Improve";
+  multithread.task = "Optimize Volume: Combine Improve";
 
 
   tbad.Start();
@@ -712,7 +712,7 @@ void MeshOptimize3d :: SplitImprove (Mesh & mesh,
   Array<double> elerrs(ne);
 
   const char * savetask = multithread.task;
-  multithread.task = "Split Improve";
+  multithread.task = "Optimize Volume: Split Improve";
 
   PrintMessage (3, "SplitImprove");
   (*testout)  << "start SplitImprove" << "\n";
@@ -826,7 +826,7 @@ void MeshOptimize3d :: SplitImproveSequential (Mesh & mesh,
   illegaltet.Clear();
 
   const char * savetask = multithread.task;
-  multithread.task = "Split Improve";
+  multithread.task = "Optimize Volume: Split Improve";
 
   PrintMessage (3, "SplitImprove");
   (*testout)  << "start SplitImprove" << "\n";
@@ -1121,7 +1121,7 @@ void MeshOptimize3d :: SwapImproveSequential (Mesh & mesh, OPTIMIZEGOAL goal,
   (*testout) << "\n" << "Start SwapImprove" << endl;
 
   const char * savetask = multithread.task;
-  multithread.task = "Swap Improve";
+  multithread.task = "Optimize Volume: Swap Improve";
   
   //  mesh.CalcSurfacesOfNode ();
     
@@ -2617,7 +2617,7 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
   (*testout) << "\n" << "Start SwapImprove" << endl;
 
   const char * savetask = multithread.task;
-  multithread.task = "Swap Improve";
+  multithread.task = "Optimize Volume: Swap Improve";
 
   INDEX_3_HASHTABLE<int> faces(mesh.GetNOpenElements()/3 + 2);
   if (goal == OPT_CONFORM)

libsrc/meshing/meshfunc.cpp

@@ -646,6 +646,8 @@ namespace netgen
   {
     static Timer t("OptimizeVolume"); RegionTimer reg(t);
     RegionTaskManager rtm(mp.parallel_meshing ? mp.nthreads : 0);
+    const char* savetask = multithread.task;
+    multithread.task = "Optimize Volume";
     
     int i;
 
@@ -663,7 +665,7 @@ namespace netgen
     */
 
     mesh3d.CalcSurfacesOfNode();
-    for (i = 1; i <= mp.optsteps3d; i++)
+    for (auto i : Range(mp.optsteps3d))
       {
 	if (multithread.terminate)
 	  break;
@@ -672,12 +674,13 @@ namespace netgen
 
 	// teterrpow = mp.opterrpow;
 	// for (size_t j = 1; j <= strlen(mp.optimize3d); j++)
-        for (size_t j = 1; j <= mp.optimize3d.length(); j++)
+        for (auto j : Range(mp.optimize3d.size()))
 	  {
+            multithread.percent = 100.* (double(j)/mp.optimize3d.size() + i)/mp.optsteps3d;
 	    if (multithread.terminate)
 	      break;
 
-	    switch (mp.optimize3d[j-1])
+	    switch (mp.optimize3d[j])
 	      {
 	      case 'c': optmesh.CombineImprove(mesh3d, OPT_REST); break;
 	      case 'd': optmesh.SplitImprove(mesh3d); break;
@@ -698,6 +701,7 @@ namespace netgen
 	MeshQuality3d (mesh3d);
       }
   
+    multithread.task = savetask;
     return MESHING3_OK;
   }
 

libsrc/meshing/meshing2.cpp

@@ -142,8 +142,8 @@ namespace netgen
   {
     p1 = ap1;
     p2 = ap2;
-    auto n1 = geo.GetNormal(gi1->trignum, p1, *gi1);
+    auto n1 = geo.GetNormal(gi1->trignum, p1, gi1);
-    auto n2 = geo.GetNormal(gi2->trignum, p2, *gi2);
+    auto n2 = geo.GetNormal(gi2->trignum, p2, gi2);
 
     ez = 0.5 * (n1+n2);
     ez.Normalize();
@@ -158,7 +158,7 @@ namespace netgen
 				     Point<2> & plainpoint, double h, int & zone)
   {
     auto& gi = geominfo.GetPGI(1);
-    auto n = geo.GetNormal(gi.trignum, locpoint, gi);
+    auto n = geo.GetNormal(gi.trignum, locpoint, &gi);
     auto p1p = locpoint - p1;
     plainpoint(0) = (p1p * ex) / h;
     plainpoint(1) = (p1p * ey) / h;
@@ -176,7 +176,7 @@ namespace netgen
   {
     locpoint = p1 + (h*plainpoint(0)) * ex + (h* plainpoint(1)) * ey;
     if (!geo.ProjectPointGI(gi.trignum, locpoint, gi))
-      geo.ProjectPoint(gi.trignum, locpoint);
+      gi = geo.ProjectPoint(gi.trignum, locpoint);
     return 0;
   }
 

libsrc/meshing/smoothing2.cpp

@@ -218,7 +218,7 @@ namespace netgen
     {
       double badness = 0;
       
-      auto n = geo.GetNormal(ld.surfi, ld.sp1, ld.gi1);
+      auto n = geo.GetNormal(ld.surfi, ld.sp1, &ld.gi1);
       Point<3> pp1 = ld.sp1 + x(0) * ld.t1 + x(1) * ld.t2;
       
       for (int j = 0; j < ld.locelements.Size(); j++)
@@ -359,7 +359,7 @@ namespace netgen
     vgrad = 0;
     double badness = 0;
 
-    auto n = geo.GetNormal(ld.surfi, ld.sp1, ld.gi1);
+    auto n = geo.GetNormal(ld.surfi, ld.sp1, &ld.gi1);
     pp1 = ld.sp1 + x(0) * ld.t1 + x(1) * ld.t2;
 
     //  meshthis -> ProjectPoint (surfi, pp1);
@@ -414,7 +414,7 @@ namespace netgen
     vgrad = 0;
     double badness = 0;
 
-    auto n = geo.GetNormal(ld.surfi, ld.sp1, ld.gi1);
+    auto n = geo.GetNormal(ld.surfi, ld.sp1, &ld.gi1);
     pp1 = ld.sp1 + x(0) * ld.t1 + x(1) * ld.t2;
 
     for (int j = 0; j < ld.locelements.Size(); j++)
@@ -577,7 +577,7 @@ namespace netgen
     vgrad = 0;
     badness = 0;
 
-    auto n = geo.GetNormal(ld.surfi, ld.sp1, ld.gi1);
+    // auto n = geo.GetNormal(ld.surfi, ld.sp1, &ld.gi1);
 
     pp1 = ld.sp1 + x(0) * ld.t1 + x(1) * ld.t2;
 
@@ -933,7 +933,7 @@ namespace netgen
 	      }
 
 
-          ld.normal = geo.GetNormal(ld.surfi, ld.sp1, ld.gi1);
+          ld.normal = geo.GetNormal(ld.surfi, ld.sp1, &ld.gi1);
 	  ld.t1 = ld.normal.GetNormal ();
 	  ld.t2 = Cross (ld.normal, ld.t1);
 	  

libsrc/meshing/smoothing3.cpp

@@ -1164,7 +1164,7 @@ void Mesh :: ImproveMesh (const CSG eometry & geometry, OPTIMIZEGOAL goal)
     }
 
   const char * savetask = multithread.task;
-  multithread.task = "Smooth Mesh";
+  multithread.task = "Optimize Volume: Smooth Mesh";
 
 
   TABLE<INDEX> surfelementsonpoint(points.Size());
@@ -1398,7 +1398,7 @@ void Mesh :: ImproveMeshSequential (const MeshingParameters & mp, OPTIMIZEGOAL g
 
 
   const char * savetask = multithread.task;
-  multithread.task = "Smooth Mesh";
+  multithread.task = "Optimize Volume: Smooth Mesh";
 
   for (PointIndex pi : points.Range())
     if ( (*this)[pi].Type() == INNERPOINT )
@@ -1524,7 +1524,7 @@ void Mesh :: ImproveMesh (const MeshingParameters & mp, OPTIMIZEGOAL goal)
     }
 
   const char * savetask = multithread.task;
-  multithread.task = "Smooth Mesh";
+  multithread.task = "Optimize Volume: Smooth Mesh";
 
   topt.Start();
   int counter = 0;
@@ -1659,7 +1659,7 @@ void Mesh :: ImproveMeshJacobian (const MeshingParameters & mp,
 
 
   const char * savetask = multithread.task;
-  multithread.task = "Smooth Mesh Jacobian";
+  multithread.task = "Optimize Volume: Smooth Mesh Jacobian";
   
   // for (PointIndex pi = points.Begin(); i < points.End(); pi++)
   for (PointIndex pi : points.Range())
@@ -1815,7 +1815,7 @@ void Mesh :: ImproveMeshJacobianOnSurface (const MeshingParameters & mp,
 
 
   const char * savetask = multithread.task;
-  multithread.task = "Smooth Mesh Jacobian";
+  multithread.task = "Optimize Volume: Smooth Mesh Jacobian";
   
   // for (PointIndex pi = points.Begin(); pi <= points.End(); pi++)
   for (PointIndex pi : points.Range())

libsrc/occ/occgeom.cpp

@@ -1034,7 +1034,7 @@ void STEP_GetEntityName(const TopoDS_Shape & theShape, STEPCAFControl_Reader * a
       SetCenter();
    }
 
-   void OCCGeometry :: ProjectPoint(int surfi, Point<3> & p) const
+   PointGeomInfo OCCGeometry :: ProjectPoint(int surfi, Point<3> & p) const
    {
       static int cnt = 0;
       if (++cnt % 1000 == 0) cout << "Project cnt = " << cnt << endl;
@@ -1048,9 +1048,12 @@ void STEP_GetEntityName(const TopoDS_Shape & theShape, STEPCAFControl_Reader * a
       suval.Coord( u, v);
       pnt = thesurf->Value( u, v );
 
-
+      PointGeomInfo gi;
+      gi.trignum = surfi;
+      gi.u = u;
+      gi.v = v;
       p = Point<3> (pnt.X(), pnt.Y(), pnt.Z());
-
+      return gi;
    }
 
   bool OCCGeometry :: ProjectPointGI(int surfind, Point<3>& p, PointGeomInfo& gi) const
@@ -1069,7 +1072,7 @@ void STEP_GetEntityName(const TopoDS_Shape & theShape, STEPCAFControl_Reader * a
   }
 
   void OCCGeometry :: ProjectPointEdge(int surfind, INDEX surfind2,
-						   Point<3> & p) const
+                                       Point<3> & p, EdgePointGeomInfo* gi) const
   {
     TopExp_Explorer exp0, exp1;
     bool done = false;
@@ -1151,26 +1154,25 @@ void STEP_GetEntityName(const TopoDS_Shape & theShape, STEPCAFControl_Reader * a
       return true;
    }
 
-  Vec<3> OCCGeometry :: GetNormal(int surfind, const Point<3> & p, const PointGeomInfo & geominfo) const
+  Vec<3> OCCGeometry :: GetNormal(int surfind, const Point<3> & p, const PointGeomInfo* geominfo) const
   {
-    gp_Pnt pnt;
+    if(geominfo)
-    gp_Vec du, dv;
+      {
+        gp_Pnt pnt;
+        gp_Vec du, dv;
 
-    Handle(Geom_Surface) occface;
+        Handle(Geom_Surface) occface;
-    occface = BRep_Tool::Surface(TopoDS::Face(fmap(surfind)));
+        occface = BRep_Tool::Surface(TopoDS::Face(fmap(surfind)));
 
-    occface->D1(geominfo.u,geominfo.v,pnt,du,dv);
+        occface->D1(geominfo->u,geominfo->v,pnt,du,dv);
 
-    auto n = Cross (Vec<3>(du.X(), du.Y(), du.Z()),
+        auto n = Cross (Vec<3>(du.X(), du.Y(), du.Z()),
-	       Vec<3>(dv.X(), dv.Y(), dv.Z()));
+                        Vec<3>(dv.X(), dv.Y(), dv.Z()));
-    n.Normalize();
+        n.Normalize();
 
-    if (fmap(surfind).Orientation() == TopAbs_REVERSED) n *= -1;
+        if (fmap(surfind).Orientation() == TopAbs_REVERSED) n *= -1;
-    return n;
+        return n;
-  }
+      }
-
-  Vec<3> OCCGeometry :: GetNormal(int surfind, const Point<3> & p) const
-  {
     Standard_Real u,v;
 
     gp_Pnt pnt(p(0), p(1), p(2));

libsrc/occ/occgeom.hpp

@@ -280,11 +280,11 @@ namespace netgen
      
     void DoArchive(Archive& ar) override;
 
-    void ProjectPoint(int surfind, Point<3> & p) const override;
+    PointGeomInfo ProjectPoint(int surfind, Point<3> & p) const override;
-    void ProjectPointEdge (int surfind, int surfind2, Point<3> & p) const override;
+    void ProjectPointEdge (int surfind, int surfind2, Point<3> & p,
+                           EdgePointGeomInfo* gi = nullptr) const override;
     bool ProjectPointGI (int surfind, Point<3> & p, PointGeomInfo & gi) const override;
-    Vec<3> GetNormal(int surfind, const Point<3> & p) const override;
+    Vec<3> GetNormal(int surfind, const Point<3> & p, const PointGeomInfo* gi) const override;
-    Vec<3> GetNormal(int surfind, const Point<3> & p, const PointGeomInfo & gi) const override;
     bool CalcPointGeomInfo(int surfind, PointGeomInfo& gi, const Point<3> & p3) const override;
 
     void PointBetweenEdge(const Point<3> & p1, const Point<3> & p2, double secpoint,

libsrc/stlgeom/stlgeom.cpp

@@ -100,14 +100,11 @@ int STLGeometry :: GenerateMesh (shared_ptr<Mesh> & mesh, MeshingParameters & mp
   return STLMeshingDummy (this, mesh, mparam, stlpar);
 }
 
-Vec<3> STLGeometry :: GetNormal(INDEX surfind, const Point<3> & p) const
+Vec<3> STLGeometry :: GetNormal(int surfind, const Point<3> & p, const PointGeomInfo* gi) const
 {
-  throw Exception("STLGeometry::GetNormal without PointGeomInfo called");
+  if(!gi)
-}
+    throw Exception("STLGeometry::GetNormal without PointGeomInfo called");
-
+  return GetChart(GetChartNr(gi->trignum)).GetNormal();
-Vec<3> STLGeometry :: GetNormal(int surfind, const Point<3> & p, const PointGeomInfo & gi) const
-{
-  return GetChart(GetChartNr(gi.trignum)).GetNormal();
 }
 
 bool STLGeometry :: CalcPointGeomInfo(int /*surfind*/, PointGeomInfo& gi, const Point<3> & p3) const
@@ -144,7 +141,7 @@ bool STLGeometry :: ProjectPointGI (int surfind, Point<3> & p, PointGeomInfo & g
   return true;
 }
 
-void STLGeometry :: ProjectPoint (INDEX surfind, Point<3> & p) const
+PointGeomInfo STLGeometry :: ProjectPoint (INDEX surfind, Point<3> & p) const
 {
   throw Exception("ProjectPoint without PointGeomInfo not implemented");
 }

libsrc/stlgeom/stlgeom.hpp

@@ -193,10 +193,9 @@ namespace netgen
     virtual void Save (string filename) const override;
 
     bool CalcPointGeomInfo(int surfind, PointGeomInfo& gi, const Point<3> & p3) const override;
-    void ProjectPoint(INDEX surfind, Point<3> & p) const override;
+    PointGeomInfo ProjectPoint(INDEX surfind, Point<3> & p) const override;
     bool ProjectPointGI (int surfind, Point<3> & p, PointGeomInfo & gi) const override;
-    Vec<3> GetNormal(int surfind, const Point<3> & p) const override;
+    Vec<3> GetNormal(int surfind, const Point<3> & p, const PointGeomInfo* gi = nullptr) const override;
-    Vec<3> GetNormal(int surfind, const Point<3> & p, const PointGeomInfo & gi) const override;
     void PointBetween(const Point<3> & p1, const Point<3> & p2,
                       double secpoint, int surfi,
                       const PointGeomInfo & gi1,