Chombo + EB + MF  3.2
IndexTMI.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 _INDEXTMI_H_
12 #define _INDEXTMI_H_
13 
14 #include <iostream>
15 using std::ostream;
16 using std::istream;
17 using std::ws;
18 
19 #include "MayDay.H"
20 #include "Misc.H"
21 #include "IndexTM.H"
22 #include "parstream.H"
23 #include <cmath>
24 #include <algorithm>
25 #include <functional>
26 #include <numeric>
27 #include "Metaprograms.H"
28 
29 #include "BaseNamespaceHeader.H"
30 
31 template<typename T, int N> ostream& operator<< (ostream & a_os,
32  const IndexTM<T,N> & a_p)
33 {
34  a_os << '(';
35 
36  for (int i = 0; i < N; ++i)
37  {
38  a_os << a_p[i];
39 
40  if (i < N-1)
41  {
42  a_os << ',';
43  }
44  }
45 
46  a_os << ')';
47 
48  if (a_os.fail())
49  {
50  MayDay::Abort("operator<<(ostream&,Index&) failed");
51  }
52 
53  return a_os;
54 }
55 
56 //
57 // Copied from <Utility.H>
58 //
59 #define CH_IGNORE_MAX 100000
60 
61 template<typename T, int N> istream& operator>> (istream & a_is,
62  IndexTM<T,N> & a_p)
63 {
64  a_is >> ws;
65 
66  char c;
67  a_is >> c;
68  a_is.putback(c);
69 
70  if (c == '(')
71  {
72  a_is.ignore(CH_IGNORE_MAX, '(') >> a_p[0];
73 
74  for (int i = 1; i < N; ++i)
75  {
76  a_is.ignore(CH_IGNORE_MAX, ',') >> a_p[i];
77  }
78 
79  a_is.ignore(CH_IGNORE_MAX, ')');
80  }
81  else if (c == '<')
82  {
83  a_is.ignore(CH_IGNORE_MAX, '<') >> a_p[0];
84 
85  for (int i = 1; i < N; ++i)
86  {
87  a_is.ignore(CH_IGNORE_MAX, ',') >> a_p[1];
88  }
89  a_is.ignore(CH_IGNORE_MAX, '>');
90  }
91  else
92  {
93  MayDay::Abort("operator>>(istream&,Index&): expected \'(\' or \'<\'");
94  }
95 
96  if (a_is.fail())
97  {
98  MayDay::Abort("operator>>(istream&,Index&) failed");
99  }
100 
101  return a_is;
102 }
103 
104 template<typename T, int N> void IndexTM<T,N>::printOn (ostream & a_os) const
105 {
106  a_os << "Index: " << *this << "\n";
107 }
108 
109 template<typename T, int N> void IndexTM<T,N>::p() const
110 {
111  pout() << *this << "\n";
112 }
113 
114 template<typename T, int N> void IndexTM<T,N>::dumpOn (ostream & a_os) const
115 {
116  a_os << "Index " << *this << "\n";
117 }
118 
119 //
120 // Static object initialization.
121 //
122 /*
123 template<typename T, int N> int IndexTM<T,N>::InitStatics()
124 {
125  IndexTM<T,N>& pz = const_cast<IndexTM<T,N>&>(IndexTM<T,N>::Zero);
126 
127  pz.setAll(0);
128 
129  IndexTM<T,N>& pu = const_cast<IndexTM<T,N>&>(IndexTM<T,N>::Unit);
130 
131  pu.setAll(1);
132 
133  // No danger of Index::Zero and Unit not having been allocated, as ARM section
134  // 3.4 says "The initialization of nonlocal static objects in a translation unit
135  // is done before the first use of any function...defined in that translation
136  // unit."
137  //
138  // Had to go through the const_cast stuff because it's nice to be able to declare
139  // Index::Zero and Index::Unit as const.
140 
141  return 0; // arbitrary
142 }
143 */
144 
145 //
146 // Inlines.
147 //
148 template<typename T, int N> inline IndexTM<T,N>::IndexTM(const T *a_a)
149  : GenericArithmeticable< T, IndexTM<T,N> >(this)
150 {
151  memcpy(m_vect, a_a, N*sizeof(T));
152 }
153 
154 template<typename T, int N> inline IndexTM<T,N>::IndexTM(const char* a_reference)
155  : GenericArithmeticable< T, IndexTM<T,N> >(this)
156 {
157  if (a_reference[0] == 48) setAll(0);
158  else if (a_reference[0] == 49) setAll(1);
159  else
160  MayDay::Error("unknown static initialization for IndexTM");
161 }
162 
163 template<typename T, int N> inline IndexTM<T,N>::IndexTM(const IndexTM<T,N> & a_iv)
164  : GenericArithmeticable< T, IndexTM<T,N> >(this)
165 {
166  memcpy(m_vect, a_iv.m_vect, N*sizeof(T));
167 }
168 
169 template<typename T, int N> inline IndexTM<T,N>& IndexTM<T,N>::operator=(const IndexTM<T,N> & a_iv)
170 {
171  memcpy(m_vect, a_iv.m_vect, N*sizeof(T));
172  return *this;
173 }
174 
175 template<typename T, int N> inline T& IndexTM<T,N>::operator[] (int a_i)
176 {
177  CH_assert(a_i >= 0 && a_i < N);
178  return m_vect[a_i];
179 }
180 
181 template<typename T, int N> inline T IndexTM<T,N>::operator[] (int a_i) const
182 {
183  CH_assert(a_i >= 0 && a_i < N);
184  return m_vect[a_i];
185 }
186 
187 template<typename T, int N> inline void IndexTM<T,N>::setVal(int a_i,
188  T a_val)
189 {
190  CH_assert(a_i >= 0 && a_i < N);
191  m_vect[a_i] = a_val;
192 }
193 
194 template<typename T, int N> inline void IndexTM<T,N>::setAll(T a_val)
195 {
196  for (int i = 0; i < N; i++)
197  {
198  m_vect[i] = a_val;
199  }
200 }
201 
202 template<typename T, int N> inline const T* IndexTM<T,N>::dataPtr() const
203 {
204  return m_vect;
205 }
206 
207 template<typename T, int N> inline T* IndexTM<T,N>::dataPtr()
208 {
209  return m_vect;
210 }
211 
212 template<typename T, int N> inline const T* IndexTM<T,N>::getVect() const
213 {
214  return m_vect;
215 }
216 
217 template<typename T, int N> inline bool IndexTM<T,N>::operator== (const IndexTM & a_p) const
218 {
219  return Metaprograms::pointwiseCompare<N,T,std::equal_to<T> >(m_vect,a_p.m_vect);
220 }
221 
222 template<typename T, int N> inline bool IndexTM<T,N>::operator!= (const IndexTM & a_p) const
223 {
224  return !(operator==(a_p));
225 }
226 
227 template<typename T, int N> inline bool IndexTM<T,N>::lexLT(const IndexTM & a_s) const
228 {
229  return Metaprograms::LexLT<N,T>()(m_vect,a_s.m_vect);
230 }
231 
232 template<typename T, int N> inline bool IndexTM<T,N>::lexGT(const IndexTM & a_s) const
233 {
234  return ! Metaprograms::LexLT<N,T>()(m_vect,a_s.m_vect);
235 }
236 
237 template<typename T, int N> inline IndexTM<T,N> IndexTM<T,N>::operator+ () const
238 {
239  return IndexTM<T,N>(*this);
240 }
241 
242 template<typename T, int N> inline IndexTM<T,N> IndexTM<T,N>::operator- () const
243 {
244  IndexTM<T,N> result(*this);
245  for (int i = 0; i < N; ++i)
246  {
247  result.m_vect[i] *= -1;
248  }
249  return result;
250 }
251 
252 template<typename T, int N> inline T IndexTM<T,N>::dotProduct(const IndexTM<T,N> & a_rhs ) const
253 {
254  return Metaprograms::InnerProduct<N,T,T,
255  std::plus<T>,
256  std::multiplies<T> >()(m_vect,
257  a_rhs.m_vect);
258 }
259 
260 template<typename T, int N> inline T IndexTM<T,N>::sum() const
261 {
263 }
264 
265 template<typename T, int N> inline T IndexTM<T,N>::product () const
266 {
268 }
269 
270 template<typename T, int N> template<typename OP> bool IndexTM<T,N>::operatorCompare(const IndexTM<T,N> & a_p,
271  const OP & a_op) const
272 {
273  return Metaprograms::pointwiseCompare<N,T,OP>(m_vect,a_p.m_vect);
274 }
275 
276 template<typename T, int N> template<typename OP> inline IndexTM<T,N>& IndexTM<T,N>::operatorOpEquals(const T & a_s,
277  const OP & a_op)
278 {
280  return *this;
281 }
282 
283 template<typename T, int N> template<typename OP> inline IndexTM<T,N>& IndexTM<T,N>::operatorOpEquals(const IndexTM<T,N> & a_p,
284  const OP & a_op)
285 {
287  return *this;
288 }
289 
290 template<typename T, int N> inline IndexTM<T,N>& IndexTM<T,N>::reciprocal()
291 {
292  std::transform(m_vect, m_vect+N, m_vect, bind1st(std::divides<T>(),T(1)));
293  return *this;
294 }
295 
296 template<typename T> static bool abscompare(const T & a_a,
297  const T & a_b)
298 {
299  return Abs(a_a) < Abs(a_b);
300 }
301 template<typename T, int N> inline int IndexTM<T,N>::minDir(bool a_doAbs) const
302 {
303  if (a_doAbs)
304  {
305  return std::min_element(m_vect, m_vect+N, std::ptr_fun(abscompare<T>)) - m_vect;
306  }
307  else
308  {
309  return std::min_element(m_vect, m_vect+N) - m_vect;
310  }
311 }
312 
313 template<typename T, int N> inline int IndexTM<T,N>::maxDir(bool a_doAbs) const
314 {
315  if (a_doAbs)
316  {
317  return std::max_element(m_vect, m_vect+N, std::ptr_fun(abscompare<T>)) - m_vect;
318  }
319  else
320  {
321  return std::max_element(m_vect, m_vect+N) - m_vect;
322  }
323 }
324 
325 template<typename T> static T ourmin(T a_a,
326  T a_b)
327 {
328  return ((a_a < a_b) ? a_a : a_b);
329 }
330 
331 template<typename T> static T ourmax(T a_a,
332  T a_b)
333 {
334  return (a_a > a_b ? a_a : a_b);
335 }
336 
337 template<typename T, int N> inline IndexTM<T,N>& IndexTM<T,N>::min(const IndexTM<T,N> & a_p)
338 {
339  std::transform(m_vect, m_vect+N, a_p.m_vect, m_vect,
340  std::ptr_fun(ourmin<T>));
341  return *this;
342 }
343 
344 template<typename T, int N> inline IndexTM<T,N>& IndexTM<T,N>::max(const IndexTM<T,N> & a_p)
345 {
346  std::transform(m_vect, m_vect+N, a_p.m_vect, m_vect,
347  std::ptr_fun(ourmax<T>));
348  return *this;
349 }
350 
351 template<typename T, int N> inline IndexTM<T,N>& IndexTM<T,N>::scale(T a_s)
352 {
353  return (*this) *= a_s;
354 }
355 
356 template<typename T, int N> inline IndexTM<T,N>& IndexTM<T,N>::reflect (T a_refIx,
357  int a_idir)
358 {
359  CH_assert(a_idir >= 0 && a_idir < N);
360  m_vect[a_idir] = -m_vect[a_idir] + 2*a_refIx;
361  return *this;
362 }
363 
364 template<typename T, int N> inline IndexTM<T,N>& IndexTM<T,N>::shift(int a_coord,
365  T a_s)
366 {
367  CH_assert(a_coord >= 0 && a_coord < N);
368  m_vect[a_coord] += a_s;
369  return *this;
370 }
371 
372 template<typename T, int N> inline IndexTM<T,N>& IndexTM<T,N>::shift(const IndexTM<T,N> & a_iv)
373 {
374  return (*this) += a_iv;
375 }
376 
377 template<typename T, int N> inline IndexTM<T,N>& IndexTM<T,N>::diagShift(T a_s)
378 {
379  return (*this) += a_s;
380 }
381 
382 template<typename T, int N> inline IndexTM<T,N> scale (const IndexTM<T,N> & a_p,
383  T a_s)
384 {
385  return a_p * a_s;
386 }
387 
388 template<typename T, int N> inline IndexTM<T,N> diagShift(const IndexTM<T,N> & a_p,
389  T a_s)
390 {
391  return a_p + a_s;
392 }
393 
394 template<typename T, int N> inline IndexTM<T,N> min(const IndexTM<T,N> & a_p1,
395  const IndexTM<T,N> & a_p2)
396 {
397  IndexTM<T,N> result(a_p1);
398  return result.min(a_p2);
399 }
400 
401 template<typename T, int N> inline IndexTM<T,N> max(const IndexTM<T,N> & a_p1,
402  const IndexTM<T,N> & a_p2)
403 {
404  IndexTM<T,N> result(a_p1);
405  return result.max(a_p2);
406 }
407 
408 template<typename T, int N> inline IndexTM<T,N> BASISV_TM(int a_dir)
409 {
410  CH_assert(a_dir >= 0 && a_dir < N);
412  tmp.dataPtr()[a_dir] = T(1);
413  return tmp;
414 }
415 
416 template<typename T, int N> inline IndexTM<T,N> reflect(const IndexTM<T,N> & a_a,
417  T a_refIx,
418  int a_idir)
419 {
420  IndexTM<T,N> result(a_a);
421  return result.reflect(a_refIx, a_idir);
422 }
423 
424 template<typename T> static T ourcoarsen(T a_a,
425  T a_b)
426 {
427  return (a_a < 0 ? T(-Abs(a_a+1))/a_b-1 : a_a/a_b);
428 }
429 
430 template<typename T, int N> inline IndexTM<T,N> coarsen(const IndexTM<T,N> & a_p,
431  T a_s)
432 {
433  IndexTM<T,N> result(a_p);
434  return result.coarsen(a_s);
435 }
436 
437 template<typename T, int N> inline IndexTM<T,N> coarsen(const IndexTM<T,N> & a_p1,
438  const IndexTM<T,N> & a_p2)
439 {
440  IndexTM<T,N> result(a_p1);
441  return result.coarsen(a_p2);
442 }
443 
444 /*
445 template<typename T, int N> IndexTM<T,N>::operator IndexTM<Real,N>()
446 {
447  IndexTM<Real,N> result;
448  for (int i = 0; i < N; ++i)
449  {
450  result.dataPtr()[i] = Real(m_vect[i]);
451  }
452  return result;
453 }
454 */
455 
456 template<typename T, int N> inline IndexTM<T,N>& IndexTM<T,N>::coarsen(T a_s)
457 {
458  CH_assert(a_s > 0);
459  std::transform(m_vect, m_vect+N, m_vect,
460  std::bind2nd(std::ptr_fun(ourcoarsen<T>),a_s));
461  return *this;
462 }
463 
464 template<typename T, int N> inline IndexTM<T,N>& IndexTM<T,N>::coarsen(const IndexTM<T,N> & a_p)
465 {
466  CH_assert(a_p > (IndexTM<T,N>::Zero));
467  std::transform(m_vect, m_vect+N, a_p.m_vect, m_vect, std::ptr_fun(ourcoarsen<T>));
468  return *this;
469 }
470 
471 template<typename T, int N> void IndexTM<T,N>::linearIn(const void* a_inBuf)
472 {
473  memcpy(m_vect, (T*)a_inBuf, N*sizeof(T));
474 }
475 
476 template<typename T, int N> void IndexTM<T,N>::linearOut(void* a_outBuf) const
477 {
478  memcpy((T*)a_outBuf, m_vect, N*sizeof(T));
479 }
480 
481 template<typename T, int N> IndexTM<T,N>::IndexTM(T a_i)
482  : GenericArithmeticable< T, IndexTM<T,N> >(this)
483 {
484  static_assert(N == 1,"N==1");
485  m_vect[0] = a_i;
486 }
487 
488 template<typename T, int N> IndexTM<T,N>::IndexTM(T a_i,
489  T a_j)
490  : GenericArithmeticable< T, IndexTM<T,N> >(this)
491 {
492  static_assert(N == 2, "N==2");
493  m_vect[0] = a_i;
494  m_vect[1] = a_j;
495 }
496 
497 template<typename T, int N> IndexTM<T,N>::IndexTM(T a_i,
498  T a_j,
499  T a_k)
500  : GenericArithmeticable< T, IndexTM<T,N> >(this)
501 {
502  static_assert(N == 3, "N==3");
503  m_vect[0] = a_i;
504  m_vect[1] = a_j;
505  m_vect[2] = a_k;
506 }
507 
508 template<typename T, int N> IndexTM<T,N>::IndexTM(T a_i,
509  T a_j,
510  T a_k,
511  T a_l)
512  : GenericArithmeticable< T, IndexTM<T,N> >(this)
513 {
514  static_assert(N == 4, "N==4");
515  m_vect[0] = a_i;
516  m_vect[1] = a_j;
517  m_vect[2] = a_k;
518  m_vect[3] = a_l;
519 }
520 
521 template<typename T, int N> IndexTM<T,N>::IndexTM(T a_i,
522  T a_j,
523  T a_k,
524  T a_l,
525  T a_m)
526  : GenericArithmeticable< T, IndexTM<T,N> >(this)
527 {
528  static_assert(N == 5, "N==5");
529  m_vect[0] = a_i;
530  m_vect[1] = a_j;
531  m_vect[2] = a_k;
532  m_vect[3] = a_l;
533  m_vect[4] = a_m;
534 }
535 
536 template<typename T, int N> IndexTM<T,N>::IndexTM(T a_i,
537  T a_j,
538  T a_k,
539  T a_l,
540  T a_m,
541  T a_n)
542  : GenericArithmeticable< T, IndexTM<T,N> >(this)
543 {
544  static_assert(N == 6, "N==6");
545  m_vect[0] = a_i;
546  m_vect[1] = a_j;
547  m_vect[2] = a_k;
548  m_vect[3] = a_l;
549  m_vect[4] = a_m;
550  m_vect[5] = a_n;
551 }
552 
553 /*
554 template<typename T, int N> const IndexTM<T,N> IndexTM<T,N>::Zero;
555 template<typename T, int N> const IndexTM<T,N> IndexTM<T,N>::Unit;
556 
557 static int s_dummyForIntVectCpp1 (IndexTM<int ,1>::InitStatics());
558 static int s_dummyForIntVectCpp2 (IndexTM<int ,2>::InitStatics());
559 static int s_dummyForIntVectCpp3 (IndexTM<int ,3>::InitStatics());
560 static int s_dummyForIntVectCpp4 (IndexTM<int ,4>::InitStatics());
561 static int s_dummyForIntVectCpp5 (IndexTM<int ,5>::InitStatics());
562 static int s_dummyForIntVectCpp6 (IndexTM<int ,6>::InitStatics());
563 static int s_dummyForRealVectCpp1(IndexTM<Real,1>::InitStatics());
564 static int s_dummyForRealVectCpp2(IndexTM<Real,2>::InitStatics());
565 static int s_dummyForRealVectCpp3(IndexTM<Real,3>::InitStatics());
566 static int s_dummyForRealVectCpp4(IndexTM<Real,4>::InitStatics());
567 static int s_dummyForRealVectCpp5(IndexTM<Real,5>::InitStatics());
568 static int s_dummyForRealVectCpp6(IndexTM<Real,6>::InitStatics());
569 */
570 
571 // If Index::Zero and Index::Unit were pointers, we wouldn't need this extra
572 // static int. But they're objects, and the danger is that the initializations
573 // right above here ("Index Index::Zero;" and "Index Index::Unit;") are hit
574 // after the last call to Index::InitStatics, and in that case the
575 // Index::Index() constructor could redefine Zero and Unit. In fact, the way
576 // things stand now, nothing bad would happen, because the Index::Index()
577 // constructor doesn't assign anything to any of the data members. But we don't
578 // want to count on that always being the case.
579 
580 #include "BaseNamespaceFooter.H"
581 
582 #endif // include guard
const T * dataPtr() const
Definition: IndexTMI.H:202
std::ostream & pout()
Use this in place of std::cout for program output.
void dumpOn(std::ostream &a_os) const
Definition: IndexTMI.H:114
istream & operator>>(istream &a_is, IndexTM< T, N > &a_p)
Definition: IndexTMI.H:61
Definition: Metaprograms.H:152
IndexTM< T, N > & operatorOpEquals(const IndexTM< T, N > &a_p, const OP &a_op)
Definition: IndexTMI.H:283
#define CH_assert(cond)
Definition: CHArray.H:37
IndexTM()
Definition: IndexTM.H:93
void setAll(T a_val)
Definition: IndexTMI.H:194
static bool abscompare(const T &a_a, const T &a_b)
Definition: IndexTMI.H:296
bool operator!=(const IndexTM &a_p) const
Definition: IndexTMI.H:222
IndexTM operator-() const
Definition: IndexTMI.H:242
static T ourmin(T a_a, T a_b)
Definition: IndexTMI.H:325
Definition: GenericArithmetic.H:54
T dotProduct(const IndexTM &a_rhs) const
Definition: IndexTMI.H:252
IndexTM< T, N > BASISV_TM(int a_dir)
Definition: IndexTMI.H:408
IndexTM & reciprocal()
Definition: IndexTMI.H:290
bool lexLT(const IndexTM &a_s) const
Definition: IndexTMI.H:227
Definition: IndexTM.H:36
IndexTM & diagShift(T a_s)
Definition: IndexTMI.H:377
void setVal(int a_i, T a_val)
Definition: IndexTMI.H:187
T sum() const
Definition: IndexTMI.H:260
T & operator[](int a_i)
Definition: IndexTMI.H:175
static T ourcoarsen(T a_a, T a_b)
Definition: IndexTMI.H:424
IndexTM & reflect(T a_refIx, int a_idir)
Definition: IndexTMI.H:356
IndexTM operator+() const
Definition: IndexTMI.H:237
bool lexGT(const IndexTM &a_s) const
Definition: IndexTMI.H:232
IndexTM & coarsen(const IndexTM &a_p)
Definition: IndexTMI.H:464
IndexTM & scale(T a_s)
Definition: IndexTMI.H:351
IndexTM & min(const IndexTM &a_p)
Definition: IndexTMI.H:337
T Abs(const T &a_a)
Definition: Misc.H:53
bool operator==(const IndexTM &a_p) const
Definition: IndexTMI.H:217
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 linearOut(void *a_outBuf) const
Definition: IndexTMI.H:476
T product() const
Definition: IndexTMI.H:265
bool operatorCompare(const IndexTM< T, N > &a_p, const OP &a_op) const
Definition: IndexTMI.H:270
#define CH_IGNORE_MAX
Definition: IndexTMI.H:59
Definition: Metaprograms.H:122
void p() const
Definition: IndexTMI.H:109
ostream & operator<<(ostream &a_os, const IndexTM< T, N > &a_p)
Definition: IndexTMI.H:31
const T * getVect() const
Definition: IndexTMI.H:212
Definition: Metaprograms.H:39
Definition: Metaprograms.H:91
void const char const int const int * N
Definition: Lapack.H:83
IndexTM & max(const IndexTM &a_p)
Definition: IndexTMI.H:344
IndexTM & shift(int a_coord, T a_s)
Definition: IndexTMI.H:364
void linearIn(const void *a_inBuf)
Definition: IndexTMI.H:471
void printOn(std::ostream &a_os) const
Definition: IndexTMI.H:104
int minDir(bool a_doAbs) const
Definition: IndexTMI.H:301
IndexTM & operator=(const IndexTM &a_rhs)
Definition: IndexTMI.H:169
static T ourmax(T a_a, T a_b)
Definition: IndexTMI.H:331
T m_vect[N]
Definition: IndexTM.H:522
static void Abort(const char *const a_msg=m_nullString)
Print out message to cerr and exit via abort() (if serial) or MPI_Abort() (if parallel).
int maxDir(bool a_doAbs) const
Definition: IndexTMI.H:313