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