easy-iso

commit bbcc969a7b81ae81862f540da06339df5c4c6f6e
parent 0fd21eba6944988e88b4fc8e9a6801b6e7ddf768
Author: Fabian Wermelinger <fabianw@mavt.ethz.ch>
Date:   Wed, 30 Nov 2016 21:14:25 +0100

added OpenMP on node level

Diffstat:
M
src/Interpolator.h
 | 
1
+
M
src/InterpolatorMPI.h
 | 
1
+
M
src/IsoExtractor.h
 | 
109
+++++++++++++++++++++++++++++++++++++++++++++++++++----------------------------

3 files changed, 73 insertions(+), 38 deletions(-)
diff --git a/src/Interpolator.h b/src/Interpolator.h
@@ -228,6 +228,7 @@ void Interpolator<TKernel>::_load_wavelet(ArgumentParser& p)
 
     Real blockdata[_BLOCKSIZE_][_BLOCKSIZE_][_BLOCKSIZE_];
     const int nBlocks = NBX * NBY * NBZ;
+#pragma omp parallel for
     for (int i = 0; i < nBlocks; ++i)
     {
         const int ix = i%NBX;
diff --git a/src/InterpolatorMPI.h b/src/InterpolatorMPI.h
@@ -354,6 +354,7 @@ void InterpolatorMPI<Tinterp>::_load_wavelet_MPI(ArgumentParser& p)
 
     Real blockdata[_BLOCKSIZE_][_BLOCKSIZE_][_BLOCKSIZE_];
     const int nBlocks = myNBX * myNBY * myNBZ;
+#pragma omp parallel for
     for (int i = 0; i < nBlocks; ++i)
     {
         const int ix = i%myNBX;
diff --git a/src/IsoExtractor.h b/src/IsoExtractor.h
@@ -422,14 +422,17 @@ protected:
     inline void _marchCube(const Real isoval,
             const Real x, const Real y, const Real z,
             const Real hx, const Real hy, const Real hz,
-            const Real Ox, const Real Oy, const Real Oz);
+            const Real Ox, const Real Oy, const Real Oz,
+            std::vector<Triangle>& tlist);
 
     inline void _marchTetrahedron(const Real isoval,
             const Real x, const Real y, const Real z,
             const Real hx, const Real hy, const Real hz,
-            const Real Ox, const Real Oy, const Real Oz);
+            const Real Ox, const Real Oy, const Real Oz,
+            std::vector<Triangle>& tlist);
 
