diff --git a/libsrc/meshing/clusters.cpp b/libsrc/meshing/clusters.cpp
index 86045061..dafe81f8 100644
--- a/libsrc/meshing/clusters.cpp
+++ b/libsrc/meshing/clusters.cpp
@@ -27,7 +27,7 @@ namespace netgen
     const MeshTopology & top = mesh.GetTopology();
 
     auto id = this->mesh.GetCommunicator().Rank();
-    auto ntasks = this->mesh.GetCommunicator().Size();
+    // auto ntasks = this->mesh.GetCommunicator().Size();
 
     bool hasedges = top.HasEdges();
     bool hasfaces = top.HasFaces();
diff --git a/libsrc/meshing/improve2gen.cpp b/libsrc/meshing/improve2gen.cpp
index 8781782d..38a15eb4 100644
--- a/libsrc/meshing/improve2gen.cpp
+++ b/libsrc/meshing/improve2gen.cpp
@@ -395,8 +395,8 @@ namespace netgen
 		SelectSurfaceOfPoint (mesh.Point(pmap.Get(1)), pgi.Get(1));
 		GetNormalVector (surfnr, mesh.Point(pmap.Get(1)), pgi.Elem(1), n);
 		  
-		for (int j = 1; j <= rule.oldels.Size(); j++)
-		  bad1 += mesh.SurfaceElement(elmap.Get(j)).CalcJacobianBadness (mesh.Points(), n);
+		for (int j = 0; j < rule.oldels.Size(); j++)
+		  bad1 += mesh[elmap[j]].CalcJacobianBadness (mesh.Points(), n);
 		  
 		// check new element:
 		for (int j = 1; j <= rule.newels.Size(); j++)
