Interpolator.h (9502B)
1 // File : Interpolator.h 2 // Date : Tue Nov 22 15:37:00 2016 3 // Author : Fabian Wermelinger 4 // Description: Data interpolator type 5 // Copyright 2016 ETH Zurich. All Rights Reserved. 6 #ifndef INTERPOLATOR_H_QCK6H0CI 7 #define INTERPOLATOR_H_QCK6H0CI 8 9 #include <cassert> 10 #include <cstdlib> 11 #include <iostream> 12 #include <cmath> 13 #include <cstring> 14 #include <string> 15 #include <functional> 16 #include <chrono> 17 #ifdef _USE_HDF_ 18 #include <hdf5.h> 19 #endif /* _USE_HDF_ */ 20 21 #include "common.h" 22 #include "ArgumentParser.h" 23 #ifdef _USE_CUBISMZ_ 24 #include "Reader_WaveletCompression.h" 25 #endif /* _USE_CUBISMZ_ */ 26 27 28 template <typename TKernel> 29 class Interpolator 30 { 31 public: 32 Interpolator() = default; 33 Interpolator(ArgumentParser& p) : 34 m_extent(p("extent").asDouble(1.0)) 35 { 36 // load the data 37 std::chrono::time_point<std::chrono::system_clock> start, end; 38 start = std::chrono::system_clock::now(); 39 if (p("load").asString("h5") == "h5") _load_h5(p); 40 else if (p("load").asString("h5") == "wavelet") _load_wavelet(p); 41 else 42 { 43 std::cerr << "ERROR: Interpolator: Undefined loader \"" << p("load").asString() << "\"" << std::endl; 44 abort(); 45 } 46 end = std::chrono::system_clock::now(); 47 std::chrono::duration<double> elapsed_seconds = end-start; 48 std::cout << "Loading data done. Elapsed time = " << elapsed_seconds.count() << "s" << std::endl; 49 50 // cell centered or nodal values? 51 if (p("data").asString("cell") == "cell") _nCells = [](const int N) { return N; }; 52 else if (p("data").asString("cell") == "node") _nCells = [](const int N) { return N-1; }; 53 else 54 { 55 std::cerr << "ERROR: Interpolator: Unknown data attribute \"" << p("data").asString() << "\"" << std::endl; 56 abort(); 57 } 58 59 const int Nmax = std::max(m_data.Nx(), std::max(m_data.Ny(), m_data.Nz())); 60 assert(Nmax > 1); 61 m_h = m_extent/_nCells(Nmax); 62 m_invh = 1.0/m_h; 63 m_extentX = m_h*_nCells(m_data.Nx()); 64 m_extentY = m_h*_nCells(m_data.Ny()); 65 m_extentZ = m_h*_nCells(m_data.Nz()); 66 67 // init coordinate transform 68 _physicalDataPos = [this](const int i) { return i*this->m_h; }; 69 m_origin_off[0] = 0.0; 70 m_origin_off[1] = 0.0; 71 m_origin_off[2] = 0.0; 72 m_isNodal = 1; 73 if (p("data").asString("cell") == "cell") 74 { 75 _physicalDataPos = [this](const int i) { return (0.5+i)*this->m_h; }; 76 m_origin_off[0] = 0.5*m_h; 77 m_origin_off[1] = 0.5*m_h; 78 m_origin_off[2] = 0.5*m_h; 79 m_isNodal = 0; 80 } 81 } 82 virtual ~Interpolator() {} 83 84 Real operator()(const Real x, const Real y, const Real z) const 85 { 86 const Real xo = x - m_origin_off[0]; 87 const Real yo = y - m_origin_off[1]; 88 const Real zo = z - m_origin_off[2]; 89 90 const int ix0 = _idx(xo); 91 const int iy0 = _idx(yo); 92 const int iz0 = _idx(zo); 93 94 Real interp = 0.0; 95 for (int iz = TKernel::start; iz < TKernel::end; ++iz) 96 for (int iy = TKernel::start; iy < TKernel::end; ++iy) 97 for (int ix = TKernel::start; ix < TKernel::end; ++ix) 98 { 99 const Real xj = _pos(ix0+ix); 100 const Real yj = _pos(iy0+iy); 101 const Real zj = _pos(iz0+iz); 102 const Real Mx = m_kernel((xo-xj)*m_invh); 103 const Real My = m_kernel((yo-yj)*m_invh); 104 const Real Mz = m_kernel((zo-zj)*m_invh); 105 interp += m_data(ix0+ix,iy0+iy,iz0+iz)*Mx*My*Mz; 106 } 107 return interp; 108 } 109 110 inline Real getH() const { return m_h; } 111 inline Real getExtent() const { return m_extent; } 112 inline Real getExtentX() const { return m_extentX; } 113 inline Real getExtentY() const { return m_extentY; } 114 inline Real getExtentZ() const { return m_extentZ; } 115 inline Real getDataPos(const int i) const { return _physicalDataPos(i); } 116 inline void setOrigin(const Real o[3]) { for (int i=0; i<3; ++i) m_origin_off[i] = o[i]; } 117 inline int getNx() const { return m_data.Nx(); } 118 inline int getNy() const { return m_data.Ny(); } 119 inline int getNz() const { return m_data.Nz(); } 120 inline int isNodal() const { return m_isNodal; } 121 122 protected: 123 Real m_extent; 124 Real m_extentX; 125 Real m_extentY; 126 Real m_extentZ; 127 Real m_origin_off[3]; 128 Real m_h, m_invh; 129 int m_isNodal; 130 131 Matrix_t m_data; 132 TKernel m_kernel; 133 134 // helper 135 std::function<int(const int)> _nCells; 136 std::function<Real(const int)> _physicalDataPos; 137 inline int _idx(const Real x) const { return std::floor(x*m_invh); } 138 inline Real _pos(const int i) const { return i*m_h; } 139 140 private: 141 // loader 142 void _load_h5(ArgumentParser& p); 143 void _load_wavelet(ArgumentParser& p); 144 }; 145 146 147 template <typename TKernel> 148 void Interpolator<TKernel>::_load_h5(ArgumentParser& parser) 149 { 150 #ifdef _USE_HDF_ 151 const std::string group = parser("hdf5_group").asString("/data"); 152 const std::string filename = parser("file").asString(""); 153 if (filename == "") 154 { 155 std::cerr << "ERROR: Interpolator: -file is not specified. No input file given" << std::endl; 156 abort(); 157 } 158 159 hid_t file_id, dataset_id, dataspace_id, file_dataspace_id; 160 hsize_t* dims; 161 hssize_t num_elem; 162 int rank, ndims, NCH; 163 int maxDim[3]; 164 Real* data; 165 166 file_id = H5Fopen(filename.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT); 167 dataset_id = H5Dopen(file_id, group.c_str(), H5P_DEFAULT); 168 file_dataspace_id = H5Dget_space(dataset_id); 169 rank = H5Sget_simple_extent_ndims(file_dataspace_id); 170 dims = new hsize_t[rank]; 171 ndims = H5Sget_simple_extent_dims(file_dataspace_id, dims, NULL); 172 num_elem = H5Sget_simple_extent_npoints(file_dataspace_id); 173 data = new Real[num_elem]; 174 maxDim[2] = dims[0]; 175 maxDim[1] = dims[1]; 176 maxDim[0] = dims[2]; 177 NCH = dims[3]; 178 dataspace_id = H5Screate_simple(rank, dims, NULL); 179 int status = H5Dread(dataset_id, HDF_PRECISION, dataspace_id, file_dataspace_id, H5P_DEFAULT, data); 180 181 /* release stuff */ 182 delete [] dims; 183 status = H5Dclose(dataset_id); 184 status = H5Sclose(dataspace_id); 185 status = H5Sclose(file_dataspace_id); 186 status = H5Fclose(file_id); 187 188 if (NCH > 1) 189 { 190 const int channel = parser("channel").asInt(0); // data channel 191 const bool magnitude = parser("magnitude").asBool(false); // vector magnitude (only if NCH == 3) 192 if (magnitude && NCH == 3) 193 { 194 int k = 0; 195 for (int i = 0; i < num_elem; i += NCH) 196 { 197 const Real a = data[i+0]; 198 const Real b = data[i+1]; 199 const Real c = data[i+2]; 200 data[k++] = std::sqrt(a*a + b*b + c*c); 201 } 202 203 } 204 else 205 { 206 assert(channel < NCH); 207 int k = 0; 208 for (int i = 0; i < num_elem; i += NCH) 209 data[k++] = data[i+channel]; 210 } 211 } 212 this->m_data = std::move(Matrix_t(maxDim[0], maxDim[1], maxDim[2], data)); 213 delete [] data; 214 #endif /* _USE_HDF_ */ 215 } 216 217 template <typename TKernel> 218 void Interpolator<TKernel>::_load_wavelet(ArgumentParser& p) 219 { 220 #ifdef _USE_CUBISMZ_ 221 const bool byte_swap = p("-swap").asBool(false); 222 const int wavelet_type = p("-wtype").asInt(1); 223 const std::string filename = p("file").asString(""); 224 if (filename == "") 225 { 226 std::cerr << "ERROR: Interpolator: -file is not specified. No input file given" << std::endl; 227 abort(); 228 } 229 230 Reader_WaveletCompression myreader(filename, byte_swap, wavelet_type); 231 myreader.load_file(); 232 const int NBX = myreader.xblocks(); 233 const int NBY = myreader.yblocks(); 234 const int NBZ = myreader.zblocks(); 235 const int maxDim[3] = {NBX*_BLOCKSIZE_, NBY*_BLOCKSIZE_, NBZ*_BLOCKSIZE_}; 236 const size_t num_elem = static_cast<size_t>(maxDim[0]) * maxDim[1] * maxDim[2]; 237 Real* const data = new Real[num_elem]; 238 239 const int nBlocks = NBX * NBY * NBZ; 240 #pragma omp parallel 241 { 242 Real blockdata[_BLOCKSIZE_][_BLOCKSIZE_][_BLOCKSIZE_]; 243 #pragma omp for 244 for (int i = 0; i < nBlocks; ++i) 245 { 246 const int ix = i%NBX; 247 const int iy = (i/NBX)%NBY; 248 const int iz = (i/(NBX*NBY))%NBZ; 249 const double zratio = myreader.load_block2(ix, iy, iz, blockdata); 250 251 for (int z=0; z < _BLOCKSIZE_; ++z) 252 for (int y=0; y < _BLOCKSIZE_; ++y) 253 { 254 assert(iy*_BLOCKSIZE_+y < maxDim[1]); 255 assert(iz*_BLOCKSIZE_+z < maxDim[2]); 256 const size_t offset = _BLOCKSIZE_*(static_cast<size_t>(ix) + NBX*(y+iy*_BLOCKSIZE_ + (z+static_cast<size_t>(iz)*_BLOCKSIZE_)*NBY*_BLOCKSIZE_)); 257 Real* const dst = data + offset; 258 memcpy(dst, &blockdata[z][y][0], _BLOCKSIZE_*sizeof(Real)); 259 } 260 } 261 } 262 263 this->m_data = std::move(Matrix_t(maxDim[0], maxDim[1], maxDim[2], data)); 264 delete [] data; 265 #else 266 fprintf(stderr, "WARNING: Executable was compiled without wavelet compressor support...\n"); 267 #endif /* _USE_CUBISMZ_ */ 268 } 269 270 #endif /* INTERPOLATOR_H_QCK6H0CI */