diff --git a/libsrc/core/array.hpp b/libsrc/core/array.hpp
index c16cc3bc..d523d96b 100644
--- a/libsrc/core/array.hpp
+++ b/libsrc/core/array.hpp
@@ -59,6 +59,26 @@ namespace ngcore
   }
 
 
+  namespace detail
+  {
+
+    // Type trait to check if a class implements a 'range_type Range()' function
+    template<typename T>
+    struct has_Range
+    {
+    private:
+      template<typename T2>
+      static constexpr auto check(T2*) ->
+        std::enable_if<std::declval<T2>().Range(), std::true_type> { std::true_type(); }
+      template<typename>
+      static constexpr std::false_type check(...);
+      using type = decltype(check<T>(nullptr)); // NOLINT
+    public:
+      NGCORE_API static constexpr bool value = type::value;
+    };
+  }
+  template<typename T>
+  constexpr bool has_range = detail::has_Range<T>::value;
 
   template <typename AO>
   class AOWrapperIterator
@@ -210,7 +230,7 @@ namespace ngcore
 
 
   template <typename  T>
-  constexpr size_t IndexBASE () { return 0; }
+  constexpr T IndexBASE () { return T(0); }
 
   
   template <class T, class IndexType = size_t> class FlatArray;
@@ -293,6 +313,11 @@ namespace ngcore
     return T_Range<T>(a,b);
   }
 
+  template<typename T>
+  NETGEN_INLINE auto Range (const T& ao)
+    -> typename std::enable_if<has_range<T>, decltype(std::declval<T>().Range())>::type
+  { return ao.Range(); }
+
   template <typename T>
   NETGEN_INLINE T_Range<T> Range_impl (T n, std::true_type)
   {
@@ -301,13 +326,15 @@ namespace ngcore
 
   template <typename T>
   NETGEN_INLINE auto Range_impl (const T & ao, std::false_type)
-    -> T_Range<decltype(ao.Size())>
+    -> T_Range<index_type<T>>
   {
-    return T_Range<decltype(ao.Size())> (0, ao.Size());
+    return T_Range<index_type<T>> (IndexBASE<index_type<T>>(),
+                                   IndexBASE<index_type<T>>() + index_type<T>(ao.Size()));
   }
 
   template <typename T>
-  auto Range(const T & x) -> decltype(Range_impl(x, std::is_integral<T>())) {
+  auto Range(const T & x)
+    -> typename std::enable_if<!has_range<T>, decltype(Range_impl(x, std::is_integral<T>()))>::type {
     return Range_impl(x, std::is_integral<T>());
   }
 
@@ -402,7 +429,7 @@ namespace ngcore
   class FlatArray : public BaseArrayObject<FlatArray<T,IndexType> >
   {
   protected:
-    static constexpr size_t BASE = IndexBASE<IndexType>();
+    static constexpr IndexType BASE = IndexBASE<IndexType>();
     /// the size
     size_t size;
     /// the data
@@ -500,9 +527,9 @@ namespace ngcore
       return data[i-BASE]; 
     }
   
-    NETGEN_INLINE T_Range<size_t> Range () const
+    NETGEN_INLINE T_Range<index_type> Range () const
     {
-      return T_Range<size_t> (BASE, Size()+BASE);
+      return T_Range<index_type> (BASE, size+BASE);
     }
     
     NETGEN_INLINE const CArray<T> Addr (size_t pos) const
diff --git a/libsrc/core/utils.hpp b/libsrc/core/utils.hpp
index 1eb6d7bd..7f38c948 100644
--- a/libsrc/core/utils.hpp
+++ b/libsrc/core/utils.hpp
@@ -126,6 +126,26 @@ namespace ngcore
       return current;
   }
 
+  namespace detail
+  {
+    template<typename T>
+    struct IndexTypeHelper
+    {
+    private:
+      template<typename T2>
+      static constexpr auto check(T2* t) -> typename T2::index_type { return *t; }
+      static constexpr auto check(...) -> decltype(std::declval<T>().Size())
+      { return decltype(std::declval<T>().Size())(); }
+    public:
+      using type = decltype(check((T*) nullptr)); // NOLINT
+    };
+  
+  } // namespace detail
+
+  // Get index type of object. If object has a typedef index_type this type is returned
+  // else decltype(obj.Size()) is returned.
+  template<typename T>
+  using index_type = typename detail::IndexTypeHelper<T>::type;
 } // namespace ngcore
 
 #endif // NETGEN_CORE_UTILS_HPP
diff --git a/libsrc/general/ngarray.hpp b/libsrc/general/ngarray.hpp
index 213a24c2..f3eea9d4 100644
--- a/libsrc/general/ngarray.hpp
+++ b/libsrc/general/ngarray.hpp
@@ -84,6 +84,7 @@ namespace netgen
     T * data;
   public:
     typedef T TELEM;
+    using index_type = TIND;
 
     /// provide size and memory
     NgFlatArray (size_t asize, T * adata) 
diff --git a/libsrc/meshing/meshtype.hpp b/libsrc/meshing/meshtype.hpp
index 3f0cf751..0d22722d 100644
--- a/libsrc/meshing/meshtype.hpp
+++ b/libsrc/meshing/meshtype.hpp
@@ -162,7 +162,7 @@ namespace netgen
     PointIndex & operator= (const PointIndex&) = default;
     PointIndex & operator= (PointIndex&&) = default;
      
