more parallel in topology

libsrc/general/parthreads.hpp

@@ -95,6 +95,12 @@ void ParallelFor( int first, int next, const TFunc & f )
 }
 
 
+  
+  template<typename T>
+  inline atomic<T> & AsAtomic (T & d)
+  {
+    return reinterpret_cast<atomic<T>&> (d);
+  }
 
   typedef void (*TaskManager)(std::function<void(int,int)>);
   typedef void (*Tracer)(string, bool);  // false .. start, true .. stop

libsrc/general/table.hpp

@@ -160,6 +160,12 @@ public:
       data[i-BASE].size++;
     }
 
+  inline void ParallelAdd (int i, const T & acont)
+    {
+      auto oldval = AsAtomic (data[i-BASE].size)++;
+      ((T*)data[i-BASE].col)[oldval] = acont;
+    }
+
   /// Inserts element acont into row i. 1-based. Does not test if already used, assumes to have mem
   inline void AddSave1 (int i, const T & acont)
     {

libsrc/include/mystdlib.h

@@ -18,7 +18,7 @@
 #include <memory>
 #include <thread>
 #include <mutex>
-
+#include <atomic>
 
 #include <new>
 #include <string>

libsrc/meshing/clusters.cpp

@@ -16,7 +16,7 @@ namespace netgen
     ;
   }
 
-  void AnisotropicClusters ::  Update(TaskManager tm)
+  void AnisotropicClusters ::  Update(TaskManager tm, Tracer tracer)
   {
     static int timer = NgProfiler::CreateTimer ("clusters");
     static int timer1 = NgProfiler::CreateTimer ("clusters1");
@@ -212,7 +212,7 @@ namespace netgen
 
     do
       {
-
+        (*tracer) ("update cluster, identify", false);
 	cnt++;
 	changed = 0;
       
@@ -332,6 +332,7 @@ namespace netgen
 	      }
 	    */
 	  }
+        (*tracer) ("update cluster, identify", true);        
       }
     while (changed);
     NgProfiler::StopTimer(timer3);

libsrc/meshing/clusters.hpp

@@ -27,7 +27,7 @@ public:
   AnisotropicClusters (const Mesh & amesh);
   ~AnisotropicClusters();
 
-  void Update(TaskManager tm = &DummyTaskManager);
+  void Update(TaskManager tm = &DummyTaskManager, Tracer trace = &DummyTracer);
 
   int GetVertexRepresentant (int vnr) const
   { return cluster_reps.Get(vnr); }

libsrc/meshing/meshclass.cpp

@@ -5773,7 +5773,7 @@ namespace netgen
   {
     topology.Update(tm, tracer);
     (*tracer)("call update clusters", false);
-    clusters->Update(tm);
+    clusters->Update(tm, tracer);
     (*tracer)("call update clusters", true);
 #ifdef PARALLEL
     if (paralleltop)

libsrc/meshing/topology.cpp

@@ -1,17 +1,9 @@
 #include <mystdlib.h>
 #include "meshing.hpp"
-#include <atomic>
 
 namespace netgen
 {
 
-  template<typename T>
-  inline atomic<T> & AsAtomic (T & d)
-  {
-    return reinterpret_cast<atomic<T>&> (d);
-  }
-
-
   
   template <class T>
   void QuickSortRec (FlatArray<T> data,
@@ -397,13 +389,26 @@ namespace netgen
        });
     
     vert2element = TABLE<ElementIndex,PointIndex::BASE> (cnt);
+    /*
     for (ElementIndex ei = 0; ei < ne; ei++)
       {
 	const Element & el = (*mesh)[ei];
 	for (int j = 0; j < el.GetNV(); j++)
 	  vert2element.AddSave (el[j], ei);
       }
-
+    */
+    ParallelForRange
+      (tm, ne,
+       [&] (size_t begin, size_t end)
+       {
+         for (ElementIndex ei = begin; ei < end; ei++)
+           {
+             const Element & el = (*mesh)[ei];
+             for (int j = 0; j < el.GetNV(); j++)
+               vert2element.ParallelAdd (el[j], ei);
+           }
+       });
+    
     cnt = 0;
     /*
     for (SurfaceElementIndex sei = 0; sei < nse; sei++)