Merge branch 'refactor_geo_meshing' into 'master'

Refactor geo meshing

See merge request jschoeberl/netgen!286

libsrc/csg/genmesh.cpp

@@ -422,7 +422,7 @@ namespace netgen
 	  geom.GetSurface((mesh.GetFaceDescriptor(k).SurfNr()));
 
 
-	Meshing2Surfaces meshing(*surf, mparam, geom.BoundingBox());
+	Meshing2Surfaces meshing(geom, *surf, mparam, geom.BoundingBox());
 	meshing.SetStartTime (starttime);
 
         double eps = 1e-8 * geom.MaxSize();

libsrc/csg/meshsurf.cpp

@@ -14,13 +14,14 @@ Meshing2Surfaces :: Meshing2Surfaces (const Surface & asurface)
   ;
 }
   */
-Meshing2Surfaces :: Meshing2Surfaces (const Surface & asurf,
+  Meshing2Surfaces :: Meshing2Surfaces (const CSGeometry& geo,
+                                        const Surface & asurf,
                                         const MeshingParameters & mp,
                                         const Box<3> & abb)
-  : Meshing2(mp, abb), surface(asurf), mparam (mp)
+    : Meshing2(geo, mp, abb), surface(asurf), mparam (mp)
-{
+  {
     ;
-}
+  }
 
 
 void Meshing2Surfaces :: DefineTransformation (const Point<3> & p1, const Point<3> & p2,

libsrc/csg/meshsurf.hpp

@@ -16,7 +16,9 @@ namespace netgen
     ///
     //  Meshing2Surfaces (const Surface & asurf);
     ///
-    Meshing2Surfaces (const Surface & asurf, const MeshingParameters & mp, 
+    Meshing2Surfaces (const CSGeometry& geo,
+                      const Surface & asurf,
+                      const MeshingParameters & mp,
 		      const Box<3> & aboundingbox);
 
   protected:

libsrc/geom2d/genmesh2d.cpp

@@ -565,7 +565,7 @@ namespace netgen
 
         mp.quad = hquad || geometry.GetDomainQuadMeshing (domnr);
 
-	Meshing2 meshing (mp, Box<3> (pmin, pmax));
+	Meshing2 meshing (geometry, mp, Box<3> (pmin, pmax));
 
 	NgArray<int, PointIndex::BASE> compress(bnp);
 	compress = -1;

libsrc/meshing/CMakeLists.txt

@@ -37,5 +37,6 @@ install(FILES
   localh.hpp meshclass.hpp meshfunc.hpp meshing2.hpp meshing3.hpp
   meshing.hpp meshtool.hpp meshtype.hpp msghandler.hpp paralleltop.hpp
   ruler2.hpp ruler3.hpp specials.hpp topology.hpp validate.hpp
+  python_mesh.hpp
   DESTINATION ${NG_INSTALL_DIR_INCLUDE}/meshing COMPONENT netgen_devel
 )

libsrc/meshing/basegeom.cpp

@@ -10,7 +10,89 @@ namespace netgen
   GeometryRegister :: ~GeometryRegister()
   { ; }
 
-  void NetgenGeometry :: OptimizeSurface(Mesh& mesh, const MeshingParameters& mparam)
+  void NetgenGeometry :: Analyse(Mesh& mesh,
+                                 const MeshingParameters& mparam) const
+  {
+    static Timer t1("SetLocalMeshsize"); RegionTimer regt(t1);
+    mesh.SetGlobalH(mparam.maxh);
+    mesh.SetMinimalH(mparam.minh);
+
+    mesh.SetLocalH(bounding_box.PMin(), bounding_box.PMax(),
+                   mparam.grading);
+
+    if(mparam.uselocalh)
+      RestrictLocalMeshsize(mesh, mparam);
+    mesh.LoadLocalMeshSize(mparam.meshsizefilename);
+  }
+
+  void NetgenGeometry :: FindEdges(Mesh& mesh,
+                                   const MeshingParameters& mparam) const
+  {
+  }
+
+  void NetgenGeometry :: MeshSurface(Mesh& mesh,
+                                     const MeshingParameters& mparam) const
+  {
+    static Timer t1("Surface Meshing"); RegionTimer regt(t1);
+
+    Array<int, PointIndex> glob2loc(mesh.GetNP());
+    for(auto k : Range(faces))
+      {
+        const auto& face = *faces[k];
+        auto bb = face.GetBoundingBox();
+        bb.Increase(bb.Diam()/10);
+        Meshing2 meshing(*this, mparam, bb);
+        glob2loc = 0;
+        int cntp = 0;
+
+        for(auto& seg : mesh.LineSegments())
+          {
+            if(seg.si == k+1)
+              {
+                for(auto j : Range(2))
+                  {
+                    auto pi = seg[j];
+                    if(glob2loc[pi] == 0)
+                      {
+                        meshing.AddPoint(mesh[pi], pi);
+                        cntp++;
+                        glob2loc[pi] = cntp;
+                      }
+                  }
+              }
+          }
+        for(auto & seg : mesh.LineSegments())
+          {
+            if(seg.si == k+1)
+              {
+                PointGeomInfo gi0, gi1;
+                gi0.trignum = gi1.trignum = k+1;
+                gi0.u = seg.epgeominfo[0].u;
+                gi0.v = seg.epgeominfo[0].v;
+                gi1.u = seg.epgeominfo[1].u;
+                gi1.v = seg.epgeominfo[1].v;
+                meshing.AddBoundaryElement(glob2loc[seg[0]],
+                                           glob2loc[seg[1]],
+                                           gi0, gi1);
+              }
+          }
+
+        // TODO Set max area 2* area of face
+
+        auto noldsurfels = mesh.GetNSE();
+
+
+        static Timer t("GenerateMesh"); RegionTimer reg(t);
+        MESHING2_RESULT res = meshing.GenerateMesh(mesh, mparam, mparam.maxh, k+1);
+
+        for(auto i : Range(noldsurfels, mesh.GetNSE()))
+          {
+            mesh.SurfaceElements()[i].SetIndex(k+1);
+          }
+      }
+  }
+
+  void NetgenGeometry :: OptimizeSurface(Mesh& mesh, const MeshingParameters& mparam) const
   {
     const auto savetask = multithread.task;
     multithread.task = "Optimizing surface";

libsrc/meshing/basegeom.hpp

@@ -12,10 +12,30 @@ struct Tcl_Interp;
 
 namespace netgen
 {
+  class GeometryEdge
+  {
+  public:
+    virtual ~GeometryEdge() {}
+  };
+
+  class GeometryFace
+  {
+  public:
+    virtual ~GeometryFace() {}
+    virtual size_t GetNBoundaries() const = 0;
+    virtual Array<GeometryEdge> GetBoundary(size_t index) 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 Box<3> GetBoundingBox() const = 0;
+  };
 
   class DLL_HEADER NetgenGeometry
   {
     unique_ptr<Refinement> ref;
+  protected:
+    Array<unique_ptr<GeometryFace>> faces;
+    Box<3> bounding_box;
   public:
     NetgenGeometry()
     {
@@ -35,10 +55,12 @@ namespace netgen
 
     virtual Mesh::GEOM_TYPE GetGeomType() const { return Mesh::NO_GEOM; }
     virtual void Analyse(Mesh& mesh,
-                         const MeshingParameters& mparam) {}
+                         const MeshingParameters& mparam) const;
-    virtual void FindEdges(Mesh& mesh, const MeshingParameters& mparam) {}
+    virtual void RestrictLocalMeshsize(Mesh& mesh,
-    virtual void MeshSurface(Mesh& mesh, const MeshingParameters& mparam) {}
+                                       const MeshingParameters& mparam) const {}
-    virtual void OptimizeSurface(Mesh& mesh, const MeshingParameters& mparam);
+    virtual void FindEdges(Mesh& mesh, const MeshingParameters& mparam) const;
+    virtual void MeshSurface(Mesh& mesh, const MeshingParameters& mparam) const;
+    virtual void OptimizeSurface(Mesh& mesh, const MeshingParameters& mparam) const;
 
     virtual void FinalizeMesh(Mesh& mesh) const {}
 

libsrc/meshing/meshing2.cpp

@@ -38,8 +38,10 @@ namespace netgen
   static Array<unique_ptr<netrule>> global_quad_rules;
 
   
-  Meshing2 :: Meshing2 (const MeshingParameters & mp, const Box<3> & aboundingbox)
+  Meshing2 :: Meshing2 (const NetgenGeometry& ageo,
-    : adfront(aboundingbox), boundingbox(aboundingbox)
+                        const MeshingParameters & mp,
+                        const Box<3> & aboundingbox)
+    : geo(ageo), adfront(aboundingbox), boundingbox(aboundingbox)
   {
     static Timer t("Mesing2::Meshing2"); RegionTimer r(t);
 
@@ -133,28 +135,37 @@ namespace netgen
   // static Vec3d ex, ey;
   // static Point3d globp1;
 
-  void Meshing2 :: DefineTransformation (const Point<3> & p1, const Point<3> & p2,
+  void Meshing2 :: DefineTransformation (const Point<3> & ap1,
-					 const PointGeomInfo * geominfo1,
+                                         const Point<3> & ap2,
-					 const PointGeomInfo * geominfo2)
+					 const PointGeomInfo * gi1,
+					 const PointGeomInfo * gi2)
   {
-    globp1 = p1;
+    p1 = ap1;
-    ex = p2 - p1;
+    p2 = ap2;
-    ex /= ex.Length();
+    auto n1 = geo.GetNormal(gi1->trignum, p1, *gi1);
-    ey.X() = -ex.Y();
+    auto n2 = geo.GetNormal(gi2->trignum, p2, *gi2);
-    ey.Y() =  ex.X();
+
-    ey.Z() = 0;
+    ez = 0.5 * (n1+n2);
+    ez.Normalize();
+    ex = (p2-p1).Normalize();
+    ez -= (ez*ex)*ex;
+    ez.Normalize();
+    ey = Cross(ez, ex);
   }
 
   void Meshing2 :: TransformToPlain (const Point<3> & locpoint, 
-				     const MultiPointGeomInfo & geominf,
+				     const MultiPointGeomInfo & geominfo,
 				     Point<2> & plainpoint, double h, int & zone)
   {
-    Vec3d p1p (globp1, locpoint);
+    auto& gi = geominfo.GetPGI(1);
+    auto n = geo.GetNormal(gi.trignum, locpoint, gi);
+    auto p1p = locpoint - p1;
+    plainpoint(0) = (p1p * ex) / h;
+    plainpoint(1) = (p1p * ey) / h;
 
-    //    p1p = locpoint - globp1;
+    if(n*ez < 0)
-    p1p /= h;
+      zone = -1;
-    plainpoint[0] = p1p * ex;
+    else
-    plainpoint[1] = p1p * ey;
       zone = 0;
   }
 
@@ -163,12 +174,9 @@ namespace netgen
 				      PointGeomInfo & gi, 
 				      double h)
   {
-    Vec3d p1p;
+    locpoint = p1 + (h*plainpoint(0)) * ex + (h* plainpoint(1)) * ey;
-    gi.trignum = 1;
+    if (!geo.ProjectPointGI(gi.trignum, locpoint, gi))
-
+      geo.ProjectPoint(gi.trignum, locpoint);
-    p1p = plainpoint[0] * ex + plainpoint[1] * ey;
-    p1p *= h;
-    locpoint = globp1 + p1p;
     return 0;
   }
 

libsrc/meshing/meshing2.hpp

@@ -41,12 +41,16 @@ class Meshing2
   ///
   double maxarea;
 
-  Vec3d ex, ey;
+  Vec3d ex, ey, ez;
-  Point3d globp1;
+  Point<3> p1, p2;
+
+  const NetgenGeometry& geo;
 
 public:
   ///
-  DLL_HEADER Meshing2 (const MeshingParameters & mp, const Box<3> & aboundingbox);
+  DLL_HEADER Meshing2 (const NetgenGeometry& geo,
+                       const MeshingParameters & mp,
+                       const Box<3> & aboundingbox);
 
   ///
   DLL_HEADER virtual ~Meshing2 ();

libsrc/occ/occgenmesh.cpp

@@ -306,7 +306,7 @@ namespace netgen
 
 
 
-  void OCCFindEdges (OCCGeometry & geom, Mesh & mesh, const MeshingParameters & mparam)
+  void OCCFindEdges (const OCCGeometry & geom, Mesh & mesh, const MeshingParameters & mparam)
   {
     static Timer t("OCCFindEdges"); RegionTimer r(t);
     static Timer tsearch("OCCFindEdges - search point");
@@ -601,7 +601,7 @@ namespace netgen
 
 
 
-  void OCCMeshSurface (OCCGeometry & geom, Mesh & mesh,
+  void OCCMeshSurface (const OCCGeometry & geom, Mesh & mesh,
                        const MeshingParameters & mparam)
   {
     static Timer t("OCCMeshSurface"); RegionTimer r(t);
@@ -651,7 +651,7 @@ namespace netgen
         
         static Timer tinit("init");
         tinit.Start();
-        Meshing2OCCSurfaces meshing(TopoDS::Face(geom.fmap(k)), bb, projecttype, mparam);
+        Meshing2OCCSurfaces meshing(geom, TopoDS::Face(geom.fmap(k)), bb, projecttype, mparam);
         tinit.Stop();
 
 

libsrc/occ/occgeom.cpp

@@ -75,19 +75,19 @@ void STEP_GetEntityName(const TopoDS_Shape & theShape, STEPCAFControl_Reader * a
 }
 
   void OCCGeometry :: Analyse(Mesh& mesh,
-                              const MeshingParameters& mparam)
+                              const MeshingParameters& mparam) const
   {
     OCCSetLocalMeshSize(*this, mesh, mparam, occparam);
   }
 
   void OCCGeometry :: FindEdges(Mesh& mesh,
-                                const MeshingParameters& mparam)
+                                const MeshingParameters& mparam) const
   {
     OCCFindEdges(*this, mesh, mparam);
   }
 
   void OCCGeometry :: MeshSurface(Mesh& mesh,
-                                  const MeshingParameters& mparam)
+                                  const MeshingParameters& mparam) const
   {
     OCCMeshSurface(*this, mesh, mparam);
   }

libsrc/occ/occgeom.hpp

@@ -224,7 +224,7 @@ namespace netgen
     Handle_XCAFDoc_ColorTool face_colours;
      
     mutable int changed;
-    NgArray<int> facemeshstatus;
+    mutable NgArray<int> facemeshstatus;
 
     // Philippose - 15/01/2009
     // Maximum mesh size for a given face
@@ -268,11 +268,11 @@ namespace netgen
     { occparam = par; }
 
     void Analyse(Mesh& mesh,
-                 const MeshingParameters& mparam) override;
+                 const MeshingParameters& mparam) const override;
     void FindEdges(Mesh& mesh,
-                   const MeshingParameters& mparam) override;
+                   const MeshingParameters& mparam) const override;
     void MeshSurface(Mesh& mesh,
-                     const MeshingParameters& mparam) override;
+                     const MeshingParameters& mparam) const override;
  
     void FinalizeMesh(Mesh& mesh) const override;
      
@@ -459,11 +459,11 @@ namespace netgen
   DLL_HEADER extern void OCCSetLocalMeshSize(const OCCGeometry & geom, Mesh & mesh, const MeshingParameters & mparam,
                                              const OCCParameters& occparam);
 
-  DLL_HEADER extern void OCCMeshSurface (OCCGeometry & geom, Mesh & mesh, const MeshingParameters & mparam);
+  DLL_HEADER extern void OCCMeshSurface (const OCCGeometry & geom, Mesh & mesh, const MeshingParameters & mparam);
 
   DLL_HEADER extern void OCCOptimizeSurface (OCCGeometry & geom, Mesh & mesh, const MeshingParameters & mparam);
 
-  DLL_HEADER extern void OCCFindEdges (OCCGeometry & geom, Mesh & mesh, const MeshingParameters & mparam);
+  DLL_HEADER extern void OCCFindEdges (const OCCGeometry & geom, Mesh & mesh, const MeshingParameters & mparam);
 }
 
 #endif

libsrc/occ/occmeshsurf.cpp

@@ -492,10 +492,12 @@ namespace netgen
   } 
 
 
-  Meshing2OCCSurfaces :: Meshing2OCCSurfaces (const TopoDS_Shape & asurf,
+  Meshing2OCCSurfaces :: Meshing2OCCSurfaces (const NetgenGeometry& geo,
+                                              const TopoDS_Shape & asurf,
 					      const Box<3> & abb, int aprojecttype,
                                               const MeshingParameters & mparam)
-    : Meshing2(mparam, Box<3>(abb.PMin(), abb.PMax())), surface(TopoDS::Face(asurf), aprojecttype)
+    : Meshing2(geo, mparam, Box<3>(abb.PMin(), abb.PMax())),
+      surface(TopoDS::Face(asurf), aprojecttype)
   {
     ;
   }

libsrc/occ/occmeshsurf.hpp

@@ -113,7 +113,8 @@ class Meshing2OCCSurfaces : public Meshing2
 
 public:
   ///
-  Meshing2OCCSurfaces (const TopoDS_Shape & asurf, const Box<3> & aboundingbox,
+  Meshing2OCCSurfaces (const NetgenGeometry& geo,
+                       const TopoDS_Shape & asurf, const Box<3> & aboundingbox,
                        int aprojecttype, const MeshingParameters & mparam);
 
   ///

libsrc/stlgeom/meshstlsurface.cpp

@@ -897,7 +897,7 @@ void STLSurfaceOptimization (STLGeometry & geom,
 
 MeshingSTLSurface :: MeshingSTLSurface (STLGeometry & ageom,
 					const MeshingParameters & mp)
-  : Meshing2(mp, ageom.GetBoundingBox()), geom(ageom)
+  : Meshing2(ageom, mp, ageom.GetBoundingBox()), geom(ageom)
 {
   ;
 }