improve STL makeatlas: searchtree, templetize searchtree

libsrc/csg/genmesh.cpp

@@ -127,7 +127,7 @@ namespace netgen
     
     Point3d pmin, pmax;
     mesh.GetBox (pmin, pmax);
-    Box3dTree segtree (pmin, pmax);
+    BoxTree<3> segtree (pmin, pmax);
     
     for (SegmentIndex si = 0; si < mesh.GetNSeg(); si++)
       {

libsrc/gprim/adtree.cpp

@@ -1794,6 +1794,243 @@ namespace netgen
 
 
 
+  
+  template <int dim>
+  T_ADTree<dim> :: T_ADTree (const float * acmin, 
+                             const float * acmax)
+    : ela(0)
+  {
+    memcpy (cmin, acmin, dim * sizeof(float));
+    memcpy (cmax, acmax, dim * sizeof(float));
+    
+    root = new T_ADTreeNode<dim>;
+    root->sep = (cmin[0] + cmax[0]) / 2;
+  }
+
+  template <int dim>  
+  T_ADTree<dim> :: ~T_ADTree ()
+  {
+    root->DeleteChilds();
+    delete root;
+  }
+
+  template <int dim>
+  void T_ADTree<dim> :: Insert (const float * p, int pi)
+  {
+    T_ADTreeNode<dim> *node(NULL);
+    T_ADTreeNode<dim> *next;
+    int dir;
+    int lr(0);
+
+    float bmin[dim];
+    float bmax[dim];
+  
+    memcpy (bmin, cmin, dim * sizeof(float));
+    memcpy (bmax, cmax, dim * sizeof(float));
+
+    next = root;
+    dir = 0;
+    while (next)
+      {
+	node = next;
+
+	if (node->pi == -1)
+	  {    
+	    memcpy (node->data, p, dim * sizeof(float));
+	    node->pi = pi;
+
+	    if (ela.Size() < pi+1)
+	      ela.SetSize (pi+1);
+	    ela[pi] = node;
+
+	    return;
+	  }
+
+	if (node->sep > p[dir])
+	  {
+	    next = node->left;
+	    bmax[dir] = node->sep;
+	    lr = 0;
+	  }
+	else
+	  {
+	    next = node->right;
+	    bmin[dir] = node->sep;
+	    lr = 1;
+	  }
+
+	dir++;
+	if (dir == dim) dir = 0;
+      }
+
+
+    next = new T_ADTreeNode<dim>;
+    memcpy (next->data, p, dim * sizeof(float));
+    next->pi = pi;
+    next->sep = (bmin[dir] + bmax[dir]) / 2;
+
+    if (ela.Size() < pi+1)
+      ela.SetSize (pi+1);
+    ela[pi] = next;
+
+    if (lr)
+      node->right = next;
+    else
+      node->left = next;
+    next -> father = node;
+
+    while (node)
+      {
+	node->nchilds++;
+	node = node->father;
+      }
+  }
+
+  template <int dim>
+  void T_ADTree<dim> :: DeleteElement (int pi)
+  {
+    T_ADTreeNode<dim> * node = ela[pi];
+
+    node->pi = -1;
+
+    node = node->father;
+    while (node)
+      {
+	node->nchilds--;
+	node = node->father;
+      }
+  }
+
+  template <int dim>
+  void T_ADTree<dim> :: PrintMemInfo (ostream & ost) const
+  {
+    ost << Elements() << " elements a " << sizeof(ADTreeNode6) 
+	<< " Bytes = "
+	<< Elements() * sizeof(T_ADTreeNode<dim>) << endl;
+    ost << "maxind = " << ela.Size() << " = " << sizeof(T_ADTreeNode<dim>*) * ela.Size() << " Bytes" << endl;
+  }
+
+
+  template <int dim>
+  class inttn {
+  public:
+    int dir;
+    T_ADTreeNode<dim> * node;
+  };
+
+
+  template <int dim>
+  void T_ADTree<dim> :: GetIntersecting (const float * bmin, 
+                                         const float * bmax,
+                                         Array<int> & pis) const
+  {
+    // static Array<inttn6> stack(10000);
+    // stack.SetSize (10000);
+    ArrayMem<inttn<dim>,10000> stack(10000);
+    pis.SetSize(0);
+
+    stack[0].node = root;
+    stack[0].dir = 0;
+    int stacks = 0;
+
+    while (stacks >= 0)
+      {
+	T_ADTreeNode<dim> * node = stack[stacks].node;
+	int dir = stack[stacks].dir; 
+
+	stacks--;
+	if (node->pi != -1)
+	  {
+            bool found = true;
+            for (int i = 0; i < dim/2; i++)
+              if (node->data[i] > bmax[i])
+                found = false;
+            for (int i = dim/2; i < dim; i++)
+              if (node->data[i] < bmin[i])
+                found = false;
+            if (found)
+              pis.Append (node->pi);            
+            /*
+	    if (node->data[0] > bmax[0] || 
+		node->data[1] > bmax[1] || 
+		node->data[2] > bmax[2] || 
+		node->data[3] < bmin[3] || 
+		node->data[4] < bmin[4] || 
+		node->data[5] < bmin[5])
+	      ;
+	    else
+              {
+                pis.Append (node->pi);
+              }
+            */
+	  }
+
+	int ndir = (dir+1) % dim;
+
+	if (node->left && bmin[dir] <= node->sep)
+	  {
+	    stacks++;
+	    stack[stacks].node = node->left;
+	    stack[stacks].dir = ndir;
+	  }
+	if (node->right && bmax[dir] >= node->sep)
+	  {
+	    stacks++;
+	    stack[stacks].node = node->right;
+	    stack[stacks].dir = ndir;
+	  }
+      }
+  }
+
+  template <int dim>
+  void T_ADTree<dim> :: PrintRec (ostream & ost, const T_ADTreeNode<dim> * node) const
+  {
+    
+    // if (node->data)     // true anyway
+      {
+	ost << node->pi << ": ";
+	ost << node->nchilds << " childs, ";
+	for (int i = 0; i < dim; i++)
+	  ost << node->data[i] << " ";
+	ost << endl;
+      }
+    if (node->left)
+      PrintRec (ost, node->left);
+    if (node->right)
+      PrintRec (ost, node->right);
+  }
+
+  template <int dim>
+  int T_ADTree<dim> :: DepthRec (const T_ADTreeNode<dim> * node) const
+  {
+    int ldepth = 0;
+    int rdepth = 0;
+
+    if (node->left)
+      ldepth = DepthRec(node->left);
+    if (node->right)
+      rdepth = DepthRec(node->right);
+    return 1 + max2 (ldepth, rdepth);
+  }
+
+  template <int dim>
+  int T_ADTree<dim> :: ElementsRec (const T_ADTreeNode<dim> * node) const
+  {
+    int els = 1;
+    if (node->left)
+      els += ElementsRec(node->left);
+    if (node->right)
+      els += ElementsRec(node->right);
+    return els;
+  }
+
+
+
+
+  
+
+
+
 #ifdef ABC
 
   /* ******************************* ADTree6F ******************************* */
@@ -2112,67 +2349,80 @@ namespace netgen
 
 
 
-
+  template <int dim>
-  Box3dTree :: Box3dTree (const Box<3> & abox)
+  BoxTree<dim> :: BoxTree (const Box<dim> & abox)
   {
     boxpmin = abox.PMin();
     boxpmax = abox.PMax();
-    float tpmin[6], tpmax[6];
+    float tpmin[2*dim], tpmax[2*dim];
-    for (int i = 0; i < 3; i++)
+    for (int i = 0; i < dim; i++)
       {
-	tpmin[i] = tpmin[i+3] = boxpmin(i);
+	tpmin[i] = tpmin[i+dim] = boxpmin(i);
-	tpmax[i] = tpmax[i+3] = boxpmax(i);
+	tpmax[i] = tpmax[i+dim] = boxpmax(i);
       }
-    tree = new ADTree6 (tpmin, tpmax);
+    tree = new T_ADTree<2*dim> (tpmin, tpmax);
   }
 
-  Box3dTree :: Box3dTree (const Point<3> & apmin, const Point<3> & apmax)
+  template <int dim>
+  BoxTree<dim> :: BoxTree (const Point<dim> & apmin, const Point<dim> & apmax)
   {
     boxpmin = apmin;
     boxpmax = apmax;
-    float tpmin[6], tpmax[6];
+    float tpmin[2*dim], tpmax[2*dim];
-    for (int i = 0; i < 3; i++)
+    for (int i = 0; i < dim; i++)
       {
-	tpmin[i] = tpmin[i+3] = boxpmin(i);
+	tpmin[i] = tpmin[i+dim] = boxpmin(i);
-	tpmax[i] = tpmax[i+3] = boxpmax(i);
+	tpmax[i] = tpmax[i+dim] = boxpmax(i);
       }
-    tree = new ADTree6 (tpmin, tpmax);
+    tree = new T_ADTree<2*dim> (tpmin, tpmax);
   }
 
-  Box3dTree :: ~Box3dTree ()
+  template <int dim>
+  BoxTree<dim> :: ~BoxTree ()
   {
     delete tree;
   }
 
-  void Box3dTree :: Insert (const Point<3> & bmin, const Point<3> & bmax, int pi)
+  template <int dim>
+  void BoxTree<dim> :: Insert (const Point<dim> & bmin, const Point<dim> & bmax, int pi)
   {
-    float tp[6];
+    float tp[2*dim];
 
-    for (int i = 0; i < 3; i++)
+    for (int i = 0; i < dim; i++)
       {
 	tp[i] = bmin(i);
-	tp[i+3] = bmax(i);
+	tp[i+dim] = bmax(i);
       }
 
     tree->Insert (tp, pi);
   }
 
-  void Box3dTree ::GetIntersecting (const Point<3> & pmin, const Point<3> & pmax, 
+  template <int dim>
+  void BoxTree<dim> ::GetIntersecting (const Point<dim> & pmin, const Point<dim> & pmax, 
                                        Array<int> & pis) const
   {
-    float tpmin[6];
+    float tpmin[2*dim];
-    float tpmax[6];
+    float tpmax[2*dim];
 
-    for (int i = 0; i < 3; i++)
+    for (int i = 0; i < dim; i++)
       {
 	tpmin[i] = boxpmin(i);
 	tpmax[i] = pmax(i);
         
-	tpmin[i+3] = pmin(i);
+	tpmin[i+dim] = pmin(i);
-	tpmax[i+3] = boxpmax(i);
+	tpmax[i+dim] = boxpmax(i);
       }
 
     tree->GetIntersecting (tpmin, tpmax, pis);
   }
 
+
+  template<> BlockAllocator T_ADTreeNode<4> :: ball(sizeof (T_ADTreeNode<4>));
+  template class T_ADTree<4>;
+  template class BoxTree<2>;
+
+  
+  template<> BlockAllocator T_ADTreeNode<6> :: ball(sizeof (T_ADTreeNode<6>));
+  template class T_ADTree<6>;
+  template class BoxTree<3>;
 }

libsrc/gprim/adtree.hpp

@@ -382,6 +382,97 @@ public:
 
 
 
+
+template <int DIM>
+class T_ADTreeNode
+{
+public:
+  T_ADTreeNode *left, *right, *father;
+  float sep;
+  float data[DIM];
+  int pi;
+  int nchilds;
+
+  T_ADTreeNode ()
+  {
+    pi = -1;
+    left = NULL;
+    right = NULL;
+    father = NULL;
+    nchilds = 0;
+  }
+  
+  void DeleteChilds ()
+  {
+    if (left)
+      {
+	left->DeleteChilds();
+	delete left;
+	left = NULL;
+      }
+    if (right)
+      {
+	right->DeleteChilds();
+	delete right;
+	right = NULL;
+      }
+  }
+  
+  // friend class T_ADTree<DIM>;
+
+  static BlockAllocator ball;
+  void * operator new(size_t)
+  {
+    return ball.Alloc();
+  }
+  
+  void operator delete (void * p)
+  {
+    ball.Free(p);
+  }
+};
+
+
+
+
+  template <int dim>
+  class T_ADTree
+  {
+    T_ADTreeNode<dim> * root;
+    float cmin[dim], cmax[dim];
+    Array<T_ADTreeNode<dim>*> ela;
+
+  public:
+    T_ADTree (const float * acmin, 
+             const float * acmax);
+    ~T_ADTree ();
+    
+    void Insert (const float * p, int pi);
+    void GetIntersecting (const float * bmin, const float * bmax,
+                          Array<int> & pis) const;
+    
+    void DeleteElement (int pi);
+    
+    
+    void Print (ostream & ost) const
+    { PrintRec (ost, root); }
+    int Depth () const
+    { return DepthRec (root); }
+    int Elements () const
+    { return ElementsRec (root); }
+    
+    void PrintRec (ostream & ost, const T_ADTreeNode<dim> * node) const;
+    int DepthRec (const T_ADTreeNode<dim> * node) const;
+    int ElementsRec (const T_ADTreeNode<dim> * node) const;
+    
+    void PrintMemInfo (ostream & ost) const;
+  };
+  
+  
+
+
+
+
 /*
 
 class ADTreeNode6F
@@ -460,26 +551,26 @@ public:
 };
 
 
-
+  template <int dim>
-class Box3dTree
+  class BoxTree
   {
-  ADTree6 * tree;
+    T_ADTree<2*dim> * tree;
-  Point<3> boxpmin, boxpmax;
+    Point<dim> boxpmin, boxpmax;
   public:
-  Box3dTree (const Box<3> & abox);
+    BoxTree (const Box<dim> & abox);
-  Box3dTree (const Point<3> & apmin, const Point<3> & apmax);
+    BoxTree (const Point<dim> & apmin, const Point<dim> & apmax);
-  ~Box3dTree ();
+    ~BoxTree ();
-  void Insert (const Point<3> & bmin, const Point<3> & bmax, int pi);
+    void Insert (const Point<dim> & bmin, const Point<dim> & bmax, int pi);
-  void Insert (const Box<3> & box, int pi)
+    void Insert (const Box<dim> & box, int pi)
     {
       Insert (box.PMin(), box.PMax(), pi);
     }
     void DeleteElement (int pi) 
     { tree->DeleteElement(pi); }
-  void GetIntersecting (const Point<3> & pmin, const Point<3> & pmax, 
+    void GetIntersecting (const Point<dim> & pmin, const Point<dim> & pmax, 
                           Array<int> & pis) const;
     
-  const ADTree6 & Tree() const { return *tree; };
+    const T_ADTree<2*dim> & Tree() const { return *tree; };
 };
 
 }

libsrc/gprim/geomobjects.hpp

@@ -232,8 +232,13 @@ namespace netgen
     enum EB_TYPE { EMPTY_BOX = 1 };
     Box ( EB_TYPE et ) 
     {
-      pmin = Point<3> (1e99, 1e99, 1e99);
+      for (int i = 0; i < D; i++)
-      pmax = Point<3> (-1e99, -1e99, -1e99);
+        {
+          pmin(i) = 1e99;
+          pmax(i) = -1e99;
+        }
+      // pmin = Point<D> (1e99, 1e99, 1e99);
+      // pmax = Point<D> (-1e99, -1e99, -1e99);
     }
 
     const Point<D> & PMin () const { return pmin; }

libsrc/meshing/adfront2.hpp

@@ -169,7 +169,7 @@ class AdFront2
   Array<FrontLine> lines;     /// front lines
 
   Box3d boundingbox;
-  Box3dTree linesearchtree;       /// search tree for lines
+  BoxTree<3> linesearchtree;       /// search tree for lines
   Point3dTree pointsearchtree;    /// search tree for points
   Point3dTree cpointsearchtree;   /// search tree for cone points (not used ???)
 

libsrc/meshing/adfront3.cpp

@@ -245,7 +245,7 @@ void AdFront3 :: CreateTrees ()
   pmin = pmin + 0.5 * (pmin - pmax);
 
   delete facetree;
-  facetree = new Box3dTree (pmin, pmax);
+  facetree = new BoxTree<3> (pmin, pmax);
   
   for (i = 1; i <= GetNF(); i++)
     {

libsrc/meshing/adfront3.hpp

@@ -212,7 +212,7 @@ class AdFront3
   Array<char> pingroup;
   
   ///
-  class Box3dTree * facetree;
+  class BoxTree<3> * facetree;
 public:
   
   ///

libsrc/meshing/delaunay.cpp

@@ -246,7 +246,7 @@ namespace netgen
   void AddDelaunayPoint (PointIndex newpi, const Point3d & newp, 
 			 Array<DelaunayTet> & tempels, 
 			 Mesh & mesh,
-			 Box3dTree & tettree, 
+			 BoxTree<3> & tettree, 
 			 MeshNB & meshnb,
 			 Array<Point<3> > & centers, Array<double> & radi2,
 			 Array<int> & connected, Array<int> & treesearch, 
@@ -680,7 +680,7 @@ namespace netgen
     pmin2 = pmin2 + 0.1 * (pmin2 - pmax2);
     pmax2 = pmax2 + 0.1 * (pmax2 - pmin2);
 
-    Box3dTree tettree(pmin2, pmax2);
+    BoxTree<3> tettree(pmin2, pmax2);
 
 
     tempels.Append (startel);
@@ -1110,7 +1110,7 @@ namespace netgen
     PrintMessage (3, "Remove intersecting");
     if (openels.Size())
       {
-	Box3dTree setree(pmin, pmax);
+	BoxTree<3> setree(pmin, pmax);
 
 	/*      
 		cout << "open elements in search tree: " << openels.Size() << endl;

libsrc/meshing/delaunay2d.cpp

@@ -285,7 +285,7 @@ namespace netgen
     Point<3> hp = mesh[startel[0]];
     hp(2) = 1;  hbox.Add (hp);
     hp(2) = -1; hbox.Add (hp);
-    Box3dTree searchtree(hbox);
+    BoxTree<3> searchtree(hbox);
 
     Array<DelaunayTrig> tempels;
     startel.CalcCenter (mesh);

libsrc/meshing/meshclass.cpp

@@ -3464,7 +3464,7 @@ namespace netgen
 
     Point3d pmin, pmax;
     GetBox (pmin, pmax);
-    Box3dTree setree(pmin, pmax);
+    BoxTree<3> setree(pmin, pmax);
     Array<int> inters;
 
     bool overlap = 0;
@@ -4282,7 +4282,7 @@ namespace netgen
                     box.Add (points[surfelements[sei].PNums()]);
                   
                   box.Increase (1.01 * box.Diam());
-                  elementsearchtree = new Box3dTree (box);
+                  elementsearchtree = new BoxTree<3> (box);
                   
                   for (SurfaceElementIndex sei = 0; sei < ne; sei++)
                     {
@@ -4297,7 +4297,7 @@ namespace netgen
                     box.Add (points[volelements[ei].PNums()]);
                   
                   box.Increase (1.01 * box.Diam());
-                  elementsearchtree = new Box3dTree (box);
+                  elementsearchtree = new BoxTree<3> (box);
                   
                   for (ElementIndex ei = 0; ei < ne; ei++)
                     {

libsrc/meshing/meshclass.hpp

@@ -99,7 +99,7 @@ namespace netgen
     int numvertices;
 
     /// geometric search tree for interval intersection search
-    Box3dTree * elementsearchtree;
+    BoxTree<3> * elementsearchtree;
     /// time stamp for tree
     mutable int elementsearchtreets;
 

libsrc/meshing/meshing2.cpp

@@ -244,7 +244,7 @@ namespace netgen
 
 
     // test for 3d overlaps
-    Box3dTree surfeltree (boundingbox.PMin(),
+    BoxTree<3> surfeltree (boundingbox.PMin(),
                            boundingbox.PMax());
 
     Array<int> intersecttrias;

libsrc/occ/occgenmesh.cpp

@@ -1181,8 +1181,8 @@ namespace netgen
 
             Array<Line> lines(sections*nedges);
 
-            Box3dTree* searchtree =
+            BoxTree<3> * searchtree =
-               new Box3dTree (bb.PMin(), bb.PMax());
+              new BoxTree<3> (bb.PMin(), bb.PMax());
 
             int nlines = 0;
             for (int i = 1; i <= nedges && !multithread.terminate; i++)

libsrc/stlgeom/stlgeom.cpp

@@ -18,6 +18,7 @@ void STLMeshing (STLGeometry & geom,
   geom.Clear();
   geom.BuildEdges();
   geom.MakeAtlas(mesh);
+  if (multithread.terminate) { return; }
   geom.CalcFaceNums();
   geom.AddFaceEdges();
   geom.LinkEdges();
@@ -48,7 +49,7 @@ void STLMeshing (STLGeometry & geom,
   meshchart = 0; // initialize all ?? JS
 
   if (geomsearchtreeon)
-    searchtree = new Box3dTree (GetBoundingBox().PMin() - Vec3d(1,1,1),
+    searchtree = new BoxTree<3> (GetBoundingBox().PMin() - Vec3d(1,1,1),
                                  GetBoundingBox().PMax() + Vec3d(1,1,1));
   else
     searchtree = NULL;
@@ -2115,7 +2116,7 @@ int STLGeometry :: CheckGeometryOverlapping()
   Point<3> pmin = geombox.PMin();
   Point<3> pmax = geombox.PMax();
 
-  Box3dTree setree(pmin, pmax);
+  BoxTree<3> setree(pmin, pmax);
 
   int oltrigs = 0;
   markedtrigs.SetSize(GetNT());

libsrc/stlgeom/stlgeomchart.cpp

@@ -24,6 +24,10 @@ void STLGeometry :: MakeAtlas(Mesh & mesh)
   int timer3 = NgProfiler::CreateTimer ("makeatlas - part 3");
   int timer4 = NgProfiler::CreateTimer ("makeatlas - part 4");
   int timer5 = NgProfiler::CreateTimer ("makeatlas - part 5");
+  int timer5a = NgProfiler::CreateTimer ("makeatlas - part 5a");
+  int timer5b = NgProfiler::CreateTimer ("makeatlas - part 5b");
+  int timer5cs = NgProfiler::CreateTimer ("makeatlas - part 5cs");
+  int timer5cl = NgProfiler::CreateTimer ("makeatlas - part 5cl");
 
   PushStatusF("Make Atlas");
 
@@ -61,6 +65,7 @@ void STLGeometry :: MakeAtlas(Mesh & mesh)
   Array<int> innerpointstochart(GetNP()); //point in chart becomes chartnum
   Array<int> chartpoints;             //point in chart becomes chartnum
   Array<int> innerchartpoints;
+  Array<Point<3>> innerchartpts;
   Array<int> dirtycharttrigs;
 
   Array<int> chartdistacttrigs (GetNT());   //outercharttrigs
@@ -131,6 +136,7 @@ void STLGeometry :: MakeAtlas(Mesh & mesh)
 
       chartpoints.SetSize(0);  
       innerchartpoints.SetSize(0);
+      innerchartpts.SetSize(0);
       chartbound.Clear();
       chartbound.SetChart(chart);
 
@@ -279,6 +285,11 @@ void STLGeometry :: MakeAtlas(Mesh & mesh)
 	    }
 	}
 
+      innerchartpts.SetSize(innerchartpoints.Size());
+      for (size_t i = 0; i < innerchartpoints.Size(); i++)
+        innerchartpts[i] = GetPoint(innerchartpoints[i]);
+      chartbound.BuildSearchTree();
+      
       NgProfiler::StopTimer (timer2);
       NgProfiler::StartTimer (timer3);
       
@@ -328,7 +339,7 @@ void STLGeometry :: MakeAtlas(Mesh & mesh)
 		    {
 		      accepted = 1;
 		      
-                      NgProfiler::StartTimer (timer4);
+                      // NgProfiler::StartTimer (timer4);
 
 		      bool isdirtytrig = false;
 		      Vec<3> gn = GetTriangle(nt).GeomNormal(points);
@@ -368,9 +379,9 @@ void STLGeometry :: MakeAtlas(Mesh & mesh)
 			    if (!accepted) break;
 			  }
 		      
-                      NgProfiler::StopTimer (timer4);		      
+                      // NgProfiler::StopTimer (timer4);		      
 
-                      NgProfiler::RegionTimer reg5(timer5);		      
+                      // NgProfiler::RegionTimer reg5(timer5);		      
 
                       
 		      // outer chart is only small environment of
@@ -378,6 +389,7 @@ void STLGeometry :: MakeAtlas(Mesh & mesh)
 		      
 		      if (accepted)
 			{
+                          // NgProfiler::StartTimer (timer5a);
 			  accepted = 0;
 			  
 			  for (int k = 1; k <= 3; k++)
@@ -387,15 +399,29 @@ void STLGeometry :: MakeAtlas(Mesh & mesh)
 				break;
 			      }
 
+                          // NgProfiler::StopTimer (timer5a);
+                          // int timer5csl = (innerchartpts.Size() < 100) ? timer5cs : timer5cl;
+                          // NgProfiler::StartTimer (timer5csl);
+			  
 			  if (!accepted)
 			    for (int k = 1; k <= 3; k++)
 			      {
 				Point<3> pt = GetPoint(ntrig.PNum(k));					  
 				double h2 = sqr(mesh.GetH(pt));
+                                /*
                                 for (int l = 1; l <= innerchartpoints.Size(); l++)
 				  {
-				    double tdist = 
+				    double tdist = Dist2(pt, GetPoint (innerchartpoints.Get(l)));
-				      Dist2(pt, GetPoint (innerchartpoints.Get(l)));
+				    if (tdist < 4 * h2)
+				      {
+					accepted = 1; 
+					break;
+				      }
+				  }
+                                */
+                                for (int l = 0; l < innerchartpts.Size(); l++)
+				  {
+				    double tdist = Dist2(pt, innerchartpts[l]);
 				    if (tdist < 4 * h2)
 				      {
 					accepted = 1; 
@@ -404,8 +430,10 @@ void STLGeometry :: MakeAtlas(Mesh & mesh)
 				  }
 				if (accepted) break;
 			      }
-			}
                           
