diff --git a/libsrc/meshing/meshclass.hpp b/libsrc/meshing/meshclass.hpp
index 287c8c0b..68d25075 100644
--- a/libsrc/meshing/meshclass.hpp
+++ b/libsrc/meshing/meshclass.hpp
@@ -388,9 +388,9 @@ namespace netgen
     auto GetNE () const { return volelements.Size(); }
 
     // [[deprecated("Use VolumeElement(ElementIndex) instead of int !")]]    
-    Element & VolumeElement(int i) { return volelements[i-1]; }
+    Element & VolumeElement(int i) { return volelements[IndexBASE<ElementIndex>()+(i-1)]; }
     // [[deprecated("Use VolumeElement(ElementIndex) instead of int !")]]        
-    const Element & VolumeElement(int i) const { return volelements[i-1]; }
+    const Element & VolumeElement(int i) const { return volelements[IndexBASE<ElementIndex>()+(i-1)]; }
     // [[deprecated("Use mesh[](VolumeElementIndex) instead !")]]
     Element & VolumeElement(ElementIndex i) { return volelements[i]; }
     // [[deprecated("Use mesh[](VolumeElementIndex) instead !")]]
diff --git a/libsrc/meshing/meshtype.hpp b/libsrc/meshing/meshtype.hpp
index 335315ef..e8c68e36 100644
--- a/libsrc/meshing/meshtype.hpp
+++ b/libsrc/meshing/meshtype.hpp
@@ -473,11 +473,16 @@ namespace netgen
   public:
     ElementIndex & operator= (const ElementIndex & ai) { i = ai.i; return *this; }
     ElementIndex & operator= (int ai) { i = ai; return *this; }
+    // private:
     constexpr operator int () const { return i; }
+  public:
     ElementIndex operator++ (int) { return ElementIndex(i++); }    
     ElementIndex operator-- (int) { return ElementIndex(i--); }
     ElementIndex & operator++ () { ++i; return *this; }
     ElementIndex & operator-- () { --i; return *this; }
+
+    int operator- (ElementIndex ei2) const { return i-ei2.i; }
+    friend constexpr ElementIndex ngcore::IndexBASE<ElementIndex>();
   };
 
   inline istream & operator>> (istream & ist, ElementIndex & pi)
@@ -485,9 +490,9 @@ namespace netgen
     int i; ist >> i; pi = i; return ist;
   }
 
-  inline ostream & operator<< (ostream & ost, const ElementIndex & pi)
+  inline ostream & operator<< (ostream & ost, const ElementIndex & ei)
   {
-    return (ost << int(pi));
+    return (ost << ei-IndexBASE<ElementIndex>());
   }
 
 
