Chombo + EB + MF  3.2
IndexedMoments.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 _INDEXEDMOMENTS_H_
12 #define _INDEXEDMOMENTS_H_
13 
14 #include "MayDay.H"
15 #include "IndexTM.H"
16 #include "Vector.H"
17 #include "SPACE.H"
18 #include "CH_EBIS_ORDER.H"
19 
20 #include "NamespaceHeader.H"
21 using std::map;
22 
23 #define CH_IM_MAX_POWER 8
24 
25 //! \class IndexedMoments
26 //! Vector-like container for multi-indexed Real values up to
27 //! some max multi-index P (sum of indicies <= P).
28 //! Layout is 0th dim first, then 1st, etc.
29 //! \tparam D The Dimension of the container
30 template <int Dim, int P>
32 {
33 public:
34  /// constructor---make statics first time called
36 
37 
38  //! Destructor
40 
41 
42  //! Retrieve the moment from the index
43  //! \params a_index The lexigraphical index to the data
44  Real& operator[](int a_i)
45  {
46  return m_moms[a_i];
47  };
48 
49  //! Retrieve the moment from the index
50  //! \params a_index The lexigraphical index to the data
51  const Real& operator[](int a_i) const
52  {
53  return m_moms[a_i];
54  };
55 
56  //! Retrieve the moment from the index
57  //! \params a_index The multi-index that's needed
58  Real& operator[](const IndexTM<int,Dim>& a_index)
59  {
60  return m_moms[indexOf(a_index)];
61  };
62 
63  //! Retrieve the moment from the index
64  //! \params a_index The multi-index that's needed
65  const Real& operator[](const IndexTM<int,Dim>& a_index) const
66  {
67  return m_moms[indexOf(a_index)];
68  };
69 
70  ///add compoenentwise
72 
73  ///multiply each component by constant
74  IndexedMoments<Dim, P>& operator*=(const Real& a_factor);
75 
76 
77  ///number of reals in the vector
78  static int size()
79  {
80  if(!s_staticsSet)
81  {
82  setStatics();
83  }
84  return s_size;
85  }
86 
87  ///for linearization
88  static size_t linearSize()
89  {
90  int numReals = size();
91  //add in the real numbers
92  size_t retval = numReals*sizeof(Real);
93 
94  //add in the m_isRegular
95  retval += sizeof(int);
96 
97  return retval;
98  }
99 
100  void linearOut(void* const a_outbuf) const
101  {
102  Real* realbuf = (Real*) a_outbuf;
103  for(int ivec = 0; ivec < size(); ivec++)
104  {
105  *realbuf = m_moms[ivec];
106  ++realbuf;
107  }
108  int* intbuf = (int *) realbuf;
109  *intbuf = m_isRegular;
110  }
111 
112  void linearIn(const void* const a_inbuf)
113  {
114  Real* realbuf = (Real*) a_inbuf;
115  for(int ivec = 0; ivec < size(); ivec++)
116  {
117  m_moms[ivec]= *realbuf;
118  ++realbuf;
119  }
120  int* intbuf = (int *) realbuf;
121  m_isRegular = *intbuf ;
122  }
123 
124  /// set to a regular IndexTM
125  void setRegular(const Real a_dx);
126 
127  //sick of misspelling this one
128  void setToRegular(const Real a_dx)
129  {
130  setRegular(a_dx);
131  }
132 
133  ///monomial powers
135  {
136  if(!s_staticsSet)
137  {
138  setStatics();
139  }
140  return s_multiIndicies;
141  }
142 
143  ///
144  bool isRegular() const
145  {
146  return (m_isRegular==1);
147  }
148 
149  /// for use with irregnode
151 
152  ///
154  {
155  if(&a_input != this)
156  {
157  m_isRegular = a_input.m_isRegular;
158  m_moms = a_input.m_moms;
159  }
160  return *this;
161  }
162 
163 
164  ///
165  /**
166  shift moment by the input vector distance.
167  this changes the current object from integral(x^p)
168  to integral((x+x0)^p), where x0 = a_distance
169  */
170  void shift(const IndexTM<Real, Dim>& a_distance);
171 
172  ///
173  void setToZero()
174  {
175  for(int ivec = 0; ivec < s_size; ivec++)
176  {
177  m_moms[ivec] = 0.;
178  }
179  }
180 
181  /// Calculate what linear index this multi-index is
182  static int indexOf(IndexTM<int,Dim> a_index);
183 
184  ///
185  static IndexTM<int,Dim> getIndex(const int& a_linearIndex)
186  {
187  return s_multiIndicies[a_linearIndex];
188  }
189 
190  /// outputs the current state to pout() (a la parstream.H)
191  void spout() const;
192 
193  ///
194  /**
195  Say <A> = sum_p(CA m^p),
196  and <B> = sum_q(CB m^q).
197 
198  This sets the current data to
199  the set of coefficents M such that
200  <AB> = sum_r(M m^r) + O(h^P+1).
201 
202  We drop all coefficents for powers s.t. p + q > P.
203  */
205  const IndexedMoments<Dim, P> & a_CB);
206 
207 
208 
209  ///divides each entry by p!
210  void divideByFactorial();
211 
212  ///multiply each entry by p!
213  void multiplyByFactorial();
214 
215 protected:
216 
217  ///
218  Real
219  getMoment(const IndexTM<int, Dim> & a_mono,
220  const map<IndexTM<int, Dim>, Real>& a_mapin) const;
221 
222  ///
223  static void setStatics();
224 
225  ///
226  static bool s_staticsSet;
227 
228  ///
229  static int s_size;
230 
231  ///
233 
234  ///
235  static void setMultiIndicies();
236 
237  ///
238  static void setSize();
239 
240 private:
241 
242  // Indicator that we contain only "full" moments
244 
245  // Our moments to store
247 
248  static const int s_max_sizes[][CH_IM_MAX_POWER+1];
249 };
250 
251 /// Calculate what linear index this multi-index is
252 ///without the order stuff
253 template<int Dim>
255  const int & a_order)
256 {
257  int retval= 0;
258  if(a_order == 1)
259  {
260  retval = IndexedMoments<Dim, 1>::indexOf(a_index);
261  }
262  else if(a_order == 2)
263  {
264  retval = IndexedMoments<Dim, 2>::indexOf(a_index);
265  }
266  else if(a_order == 3)
267  {
268  retval = IndexedMoments<Dim, 3>::indexOf(a_index);
269  }
270  else if(a_order == 4)
271  {
272  retval = IndexedMoments<Dim, 4>::indexOf(a_index);
273  }
274  else
275  {
276  MayDay::Error("need to add more cases to getIndMomLinearIndex");
277  }
278  return retval;
279 }
280 
281 ///
282 template<int Dim>
283 const IndexTM<int,Dim>
284 getIndMomMultiIndex(const int & a_index,
285  const int & a_order)
286 {
287  IndexTM<int,Dim> retval;
288  if(a_order == 1)
289  {
290  retval = IndexedMoments<Dim,1>::getIndex(a_index);
291  }
292  else if(a_order == 2)
293  {
294  retval = IndexedMoments<Dim,2>::getIndex(a_index);
295  }
296  else if(a_order == 3)
297  {
298  retval = IndexedMoments<Dim,3>::getIndex(a_index);
299  }
300  else if(a_order == 4)
301  {
302  retval = IndexedMoments<Dim,4>::getIndex(a_index);
303  }
304  else
305  {
306  MayDay::Error("need to add more cases to getIndMomMultiIndex");
307  }
308  return retval;
309 }
310 
311 ///
312 template<int Dim>
313 int
314 getIndMomSize(const int & a_order)
315 {
316  int retval = 0;
317  if(a_order == 1)
318  {
319  retval = IndexedMoments<Dim,1>::size();
320  }
321  else if(a_order == 2)
322  {
323  retval = IndexedMoments<Dim,2>::size();
324  }
325  else if(a_order == 3)
326  {
327  retval = IndexedMoments<Dim,3>::size();
328  }
329  else if(a_order == 4)
330  {
331  retval = IndexedMoments<Dim,4>::size();
332  }
333  else
334  {
335  MayDay::Error("need to add more cases to getIndMomSize");
336  }
337  return retval;
338 }
339 
340 ///to see if all powers of p are even
341 template<int Dim>
342 bool allEven(const IndexTM<int, Dim>& a_p)
343 {
344  bool retval = true;
345  for(int idir = 0; idir < Dim; idir++)
346  {
347  if(a_p[idir]%2 != 0)
348  {
349  retval = false;
350  break;
351  }
352  }
353  return retval;
354 }
355 
356 /// computes x^p
357 template<int Dim>
358 Real
359 POW(const Real& a_x, const IndexTM<int, Dim> a_p)
360 {
361  Real retval = 1;
362  for(int idir = 0; idir < Dim; idir++)
363  {
364  for(int iexp = 0; iexp < a_p[idir]; iexp++)
365  {
366  retval *= a_x;
367  }
368  }
369  return retval;
370 }
371 ///
372 /**
373  Moments are centered at the center of the cell.
374  For each of these moments I shift them to the lowest
375  corner of the cell, where I know what the bounds of the integrals
376  is (lower bound always zero, upper bound = dx^d dx^px dx^py dx^pz
377  If the shifted moment is out of bounds, I bound it.
378  The tolerance is about verbosity. If the moment is outside the tolerance
379  then it gets included into a_bogusPowers.
380  EBMoments do not get checked for maxvals.
381 **/
382 template<int Dim, int ORDER>
383 void
385  Vector<IndexTM<int, Dim> > & a_bogusPowers,
386  const Real & a_dx,
387  const Real & a_tolerance,
388  const bool & a_ebMoment,
389  const bool & a_bindMoments);
390 
391 ///
392 /**
393  return true if all of a_index >= 0, false otherwise
394  */
395 template <int Dim>
396 bool allPositive(const IndexTM<int, Dim>& a_index);
397 
398 /**/
399 ///template specializations for debugger
400 /**
401 template< >
402 void
403 checkMoments(IndexedMoments<SpaceDim, CH_EBIS_ORDER> & a_moments,
404  Vector<IndexTM<int, SpaceDim> > & a_bogusPowers,
405  const Real & a_dx,
406  const Real & a_tolerance,
407  const bool & a_ebMoment,
408  const bool & a_bindMoments);
409 
410 
411 template< >
412 void
413 checkMoments(IndexedMoments<SpaceDim-1, CH_EBIS_ORDER> & a_moments,
414  Vector<IndexTM<int, SpaceDim-1> > & a_bogusPowers,
415  const Real & a_dx,
416  const Real & a_tolerance,
417  const bool & a_ebMoment,
418  const bool & a_bindMoments);
419 **/
420 
421 #include "NamespaceFooter.H"
422 
423 #include "IndexedMomentsImplem.H"
424 
425 #endif
const Real & operator[](int a_i) const
Definition: IndexedMoments.H:51
Real & operator[](int a_i)
Definition: IndexedMoments.H:44
static int indexOf(IndexTM< int, Dim > a_index)
Calculate what linear index this multi-index is.
Definition: IndexedMomentsImplem.H:303
void linearOut(void *const a_outbuf) const
Definition: IndexedMoments.H:100
int getIndMomSize(const int &a_order)
Definition: IndexedMoments.H:314
Definition: IndexedMoments.H:31
int getIndMomLinearIndex(const IndexTM< int, Dim > &a_index, const int &a_order)
Definition: IndexedMoments.H:254
void setRegular(const Real a_dx)
set to a regular IndexTM
Definition: IndexedMomentsImplem.H:348
static int size()
number of reals in the vector
Definition: IndexedMoments.H:78
void setToZero()
Definition: IndexedMoments.H:173
void shift(const IndexTM< Real, Dim > &a_distance)
shift moment by the input distance.
Definition: IndexedMomentsImplem.H:222
IndexedMoments()
constructor—make statics first time called
Definition: IndexedMomentsImplem.H:333
Vector< Real > m_moms
Definition: IndexedMoments.H:246
one dimensional dynamic array
Definition: Vector.H:53
static int s_size
Definition: IndexedMoments.H:229
static const Vector< IndexTM< int, Dim > > & getMonomialPowers()
monomial powers
Definition: IndexedMoments.H:134
static size_t linearSize()
for linearization
Definition: IndexedMoments.H:88
static Vector< IndexTM< int, Dim > > s_multiIndicies
Definition: IndexedMoments.H:232
int m_isRegular
Definition: IndexedMoments.H:243
bool isRegular() const
Definition: IndexedMoments.H:144
Definition: IndexTM.H:36
IndexedMoments< Dim, P > & operator*=(const Real &a_factor)
multiply each component by constant
Definition: IndexedMomentsImplem.H:198
IndexedMoments< Dim, P > & operator=(const map< IndexTM< int, Dim >, Real > &a_mapin)
for use with irregnode
Definition: IndexedMomentsImplem.H:261
static void setStatics()
Definition: IndexedMomentsImplem.H:323
void setToRegular(const Real a_dx)
Definition: IndexedMoments.H:128
void multiplyByFactorial()
multiply each entry by p!
Definition: IndexedMomentsImplem.H:459
IndexedMoments< Dim, P > & operator=(const IndexedMoments< Dim, P > &a_input)
Definition: IndexedMoments.H:153
~IndexedMoments()
Destructor.
Definition: IndexedMoments.H:39
Real & operator[](const IndexTM< int, Dim > &a_index)
Definition: IndexedMoments.H:58
double Real
Definition: REAL.H:33
void setToTruncatedMultiply(const IndexedMoments< Dim, P > &a_CA, const IndexedMoments< Dim, P > &a_CB)
Definition: IndexedMomentsImplem.H:418
Real POW(const Real &a_x, const IndexTM< int, Dim > a_p)
computes x^p
Definition: IndexedMoments.H:359
static void Error(const char *const a_msg=m_nullString, int m_exitCode=CH_DEFAULT_ERROR_CODE)
Print out message to cerr and exit with the specified exit code.
void linearIn(const void *const a_inbuf)
Definition: IndexedMoments.H:112
void spout() const
outputs the current state to pout() (a la parstream.H)
Definition: IndexedMomentsImplem.H:276
static IndexTM< int, Dim > getIndex(const int &a_linearIndex)
Definition: IndexedMoments.H:185
static void setSize()
Definition: IndexedMomentsImplem.H:383
static bool s_staticsSet
Definition: IndexedMoments.H:226
void checkMoments(IndexedMoments< Dim, ORDER > &a_moments, Vector< IndexTM< int, Dim > > &a_bogusPowers, const Real &a_dx, const Real &a_tolerance, const bool &a_ebMoment, const bool &a_bindMoments)
Definition: IndexedMomentsImplem.H:40
IndexedMoments< Dim, P > & operator+=(const IndexedMoments< Dim, P > &increment)
add compoenentwise
Definition: IndexedMomentsImplem.H:210
static void setMultiIndicies()
Definition: IndexedMomentsImplem.H:397
bool allPositive(const IndexTM< int, Dim > &a_index)
Definition: IndexedMomentsImplem.H:286
bool allEven(const IndexTM< int, Dim > &a_p)
to see if all powers of p are even
Definition: IndexedMoments.H:342
const IndexTM< int, Dim > getIndMomMultiIndex(const int &a_index, const int &a_order)
Definition: IndexedMoments.H:284
void divideByFactorial()
divides each entry by p!
Definition: IndexedMomentsImplem.H:446
static const int s_max_sizes[][CH_IM_MAX_POWER+1]
Definition: IndexedMoments.H:248
#define CH_IM_MAX_POWER
Definition: IndexedMoments.H:23
const Real & operator[](const IndexTM< int, Dim > &a_index) const
Definition: IndexedMoments.H:65
Real getMoment(const IndexTM< int, Dim > &a_mono, const map< IndexTM< int, Dim >, Real > &a_mapin) const
Definition: IndexedMomentsImplem.H:167