MeshMap: Grid(MPI).h allocate their own copy of MeshMap

include/Cubism/Grid.h

@@ -39,7 +39,6 @@ protected:
     const double maxextent;
     const unsigned int N, NX, NY, NZ;
 
-    const bool m_own_mesh_maps;
     std::vector<MeshMap<Block>*> m_mesh_maps;
 
     void _dealloc()
@@ -48,13 +47,10 @@ protected:
 
         alloc.deallocate(m_blocks, N);
 
-        if (m_own_mesh_maps)
+        for (size_t i = 0; i < m_mesh_maps.size(); ++i)
         {
-            for (size_t i = 0; i < m_mesh_maps.size(); ++i)
-            {
-                delete m_mesh_maps[i];
-                m_mesh_maps[i] = NULL;
-            }
+            delete m_mesh_maps[i];
+            m_mesh_maps[i] = NULL;
         }
     }
 
@@ -103,8 +99,7 @@ public:
     typedef typename Block::RealType Real;  // Block MUST provide `RealType`.
 
     Grid(const unsigned int _NX, const unsigned int _NY = 1, const unsigned int _NZ = 1, const double _maxextent = 1) :
-        m_blocks(NULL), maxextent(_maxextent), N(_NX*_NY*_NZ), NX(_NX), NY(_NY), NZ(_NZ),
-        m_own_mesh_maps(true)
+        m_blocks(NULL), maxextent(_maxextent), N(_NX*_NY*_NZ), NX(_NX), NY(_NY), NZ(_NZ)
     {
         _alloc();
 
@@ -133,19 +128,20 @@ public:
                 }
     }
 
-    Grid(MeshMap<Block>* const mapX, MeshMap<Block>* const mapY, MeshMap<Block>* const mapZ,
-            const int _NX, const int _NY=1, const int _NZ=1) :
+    Grid(const MeshMap<Block>* const mapX,
+         const MeshMap<Block>* const mapY,
+         const MeshMap<Block>* const mapZ,
+         const int _NX, const int _NY=1, const int _NZ=1) :
         m_blocks(NULL),
         maxextent(-1.0), // not used
         N(_NX*_NY*_NZ),
-        NX(_NX), NY(_NY), NZ(_NZ),
-        m_own_mesh_maps(false)
+        NX(_NX), NY(_NY), NZ(_NZ)
     {
         _alloc();
 
-        m_mesh_maps.push_back(mapX);
-        m_mesh_maps.push_back(mapY);
-        m_mesh_maps.push_back(mapZ);
+        m_mesh_maps.push_back(new MeshMap<Block>(*mapX));
+        m_mesh_maps.push_back(new MeshMap<Block>(*mapY));
+        m_mesh_maps.push_back(new MeshMap<Block>(*mapZ));
 
         for(unsigned int iz=0; iz<NZ; iz++)
             for(unsigned int iy=0; iy<NY; iy++)
@@ -154,7 +150,11 @@ public:
                     const long long blockID = _encode(ix, iy, iz);
                     const int idx[3] = {(int)ix, (int)iy, (int)iz};
 
-                    m_vInfo.push_back(BlockInfo(blockID, idx, mapX, mapY, mapZ, _linaccess(blockID)));
+                    m_vInfo.push_back(BlockInfo(blockID, idx,
+                                                m_mesh_maps[0], // mmapX
+                                                m_mesh_maps[1], // mmapY
+                                                m_mesh_maps[2], // mmapZ
+                                                _linaccess(blockID)));
                 }
     }
 

include/Cubism/GridMPI.h

@@ -160,14 +160,16 @@ public:
     }
 
 
