diff --git a/libsrc/geom2d/csg2d.cpp b/libsrc/geom2d/csg2d.cpp
index f97a0633..09f743b7 100644
--- a/libsrc/geom2d/csg2d.cpp
+++ b/libsrc/geom2d/csg2d.cpp
@@ -250,13 +250,20 @@ IntersectionType IntersectSplineSegment( const Spline & s, const Point<2> & r0,
   double det =  b_*b_ - 4*a_*c_;
   if(det<0.0)
     return NO_INTERSECTION;
-  double sqrt_det =  sqrt(det);
-  double t1 = 1.0/(2*a_) * (-b_ + sqrt_det);
-  double t2 = 1.0/(2*a_) * (-b_ - sqrt_det);
 
-  double t = min(t1,t2);
-  if(t<alpha)
-    t = max(t1,t2);
+  double t;
+
+  if(fabs(a_)>EPSILON)
+  {
+    double sqrt_det =  sqrt(det);
+    double t1 = 1.0/(2*a_) * (-b_ + sqrt_det);
+    double t2 = 1.0/(2*a_) * (-b_ - sqrt_det);
+    t = min(t1, t2);
+    if(t<alpha)
+      t = max(t1,t2);
+  }
+  else // degenerate quadratic equation
+    t = -c_/b_;
 
   if(t+EPSILON<alpha)
     return NO_INTERSECTION;
@@ -288,10 +295,17 @@ IntersectionType IntersectSplineSegment1( const Spline & s, const Point<2> & r0,
   double det =  b_*b_ - 4*a_*c_;
   if(det<0.0)
     return NO_INTERSECTION;
+
   double sqrt_det =  sqrt(det);
   double vbeta[2];
-  vbeta[0] = 1.0/(2*a_) * (-b_ + sqrt_det);
-  vbeta[1] = 1.0/(2*a_) * (-b_ - sqrt_det);
+
+  if(fabs(a_)>EPSILON)
+  {
+    vbeta[0] = 1.0/(2*a_) * (-b_ + sqrt_det);
+    vbeta[1] = 1.0/(2*a_) * (-b_ - sqrt_det);
+  }
+  else // degenrate quadratic equation
+    vbeta[0] = vbeta[1] = -c_/b_;
 
   int dim = fabs(vr[0]) > fabs(vr[1]) ? 0 : 1;
   double valpha[2];
diff --git a/libsrc/meshing/meshclass.hpp b/libsrc/meshing/meshclass.hpp
index e3b45b9b..1de31c97 100644
--- a/libsrc/meshing/meshclass.hpp
+++ b/libsrc/meshing/meshclass.hpp
@@ -38,7 +38,7 @@ namespace netgen
     NgMPI_Comm comm;
     
     /// line-segments at edges
-    Array<Segment> segments;
+    Array<Segment, SegmentIndex> segments;
     /// surface elements, 2d-inner elements
     Array<Element2d, SurfaceElementIndex> surfelements;
     /// volume elements
diff --git a/libsrc/meshing/meshtype.hpp b/libsrc/meshing/meshtype.hpp
index 2f49ad70..8644b6aa 100644
--- a/libsrc/meshing/meshtype.hpp
+++ b/libsrc/meshing/meshtype.hpp
@@ -306,13 +306,17 @@ namespace netgen
     SegmentIndex & operator= (const SegmentIndex & ai) 
     { i = ai.i; return *this; }
     SegmentIndex & operator= (int ai) { i = ai; return *this; }
-    operator int () const { return i; }
+    constexpr operator int () const { return i; }
     SegmentIndex& operator++ () { ++i; return *this; }
     SegmentIndex& operator-- () { --i; return *this; }
     SegmentIndex operator++ (int) { return i++; }
     SegmentIndex operator-- (int) { return i--; }
   };
 
+  inline void SetInvalid (SegmentIndex & id) { id = -1; }
+  inline bool IsInvalid (SegmentIndex & id) { return id == -1; }
+
+
   inline istream & operator>> (istream & ist, SegmentIndex & pi)
   {
     int i; ist >> i; pi = i; return ist;
diff --git a/libsrc/meshing/python_mesh.cpp b/libsrc/meshing/python_mesh.cpp
index 9987ed0f..caee2d35 100644
--- a/libsrc/meshing/python_mesh.cpp
+++ b/libsrc/meshing/python_mesh.cpp
@@ -603,7 +603,7 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
 
   ExportArray<Element,size_t>(m);
   ExportArray<Element2d,SurfaceElementIndex>(m);
-  ExportArray<Segment,size_t>(m);
+  ExportArray<Segment,SegmentIndex>(m);
   ExportArray<Element0d>(m);
   ExportArray<MeshPoint,PointIndex>(m);
   ExportArray<FaceDescriptor>(m);
@@ -823,7 +823,7 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
          py::return_value_policy::reference)
 
     .def("Elements1D", 
-         static_cast<Array<Segment>&(Mesh::*)()> (&Mesh::LineSegments),
+         static_cast<Array<Segment, SegmentIndex>&(Mesh::*)()> (&Mesh::LineSegments),
          py::return_value_policy::reference)
 
     .def("Elements0D", FunctionPointer([] (Mesh & self) -> Array<Element0d>&