Chombo + EB + MF  3.2
BinFabImplem.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 _BINFABIMPLEM_H_
12 #define _BINFABIMPLEM_H_
13 
14 #include <cmath>
15 
16 #include "BoxIterator.H"
17 #include "parstream.H"
18 #include "NamespaceHeader.H"
19 
20 //
21 // Implementation.=====================================
22 //
23 
24 template <class T>
26  : BaseFab<List<T> >(),
27  //[NOTE: dont use RealVect::Zero because of unreliable order of static initializations. <dbs>]
28  m_origin(D_DECL6(0,0,0,0,0,0)),
29  m_mesh_spacing(D_DECL6(0,0,0,0,0,0))
30 {
31 }
32 
33 template <class T>
34 BinFab<T>::BinFab(const Box& a_domain,
35  const RealVect& a_meshSpacing,
36  const RealVect& a_origin)
37 {
38  define(a_domain, a_meshSpacing, a_origin);
39 }
40 
41 #ifdef __IBMCPP__
42 //[NOTE: this is here only because the IBM xlC compiler tries to compile
43 // DefaultDataFactory<BinFab<BinItem>> in particleTest.cpp, even
44 // though it isn't used anywhere. Stupid compiler. <dbs> July2005]
45 template <class T>
46 BinFab<T>::BinFab(const Box& a_domain, int a_ncomps)
47 {
48  MayDay::Abort("BinFab(Box,int): Dont use this constructor");
49 }
50 #endif
51 
52 template <class T>
54 {
55  this->undefine();
56 }
57 
58 template <class T>
59 void BinFab<T>::define(const Box& a_domain,
60  const RealVect& a_meshSpacing,
61  const RealVect& a_origin)
62 {
63  // BaseFab vars
64  this->m_domain = a_domain;
65  this->m_nvar = 1;
66  this->m_numpts = this->m_domain.numPts();
67  BaseFab<List<T> >::define();
68  // BinFab vars
69  m_mesh_spacing = a_meshSpacing;
70  m_origin = a_origin;
71 }
72 
73 // set the mesh spacing
74 //[NOTE: this is a vector even though Chombo doesn't support isotropic grids yet.]
75 template <class T>
76 void BinFab<T>::meshSpacing(const RealVect & a_meshSpacing)
77 {
78  m_mesh_spacing = a_meshSpacing;
79 }
80 
81 // retrieve a copy of the mesh spacing
82 template <class T>
84 {
85  return m_mesh_spacing;
86 }
87 
88 // set the coordinate origin;
89 // i.e. the physical location of the lower corner of cell IntVect::Zero
90 template <class T>
91 void BinFab<T>::origin(const RealVect& a_origin)
92 {
93  m_origin = a_origin;
94 }
95 
96 // retrieve a copy of the coordinate origin
97 template <class T>
99 {
100  return m_origin;
101 }
102 
103 ////////////////////////////////////////////////////////////////
104 
105 /// Insert a single item into this BinFab in the appropriate cell.
106 template <class T>
107 void BinFab<T>::addItem(const T& a_item)
108 {
109  int comp = 0;
110  IntVect binLoc = locateBin(a_item);
111  if (this->m_domain.contains(binLoc))
112  {
113  this->operator()(binLoc,comp).append(a_item);
114  }
115 }
116 
117 /// Insert a single item into this BinFab in the specified cell.
118 template <class T>
119 void BinFab<T>::addItem(const T& a_item, const IntVect& a_binLoc)
120 {
121  int comp = 0;
122  if (this->m_domain.contains(a_binLoc))
123  {
124  this->operator()(a_binLoc,comp).append(a_item);
125  }
126 }
127 
128 /// Insert multiple items into this BinFab.
129 template <class T>
130 void BinFab<T>::addItems(const List<T>& a_List)
131 {
132  if (a_List.isNotEmpty())
133  {
134  // loop through items in List, and add to appropriate locations
135  ListIterator<T> lit(a_List);
136  int comp = 0;
137  IntVect binLoc;
138 
139  for (lit.rewind(); lit; ++lit)
140  {
141  const T& thisItem = a_List[lit];
142  binLoc = locateBin(thisItem);
143 
144  if (this->m_domain.contains(binLoc))
145  {
146  this->operator()(binLoc,comp).append(thisItem);
147  }
148  }
149  }
150 }
151 
152 ////////////////////////////////////////////////////////////////
153 
154 /// Insert multiple items into BinFab and remove from source list.
155 // [NOTE: this is mostly the same as the previous version, except for the inner loop]
156 template <class T>
158 {
159  if (a_List.isNotEmpty())
160  {
161  // loop through items in List, add to appropriate bins,
162  // and remove from List
163  ListIterator<T> lit(a_List);
164  int comp = 0;
165  IntVect binLoc;
166 
167  for (lit.rewind(); lit; )
168  {
169  const T& thisItem = a_List[lit];
170  binLoc = locateBin(thisItem);
171  if (this->m_domain.contains(binLoc)) //[NOTE: BaseFab domain, not problem domain]
172  {
173  //this->operator()(binLoc,comp).append(thisItem);
174  //a_List.remove(lit);
175  this->operator()(binLoc,comp).transfer(lit);
176  }
177  else
178  {
179  // removing a list element has the side-effect of moving the
180  // iterator forward, so only increment explicitly if not removed
181  ++lit ;
182  }
183  }
184  }
185 }
186 
187 /// Move or copy items from arg list to *this, depending on in or not-in the given Box.
188 // [NOTE: this is mostly the same as the previous version, except for the inner loop]
189 template <class T>
190 void BinFab<T>::addItemsDestructive(List<T>& a_List, const Box& a_valid ,bool a_in)
191 {
192  CH_assert( a_valid.cellCentered() );
193  if (a_List.isNotEmpty())
194  {
195  // loop through items in List, add to appropriate bins,
196  // and remove from List
197  ListIterator<T> lit(a_List);
198  int comp = 0;
199  IntVect binLoc;
200 
201  for (lit.rewind(); lit; )
202  {
203  const T& thisItem = a_List[lit];
204  binLoc = locateBin(thisItem);
205  if (this->m_domain.contains(binLoc)) //[NOTE: BaseFab domain, not problem domain]
206  {
207 
208  // If the a_in flag is true, then delete the particle from the
209  // List if it is in the valid Box; if the a_in flag is false,
210  // delete the particle if it is NOT in the Box.
211  // Transfering a list element has the side-effect of moving the
212  // iterator forward, so explicitly increment the iterator only
213  // if the particle was not transfered. Wish C++ have a nice XOR function.
214  if ( a_valid.contains(binLoc) )
215  {
216  if ( a_in )
217  {
218  //this->operator()(binLoc,comp).append(thisItem);
219  //a_List.remove(lit);
220  this->operator()(binLoc,comp).transfer(lit);
221  }
222  else
223  {
224  this->operator()(binLoc,comp).append(thisItem);
225  ++lit ;
226  }
227  }
228  else
229  {
230  if ( a_in )
231  {
232  this->operator()(binLoc,comp).append(thisItem);
233  ++lit;
234  }
235  else
236  {
237  //this->operator()(binLoc,comp).append(thisItem);
238  //a_List.remove(lit);
239  this->operator()(binLoc,comp).transfer(lit);
240  }
241  }
242 
243  }
244  else
245  {
246  // not added, so just increment the iterator
247  ++lit ;
248  }
249  }
250  }
251 }
252 
253 ////////////////////////////////////////////////////////////////
254 
255 // Look at every item in the BinFab and move to a different cell
256 // if it is in the wrong one. Useful after computing new physical
257 // locations of particles, e.g. moving in time.
258 template <class T>
260 {
261  // loop through each bin and determine if each binItem is
262  // in the correct bin.
263  BoxIterator bit(this->m_domain);
264  int comp = 0;
265  IntVect binLoc;
266 // int lost = 0;
267 
268  for (BoxIterator bit(this->m_domain); bit.ok(); ++bit)
269  {
270  const IntVect thisIV = bit();
271  List<T>& thisList = this->operator()(thisIV,comp);
272 
273  if (thisList.isNotEmpty())
274  {
275  // now index through the List and figure out which bin each
276  // particle should be in
277  for (ListIterator<T> thisLit(thisList); thisLit;)
278  {
279  T& thisItem = thisList[thisLit];
280  binLoc = locateBin(thisItem);
281 
282  if (binLoc != thisIV)
283  {
284  // binItem needs to be moved.
285  // note that if binItem moves outside of domain, it is lost
286  if (this->m_domain.contains(binLoc))
287  {
288  // this->operator()(binLoc,comp).append(thisItem);
289  // thisList.remove(thisLit);
290  this->operator()(binLoc,comp).transfer(thisLit);
291  }
292  else
293  {
294 // ++lost ;
295  thisList.remove(thisLit);
296  }
297  }
298  else
299  {
300  // if this list element is not moved, increment the iterator
301  ++thisLit;
302  }
303  } // end loop over items in this List of binItems
304  } // end if this list has items
305  } // end loop over bins
306 // return lost;
307 }
308 
309 // Sort all the items into the appropriate bins and return items that are in (or not in) the \p Box.
310 // The returned items are appended to the list, so the caller should clear() first if apprpriate.
311 // Same as reBin(void) except it keeps track of the lost items.
312 //[NOTE: items that are moved forward in the box wil be looked at twice. <dbs>]
313 template <class T>
314 void BinFab<T>::reBin(List<T>& a_lost, const Box & a_valid, bool a_in)
315 {
316 //XXX What to do if a_valid is not contained in m_domain?? <dbs>
317  int comp = 0; //assume only one component
318  IntVect binLoc;
319  Box valid = a_valid.isEmpty() ? this->m_domain : a_valid ;
320 
321  // Loop through _all_ the bins in this BinFab, regardless of the valid Box
322  for (BoxIterator bit(this->m_domain); bit.ok(); ++bit)
323  {
324  const IntVect thisIV = bit();
325  List<T>& thisList = this->operator()(thisIV,comp);
326  if (thisList.isNotEmpty())
327  {
328  // Compute the proper cell for each element in this list
329  for (ListIterator<T> thisLit(thisList); thisLit;) //[NOTE: no auto-increment operation]
330  {
331  T& thisItem = thisList[thisLit];
332  binLoc = locateBin(thisItem);
333  if (binLoc != thisIV)
334  {
335  // move it or return it, as appropriate
336  if (valid.contains(binLoc))
337  { // new location is in the arg box
338  if (a_in) a_lost.transfer(thisLit) ;//return it
339  else this->operator()(binLoc,comp).transfer(thisLit); //move it
340  }
341  else if (this->m_domain.contains(binLoc))
342  { // new location not in the arg box but in the domain
343  this->operator()(binLoc,comp).transfer(thisLit); //move it
344  }
345  else
346  { // new location is not in the domain
347  if (!a_in) a_lost.transfer(thisLit); //return it
348  else thisList.remove(thisLit); // and not in the box, so throw it away
349  }
350  //
351  }
352  else
353  {
354  // increment the iterator only if the item stays in place
355  ++thisLit;
356  }
357  } // end loop over items in this List of binItems
358  } // end if this list has items
359  } // end loop over bins
360 }
361 
362 ////////////////////////////////////////////////////////////////
363 
364 // Move items from a_srcBox in the source BinFab to a_destBox in *this.
365 // This is what ::copy() would do if it allowed a non-const source.
366 //[NOTE: if the boxes are different, the items will end up in cells with
367 // different IntVects from the ones they started in. <dbs>]
368 template <class T>
370  const Box& a_srcBox,
371  const Box& a_destBox)
372 {
373  CH_assert(a_src.box().contains(a_srcBox));
374  CH_assert(a_destBox.sameSize(a_srcBox));
375  // Ignore components; BinFab assumes only one.
376  CH_assert( a_src.nComp() == 1 );
377  const int comp = 0, ncomps = 1;
378  //XXX stupid lack of operator==
379  //XXX CH_assert( a_comps == Interval(0,0) );
380  // loop through the cells in the intersection of the two boxes
381  ForAllThisBNNXCBN(List<T>,a_destBox,comp,ncomps,a_src,a_srcBox,comp)
382  {
383  if (a_srcR.isNotEmpty())
384  {
385  // cast away const-ness that is assumed by the ForAll...
386  List<T>& src = (List<T>&)a_srcR ;
387 //XXX debug
388 //XXX pout() << "dbg: transfer: i=" << iR << " has " << a_srcR.length() << endl ;
389  thisR.catenate(src);
390  }
391  } EndForTX
392 }
393 
394 ////////////////////////////////////////////////////////////////
395 
396 // Deletes all items from this BinFab and resets it to the undefined state.
397 template <class T>
399 {
400  this->undefine();
401  this->m_domain = Box();
402  this->m_nvar = 0;
403  this->m_numpts = 0;
404  m_origin = RealVect::Zero;
405  m_mesh_spacing = RealVect::Zero;
406 }
407 
408 ////////////////////////////////////////////////////////////////
409 
410 // Counts the number of items in the cells of a_box.
411 template <class T>
412 int BinFab<T>::numItems(const Box& a_box) const
413 {
414  // Ignore components; BinFab assumes only one.
415  const int comp = 0 ;
416  // default is an empty box, so use this BinFab's box instead
417  // [NOTE: this leaves the caller no way to ask for an empty box;
418  // it must avoid calling this routine in that case. Sorry.]
419  Box ibox = a_box;
420  if (ibox.isEmpty())
421  {
422  ibox = this->box();
423  }
424  // Loop over all bins in the BinFab
425  int totalSize = 0;
426  for (BoxIterator bi(ibox); bi.ok(); ++bi)
427  {
428  totalSize += this->operator()(bi(),comp).length();
429  }
430 //XXX -- debug
431 //XXX pout() << "dbg: BinFab::numItems is " << totalSize << " for box " << ibox << endl ;
432  return totalSize;
433 }
434 
435 ////////////////////////////////////////////////////////////////
436 
437 // Computes the linear size (in bytes) of the items in the cells of a_box.
438 template <class T>
439 size_t BinFab<T>::size(const Box& a_box, const Interval& a_comps) const
440 {
441  // Ignore components; BinFab assumes only one.
442  const int ncomps = 1 ;
443  //XXX stupid lack of operator==
444  //XXX CH_assert( a_comps == Interval(0,0) );
445  // Count the number of items in the specified Box.
446  size_t totalSize = numItems(a_box);
447  // Compute size as the number of items times the size of the items
448  {
449  // Compute the size of an item, assuming all are the same size.
450  int sizeOfT;
451  T tmp;
452  sizeOfT = tmp.size();
453  // totalSize has total number of items, so now compute total size
454  totalSize *= sizeOfT;
455  }
456  // Will also need to include an integer for each list to specify how
457  // many items are in each list.
458  //[NOTE: Is this the right thing to do? An alternate approach would be to
459  // simply unpack one List<T> for all the items and then call addItems()
460  // with that list. But that probably will be slower. <dfmartin>]
461  size_t numBins = a_box.numPts() * ncomps ;
462  totalSize += numBins*sizeof(int);
463  return totalSize;
464 }
465 
466 ////////////////////////////////////////////////////////////////
467 
468 // Linearizes the items in the cells of a_box into a_buf.
469 template <class T>
470 void BinFab<T>::linearOut(void* a_buf, const Box& a_box, const Interval& a_comps) const
471 {
472  // Ignore components; BinFab assumes only one.
473  const int comp = 0 ;
474  //XXX stupid lack of operator==
475  //XXX CH_assert( a_comps == Interval(0,0) );
476  // make a usable pointer to the buffer
477  char* buf_ch = (char*)a_buf;
478  // All we do in this function is loop over the Box
479  // and call linearOut() on the List elements.
480  for (BoxIterator bit(a_box); bit.ok(); ++bit)
481  {
482  const List<T>& thisList = this->operator()(bit(), comp);
483  // write the number of list items first
484  int *intBuffer = (int*)buf_ch;
485  *intBuffer = thisList.length();
486  buf_ch += sizeof(int);
487  // now loop over the items in the list and write each one
488  for (ListIterator<T> lit(thisList); lit.ok(); ++lit)
489  {
490  thisList[lit].linearOut(buf_ch);
491  buf_ch += thisList[lit].size();
492  }
493  }
494 }
495 
496 ////////////////////////////////////////////////////////////////
497 
498 // Same as linearOut() except it clear()s the List after writing it.
499 template <class T>
500 void
501 BinFab<T>::linearOutDestructive(void* a_buf, const Box& a_box, const Interval& a_comps)
502 {
503  const int comp = 0 ;
504  //XXX stupid lack of operator==
505  //XXX CH_assert( a_comps == Interval(0,0) );
506  char* buf_ch = (char*)a_buf;
507  for (BoxIterator bit(a_box); bit.ok(); ++bit)
508  {
509  List<T>& thisList = this->operator()(bit(), comp);
510  int *intBuffer = (int*)buf_ch;
511  *intBuffer = thisList.length();
512  buf_ch += sizeof(int);
513  for (ListIterator<T> lit(thisList); lit.ok(); ++lit)
514  {
515  thisList[lit].linearOut(buf_ch);
516  buf_ch += thisList[lit].size();
517  }
518  thisList.clear();
519  }
520 }
521 
522 ////////////////////////////////////////////////////////////////
523 
524 // Fills this BinFab from the linearization in a_buf of the cells in a_box.
525 template <class T>
526 void
527 BinFab<T>::linearIn(void* a_buf, const Box& a_box, const Interval& a_comps)
528 {
529  // Should be just the inverse of linearOut().
530  // Ignore components; BinFab assumes only one.
531  const int comp = 0 ;
532  //XXX stupid lack of operator==
533  //XXX CH_assert( a_comps == Interval(0,0) );
534  // make a usable pointer to the buffer
535  char *buf_ch = (char*)a_buf;
536  // All we do in this function is loop over the Box
537  // and call linearOut() on the List elements.
538  for (BoxIterator bit(a_box); bit.ok(); ++bit)
539  {
540  List<T>& thisList = this->operator()(bit(), comp);
541  // clear the list to prevent accumulation of binItems
542  thisList.clear();
543 
544  // first extract the length of the list
545  int numItems = *((int*)buf_ch);
546  buf_ch += sizeof(int);
547  // loop over the list items and make the list
548  for (int n=0; n<numItems; ++n)
549  {
550  T thisBinItem;
551  thisBinItem.linearIn( buf_ch );
552  thisList.append(thisBinItem);
553  buf_ch += thisBinItem.size();
554  }
555  }
556 }
557 
558 //X /*******************************************************************************
559 //X This is an example that I want to keep here to remind that the loops over boxes
560 //X should be unrolled through the BaseFabMacro. <fm>
561 //X *******************************************************************************
562 //X template <class T>
563 //X inline
564 //X void
565 //X BaseFab<T>::linearIn(void* buf, const Box& R, const Interval& comps)
566 //X {
567 //X T* buffer = (T*)buf;
568 //X ForAllThisBNN(T,R,comps.begin(), comps.size())
569 //X {
570 //X thisR = *buffer;
571 //X ++buffer;
572 //X } EndFor;
573 //X }
574 //X */
575 
576 ////////////////////////////////////////////////////////////////
577 
578 // Compute the IntVect of the cell containing a_binItem.
579 template <class T>
580 inline
581 IntVect
582 BinFab<T>::locateBin(const T& a_binItem) const
583 {
584  IntVect binLoc;
585  RealVect thisPos = a_binItem.position();
586  thisPos -= m_origin;
587  thisPos /= m_mesh_spacing;
588 
589  for ( int d=0 ; d<SpaceDim ; ++d )
590  {
591  binLoc[d] = (int)floor(thisPos[d]);
592  }
593  return binLoc;
594 }
595 
596 #include "NamespaceFooter.H"
597 #endif
#define D_DECL6(a, b, c, d, e, f)
Definition: CHArray.H:39
void catenate(List< T > &src)
Appends a copy of all items in List<T> src to this List<T>.
Definition: ListImplem.H:211
bool ok()
Definition: BoxIterator.H:281
IntVect locateBin(const RealVect a_x, const RealVect a_dx, const RealVect a_origin)
compute the cell index containing the physical position of the item
Definition: ParticleData.H:47
virtual void addItem(const T &a_item)
Copy one item into the appropriate bin.
Definition: BinFabImplem.H:107
virtual int numItems(const Box &a_box=Box()) const
Count the number of items in this BinFab in the specified Box.
Definition: BinFabImplem.H:412
virtual RealVect meshSpacing() const
Retrieve the mesh cell size, in physical coordinates.
Definition: BinFabImplem.H:83
#define CH_assert(cond)
Definition: CHArray.H:37
bool isNotEmpty() const
Returns true if the List<T> is not empty.
Definition: List.H:627
void transfer(ListIterator< T > &lit)
Transfer the object pointed to by lit from the List<T> lit is associated with to this one...
Definition: ListImplem.H:317
void remove(const T &value)
Removes all objects in the List<T> equal to value.
Definition: ListImplem.H:356
void append(const T &value)
Adds a copy of the value to the end of the List<T>.
Definition: List.H:582
void clear()
Removes all objects from the List<T>.
Definition: ListImplem.H:263
IntVect size() const
Should be inherited from BaseFab but it isn't, sigh; code is the same.
Definition: BinFab.H:160
void rewind()
Reset this ListIterator<T> to point to the first element in the List<T>.
Definition: List.H:422
iterates through the IntVects of a Box
Definition: BoxIterator.H:37
virtual void reBin()
Sort all the items into the appropriate bins, losing items not in the FAB's domain.
Definition: BinFabImplem.H:259
const int SpaceDim
Definition: SPACE.H:38
virtual void linearIn(void *a_buf, const Box &a_box, const Interval &a_comps=Interval(0, 0))
Extract a serialized (packed) representation from a_buf.
Definition: BinFabImplem.H:527
int nComp() const
{ accessors}
Definition: EBInterface.H:45
virtual void linearOut(void *a_buf, const Box &a_box, const Interval &a_comps) const
Write a serialized (packed) representation into a_buf.
Definition: BinFabImplem.H:470
static const RealVect Zero
Definition: RealVect.H:421
Structure for passing component ranges in code.
Definition: Interval.H:23
A Doubly-Linked List Class.
Definition: List.H:21
virtual IntVect locateBin(const T &a_item) const
compute the cell index containing the physical position of the item
Definition: BinFabImplem.H:582
Iterator over a List.
Definition: List.H:20
bool cellCentered() const
Definition: Box.H:1939
virtual ~BinFab()
Destructor.
Definition: BinFabImplem.H:53
int length() const
Returns the number of objects in the List<T>.
Definition: ListImplem.H:56
bool isEmpty() const
{ Comparison Functions}
Definition: Box.H:1862
Definition: BaseFab.H:81
A Rectangular Domain on an Integer Lattice.
Definition: Box.H:469
A Real vector in SpaceDim-dimensional space.
Definition: RealVect.H:41
virtual void linearOutDestructive(void *buf, const Box &a_box, const Interval &a_comps=Interval(0, 0))
Write a serialized (packed) representation into a_buf and remove from *this.
Definition: BinFabImplem.H:501
virtual void clear()
Delete all the items in this BinFab and reset it to the null Box.
Definition: BinFabImplem.H:398
BinFab()
Null constructor, copy constructor and operator= can be compiler defined.
Definition: BinFabImplem.H:25
size_t numPts() const
bool sameSize(const Box &b) const
Definition: Box.H:1912
bool contains(const IntVect &p) const
Definition: Box.H:1887
An integer Vector in SpaceDim-dimensional space.
Definition: CHArray.H:42
bool ok() const
Return true if the iterator is not past the end of the list.
Definition: List.H:465
Base class for particle data on a Box.
Definition: BinFab.H:30
RealVect origin() const
Retrieve the physical coordinates of the lower corner of cell zero.
Definition: BinFabImplem.H:98
virtual void addItemsDestructive(List< T > &a_list)
Copy items into the appropriate bins and remove them from the List.
Definition: BinFabImplem.H:157
virtual void transfer(BinFab< T > &a_src, const Box &a_srcBox, const Box &a_destBox)
Move particles from a BinFab to the overlapping cells in this one.
Definition: BinFabImplem.H:369
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).
virtual void addItems(const List< T > &a_list)
Copy items into the appropriate bins.
Definition: BinFabImplem.H:130
const Box & box() const