-    inline void _evalTetrahedron(const Real coords[4][3], const Real val[4], const Real isoval);
+    inline void _evalTetrahedron(const Real coords[4][3], const Real val[4],
+            const Real isoval, std::vector<Triangle>& tlist);
 
     void _writePLY(const std::string basename, const std::vector<Triangle>& tList)
     {
@@ -516,39 +519,65 @@ void IsoExtractor<Tkernel,Tinterpolator>::_extractIso(const Real isoval,
         const Real chx, const Real chy, const Real chz,
         const Real Ox, const Real Oy, const Real Oz, const bool verbose)
 {
-    const int speculate = m_parser("speculate_triangle").asInt(10000);
-    m_triList.clear();
-    m_triList.reserve(speculate);
-    if (m_parser("march").asString("cube") == "cube")
-    {
-        for (int k = 0; k < Nz; ++k)
-            for (int j = 0; j < Ny; ++j)
-                for (int i = 0; i < Nx; ++i)
-                {
-                    const Real x = i*chx;
-                    const Real y = j*chy;
-                    const Real z = k*chz;
-                    _marchCube(isoval, x, y, z, chx, chy, chz, Ox, Oy, Oz);
-                }
-    }
-    else if (m_parser("march").asString("cube") == "tetrahedron")
-    {
-        for (int k = 0; k < Nz; ++k)
-            for (int j = 0; j < Ny; ++j)
-                for (int i = 0; i < Nx; ++i)
-                {
-                    const Real x = i*chx;
-                    const Real y = j*chy;
-                    const Real z = k*chz;
-                    _marchTetrahedron(isoval, x, y, z, chx, chy, chz, Ox, Oy, Oz);
-                }
-    }
-    else
+    const int speculate = m_parser("speculate_triangle").asInt(150000);
+    const std::string algorithm = m_parser("march").asString("cube");
+    if (!(algorithm == "cube" || algorithm == "tetrahedron"))
     {
         if (verbose)
-            std::cerr << "ERROR: IsoExtractor: Unknown algorithm " << m_parser("march").asString() << std::endl;
+            std::cerr << "ERROR: IsoExtractor: Unknown algorithm " << algorithm << std::endl;
         abort();
     }
+
+    const size_t nCubes = static_cast<size_t>(Nx) * Ny * Nz;
+    size_t totalSize = 0;
+#pragma omp parallel
+    {
+        std::vector<Triangle> myTriList;
+        myTriList.reserve(speculate);
+
+        if (algorithm == "cube")
+        {
+#pragma omp for
+            for (size_t i = 0; i < nCubes; ++i)
+            {
+                const size_t ix = i%Nx;
+                const size_t iy = (i/Nx)%Ny;
+                const size_t iz = (i/(Nx*Ny))%Nz;
+                const Real x = ix*chx;
+                const Real y = iy*chy;
+                const Real z = iz*chz;
+                _marchCube(isoval, x, y, z, chx, chy, chz, Ox, Oy, Oz, myTriList);
+            }
+        }
+        else if (algorithm == "tetrahedron")
+        {
+#pragma omp for
+            for (size_t i = 0; i < nCubes; ++i)
+            {
+                const size_t ix = i%Nx;
+                const size_t iy = (i/Nx)%Ny;
+                const size_t iz = (i/(Nx*Ny))%Nz;
+                const Real x = ix*chx;
+                const Real y = iy*chy;
+                const Real z = iz*chz;
+                _marchTetrahedron(isoval, x, y, z, chx, chy, chz, Ox, Oy, Oz, myTriList);
+            }
+        }
+#pragma omp critical
+        {
+            totalSize += myTriList.size();
+        }
+#pragma omp master
+        {
+            m_triList.clear();
+            m_triList.reserve(totalSize);
+        }
+#pragma omp barrier
+#pragma omp critical
+        {
+            m_triList.insert(m_triList.end(), myTriList.begin(), myTriList.end());
+        }
+    }
 }
 
 
@@ -556,7 +585,8 @@ template <typename Tkernel, template<typename> class Tinterpolator>
 void IsoExtractor<Tkernel,Tinterpolator>::_marchCube(const Real isoval,
         const Real x, const Real y, const Real z,
         const Real hx, const Real hy, const Real hz,
-        const Real Ox, const Real Oy, const Real Oz)
+        const Real Ox, const Real Oy, const Real Oz,
+        std::vector<Triangle>& tlist)
 {
     Real valCubeVerts[8];
     Real coordEdgeVert[12][3];
@@ -608,7 +638,7 @@ void IsoExtractor<Tkernel,Tinterpolator>::_marchCube(const Real isoval,
             tri.vertices[j][1] = coordEdgeVert[k][1];
             tri.vertices[j][2] = coordEdgeVert[k][2];
         }
-        m_triList.push_back(tri);
+        tlist.push_back(tri);
     }
 }
 
@@ -617,7 +647,8 @@ template <typename Tkernel, template<typename> class Tinterpolator>
 void IsoExtractor<Tkernel,Tinterpolator>::_marchTetrahedron(const Real isoval,
         const Real x, const Real y, const Real z,
         const Real hx, const Real hy, const Real hz,
-        const Real Ox, const Real Oy, const Real Oz)
+        const Real Ox, const Real Oy, const Real Oz,
+        std::vector<Triangle>& tlist)
 {
     Real valCubeVerts[8];
     Real coordCubeVerts[8][3];
@@ -654,13 +685,15 @@ void IsoExtractor<Tkernel,Tinterpolator>::_marchTetrahedron(const Real isoval,
             coordTetraVerts[j][2] = coordCubeVerts[vertID][2];
             valTetraVerts[j] = valCubeVerts[vertID];
         }
-        _evalTetrahedron(coordTetraVerts, valTetraVerts, isoval);
+        _evalTetrahedron(coordTetraVerts, valTetraVerts, isoval, tlist);
     }
 }
 
 
 template <typename Tkernel, template<typename> class Tinterpolator>
-void IsoExtractor<Tkernel,Tinterpolator>::_evalTetrahedron(const Real coords[4][3], const Real val[4], const Real isoval)
+void IsoExtractor<Tkernel,Tinterpolator>::_evalTetrahedron(const Real coords[4][3], const Real val[4],
+        const Real isoval,
+        std::vector<Triangle>& tlist)
 {
     Real coordEdgeVert[6][3];
 
@@ -709,7 +742,7 @@ void IsoExtractor<Tkernel,Tinterpolator>::_evalTetrahedron(const Real coords[4][
             tri.vertices[j][1] = coordEdgeVert[k][1];
             tri.vertices[j][2] = coordEdgeVert[k][2];
         }
-        m_triList.push_back(tri);
+        tlist.push_back(tri);
     }
 }