diff --git a/libsrc/interface/writeuser.cpp b/libsrc/interface/writeuser.cpp
index 5ad3c674..522aebe7 100644
--- a/libsrc/interface/writeuser.cpp
+++ b/libsrc/interface/writeuser.cpp
@@ -427,7 +427,7 @@ void WriteSTLExtFormat (const Mesh & mesh,
 
       for(int faceNr = 1;faceNr <= faceBCMapping.EntrySize(bcInd); faceNr++)
       {
-          NgArray<SurfaceElementIndex> faceSei;
+        Array<SurfaceElementIndex> faceSei;
           mesh.GetSurfaceElementsOfFace(faceBCMapping.Get(bcInd,faceNr),faceSei);
 
           for (int i = 0; i < faceSei.Size(); i++)
diff --git a/libsrc/meshing/bcfunctions.cpp b/libsrc/meshing/bcfunctions.cpp
index 1bede57c..d430a3a7 100644
--- a/libsrc/meshing/bcfunctions.cpp
+++ b/libsrc/meshing/bcfunctions.cpp
@@ -341,9 +341,9 @@ namespace netgen
       // Extract the number of surface elements having a given colour
       // And save this number into an array for later sorting
       for(int face_index = 1; face_index <= nfd; face_index++)
-      {
-         NgArray<SurfaceElementIndex> se_face;
-
+        {
+          Array<SurfaceElementIndex> se_face;
+          
          mesh.GetSurfaceElementsOfFace(face_index, se_face);
 
          Vec3d face_colour;
diff --git a/libsrc/meshing/improve2.cpp b/libsrc/meshing/improve2.cpp
index 74a322a3..7afc65dc 100644
--- a/libsrc/meshing/improve2.cpp
+++ b/libsrc/meshing/improve2.cpp
@@ -74,16 +74,24 @@ namespace netgen
     NgProfiler::StartTimer (timerstart);
 
 
-    NgArray<SurfaceElementIndex> seia;
+    Array<SurfaceElementIndex> seia;
     mesh.GetSurfaceElementsOfFace (faceindex, seia);
 
+    /*
     for (int i = 0; i < seia.Size(); i++)
       if (mesh[seia[i]].GetNP() != 3)
 	{
 	  GenericImprove (mesh);
 	  return;
 	}
-
+    */
+    for (SurfaceElementIndex sei : seia)
+      if (mesh[sei].GetNP() != 3)
+	{
+	  GenericImprove (mesh);
+	  return;
+	}
+      
     int surfnr = mesh.GetFaceDescriptor (faceindex).SurfNr();
 
     NgArray<Neighbour> neighbors(mesh.GetNSE());
@@ -356,8 +364,8 @@ namespace netgen
 			Element2d sw2 (pi3, pi4, pi2);
 
 			int legal1 = 
-			  mesh.LegalTrig (mesh.SurfaceElement (t1)) + 
-			  mesh.LegalTrig (mesh.SurfaceElement (t2));
+			  mesh.LegalTrig (mesh[t1]) + 
+			  mesh.LegalTrig (mesh[t2]);
 			int legal2 = 
 			  mesh.LegalTrig (sw1) + mesh.LegalTrig (sw2);
 
@@ -438,35 +446,24 @@ namespace netgen
     static int timerstart1 = NgProfiler::CreateTimer ("Combineimprove 2D start1");
     NgProfiler::StartTimer  (timerstart1);
 
-
-
-    // int i, j, k, l;
-    // PointIndex pi;
-    // SurfaceElementIndex sei;
-
-
-    NgArray<SurfaceElementIndex> seia;
+    
+    Array<SurfaceElementIndex> seia;
     mesh.GetSurfaceElementsOfFace (faceindex, seia);
 
 
-    for (int i = 0; i < seia.Size(); i++)
-      if (mesh[seia[i]].GetNP() != 3)
+    for (SurfaceElementIndex sei : seia)
+      if (mesh[sei].GetNP() != 3)
 	return;
 
 
-
     int surfnr = 0;
     if (faceindex)
       surfnr = mesh.GetFaceDescriptor (faceindex).SurfNr();
 
 
-    // PointIndex pi1, pi2;
-    // MeshPoint p1, p2, pnew;
-    double bad1, bad2;
     Vec<3> nv;
 
     int np = mesh.GetNP();
-    //int nse = mesh.GetNSE();
 
     TABLE<SurfaceElementIndex,PointIndex::BASE> elementsonnode(np); 
     NgArray<SurfaceElementIndex> hasonepi, hasbothpi;
@@ -639,8 +636,9 @@ namespace netgen
 		  hasonepi.Append (elementsonnode[pi2][k]);
 	      } 
 
-	    bad1 = 0;
+	    double bad1 = 0;
 	    int illegal1 = 0, illegal2 = 0;
+            
 	    for (int k = 0; k < hasonepi.Size(); k++)
 	      {
 		const Element2d & el = mesh[hasonepi[k]];
@@ -665,7 +663,7 @@ namespace netgen
 	    mesh[pi1] = pnew;
 	    mesh[pi2] = pnew;
 
-	    bad2 = 0;
+	    double bad2 = 0;
 	    for (int k = 0; k < hasonepi.Size(); k++)
 	      {
 		Element2d & el = mesh[hasonepi[k]];
@@ -753,24 +751,17 @@ namespace netgen
                 for (SurfaceElementIndex sei2 : elementsonnode[pi2])
 		  {
                     Element2d & el = mesh[sei2];
-		    if(el.IsDeleted()) continue;
-		    elementsonnode.Add (pi1, sei2);
-
-                    /*
-		    bool haspi1 = 0;
-		    for (int l = 0; l < el.GetNP(); l++)
-		      if (el[l] == pi1)
-			haspi1 = true;
-		    if (haspi1) continue;
-                    */
+		    if (el.IsDeleted()) continue;
                     if (el.PNums().Contains(pi1)) continue;
+
+		    elementsonnode.Add (pi1, sei2);
                     
-		    for (int l = 0; l < el.GetNP(); l++)
+		    for (auto l : Range(el.GetNP()))
 		      {
 			if (el[l] == pi2)
 			  {
 			    el[l] = pi1;
-			    el.GeomInfoPi (l+1) = gi;
+			    el.GeomInfo()[l] = gi;
 			  }
 
 			fixed[el[l]] = true;
diff --git a/libsrc/meshing/improve2gen.cpp b/libsrc/meshing/improve2gen.cpp
index ab62bc79..8781782d 100644
--- a/libsrc/meshing/improve2gen.cpp
+++ b/libsrc/meshing/improve2gen.cpp
@@ -430,7 +430,7 @@ namespace netgen
 		  }
 		  
 		for (int j = 0; j < rule.oldels.Size(); j++)
-		  mesh.DeleteSurfaceElement ( elmap[j] );
+		  mesh.Delete (elmap[j]);
 
 		for (int j = 1; j <= pmap.Size(); j++)
 		  nelonnode[pmap.Get(j)] += rule.incelsonnode.Get(j);
diff --git a/libsrc/meshing/improve3.hpp b/libsrc/meshing/improve3.hpp
index 925d12d3..6a94af8c 100644
--- a/libsrc/meshing/improve3.hpp
+++ b/libsrc/meshing/improve3.hpp
@@ -3,7 +3,7 @@
 
 
 extern double CalcTotalBad (const Mesh::T_POINTS & points, 
-			    const Mesh::T_VOLELEMENTS & elements,
+			    const NgArray<Element> & elements,
 			    const MeshingParameters & mp);
 
 
@@ -33,7 +33,7 @@ public:
 
 
   double CalcTotalBad (const Mesh::T_POINTS & points, 
-		       const Mesh::T_VOLELEMENTS & elements)
+		       const NgArray<Element> & elements)
   {
     return netgen::CalcTotalBad (points, elements, mp);
   }
@@ -100,7 +100,7 @@ class JacobianPointFunction : public MinFunction
 {
 public:
   Mesh::T_POINTS & points;
-  const Mesh::T_VOLELEMENTS & elements;
+  const NgArray<Element> & elements;
   TABLE<INDEX> elementsonpoint;
   PointIndex actpind;
 
@@ -109,7 +109,7 @@ public:
   
 public:
   JacobianPointFunction (Mesh::T_POINTS & apoints, 
-			 const Mesh::T_VOLELEMENTS & aelements);
+			 const NgArray<Element> & aelements);
   virtual ~JacobianPointFunction () { ; }
   virtual void SetPointIndex (PointIndex aactpind);
   virtual double Func (const Vector & x) const;
diff --git a/libsrc/meshing/meshclass.cpp b/libsrc/meshing/meshclass.cpp
index 08277f12..b00a630e 100644
--- a/libsrc/meshing/meshclass.cpp
+++ b/libsrc/meshing/meshclass.cpp
@@ -2467,7 +2467,7 @@ namespace netgen
 
     for (int i = 1; i <= GetNSE(); i++)
       {
-        Element2d & sel = surfelements.Elem(i);
+        Element2d & sel = surfelements[i-1];
         bool remove = false;
         for (int j = 1; j <= sel.GetNP(); j++)
           if (frontpoints.Test(sel.PNum(j)))
@@ -2478,9 +2478,9 @@ namespace netgen
 
     for (int i = surfelements.Size(); i >= 1; i--)
       {
-        if (!surfelements.Elem(i).PNum(1).IsValid())
+        if (!surfelements[i-1].PNum(1).IsValid())
           {
-            surfelements.Elem(i) = surfelements.Last();
+            surfelements[i-1] = surfelements.Last();
             surfelements.DeleteLast();
           }
       }
@@ -3270,8 +3270,12 @@ namespace netgen
     for (int i = 0; i < surfelements.Size(); i++)
       if (surfelements[i].IsDeleted())
         {
+          surfelements[i] = surfelements.Last();
+          surfelements.DeleteLast();
+          /*
           surfelements.Delete(i);
           i--;
+          */
         }
 
     for (int i = 0; i < segments.Size(); i++)
@@ -3907,7 +3911,7 @@ namespace netgen
             for (i = 1; i <= nse; i++)
               if (!used.Test(i))
                 {
-                  Element2d & el = surfelements.Elem(i);
+                  Element2d & el = surfelements[i-1];
                   int found = 0, foundrev = 0;
                   for (j = 1; j <= 3; j++)
                     {
@@ -5372,7 +5376,7 @@ namespace netgen
     int np = GetNP();
 
     BitArray usedp(np);
-    NgArray<SurfaceElementIndex> els_of_face;
+    Array<SurfaceElementIndex> els_of_face;
 
     fdi = 1;
     while (fdi <= GetNFD())
@@ -5546,7 +5550,7 @@ namespace netgen
       }
   }
 
-  void Mesh :: GetSurfaceElementsOfFace (int facenr, NgArray<SurfaceElementIndex> & sei) const
+  void Mesh :: GetSurfaceElementsOfFace (int facenr, Array<SurfaceElementIndex> & sei) const
   {
     static int timer = NgProfiler::CreateTimer ("GetSurfaceElementsOfFace");
     NgProfiler::RegionTimer reg (timer);
diff --git a/libsrc/meshing/meshclass.hpp b/libsrc/meshing/meshclass.hpp
index 09d17753..9e31bd0f 100644
--- a/libsrc/meshing/meshclass.hpp
+++ b/libsrc/meshing/meshclass.hpp
@@ -24,9 +24,9 @@ namespace netgen
   {
   public:
     typedef ::netgen::T_POINTS T_POINTS;
-    typedef NgArray<Element, 0, size_t> T_VOLELEMENTS;
+    // typedef NgArray<Element, 0, size_t> T_VOLELEMENTS;
     // typedef NgArray<Element2d, 0, SurfaceElementIndex> T_SURFELEMENTS;
-    typedef NgArray<Element2d, 0, size_t> T_SURFELEMENTS;
+    // typedef NgArray<Element2d, 0, size_t> T_SURFELEMENTS;
 
   private:
     /// point coordinates
@@ -38,9 +38,9 @@ namespace netgen
     /// line-segments at edges
     NgArray<Segment, 0, size_t> segments;
     /// surface elements, 2d-inner elements
-    T_SURFELEMENTS surfelements;
+    Array<Element2d> surfelements;
     /// volume elements
-    T_VOLELEMENTS volelements;
+    NgArray<Element> volelements;
     /// points will be fixed forever
     NgArray<PointIndex> lockedpoints;
 
@@ -263,10 +263,6 @@ namespace netgen
     const Segment & operator[] (SegmentIndex si) const { return segments[si]; }
     Segment & operator[] (SegmentIndex si) { return segments[si]; }
 
-    /*
-    const NgArray<Segment> & LineSegments() const { return segments; }
-    NgArray<Segment> & LineSegments() { return segments; }
-    */
     const auto & LineSegments() const { return segments; }
     auto & LineSegments() { return segments; }
     
@@ -276,30 +272,46 @@ namespace netgen
     // write to pre-allocated container, thread-safe
     DLL_HEADER void SetSurfaceElement (SurfaceElementIndex sei, const Element2d & el);
     
-    // [[deprecated("Use DeleteSurfaceElement(SurfaceElementIndex) instead of int !")]]
+    [[deprecated("Use Delete(SurfaceElementIndex) instead of int !")]]
     void DeleteSurfaceElement (int eli)
-    { 
+    {
+      /*
       surfelements.Elem(eli).Delete();
       surfelements.Elem(eli).PNum(1).Invalidate();
       surfelements.Elem(eli).PNum(2).Invalidate();
       surfelements.Elem(eli).PNum(3).Invalidate();
+      */
+      surfelements[eli-1].Delete();
+      surfelements[eli-1].PNum(1).Invalidate();
+      surfelements[eli-1].PNum(2).Invalidate();
+      surfelements[eli-1].PNum(3).Invalidate();
       timestamp = NextTimeStamp();
     }
 
+    [[deprecated("Use Delete(SurfaceElementIndex) instead !")]]        
     void DeleteSurfaceElement (SurfaceElementIndex eli)
     {
       for (auto & p : surfelements[eli].PNums()) p.Invalidate();
       surfelements[eli].Delete();
       timestamp = NextTimeStamp();
     }
+    
+    void Delete (SurfaceElementIndex eli)
+    {
+      for (auto & p : surfelements[eli].PNums()) p.Invalidate();
+      surfelements[eli].Delete();
+      timestamp = NextTimeStamp();
+    }
 
     auto GetNSE () const { return surfelements.Size(); }
 
     // [[deprecated("Use SurfaceElement(SurfaceElementIndex) instead of int !")]]    
-    Element2d & SurfaceElement(int i) { return surfelements.Elem(i); }
+    Element2d & SurfaceElement(int i) { return surfelements[i-1]; }
     // [[deprecated("Use SurfaceElement(SurfaceElementIndex) instead of int !")]]        
-    const Element2d & SurfaceElement(int i) const { return surfelements.Get(i); }
+    const Element2d & SurfaceElement(int i) const { return surfelements[i-1]; }
+    [[deprecated("Use mesh[](SurfaceElementIndex) instead !")]]
     Element2d & SurfaceElement(SurfaceElementIndex i) { return surfelements[i]; }
+    [[deprecated("Use mesh[](SurfaceElementIndex) instead !")]]
     const Element2d & SurfaceElement(SurfaceElementIndex i) const { return surfelements[i]; }
 
     const Element2d & operator[] (SurfaceElementIndex ei) const
@@ -307,12 +319,12 @@ namespace netgen
     Element2d & operator[] (SurfaceElementIndex ei)
     { return surfelements[ei]; }
 
-    const T_SURFELEMENTS & SurfaceElements() const { return surfelements; }
-    T_SURFELEMENTS & SurfaceElements() { return surfelements; }
+    const auto & SurfaceElements() const { return surfelements; }
+    auto & SurfaceElements() { return surfelements; }
 
   
     DLL_HEADER void RebuildSurfaceElementLists ();
-    DLL_HEADER void GetSurfaceElementsOfFace (int facenr, NgArray<SurfaceElementIndex> & sei) const;
+    DLL_HEADER void GetSurfaceElementsOfFace (int facenr, Array<SurfaceElementIndex> & sei) const;
 
     DLL_HEADER ElementIndex AddVolumeElement (const Element & el);
     // write to pre-allocated container, thread-safe
@@ -324,15 +336,13 @@ namespace netgen
     Element & VolumeElement(int i) { return volelements.Elem(i); }
     // [[deprecated("Use VolumeElement(ElementIndex) instead of int !")]]        
     const Element & VolumeElement(int i) const { return volelements.Get(i); }
+    [[deprecated("Use mesh[](VolumeElementIndex) instead !")]]
     Element & VolumeElement(ElementIndex i) { return volelements[i]; }
+    [[deprecated("Use mesh[](VolumeElementIndex) instead !")]]
     const Element & VolumeElement(ElementIndex i) const { return volelements[i]; }
 
-    const Element & operator[] (ElementIndex ei) const 
-    { return volelements[ei]; }
-    Element & operator[] (ElementIndex ei)
-    { return volelements[ei]; }
-
-
+    const Element & operator[] (ElementIndex ei) const { return volelements[ei]; }
+    Element & operator[] (ElementIndex ei) { return volelements[ei]; }
 
     ELEMENTTYPE ElementType (ElementIndex i) const 
     { return (volelements[i].flags.fixed) ? FIXEDELEMENT : FREEELEMENT; }
diff --git a/libsrc/meshing/meshing2.cpp b/libsrc/meshing/meshing2.cpp
index 01303834..7780b545 100644
--- a/libsrc/meshing/meshing2.cpp
+++ b/libsrc/meshing/meshing2.cpp
@@ -305,7 +305,7 @@ namespace netgen
 	  }
       }
     */
-    NgArray<SurfaceElementIndex> seia;
+    Array<SurfaceElementIndex> seia;
     mesh.GetSurfaceElementsOfFace (facenr, seia);
     for (int i = 0; i < seia.Size(); i++)
       {
diff --git a/libsrc/meshing/parallelmesh.cpp b/libsrc/meshing/parallelmesh.cpp
index b7930162..bfb7e3a7 100644
--- a/libsrc/meshing/parallelmesh.cpp
+++ b/libsrc/meshing/parallelmesh.cpp
@@ -749,9 +749,9 @@ namespace netgen
 
     auto & self = const_cast<Mesh&>(*this);
     self.points = T_POINTS(0);
-    self.surfelements = T_SURFELEMENTS(0);
-    self.volelements = T_VOLELEMENTS(0);
-    self.segments = NgArray<Segment, 0, size_t>(0);
+    self.surfelements = Array<Element2d>(0);
+    self.volelements = NgArray<Element>(0);
+    self.segments = NgArray<Segment>(0);
     self.lockedpoints = NgArray<PointIndex>(0);
     auto cleanup_ptr = [](auto & ptr) {
       if (ptr != nullptr) {
diff --git a/libsrc/meshing/python_mesh.cpp b/libsrc/meshing/python_mesh.cpp
index 16f1b269..4646e3b4 100644
--- a/libsrc/meshing/python_mesh.cpp
+++ b/libsrc/meshing/python_mesh.cpp
@@ -740,11 +740,11 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
     .def_property("dim", &Mesh::GetDimension, &Mesh::SetDimension)
 
     .def("Elements3D", 
-         static_cast<NgArray<Element,0,size_t>&(Mesh::*)()> (&Mesh::VolumeElements),
+         static_cast<NgArray<Element>&(Mesh::*)()> (&Mesh::VolumeElements),
          py::return_value_policy::reference)
 
     .def("Elements2D", 
-         static_cast<NgArray<Element2d,0,size_t>&(Mesh::*)()> (&Mesh::SurfaceElements),
+         static_cast<Array<Element2d>&(Mesh::*)()> (&Mesh::SurfaceElements),
          py::return_value_policy::reference)
 
     .def("Elements1D", 
@@ -813,7 +813,7 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
     .def ("DeleteSurfaceElement",
           FunctionPointer ([](Mesh & self, SurfaceElementIndex i)
                            {
-                             return self.DeleteSurfaceElement (i);
+                             return self.Delete(i);
                            }))
     
     .def ("Compress", FunctionPointer ([](Mesh & self)
diff --git a/libsrc/meshing/smoothing2.5.cpp b/libsrc/meshing/smoothing2.5.cpp
index 351771e6..c20daf24 100644
--- a/libsrc/meshing/smoothing2.5.cpp
+++ b/libsrc/meshing/smoothing2.5.cpp
@@ -51,7 +51,7 @@ namespace netgen
     int i, j, k;
     SurfaceElementIndex sei;
 
-    NgArray<SurfaceElementIndex> seia;
+    Array<SurfaceElementIndex> seia;
     mesh.GetSurfaceElementsOfFace (faceindex, seia);
 
     /*
diff --git a/libsrc/meshing/smoothing2.cpp b/libsrc/meshing/smoothing2.cpp
index ec245b16..5b91dc3a 100644
--- a/libsrc/meshing/smoothing2.cpp
+++ b/libsrc/meshing/smoothing2.cpp
@@ -741,13 +741,13 @@ namespace netgen
     Opti2dLocalData ld;
 
 
-    NgArray<SurfaceElementIndex> seia;
+    Array<SurfaceElementIndex> seia;
     mesh.GetSurfaceElementsOfFace (faceindex, seia);
     bool mixed = 0;
-    for (int i = 0; i < seia.Size(); i++)
-      if (mesh[seia[i]].GetNP() != 3)
+    for (auto sei : seia)
+      if (mesh[sei].GetNP() != 3)
 	{
-	  mixed = 1;
+	  mixed = true;
 	  break;
 	}
 
diff --git a/libsrc/meshing/smoothing3.cpp b/libsrc/meshing/smoothing3.cpp
index 18e7077f..f84f728c 100644
--- a/libsrc/meshing/smoothing3.cpp
+++ b/libsrc/meshing/smoothing3.cpp
@@ -300,7 +300,7 @@ namespace netgen
   {
   public:
     Mesh::T_POINTS & points;
-    const Mesh::T_VOLELEMENTS & elements;
+    const NgArray<Element> & elements;
     TABLE<int,PointIndex::BASE> elementsonpoint;
     const MeshingParameters & mp;
     PointIndex actpind;
@@ -308,7 +308,7 @@ namespace netgen
   
   public:
     PointFunction (Mesh::T_POINTS & apoints, 
-		   const Mesh::T_VOLELEMENTS & aelements,
+		   const NgArray<Element> & aelements,
 		   const MeshingParameters & amp);
     virtual ~PointFunction () { ; }
     virtual void SetPointIndex (PointIndex aactpind);
@@ -323,7 +323,7 @@ namespace netgen
 
 
   PointFunction :: PointFunction (Mesh::T_POINTS & apoints, 
-				  const Mesh::T_VOLELEMENTS & aelements,
+				  const NgArray<Element> & aelements,
 				  const MeshingParameters & amp)
     : points(apoints), elements(aelements), elementsonpoint(apoints.Size()), mp(amp)
   {
@@ -477,7 +477,7 @@ namespace netgen
     DenseMatrix m;
   public:
     CheapPointFunction (Mesh::T_POINTS & apoints, 
-			const Mesh::T_VOLELEMENTS & aelements,
+			const NgArray<Element> & aelements,
 			const MeshingParameters & amp);
     virtual void SetPointIndex (PointIndex aactpind);
     virtual double PointFunctionValue (const Point<3> & pp) const;
@@ -486,7 +486,7 @@ namespace netgen
 
 
   CheapPointFunction :: CheapPointFunction (Mesh::T_POINTS & apoints, 
-					    const Mesh::T_VOLELEMENTS & aelements,
+					    const NgArray<Element> & aelements,
 					    const MeshingParameters & amp)
     : PointFunction (apoints, aelements, amp)
   {
@@ -920,7 +920,7 @@ double Opti3EdgeMinFunction :: FuncGrad (const Vector & x, Vector & grad) const
 
 
 double CalcTotalBad (const Mesh::T_POINTS & points, 
-		     const Mesh::T_VOLELEMENTS & elements,
+		     const NgArray<Element> & elements,
 		     const MeshingParameters & mp)
 {
   static Timer t("CalcTotalBad"); RegionTimer reg(t);
@@ -985,13 +985,13 @@ int WrongOrientation (const Mesh::T_POINTS & points, const Element & el)
 // {
 // public:
 //   Mesh::T_POINTS & points;
-//   const Mesh::T_VOLELEMENTS & elements;
+//   const NgArray<Element> & elements;
 //   TABLE<INDEX> elementsonpoint;
 //   PointIndex actpind;
   
 // public:
 //   JacobianPointFunction (Mesh::T_POINTS & apoints, 
-// 			 const Mesh::T_VOLELEMENTS & aelements);
+// 			 const NgArray<Element> & aelements);
   
 //   virtual void SetPointIndex (PointIndex aactpind);
 //   virtual double Func (const Vector & x) const;
@@ -1002,7 +1002,7 @@ int WrongOrientation (const Mesh::T_POINTS & points, const Element & el)
 
 JacobianPointFunction :: 
 JacobianPointFunction (Mesh::T_POINTS & apoints, 
-		       const Mesh::T_VOLELEMENTS & aelements)
+		       const NgArray<Element> & aelements)
   : points(apoints), elements(aelements), elementsonpoint(apoints.Size())
 {
   INDEX i;
diff --git a/libsrc/meshing/topology.cpp b/libsrc/meshing/topology.cpp
index 850ab808..7f04d2bc 100644
--- a/libsrc/meshing/topology.cpp
+++ b/libsrc/meshing/topology.cpp
@@ -219,7 +219,7 @@ namespace netgen
     
     for (SurfaceElementIndex elnr : top.GetVertexSurfaceElements(v))          
       {
-        const Element2d & el = mesh.SurfaceElement (elnr);
+        const Element2d & el = mesh[elnr];
         
         const ELEMENT_FACE * elfaces = MeshTopology::GetFaces1 (el.GetType());
 	
diff --git a/libsrc/occ/occgenmesh.cpp b/libsrc/occ/occgenmesh.cpp
index 28176009..30e42cc6 100644
--- a/libsrc/occ/occgenmesh.cpp
+++ b/libsrc/occ/occgenmesh.cpp
@@ -831,8 +831,8 @@ namespace netgen
           {
             if (notrys == 1)
               {
-                for (int i = noldsurfel+1; i <= mesh.GetNSE(); i++)
-                  mesh.DeleteSurfaceElement (i);
+                for (SurfaceElementIndex sei = noldsurfel; sei < mesh.GetNSE(); sei++)
+                  mesh.Delete(sei);
 
                 mesh.Compress();
 
diff --git a/libsrc/occ/occpkg.cpp b/libsrc/occ/occpkg.cpp
index 51673451..0c6531ea 100644
--- a/libsrc/occ/occpkg.cpp
+++ b/libsrc/occ/occpkg.cpp
@@ -711,7 +711,7 @@ namespace netgen
 
         for(int i = 1; i <= mesh->GetNFD(); i++)
         {
-           NgArray<SurfaceElementIndex> surfElems;
+          Array<SurfaceElementIndex> surfElems;
            mesh->GetSurfaceElementsOfFace(i,surfElems);
 
            if(ColourMatch(face_colours[colourind],mesh->GetFaceDescriptor(i).SurfColour()))
@@ -735,7 +735,7 @@ namespace netgen
 
         for(int i = 1; i <= mesh->GetNFD(); i++)
         {
-           NgArray<SurfaceElementIndex> surfElems;
+           Array<SurfaceElementIndex> surfElems;
            mesh->GetSurfaceElementsOfFace(i,surfElems);
 
            if(ColourMatch(face_colours[colourind],mesh->GetFaceDescriptor(i).SurfColour()))
@@ -759,7 +759,7 @@ namespace netgen
 
         for(int i = 1; i <= mesh->GetNFD(); i++)
         {
-           NgArray<SurfaceElementIndex> surfElems;
+           Array<SurfaceElementIndex> surfElems;
            mesh->GetSurfaceElementsOfFace(i,surfElems);
 
            if(ColourMatch(face_colours[colourind],mesh->GetFaceDescriptor(i).SurfColour()))
@@ -783,14 +783,14 @@ namespace netgen
 
      if(strcmp(argv[1], "hideonly") == 0)
      {
-        NgArray<Vec3d> face_colours;
+       NgArray<Vec3d> face_colours;
         GetFaceColours(*mesh,face_colours);
 
         int colourind = atoi (argv[2]);
 
         for(int i = 1; i <= mesh->GetNFD(); i++)
         {
-           NgArray<SurfaceElementIndex> surfElems;
+           Array<SurfaceElementIndex> surfElems;
            mesh->GetSurfaceElementsOfFace(i,surfElems);
 
            if(ColourMatch(face_colours[colourind],mesh->GetFaceDescriptor(i).SurfColour()))
diff --git a/libsrc/visualization/vsfieldlines.cpp b/libsrc/visualization/vsfieldlines.cpp
index dce298bd..def0b5f9 100644
--- a/libsrc/visualization/vsfieldlines.cpp
+++ b/libsrc/visualization/vsfieldlines.cpp
@@ -576,7 +576,7 @@ namespace netgen
     shared_ptr<Mesh> mesh = GetMesh();
     if (!mesh) return;
 
-    NgArray<SurfaceElementIndex> elements_2d;
+    Array<SurfaceElementIndex> elements_2d;
     
     //cout << "fieldlines_startface " << fieldlines_startface << endl;
     mesh->GetSurfaceElementsOfFace(fieldlines_startface,elements_2d);
@@ -592,7 +592,7 @@ namespace netgen
 	int i;
     for(i=0; i<elements_2d.Size(); i++)
       {
-	const Element2d & elem = mesh->SurfaceElement(elements_2d[i]);
+	const Element2d & elem = (*mesh)[elements_2d[i]];
 	
 	v1 = mesh->Point(elem[1]) - mesh->Point(elem[0]);
 	v2 = mesh->Point(elem[2]) - mesh->Point(elem[0]);
@@ -613,7 +613,7 @@ namespace netgen
     
     while(startpointsp < startpoints.Size())
       {
-	const Element2d & elem = mesh->SurfaceElement(elements_2d[i]);
+	const Element2d & elem = (*mesh)[elements_2d[i]];
 	
 	int numtri = (elem.GetNV() == 3) ? 1 : 2;
 	
diff --git a/libsrc/visualization/vsmesh.cpp b/libsrc/visualization/vsmesh.cpp
index b784812f..64810b00 100644
--- a/libsrc/visualization/vsmesh.cpp
+++ b/libsrc/visualization/vsmesh.cpp
@@ -1005,7 +1005,7 @@ namespace netgen
 
     int hoplotn = 1 << vispar.subdivisions;
     
-    NgArray<SurfaceElementIndex> seia;
+    Array<SurfaceElementIndex> seia;
 
 
     for (int faceindex = 1; faceindex <= mesh->GetNFD(); faceindex++)