diff --git a/libsrc/meshing/basegeom.cpp b/libsrc/meshing/basegeom.cpp
index 429a5d58..d7489a81 100644
--- a/libsrc/meshing/basegeom.cpp
+++ b/libsrc/meshing/basegeom.cpp
@@ -11,7 +11,7 @@ namespace netgen
   { ; }
 
   void NetgenGeometry :: Analyse(Mesh& mesh,
-                                 const MeshingParameters& mparam)
+                                 const MeshingParameters& mparam) const
   {
     static Timer t1("SetLocalMeshsize"); RegionTimer regt(t1);
     mesh.SetGlobalH(mparam.maxh);
@@ -25,7 +25,74 @@ namespace netgen
     mesh.LoadLocalMeshSize(mparam.meshsizefilename);
   }
 
-  void NetgenGeometry :: OptimizeSurface(Mesh& mesh, const MeshingParameters& mparam)
+  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";
diff --git a/libsrc/meshing/basegeom.hpp b/libsrc/meshing/basegeom.hpp
index 86d19546..4efe2cbb 100644
--- a/libsrc/meshing/basegeom.hpp
+++ b/libsrc/meshing/basegeom.hpp
@@ -12,11 +12,29 @@ 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()
@@ -37,12 +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 RestrictLocalMeshsize(Mesh& mesh,
                                        const MeshingParameters& mparam) const {}
-    virtual void FindEdges(Mesh& mesh, const MeshingParameters& mparam) {}
-    virtual void MeshSurface(Mesh& mesh, const MeshingParameters& mparam) {}
-    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 {}
 
diff --git a/libsrc/occ/occgenmesh.cpp b/libsrc/occ/occgenmesh.cpp
index bc63203b..c6732a9d 100644
--- a/libsrc/occ/occgenmesh.cpp
+++ b/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);
diff --git a/libsrc/occ/occgeom.cpp b/libsrc/occ/occgeom.cpp
index bb95ee88..96434678 100644
--- a/libsrc/occ/occgeom.cpp
+++ b/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);
   }
diff --git a/libsrc/occ/occgeom.hpp b/libsrc/occ/occgeom.hpp
index 25369ac7..09561b29 100644
--- a/libsrc/occ/occgeom.hpp
+++ b/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