get rid of OpenMP

CMakeLists.txt

@@ -144,12 +144,6 @@ if(APPLE)
 #     set(CMAKE_SHARED_LIBRARY_SUFFIX ".so")
 endif(APPLE)
 
-#######################################################################
-if(NOT INTEL_MIC AND NOT WIN32 AND NOT APPLE)
-    find_package(OpenMP)
-    set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")
-endif(NOT INTEL_MIC AND NOT WIN32 AND NOT APPLE)
-
 #######################################################################
 find_package(ZLIB REQUIRED)
 include_directories(${ZLIB_INCLUDE_DIRS})

libsrc/general/optmem.cpp

@@ -14,6 +14,7 @@
 
 namespace netgen
 {
+  static mutex block_allocator_mutex;
 
   BlockAllocator :: BlockAllocator (unsigned asize, unsigned ablocks)
     : bablocks (0)
@@ -36,8 +37,8 @@ namespace netgen
   void * BlockAllocator :: Alloc ()
   {
     void * p;
-#pragma omp critical (BlockAllocator)
     {
+      lock_guard<mutex> guard(block_allocator_mutex); 
       //  return new char[size];
       if (!freelist)
         {
@@ -60,8 +61,8 @@ namespace netgen
 
   void BlockAllocator :: Free (void * p)
   {
-#pragma omp critical (BlockAllocator)
     {
+      lock_guard<mutex> guard(block_allocator_mutex); 
       if (bablocks.Size())
         {
           *(void**)p = freelist;

libsrc/general/parthreads.hpp

@@ -187,6 +187,27 @@ public:
 
 #endif
 
+// Simple ParallelFor function to replace OpenMP
+template<typename TFunc>
+void ParallelFor( int first, int next, const TFunc & f )
+{
+  int nthreads = thread::hardware_concurrency();
+  thread * threads = new thread[nthreads];
+  for (int i=0; i<nthreads; i++)
+    {
+      threads[i] = std::thread( [&] ()
+        {
+          int myfirst = first + (next-first)*i/nthreads;
+          int mynext = first + (next-first)*(i+1)/nthreads;
+          f(myfirst, mynext);
+        });
+    }
+
+  for (int i=0; i<nthreads; i++)
+    threads[i].join();
+  delete [] threads;
+}
+
 }
 
 #endif

libsrc/include/mystdlib.h

@@ -21,6 +21,8 @@
 #include <climits>
 #include <algorithm>
 #include <memory>
+#include <thread>
+#include <mutex>
 
 
 #include <new>
@@ -35,9 +37,6 @@
 #include <unistd.h>  // for usleep (only for parallel)
 #endif
 
-#ifdef _OPENMP
-#include <omp.h>
-#endif
 
 
 /*

libsrc/interface/nginterface.cpp

@@ -2363,8 +2363,6 @@ void LinkFunction ()
 
 void Ng_TclCmd(string cmd)
 {
-#pragma omp critical(tcltodo)
+  lock_guard<mutex> guard(tcl_todo_mutex);
-  {
+  *(multithread.tcl_todo) += cmd;
-    *(multithread.tcl_todo) += cmd;
-  }
 }

libsrc/meshing/global.cpp

@@ -68,6 +68,8 @@ namespace netgen
 
   Array<int> tets_in_qualclass;
 
+  mutex tcl_todo_mutex;
+
   int h_argc = 0;
   char ** h_argv = NULL;
 

libsrc/meshing/global.hpp

@@ -29,6 +29,8 @@ namespace netgen
 
   DLL_HEADER extern Array<int> tets_in_qualclass;
 
+  DLL_HEADER extern mutex tcl_todo_mutex;
+
   class multithreadt
   {
   public:

libsrc/meshing/meshclass.cpp

@@ -4,6 +4,8 @@
 namespace netgen
 {
 
+  static mutex buildsearchtree_mutex;
+
   Mesh :: Mesh ()
     : surfarea(*this)
   {
@@ -4101,8 +4103,8 @@ namespace netgen
   {
     if (elementsearchtreets == GetTimeStamp()) return;
 
-#pragma omp critical (buildsearchtree)
     {
+      std::lock_guard<std::mutex> guard(buildsearchtree_mutex);
       if (elementsearchtreets != GetTimeStamp())
         {
           NgLock lock(mutex);

libsrc/stlgeom/stlgeom.cpp

@@ -2138,13 +2138,13 @@ int STLGeometry :: CheckGeometryOverlapping()
     }
 
 
-#pragma omp parallel
   {
     Array<int> inters;
+    mutex inters_mutex;
     
-#pragma omp for
+    ParallelFor( 1, GetNT()+1, [&] (int first, int next)
-    for (int i = 1; i <= GetNT(); i++)
+    {
-      {
+      for (int i=first; i<next; i++) {
 	const STLTriangle & tri = GetTriangle(i);
 	
 	Point<3> tpmin = tri.box.PMin();
@@ -2176,7 +2176,7 @@ int STLGeometry :: CheckGeometryOverlapping()
 
 	    if (IntersectTriangleTriangle (&trip1[0], &trip2[0]))
 	      {
-#pragma omp critical
+                lock_guard<mutex> guard(inters_mutex);
 		{
 		  oltrigs++;
 		  PrintMessage(5,"Intersecting Triangles: trig ",i," with ",inters.Get(j),"!");
@@ -2186,6 +2186,7 @@ int STLGeometry :: CheckGeometryOverlapping()
 	      }
 	  }
       }
+    });
   }
   PrintMessage(3,"Check overlapping geometry ... ", oltrigs, " triangles overlap");
   return oltrigs;

libsrc/visualization/vssolution.cpp

@@ -2474,27 +2474,36 @@ namespace netgen
               NgProfiler::RegionTimer reg1 (timer1);
               
               int ne = mesh->GetNE();
-              double hminv = numeric_limits<double>::max();
-              double hmaxv = -numeric_limits<double>::max();
-              bool hhasit = false;
-#if defined _OPENMP && _OPENMP >= 201107
-#pragma omp parallel for reduction (max : hmaxv) reduction (min : hminv) reduction (|| : hhasit)
-#endif
-              for (int i = 0; i < ne; i++) 
-                {
-                  double val;
-                  bool considerElem = GetValue (sol, i, 0.333, 0.333, 0.333, comp, val);
-                  if (considerElem)
-                    {
-                      if (val > hmaxv) hmaxv = val;
-                      if (val < hminv) hminv = val;
-                      hhasit = true;
-                    }
-                }
 
-              minv = min(minv, hminv);
+              mutex min_mutex;
-              maxv = max(maxv, hmaxv);
+              mutex max_mutex;
-              hasit |= hhasit;
+
+              ParallelFor(0, ne, [&] (int first, int next)
+                {
+                  double minv_local = numeric_limits<double>::max();
+                  double maxv_local = -numeric_limits<double>::max();
+                  for (int i=first; i<next; i++)
+                    {
+                      double val;
+                      bool considerElem = GetValue (sol, i, 0.333, 0.333, 0.333, comp, val);
+                      if (considerElem)
+                        {
+                          if (val > maxv_local) maxv_local = val;
+                          if (val < minv_local) minv_local = val;
+                          hasit = true;
+                        }
+                    }
+                  if(minv_local < minv)
+                    {
+                      lock_guard<mutex> guard(min_mutex);
+                      minv = minv_local;
+                    }
+                  if(maxv_local > maxv)
+                    {
+                      lock_guard<mutex> guard(max_mutex);
+                      maxv = maxv_local;
+                    }
+                });
             }
           
           if (sol->draw_surface)

ng/ngpkg.cpp

@@ -710,8 +710,8 @@ namespace netgen
       strcat (lstring, " 0");
     Tcl_SetVar  (interp, "::status_tetqualclasses", lstring, 0);
 
-#pragma omp critical(tcltodo)
     {
+      lock_guard<mutex> guard(tcl_todo_mutex);
       if (multithread.tcl_todo->length())
         {
           Tcl_Eval (interp, multithread.tcl_todo->c_str());