+                          // NgProfiler::StopTimer (timer5csl);
+			}
+                      // NgProfiler::StartTimer (timer5b);
 		      
 		      if (accepted)
 			{
@@ -427,11 +455,14 @@ void STLGeometry :: MakeAtlas(Mesh & mesh)
 				}
 			    }
 			}
+                      // NgProfiler::StopTimer (timer5b);
 		    }	       
 		}
 	    }
 	}            
 
+      chartbound.DeleteSearchTree();
+
       NgProfiler::StopTimer (timer3);
 
       //end of while loop for outer chart

libsrc/stlgeom/stlgeommesh.cpp

@@ -930,7 +930,7 @@ void STLGeometry :: RestrictLocalH(class Mesh & mesh, double gh)
       double mindist = 1E50;
       
       PrintMessage(7,"build search tree...");
-      Box3dTree* lsearchtree = new Box3dTree (GetBoundingBox().PMin() - Vec3d(1,1,1),
+      BoxTree<3> * lsearchtree = new BoxTree<3> (GetBoundingBox().PMin() - Vec3d(1,1,1),
                                                  GetBoundingBox().PMax() + Vec3d(1,1,1));
       
       Array<Point3d> pmins(GetNLines());

libsrc/stlgeom/stltool.cpp

@@ -617,7 +617,7 @@ STLChart :: STLChart(STLGeometry * ageometry)
   geometry = ageometry;
 
   if ( stlparam.usesearchtree == 1)
-    searchtree = new Box3dTree (geometry->GetBoundingBox().PMin() - Vec3d(1,1,1),
+    searchtree = new BoxTree<3> (geometry->GetBoundingBox().PMin() - Vec3d(1,1,1),
                                  geometry->GetBoundingBox().PMax() + Vec3d(1,1,1));
   else
     searchtree = NULL;
@@ -745,7 +745,7 @@ void STLChart :: DelChartTrigs(const Array<int>& trigs)
     {
       PrintMessage(7, "Warning: unsecure routine due to first use of searchtrees!!!");
       //bould new searchtree!!!
-      searchtree = new Box3dTree (geometry->GetBoundingBox().PMin() - Vec3d(1,1,1),
+      searchtree = new BoxTree<3> (geometry->GetBoundingBox().PMin() - Vec3d(1,1,1),
                                    geometry->GetBoundingBox().PMax() + Vec3d(1,1,1));
 
       for (int i = 1; i <= charttrigs->Size(); i++)
@@ -1088,14 +1088,49 @@ int STLBoundary :: TestSeg(const Point<3>& p1, const Point<3> & p2, const Vec<3>
   //  return (maxvalnew < eps);
 }
 
+void STLBoundary :: BuildSearchTree()
+{
+  static int timer = NgProfiler::CreateTimer ("BuildSearchTree");
+  NgProfiler::RegionTimer reg(timer);
   
+  delete searchtree;
+
+  Box<2> box2d(Box<2>::EMPTY_BOX);  
+
+  int nseg = NOSegments();  
+  for (int j = 1; j <= nseg; j++)
+    {
+      const STLBoundarySeg & seg = GetSegment(j);
+      if (seg.IsSmoothEdge()) continue;
+      box2d.Add(seg.BoundingBox().PMin());
+      box2d.Add(seg.BoundingBox().PMax());
+    }
+
+  searchtree = new BoxTree<2> (box2d);
+
+  for (int j = 1; j <= nseg; j++)
+    {
+      const STLBoundarySeg & seg = GetSegment(j);
+      if (seg.IsSmoothEdge()) continue;
+      searchtree -> Insert (seg.BoundingBox(), j);
+    }  
+}
+
+void STLBoundary :: DeleteSearchTree()
+{
+  static int timer = NgProfiler::CreateTimer ("DeleteSearchTree");
+  NgProfiler::RegionTimer reg(timer);
+  
+  delete searchtree;
+  searchtree = nullptr;
+}
 
 // checks, whether 2d projection intersects
 int STLBoundary :: TestSegChartNV(const Point3d & p1, const Point3d& p2, 
 				  const Vec3d& sn)
 {
-  int timer = NgProfiler::CreateTimer ("TestSegChartNV");
+  static int timerquick = NgProfiler::CreateTimer ("TestSegChartNV-searchtree");
-  NgProfiler::StartTimer (timer);
+  static int timer = NgProfiler::CreateTimer ("TestSegChartNV");
 
   int nseg = NOSegments();
  
@@ -1120,6 +1155,40 @@ int STLBoundary :: TestSegChartNV(const Point3d & p1, const Point3d& p2,
   cout << "avg nseg = " << double(totnseg)/cnt << endl;
   */
 
+  if (searchtree)
+    {
+      // NgProfiler::RegionTimer reg(timerquick);      
+      
+      ArrayMem<int,100> pis;
+      searchtree -> GetIntersecting (box2d.PMin(), box2d.PMax(), pis);
+      
+      for (int j : pis)
+        {
+          const STLBoundarySeg & seg = GetSegment(j);
+          
+          if (seg.IsSmoothEdge()) continue;
+          if (!box2d.Intersect (seg.BoundingBox())) continue;
+          
+          const Point<2> & sp1 = seg.P2D1();
+          const Point<2> & sp2 = seg.P2D2();
+          
+          Line2d l2 (sp1, sp2);
+          double lam1, lam2;
+          
+          int err = CrossPointBarycentric (l1, l2, lam1, lam2);
+          
+          if (!err && lam1 > eps && lam1 < 1-eps &&
+              lam2 > eps && lam2 < 1-eps)
+            {
+              ok = false;
+              break;
+            }
+        }
+    }
+  
+  else
+    {      
+      // NgProfiler::RegionTimer reg(timer);      
     for (int j = 1; j <= nseg; j++)
       {
         const STLBoundarySeg & seg = GetSegment(j);
@@ -1143,7 +1212,7 @@ int STLBoundary :: TestSegChartNV(const Point3d & p1, const Point3d& p2,
           }
       }
     
-  NgProfiler::StopTimer (timer);
+    }
 
   return ok;
 }

libsrc/stlgeom/stltool.hpp

@@ -45,7 +45,7 @@ private:
   STLGeometry * geometry;
   Array<int>* charttrigs; // trigs which only belong to this chart
   Array<int>* outertrigs; // trigs which belong to other charts
-  Box3dTree * searchtree; // ADT containing outer trigs
+  BoxTree<3> * searchtree; // ADT containing outer trigs
 
   Array<twoint>* olimit; //outer limit of outer chart
   Array<twoint>* ilimit; //outer limit of inner chart
@@ -160,9 +160,10 @@ private:
   STLGeometry * geometry;
   const STLChart * chart;
   Array<STLBoundarySeg> boundary;
+  BoxTree<2> * searchtree = nullptr;
 public:
   STLBoundary(STLGeometry * ageometry);
-  // : boundary() {};
+  ~STLBoundary() { delete searchtree; }
 
   void Clear() {boundary.SetSize(0);};
   void SetChart (const STLChart * achart) { chart = achart; }
@@ -175,6 +176,8 @@ public:
   int NOSegments() {return boundary.Size();};
   const STLBoundarySeg & GetSegment(int i) {return boundary.Get(i);}
 
+  void BuildSearchTree();
+  void DeleteSearchTree();
   int TestSeg(const Point<3> & p1, const Point<3> & p2, const Vec<3> & sn, 
 	      double sinchartangle, int divisions, Array<Point<3> >& points,
 	      double eps);

libsrc/stlgeom/stltopology.hpp

@@ -252,7 +252,7 @@ protected:
   
   // searchtree for trigs and points
 
-  Box3dTree * searchtree; // ADT
+  BoxTree<3> * searchtree; // ADT
   Point3dTree * pointtree;
 
   Box<3> boundingbox;