diff --git a/libsrc/meshing/topology.cpp b/libsrc/meshing/topology.cpp
index eb12900c..ce4f9ce7 100644
--- a/libsrc/meshing/topology.cpp
+++ b/libsrc/meshing/topology.cpp
@@ -1807,29 +1807,36 @@ namespace netgen
   void MeshTopology :: GetElementEdges (int elnr, NgArray<int> & eledges) const
   {
     int ned = GetNEdges (mesh->VolumeElement(elnr).GetType());
+    ElementIndex ei = IndexBASE<ElementIndex>() +(elnr-1);
     eledges.SetSize (ned);
     for (int i = 0; i < ned; i++)
-      eledges[i] = edges.Get(elnr)[i]+1;
-      // eledges[i] = abs (edges.Get(elnr)[i]);
+      // eledges[i] = edges.Get(elnr)[i]+1;
+      eledges[i] = edges[ei][i]+1;
   }
 
   void MeshTopology :: GetElementFaces (int elnr, NgArray<int> & elfaces) const
   {
     int nfa = GetNFaces (mesh->VolumeElement(elnr).GetType());
+    ElementIndex ei = IndexBASE<ElementIndex>() +(elnr-1);
+    
     elfaces.SetSize (nfa);
 
     for (auto i : Range(nfa))
-      elfaces[i] = faces.Get(elnr)[i]+1;
+      // elfaces[i] = faces.Get(elnr)[i]+1;
+      elfaces[i] = faces[ei][i]+1;
   }
 
   
   void MeshTopology :: GetElementFaces (int elnr, NgArray<int> & elfaces, bool withorientation) const
   {
     int nfa = GetNFaces (mesh->VolumeElement(elnr).GetType());
+    ElementIndex ei = IndexBASE<ElementIndex>() +(elnr-1);
+    
     elfaces.SetSize (nfa);
 
     for (auto i : Range(nfa))
-      elfaces[i] = faces.Get(elnr)[i]+1;
+      // elfaces[i] = faces.Get(elnr)[i]+1;
+      elfaces[i] = faces[ei][i]+1;
     
     if(withorientation)
     {
@@ -1875,7 +1882,7 @@ namespace netgen
   int MeshTopology :: GetElementEdges (int elnr, int * eledges, int * orient) const
   {
     //  int ned = GetNEdges (mesh.VolumeElement(elnr).GetType());
-    
+    ElementIndex ei = IndexBASE<ElementIndex>() +(elnr-1);
     if (mesh->GetDimension()==3 || 1)
       {
         if (orient)
@@ -1887,8 +1894,12 @@ namespace netgen
 		eledges[i] = abs (edges.Get(elnr)[i]);
 		orient[i] = (edges.Get(elnr)[i] > 0 ) ? 1 : -1;
                 */
-                if (edges.Get(elnr)[i] == -1) return i;
-                eledges[i] = edges.Get(elnr)[i]+1;
+                
+                // if (edges.Get(elnr)[i] == -1) return i;
+                // eledges[i] = edges[ei].Get(elnr)[i]+1;
+                if (edges[ei][i] == -1) return i;
+                eledges[i] = edges[ei][i]+1;
+                
 		// orient[i] = edges.Get(elnr)[i].orient ? -1 : 1;
                 orient[i] = GetElementEdgeOrientation(elnr, i) ? -1 : 1;
 	      }
@@ -1899,8 +1910,11 @@ namespace netgen
 	      {
 		// if (!edges.Get(elnr)[i]) return i;
 		// eledges[i] = abs (edges.Get(elnr)[i]);
-                if (edges.Get(elnr)[i] == -1) return i;
-                eledges[i] = edges.Get(elnr)[i]+1;
+                
+                // if (edges.Get(elnr)[i] == -1) return i;
+                //eledges[i] = edges.Get(elnr)[i]+1;
+                if (edges[ei][i] == -1) return i;
+                eledges[i] = edges[ei][i]+1;
 
 	      }
 	  }
@@ -1933,6 +1947,8 @@ namespace netgen
 
   int MeshTopology :: GetElementFaces (int elnr, int * elfaces, int * orient) const
   {
+    ElementIndex ei = IndexBASE<ElementIndex>() +(elnr-1);
+    
     //  int nfa = GetNFaces (mesh.VolumeElement(elnr).GetType());
     if (orient)
       {
@@ -1943,8 +1959,10 @@ namespace netgen
 	    elfaces[i] = (faces.Get(elnr)[i]-1) / 8 + 1;
 	    orient[i] = (faces.Get(elnr)[i]-1) % 8;
             */
-	    if (faces.Get(elnr)[i] == -1) return i;
-	    elfaces[i] = faces.Get(elnr)[i]+1;
+	    // if (faces.Get(elnr)[i] == -1) return i;
+	    // elfaces[i] = faces.Get(elnr)[i]+1;
+	    if (faces[ei][i] == -1) return i;
+	    elfaces[i] = faces[ei][i]+1;
 	    // orient[i] = faces.Get(elnr)[i].forient;
             orient[i] = GetElementFaceOrientation (elnr, i);
 	  }
@@ -1955,8 +1973,11 @@ namespace netgen
 	  {
 	    // if (!faces.Get(elnr)[i]) return i;
 	    // elfaces[i] = (faces.Get(elnr)[i]-1) / 8 + 1;
-	    if (faces.Get(elnr)[i] == -1) return i;
-	    elfaces[i] = faces.Get(elnr)[i]+1;
+
+	    // if (faces.Get(elnr)[i] == -1) return i;
+	    // elfaces[i] = faces.Get(elnr)[i]+1;
+	    if (faces[ei][i] == -1) return i;
+	    elfaces[i] = faces[ei][i]+1;
 	  }
       }
     return 6;
@@ -1966,9 +1987,11 @@ namespace netgen
   void MeshTopology :: GetSurfaceElementEdges (int elnr, NgArray<int> & eledges) const
   {
     int ned = GetNEdges (mesh->SurfaceElement(elnr).GetType());
+    SurfaceElementIndex sei = IndexBASE<SurfaceElementIndex>() +(elnr-1);    
     eledges.SetSize (ned);
     for (int i = 0; i < ned; i++)
-      eledges[i] = surfedges.Get(elnr)[i]+1;
+      // eledges[i] = surfedges.Get(elnr)[i]+1;
+      eledges[i] = surfedges[sei][i]+1;
   }
 
   void MeshTopology :: GetEdges (SurfaceElementIndex elnr, NgArray<int> & eledges) const
@@ -2030,6 +2053,7 @@ namespace netgen
 
   int MeshTopology :: GetSurfaceElementEdges (int elnr, int * eledges, int * orient) const
   {
+    SurfaceElementIndex sei = IndexBASE<SurfaceElementIndex>() +(elnr-1);        
     int i;
     if (mesh->GetDimension() == 3 || 1)
       {
@@ -2042,8 +2066,10 @@ namespace netgen
 		eledges[i] = abs (surfedges.Get(elnr)[i]);
 		orient[i] = (surfedges.Get(elnr)[i] > 0 ) ? 1 : -1;
                 */
-		if (surfedges.Get(elnr)[i] == -1) return i;
-		eledges[i] = surfedges.Get(elnr)[i]+1;
+		// if (surfedges.Get(elnr)[i] == -1) return i;
+		// eledges[i] = surfedges.Get(elnr)[i]+1;
+		if (surfedges[sei][i] == -1) return i;
+		eledges[i] = surfedges[sei][i]+1;
 		// orient[i] = (surfedges.Get(elnr)[i].orient) ? -1 : 1;
                 // orient[i] = GetSurfaceElementEdgeOrientation(elnr, i) ? -1 : 1;
                 orient[i] = 1;
@@ -2058,8 +2084,10 @@ namespace netgen
 		if (!surfedges.Get(elnr)[i]) return i;
 		eledges[i] = abs (surfedges.Get(elnr)[i]);
                 */
-		if (surfedges.Get(elnr)[i] == -1) return i;
-		eledges[i] = surfedges.Get(elnr)[i]+1;
+		// if (surfedges.Get(elnr)[i] == -1) return i;
+		// eledges[i] = surfedges.Get(elnr)[i]+1;
+		if (surfedges[sei][i] == -1) return i;
+		eledges[i] = surfedges[sei][i]+1;
 	      }
 	  }
 	return 4;
diff --git a/libsrc/meshing/topology.hpp b/libsrc/meshing/topology.hpp
index b8a00268..fcb1b967 100644
--- a/libsrc/meshing/topology.hpp
+++ b/libsrc/meshing/topology.hpp
@@ -32,9 +32,9 @@ class MeshTopology
   Array<std::array<PointIndex,2>> edge2vert;
   Array<std::array<PointIndex,4>> face2vert;
 
-  NgArray<std::array<T_EDGE,12>> edges;
-  NgArray<std::array<T_FACE,6>> faces;
-  NgArray<std::array<T_EDGE,4>> surfedges;
+  Array<std::array<T_EDGE,12>, ElementIndex> edges;
+  Array<std::array<T_FACE,6>, ElementIndex> faces;
+  Array<std::array<T_EDGE,4>, SurfaceElementIndex> surfedges;
   
   NgArray<T_EDGE> segedges;
   NgArray<T_FACE> surffaces;
@@ -170,11 +170,12 @@ public:
 
   int GetSurfaceElementEdges (int elnr, int * edges, int * orient) const;
 
-  const T_EDGE * GetElementEdgesPtr (int elnr) const { return &edges[elnr][0]; }
+  [[deprecated("use GetEdges(ElementIndex) instead")]]
+  const T_EDGE * GetElementEdgesPtr (int elnr) const { return &edges[IndexBASE<ElementIndex>()+elnr][0]; }
   const T_EDGE * GetSurfaceElementEdgesPtr (int selnr) const { return &surfedges[selnr][0]; }
   const T_EDGE * GetSegmentElementEdgesPtr (int selnr) const { return &segedges[selnr]; }
 
-  const T_FACE * GetElementFacesPtr (int elnr) const { return &faces[elnr][0]; }
+  const T_FACE * GetElementFacesPtr (int elnr) const { return &faces[IndexBASE<ElementIndex>()+elnr][0]; }
   const T_FACE * GetSurfaceElementFacesPtr (int selnr) const { return &surffaces[selnr]; }