1//
2// Created by tommaso on 17/06/19.
3//
4
5#ifndef OPENFPM_PDATA_BLOCKGEOMETRY_HPP
6#define OPENFPM_PDATA_BLOCKGEOMETRY_HPP
7
8#include <boost/mpl/size_t.hpp>
9#include <cstring>
10#include <Grid/grid_sm.hpp>
11#include "SparseGridGpu/TemplateUtils/mathUtils.hpp"
12
13/**
14 * This class provides an interface to linearization of coordinates and viceversa when blocks are involved.
15 * This can be seen as a lightweight version of grid_sm, with just LinId and InvLinId methods, but
16 * tuned for blocked data.
17 */
18template<unsigned int dim, unsigned int blockEdgeSize, typename indexT = long int>
19class grid_smb
20{
21private:
22
23 indexT blockSz[dim];
24 indexT sz[dim];
25
26protected:
27
28 constexpr static indexT blockSize = IntPow<blockEdgeSize, dim>::value;
29
30public:
31
32 grid_smb() {}
33
34 __host__ __device__ grid_smb(const size_t (& sz)[dim])
35 {
36 for (int d=0; d<dim; ++d)
37 {
38 this->sz[d] = sz[d];
39 blockSz[d] = sz[d] / blockEdgeSize + ((sz[d] % blockEdgeSize) != 0);
40 }
41 }
42
43 __host__ __device__ grid_smb(const indexT (& sz)[dim])
44 {
45 for (int d=0; d<dim; ++d)
46 {
47 this->sz[d] = sz[d];
48 blockSz[d] = sz[d] / blockEdgeSize + ((sz[d] % blockEdgeSize) != 0);
49 }
50 }
51
52 __host__ __device__ grid_smb(const size_t domainBlockEdgeSize)
53 {
54 for (int i = 0; i < dim; ++i)
55 {
56 blockSz[i] = (indexT)domainBlockEdgeSize;
57 sz[i] = (indexT)domainBlockEdgeSize * blockEdgeSize;
58 }
59 }
60
61 template<typename T>
62 __host__ __device__ grid_smb(const grid_sm<dim, T> blockGrid)
63 {
64 for (int i = 0; i < dim; ++i)
65 {
66 blockSz[i] = blockGrid.size(i);
67 sz[i] = blockGrid.size(i) * blockEdgeSize;
68 }
69 }
70
71#ifdef __NVCC__
72 //Constructors from dim3 and uint3 objects
73 __host__ __device__ grid_smb(const dim3 blockDimensions)
74 {
75 unsigned int i = 0;
76 assert(dim <= 3);
77 blockSz[i] = blockDimensions.x;
78 sz[i] = blockSz[i] * blockEdgeSize;
79 if (dim > 1)
80 {
81 ++i;
82 blockSz[i] = blockDimensions.y;
83 sz[i] = blockSz[i] * blockEdgeSize;
84 if (dim > 2)
85 {
86 ++i;
87 blockSz[i] = blockDimensions.z;
88 sz[i] = blockSz[i] * blockEdgeSize;
89 }
90 }
91 }
92
93
94#endif // __NVCC__
95
96 __host__ __device__ grid_smb(const grid_smb<dim, blockEdgeSize> &other)
97 {
98 for (indexT i = 0 ; i < dim ; i++)
99 {
100 blockSz[i] = other.blockSz[i];
101 sz[i] = other.sz[i];
102 }
103 }
104
105 __host__ __device__ grid_smb &operator=(const grid_smb<dim, blockEdgeSize> &other)
106 {
107 for (indexT i = 0 ; i < dim ; i++)
108 {
109 blockSz[i] = other.blockSz[i];
110 sz[i] = other.sz[i];
111 }
112 return *this;
113 }
114
115 __host__ __device__ const indexT (& getSize() const)[dim]
116 {
117 return sz;
118 }
119
120 __host__ __device__ const indexT & size(int i) const
121 {
122 return sz[i];
123 }
124
125 /*! \brief Linearize the coordinate index
126 *
127 * The linearization is given by getting the block indexes and the local coordinate indexes
128 *
129 * Linearize the block index (blockLinId), linearize the local index (localLinId) and return
130 * blockLinId * blockSize + localLinId
131 *
132 * \param coord coordinates
133 *
134 * \return linearized index
135 *
136 */
137 template<typename indexT_>
138 inline __host__ __device__ indexT LinId(const grid_key_dx<dim, indexT_> coord) const
139 {
140 //todo: Check (in debug mode only) that the coordinates passed here are valid and not overflowing dimensions (???)
141 indexT blockLinId = coord.get(dim - 1) / blockEdgeSize;
142 indexT localLinId = coord.get(dim - 1) % blockEdgeSize;
143 for (int d = dim - 2; d >= 0; --d)
144 {
145 blockLinId *= blockSz[d];
146 localLinId *= blockEdgeSize;
147 blockLinId += coord.get(d) / blockEdgeSize;
148 localLinId += coord.get(d) % blockEdgeSize;
149 }
150 return blockLinId * blockSize + localLinId;
151 }
152
153 inline __host__ __device__ grid_key_dx<dim, int> InvLinId(const indexT linId) const
154 {
155 indexT blockLinId = linId / blockSize;
156 indexT localLinId = linId % blockSize;
157 return InvLinId(blockLinId, localLinId);
158 }
159
160 /*! brief Given a local offset it provide the the position in coordinates
161 *
162 * \param linId point
163 *
164 * \return the point in coordinates
165 *
166 */
167 inline __host__ __device__ grid_key_dx<dim, int> LocalInvLinId(unsigned int localLinId) const
168 {
169 grid_key_dx<dim, int> coord;
170 for (int d = 0; d < dim; ++d)
171 {
172 auto c = localLinId % blockEdgeSize;
173 coord.set_d(d, c);
174 localLinId /= blockEdgeSize;
175 }
176 return coord;
177 }
178
179 /*! \brief Invert from the linearized block id + local id to the position of the point in coordinates
180 *
181 * From the point coordinated x,y,z you can get the block coordinated block_x,block_y,block_z
182 * and the local coordinates inside the chunk loc_x,loc_y,loc_z
183 *
184 * linearizing block coordinated you get blockId and linearizing the local coordinated you get localLinId
185 *
186 * each point in a sparse grid is identified by the formula blockId*blockSize + localLinId
187 *
188 * This function invert such formula.
189 *
190 * \param blockLinId is blockId
191 *
192 * \param localLinId is localLinId in the formula
193 *
194 *
195 * \return the spatial coordinates of the point
196 *
197 */
198 inline __host__ __device__ grid_key_dx<dim, int> InvLinId(indexT blockLinId, indexT localLinId) const
199 {
200 grid_key_dx<dim, int> coord;
201 for (int d = 0; d < dim; ++d)
202 {
203 auto c = blockLinId % blockSz[d];
204 c *= blockEdgeSize;
205 c += localLinId % blockEdgeSize;
206 coord.set_d(d, c);
207 blockLinId /= blockSz[d];
208 localLinId /= blockEdgeSize;
209 }
210 return coord;
211 }
212
213 // Now methods to handle blockGrid coordinates (e.g. to load neighbouring blocks)
214 template<typename indexT_>
215 inline __host__ __device__ indexT BlockLinId(const grid_key_dx<dim, indexT_> & blockCoord) const
216 {
217 indexT blockLinId = blockCoord.get(dim - 1);
218 if (blockLinId >= blockSz[dim-1])
219 {return -1;}
220
221 for (int d = dim - 2; d >= 0; --d)
222 {
223 blockLinId *= blockSz[d];
224 indexT cur = blockCoord.get(d);
225 if (cur >= blockSz[d])
226 {
227 return -1;
228 }
229 blockLinId += cur;
230 }
231 return blockLinId;
232 }
233
234 // Now methods to handle blockGrid coordinates (e.g. to load neighbouring blocks)
235 template<typename indexT_>
236 inline __host__ __device__ grid_key_dx<dim,indexT> getGlobalCoord(const grid_key_dx<dim, indexT_> & blockCoord, unsigned int offset) const
237 {
238 grid_key_dx<dim,indexT> k;
239
240 for (unsigned int i = 0 ; i < dim ; i++)
241 {
242 k.set_d(i,blockCoord.get(i)*blockEdgeSize + offset%blockEdgeSize);
243 offset /= blockEdgeSize;
244 }
245 return k;
246 }
247
248 inline __host__ __device__ grid_key_dx<dim, int> BlockInvLinId(indexT blockLinId) const
249 {
250 grid_key_dx<dim, int> blockCoord;
251 for (int d = 0; d < dim; ++d)
252 {
253 auto c = blockLinId % blockSz[d];
254 blockCoord.set_d(d, c);
255 blockLinId /= blockSz[d];
256 }
257 return blockCoord;
258 }
259
260 inline indexT size_blocks() const
261 {
262 indexT sz = 1;
263
264 for (indexT i = 0 ; i < dim ; i++)
265 {sz *= blockSz[i];}
266
267 return sz;
268 }
269
270 inline indexT size() const
271 {
272 indexT sz = 1;
273
274 for (indexT i = 0 ; i < dim ; i++)
275 {sz *= this->sz[i];}
276
277 return sz;
278 }
279
280 __host__ __device__ inline indexT getBlockSize() const
281 {
282 return blockSize;
283 }
284
285 __host__ __device__ inline void swap(grid_smb<dim, blockEdgeSize, indexT> &other)
286 {
287 indexT blockSz_tmp[dim];
288 indexT sz_tmp[dim];
289
290 for (indexT i = 0 ; i < dim ; i++)
291 {
292 blockSz_tmp[i] = blockSz[i];
293 blockSz[i] = other.blockSz[i];
294 other.blockSz[i] = blockSz_tmp[i];
295
296 sz_tmp[i] = sz[i];
297 sz[i] = other.sz[i];
298 other.sz[i] = sz_tmp[i];
299 }
300 }
301};
302
303#endif //OPENFPM_PDATA_BLOCKGEOMETRY_HPP
304