diff --git a/libsrc/meshing/improve3.cpp b/libsrc/meshing/improve3.cpp
index 5774caf9..b8cba78f 100644
--- a/libsrc/meshing/improve3.cpp
+++ b/libsrc/meshing/improve3.cpp
@@ -44,7 +44,7 @@ void MeshOptimize3d :: CombineImprove (Mesh & mesh,
   double totalbad = 0;
   for (ElementIndex ei = 0; ei < ne; ei++)
     {
-      if(mesh.GetDimension()==3 && mp.only3D_domain_nr && mp.only3D_domain_nr != mesh.VolumeElement(ei).GetIndex())
+      if(mesh.GetDimension()==3 && mp.only3D_domain_nr && mp.only3D_domain_nr != mesh[ei].GetIndex())
 	continue;
       double elerr = CalcBad (mesh.Points(), mesh[ei], 0);
       totalbad += elerr;
diff --git a/libsrc/meshing/improve3.hpp b/libsrc/meshing/improve3.hpp
index 6a94af8c..c3baf251 100644
--- a/libsrc/meshing/improve3.hpp
+++ b/libsrc/meshing/improve3.hpp
@@ -3,7 +3,7 @@
 
 
 extern double CalcTotalBad (const Mesh::T_POINTS & points, 
-			    const NgArray<Element> & elements,
+			    const Array<Element> & elements,
 			    const MeshingParameters & mp);
 
 
@@ -33,7 +33,7 @@ public:
 
 
   double CalcTotalBad (const Mesh::T_POINTS & points, 
-		       const NgArray<Element> & elements)
+		       const Array<Element> & elements)
   {
     return netgen::CalcTotalBad (points, elements, mp);
   }
@@ -100,7 +100,7 @@ class JacobianPointFunction : public MinFunction
 {
 public:
   Mesh::T_POINTS & points;
-  const NgArray<Element> & elements;
+  const Array<Element> & elements;
   TABLE<INDEX> elementsonpoint;
   PointIndex actpind;
 
@@ -109,7 +109,7 @@ public:
   
 public:
   JacobianPointFunction (Mesh::T_POINTS & apoints, 
-			 const NgArray<Element> & aelements);
+			 const Array<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 1e111a7d..5f364cee 100644
--- a/libsrc/meshing/meshclass.cpp
+++ b/libsrc/meshing/meshclass.cpp
@@ -372,7 +372,7 @@ namespace netgen
     */
 
     volelements[ei]  = el;
-    volelements.Last().flags.illegal_valid = 0;
+    volelements[ei].flags.illegal_valid = 0;
   }
 
 
@@ -3222,7 +3222,7 @@ namespace netgen
 
   double Mesh :: ElementError (int eli, const MeshingParameters & mp) const
   {
-    const Element & el = volelements.Get(eli);
+    const Element & el = volelements[eli-1];
     return CalcTetBadness (points.Get(el[0]), points.Get(el[1]),
                            points.Get(el[2]), points.Get(el[3]), -1, mp);
   }
@@ -3262,7 +3262,7 @@ namespace netgen
       if (volelements[i][0] <= PointIndex::BASE-1 ||
           volelements[i].IsDeleted())
         {
-          volelements.Delete(i);
+          volelements.DeleteElement(i);
           i--;
         }
 
@@ -3277,13 +3277,13 @@ namespace netgen
     for (int i = 0; i < segments.Size(); i++)
       if (segments[i][0] <= PointIndex::BASE-1)
         {
-          segments.Delete(i);
+          segments.DeleteElement(i);
           i--;
         }
 
     for(int i=0; i < segments.Size(); i++)
       if(segments[i].edgenr < 0)
-          segments.Delete(i--);
+          segments.DeleteElement(i--);
 
     pused.Clear();
     for (int i = 0; i < volelements.Size(); i++)
diff --git a/libsrc/meshing/meshclass.hpp b/libsrc/meshing/meshclass.hpp
index 9e31bd0f..7eb4f849 100644
--- a/libsrc/meshing/meshclass.hpp
+++ b/libsrc/meshing/meshclass.hpp
@@ -36,13 +36,13 @@ namespace netgen
     NgMPI_Comm comm;
     
     /// line-segments at edges
-    NgArray<Segment, 0, size_t> segments;
+    Array<Segment> segments;
     /// surface elements, 2d-inner elements
     Array<Element2d> surfelements;
     /// volume elements
-    NgArray<Element> volelements;
+    Array<Element> volelements;
     /// points will be fixed forever
-    NgArray<PointIndex> lockedpoints;
+    Array<PointIndex> lockedpoints;
 
 
     /// surface indices at boundary nodes
@@ -242,8 +242,8 @@ namespace netgen
     DLL_HEADER SegmentIndex AddSegment (const Segment & s);
     void DeleteSegment (int segnr)
     {
-      segments.Elem(segnr)[0].Invalidate();
-      segments.Elem(segnr)[1].Invalidate();
+      segments[segnr-1][0].Invalidate();
+      segments[segnr-1][1].Invalidate();
     }
     /*
     void FullDeleteSegment (int segnr)  // von wem ist das ???
@@ -254,9 +254,9 @@ namespace netgen
 
     int GetNSeg () const { return segments.Size(); }
     // [[deprecated("Use LineSegment(SegmentIndex) instead of int !")]]                
-    Segment & LineSegment(int i) { return segments.Elem(i); }
+    Segment & LineSegment(int i) { return segments[i-1]; }
     // [[deprecated("Use LineSegment(SegmentIndex) instead of int !")]]                    
-    const Segment & LineSegment(int i) const { return segments.Get(i); }
+    const Segment & LineSegment(int i) const { return segments[i-1]; }
 
     Segment & LineSegment(SegmentIndex si) { return segments[si]; }
     const Segment & LineSegment(SegmentIndex si) const { return segments[si]; }
@@ -266,7 +266,7 @@ namespace netgen
     const auto & LineSegments() const { return segments; }
     auto & LineSegments() { return segments; }
     
-    NgArray<Element0d> pointelements;  // only via python interface
+    Array<Element0d> pointelements;  // only via python interface
 
     DLL_HEADER SurfaceElementIndex AddSurfaceElement (const Element2d & el);
     // write to pre-allocated container, thread-safe
@@ -282,32 +282,34 @@ namespace netgen
       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();
+      // 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();
+      // 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 !")]]    
+    [[deprecated("Use SurfaceElement(SurfaceElementIndex) instead of int !")]]    
     Element2d & SurfaceElement(int i) { return surfelements[i-1]; }
-    // [[deprecated("Use SurfaceElement(SurfaceElementIndex) instead of int !")]]        
+    [[deprecated("Use SurfaceElement(SurfaceElementIndex) instead of int !")]]        
     const Element2d & SurfaceElement(int i) const { return surfelements[i-1]; }
     [[deprecated("Use mesh[](SurfaceElementIndex) instead !")]]
     Element2d & SurfaceElement(SurfaceElementIndex i) { return surfelements[i]; }
@@ -333,9 +335,9 @@ namespace netgen
     auto GetNE () const { return volelements.Size(); }
 
     // [[deprecated("Use VolumeElement(ElementIndex) instead of int !")]]    
-    Element & VolumeElement(int i) { return volelements.Elem(i); }
+    Element & VolumeElement(int i) { return volelements[i-1]; }
     // [[deprecated("Use VolumeElement(ElementIndex) instead of int !")]]        
-    const Element & VolumeElement(int i) const { return volelements.Get(i); }
+    const Element & VolumeElement(int i) const { return volelements[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 18be4bef..c37db509 100644
--- a/libsrc/meshing/meshtype.hpp
+++ b/libsrc/meshing/meshtype.hpp
@@ -602,7 +602,7 @@ namespace netgen
     void Delete ()
     {
       deleted = 1;
-      for (PointIndex & p : pnum) p.Invalidate(); 
+      // for (PointIndex & p : pnum) p.Invalidate(); 
     }
     
     bool IsDeleted () const 
diff --git a/libsrc/meshing/parallelmesh.cpp b/libsrc/meshing/parallelmesh.cpp
index bfb7e3a7..85deabc0 100644
--- a/libsrc/meshing/parallelmesh.cpp
+++ b/libsrc/meshing/parallelmesh.cpp
@@ -750,9 +750,9 @@ namespace netgen
     auto & self = const_cast<Mesh&>(*this);
     self.points = T_POINTS(0);
     self.surfelements = Array<Element2d>(0);
-    self.volelements = NgArray<Element>(0);
-    self.segments = NgArray<Segment>(0);
-    self.lockedpoints = NgArray<PointIndex>(0);
+    self.volelements = Array<Element>(0);
+    self.segments = Array<Segment>(0);
+    self.lockedpoints = Array<PointIndex>(0);
     auto cleanup_ptr = [](auto & ptr) {
       if (ptr != nullptr) {
 	delete ptr;
diff --git a/libsrc/meshing/python_mesh.cpp b/libsrc/meshing/python_mesh.cpp
index e72dc282..11b4d85e 100644
--- a/libsrc/meshing/python_mesh.cpp
+++ b/libsrc/meshing/python_mesh.cpp
@@ -741,7 +741,7 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
     .def_property("dim", &Mesh::GetDimension, &Mesh::SetDimension)
 
     .def("Elements3D", 
-         static_cast<NgArray<Element>&(Mesh::*)()> (&Mesh::VolumeElements),
+         static_cast<Array<Element>&(Mesh::*)()> (&Mesh::VolumeElements),
          py::return_value_policy::reference)
 
     .def("Elements2D", 
@@ -749,10 +749,10 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
          py::return_value_policy::reference)
 
     .def("Elements1D", 
-         static_cast<NgArray<Segment,0,size_t>&(Mesh::*)()> (&Mesh::LineSegments),
+         static_cast<Array<Segment>&(Mesh::*)()> (&Mesh::LineSegments),
          py::return_value_policy::reference)
 
-    .def("Elements0D", FunctionPointer([] (Mesh & self) -> NgArray<Element0d>&
+    .def("Elements0D", FunctionPointer([] (Mesh & self) -> Array<Element0d>&
                                        {
                                          return self.pointelements;
                                        } ),
diff --git a/libsrc/meshing/refine.cpp b/libsrc/meshing/refine.cpp
index 7b0e2c3c..97c26d73 100644
--- a/libsrc/meshing/refine.cpp
+++ b/libsrc/meshing/refine.cpp
@@ -185,7 +185,7 @@ namespace netgen
     int oldnf = mesh.GetNSE();
     for (SurfaceElementIndex sei = 0; sei < oldnf; sei++)
       {
-	const Element2d & el = mesh.SurfaceElement(sei);
+	const Element2d & el = mesh[sei];
 
 	switch (el.GetType())
 	  {
@@ -265,7 +265,7 @@ namespace netgen
 		  nel.SetIndex(ind);
 
 		  if (j == 0)
-		    mesh.SurfaceElement(sei) = nel;
+		    mesh[sei] = nel;
 		  else
 		    mesh.AddSurfaceElement(nel);
 		}
@@ -343,7 +343,7 @@ namespace netgen
 		  nel.SetIndex(ind);
 
 		  if (j == 0)
-		    mesh.SurfaceElement(sei) = nel;
+		    mesh[sei] = nel;
 		  else
 		    mesh.AddSurfaceElement(nel);
 		}
diff --git a/libsrc/meshing/secondorder.cpp b/libsrc/meshing/secondorder.cpp
index 25a27dab..cc694136 100644
--- a/libsrc/meshing/secondorder.cpp
+++ b/libsrc/meshing/secondorder.cpp
@@ -115,10 +115,9 @@ namespace netgen
     // refine surface elements
     for (SurfaceElementIndex sei = 0; sei < mesh.GetNSE(); sei++)
       {
-	int j;
-	const Element2d & el = mesh.SurfaceElement(sei);
+	const Element2d & el = mesh[sei];
 
-	int onp(0);
+	int onp = 0;
       
 	Element2d newel(TRIG);
 	newel.SetIndex (el.GetIndex());
@@ -168,11 +167,11 @@ namespace netgen
 	    PrintSysError ("Unhandled element in secondorder:", int(el.GetType()));
 	  }
 
-	for (j = 0; j < onp; j++)
+	for (int j = 0; j < onp; j++)
 	  newel[j] = el[j];
       
 	int nnp = newel.GetNP();
-	for (j = 0; j < nnp-onp; j++)
+	for (int j = 0; j < nnp-onp; j++)
 	  {
 	    int pi1 = newel[betw[j][0]];
 	    int pi2 = newel[betw[j][1]];
@@ -198,7 +197,7 @@ namespace netgen
 	      }
 	  }
       
-	mesh.SurfaceElement(sei) = newel;
+	mesh[sei] = newel;
       }
 
  
diff --git a/libsrc/meshing/smoothing3.cpp b/libsrc/meshing/smoothing3.cpp
index f84f728c..5bc028c6 100644
--- a/libsrc/meshing/smoothing3.cpp
+++ b/libsrc/meshing/smoothing3.cpp
@@ -300,7 +300,7 @@ namespace netgen
   {
   public:
     Mesh::T_POINTS & points;
-    const NgArray<Element> & elements;
+    const Array<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 NgArray<Element> & aelements,
+		   const Array<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 NgArray<Element> & aelements,
+				  const Array<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 NgArray<Element> & aelements,
+			const Array<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 NgArray<Element> & aelements,
+					    const Array<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 NgArray<Element> & elements,
+		     const Array<Element> & elements,
 		     const MeshingParameters & mp)
 {
   static Timer t("CalcTotalBad"); RegionTimer reg(t);
@@ -933,9 +933,9 @@ double CalcTotalBad (const Mesh::T_POINTS & points,
 
   double teterrpow = mp.opterrpow;
 
-  for (int i = 1; i <= elements.Size(); i++)
+  for (int i = 0; i < elements.Size(); i++)
     {
-      elbad = pow (max2(CalcBad (points, elements.Get(i), 0, mp),1e-10),
+      elbad = pow (max2(CalcBad (points, elements[i], 0, mp),1e-10),
 		   1/teterrpow);
 
       int qualclass = int (20 / elbad + 1);
@@ -1002,17 +1002,12 @@ int WrongOrientation (const Mesh::T_POINTS & points, const Element & el)
 
 JacobianPointFunction :: 
 JacobianPointFunction (Mesh::T_POINTS & apoints, 
-		       const NgArray<Element> & aelements)
+		       const Array<Element> & aelements)
   : points(apoints), elements(aelements), elementsonpoint(apoints.Size())
 {
-  INDEX i;
-  int j;
-  
-  for (i = 1; i <= elements.Size(); i++)
-    {
-      for (j = 1; j <= elements.Get(i).NP(); j++)
-	elementsonpoint.Add1 (elements.Get(i).PNum(j), i);  
-    }
+  for (int i = 0; i < elements.Size(); i++)
+      for (int j = 1; j <= elements[i].NP(); j++)
+          elementsonpoint.Add1 (elements[i].PNum(j), i+1);
 
   onplane = false;
 }
@@ -1039,7 +1034,7 @@ double JacobianPointFunction :: Func (const Vector & v) const
   for (j = 1; j <= elementsonpoint.EntrySize(actpind); j++)
     {
       int eli = elementsonpoint.Get(actpind, j);
-      badness += elements.Get(eli).CalcJacobianBadness (points);
+      badness += elements[eli-1].CalcJacobianBadness (points);
     }
   
   points.Elem(actpind) = hp; 
@@ -1072,7 +1067,7 @@ FuncGrad (const Vector & x, Vector & g) const
   for (j = 1; j <= elementsonpoint.EntrySize(actpind); j++)
     {
       int eli = elementsonpoint.Get(actpind, j);
-      const Element & el = elements.Get(eli);
+      const Element & el = elements[eli-1];
 
       lpi = 0;
       for (k = 1; k <= el.GetNP(); k++)
@@ -1080,7 +1075,7 @@ FuncGrad (const Vector & x, Vector & g) const
 	  lpi = k;
       if (!lpi) cerr << "loc point not found" << endl;
 
-      badness += elements.Get(eli).
+      badness += elements[eli-1].
 	CalcJacobianBadnessGradient (points, lpi, hderiv);
 
       for(k=0; k<3; k++)
@@ -1145,7 +1140,7 @@ FuncDeriv (const Vector & x, const Vector & dir, double & deriv) const
   for (j = 1; j <= elementsonpoint.EntrySize(actpind); j++)
     {
       int eli = elementsonpoint.Get(actpind, j);
-      const Element & el = elements.Get(eli);
+      const Element & el = elements[eli-1];
 
       lpi = 0;
       for (k = 1; k <= el.GetNP(); k++)
@@ -1153,7 +1148,7 @@ FuncDeriv (const Vector & x, const Vector & dir, double & deriv) const
 	  lpi = k;
       if (!lpi) cerr << "loc point not found" << endl;
 
-      badness += elements.Get(eli).
+      badness += elements[eli-1].
 	CalcJacobianBadnessDirDeriv (points, lpi, vdir, hderiv);
       deriv += hderiv;
     }
diff --git a/libsrc/occ/occgenmesh.cpp b/libsrc/occ/occgenmesh.cpp
index 30e42cc6..6ba92270 100644
--- a/libsrc/occ/occgenmesh.cpp
+++ b/libsrc/occ/occgenmesh.cpp
@@ -858,8 +858,8 @@ namespace netgen
 
         notrys = 1;
 
-        for (int i = oldnf+1; i <= mesh.GetNSE(); i++)
-          mesh.SurfaceElement(i).SetIndex (k);
+        for (SurfaceElementIndex sei = oldnf; sei < mesh.GetNSE(); sei++)
+          mesh[sei].SetIndex (k);
       }
 
     //      ofstream problemfile("occmesh.rep");
diff --git a/libsrc/visualization/meshdoc.cpp b/libsrc/visualization/meshdoc.cpp
index 58bcef75..d0993272 100644
--- a/libsrc/visualization/meshdoc.cpp
+++ b/libsrc/visualization/meshdoc.cpp
@@ -119,9 +119,6 @@ void VisualSceneMeshDoctor :: DrawScene ()
 
 void VisualSceneMeshDoctor :: BuildScene (int zoomall)
 {
-  int i, j;
- 
-  
   if (zoomall)
     {
       Point3d pmin, pmax;
@@ -159,15 +156,16 @@ void VisualSceneMeshDoctor :: BuildScene (int zoomall)
   
   glDisable (GL_COLOR_MATERIAL);
     
-  for (i = 1; i <= mesh->GetNSE(); i++)
+  for (int i = 1; i <= mesh->GetNSE(); i++)
     {
       glLoadName (i);
 
       // copy to be thread-safe
-      Element2d el = mesh->SurfaceElement (i);
+      // Element2d el = mesh->SurfaceElement (i);
+      Element2d el = (*mesh)[SurfaceElementIndex(i-1)];
 
       int drawel = 1;
-      for (j = 1; j <= el.GetNP(); j++)
+      for (int j = 1; j <= el.GetNP(); j++)
 	{
 	  if (!el.PNum(j))
 	    drawel = 0;
@@ -245,7 +243,7 @@ void VisualSceneMeshDoctor :: BuildScene (int zoomall)
 	    { 3, 5, 4 },
 	    { 4, 5, 6 } };
 
-	  for (j = 0; j < 4; j++)
+	  for (int j = 0; j < 4; j++)
 	    {
 	      const Point3d & lp1 = mesh->Point (el.PNum(trigs[j][0]));
 	      const Point3d & lp2 = mesh->Point (el.PNum(trigs[j][1]));
@@ -275,12 +273,12 @@ void VisualSceneMeshDoctor :: BuildScene (int zoomall)
   glColor3f (0.0f, 0.0f, 0.0f);
   glEnable (GL_COLOR_MATERIAL);
   
-  for (i = 1; i <= mesh->GetNSE(); i++)
+  for (int i = 1; i <= mesh->GetNSE(); i++)
     {
-      Element2d el = mesh->SurfaceElement(i);
+      Element2d el = (*mesh)[SurfaceElementIndex(i-1)];
 
       int drawel = 1;
-      for (j = 1; j <= el.GetNP(); j++)
+      for (int j = 1; j <= el.GetNP(); j++)
 	{
 	  if (!el.PNum(j))
 	    drawel = 0;
@@ -372,7 +370,7 @@ void VisualSceneMeshDoctor :: BuildScene (int zoomall)
 
   glLineWidth (2.0f);
 
-  for (i = 1; i <= mesh->GetNSeg(); i++)
+  for (int i = 1; i <= mesh->GetNSeg(); i++)
     {
       const Segment & seg = mesh->LineSegment(i);
       const Point3d & p1 = mesh->Point(seg[0]);
@@ -498,14 +496,14 @@ void VisualSceneMeshDoctor :: SetMarkEdgeDist (int dist)
 void VisualSceneMeshDoctor :: ClickElement (int elnr)
 {
   selelement = elnr;
-
   int oldlocpi = locpi;
   locpi = locpi % 3 + 1;
   
   if (selelement > 0 && selelement <= mesh->GetNSE())
     {
-      selpoint = mesh->SurfaceElement(selelement).PNum(locpi);
-      selpoint2 = mesh->SurfaceElement(selelement).PNum(oldlocpi);
+      SurfaceElementIndex sei(elnr-1);
+      selpoint = (*mesh)[sei].PNum(locpi);
+      selpoint2 = (*mesh)[sei].PNum(oldlocpi);
       cout << "selpts = " << selpoint << ", " << selpoint2 << endl;
     }
 
diff --git a/libsrc/visualization/vsmesh.cpp b/libsrc/visualization/vsmesh.cpp
index 64810b00..c3949e03 100644
--- a/libsrc/visualization/vsmesh.cpp
+++ b/libsrc/visualization/vsmesh.cpp
@@ -665,10 +665,10 @@ namespace netgen
 
 	for (ElementIndex ei : mesh->VolumeElements().Range())
 	  {
-            if (mesh->VolumeElement(ei).flags.badel)
+            if ((*mesh)[ei].flags.badel)
 	      {
 		// copy to be thread-safe
-		Element el = mesh->VolumeElement (ei);
+		Element el = (*mesh)[ei];
 		if ( (el.GetNP() == 4) || (el.GetNP() == 10))
 		  {
 		    glBegin (GL_LINES);
@@ -749,10 +749,12 @@ namespace netgen
 
         for (SurfaceElementIndex sei : mesh->SurfaceElements().Range())
 	  {
-            Element2d el = mesh->SurfaceElement(sei); // copy to be thread-safe
+            Element2d el = (*mesh)[sei]; // copy to be thread-safe
             if (!el.BadElement())
 	      continue;
 
+            if (el.IsDeleted()) continue;
+            
             bool drawel = true;
             for (int j = 1; j <= el.GetNP(); j++)
 	      if (!el.PNum(j).IsValid())