-    GridMPI(MeshMap<Block>* const mapX, MeshMap<Block>* const mapY, MeshMap<Block>* const mapZ,
+    GridMPI(const MeshMap<Block>* const mapX,
+            const MeshMap<Block>* const mapY,
+            const MeshMap<Block>* const mapZ,
             const int npeX, const int npeY, const int npeZ,
             const int nX=0, const int nY=0, const int nZ=0,
             const MPI_Comm comm = MPI_COMM_WORLD): TGrid(mapX,mapY,mapZ,nX,nY,nZ), timestamp(0), worldcomm(comm)
     {
-        assert(mapX->nblocks() == static_cast<size_t>(npeX*nX));
-        assert(mapY->nblocks() == static_cast<size_t>(npeY*nY));
-        assert(mapZ->nblocks() == static_cast<size_t>(npeZ*nZ));
+        assert(this->m_mesh_maps[0]->nblocks() == static_cast<size_t>(npeX*nX));
+        assert(this->m_mesh_maps[1]->nblocks() == static_cast<size_t>(npeY*nY));
+        assert(this->m_mesh_maps[2]->nblocks() == static_cast<size_t>(npeZ*nZ));
 
         blocksize[0] = Block::sizeX;
         blocksize[1] = Block::sizeY;
@@ -196,7 +198,6 @@ public:
 
         const std::vector<BlockInfo> vInfo = TGrid::getBlocksInfo();
 
-        MeshMap<Block>* const ptr_map[3] = {mapX, mapY, mapZ};
         for(size_t i=0; i<vInfo.size(); ++i)
         {
             BlockInfo info = vInfo[i];
@@ -205,11 +206,11 @@ public:
             {
                 info.index[j] += mypeindex[j]*mybpd[j];
 
-                info.origin[j] = ptr_map[j]->block_origin(info.index[j]);
+                info.origin[j] = this->m_mesh_maps[j]->block_origin(info.index[j]);
 
-                info.block_extent[j] = ptr_map[j]->block_width(info.index[j]);
+                info.block_extent[j] = this->m_mesh_maps[j]->block_width(info.index[j]);
 
-                info.ptr_grid_spacing[j] = ptr_map[j]->get_grid_spacing(info.index[j]);
+                info.ptr_grid_spacing[j] = this->m_mesh_maps[j]->get_grid_spacing(info.index[j]);
             }
 
             cached_blockinfo.push_back(info);

include/Cubism/MeshMap.h

@@ -12,6 +12,7 @@
 #include <cassert>
 #include <cstddef>
 #include <string>
+#include <cstring>
 
 #include "Cubism/Common.h"
 
@@ -60,7 +61,6 @@ template <typename TBlock>
 class MeshMap
 {
 public:
-    // constructor used to assume uniform cells in all directions!FIX EVERYWHERE
     MeshMap(const double xS, const double xE, const unsigned int Nblocks,
         const int nElemPerBlock = TBlock::sizeX) :
         m_xS(xS), m_xE(xE), m_extent(xE-xS), m_Nblocks(Nblocks),
@@ -68,6 +68,20 @@ public:
         m_uniform(true), m_initialized(false)
     {}
 
+    MeshMap(const MeshMap& c):
+        m_xS(c.m_xS), m_xE(c.m_xE), m_extent(c.m_extent), m_Nblocks(c.m_Nblocks),
+        m_Ncells(c.m_Ncells), m_blockSize(c.m_blockSize),
+        m_uniform(c.m_uniform), m_initialized(c.m_initialized),
+        m_kernel_name(c.m_kernel_name)
+    {
+        if (m_initialized)
+        {
+            _alloc();
+            std::memcpy(m_grid_spacing, c.m_grid_spacing, m_Ncells*sizeof(double));
+            std::memcpy(m_block_spacing, c.m_block_spacing, m_Nblocks*sizeof(double));
+        }
+    }
+
     ~MeshMap()
     {
         if (m_initialized)
@@ -77,6 +91,8 @@ public:
         }
     }
 
+    MeshMap& operator=(const MeshMap& rhs) = delete;
+
     void init(const MeshDensity* const kernel, const unsigned int ghostS=0, const unsigned int ghostE=0, double* const ghost_spacing=NULL)
     {
         _alloc();