Chombo + EB + MF  3.2
FourthOrderMappedFineInterpSup.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 _FOURTHORDERMAPPEDFINEINTERPSUP_H_
12 #define _FOURTHORDERMAPPEDFINEINTERPSUP_H_
13 
14 
15 /******************************************************************************/
16 /**
17  * \file
18  *
19  * \brief Supporting classes and routines for FourthOrderMappedFineInterp
20  *
21  *//***************************************************************************/
22 
23 #include <cstring>
24 
25 #include "NamespaceHeader.H"
26 
27 
28 /*******************************************************************************
29  *
30  * Routines
31  *
32  ******************************************************************************/
33 
34 /*--------------------------------------------------------------------*/
35 /// Calculate a binomial coefficient
36 /** \param[in] n
37  * \param[in] k
38  * \return Binomial coefficient
39  *//*-----------------------------------------------------------------*/
40 
41 inline int
42 binomial(const int n, int k)
43 {
44  CH_assert(k >= 0);
45  CH_assert(n >= k);
46  if (k == 0) return 1;
47  int cnum = n - k + 1;
48  int cden = k;
49  while (--k)
50  {
51  cnum *= n - k + 1;
52  cden *= k;
53  }
54  return cnum/cden;
55 }
56 
57 /*--------------------------------------------------------------------*/
58 /// Find the sequential index of a power
59 /** The powers are indexed as follows (i.e for \f$d=3\f$ dimensions) \verbatim
60  * int idx = 0
61  * for (px = 0; px <= m; ++px)
62  * for (py = 0; px+py <= m; ++py)
63  * for (pz = 0; px+py+pz <= m; ++py)
64  * ++idx \endverbatim
65  * where \f$m\f$ is the degree of the polynomial and we wish to find
66  * idx. To compute the sequential index of any given power, we can
67  * use the relations
68  * \f[
69  * \mbox{number of powers} = {m + d \choose d}
70  * \f]
71  * and
72  * \f[
73  * \sum_{j=k}^n {j\choose k} = { n+1 \choose k+1 }
74  * \f]
75  * With these, the amount to add to the sequential index for a power
76  * at a spatial index is total number of powers remaining at this
77  * spatial index (remainder of \f$m\f$) minus the number of powers
78  * not used at this spatial index. E.g, if \f$m=3\f$, \f$d=3\f$,
79  * and \f$p_x=1\f$, there are 2 powers left for the remaining 2
80  * dimensions,
81  * \f[
82  * {2+2 \choose 2}\,.
83  * \f]
84  * The increment to the sequential index is then
85  * \f[
86  * {2+3+1 \choose 2+1} - {2+2+1 \choose 2+1}
87  * \f]
88  * In general, this can be written for direction index \f$i\f$, in
89  * \f$d\f$ space dimensions, with \f$m_i\f$ giving the remaining
90  * available powers at \f$i\f$, and \f$p_i\f$ giving the power used
91  * at index \f$i\f$ as
92  * \f[
93  * {d-i+m_i \choose d-i} - {d-i+m_i-p_i \choose d-i } \,.
94  * \f]
95  *
96  * \param[in] a_m Degree of the polynomial
97  * \param[in] a_p Power for each direction
98  * \return Sequential index
99  *//*-----------------------------------------------------------------*/
100 
101 inline int
102 powerIndex(int a_m, IntVect a_p)
103 {
104  int idx = 0;
105  for (int iDir = 0; iDir != SpaceDim; ++iDir)
106  {
107  const int remdir = SpaceDim - iDir;
108  idx += binomial(remdir + a_m, remdir);
109  a_m -= a_p[iDir];
110  idx -= binomial(remdir + a_m, remdir);
111  }
112  return idx;
113 }
114 
115 
116 /*******************************************************************************
117  */
118 /// Coordinate transformations
119 /**
120  * Provides for coordinate transformations. The transformation is input in
121  * compact form where a vector of size SpaceDim describes for each direction,
122  * the corresponding direction (including sign) in transformed space.
123  * Directions are number from 1, e.g., [-2, -1, -3] transforms
124  * i -> -j, j -> -i, and k -> -k. For transformation of IntVects, all
125  * coordinate systems are "fixed" at cell IntVect::Zero such that the
126  * transformation of IntVect::Zero is also IntVect::Zero. See
127  * http://www.grc.nasa.gov/WWW/cgns/sids/cnct.html#GridConnectivity1to1
128  * for a more detailed explanation on transformations.
129  *
130  * \note
131  * <ul>
132  * <li> The transformation description can be passed to Chombo Fortran as
133  * a 1-D array
134  * \verbatim CHF_CONST_I1D(ct.xctm(), 2*SpaceDim) \endverbatim
135  * and then used with a routine that duplicates
136  * CoordTransform::transform.
137  * </ul>
138  *
139  ******************************************************************************/
140 
142 {
143 public:
144 
145  /// Default constructor
146  CoordTransform();
147 
148  /// Constructor taking a compact description of the transformation
149  CoordTransform(const int *const a_ctm);
150 
151  /// Copy constructor
152  CoordTransform(const CoordTransform &a_ct);
153 
154  /// Assignment
156 
157  /// Defines the transformation by expanding a compact description
158  void expand(const int *const a_ctm);
159 
160  /// Transforms an IntVect to the transformed space
161  IntVect transform(const IntVect& a_iv) const;
162 
163  /// Returns the transformed direction associated with direction 'i'
164  int indexAt(const int a_i) const;
165 
166  /// Returns the sign of the transformation associated with direction 'i'
167  int signAt(const int a_i) const;
168 
169  /// Returns the expanded transformation description
170  const int* xctm() const;
171 
172  /// Reverse the transformation to the other way
173  void reverse();
174 
175 private:
176  /// Sign function
177  int sign(const int a_i) const
178  {
179  return ( a_i < 0 ) ? -1 : 1;
180  }
181 
182  int m_xctm[2*SpaceDim]; ///< An expanded description of the
183  ///< transformation. For direction i:
184  ///< m_xctm[2*i] = new direction
185  ///< m_xctm[2*i+1] = sign
186 };
187 
188 
189 /*******************************************************************************
190  *
191  * Class CoordTransform: inline member definitions
192  *
193  ******************************************************************************/
194 
195 /*--------------------------------------------------------------------*/
196 // Default sets m_xctm everywhere to zero.
197 /** \param[in] a_ctm Compact description of the transformation
198  *//*-----------------------------------------------------------------*/
199 
200 inline
202 {
203  std::memset(m_xctm, 0, 2*SpaceDim*sizeof(int));
204 }
205 
206 /*--------------------------------------------------------------------*/
207 // Constructor taking a compact description of the transformation
208 /** \param[in] a_ctm Compact description of the transformation
209  *//*-----------------------------------------------------------------*/
210 
211 inline
212 CoordTransform::CoordTransform(const int *const a_ctm)
213 {
214  expand(a_ctm);
215 }
216 
217 /*--------------------------------------------------------------------*/
218 // Copy constructor
219 /** \param[in] a_ct CoordTransform to copy
220  *//*-----------------------------------------------------------------*/
221 
222 inline
224 {
225  std::memcpy(m_xctm, a_ct.m_xctm, 2*SpaceDim*sizeof(int));
226 }
227 
228 /*--------------------------------------------------------------------*/
229 // Assignment
230 /** \param[in] a_ct CoordTransform to copy
231  * \return Updated with assignment
232  *//*-----------------------------------------------------------------*/
233 
234 inline CoordTransform &
236 {
237  if (&a_ct != this)
238  {
239  std::memcpy(m_xctm, a_ct.m_xctm, 2*SpaceDim*sizeof(int));
240  }
241  return *this;
242 }
243 
244 /*--------------------------------------------------------------------*/
245 // Defines the transformation by expanding a compact description
246 /** \param[in] a_ctm Compact description of the transformation
247  *//*-----------------------------------------------------------------*/
248 
249 inline void
250 CoordTransform::expand(const int *const a_ctm)
251 {
252  for (int i = 0; i != SpaceDim; ++i)
253  {
254  m_xctm[2*i] = std::abs(a_ctm[i]) - 1;
255  m_xctm[2*i+1] = sign(a_ctm[i]);
256  }
257 }
258 
259 /*--------------------------------------------------------------------*/
260 // Transforms an IntVect to the transformed space
261 /** \param[in] a_iv IntVect to transform
262  * \return Tranformed IntVect
263  *//*-----------------------------------------------------------------*/
264 
265 inline IntVect
267 {
268  IntVect tiv;
269  D_EXPR6(tiv[m_xctm[ 0]] = m_xctm[ 1]*a_iv[0],
270  tiv[m_xctm[ 2]] = m_xctm[ 3]*a_iv[1],
271  tiv[m_xctm[ 4]] = m_xctm[ 5]*a_iv[2],
272  tiv[m_xctm[ 6]] = m_xctm[ 7]*a_iv[3],
273  tiv[m_xctm[ 8]] = m_xctm[ 9]*a_iv[4],
274  tiv[m_xctm[10]] = m_xctm[11]*a_iv[5]);
275  return tiv;
276 }
277 
278 /*--------------------------------------------------------------------*/
279 // Returns the transformed direction associated with direction 'i'
280 /** \param[in] a_i Input direction
281  * \return Transformed direction
282  *//*-----------------------------------------------------------------*/
283 
284 inline int
285 CoordTransform::indexAt(const int a_i) const
286 {
287  return m_xctm[2*a_i];
288 }
289 
290 /*--------------------------------------------------------------------*/
291 // Returns the sign of the transformation associated with direction
292 // 'i'
293 /** \param[in] a_i Input direction
294  * \return Sign associated with transformation
295  *//*-----------------------------------------------------------------*/
296 
297 inline int
298 CoordTransform::signAt(const int a_i) const
299 {
300  return m_xctm[2*a_i + 1];
301 }
302 
303 /*--------------------------------------------------------------------*/
304 // Returns the expanded transformation description
305 /** \return Exanded transformation decsription
306  *//*-----------------------------------------------------------------*/
307 
308 inline const int*
310 {
311  return m_xctm;
312 }
313 
314 /*--------------------------------------------------------------------*/
315 // Reverse the transformation to the other way
316 /*--------------------------------------------------------------------*/
317 
318 inline void
320 {
321  int origxctm[2*SpaceDim];
322  std::memcpy(origxctm, m_xctm, 2*SpaceDim*sizeof(int));
323  for (int i = 0; i != SpaceDim; ++i)
324  {
325  const int j = origxctm[2*i];
326  m_xctm[2*j] = i;
327  m_xctm[2*j+1] = origxctm[2*i+1];
328  }
329 }
330 
331 #include "NamespaceFooter.H"
332 #endif
int m_xctm[2 *SpaceDim]
Definition: FourthOrderMappedFineInterpSup.H:182
#define CH_assert(cond)
Definition: CHArray.H:37
IntVect transform(const IntVect &a_iv) const
Transforms an IntVect to the transformed space.
Definition: FourthOrderMappedFineInterpSup.H:266
int signAt(const int a_i) const
Returns the sign of the transformation associated with direction 'i'.
Definition: FourthOrderMappedFineInterpSup.H:298
int binomial(const int n, int k)
Calculate a binomial coefficient.
Definition: FourthOrderMappedFineInterpSup.H:42
CoordTransform()
Default constructor.
Definition: FourthOrderMappedFineInterpSup.H:201
Coordinate transformations.
Definition: FourthOrderMappedFineInterpSup.H:141
const int * xctm() const
Returns the expanded transformation description.
Definition: FourthOrderMappedFineInterpSup.H:309
void reverse()
Reverse the transformation to the other way.
Definition: FourthOrderMappedFineInterpSup.H:319
const int SpaceDim
Definition: SPACE.H:38
CoordTransform & operator=(const CoordTransform &a_ct)
Assignment.
Definition: FourthOrderMappedFineInterpSup.H:235
int sign(const int a_i) const
Sign function.
Definition: FourthOrderMappedFineInterpSup.H:177
void expand(const int *const a_ctm)
Defines the transformation by expanding a compact description.
Definition: FourthOrderMappedFineInterpSup.H:250
int powerIndex(int a_m, IntVect a_p)
Find the sequential index of a power.
Definition: FourthOrderMappedFineInterpSup.H:102
An integer Vector in SpaceDim-dimensional space.
Definition: CHArray.H:42
int indexAt(const int a_i) const
Returns the transformed direction associated with direction 'i'.
Definition: FourthOrderMappedFineInterpSup.H:285