Merge branch 'constexpr_experiments' into 'master'

Constexpr experiments

See merge request ngsolve/netgen!692

libsrc/core/array.hpp

@@ -216,6 +216,10 @@ namespace ngcore
   template <typename  T>
   constexpr T IndexBASE () { return T(0); }
 
+  template <typename  T>
+  constexpr T IndexBASE (T ind) { return IndexBASE<T>(); }
+
+  
 
   class IndexFromEnd
   {

libsrc/interface/writegmsh2.cpp

@@ -58,7 +58,7 @@ namespace netgen
       int np = mesh.GetNP();  /// number of points in mesh
       int ne = mesh.GetNE();  /// number of 3D elements in mesh
       int nse = mesh.GetNSE();  /// number of surface elements (BC)
-      int i, j, k, l;
+      // int i, j, k, l;
 
 
       /*
@@ -66,8 +66,8 @@ namespace netgen
       */
 
       if ((ne > 0)
-         && (mesh.VolumeElement(1).GetNP() <= 10)
+          && (mesh.VolumeElements().First().GetNP() <= 10)
-         && (mesh.SurfaceElement(1).GetNP() <= 6))
+          && (mesh.SurfaceElements().First().GetNP() <= 6))
       {
          cout << "Write GMSH v2.xx Format \n";
          cout << "The GMSH v2.xx export is currently available for elements upto 2nd Order\n" << endl;
@@ -86,7 +86,7 @@ namespace netgen
          outfile << "$Nodes\n";
          outfile << np << "\n";
 
-         for (i = 1; i <= np; i++)
+         for (int i = 1; i <= np; i++)
          {
             const Point3d & p = mesh.Point(i);
             outfile << i << " "; /// node number
@@ -101,11 +101,11 @@ namespace netgen
          outfile << "$Elements\n";
          outfile << ne + nse << "\n";  ////  number of elements + number of surfaces BC
 
-         for (i = 1; i <= nse; i++)
+         for (auto sei : Range(mesh.SurfaceElements()))
          {
             int elType = 0;
 
-            Element2d el = mesh.SurfaceElement(i);
+            Element2d el = mesh[sei]; // .SurfaceElement(i);
             if(invertsurf) el.Invert();
             
             if(el.GetNP() == 3) elType = GMSH_TRIG;	//// GMSH Type for a 3 node triangle
@@ -116,7 +116,7 @@ namespace netgen
                return;
             }
 
-            outfile << i;
+            outfile << sei-IndexBASE(sei)+1;  
             outfile << " ";
             outfile << elType;
             outfile << " ";
@@ -125,7 +125,7 @@ namespace netgen
             outfile << mesh.GetFaceDescriptor (el.GetIndex()).BCProperty() << " ";
             /// that means that physical entity = elementary entity (arbitrary approach)
             outfile << mesh.GetFaceDescriptor (el.GetIndex()).BCProperty() << " ";
-            for (j = 1; j <= el.GetNP(); j++)
+            for (int j = 1; j <= el.GetNP(); j++)
             {
                outfile << " ";
                outfile << el.PNum(triGmsh[j]);
@@ -133,12 +133,12 @@ namespace netgen
             outfile << "\n";
          }
 
-
+         for (ElementIndex ei : Range(mesh.VolumeElements()))
-         for (i = 1; i <= ne; i++)
          {
+           int i = ei-IndexBASE(ei)+1;
             int elType = 0;
 
-            Element el = mesh.VolumeElement(i);
+            Element el = mesh[ei];
             if (inverttets) el.Invert();
 
             if(el.GetNP() == 4) elType = GMSH_TET;    //// GMSH Element type for 4 node tetrahedron
@@ -160,7 +160,7 @@ namespace netgen
             outfile << " ";
             outfile << 100000 + el.GetIndex();   /// volume number
             outfile << " ";
-            for (j = 1; j <= el.GetNP(); j++)
+            for (int j = 1; j <= el.GetNP(); j++)
             {
                outfile << " ";
                outfile << el.PNum(tetGmsh[j]);
@@ -193,7 +193,7 @@ namespace netgen
          outfile << "$Nodes\n";
          outfile << np << "\n";
 
-         for (i = 1; i <= np; i++)
+         for (int i = 1; i <= np; i++)
          {
             const Point3d & p = mesh.Point(i);
             outfile << i << " "; /// node number
@@ -207,7 +207,7 @@ namespace netgen
          outfile << "$Elements\n";
          outfile << nse << "\n";
 
-         for (k = 1; k <= nse; k++)
+         for (int k = 1; k <= nse; k++)
          {
             int elType = 0;
 
@@ -232,7 +232,7 @@ namespace netgen
             outfile << mesh.GetFaceDescriptor (el.GetIndex()).BCProperty() << " ";
             /// that means that physical entity = elementary entity (arbitrary approach)
             outfile << mesh.GetFaceDescriptor (el.GetIndex()).BCProperty() << " ";
-            for (l = 1; l <= el.GetNP(); l++)
+            for (int l = 1; l <= el.GetNP(); l++)
             {
                outfile << " ";
                if((elType == GMSH_TRIG) || (elType == GMSH_TRIG6))

libsrc/meshing/bisect.cpp

@@ -5,10 +5,8 @@
 
 #include "meshing.hpp"  // quickfix for parallel 
 
-
 #define noDEBUG
 
-
 namespace netgen
 {
   class MarkedTet

libsrc/meshing/clusters.cpp

@@ -143,19 +143,21 @@ namespace netgen
 	  cluster_reps.Elem(nnums[j]) = nnums[j];
       }
     */
+
+    
     ngcore::ParallelForRange
       (mesh.SurfaceElements().Range(),
        [&] (auto myrange)
        {
          NgArrayMem<int,9> nnums; // , ednums;
-         for (int i_ : myrange)
+         for (SurfaceElementIndex i_ : myrange)
            {
              int i = i_+1;
-             const Element2d & el = mesh.SurfaceElement(i);
+             const Element2d & el = mesh[i_]; // .SurfaceElement(i);
              ELEMENT_TYPE typ = el.GetType();
              
              // top.GetSurfaceElementEdges (i, ednums);
-             auto ednums = top.GetEdges (SurfaceElementIndex(i_));
+             auto ednums = top.GetEdges (i_);
              // cout << "ednums = " << ednums << endl;
              
              int fanum = top.GetSurfaceElementFace (i);             

libsrc/meshing/improve3.cpp

@@ -357,7 +357,7 @@ void MeshOptimize3d :: CombineImprove ()
 {
   static Timer t("MeshOptimize3d::CombineImprove"); RegionTimer reg(t);
   static Timer topt("Optimize");
-  static Timer tsearch("Search");
+  static Timer tsearch("Search-combine");
   static Timer tbuild_elements_table("Build elements table");
 
   mesh.BuildBoundaryEdges(false);
@@ -613,7 +613,7 @@ void MeshOptimize3d :: SplitImprove ()
 {
   static Timer t("MeshOptimize3d::SplitImprove"); RegionTimer reg(t);
   static Timer topt("Optimize");
-  static Timer tsearch("Search");
+  static Timer tsearch("Search-split");
 
   // int np = mesh.GetNP();
   int ne = mesh.GetNE();

libsrc/meshing/meshclass.cpp

@@ -3051,6 +3051,7 @@ namespace netgen
     // bool buggy = false;
     // ofstream bout("buggy.out");
 
+
     for (int i = 1; i <= GetNSE(); i++)
       {
         const Element2d & el = SurfaceElement(i);
@@ -3060,10 +3061,10 @@ namespace netgen
           {
             for (int j = 1; j <= el.GetNP(); j++)
               {
-                INDEX_3 seg (el.PNumMod(j), el.PNumMod(j+1), el.GetIndex());
+                PointIndices<3> seg (el.PNumMod(j), el.PNumMod(j+1), el.GetIndex());
                 // int data;
 
-                if (seg.I1() < PointIndex::BASE || seg.I2() < PointIndex::BASE)
+                if (!seg.I1().IsValid() || !seg.I2().IsValid())
                   cerr << "seg = " << seg << endl;
 
                 if (faceht.Used(seg))
@@ -7299,6 +7300,7 @@ namespace netgen
       {
         int eli0, eli1;
         GetTopology().GetSurface2VolumeElement(sei+1, eli0, eli1);
+        // auto [ei0,ei1] = GetTopology().GetSurface2VolumeElement(sei); // the way to go
         auto & sel = (*this)[sei];
         int face = sel.GetIndex();
         int domin = VolumeElement(eli0).GetIndex();

libsrc/meshing/meshing3.cpp

@@ -9,8 +9,6 @@ double minwithoutother;
 
 
 
-
-
 MeshingStat3d :: MeshingStat3d ()
 {
   cntsucc = cnttrials = cntelem = qualclass = 0;

libsrc/meshing/meshtype.hpp

@@ -27,7 +27,6 @@ namespace netgen
   */
 
 
-
   enum ELEMENT_TYPE : unsigned char { 
     SEGMENT = 1, SEGMENT3 = 2,
     TRIG = 10, QUAD=11, TRIG6 = 12, QUAD6 = 13, QUAD8 = 14,
@@ -149,129 +148,108 @@ namespace netgen
   }
 
 
-
+  template <typename T, typename TIndex, int Base>
-
+  class Index
-  /*
-  class PointIndex
   {
-    int i;
+  public:
+
+    T i;
+    static constexpr int BASE = Base;
+    
   public:
     class t_invalid { public: constexpr t_invalid() = default; };
     static constexpr t_invalid INVALID{};
     
-    PointIndex () = default;
+    constexpr Index () = default;
-    PointIndex (const PointIndex&) = default;
+    constexpr Index (const Index& i2) = default;
-    PointIndex (PointIndex &&) = default;
+    constexpr Index (Index &&) = default;
-    PointIndex & operator= (const PointIndex&) = default;
+    Index & operator= (const Index&) = default;
-    PointIndex & operator= (PointIndex&&) = default;
+    Index & operator= (Index&&) = default;
 
-    constexpr PointIndex (int ai) : i(ai)
+    // private:
+    constexpr Index (T ai) : i(ai)
     {
 #ifdef DEBUG
-      if (ai < PointIndex::BASE)
+      if (ai < Base)
-        cout << "illegal PointIndex, use PointIndex::INVALID instead" << endl;
+        cout << "illegal Index, use Index::INVALID instead" << endl;
-        // throw Exception("illegal PointIndex, use PointIndex::INVALID instead");
 #endif
     }
-    constexpr PointIndex (t_invalid inv) : i(PointIndex::BASE-1) { ; }
-    // PointIndex & operator= (const PointIndex &ai) { i = ai.i; return *this; }
-    constexpr operator const int& () const { return i; }
-    explicit constexpr operator int& () { 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; }
-    PointIndex operator-- () { i--; return *this; }
-    PointIndex operator+= (int add) { i += add; return *this; }
-    void Invalidate() { i = PointIndex::BASE-1; }
-    bool IsValid() const { return i != PointIndex::BASE-1; }
-#ifdef BASE0
-    static constexpr size_t BASE = 0;
-#else
-    static constexpr size_t BASE = 1;
-#endif  
 
-    void DoArchive (Archive & ar) { ar & i; }
+
-  }
+    // friend constexpr netgen::TIndex ngcore::IndexBASE<netgen::TIndex> ();
+    // template <int N> friend class PointIndices;
+
+    /*
+    friend auto operator+ (Index, int) -> TIndex;
+    friend TIndex operator+ (Index, size_t);    
+    friend TIndex operator+ (int, Index);
+    friend TIndex operator+ (size_t, Index);
+    friend constexpr TIndex operator- (Index, int);
+    friend int operator- (Index, Index);
+    friend bool operator< (Index a, Index b);
+    friend bool operator> (Index a, Index b);
+    friend bool operator>= (Index a, Index b);
+    friend bool operator<= (Index a, Index b);
+    friend bool operator== (Index a, Index b);
+    friend bool operator!= (Index a, Index b);
     */
     
-
-  // #define BASE0
-  
-  class PointIndex
-  {
-    int i;
   public:
-    class t_invalid { public: constexpr t_invalid() = default; };
+    constexpr Index (t_invalid inv) : i(long(BASE)-1) { ; }
-    static constexpr t_invalid INVALID{};
+    // TIndex & operator= (const TIndex &ai) { i = ai.i; return *this; }
-    
-    PointIndex () = default;
-    PointIndex (const PointIndex&) = default;
-    PointIndex (PointIndex &&) = default;
-    PointIndex & operator= (const PointIndex&) = default;
-    PointIndex & operator= (PointIndex&&) = default;
-
     // private:
-    constexpr PointIndex (int ai) : i(ai)
+    constexpr operator T () const { return i; }
-    {
+    explicit constexpr operator T& () { return i; }
-#ifdef DEBUG
-      if (ai < PointIndex::BASE)
-        cout << "illegal PointIndex, use PointIndex::INVALID instead" << endl;
-        // throw Exception("illegal PointIndex, use PointIndex::INVALID instead");
-#endif
-    }
-
-    friend constexpr netgen::PointIndex ngcore::IndexBASE<netgen::PointIndex> ();
-    template <int N> friend class PointIndices;
-    
-    friend constexpr PointIndex operator+ (PointIndex, int);
-    friend constexpr PointIndex operator+ (PointIndex, size_t);    
-    friend constexpr PointIndex operator+ (int, PointIndex);
-    friend constexpr PointIndex operator+ (size_t, PointIndex);
-    friend constexpr PointIndex operator- (PointIndex, int);
-    friend constexpr int operator- (PointIndex, PointIndex);
-    friend bool operator< (PointIndex a, PointIndex b);
-    friend bool operator> (PointIndex a, PointIndex b);
-    friend bool operator>= (PointIndex a, PointIndex b);
-    friend bool operator<= (PointIndex a, PointIndex b);
-    friend constexpr bool operator== (PointIndex a, PointIndex b);
-    friend constexpr bool operator!= (PointIndex a, PointIndex b);
-    
   public:
-    constexpr PointIndex (t_invalid inv) : i(long(PointIndex::BASE)-1) { ; }
+    // constexpr operator TIndex() const { return TIndex(i); }
-    // PointIndex & operator= (const PointIndex &ai) { i = ai.i; return *this; }
+    // operator TIndex&() { return static_cast<TIndex&>(*this); }    
-    // private:
+
-    constexpr operator const int& () const { return i; }
+    TIndex operator++ (int) { TIndex hi{*this}; i++; return hi; }
-    explicit constexpr operator int& () { return i; }
+    TIndex operator-- (int) { TIndex hi(*this); i--; return hi; }
-  public:
+    TIndex & operator++ () { i++; return static_cast<TIndex&>(*this); }
-    PointIndex operator++ (int) { PointIndex hi(*this); i++; return hi; }
+    TIndex & operator-- () { i--; return static_cast<TIndex&>(*this); }
-    PointIndex operator-- (int) { PointIndex hi(*this); i--; return hi; }
+    constexpr TIndex operator+= (T add) { i += add; return TIndex{*this}; }
-    PointIndex & operator++ () { i++; return *this; }
+    void Invalidate() { i = long(TIndex::BASE)-1; }
-    PointIndex operator-- () { i--; return *this; }
+    bool IsValid() const { return i+1 != TIndex::BASE; }
-    PointIndex operator+= (int add) { i += add; return *this; }
-    void Invalidate() { i = long(PointIndex::BASE)-1; }
-    bool IsValid() const { return i+1 != PointIndex::BASE; }
     // operator bool() const { return IsValid(); }
-#ifdef BASE0
-    static constexpr size_t BASE = 0;
-#else
-    static constexpr size_t BASE = 1;
-#endif  
 
     void DoArchive (Archive & ar) { ar & i; }
   };
 
-  constexpr inline PointIndex operator+ (PointIndex pi, int i) { return PointIndex(pi.i+i); }
+
-  constexpr inline PointIndex operator+ (PointIndex pi, size_t i) { return PointIndex(pi.i+i); }  
+  template <typename T, typename TIndex, int Base>  
-  constexpr inline PointIndex operator+ (int i, PointIndex pi) { return PointIndex(pi.i+i); }
+  constexpr auto operator+ (Index<T,TIndex,Base> ind, int i) { return TIndex(ind.i+i); }
-  constexpr inline PointIndex operator+ (size_t i, PointIndex pi) { return PointIndex(pi.i+i); }
+  template <typename T, typename TIndex, int Base>    
-  constexpr inline PointIndex operator- (PointIndex pi, int i) { return PointIndex(pi.i-i); }
+  constexpr TIndex operator+ (Index<T,TIndex,Base> pi, size_t i) { return TIndex(pi.i+i); }
-  constexpr inline int operator- (PointIndex pa, PointIndex pb) { return pa.i-pb.i; }
+  template <typename T, typename TIndex, int Base>    
-  inline bool operator< (PointIndex a, PointIndex b) { return a.i-b.i < 0; }
+  constexpr TIndex operator+ (int i, Index<T,TIndex,Base> pi) { return TIndex(pi.i+i); }
-  inline bool operator> (PointIndex a, PointIndex b) { return a.i-b.i > 0; }
+  template <typename T, typename TIndex, int Base>    
-  inline bool operator>= (PointIndex a, PointIndex b) { return a.i-b.i >= 0; }
+  inline TIndex operator+ (size_t i, Index<T,TIndex,Base> pi) { return TIndex(pi.i+i); }
-  inline bool operator<= (PointIndex a, PointIndex b) { return a.i-b.i <= 0; }
+  template <typename T, typename TIndex, int Base>    
-  inline constexpr bool operator== (PointIndex a, PointIndex b) { return a.i == b.i; }
+  constexpr inline auto operator- (Index<T,TIndex,Base> pi, int i) -> TIndex { return TIndex(pi.i-i); }
-  inline constexpr bool operator!= (PointIndex a, PointIndex b) { return a.i != b.i; }
+  template <typename T, typename TIndex, int Base>    
+  constexpr inline int operator- (Index<T,TIndex,Base> pa, Index<T,TIndex,Base> pb) { return pa.i-pb.i; }
+  template <typename T, typename TIndex, int Base>      
+  inline bool operator< (Index<T,TIndex,Base> a, Index<T,TIndex,Base> b) { return a.i-b.i < 0; }
+  template <typename T, typename TIndex, int Base>      
+  inline bool operator> (Index<T,TIndex,Base> a, Index<T,TIndex,Base> b) { return a.i-b.i > 0; }
+  template <typename T, typename TIndex, int Base>      
+  inline bool operator>= (Index<T,TIndex,Base> a, Index<T,TIndex,Base> b) { return a.i-b.i >= 0; }
+  template <typename T, typename TIndex, int Base>      
+  inline bool operator<= (Index<T,TIndex,Base> a, Index<T,TIndex,Base> b) { return a.i-b.i <= 0; }
+  template <typename T, typename TIndex, int Base>      
+  inline bool operator== (Index<T,TIndex,Base> a, Index<T,TIndex,Base> b) { return a.i == b.i; }
+  template <typename T, typename TIndex, int Base>      
+  inline bool operator!= (Index<T,TIndex,Base> a, Index<T,TIndex,Base> b) { return a.i != b.i; }
+  
+
+
+  class PointIndex : public Index<int,PointIndex,1>
+  {
+  public:
+    using Index::Index;
+  };
+
 }
 
 namespace ngcore
@@ -288,16 +266,15 @@ namespace netgen
   {
     // int i; ist >> i; pi = PointIndex(i); return ist;
     int i; ist >> i;
-    // pi = PointIndex(i);
     pi = IndexBASE<PointIndex>()+i-1;
     return ist;
   }
 
   inline ostream & operator<< (ostream & ost, const PointIndex & pi)
   {
+    // return (ost << int(pi));
     int intpi = pi - IndexBASE<PointIndex>() + 1;
     return (ost << intpi);    
-    // return (ost << int(pi));
   }
 
   template <int N> class PointIndices;
@@ -463,23 +440,10 @@ namespace std
 namespace netgen
 {
 
-  class ElementIndex
+  class ElementIndex : public Index<int,ElementIndex,0>
   {
-    int i;
   public:
-    ElementIndex () = default;
+    using Index<int,ElementIndex,0>::Index;
-    constexpr /* explicit */ ElementIndex (int ai) : i(ai) { ; }
-    ElementIndex & operator= (const ElementIndex & ai) { i = ai.i; return *this; }
-    ElementIndex & operator= (int ai) { i = ai; return *this; }
-
-    constexpr /* explicit */  operator int () const { return i; }
-    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)
@@ -492,45 +456,33 @@ namespace netgen
     return (ost << ei-IndexBASE<ElementIndex>());
   }
 
-  inline ElementIndex operator+ (ElementIndex ei, int i) { return ElementIndex { int(ei) + i }; }
-  inline ElementIndex operator+ (size_t s, ElementIndex ei) { return ElementIndex(int(ei) + s); }  
-  inline ElementIndex operator+ (ElementIndex ei, size_t s) { return ElementIndex(int(ei) + s); }  
-  inline bool operator== (ElementIndex ei1, ElementIndex ei2) { return int(ei1) == int(ei2); };
-  inline bool operator!= (ElementIndex ei1, ElementIndex ei2) { return int(ei1) != int(ei2); };
-  inline bool operator< (ElementIndex ei1, ElementIndex ei2) { return int(ei1) < int(ei2); };  
-  inline bool operator> (ElementIndex ei1, ElementIndex ei2) { return int(ei1) > int(ei2); };  
-  inline bool operator>= (ElementIndex ei1, ElementIndex ei2) { return int(ei1) >= int(ei2); };
-  inline bool operator<= (ElementIndex ei1, ElementIndex ei2) { return int(ei1) <= int(ei2); };
-  // these should not be needed:
   
-  inline bool operator== (ElementIndex ei1, int ei2) { return int(ei1) == int(ei2); };  
+  // these should not be needed soon
-  inline bool operator< (size_t s, ElementIndex ei2) { return int(s) < int(ei2); };    
+  inline bool operator== (Index<int,ElementIndex,0> ei1, int ei2) { return int(ei1) == int(ei2); };  
-  inline bool operator< (ElementIndex ei1, size_t s) { return int(ei1) < int(s); };   // should not need
+  inline bool operator< (size_t s, Index<int,ElementIndex,0> ei2) { return int(s) < int(ei2); };    
-  inline bool operator< (ElementIndex ei1, int s) { return int(ei1) < int(s); };   // should not need
+  inline bool operator< ( Index<int,ElementIndex,0> ei1, size_t s) { return int(ei1) < int(s); };   // should not need
-  inline bool operator>= (size_t s, ElementIndex ei2) { return int(s) >= int(ei2); };
+  inline bool operator< ( Index<int,ElementIndex,0> ei1, int s) { return int(ei1) < int(s); };   // should not need
+  inline bool operator>= (size_t s,  Index<int,ElementIndex,0> ei2) { return int(s) >= int(ei2); };
 
 
-  class SurfaceElementIndex
+  class SurfaceElementIndex : public Index<int,SurfaceElementIndex,0>
   {
-    int i;
   public:
-    SurfaceElementIndex () = default;
+    using Index::Index;
-    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; }
-    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; }
-    SurfaceElementIndex & operator-- () { --i; return *this; }
-    SurfaceElementIndex & operator+= (int inc) { i+=inc; return *this; }
-    void DoArchive (Archive & ar) { ar & i; }
   };
 
+  
+  // these should not be needed soon
+  inline bool operator== (Index<int, SurfaceElementIndex,0> ei1, int ei2) { return int(ei1) == int(ei2); };
+  inline bool operator== (int ei2, Index<int, SurfaceElementIndex,0> ei1) { return int(ei1) == int(ei2); };
+  inline bool operator!= (Index<int, SurfaceElementIndex,0> ei1, int ei2) { return int(ei1) != int(ei2); };    
+  inline bool operator< (size_t s, Index<int, SurfaceElementIndex,0> ei2) { return int(s) < int(ei2); };    
+  inline bool operator< (Index<int, SurfaceElementIndex,0> ei1, size_t s) { return int(ei1) < int(s); };   // should not need
+  inline bool operator< (Index<int, SurfaceElementIndex,0> ei1, int s) { return int(ei1) < int(s); };   // should not need
+  inline bool operator>= (size_t s, Index<int, SurfaceElementIndex,0> ei2) { return int(s) >= int(ei2); };
+  inline bool operator>= (Index<int, SurfaceElementIndex,0> ei1, int s) { return int(ei1) >= int(s); };
+
+  
   inline void SetInvalid (SurfaceElementIndex & id) { id = -1; }
   inline bool IsInvalid (SurfaceElementIndex & id) { return id == -1; }
 
@@ -544,20 +496,11 @@ namespace netgen
     return (ost << int(pi));
   }
 
-  class SegmentIndex
+
+  class SegmentIndex : public Index<int,SegmentIndex,0>
   {
-    int i;
   public:
-    SegmentIndex () = default;
+    using Index::Index;
-    constexpr SegmentIndex (int ai) : i(ai) { ; }
-    SegmentIndex & operator= (const SegmentIndex & ai) 
-    { i = ai.i; return *this; }
-    SegmentIndex & operator= (int ai) { i = ai; return *this; }
-    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; }

libsrc/meshing/parallelmesh.cpp

@@ -593,17 +593,16 @@ namespace netgen
 
     DynamicTable<int> elementarrays(elarraysize);
     
-    for (int ei = 1; ei <= GetNE(); ei++)
+    for (ElementIndex ei : VolumeElements().Range())    
       {
 	const Element & el = VolumeElement (ei);
-	// int dest = el.GetPartition();
+        int dest = vol_partition[ei];
-        int dest = vol_partition[ei-1];
         
-	elementarrays.Add (dest, ei);
+	elementarrays.Add (dest, int(ei+1));
 	elementarrays.Add (dest, el.GetIndex());
 	elementarrays.Add (dest, el.GetNP());
-	for (int i = 0; i < el.GetNP(); i++)
+        for (PointIndex pi : el.PNums())
-	  elementarrays.Add (dest, el[i]);
+	  elementarrays.Add (dest, pi);
       }
     tbuildelementtable.Stop();
     

libsrc/meshing/python_mesh.cpp

@@ -268,7 +268,7 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
     .def(py::init<int>())
     .def("__repr__", &ToString<PointIndex>)
     .def("__str__", &ToString<PointIndex>)
-    .def_property_readonly("nr", &PointIndex::operator const int&)
+    .def_property_readonly("nr", &PointIndex::operator int)
     .def("__eq__" , FunctionPointer( [](PointIndex &self, PointIndex &other)
                   { return static_cast<int>(self)==static_cast<int>(other); }) )
     .def("__hash__" , FunctionPointer( [](PointIndex &self ) { return static_cast<int>(self); }) )

libsrc/meshing/refine.cpp

@@ -865,9 +865,10 @@ namespace netgen
 		for (int k = 1; k <= 3; k++)
 		  {
 		    fhelp.Clear();
-		    for (int i = 1; i <= mesh.GetNE(); i++)
+		    // for (int i = 1; i <= mesh.GetNE(); i++) 
+                    for (const Element & el : mesh.VolumeElements())
 		      {
-			const Element & el = mesh.VolumeElement(i);
+			// const Element & el = mesh.VolumeElement(i);
 			int freeel = 0;
 			for (int j = 1; j <= el.GetNP(); j++)
 			  if (free.Test(el.PNum(j)))
@@ -897,19 +898,19 @@ namespace netgen
 
 
 		wrongels = 0;
-		for (int i = 1; i <= mesh.GetNE(); i++)
+                for (ElementIndex ei : mesh.VolumeElements().Range())
 		  {
-		    if (mesh.VolumeElement(i).Volume(mesh.Points()) < 0)
+		    if (mesh.VolumeElement(ei).Volume(mesh.Points()) < 0)
 		      {
 			wrongels++;
-			mesh.VolumeElement(i).Flags().badel = 1;
+			mesh.VolumeElement(ei).Flags().badel = 1;
 			(*testout) << "wrong el: ";
 			for (int j = 1; j <= 4; j++)
-			  (*testout) << mesh.VolumeElement(i).PNum(j) << " ";
+			  (*testout) << mesh.VolumeElement(ei).PNum(j) << " ";
 			(*testout) << endl;
 		      }
 		    else
-		      mesh.VolumeElement(i).Flags().badel = 0;
+		      mesh.VolumeElement(ei).Flags().badel = 0;
 		  }
 		cout << "wrongels = " << wrongels << endl;
 	      }

libsrc/meshing/topology.hpp

@@ -185,6 +185,13 @@ public:
     elnr2 = surf2volelement.Get(selnr)[1];
   }
 
+  std::array<ElementIndex,2> GetSurface2VolumeElement (SurfaceElementIndex sei) 
+  {
+    return { ElementIndex( surf2volelement.Get(sei+1)[0] - 1),
+             ElementIndex( surf2volelement.Get(sei+1)[1] - 1) };
+  }
+
+  
   int GetFace2SurfaceElement (int fnr) const { return face2surfel[fnr-1]; }
 
   SegmentIndex GetSegmentOfEdge(int edgenr) const { return edge2segment[edgenr-1]; }

libsrc/occ/occmeshsurf.cpp

@@ -366,8 +366,11 @@ namespace netgen
 
   void OCCSurface :: Project (Point<3> & ap, PointGeomInfo & gi)
   {
-    // static Timer t("OccSurface::Project"); RegionTimer reg(t);
+    static Timer t("OccSurface::Project"); RegionTimer reg(t);
-    // static Timer t2("OccSurface::Project actual"); 
+    static Timer tanal("OccSurface::Project analysis"); 
+    static Timer ttol("OccSurface::Project approximation"); 
+
+    static Timer t2("OccSurface::Project actual"); 
 
 
     // try Newton's method ...
@@ -473,13 +476,17 @@ namespace netgen
 
     // double u,v;
     // JS : shouldn't we move these 2 lines to the constructor ?
+    // tanal.Start();
     Handle( ShapeAnalysis_Surface ) su = new ShapeAnalysis_Surface( occface );
+    // ShapeAnalysis_Surface su( occface );    
+    // tanal.Stop();
+    ttol.Start();
     auto toltool =  BRep_Tool::Tolerance( topods_face );
-
+    ttol.Stop();
     // gp_Pnt2d suval = su->ValueOfUV ( pnt, toltool);
-    // t2.Start();
+    t2.Start();
     gp_Pnt2d suval = su->NextValueOfUV (gp_Pnt2d(u,v), pnt, toltool);
-    // t2.Stop();
+    t2.Stop();
     suval.Coord( u, v);
     pnt = occface->Value( u, v );
     

libsrc/occ/occmeshsurf.hpp

@@ -5,10 +5,13 @@
 
 #include "occgeom.hpp"
 #include "mydefs.hpp"
-
+#include <BRep_Tool.hxx>
 #include <TopoDS_Face.hxx>
 #include <Geom_Surface.hxx>
 #include <ShapeAnalysis.hxx>
+#include <ShapeAnalysis_Surface.hxx>
+#include <GeomLProp_SLProps.hxx>
+
 
 #define PARAMETERSPACE -1
 #define PLANESPACE     1
@@ -30,7 +33,8 @@ public:
   Handle(Geom_Surface) occface;
   TopAbs_Orientation orient;
   int projecttype;
-
+  ShapeAnalysis_Surface su;
+  Standard_Real toltool;
 protected:
   Point<3> p1;
   Point<3> p2;
@@ -60,10 +64,15 @@ protected:
 
 public:
   OCCSurface (const TopoDS_Face & aface, int aprojecttype)
+    : topods_face(aface),
+      occface(BRep_Tool::Surface(topods_face)),
+      su( occface ),
+      toltool(BRep_Tool::Tolerance(topods_face))
+      
   {
     static Timer t("occurface ctor"); RegionTimer r(t);
     topods_face = aface;
-    occface = BRep_Tool::Surface(topods_face);
+    // occface = BRep_Tool::Surface(topods_face);
     orient = topods_face.Orientation();
     projecttype = aprojecttype;
     ShapeAnalysis::GetFaceUVBounds (topods_face, umin, umax, vmin, vmax);
@@ -72,6 +81,8 @@ public:
     umax += fabs(umax-umin)/100.0;
     vmax += fabs(vmax-vmin)/100.0;
     // projecttype = PLANESPACE;
+
+    // su = ShapeAnalysis_Surface( occface );        
     /*
     TopExp_Explorer exp1;
     exp1.Init (topods_face, TopAbs_WIRE);

libsrc/occ/occpkg.cpp

@@ -817,20 +817,28 @@ namespace netgen
 
      if(strcmp(argv[1], "showall") == 0)
      {
+       /*
         for(int i = 1; i <= mesh->GetNSE(); i++)
         {
-           mesh->SurfaceElement(i).Visible(1);
+          mesh->SurfaceElement(i).Visible(1);
         }
+       */
+       for (auto & el : mesh->SurfaceElements())
+         el.Visible(1);
        
-        mesh->SetNextTimeStamp();
+       mesh->SetNextTimeStamp();
      }
 
      if(strcmp(argv[1], "hideall") == 0)
      {
+       /*
         for(int i = 1; i <= mesh->GetNSE(); i++)
         {
-           mesh->SurfaceElement(i).Visible(0);
+          mesh->SurfaceElement(i).Visible(0);
         }
+       */
+       for (auto & el : mesh->SurfaceElements())
+         el.Visible(0);
        
         mesh->SetNextTimeStamp();
      }

nglib/nglib.cpp

@@ -270,7 +270,7 @@ namespace nglib
    NGLIB_API Ng_Surface_Element_Type 
       Ng_GetSurfaceElement (Ng_Mesh * mesh, int num, int * pi)
    {
-      const Element2d & el = ((Mesh*)mesh)->SurfaceElement(num);
+     const Element2d & el = ((Mesh*)mesh)->SurfaceElement(SurfaceElementIndex(num-1));
       for (int i = 1; i <= el.GetNP(); i++)
          pi[i-1] = el.PNum(i);
       Ng_Surface_Element_Type et;
@@ -301,7 +301,7 @@ namespace nglib
    NGLIB_API Ng_Volume_Element_Type
       Ng_GetVolumeElement (Ng_Mesh * mesh, int num, int * pi)
    {
-      const Element & el = ((Mesh*)mesh)->VolumeElement(num);
+     const Element & el = ((Mesh*)mesh)->VolumeElement(ElementIndex(num-1));
       for (int i = 1; i <= el.GetNP(); i++)
          pi[i-1] = el.PNum(i);
       Ng_Volume_Element_Type et;
@@ -439,7 +439,7 @@ namespace nglib
    NGLIB_API Ng_Surface_Element_Type
       Ng_GetElement_2D (Ng_Mesh * mesh, int num, int * pi, int * matnum)
    {
-      const Element2d & el = ((Mesh*)mesh)->SurfaceElement(num);
+     const Element2d & el = ((Mesh*)mesh)->SurfaceElement(SurfaceElementIndex(num-1));
       for (int i = 1; i <= el.GetNP(); i++)
          pi[i-1] = el.PNum(i);