diff --git a/libsrc/core/array.hpp b/libsrc/core/array.hpp
index 00859e55..35d1ab30 100644
--- a/libsrc/core/array.hpp
+++ b/libsrc/core/array.hpp
@@ -216,6 +216,10 @@ namespace ngcore
   template <typename  T>
   constexpr T IndexBASE () { return T(0); }
 
+  template <typename  T>
+  constexpr T IndexBASE (T ind) { return IndexBASE<T>(); }
+
+  
 
   class IndexFromEnd
   {
diff --git a/libsrc/interface/writegmsh2.cpp b/libsrc/interface/writegmsh2.cpp
index 845556bd..74843195 100644
--- a/libsrc/interface/writegmsh2.cpp
+++ b/libsrc/interface/writegmsh2.cpp
@@ -58,7 +58,7 @@ namespace netgen
       int np = mesh.GetNP();  /// number of points in mesh
       int ne = mesh.GetNE();  /// number of 3D elements in mesh
       int nse = mesh.GetNSE();  /// number of surface elements (BC)
-      int i, j, k, l;
+      // int i, j, k, l;
 
 
       /*
@@ -66,8 +66,8 @@ namespace netgen
       */
 
       if ((ne > 0)
-         && (mesh.VolumeElement(1).GetNP() <= 10)
-         && (mesh.SurfaceElement(1).GetNP() <= 6))
+          && (mesh.VolumeElements().First().GetNP() <= 10)
+          && (mesh.SurfaceElements().First().GetNP() <= 6))
       {
          cout << "Write GMSH v2.xx Format \n";
          cout << "The GMSH v2.xx export is currently available for elements upto 2nd Order\n" << endl;
@@ -86,7 +86,7 @@ namespace netgen
          outfile << "$Nodes\n";
          outfile << np << "\n";
 
-         for (i = 1; i <= np; i++)
+         for (int i = 1; i <= np; i++)
          {
             const Point3d & p = mesh.Point(i);
             outfile << i << " "; /// node number
@@ -101,13 +101,13 @@ namespace netgen
          outfile << "$Elements\n";
          outfile << ne + nse << "\n";  ////  number of elements + number of surfaces BC
 
-         for (i = 1; i <= nse; i++)
+         for (auto sei : Range(mesh.SurfaceElements()))
          {
             int elType = 0;
 
-            Element2d el = mesh.SurfaceElement(i);
+            Element2d el = mesh[sei]; // .SurfaceElement(i);
             if(invertsurf) el.Invert();
-
+            
             if(el.GetNP() == 3) elType = GMSH_TRIG;	//// GMSH Type for a 3 node triangle
             if(el.GetNP() == 6) elType = GMSH_TRIG6;  //// GMSH Type for a 6 node triangle
             if(elType == 0)
@@ -116,7 +116,7 @@ namespace netgen
                return;
             }
 
-            outfile << i;
+            outfile << sei-IndexBASE(sei)+1;  
             outfile << " ";
             outfile << elType;
             outfile << " ";
@@ -125,7 +125,7 @@ namespace netgen
             outfile << mesh.GetFaceDescriptor (el.GetIndex()).BCProperty() << " ";
             /// that means that physical entity = elementary entity (arbitrary approach)
             outfile << mesh.GetFaceDescriptor (el.GetIndex()).BCProperty() << " ";
-            for (j = 1; j <= el.GetNP(); j++)
+            for (int j = 1; j <= el.GetNP(); j++)
             {
                outfile << " ";
                outfile << el.PNum(triGmsh[j]);
@@ -133,12 +133,12 @@ namespace netgen
             outfile << "\n";
          }
 
-
-         for (i = 1; i <= ne; i++)
+         for (ElementIndex ei : Range(mesh.VolumeElements()))
          {
+           int i = ei-IndexBASE(ei)+1;
             int elType = 0;
 
-            Element el = mesh.VolumeElement(i);
+            Element el = mesh[ei];
             if (inverttets) el.Invert();
 
             if(el.GetNP() == 4) elType = GMSH_TET;    //// GMSH Element type for 4 node tetrahedron
@@ -160,7 +160,7 @@ namespace netgen
             outfile << " ";
             outfile << 100000 + el.GetIndex();   /// volume number
             outfile << " ";
-            for (j = 1; j <= el.GetNP(); j++)
+            for (int j = 1; j <= el.GetNP(); j++)
             {
                outfile << " ";
                outfile << el.PNum(tetGmsh[j]);
@@ -193,7 +193,7 @@ namespace netgen
          outfile << "$Nodes\n";
          outfile << np << "\n";
 
-         for (i = 1; i <= np; i++)
+         for (int i = 1; i <= np; i++)
          {
             const Point3d & p = mesh.Point(i);
             outfile << i << " "; /// node number
@@ -207,7 +207,7 @@ namespace netgen
          outfile << "$Elements\n";
          outfile << nse << "\n";
 
-         for (k = 1; k <= nse; k++)
+         for (int k = 1; k <= nse; k++)
          {
             int elType = 0;
 
@@ -232,7 +232,7 @@ namespace netgen
             outfile << mesh.GetFaceDescriptor (el.GetIndex()).BCProperty() << " ";
             /// that means that physical entity = elementary entity (arbitrary approach)
             outfile << mesh.GetFaceDescriptor (el.GetIndex()).BCProperty() << " ";
-            for (l = 1; l <= el.GetNP(); l++)
+            for (int l = 1; l <= el.GetNP(); l++)
             {
                outfile << " ";
                if((elType == GMSH_TRIG) || (elType == GMSH_TRIG6))
diff --git a/libsrc/meshing/bisect.cpp b/libsrc/meshing/bisect.cpp
index 8043cb87..2870bb1d 100644
--- a/libsrc/meshing/bisect.cpp
+++ b/libsrc/meshing/bisect.cpp
@@ -5,10 +5,8 @@
 
 #include "meshing.hpp"  // quickfix for parallel 
 
-
 #define noDEBUG
 
-
 namespace netgen
 {
   class MarkedTet
diff --git a/libsrc/meshing/clusters.cpp b/libsrc/meshing/clusters.cpp
index 2cba410c..c71bc733 100644
--- a/libsrc/meshing/clusters.cpp
+++ b/libsrc/meshing/clusters.cpp
@@ -143,19 +143,21 @@ namespace netgen
 	  cluster_reps.Elem(nnums[j]) = nnums[j];
       }
     */
+
+    
     ngcore::ParallelForRange
       (mesh.SurfaceElements().Range(),
        [&] (auto myrange)
        {
          NgArrayMem<int,9> nnums; // , ednums;
-         for (int i_ : myrange)
+         for (SurfaceElementIndex i_ : myrange)
            {
              int i = i_+1;
-             const Element2d & el = mesh.SurfaceElement(i);
+             const Element2d & el = mesh[i_]; // .SurfaceElement(i);
              ELEMENT_TYPE typ = el.GetType();
              
              // top.GetSurfaceElementEdges (i, ednums);
-             auto ednums = top.GetEdges (SurfaceElementIndex(i_));
+             auto ednums = top.GetEdges (i_);
              // cout << "ednums = " << ednums << endl;
              
              int fanum = top.GetSurfaceElementFace (i);             
diff --git a/libsrc/meshing/improve3.cpp b/libsrc/meshing/improve3.cpp
index a0811112..a021f8bc 100644
--- a/libsrc/meshing/improve3.cpp
+++ b/libsrc/meshing/improve3.cpp
@@ -357,7 +357,7 @@ void MeshOptimize3d :: CombineImprove ()
 {
   static Timer t("MeshOptimize3d::CombineImprove"); RegionTimer reg(t);
   static Timer topt("Optimize");
-  static Timer tsearch("Search");
+  static Timer tsearch("Search-combine");
   static Timer tbuild_elements_table("Build elements table");
 
   mesh.BuildBoundaryEdges(false);
@@ -613,7 +613,7 @@ void MeshOptimize3d :: SplitImprove ()
 {
   static Timer t("MeshOptimize3d::SplitImprove"); RegionTimer reg(t);
   static Timer topt("Optimize");
-  static Timer tsearch("Search");
+  static Timer tsearch("Search-split");
 
   // int np = mesh.GetNP();
   int ne = mesh.GetNE();
diff --git a/libsrc/meshing/meshclass.cpp b/libsrc/meshing/meshclass.cpp
index be26bf80..16ee4a0a 100644
--- a/libsrc/meshing/meshclass.cpp
+++ b/libsrc/meshing/meshclass.cpp
@@ -3051,6 +3051,7 @@ namespace netgen
     // bool buggy = false;
     // ofstream bout("buggy.out");
 
+
     for (int i = 1; i <= GetNSE(); i++)
       {
         const Element2d & el = SurfaceElement(i);
@@ -3060,10 +3061,10 @@ namespace netgen
           {
             for (int j = 1; j <= el.GetNP(); j++)
               {
-                INDEX_3 seg (el.PNumMod(j), el.PNumMod(j+1), el.GetIndex());
+                PointIndices<3> seg (el.PNumMod(j), el.PNumMod(j+1), el.GetIndex());
                 // int data;
 
-                if (seg.I1() < PointIndex::BASE || seg.I2() < PointIndex::BASE)
+                if (!seg.I1().IsValid() || !seg.I2().IsValid())
                   cerr << "seg = " << seg << endl;
 
                 if (faceht.Used(seg))
@@ -7299,6 +7300,7 @@ namespace netgen
       {
         int eli0, eli1;
         GetTopology().GetSurface2VolumeElement(sei+1, eli0, eli1);
+        // auto [ei0,ei1] = GetTopology().GetSurface2VolumeElement(sei); // the way to go
         auto & sel = (*this)[sei];
         int face = sel.GetIndex();
         int domin = VolumeElement(eli0).GetIndex();
diff --git a/libsrc/meshing/meshing3.cpp b/libsrc/meshing/meshing3.cpp
index 57e19d49..0ed40bb1 100644
--- a/libsrc/meshing/meshing3.cpp
+++ b/libsrc/meshing/meshing3.cpp
@@ -9,8 +9,6 @@ double minwithoutother;
 
 
 
-
-
 MeshingStat3d :: MeshingStat3d ()
 {
   cntsucc = cnttrials = cntelem = qualclass = 0;
diff --git a/libsrc/meshing/meshtype.hpp b/libsrc/meshing/meshtype.hpp
index 843ba9ca..1fa2ed47 100644
--- a/libsrc/meshing/meshtype.hpp
+++ b/libsrc/meshing/meshtype.hpp
@@ -27,7 +27,6 @@ namespace netgen
   */
 
 
-
   enum ELEMENT_TYPE : unsigned char { 
     SEGMENT = 1, SEGMENT3 = 2,
     TRIG = 10, QUAD=11, TRIG6 = 12, QUAD6 = 13, QUAD8 = 14,
@@ -149,8 +148,107 @@ namespace netgen
   }
 
 
+  template <typename T, typename TIndex, int Base>
+  class Index
+  {
+  public:
+
+    T i;
+    static constexpr int BASE = Base;
+    
+  public:
+    class t_invalid { public: constexpr t_invalid() = default; };
+    static constexpr t_invalid INVALID{};
+    
+    constexpr Index () = default;
+    constexpr Index (const Index&) = default;
+    constexpr Index (Index &&) = default;
+    Index & operator= (const Index&) = default;
+    Index & operator= (Index&&) = default;
+
+    // private:
+    constexpr Index (int ai) : i(ai)
+    {
+#ifdef DEBUG
+      if (ai < Base)
+        cout << "illegal Index, use Index::INVALID instead" << endl;
+#endif
+    }
 
 
+    // friend constexpr netgen::TIndex ngcore::IndexBASE<netgen::TIndex> ();
+    // template <int N> friend class PointIndices;
+
+    /*
+    friend auto operator+ (Index, int) -> TIndex;
+    friend TIndex operator+ (Index, size_t);    
+    friend TIndex operator+ (int, Index);
+    friend TIndex operator+ (size_t, Index);
+    friend constexpr TIndex operator- (Index, int);
+    friend int operator- (Index, Index);
+    friend bool operator< (Index a, Index b);
+    friend bool operator> (Index a, Index b);
+    friend bool operator>= (Index a, Index b);
+    friend bool operator<= (Index a, Index b);
+    friend bool operator== (Index a, Index b);
+    friend bool operator!= (Index a, Index b);
+    */
+    
+  public:
+    constexpr Index (t_invalid inv) : i(long(BASE)-1) { ; }
+    // TIndex & operator= (const TIndex &ai) { i = ai.i; return *this; }
+    // private:
+    constexpr operator const int& () const { return i; }
+    explicit constexpr operator int& () { return i; }
+  public:
+    constexpr operator TIndex() const { return TIndex(i); }
+    constexpr operator TIndex&() { return static_cast<TIndex&>(*this); }    
+    TIndex operator++ (int) { TIndex hi(*this); i++; return hi; }
+    TIndex operator-- (int) { TIndex hi(*this); i--; return hi; }
+    TIndex & operator++ () { i++; return *this; }
+    TIndex operator-- () { i--; return *this; }
+    TIndex operator+= (int add) { i += add; return *this; }
+    void Invalidate() { i = long(TIndex::BASE)-1; }
+    bool IsValid() const { return i+1 != TIndex::BASE; }
+    // operator bool() const { return IsValid(); }
+
+    void DoArchive (Archive & ar) { ar & i; }
+  };
+
+
+  /*
+  template <typename T, typename TIndex, int Base>  
+  auto operator+ (Index<T,TIndex,Base> pi, int i) -> TIndex { return TIndex(pi.i+i); }
+  template <typename T, typename TIndex, int Base>    
+  inline TIndex operator+ (Index<T,TIndex,Base> pi, size_t i) { return TIndex(pi.i+i); }
+  template <typename T, typename TIndex, int Base>    
+  inline TIndex operator+ (int i, Index<T,TIndex,Base> pi) { return TIndex(pi.i+i); }
+  template <typename T, typename TIndex, int Base>    
+  inline TIndex operator+ (size_t i, Index<T,TIndex,Base> pi) { return TIndex(pi.i+i); }
+  
+  template <typename T, typename TIndex, int Base>    
+  constexpr inline auto operator- (Index<T,TIndex,Base> pi, int i) -> TIndex { return TIndex(pi.i-i); }
+  
+  template <typename T, typename TIndex, int Base>    
+  inline int operator- (Index<T,TIndex,Base> pa, Index<T,TIndex,Base> pb) { return pa.i-pb.i; }
+  template <typename T, typename TIndex, int Base>      
+  inline bool operator< (Index<T,TIndex,Base> a, Index<T,TIndex,Base> b) { return a.i-b.i < 0; }
+  template <typename T, typename TIndex, int Base>      
+  inline bool operator> (Index<T,TIndex,Base> a, Index<T,TIndex,Base> b) { return a.i-b.i > 0; }
+  template <typename T, typename TIndex, int Base>      
+  inline bool operator>= (Index<T,TIndex,Base> a, Index<T,TIndex,Base> b) { return a.i-b.i >= 0; }
+  template <typename T, typename TIndex, int Base>      
+  inline bool operator<= (Index<T,TIndex,Base> a, Index<T,TIndex,Base> b) { return a.i-b.i <= 0; }
+  template <typename T, typename TIndex, int Base>      
+  inline bool operator== (Index<T,TIndex,Base> a, Index<T,TIndex,Base> b) { return a.i == b.i; }
+  template <typename T, typename TIndex, int Base>      
+  inline bool operator!= (Index<T,TIndex,Base> a, Index<T,TIndex,Base> b) { return a.i != b.i; }
+  */
+  
+
+  
+
+    
   /*
   class PointIndex
   {
@@ -196,7 +294,9 @@ namespace netgen
 
 
   // #define BASE0
-  
+
+
+  /*
   class PointIndex
   {
     int i;
@@ -260,6 +360,14 @@ namespace netgen
     void DoArchive (Archive & ar) { ar & i; }
   };
 
+  */
+
+  class PointIndex : public Index<int,PointIndex,1>
+  {
+  public:
+    using Index::Index;
+  };
+    
   constexpr inline PointIndex operator+ (PointIndex pi, int i) { return PointIndex(pi.i+i); }
   constexpr inline PointIndex operator+ (PointIndex pi, size_t i) { return PointIndex(pi.i+i); }  
   constexpr inline PointIndex operator+ (int i, PointIndex pi) { return PointIndex(pi.i+i); }
@@ -501,8 +609,8 @@ namespace netgen
   inline bool operator> (ElementIndex ei1, ElementIndex ei2) { return int(ei1) > int(ei2); };  
   inline bool operator>= (ElementIndex ei1, ElementIndex ei2) { return int(ei1) >= int(ei2); };
   inline bool operator<= (ElementIndex ei1, ElementIndex ei2) { return int(ei1) <= int(ei2); };
-  // these should not be needed:
 
+  // these should not be needed:
   inline bool operator== (ElementIndex ei1, int ei2) { return int(ei1) == int(ei2); };  
   inline bool operator< (size_t s, ElementIndex ei2) { return int(s) < int(ei2); };    
   inline bool operator< (ElementIndex ei1, size_t s) { return int(ei1) < int(s); };   // should not need
@@ -530,6 +638,27 @@ namespace netgen
     SurfaceElementIndex & operator+= (int inc) { i+=inc; return *this; }
     void DoArchive (Archive & ar) { ar & i; }
   };
+
+  inline SurfaceElementIndex operator+ (SurfaceElementIndex ei, int i) { return SurfaceElementIndex { int(ei) + i }; }
+  inline SurfaceElementIndex operator+ (size_t s, SurfaceElementIndex ei) { return SurfaceElementIndex(int(ei) + s); }  
+  inline SurfaceElementIndex operator+ (SurfaceElementIndex ei, size_t s) { return SurfaceElementIndex(int(ei) + s); }  
+  inline bool operator== (SurfaceElementIndex ei1, SurfaceElementIndex ei2) { return int(ei1) == int(ei2); };
+  inline bool operator!= (SurfaceElementIndex ei1, SurfaceElementIndex ei2) { return int(ei1) != int(ei2); };
+  inline bool operator< (SurfaceElementIndex ei1, SurfaceElementIndex ei2) { return int(ei1) < int(ei2); };  
+  inline bool operator> (SurfaceElementIndex ei1, SurfaceElementIndex ei2) { return int(ei1) > int(ei2); };  
+  inline bool operator>= (SurfaceElementIndex ei1, SurfaceElementIndex ei2) { return int(ei1) >= int(ei2); };
+  inline bool operator<= (SurfaceElementIndex ei1, SurfaceElementIndex ei2) { return int(ei1) <= int(ei2); };
+
+  // these should not be needed:
+  inline bool operator== (SurfaceElementIndex ei1, int ei2) { return int(ei1) == int(ei2); };
+  inline bool operator== (int ei2, SurfaceElementIndex ei1) { return int(ei1) == int(ei2); };
+  inline bool operator!= (SurfaceElementIndex ei1, int ei2) { return int(ei1) != int(ei2); };    
+  inline bool operator< (size_t s, SurfaceElementIndex ei2) { return int(s) < int(ei2); };    
+  inline bool operator< (SurfaceElementIndex ei1, size_t s) { return int(ei1) < int(s); };   // should not need
+  inline bool operator< (SurfaceElementIndex ei1, int s) { return int(ei1) < int(s); };   // should not need
+  inline bool operator>= (size_t s, SurfaceElementIndex ei2) { return int(s) >= int(ei2); };
+  inline bool operator>= (SurfaceElementIndex ei1, int s) { return int(ei1) >= int(s); };
+
   
   inline void SetInvalid (SurfaceElementIndex & id) { id = -1; }
   inline bool IsInvalid (SurfaceElementIndex & id) { return id == -1; }
diff --git a/libsrc/meshing/parallelmesh.cpp b/libsrc/meshing/parallelmesh.cpp
index 4d7b9a1d..7745f33a 100644
--- a/libsrc/meshing/parallelmesh.cpp
+++ b/libsrc/meshing/parallelmesh.cpp
@@ -592,18 +592,17 @@ namespace netgen
       }
 
     DynamicTable<int> elementarrays(elarraysize);
-
-    for (int ei = 1; ei <= GetNE(); ei++)
+    
+    for (ElementIndex ei : VolumeElements().Range())    
       {
 	const Element & el = VolumeElement (ei);
-	// int dest = el.GetPartition();
-        int dest = vol_partition[ei-1];
+        int dest = vol_partition[ei];
         
-	elementarrays.Add (dest, ei);
+	elementarrays.Add (dest, int(ei+1));
 	elementarrays.Add (dest, el.GetIndex());
 	elementarrays.Add (dest, el.GetNP());
-	for (int i = 0; i < el.GetNP(); i++)
-	  elementarrays.Add (dest, el[i]);
+        for (PointIndex pi : el.PNums())
+	  elementarrays.Add (dest, pi);
       }
     tbuildelementtable.Stop();
     
diff --git a/libsrc/meshing/refine.cpp b/libsrc/meshing/refine.cpp
index 9c5f5556..dbbfec10 100644
--- a/libsrc/meshing/refine.cpp
+++ b/libsrc/meshing/refine.cpp
@@ -865,9 +865,10 @@ namespace netgen
 		for (int k = 1; k <= 3; k++)
 		  {
 		    fhelp.Clear();
-		    for (int i = 1; i <= mesh.GetNE(); i++)
+		    // for (int i = 1; i <= mesh.GetNE(); i++) 
+                    for (const Element & el : mesh.VolumeElements())
 		      {
-			const Element & el = mesh.VolumeElement(i);
+			// const Element & el = mesh.VolumeElement(i);
 			int freeel = 0;
 			for (int j = 1; j <= el.GetNP(); j++)
 			  if (free.Test(el.PNum(j)))
@@ -897,19 +898,19 @@ namespace netgen
 
 
 		wrongels = 0;
-		for (int i = 1; i <= mesh.GetNE(); i++)
+                for (ElementIndex ei : mesh.VolumeElements().Range())
 		  {
-		    if (mesh.VolumeElement(i).Volume(mesh.Points()) < 0)
+		    if (mesh.VolumeElement(ei).Volume(mesh.Points()) < 0)
 		      {
 			wrongels++;
-			mesh.VolumeElement(i).Flags().badel = 1;
+			mesh.VolumeElement(ei).Flags().badel = 1;
 			(*testout) << "wrong el: ";
 			for (int j = 1; j <= 4; j++)
-			  (*testout) << mesh.VolumeElement(i).PNum(j) << " ";
+			  (*testout) << mesh.VolumeElement(ei).PNum(j) << " ";
 			(*testout) << endl;
 		      }
 		    else
-		      mesh.VolumeElement(i).Flags().badel = 0;
+		      mesh.VolumeElement(ei).Flags().badel = 0;
 		  }
 		cout << "wrongels = " << wrongels << endl;
 	      }
diff --git a/libsrc/meshing/topology.hpp b/libsrc/meshing/topology.hpp
index fcb1b967..782d1ba6 100644
--- a/libsrc/meshing/topology.hpp
+++ b/libsrc/meshing/topology.hpp
@@ -185,6 +185,13 @@ public:
     elnr2 = surf2volelement.Get(selnr)[1];
   }
 
+  std::array<ElementIndex,2> GetSurface2VolumeElement (SurfaceElementIndex sei) 
+  {
+    return { ElementIndex( surf2volelement.Get(sei+1)[0] - 1),
+             ElementIndex( surf2volelement.Get(sei+1)[1] - 1) };
+  }
+
+  
   int GetFace2SurfaceElement (int fnr) const { return face2surfel[fnr-1]; }
 
   SegmentIndex GetSegmentOfEdge(int edgenr) const { return edge2segment[edgenr-1]; }
diff --git a/libsrc/occ/occmeshsurf.cpp b/libsrc/occ/occmeshsurf.cpp
index 37f23cfc..16a0e427 100644
--- a/libsrc/occ/occmeshsurf.cpp
+++ b/libsrc/occ/occmeshsurf.cpp
@@ -366,8 +366,11 @@ namespace netgen
 
   void OCCSurface :: Project (Point<3> & ap, PointGeomInfo & gi)
   {
-    // static Timer t("OccSurface::Project"); RegionTimer reg(t);
-    // static Timer t2("OccSurface::Project actual"); 
+    static Timer t("OccSurface::Project"); RegionTimer reg(t);
+    static Timer tanal("OccSurface::Project analysis"); 
+    static Timer ttol("OccSurface::Project approximation"); 
+
+    static Timer t2("OccSurface::Project actual"); 
 
 
     // try Newton's method ...
@@ -472,14 +475,18 @@ namespace netgen
     */
 
     // double u,v;
-    // JS : shouldn't we move these 2 lines to the constructor ? 
+    // JS : shouldn't we move these 2 lines to the constructor ?
+    // tanal.Start();
     Handle( ShapeAnalysis_Surface ) su = new ShapeAnalysis_Surface( occface );
+    // ShapeAnalysis_Surface su( occface );    
+    // tanal.Stop();
+    ttol.Start();
     auto toltool =  BRep_Tool::Tolerance( topods_face );
-
+    ttol.Stop();
     // gp_Pnt2d suval = su->ValueOfUV ( pnt, toltool);
-    // t2.Start();
+    t2.Start();
     gp_Pnt2d suval = su->NextValueOfUV (gp_Pnt2d(u,v), pnt, toltool);
-    // t2.Stop();
+    t2.Stop();
     suval.Coord( u, v);
     pnt = occface->Value( u, v );
     
diff --git a/libsrc/occ/occmeshsurf.hpp b/libsrc/occ/occmeshsurf.hpp
index b4c4b06e..c045da69 100644
--- a/libsrc/occ/occmeshsurf.hpp
+++ b/libsrc/occ/occmeshsurf.hpp
@@ -5,10 +5,13 @@
 
 #include "occgeom.hpp"
 #include "mydefs.hpp"
-
+#include <BRep_Tool.hxx>
 #include <TopoDS_Face.hxx>
 #include <Geom_Surface.hxx>
 #include <ShapeAnalysis.hxx>
+#include <ShapeAnalysis_Surface.hxx>
+#include <GeomLProp_SLProps.hxx>
+
 
 #define PARAMETERSPACE -1
 #define PLANESPACE     1
@@ -30,7 +33,8 @@ public:
   Handle(Geom_Surface) occface;
   TopAbs_Orientation orient;
   int projecttype;
-
+  ShapeAnalysis_Surface su;
+  Standard_Real toltool;
 protected:
   Point<3> p1;
   Point<3> p2;
@@ -60,10 +64,15 @@ protected:
 
 public:
   OCCSurface (const TopoDS_Face & aface, int aprojecttype)
+    : topods_face(aface),
+      occface(BRep_Tool::Surface(topods_face)),
+      su( occface ),
+      toltool(BRep_Tool::Tolerance(topods_face))
+      
   {
     static Timer t("occurface ctor"); RegionTimer r(t);
     topods_face = aface;
-    occface = BRep_Tool::Surface(topods_face);
+    // occface = BRep_Tool::Surface(topods_face);
     orient = topods_face.Orientation();
     projecttype = aprojecttype;
     ShapeAnalysis::GetFaceUVBounds (topods_face, umin, umax, vmin, vmax);
@@ -72,6 +81,8 @@ public:
     umax += fabs(umax-umin)/100.0;
     vmax += fabs(vmax-vmin)/100.0;
     // projecttype = PLANESPACE;
+
+    // su = ShapeAnalysis_Surface( occface );        
     /*
     TopExp_Explorer exp1;
     exp1.Init (topods_face, TopAbs_WIRE);
diff --git a/libsrc/occ/occpkg.cpp b/libsrc/occ/occpkg.cpp
index ecc11af8..57b1b78e 100644
--- a/libsrc/occ/occpkg.cpp
+++ b/libsrc/occ/occpkg.cpp
@@ -817,21 +817,29 @@ namespace netgen
 
      if(strcmp(argv[1], "showall") == 0)
      {
+       /*
         for(int i = 1; i <= mesh->GetNSE(); i++)
         {
-           mesh->SurfaceElement(i).Visible(1);
+          mesh->SurfaceElement(i).Visible(1);
         }
-
-        mesh->SetNextTimeStamp();
+       */
+       for (auto & el : mesh->SurfaceElements())
+         el.Visible(1);
+       
+       mesh->SetNextTimeStamp();
      }
 
      if(strcmp(argv[1], "hideall") == 0)
      {
+       /*
         for(int i = 1; i <= mesh->GetNSE(); i++)
         {
-           mesh->SurfaceElement(i).Visible(0);
+          mesh->SurfaceElement(i).Visible(0);
         }
-
+       */
+       for (auto & el : mesh->SurfaceElements())
+         el.Visible(0);
+       
         mesh->SetNextTimeStamp();
      }
 
diff --git a/nglib/nglib.cpp b/nglib/nglib.cpp
index 11a43ddb..a73fb324 100644
--- a/nglib/nglib.cpp
+++ b/nglib/nglib.cpp
@@ -270,7 +270,7 @@ namespace nglib
    NGLIB_API Ng_Surface_Element_Type 
       Ng_GetSurfaceElement (Ng_Mesh * mesh, int num, int * pi)
    {
-      const Element2d & el = ((Mesh*)mesh)->SurfaceElement(num);
+     const Element2d & el = ((Mesh*)mesh)->SurfaceElement(SurfaceElementIndex(num-1));
       for (int i = 1; i <= el.GetNP(); i++)
          pi[i-1] = el.PNum(i);
       Ng_Surface_Element_Type et;
@@ -301,7 +301,7 @@ namespace nglib
    NGLIB_API Ng_Volume_Element_Type
       Ng_GetVolumeElement (Ng_Mesh * mesh, int num, int * pi)
    {
-      const Element & el = ((Mesh*)mesh)->VolumeElement(num);
+     const Element & el = ((Mesh*)mesh)->VolumeElement(ElementIndex(num-1));
       for (int i = 1; i <= el.GetNP(); i++)
          pi[i-1] = el.PNum(i);
       Ng_Volume_Element_Type et;
@@ -439,7 +439,7 @@ namespace nglib
    NGLIB_API Ng_Surface_Element_Type
       Ng_GetElement_2D (Ng_Mesh * mesh, int num, int * pi, int * matnum)
    {
-      const Element2d & el = ((Mesh*)mesh)->SurfaceElement(num);
+     const Element2d & el = ((Mesh*)mesh)->SurfaceElement(SurfaceElementIndex(num-1));
       for (int i = 1; i <= el.GetNP(); i++)
          pi[i-1] = el.PNum(i);