-    PointIndex (int ai) : i(ai)
+    constexpr PointIndex (int ai) : i(ai)
     {
 #ifdef DEBUG
       if (ai < PointIndex::BASE)
@@ -172,7 +172,7 @@ namespace netgen
     }
     constexpr PointIndex (t_invalid inv) : i(PointIndex::BASE-1) { ; }
     // PointIndex & operator= (const PointIndex &ai) { i = ai.i; return *this; }
-    operator int () const { return i; }
+    constexpr operator int () const { return i; }
     PointIndex operator++ (int) { PointIndex hi(*this); i++; return hi; }
     PointIndex operator-- (int) { PointIndex hi(*this); i--; return hi; }
     PointIndex operator++ () { i++; return *this; }
@@ -180,9 +180,9 @@ namespace netgen
     void Invalidate() { i = PointIndex::BASE-1; }
     bool IsValid() const { return i != PointIndex::BASE-1; }
 #ifdef BASE0
-    enum { BASE = 0 };
+    static constexpr size_t BASE = 0;
 #else
-    enum { BASE = 1 };
+    static constexpr size_t BASE = 1;
 #endif  
 
     void DoArchive (Archive & ar) { ar & i; }
@@ -193,7 +193,7 @@ namespace netgen
 namespace ngcore
 {
   template<> 
-  constexpr size_t IndexBASE<netgen::PointIndex> () { return netgen::PointIndex::BASE; }
+  constexpr netgen::PointIndex IndexBASE<netgen::PointIndex> () { return netgen::PointIndex(netgen::PointIndex::BASE); }
 }
 
 namespace netgen
@@ -255,14 +255,14 @@ namespace netgen
     int i;
   public:
     SurfaceElementIndex () = default;
-    SurfaceElementIndex (int ai) : i(ai) { ; }
+    constexpr SurfaceElementIndex (int ai) : i(ai) { ; }
     /*
     SurfaceElementIndex & operator= (const SurfaceElementIndex & ai) 
     { i = ai.i; return *this; }
     */
     SurfaceElementIndex & operator= (const SurfaceElementIndex & ai) = default;
     SurfaceElementIndex & operator= (int ai) { i = ai; return *this; }
-    operator int () const { return i; }
+    constexpr operator int () const { return i; }
     SurfaceElementIndex operator++ (int) { SurfaceElementIndex hi(*this); i++; return hi; }
     SurfaceElementIndex operator-- (int) { SurfaceElementIndex hi(*this); i--; return hi; }
     SurfaceElementIndex & operator++ () { ++i; return *this; }
diff --git a/tests/catch/CMakeLists.txt b/tests/catch/CMakeLists.txt
index 3b727af8..29f234ff 100644
--- a/tests/catch/CMakeLists.txt
+++ b/tests/catch/CMakeLists.txt
@@ -26,6 +26,7 @@ macro(add_unit_test name sources)
 endmacro()
 
 add_unit_test(archive archive.cpp)
+add_unit_test(array array.cpp)
 add_unit_test(symboltable symboltable.cpp)
 add_unit_test(utils utils.cpp)
 add_unit_test(version version.cpp)
diff --git a/tests/catch/array.cpp b/tests/catch/array.cpp
index 471806c0..e7564324 100644
--- a/tests/catch/array.cpp
+++ b/tests/catch/array.cpp
@@ -1,9 +1,11 @@
 
 #include "catch.hpp"
-#include <array.hpp>
+#include <core/array.hpp>
 using namespace ngcore;
 using namespace std;
 
+#include "meshing.hpp"
+
 
 TEST_CASE("Array")
 {
@@ -15,13 +17,14 @@ TEST_CASE("Array")
 #endif // DEBUG
   Array<double> a_initlst = { 1., 2., 3.};
   CHECK(a_initlst[1] == 2.);
-  CHECK(a_initlst.size() == 3);
+  CHECK(a_initlst.Size() == 3);
   FlatArray fa_a = a_initlst;
   CHECK(typeid(fa_a) == typeid(FlatArray<double>));
-  CHECK(fa_a.size() == 3);
-  CHECK(a.Last() == 3.);
-  a.DeleteLast();
-  CHECK(a.Last() == 2. && a.Size() == 2);
+  CHECK(fa_a.Size() == 3);
+  CHECK(fa_a.Last() == 3.);
+  a_initlst.DeleteLast();
+  CHECK(a_initlst.Last() == 2.);
+  CHECK(a_initlst.Size() == 2);
 #ifdef DEBUG
   CHECK_THROWS_AS(fa_a[5], RangeException);
 #endif // DEBUG
@@ -34,4 +37,25 @@ TEST_CASE("Array")
     CHECK(val == 2);
   }
   CHECK(count == 4);
+
+  // range tests
+  CHECK(typeid(array.Range()) == typeid(T_Range<size_t>));
+  Array<int, int> intarray;
+  CHECK(typeid(intarray.Range()) == typeid(T_Range<int>));
+  CHECK(typeid(Range(intarray)) == typeid(T_Range<int>));
+  int i = 0;
+  for(auto j : Range(b))
+    CHECK(j == i++);
+  i = 0;
+  for(auto j : b.Range())
+    CHECK(j == i++);
+
+  // pointindex is still 1 based
+  Array<double, netgen::PointIndex> piarray(2);
+  i = 1;
+  for(auto j : Range(piarray))
+    CHECK(j == i++);
+  i = 1;
+  for(auto j : piarray.Range())
+    CHECK(j == i++);
 }