quicksort for CSG - special point sorting

libsrc/csg/edgeflw.cpp

@@ -33,10 +33,13 @@ namespace netgen
 
   void EdgeCalculation :: Calc(double h, Mesh & mesh)
   {
+    static int timer = NgProfiler::CreateTimer ("CSG: mesh edges");
+    NgProfiler::RegionTimer reg (timer);
+
+
     PrintMessage (1, "Find edges");
     PushStatus ("Find edges");
 
-
     for (int i = 1; i <= mesh.GetNP(); i++)
       meshpoint_tree->Insert (mesh.Point(i), i);
 

libsrc/csg/specpoin.cpp

@@ -1365,6 +1365,10 @@ namespace netgen
 			Array<MeshPoint> & apoints, 
 			Array<SpecialPoint> & specpoints)
   {
+    static int timer = NgProfiler::CreateTimer ("CSG: analyze special points");
+    NgProfiler::RegionTimer reg (timer);
+
+
     Array<int> surfind, rep_surfind, surfind2, rep_surfind2, surfind3;
 
     Array<Vec<3> > normalvecs;
@@ -1381,15 +1385,28 @@ namespace netgen
 
     if (!apoints.Size()) return;
 
-#define VERTSORT
 
-#ifdef VERTSORT
+    {
+      /*
+	sort points in the (arbitrary) direction dir
+	important for periodic boundaries: 
+	corner points on the left and the right boundary come in the same ordering
+      */
       Vec<3> dir(1.2, 1.7, 0.9);
+      
+      Array<double> coord(apoints.Size());
+      for (int i = 0; i < apoints.Size(); i++)
+	coord[i] = dir * Vec<3> (apoints[i]);
+      
+      QuickSort (coord, apoints);
+      
+      /*
 	for (int i = 0; i < apoints.Size(); i++)
 	for (int j = 0; j < apoints.Size()-1; j++)
 	if ( (dir * Vec<3> (apoints[j])) > (dir * Vec<3> (apoints[j+1])))
 	swap (apoints[j], apoints[j+1]);
-#endif
+      */
+    }
 
 
 

libsrc/gprim/adtree.cpp

@@ -2065,13 +2065,13 @@ namespace netgen
 
 
 
-  Point3dTree :: Point3dTree (const Point3d & pmin, const Point3d & pmax)
+  Point3dTree :: Point3dTree (const Point<3> & pmin, const Point<3> & pmax)
   {
     float pmi[3], pma[3];
     for (int i = 0; i < 3; i++)
       {
-	pmi[i] = pmin.X(i+1);
+	pmi[i] = pmin(i);
-	pma[i] = pmax.X(i+1);
+	pma[i] = pmax(i);
       }
     tree = new ADTree3 (pmi, pma);
   }
@@ -2083,23 +2083,23 @@ namespace netgen
 
 
 
-  void Point3dTree :: Insert (const Point3d & p, int pi)
+  void Point3dTree :: Insert (const Point<3> & p, int pi)
   {
-    static float pd[3];
+    float pd[3];
-    pd[0] = p.X();
+    pd[0] = p(0);
-    pd[1] = p.Y();
+    pd[1] = p(1);
-    pd[2] = p.Z();
+    pd[2] = p(2);
     tree->Insert (pd, pi);
   }
 
-  void Point3dTree :: GetIntersecting (const Point3d & pmin, const Point3d & pmax, 
+  void Point3dTree :: GetIntersecting (const Point<3> & pmin, const Point<3> & pmax, 
 				       Array<int> & pis) const
   {
     float pmi[3], pma[3];
     for (int i = 0; i < 3; i++)
       {
-	pmi[i] = pmin.X(i+1);
+	pmi[i] = pmin(i);
-	pma[i] = pmax.X(i+1);
+	pma[i] = pmax(i);
       }
     tree->GetIntersecting (pmi, pma, pis);
   }
@@ -2113,15 +2113,15 @@ namespace netgen
 
 
 
-  Box3dTree :: Box3dTree (const Point3d & apmin, const Point3d & apmax)
+  Box3dTree :: Box3dTree (const Point<3> & apmin, const Point<3> & apmax)
   {
     boxpmin = apmin;
     boxpmax = apmax;
     float tpmin[6], tpmax[6];
     for (int i = 0; i < 3; i++)
       {
-	tpmin[i] = tpmin[i+3] = boxpmin.X(i+1);
+	tpmin[i] = tpmin[i+3] = boxpmin(i);
-	tpmax[i] = tpmax[i+3] = boxpmax.X(i+1);
+	tpmax[i] = tpmax[i+3] = boxpmax(i);
       }
     tree = new ADTree6 (tpmin, tpmax);
   }
@@ -2131,20 +2131,20 @@ namespace netgen
     delete tree;
   }
 
-  void Box3dTree :: Insert (const Point3d & bmin, const Point3d & bmax, int pi)
+  void Box3dTree :: Insert (const Point<3> & bmin, const Point<3> & bmax, int pi)
   {
-    static float tp[6];
+    float tp[6];
 
     for (int i = 0; i < 3; i++)
       {
-	tp[i] = bmin.X(i+1);
+	tp[i] = bmin(i);
-	tp[i+3] = bmax.X(i+1);
+	tp[i+3] = bmax(i);
       }
 
     tree->Insert (tp, pi);
   }
 
-  void Box3dTree ::GetIntersecting (const Point3d & pmin, const Point3d & pmax, 
+  void Box3dTree ::GetIntersecting (const Point<3> & pmin, const Point<3> & pmax, 
 				    Array<int> & pis) const
   {
     float tpmin[6];
@@ -2152,11 +2152,11 @@ namespace netgen
 
     for (int i = 0; i < 3; i++)
       {
-	tpmin[i] = boxpmin.X(i+1);
+	tpmin[i] = boxpmin(i);
-	tpmax[i] = pmax.X(i+1);
+	tpmax[i] = pmax(i);
       
-	tpmin[i+3] = pmin.X(i+1);
+	tpmin[i+3] = pmin(i);
-	tpmax[i+3] = boxpmax.X(i+1);
+	tpmax[i+3] = boxpmax(i);
       }
 
     tree->GetIntersecting (tpmin, tpmax, pis);

libsrc/gprim/adtree.hpp

@@ -447,12 +447,12 @@ class Point3dTree
   ADTree3 * tree;
 
 public:
-  Point3dTree (const Point3d & pmin, const Point3d & pmax);
+  Point3dTree (const Point<3> & pmin, const Point<3> & pmax);
   ~Point3dTree ();
-  void Insert (const Point3d & p, int pi);
+  void Insert (const Point<3> & p, int pi);
   void DeleteElement (int pi) 
     { tree->DeleteElement(pi); }
-  void GetIntersecting (const Point3d & pmin, const Point3d & pmax, 
+  void GetIntersecting (const Point<3> & pmin, const Point<3> & pmax, 
 			Array<int> & pis) const;
   const ADTree3 & Tree() const { return *tree; };
 };
@@ -462,18 +462,18 @@ public:
 class Box3dTree
 {
   ADTree6 * tree;
-  Point3d boxpmin, boxpmax;
+  Point<3> boxpmin, boxpmax;
 public:
-  Box3dTree (const Point3d & apmin, const Point3d & apmax);
+  Box3dTree (const Point<3> & apmin, const Point<3> & apmax);
   ~Box3dTree ();
-  void Insert (const Point3d & bmin, const Point3d & bmax, int pi);
+  void Insert (const Point<3> & bmin, const Point<3> & bmax, int pi);
   void Insert (const Box<3> & box, int pi)
   {
     Insert (box.PMin(), box.PMax(), pi);
   }
   void DeleteElement (int pi) 
     { tree->DeleteElement(pi); }
-  void GetIntersecting (const Point3d & pmin, const Point3d & pmax, 
+  void GetIntersecting (const Point<3> & pmin, const Point<3> & pmax, 
 			Array<int> & pis) const;
 
   const ADTree6 & Tree() const { return *tree; };

libsrc/include/nginterface_v2.hpp

@@ -138,7 +138,7 @@ namespace netgen
   DLL_HEADER int Ng_GetElementIndex (int nr);
 
   
-  
+  /*  
   /// Curved Elements:
   /// xi..... DIM_EL local coordinates
   /// sxi ... step xi
@@ -149,7 +149,7 @@ namespace netgen
                                                  const double * xi, int sxi,
                                                  double * x, int sx,
                                                  double * dxdxi, int sdxdxi);
-  
+  */
   
 }
 #endif

libsrc/occ/occgeom.cpp

@@ -1228,16 +1228,16 @@ namespace netgen
 
 
 
-   char * shapesname[] =
+  const char * shapesname[] =
    {" ", "CompSolids", "Solids", "Shells",
 
    "Faces", "Wires", "Edges", "Vertices"};
 
-   char * shapename[] =
+  const char * shapename[] =
    {" ", "CompSolid", "Solid", "Shell",
    "Face", "Wire", "Edge", "Vertex"};
 
-   char * orientationstring[] =
+  const char * orientationstring[] =
      {"+", "-"};
 
 
@@ -1253,7 +1253,7 @@ namespace netgen
 
       TopExp_Explorer e;
       int count = 0;
-      int count2;
+      int count2 = 0;
 
       if (isfree)
          e.Init(sh, l, TopAbs_ShapeEnum(l-1));
@@ -1282,6 +1282,8 @@ namespace netgen
 	     count2 = emap.FindIndex(TopoDS::Edge(e.Current())); break;
 	   case TopAbs_VERTEX:
 	     count2 = vmap.FindIndex(TopoDS::Vertex(e.Current())); break;
+	   default:
+	     cout << "RecursiveTopologyTree: Case " << e.Current().ShapeType() << " not handeled" << endl;
          }
 
          int nrsubshapes = 0;
@@ -1349,7 +1351,7 @@ namespace netgen
       int stretchedpinfaces = 0;
       int smoothpinfaces = 0;
       int twistedfaces = 0;
-      int edgessamebutnotidentified = 0;
+      // int edgessamebutnotidentified = 0;
 
       cout << "checking faces ... " << flush;
 
@@ -1453,8 +1455,8 @@ namespace netgen
       cout << "done" << endl;
       cout << "checking edges ... " << flush;
 
-      double dmax;
+      // double dmax;
-      int cnt = 0;
+      // int cnt = 0;
       Array <double> edgeLengths;
       Array <int> order;
       edgeLengths.SetSize (emap.Extent());