Chombo + EB + MF  3.2
ShapeArray.H
Go to the documentation of this file.
1 #ifdef CH_LANG_CC
2 /*
3  * _______ __
4  * / ___/ / ___ __ _ / / ___
5  * / /__/ _ \/ _ \/ V \/ _ \/ _ \
6  * \___/_//_/\___/_/_/_/_.__/\___/
7  * Please refer to Copyright.txt, in Chombo's root directory.
8  */
9 #endif
10 
11 #ifndef _SHAPEARRAY_H_
12 #define _SHAPEARRAY_H_
13 
14 #include <cstddef>
15 #include <tuple>
16 #include "IntVect.H"
17 
18 #include "NamespaceHeader.H"
19 
20 #if CXXSTD>=14
21 namespace shape
22 {
23 
24 /// Returns length of dimension D from a tuple of lengths
25 /** If highest dimension, returns 0
26  * \tparam IIx Indexing type
27  * \tparam D Dimension to query length
28  * \tparam Tuple Tuple type
29  * \param[in] a_lens Tuple of all dimension lengths
30  * \return D = rank : 0
31  * D < rank : length of the dimension
32  */
33 template <typename IIx, std::size_t D, typename Tuple>
34 constexpr IIx getLen(const Tuple& a_lens) noexcept
35 {
36  static_assert(std::is_integral<
37  typename std::tuple_element<D, Tuple>::type
38  >::value,
39  "ShapeArray requires an integer dimension");
40  return (D == std::tuple_size<Tuple>::value - 1) ? 0 : std::get<D>(a_lens);
41 }
42 
43 
44 /*******************************************************************************
45  */
46 /// Shapes linear memory according to some rank and dimensions and gives access
47 /**
48  * Template recursion is used to store the size of each dimension and at D==1,
49  * a pointer to the data is also stored. Access is provided through
50  * multidimensional notation x[i2][i1][i0] for rank = 3. C (row) ordering is
51  * assumed, where i0 (D==1) is unit stride. The size of the highest dimension
52  * (D==rank) is not required and this space is reused during access.
53  *
54  * 3-D (rank==3) example:
55  * We want an array of size i=4 x j=3 x k=2 where i is unit stride.
56  * Construct with a pointer to linear memory and the size of each dimension.
57  * \code{.cpp}
58  * double *data = new double[4*3*2];
59  * ShapeArray<double, 3> array3d(data, 4, 3, 2);
60  * \endcode
61  * Alternatively, use the helper function where template parameters are
62  * deduced.
63  * \code{.cpp}
64  * auto array3d = make_MD_array(data, 4, 3, 2);
65  * \endcode
66  * Access the array using operator[] for each dimension
67  * \code{.cpp}
68  * for (int k = 0; k != 2; ++k)
69  * for (int j = 0; j != 3; ++j)
70  * for (int i = 0; i != 4; ++i)
71  * {
72  * assert(&array3d[k][j][i] == &data[(k*3 + j)*4 + i]);
73  * }
74  * \endcode
75  *
76  * Note that operator[] returns a ShapeArray object of rank D-1. When D=1,
77  * the data is returned. During access, the size of the highest dimension
78  * is modified to tabulate the math on the right-hand-side of the equality
79  * above; lower dimensions are propagated.
80  *
81  * Performance:
82  * Compilers are very effective at optimizing the access. Constructing the
83  * ShapeArray is costly, so you would not want to do that for every access.
84  * However access is equivalent in performance to built-in Fortran arrays or
85  * in C, casting to pointers-to-VLAs to make use of built-in C-style
86  * multidimensional indexing.
87  *
88  * Non-zero lower bounds:
89  * The indexing assumes a zero-based lower bound. If otherwise, you must
90  * offset the array data. For a Chombo box, this can be done with.
91  * \code{.cpp}
92  * double *offsetData = data + box.index(IntVect_zero);
93  * \endcode
94  * If you ever shift the box, make sure you reset the data with a new offset.
95  *
96  * \tparam T Type of data in the array
97  * \tparam D Rank of the array (during recursion, this is the
98  * current dimension)
99  * \tparam IIx Type of integer for indexing the array (default
100  * std::ptrdiff_t)
101  *
102  *//*+*************************************************************************/
103 
104 template <typename T, std::size_t D, typename IIx = std::ptrdiff_t>
105 struct ShapeArray : ShapeArray<T, D-1, IIx>
106 {
107  /// Default constructor
108  ShapeArray() noexcept
109  :
110  ShapeArray<T, D-1, IIx>()
111  { }
112 
113  /// Construct from integer sequence
114  /** \tparam Is... Integer sequence of array dimensions
115  * \param[in] a_data Linear memory to reshaped
116  * \param[in] a_lens Size of each dimension
117  */
118  template <typename... Is>
119  ShapeArray(T* a_data, Is... a_lens) noexcept
120  :
121  ShapeArray<T, D-1, IIx>(a_data, a_lens...),
122  m_thisDlen(getLen<IIx, D-1>(std::make_tuple(a_lens...)))
123  { }
124 
125  /// Construct from a static vector. The rank of the array must be SpaceDim+1
126  /** \tparam T2 Type of data in the vector
127  * \tparam D2 Size of the vector (must be >= SpaceDim)
128  * \param[in] a_data Linear memory to reshaped
129  * \param[in] a_lens Size of each dimension
130  */
131  // template <typename T2, stc::array_size_type D2>
132  // ShapeArray(T* a_data, const stc::Vector<T2, D2>& a_lens) noexcept
133  // :
134  // ShapeArray<T, D-1, IIx>(a_data, a_lens),
135  // m_thisDlen((D-1 == SpaceDim) ? 0 : a_lens[D-1])
136  // { }
137  ShapeArray(T* a_data, const IntVect& a_lens) noexcept
138  :
139  ShapeArray<T, D-1, IIx>(a_data, a_lens),
140  m_thisDlen((D-1 == SpaceDim) ? 0 : a_lens[D-1])
141  { }
142 
143  /// Copy constructor
144  ShapeArray(const ShapeArray& a_other) noexcept
145  :
146  ShapeArray<T, D-1, IIx>((ShapeArray<T, D-1, IIx>)a_other),
147  m_thisDlen(a_other.m_thisDlen)
148  { }
149 
150  /// Assignment operator
151  ShapeArray& operator=(const ShapeArray& a_other) = default;
152 
153  /// Define from integer sequence
154  /** \tparam Is... Integer sequence of array dimensions
155  * \param[in] a_data Linear memory to reshaped
156  * \param[in] a_lens Size of each dimension
157  */
158  template <typename... Is>
159  void define(T* a_data, Is... a_lens) noexcept
160  {
161  ShapeArray<T, D-1, IIx>::define(a_data, a_lens...);
162  m_thisDlen = getLen<IIx, D-1>(std::make_tuple(a_lens...));
163  }
164 
165  /// Define from a static vector. The rank of the array must be SpaceDim+1
166  /** \tparam T2 Type of data in the vector
167  * \tparam D2 Size of the vector (must be >= SpaceDim)
168  * \param[in] a_data Linear memory to reshaped
169  * \param[in] a_lens Size of each dimension
170  */
171  // template <typename T2, stc::array_size_type D2>
172  // void define(T* a_data, const stc::Vector<T2, D2>& a_lens) noexcept
173  // {
174  // ShapeArray<T, D-1, IIx>::define(a_data, a_lens);
175  // m_thisDlen = (D-1 == SpaceDim) ? 0 : a_lens[D-1];
176  // }
177  void define(T* a_data, const IntVect& a_lens) noexcept
178  {
179  ShapeArray<T, D-1, IIx>::define(a_data, a_lens);
180  m_thisDlen = (D-1 == SpaceDim) ? 0 : a_lens[D-1];
181  }
182 
183  /// Reset only the data
184  /** This is often used if the bounds change (e.g., a box is shifted)
185  * \param[in] a_data Linear memory to reshaped
186  */
187  void resetData(T* a_data) noexcept
188  {
189  ShapeArray<T, 1, IIx>::resetData(a_data);
190  }
191 
192  /// Make it obvious that the array is not defined
193  void clear() noexcept
194  {
195  ShapeArray<T, D-1, IIx>::clear();
196  m_thisDlen = 0;
197  }
198 
199  /// Access the data
200  /** \param[in] a_idx Index for this dimension
201  * \return A temporary ShapeArray of one lower dimension
202  */
203  auto operator[](const IIx a_idx) const noexcept
204  {
205  return ShapeArray<T, D-1, IIx>(m_thisDlen + a_idx,
206  // propagate the original array
207  (ShapeArray<T, D-1, IIx>)(*this));
208  }
209 
210  // Higher dimension object can access our data
211  friend struct ShapeArray<T, D+1, IIx>;
212 
213 protected:
214 
215  /// Access builder
216  /** Higher-dimension indexing is stored in m_thisDlen. Lower dimensions are
217  * propagated.
218  */
219  ShapeArray(const IIx a_idx, const ShapeArray& a_other) noexcept
220  :
221  ShapeArray<T, D-1, IIx>((ShapeArray<T, D-1, IIx>)a_other),
222  m_thisDlen(a_idx*a_other.m_thisDlen)
223  { }
224 
225  IIx m_thisDlen; ///< Length of this dimension (or storage
226  ///< of index math if this is a temporary
227  ///< during access)
228 };
229 
230 // Specialization for D==1
231 template <typename T, typename IIx>
232 struct ShapeArray<T, 1, IIx>
233 {
234  ShapeArray() noexcept
235  :
236  m_data(nullptr)
237  { }
238  template <typename... Is>
239  ShapeArray(T* a_data, Is... a_lens) noexcept
240  :
241  m_thisDlen(getLen<IIx, 0>(std::make_tuple(a_lens...))),
242  m_data(a_data)
243  { }
244  // template <typename T2, stc::array_size_type D2>
245  // ShapeArray(T* a_data, const stc::Vector<T2, D2>& a_lens) noexcept
246  // :
247  // m_thisDlen(a_lens[0]),
248  // m_data(a_data)
249  // {
250  // static_assert(D2 >= SpaceDim, "ShapeArray constructed from static vector "
251  // "with insufficient size");
252  // }
253  ShapeArray(T* a_data, const IntVect& a_lens) noexcept
254  :
255  m_thisDlen(a_lens[0]),
256  m_data(a_data)
257  { }
258  ShapeArray(const ShapeArray& a_other) noexcept
259  :
260  m_thisDlen(a_other.m_thisDlen),
261  m_data(a_other.m_data)
262  { }
263  ShapeArray& operator=(const ShapeArray& a_other) = default;
264  template <typename... Is>
265  void define(T* a_data, Is... a_lens) noexcept
266  {
267  m_thisDlen = getLen<IIx, 0>(std::make_tuple(a_lens...));
268  m_data = a_data;
269  }
270  // template <typename T2, stc::array_size_type D2>
271  // void define(T* a_data, const stc::Vector<T2, D2>& a_lens) noexcept
272  // {
273  // static_assert(D2 >= SpaceDim, "ShapeArray defined from static vector "
274  // "with insufficient size");
275  // m_thisDlen = a_lens[0];
276  // m_data = a_data;
277  // }
278  void define(T* a_data, const IntVect& a_lens) noexcept
279  {
280  m_thisDlen = a_lens[0];
281  m_data = a_data;
282  }
283  void resetData(T* a_data) noexcept
284  {
285  m_data = a_data;
286  }
287  void clear() noexcept
288  {
289  m_thisDlen = 0;
290  m_data = nullptr;
291  }
292 
293  /// Access the data
294  /** \param[in] a_idx Index for this dimension
295  * \return Modifiable data at requested element
296  */
297  T& operator[](const IIx a_idx) const noexcept
298  {
299  return *(m_data + m_thisDlen + a_idx);
300  }
301 
302  /// Get pointer to the 1-D data
303  T* data() const noexcept
304  {
305  return m_data + m_thisDlen;
306  }
307 
308  // Higher dimension object can access our data
309  friend struct ShapeArray<T, 2, IIx>;
310 
311 protected:
312 
313  /// Access builder
314  /** Higher-dimension indexing is stored in m_thisDlen.
315  */
316  ShapeArray(const IIx a_idx, const ShapeArray& a_other) noexcept
317  :
318  m_thisDlen(a_idx*a_other.m_thisDlen),
319  m_data(a_other.m_data)
320  { }
321 
322  IIx m_thisDlen; ///< Length of this dimension
323  T* m_data; ///< Linear data
324 };
325 
326 /// Helper function that builds a ShapeArray deducing all template parameters
327 /** \param[in] a_data Linear memory to reshaped
328  * \param[in] a_lens Size of each dimension
329  * \return ShapeArray with type given by a_data and rank defined by the
330  * number of dimension lengths.
331  */
332 template <typename T, typename... Is>
333 inline auto make_MD_array(T* a_data, Is... a_lens) noexcept
334 {
335  return ShapeArray<T, sizeof...(a_lens), int>(a_data, a_lens...);
336 }
337 
338 } // namespace shape
339 #endif
340 
341 #include "NamespaceFooter.H"
342 
343 #endif /* ! defined _SHAPEARRAY_H_ */
const int SpaceDim
Definition: SPACE.H:38
Ordered Tuples for Types T.
Definition: Tuple.H:30
An integer Vector in SpaceDim-dimensional space.
Definition: CHArray.H:42