Chombo + EB  3.2
IndicesTransformation.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 _INDICESTRANSFORMATION_H_
12 #define _INDICESTRANSFORMATION_H_
13 
14 #include <iostream>
15 #include "IntVect.H"
16 #include "RealVect.H"
17 #include "BaseFab.H"
18 #include "Tuple.H"
19 #include "IndicesFunctions.H"
20 
21 #include "NamespaceHeader.H"
22 
23 // ---------------------------------------------------------
24 /// Class to describe transformation of indices from one block to another.
26 {
27 public:
28 
29  /// null constructor leaves object in undefined state.
31 
32  ///
33  /**
34  Constructor.
35 
36  IntVect pOld: indices of a cell in old index space
37  IntVect pNew: indices of the same cell in new index space
38 
39  Then for each direction idir:
40  pNew[idir] = a_sign[idir]*pOld[a_permutation[idir]] + a_translation[idir].
41  */
42  IndicesTransformation(const IntVect& a_permutation,
43  const IntVect& a_sign,
44  const IntVect& a_translation);
45 
46  ///
47  /*
48  Same as constructor.
49  */
50  void define(const IntVect& a_permutation,
51  const IntVect& a_sign,
52  const IntVect& a_translation);
53 
54  ///
55  /*
56  Define an IndicesTransformation that is a swap of two indices.
57  Is Identity if the indices are equal.
58  */
59  void defineFromSwap(int a_ind1,
60  int a_ind2);
61 
62  ///
63  /*
64  Define an IndicesTransformation that is just a translation.
65  */
66  void defineFromTranslation(const IntVect& a_translation);
67 
68  ///
69  /*
70  Define an IndicesTransformation
71  such that the NODE a_pivotOld maps to the NODE a_pivotNew
72  and for each idir, [idir] maps to [a_permutation[idir]]
73  with sign a_sign[idir].
74  */
75  void defineFromPivot(const IntVect& a_pivotOld,
76  const IntVect& a_pivotNew,
77  const IntVect& a_permutation,
78  const IntVect& a_sign);
79 
80  ///
81  /*
82  Define an IndicesTransformation
83  such that a_srcBox's face in dimension a_srcDim side a_srcSide
84  maps to a_dstBox's face in dimension a_dstDim side a_dstSide
85  and for each idir, a_sign[idir] is -1 or +1 according to whether
86  source direction idir is flipped or not in the opposite block.
87  */
88  void defineFromFaces(const Box& a_srcBox,
89  int a_srcDim,
90  Side::LoHiSide a_srcSide,
91  const Box& a_dstBox,
92  int a_dstDim,
93  Side::LoHiSide a_dstSide,
94  const IntVect& a_sign = IntVect::Unit);
95 
96  ///
97  /*
98  Is this transformation defined?
99  */
100  bool isDefined() const;
101 
102  bool operator==(const IndicesTransformation& a_itOther) const;
103 
104  bool operator!=(const IndicesTransformation& a_itOther) const;
105 
106  friend std::ostream& operator<< (std::ostream& a_os,
107  const IndicesTransformation& a_it);
108 
109  ///
110  /**
111  IntVect pOld: indices of a cell in old index space
112  IntVect pNew: indices of the same cell in new index space
113 
114  Then for each direction idir:
115  pNew[idir] = m_sign[idir]*pOld[m_permutation[idir]] + m_translation[idir]
116  and hence:
117  pOld[m_permutation[idir]] = m_sign[idir]*pNew[idir]
118  - m_sign[idir]*m_translation[idir]
119  */
120  IntVect transformFwd(const IntVect& a_ivOld) const;
121 
122  IntVect transformBack(const IntVect& a_ivNew) const;
123 
124  IntVect transform(const IntVect& a_iv, bool a_forward = true) const;
125 
126  ///
127  /**
128  Transform NODE indices.
129  */
130  IntVect transformNode(const IntVect& a_iv) const;
131 
132  ///
133  /**
134  Transform coordinates in mapped space, not an index.
135  */
136  RealVect transformMapped(const RealVect& a_pointOld,
137  const RealVect& a_dxOld,
138  const RealVect& a_dxNew) const;
139 
140  ///
141  /**
142  Transform a vector, not an index.
143  There may be permutation of indices and sign change, but no translation.
144  */
145  IntVect transformVectorFwd(const IntVect& a_vecOld) const;
146 
147  IntVect transformVectorBack(const IntVect& a_vecNew) const;
148 
149  IntVect transformVector(const IntVect& a_vec, bool a_forward = true) const;
150 
151  ///
152  /**
153  Transform indices of either cells or nodes or faces or edges.
154  */
155  IntVect transformWithType(const IntVect& a_iv,
156  const IntVect& a_tp,
157  bool a_forward = true) const;
158 
159  ///
160  /**
161  Transform the type of Box:
162  no change if a_tp is the cell type or the node type,
163  but if a face or edge type, then the transformed type may differ.
164  */
165  IntVect transformType(const IntVect& a_tp, bool a_forward = true) const;
166 
167  ///
168  /**
169  Transform a whole box.
170  */
171  Box transformFwd(const Box& a_bxOld) const;
172 
173  Box transformBack(const Box& a_bxNew) const;
174 
175  Box transform(const Box& a_bx, bool a_forward = true) const;
176 
177  ///
178  /**
179  Apply forward transformation on a Box within a BaseFab.
180  */
181  template<typename T> void transformFwd(BaseFab<T>& a_dstFab,
182  const BaseFab<T>& a_srcFab,
183  const Box& a_srcBox,
184  const Interval& a_dstIntvl,
185  const Interval& a_srcIntvl) const
186  {
187  transform(a_dstFab, a_srcFab, a_srcBox,
188  a_dstIntvl, a_srcIntvl, true);
189  }
190 
191  ///
192  /**
193  Apply inverse transformation on a Box within a BaseFab.
194  */
195  template<typename T> void transformBack(BaseFab<T>& a_dstFab,
196  const BaseFab<T>& a_srcFab,
197  const Box& a_srcBox,
198  const Interval& a_dstIntvl,
199  const Interval& a_srcIntvl) const
200  {
201  transform(a_dstFab, a_srcFab, a_srcBox,
202  a_dstIntvl, a_srcIntvl, false);
203  }
204 
205  ///
206  /**
207  Apply forward or inverse transformation on a Box within a BaseFab.
208  */
209  template<typename T> void transform(BaseFab<T>& a_dstFab,
210  const BaseFab<T>& a_srcFab,
211  const Box& a_srcBox,
212  const Interval& a_dstIntvl,
213  const Interval& a_srcIntvl,
214  bool a_forward = true) const
215  {
216 #if 0
217  CH_assert(a_srcFab.box().contains(a_srcBox));
218  Box transformedBox = transform(a_srcBox, a_forward);
219  CH_assert(a_dstFab.box().contains(transformedBox));
220  IntVect tp = a_srcBox.type();
221  int dstComp = a_dstIntvl.begin();
222  int srcComp = a_srcIntvl.begin();
223  int srcCompHi = a_srcIntvl.end();
224  for (; srcComp <= srcCompHi; ++srcComp, ++dstComp)
225  {
226  for (BoxIterator bit(a_srcBox); bit.ok(); ++bit)
227  {
228  IntVect ivSrc = bit();
229  IntVect ivDst = transformWithType(ivSrc, tp, a_forward);
230  a_dstFab(ivDst, dstComp) =
231  a_srcFab(ivSrc, srcComp);
232  }
233  }
234 #else
235  // const IndicesTransformation itThis = (a_forward) ? *this : this->inverse();
236 
237  if (a_dstIntvl.size() == 0)
238  { // No components to set.
239  return;
240  }
241  if (a_srcBox.isEmpty())
242  { // Nothing to set.
243  return;
244  }
245  CH_assert(a_dstIntvl.size() == a_srcIntvl.size());
246  // dstIV[idir] = sign[idir]*srcIV[perm[idir]] + translation[idir]
247  const Box& dstAllBox = a_dstFab.box();
248  const Box& srcAllBox = a_srcFab.box();
249  const IntVect& srcAllDims = srcAllBox.size();
250  CH_assert(srcAllBox.contains(a_srcBox));
251  const IntVect& srcType = a_srcBox.type();
252  // Box dstSetBox = boxtransform(a_srcBox, perm, sign, translation);
253  Box dstSetBox = transform(a_srcBox, a_forward);
254  // Check that dstSetBox, the transformation of a_srcBox, is in dstAllBox.
255  CH_assert(dstAllBox.contains(dstSetBox));
256  const IntVect& srcSetLo = a_srcBox.smallEnd();
257  const IntVect& srcSetDims = a_srcBox.size();
258  IntVect dstSetLo = transformWithType(srcSetLo, srcType, a_forward);
259  const size_t dstIndexLo =
260  FortranArrayIndex(dstSetLo, dstAllBox);
261  const size_t srcIndexLo =
262  FortranArrayIndex(srcSetLo, srcAllBox);
263  // Not unsigned, because increments can be negative.
264  Tuple<long long int, SpaceDim> srcInc, dstInc, srcAllInc;
265  for (int idir = 0; idir < SpaceDim; idir++)
266  { // What happens if you move by BASISV(idir) from origin of set box?
267  IntVect src1dir = srcSetLo + BASISV(idir);
268  IntVect dst1dir = transformWithType(src1dir, srcType, a_forward);
269  if ( ! (dstAllBox.contains(dst1dir)) )
270  {
271  // FortranArrayIndex(dst1dir, dstAllBox) will be negative.
272  dstInc[idir] = 0;
273  }
274  else
275  {
276  size_t dstIndex1dirAll =
277  FortranArrayIndex(dst1dir, dstAllBox);
278  // Need to be careful when subtracting unsigned ints.
279  if (dstIndex1dirAll >= dstIndexLo)
280  {
281  dstInc[idir] = dstIndex1dirAll - dstIndexLo;
282  }
283  else
284  { // dstIndexLo > dstIndex1dirAll
285  dstInc[idir] = -(dstIndexLo - dstIndex1dirAll);
286  }
287  }
288  // same as FortranArrayIndex(src1dir, srcAllBox) - srcIndexLo;
289  srcInc[idir] =
290  FortranArrayIndex(BASISV(idir), srcAllDims);
291  srcAllInc[idir] = srcInc[idir] * srcSetDims[idir];
292  }
293 
294  Tuple<size_t, SpaceDim> srcBegin, dstBegin, srcEnd;
295  srcBegin[CH_SPACEDIM-1] = srcIndexLo;
296  dstBegin[CH_SPACEDIM-1] = dstIndexLo;
297 
298  size_t srcIndex, dstIndex;
299 
300  for (int srcComp = a_srcIntvl.begin(), dstComp = a_dstIntvl.begin();
301  srcComp <= a_srcIntvl.end();
302  ++srcComp, ++dstComp)
303  {
304  Interval src1Intvl(srcComp, srcComp);
305  Interval dst1Intvl(dstComp, dstComp);
306  BaseFab<T> dst1Fab(dst1Intvl, a_dstFab);
307  const BaseFab<T> src1Fab(src1Intvl, const_cast<BaseFab<T>&>(a_srcFab));
308 
309  const T* src1Array = src1Fab.dataPtr();
310  T* dst1Array = dst1Fab.dataPtr();
311 
312 #if CH_SPACEDIM>5
313  srcEnd[5] = srcBegin[5] + srcAllInc[5];
314  for (srcBegin[4] = srcBegin[5], dstBegin[4] = dstBegin[5];
315  srcBegin[4] < srcEnd[5];
316  srcBegin[4] += srcInc[5], dstBegin[4] += dstInc[5])
317  {
318 #endif
319 #if CH_SPACEDIM>4
320  srcEnd[4] = srcBegin[4] + srcAllInc[4];
321  for (srcBegin[3] = srcBegin[4], dstBegin[3] = dstBegin[4];
322  srcBegin[3] < srcEnd[4];
323  srcBegin[3] += srcInc[4], dstBegin[3] += dstInc[4])
324  {
325 #endif
326 #if CH_SPACEDIM>3
327  srcEnd[3] = srcBegin[3] + srcAllInc[3];
328  for (srcBegin[2] = srcBegin[3], dstBegin[2] = dstBegin[3];
329  srcBegin[2] < srcEnd[3];
330  srcBegin[2] += srcInc[3], dstBegin[2] += dstInc[3])
331  {
332 #endif
333 #if CH_SPACEDIM>2
334  srcEnd[2] = srcBegin[2] + srcAllInc[2];
335  for (srcBegin[1] = srcBegin[2], dstBegin[1] = dstBegin[2];
336  srcBegin[1] < srcEnd[2];
337  srcBegin[1] += srcInc[2], dstBegin[1] += dstInc[2])
338  {
339 #endif
340 #if CH_SPACEDIM>1
341  srcEnd[1] = srcBegin[1] + srcAllInc[1];
342  for (srcBegin[0] = srcBegin[1], dstBegin[0] = dstBegin[1];
343  srcBegin[0] < srcEnd[1];
344  srcBegin[0] += srcInc[1], dstBegin[0] += dstInc[1])
345  {
346 #endif
347  srcEnd[0] = srcBegin[0] + srcAllInc[0];
348  for (srcIndex = srcBegin[0], dstIndex = dstBegin[0];
349  srcIndex < srcEnd[0];
350  srcIndex += srcInc[0], dstIndex += dstInc[0])
351  {
352 #ifndef NDEBUG
353  // No need to check dstIndex >=0, srcIndex >= 0,
354  // because they are unsigned ints.
355  CH_assert(srcIndex < srcAllBox.numPts());
356  CH_assert(dstIndex < dstAllBox.numPts());
357 #endif
358  dst1Array[dstIndex] = src1Array[srcIndex];
359  }
360 #if CH_SPACEDIM>1
361  }
362 #endif
363 #if CH_SPACEDIM>2
364  }
365 #endif
366 #if CH_SPACEDIM>3
367  }
368 #endif
369 #if CH_SPACEDIM>4
370  }
371 #endif
372 #if CH_SPACEDIM>5
373  }
374 #endif
375  } // loop over components
376 #endif
377  }
378 
379  ///
380  /**
381  Return the inverse transformation.
382  */
384 
385  ///
386  /**
387  Return the composite transformation:
388  that is, this transformation followed by a_next.
389  */
391 
392  ///
393  /**
394  Return this transformation with index spaces coarsened by a_refRatio.
395  */
396  IndicesTransformation coarsen(int a_refRatio) const;
397 
398  ///
399  /**
400  Return this transformation with index spaces refined by a_refRatio.
401  */
402  IndicesTransformation refine(int a_refRatio) const;
403 
404  static const IntVect NoPermutation;
405 
406  ///
407  /** The identity transformation, which has no permutation,
408  no sign change, and no translation.
409  */
411 
412  ///
413  /** The undefined transformation.
414  */
416 
418  {
419  return m_permutation;
420  }
421  IntVect getSign() const
422  {
423  return m_sign;
424  }
426  {
427  return m_translation;
428  }
429 
430  ///
431  /** Initializes IndicesTransformation::Identity .
432  */
433  static int InitStatics();
434 
435 protected:
436 
437  // Zero if undefined
439 
440  // Zero if undefined
442 
443  // Zero if undefined
445 };
446 
447 //
448 //// Static initialization. Gotta make sure there are no static object
449 //// definitions above here (except possibly stuff in headers). Otherwise,
450 //// the danger is that some static object's constructor calls
451 //// IndicesTransformation::Identity, the very things the following
452 //// definition is supposed to initialize.
453 //// (I got this from IntVect.H)
455 
456 #include "NamespaceFooter.H"
457 
458 #endif // include guard
bool ok()
Definition: BoxIterator.H:281
bool isDefined() const
IntVect transformBack(const IntVect &a_ivNew) const
IntVect getTranslation() const
Definition: IndicesTransformation.H:425
IndicesTransformation()
null constructor leaves object in undefined state.
IntVect transformNode(const IntVect &a_iv) const
#define CH_SPACEDIM
Definition: SPACE.H:51
#define CH_assert(cond)
Definition: CHArray.H:37
IntVect getPermutation() const
Definition: IndicesTransformation.H:417
void transform(BaseFab< T > &a_dstFab, const BaseFab< T > &a_srcFab, const Box &a_srcBox, const Interval &a_dstIntvl, const Interval &a_srcIntvl, bool a_forward=true) const
Definition: IndicesTransformation.H:209
IntVect transformFwd(const IntVect &a_ivOld) const
IntVect m_translation
Definition: IndicesTransformation.H:444
IndicesTransformation compose(const IndicesTransformation &a_next) const
IntVect getSign() const
Definition: IndicesTransformation.H:421
IntVect size() const
size functions
Definition: Box.H:1803
IntVect BASISV(int dir)
Definition: IntVect.H:1257
void defineFromFaces(const Box &a_srcBox, int a_srcDim, Side::LoHiSide a_srcSide, const Box &a_dstBox, int a_dstDim, Side::LoHiSide a_dstSide, const IntVect &a_sign=IntVect::Unit)
iterates through the IntVects of a Box
Definition: BoxIterator.H:37
IntVect transformType(const IntVect &a_tp, bool a_forward=true) const
void transformBack(BaseFab< T > &a_dstFab, const BaseFab< T > &a_srcFab, const Box &a_srcBox, const Interval &a_dstIntvl, const Interval &a_srcIntvl) const
Definition: IndicesTransformation.H:195
int size() const
Definition: Interval.H:75
void transformFwd(BaseFab< T > &a_dstFab, const BaseFab< T > &a_srcFab, const Box &a_srcBox, const Interval &a_dstIntvl, const Interval &a_srcIntvl) const
Definition: IndicesTransformation.H:181
IntVect transformWithType(const IntVect &a_iv, const IntVect &a_tp, bool a_forward=true) const
const int SpaceDim
Definition: SPACE.H:38
IndicesTransformation refine(int a_refRatio) const
static const IndicesTransformation Undefined
Definition: IndicesTransformation.H:415
size_t FortranArrayIndex(IntVect a_iv, IntVect a_dims)
Definition: IndicesFunctions.H:18
Structure for passing component ranges in code.
Definition: Interval.H:23
void defineFromPivot(const IntVect &a_pivotOld, const IntVect &a_pivotNew, const IntVect &a_permutation, const IntVect &a_sign)
static const IntVect Unit
Definition: IntVect.H:659
RealVect transformMapped(const RealVect &a_pointOld, const RealVect &a_dxOld, const RealVect &a_dxNew) const
static const IntVect NoPermutation
Definition: IndicesTransformation.H:404
void defineFromSwap(int a_ind1, int a_ind2)
Ordered Tuples for Types T.
Definition: Tuple.H:30
IntVect transformVector(const IntVect &a_vec, bool a_forward=true) const
LoHiSide
Definition: LoHiSide.H:27
static int InitStatics()
bool operator!=(const IndicesTransformation &a_itOther) const
int begin() const
Definition: Interval.H:97
bool isEmpty() const
{ Comparison Functions}
Definition: Box.H:1846
void defineFromTranslation(const IntVect &a_translation)
Definition: BaseFab.H:76
void define(const IntVect &a_permutation, const IntVect &a_sign, const IntVect &a_translation)
T * dataPtr(int a_n=0)
Definition: BaseFabImplem.H:328
IntVect type() const
Definition: Box.H:1789
static const IndicesTransformation Identity
Definition: IndicesTransformation.H:410
A Rectangular Domain on an Integer Lattice.
Definition: Box.H:465
A Real vector in SpaceDim-dimensional space.
Definition: RealVect.H:41
size_t numPts() const
IndicesTransformation inverse() const
bool contains(const IntVect &p) const
Definition: Box.H:1871
An integer Vector in SpaceDim-dimensional space.
Definition: CHArray.H:42
friend std::ostream & operator<<(std::ostream &a_os, const IndicesTransformation &a_it)
static int s_dummyForIndicesTransformationH
Definition: IndicesTransformation.H:454
Class to describe transformation of indices from one block to another.
Definition: IndicesTransformation.H:25
IntVect m_permutation
Definition: IndicesTransformation.H:438
IndicesTransformation coarsen(int a_refRatio) const
IntVect transformVectorBack(const IntVect &a_vecNew) const
int end() const
Definition: Interval.H:102
IntVect transform(const IntVect &a_iv, bool a_forward=true) const
bool operator==(const IndicesTransformation &a_itOther) const
const IntVect & smallEnd() const
{ Accessors}
Definition: Box.H:1754
IntVect m_sign
Definition: IndicesTransformation.H:441
IntVect transformVectorFwd(const IntVect &a_vecOld) const
const Box & box() const
Definition: BaseFabImplem.H:226