Chombo + EB  3.2
IVSFABI.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 _IVSFABI_H_
12 #define _IVSFABI_H_
13 
14 #include "MayDay.H"
15 #include "IntVectSet.H"
16 #include "parstream.H"
17 #include "SPMD.H"
18 #include "DebugOut.H"
19 #include "SPMD.H"
20 #include "NamespaceHeader.H"
21 
22 template <class T>
23 bool IVSFAB<T>::s_verbose = false;
24 
25 template <class T>
26 void
27 IVSFAB<T>::setVerbose(bool a_verbose)
28 {
29  s_verbose = a_verbose;
30 }
31 /*************/
32 template <class T> inline
33 const IntVectSet&
34 IVSFAB<T>::getIVS() const
35 {
36  return m_ivs;
37 }
38 /******************/
39 template <class T> inline
40 IVSFAB<T>::IVSFAB():m_dataPtr(0),m_nComp(0),m_nIvs(0),m_loVect(IntVect::Unit), m_hiVect(IntVect::Zero), m_isDefined(false)
41 {
42 
43 }
44 /******************/
45 template <class T> inline
46 IVSFAB<T>::~IVSFAB()
47 {
48  clear();
49 }
50 /******************/
51 template <class T> inline
52 IVSFAB<T>::IVSFAB(const IntVectSet& a_ivsin,
53  const int& a_nvarin)
54 {
55  define( a_ivsin, a_nvarin );
56 }
57 
58 
59 /******************/
60 template<class T> inline
61 void IVSFAB<T>::define(const IntVectSet& a_ivsin,
62  const int& a_nvarin )
63 {
64  m_dataPtr = NULL;
65  m_nComp = a_nvarin;
66  m_nIvs = 0;
67  m_loVect = IntVect::Zero;
68  m_hiVect = IntVect::Unit;
69  m_ivmap.define( a_ivsin.minBox(), a_nvarin );
70  m_ivs = a_ivsin;
71 
72  CH_assert(a_nvarin > 0);
73  if (!a_ivsin.isEmpty())
74  {
75  IVSIterator ivsit(a_ivsin);
76  for (ivsit.reset(); ivsit.ok(); ++ivsit)
77  {
78  const IntVect& iv = ivsit();
79  m_ivmap(iv,0) = m_nIvs;
80  m_nIvs++;
81  }
82 
83  if ( (m_nIvs > 0) && m_nComp > 0)
84  {
85  m_dataPtr = new T[m_nComp*m_nIvs];
86  }
87  }
88  if (m_nIvs > 0)
89  {
90  //set up face low and high vectors
91  for (int dir = 0; dir < CH_SPACEDIM; ++dir)
92  {
93  m_loVect[dir] = 1;
94  m_hiVect[dir] = 1;
95  }
96  m_hiVect[0] = m_nIvs;
97 
98  }
99  else
100  {
101  m_loVect = IntVect::Unit;
102  m_hiVect = IntVect::Zero;
103  m_dataPtr = NULL;
104  }
105  m_isDefined=true;
106 
107 }
108 /******************/
109 template <class T> inline
110 void
111 IVSFAB<T>::setVal(const T& a_value)
112 {
113  CH_assert(isDefined());
114  for (int ivec = 0; ivec < m_nIvs*m_nComp; ivec++)
115  m_dataPtr[ivec] = a_value;
116 }
117 /******************/
118 template <class T> inline
119 void
120 IVSFAB<T>::copy(const Box& a_fromBox,
121  const Interval& a_dstInterval,
122  const Box& a_toBox,
123  const IVSFAB<T>& a_src,
124  const Interval& a_srcInterval)
125 {
126  IntVect shift = a_toBox.smallEnd() - a_fromBox.smallEnd();
127  IntVect shift2 = a_toBox.bigEnd() - a_fromBox.bigEnd();
128  CH_assert(shift==shift2);
129  CH_assert(isDefined());
130  CH_assert(a_src.isDefined());
131  CH_assert(a_srcInterval.size() == a_dstInterval.size());
132  CH_assert(a_dstInterval.begin() >= 0);
133  CH_assert(a_srcInterval.begin() >= 0);
134  CH_assert(a_dstInterval.end() < m_nComp);
135  CH_assert(a_srcInterval.end() < a_src.m_nComp);
136 
137  if ((!m_ivs.isEmpty()) && (!a_src.m_ivs.isEmpty()))
138  {
139  // my points in toBox
140  IntVectSet ivsIntersect = m_ivs & a_toBox;
141  // his points in fromBox
142  IntVectSet set2 = a_src. m_ivs & a_fromBox;
143  // his points mapped to my Box
144  set2.shift(shift);
145  // overlap
146  ivsIntersect &= set2;
147 
148  IVSFAB<T>& thisFAB = *this;
149  int compSize = a_srcInterval.size();
150  for (IVSIterator ivit(ivsIntersect); ivit.ok(); ++ivit)
151  {
152  const IntVect& iv = ivit();
153  IntVect iv2 = iv - shift;
154  for (int icomp = 0; icomp < compSize; icomp++)
155  {
156  int isrccomp = a_srcInterval.begin() + icomp;
157  int idstcomp = a_dstInterval.begin() + icomp;
158  thisFAB(iv, idstcomp) = a_src(iv2, isrccomp);
159  }
160  }
161  }
162 }
163 /******************/
164 //template specializations for real
165 template <> inline
166 int IVSFAB<Real>::size(const Box& a_region,
167  const Interval& a_comps) const
168 {
169  // size of actual data plus size of count ..
170  // specialization for T=Real.
171  IntVectSet subset = m_ivs & a_region;
172  int numberThings = subset.numPts();
173  int retval = numberThings*sizeof(Real)*a_comps.size() + sizeof(int);
174  retval += SpaceDim*sizeof(int)*numberThings; // for intvects
175  retval += 2*SpaceDim*sizeof(Real); // box done as reals for some reason
176  return retval;
177 }
178 
179 /********************/
180 template <> inline void IVSFAB<Real>::linearIn(void* a_buf, const Box& a_region, const Interval& a_comps)
181 {
182  // message includes count
183  IntVectSet subset = m_ivs & a_region;
184 
185  // input the box
186  IntVect ibl, ibh;
187  Real* rptr = (Real*)a_buf;
188  for ( int idir=0; idir<SpaceDim; idir++ ) {
189  ibl[idir] = *rptr++;
190  ibh[idir] = *rptr++;
191  }
192 
193  // compute the offset (for periodic case)
194  IntVect shift = a_region.smallEnd() - ibl;
195  IntVect shift2 = a_region.bigEnd() - ibh;
196  CH_assert(shift==shift2);
197 
198  int* iptr = (int*)rptr;
199  int count_received = *iptr++;
200  rptr = (Real*)iptr;
201  for ( int i=0; i<count_received; i++ )
202  {
203  IntVect iv2;
204  iptr = (int*)rptr;
205  for ( int idir=0; idir<SpaceDim; idir++ ) iv2[idir] = *iptr++;
206  rptr = (Real*)iptr;
207 
208  // perform shift
209  iv2 += shift;
210 
211  for ( int c=a_comps.begin(); c<=a_comps.end(); c++ )
212  {
213  int iloc = getIndex(iv2, c);
214  m_dataPtr[iloc] = *rptr++;
215  }
216  }
217 }
218 
219 /********************/
220 template <> inline void IVSFAB<Real>::linearOut(void* a_buf,
221  const Box& a_region,
222  const Interval& a_comps) const
223 {
224  IntVectSet subset = m_ivs & a_region;
225 
226  Real* rptr = (Real*)a_buf;
227  for ( int idir=0; idir<SpaceDim; idir++ ) {
228  *rptr++ = a_region.smallEnd(idir);
229  *rptr++ = a_region.bigEnd(idir);
230  }
231 
232  int* iptr = (int*)rptr;
233  int count = subset.numPts();
234  *iptr++ = count;
235  rptr = (Real*)iptr;
236 
237  for ( IVSIterator ivsit(subset); ivsit.ok(); ++ivsit )
238  {
239  const IntVect& iv = ivsit();
240  iptr = (int*)rptr;
241  for ( int idir=0; idir<SpaceDim; idir++ ) *iptr++ = iv[idir];
242  rptr = (Real*)iptr;
243  for ( int c=a_comps.begin(); c<=a_comps.end(); c++ )
244  {
245  int iloc = getIndex(iv, c);
246  *rptr++ = m_dataPtr[iloc];
247  }
248  }
249  CH_assert( ((char*)rptr) - ((char*)a_buf) == size(a_region,a_comps) );
250 }
251 
252 
253 /******************/
254 //template specializations for int
255 template <> inline
256 int IVSFAB<int>::size(const Box& a_region,
257  const Interval& a_comps) const
258 {
259  // size of actual data plus size of count ..
260  // specialization for T=int
261  IntVectSet subset = m_ivs & a_region;
262  int numberThings = subset.numPts();
263  int retval = numberThings*sizeof(int)*a_comps.size() + sizeof(int);
264  retval += SpaceDim*sizeof(int)*numberThings; // for intvects
265  retval += 2*SpaceDim*sizeof(int); // box
266  return retval;
267 }
268 
269 
270 /********************/
271 template <> inline void IVSFAB<int>::linearOut(void* a_buf,
272  const Box& a_region,
273  const Interval& a_comps) const
274 {
275  IntVectSet subset = m_ivs & a_region;
276 
277  int* rptr = (int*)a_buf;
278  for ( int idir=0; idir<SpaceDim; idir++ )
279  {
280  *rptr++ = a_region.smallEnd(idir);
281  *rptr++ = a_region.bigEnd(idir);
282  }
283 
284  int* iptr = (int*)rptr;
285  int count = subset.numPts();
286  *iptr++ = count;
287  rptr = (int*)iptr;
288 
289  for ( IVSIterator ivsit(subset); ivsit.ok(); ++ivsit )
290  {
291  const IntVect& iv = ivsit();
292  iptr = (int*)rptr;
293  for ( int idir=0; idir<SpaceDim; idir++ ) *iptr++ = iv[idir];
294  rptr = (int*)iptr;
295  for ( int c=a_comps.begin(); c<=a_comps.end(); c++ )
296  {
297  int iloc = getIndex(iv, c);
298  *rptr++ = m_dataPtr[iloc];
299  }
300  }
301  // CH_assert( ((char*)rptr) - ((char*)a_buf) == size(a_region,a_comps) );
302 }
303 /********************/
304 template <> inline void IVSFAB<int>::linearIn(void* a_buf, const Box& a_region, const Interval& a_comps)
305 {
306  // message includes count
307  IntVectSet subset = m_ivs & a_region;
308 
309  // input the box
310  IntVect ibl, ibh;
311  int* rptr = (int*)a_buf;
312  for ( int idir=0; idir<SpaceDim; idir++ )
313  {
314  ibl[idir] = *rptr++;
315  ibh[idir] = *rptr++;
316  }
317 
318  // compute the offset (for periodic case)
319  IntVect shift = a_region.smallEnd() - ibl;
320  IntVect shift2 = a_region.bigEnd() - ibh;
321  CH_assert(shift==shift2);
322 
323  int* iptr = (int*)rptr;
324  int count_received = *iptr++;
325  rptr = (int*)iptr;
326  for ( int i=0; i<count_received; i++ )
327  {
328  IntVect iv2;
329  iptr = (int*)rptr;
330  for ( int idir=0; idir<SpaceDim; idir++ ) iv2[idir] = *iptr++;
331  rptr = (int*)iptr;
332 
333  // perform shift
334  iv2 += shift;
335 
336  for ( int c=a_comps.begin(); c<=a_comps.end(); c++ )
337  {
338  int iloc = getIndex(iv2, c);
339  m_dataPtr[iloc] = *rptr++;
340  }
341  }
342 }
343 ////////////////////////////////////////
344 template <class T> inline
345 int IVSFAB<T>::size(const Box& a_region,
346  const Interval& a_comps) const
347 {
348  IntVectSet subset = m_ivs & a_region;
349  int numberThings = subset.numPts();
350 
351 // can't do this for objects like vectors ...
352 // T dummy;
353 // int sizeofThing = dummy.linearSize();
354 // int retval = numberThings*sizeofThing + sizeof(int);
355 
356  int retval=sizeof(int); // count
357  for( int icomp=a_comps.begin(); icomp<=a_comps.end(); icomp++ ) {
358  for( IVSIterator ivsit(subset); ivsit.ok(); ++ivsit ) {
359  const IntVect& iv = ivsit();
360  int iloc = getIndex(iv,icomp);
361  retval += m_dataPtr[iloc].linearSize(); // size of each object
362  }
363  }
364  retval += SpaceDim*sizeof(int)*numberThings; // for intvects
365  retval += 2*SpaceDim*sizeof(int); // box
366  return retval;
367 }
368 ////////////////////
369 template <class T> inline
370 void IVSFAB<T>::linearOut(void* a_buf,
371  const Box& a_region,
372  const Interval& a_comps) const
373 {
374  // Note: it WAS assumed that there is one component.
375  // and all "T" have equal size.
376  // message includes count
377  IntVectSet subset = m_ivs & a_region;
378 
379 
380  Real* rptr = (Real*)a_buf;
381  for ( int idir=0; idir<SpaceDim; idir++ ) {
382  *rptr++ = a_region.smallEnd(idir);
383  *rptr++ = a_region.bigEnd(idir);
384  }
385 
386  int* iptr = (int*)rptr;
387 
388  int count = subset.numPts();
389  *iptr++ = count;
390 // cannot count on T having constant size
391 // T dummy;
392 // int Tsz = dummy.linearSize();
393  char* cptr = (char*)iptr;
394  for ( IVSIterator ivsit(subset); ivsit.ok(); ++ivsit )
395  {
396  const IntVect& iv = ivsit();
397 
398  iptr = (int*)cptr;
399  for ( int idir=0; idir<SpaceDim; idir++ ) *iptr++ = iv[idir];
400  cptr = (char*)iptr;
401 
402  for( int c=a_comps.begin(); c<=a_comps.end(); c++ ) {
403  int iloc = getIndex(iv, c); // assume 1 component ...
404  m_dataPtr[iloc].linearOut( (void*)cptr );
405  cptr += m_dataPtr[iloc].linearSize();
406  }
407  }
408 if( cptr - ((char*)a_buf) != size(a_region,a_comps) ) {
409  pout() << " region " << a_region << endl;
410  pout() << " comps " << a_comps.begin() << " to " << a_comps.end() << endl;
411  pout() << " size " << size(a_region,a_comps) << endl;
412  pout() << " ptrdiff " << (cptr - ((char*)a_buf)) << endl;
413 }
414  CH_assert( cptr - ((char*)a_buf) == size(a_region,a_comps) );
415 }
416 
417 template <class T> inline
418 void IVSFAB<T>::linearIn(void* a_buf, const Box& a_region, const Interval& a_comps)
419 {
420  // Note: it is assumed that there is one component.
421  // and all "T" have same size
422  // message includes count
423  IntVectSet subset = m_ivs & a_region;
424 
425  // input the box
426  IntVect ibl, ibh;
427  Real* rptr = (Real*)a_buf;
428  for ( int idir=0; idir<SpaceDim; idir++ ) {
429  ibl[idir] = *rptr++;
430  ibh[idir] = *rptr++;
431  }
432 
433  // compute the offset (for periodic case)
434  IntVect shift = a_region.smallEnd() - ibl;
435  IntVect shift2 = a_region.bigEnd() - ibh;
436  CH_assert(shift==shift2);
437 
438  int* iptr = (int*)rptr;
439 
440  int count_received = *iptr++;
441 // cannot count on T having constant size
442 // T dummy;
443 // int Tsz = dummy.linearSize();
444  char* cptr = (char*)iptr;
445  for ( int i=0; i<count_received; i++ )
446  {
447  IntVect iv2;
448  iptr = (int*)cptr;
449  for ( int idir=0; idir<SpaceDim; idir++ ) iv2[idir] = *iptr++;
450  cptr = (char*)iptr;
451 
452  // perform shift
453  iv2 += shift;
454 
455  for( int c=a_comps.begin(); c<=a_comps.end(); c++ ) {
456  int iloc = getIndex(iv2, c);
457  m_dataPtr[iloc].linearIn( (void*)cptr );
458 // cptr += Tsz;
459  cptr += m_dataPtr[iloc].linearSize();
460  }
461  }
462 }
463 
464 template <class T> inline
465 int
466 IVSFAB<T>::getIndex(const IntVect& a_iv, const int& a_comp) const
467 {
468  CH_assert(isDefined());
469  CH_assert(m_ivs.contains(a_iv));
470  CH_assert((a_comp >= 0) && (a_comp < m_nComp));
471 
472  int ioffset = m_ivmap(a_iv, 0);
473  CH_assert(ioffset >= 0);
474  CH_assert(ioffset < m_nIvs);
475  //now add offset from componentnitude
476  ioffset += m_nIvs*a_comp;
477  return ioffset;
478 }
479 /********************/
480 template <class T> inline
481 void
482 IVSFAB<T>::clear()
483 {
484  m_nComp = 0;
485  m_nIvs = 0;
486  m_ivs.makeEmpty();
487  m_ivmap.clear();
488  if (m_dataPtr != NULL)
489  {
490  delete[] m_dataPtr;
491  m_dataPtr = NULL;
492  }
493  m_isDefined = false;
494 }
495 /*************************/
496 template <class T> inline
497 bool
498 IVSFAB<T>::isDefined() const
499 {
500  return (m_isDefined);
501 }
502 /*************************/
503 template <class T> inline
504 int
505 IVSFAB<T>::numIvs() const
506 {
507  return m_nIvs;
508 }
509 /*************************/
510 template <class T> inline
511 int
512 IVSFAB<T>::nComp() const
513 {
514  return m_nComp;
515 }
516 /*************************/
517 template <class T> inline
518 T&
519 IVSFAB<T>::operator() (const IntVect& a_ndin,
520  const int& a_comp)
521 {
522  CH_assert(isDefined());
523  CH_assert(a_comp >= 0);
524  CH_assert(a_comp < m_nComp);
525  int iloc = getIndex(a_ndin, a_comp);
526  return(m_dataPtr[iloc]);
527 }
528 /**************************/
529 template <class T> inline
530 const T&
531 IVSFAB<T>::operator() (const IntVect& a_ndin,
532  const int& a_comp) const
533 {
534  CH_assert(isDefined());
535  CH_assert(a_comp >= 0);
536  CH_assert(a_comp < m_nComp);
537  int iloc = getIndex(a_ndin, a_comp);
538  return(m_dataPtr[iloc]);
539 }
540 /******************/
541 template <class T> inline
542 const T*
543 IVSFAB<T>::dataPtr(const int& a_comp) const
544 {
545  CH_assert(isDefined());
546  CH_assert(a_comp >= 0);
547  CH_assert(a_comp < m_nComp);
548  return m_dataPtr + a_comp*m_nIvs;
549 }
550 /******************/
551 template <class T> inline
552 T*
553 IVSFAB<T>::dataPtr(const int& a_comp)
554 {
555  CH_assert(isDefined());
556  CH_assert(a_comp >= 0);
557  CH_assert(a_comp < m_nComp);
558  return m_dataPtr + a_comp*m_nIvs;
559 }
560 /******************/
561 template <class T> inline
562 const int*
563 IVSFAB<T>::loVect() const
564 {
565  return m_loVect.getVect();
566 }
567 /******************/
568 template <class T> inline
569 const int*
570 IVSFAB<T>::hiVect() const
571 {
572  return m_hiVect.getVect();
573 }
574 /******************/
575 /*
576 template <class T> inline
577 void
578 IVSFAB<T>::setDefaultValues()
579 {
580  m_isDefined = false;
581  m_dataPtr = NULL;
582  m_nIvs = 0;
583  m_nComp = 0;
584  m_loVect = IntVect::Unit;
585  m_hiVect = IntVect::Zero;
586 }
587 */
588 #include "NamespaceFooter.H"
589 #endif
pout
std::ostream & pout()
Use this in place of std::cout for program output.
IVSFAB::setVerbose
static void setVerbose(bool a_verbose)
Definition: IVSFABI.H:27
IntVectSet
An irregular domain on an integer lattice.
Definition: IntVectSet.H:44
CH_SPACEDIM
#define CH_SPACEDIM
Definition: SPACE.H:51
CH_assert
#define CH_assert(cond)
Definition: CHArray.H:37
IVSFAB::copy
void copy(const Box &a_fromBox, const Interval &a_destInterval, const Box &a_toBox, const IVSFAB< T > &a_src, const Interval &a_srcInterval)
Definition: IVSFABI.H:120
IVSFAB::define
void define(const IntVectSet &a_region, const int &a_nvarin)
Definition: IVSFABI.H:61
IVSFAB::m_ivs
IntVectSet m_ivs
Definition: IVSFAB.H:190
IVSIterator::ok
bool ok() const
returns true if this iterator is still in its IntVectSet
Definition: IntVectSet.H:711
IVSFAB::linearIn
void linearIn(void *buf, const Box &R, const Interval &comps)
Definition: IVSFABI.H:418
IVSFAB::nComp
int nComp() const
Definition: IVSFABI.H:512
IVSFAB
Definition: IVSFAB.H:31
IVSFAB::m_nComp
int m_nComp
Definition: IVSFAB.H:178
IVSFAB::dataPtr
const T * dataPtr(const int &a_comp) const
Definition: IVSFABI.H:543
IVSFAB::isDefined
bool isDefined() const
Definition: IVSFABI.H:498
IVSIterator::reset
void reset()
same as begin()
Definition: IntVectSet.H:722
IntVectSet::numPts
int numPts() const
Returns the number of IntVects in this IntVectSet.
Interval::size
int size() const
Definition: Interval.H:75
IVSFAB::setVal
void setVal(const T &value)
Definition: IVSFABI.H:111
SpaceDim
const int SpaceDim
Definition: SPACE.H:38
IntVectSet::minBox
const Box & minBox() const
Returns the minimum enclosing box of this IntVectSet.
IVSFAB::size
int size(const Box &R, const Interval &comps) const
Definition: IVSFABI.H:345
Interval
Structure for passing component ranges in code.
Definition: Interval.H:23
IntVect::Unit
static const IntVect Unit
Definition: IntVect.H:659
IVSFAB::~IVSFAB
~IVSFAB()
Definition: IVSFABI.H:46
Real
double Real
Definition: REAL.H:33
IntVectSet::isEmpty
bool isEmpty() const
Returns true if no IntVects are in this IntVectSet.
Box::bigEnd
const IntVect & bigEnd() const
Definition: Box.H:1768
IVSFAB::operator()
T & operator()(const IntVect &a_iv, const int &varlocin)
Definition: IVSFABI.H:519
Interval::begin
int begin() const
Definition: Interval.H:97
IntVect::Zero
static const IntVect Zero
Definition: IntVect.H:654
Box
A Rectangular Domain on an Integer Lattice.
Definition: Box.H:465
IVSFAB::linearOut
void linearOut(void *buf, const Box &R, const Interval &comps) const
Definition: IVSFABI.H:370
IVSFAB::loVect
const int * loVect() const
Definition: IVSFABI.H:563
IntVect
An integer Vector in SpaceDim-dimensional space.
Definition: CHArray.H:42
IVSFAB::numIvs
int numIvs() const
Definition: IVSFABI.H:505
Interval::end
int end() const
Definition: Interval.H:102
IVSIterator
Iterator for an IntVectSet.
Definition: IntVectSet.H:640
IVSFAB::clear
void clear()
Definition: IVSFABI.H:482
Box::smallEnd
const IntVect & smallEnd() const
{ Accessors}
Definition: Box.H:1754
IntVectSet::shift
void shift(const IntVect &iv)
Increment all the IntVects in this IntVectSet by the IntVect iv.
IVSFAB::getIndex
int getIndex(const IntVect &a_iv, const int &a_comp) const
Definition: IVSFABI.H:466
IVSFAB::hiVect
const int * hiVect() const
Definition: IVSFABI.H:570
IVSFAB::IVSFAB
IVSFAB()
Definition: IVSFABI.H:40
IVSFAB::getIVS
const IntVectSet & getIVS() const
Definition: IVSFABI.H:34
Generated on Sat Nov 28 2020 03:57:10 for Chombo + EB by   doxygen 1.8.6