Merge branch 'master' into stl_to_python

libsrc/meshing/meshclass.cpp

@@ -2148,7 +2148,7 @@ namespace netgen
                 INDEX_2 seg (el.PNumMod(j), el.PNumMod(j+1));
                 INDEX_2 data;
 
-                if (seg.I1() <= 0 || seg.I2() <= 0)
+                if (seg.I1() < PointIndex::BASE || seg.I2() < PointIndex::BASE)
                   cerr << "seg = " << seg << endl;
 
                 if (faceht.Used(seg))
@@ -2355,34 +2355,33 @@ namespace netgen
 
   void Mesh :: RemoveOneLayerSurfaceElements ()
   {
-    int i, j;
     int np = GetNP();
 
     FindOpenSegments();
-    BitArray frontpoints(np);
+    BitArray frontpoints(np+1);  // for 0- and 1-based
-
     frontpoints.Clear();
-    for (i = 1; i <= GetNOpenSegments(); i++)
+    
+    for (int i = 1; i <= GetNOpenSegments(); i++)
       {
         const Segment & seg = GetOpenSegment(i);
         frontpoints.Set (seg[0]);
         frontpoints.Set (seg[1]);
       }
 
-    for (i = 1; i <= GetNSE(); i++)
+    for (int i = 1; i <= GetNSE(); i++)
       {
         Element2d & sel = surfelements.Elem(i);
-        int remove = 0;
+        bool remove = false;
-        for (j = 1; j <= sel.GetNP(); j++)
+        for (int j = 1; j <= sel.GetNP(); j++)
           if (frontpoints.Test(sel.PNum(j)))
-            remove = 1;
+            remove = true;
         if (remove)
-          sel.PNum(1) = 0;
+          sel.PNum(1).Invalidate();
       }
 
-    for (i = surfelements.Size(); i >= 1; i--)
+    for (int i = surfelements.Size(); i >= 1; i--)
       {
-        if (surfelements.Elem(i).PNum(1) == 0)
+        if (!surfelements.Elem(i).PNum(1).IsValid())
           {
             surfelements.Elem(i) = surfelements.Last();
             surfelements.DeleteLast();

libsrc/meshing/meshclass.hpp

@@ -170,7 +170,7 @@ namespace netgen
     /// number of refinement levels
     int mglevels;
     /// refinement hierarchy
-    Array<INDEX_2,PointIndex::BASE> mlbetweennodes;
+    Array<PointIndices<2>,PointIndex::BASE> mlbetweennodes;
     /// parent element of volume element
     Array<int> mlparentelement;
     /// parent element of surface element

libsrc/meshing/meshtype.hpp

@@ -144,7 +144,18 @@ namespace netgen
     return (ost << int(pi));
   }
 
-
+  template <int N> class PointIndices;
+  template <> class PointIndices<2> : public INDEX_2
+  {
+  public:
+    PointIndices () = default;
+    PointIndices (INDEX_2 i2) : INDEX_2(i2) { ; }
+    PointIndices (PointIndex i1, PointIndex i2) : INDEX_2(i1,i2) { ; } 
+    PointIndex operator[] (int i) const { return PointIndex(INDEX_2::operator[](i)); }
+    PointIndex & operator[] (int i) { return reinterpret_cast<PointIndex&>(INDEX_2::operator[](i)); }
+    static PointIndices Sort(PointIndex i1, PointIndex i2) { return INDEX_2::Sort(i1, i2); } 
+  };
+  
 
 
   class ElementIndex

libsrc/meshing/refine.cpp

@@ -54,14 +54,13 @@ namespace netgen
 	  case TRIG:
 	  case TRIG6:
 	    {
-	      static int betw[3][3] =
+              static int betw[3][3] =
-		{ { 2, 3, 4 },
+		{ { 1, 2, 3 },
-		  { 1, 3, 5 },
+		  { 0, 2, 4 },
-		  { 1, 2, 6 } };
+		  { 0, 1, 5 } };
-
               for (int j = 0; j < 3; j++)
                 {
-                  INDEX_2 i2 = INDEX_2::Sort(el.PNum(betw[j][0]),el.PNum(betw[j][1]));
+                  auto i2 = PointIndices<2>::Sort(el[betw[j][0]],el[betw[j][1]]);
                   if (!between.Used(i2))
                     {
                       between.Set (i2, 0);          
@@ -771,12 +770,9 @@ namespace netgen
 
 	BitArray boundp(np);
 	boundp.Clear();
-	for (int i = 1; i <= mesh.GetNSE(); i++)
+	for (auto & sel : mesh.SurfaceElements())
-	  {
+          for (auto pi : sel.PNums())
-	    const Element2d & sel = mesh.SurfaceElement(i);
+            boundp.Set(pi);
-	    for (int j = 1; j <= sel.GetNP(); j++)
-	      boundp.Set(sel.PNum(j));
-	  }
 
 
 	double lam = 0.5;
@@ -835,11 +831,9 @@ namespace netgen
 		    (*testout) << "p " << i << endl;
 
 		(*testout) << "surf points: " << endl;
-		for (int i = 1; i <= mesh.GetNSE(); i++)
+		for (auto & sel : mesh.SurfaceElements())
-		  for (int j = 1; j <= 3; j++)
+		  for (auto pi : sel.PNums())
-		    (*testout) << mesh.SurfaceElement(i).PNum(j) << endl;
+		    (*testout) << pi << endl;
-		  
-
 
 		mesh.CalcSurfacesOfNode();
 		free.Invert();

libsrc/stlgeom/meshstlsurface.cpp

@@ -121,8 +121,8 @@ static void STLFindEdges (STLGeometry & geom,
 	  */
 	  Point3d hp, hp2;
 	  Segment seg;
-	  seg[0] = p1;
+	  seg[0] = p1 + PointIndex::BASE-1;
-	  seg[1] = p2;
+	  seg[1] = p2 + PointIndex::BASE-1;
 	  seg.si = geom.GetTriangle(trig1).GetFaceNum();
 	  seg.edgenr = i;
 
@@ -177,8 +177,8 @@ static void STLFindEdges (STLGeometry & geom,
 
 
 	  Segment seg2;
-	  seg2[0] = p2;
+	  seg2[0] = p2 + PointIndex::BASE-1;;
-	  seg2[1] = p1;
+	  seg2[1] = p1 + PointIndex::BASE-1;;
 	  seg2.si = geom.GetTriangle(trig2).GetFaceNum();
 	  seg2.edgenr = i;