Chombo + EB: BaseIFFABI.H Source File

00001 #ifdef CH_LANG_CC
00002 /*
00003 *      _______              __
00004 *     / ___/ /  ___  __ _  / /  ___
00005 *    / /__/ _ \/ _ \/  V \/ _ \/ _ \
00006 *    \___/_//_/\___/_/_/_/_.__/\___/
00007 *    Please refer to Copyright.txt, in Chombo's root directory.
00008 */
00009 #endif
00010 
00011 #ifndef _BASEIFFABI_H_
00012 #define _BASEIFFABI_H_
00013 #include "MayDay.H"
00014 #include "FaceIterator.H"
00015 #include <typeinfo>
00016 #include "NamespaceHeader.H"
00017 
00018 template <class T>
00019 bool BaseIFFAB<T>::s_verbose  = false;
00020 /******************/
00021 template <class T> inline   void
00022 BaseIFFAB<T>::setVerbose(bool a_verbose)
00023 {
00024   s_verbose = a_verbose;
00025 }
00026 /******************/
00027 template <class T> inline
00028 const EBGraph&
00029 BaseIFFAB<T>::getEBGraph() const
00030 {
00031   return m_ebgraph;
00032 }
00033 /******************/
00034 template <class T> inline
00035 BaseIFFAB<T>::BaseIFFAB()
00036 {
00037   setDefaultValues();
00038 }
00039 /******************/
00040 template <class T> inline
00041 BaseIFFAB<T>::~BaseIFFAB()
00042 {
00043   clear();
00044 }
00045 /******************/
00046 template <class T> inline
00047 BaseIFFAB<T>::BaseIFFAB(const IntVectSet& a_ivsin,
00048                         const EBGraph& a_ebgraph,
00049                         const int& a_direction,
00050                         const int& a_nvarin)
00051 {
00052   setDefaultValues();
00053   define(a_ivsin, a_ebgraph, a_direction, a_nvarin);
00054 }
00055 
00056 template <class T> Arena* BaseIFFAB<T>::s_Arena = NULL;
00057 
00058 /******************/
00059 template <class T> inline
00060 void
00061 BaseIFFAB<T>::define(const IntVectSet& a_ivsin,
00062                      const EBGraph& a_ebgraph,
00063                      const int& a_direction,
00064                      const int& a_nvarin)
00065 {
00066   clear();
00067   m_isDefined = true;
00068   CH_assert(a_nvarin > 0);
00069   CH_assert((a_direction >= 0) && (a_direction < SpaceDim));
00070   const ProblemDomain& domain = a_ebgraph.getDomain();
00071   m_direction = a_direction;
00072   m_ivs = a_ivsin;
00073   m_ebgraph = a_ebgraph;
00074   m_nComp = a_nvarin;
00075   Box minbox = a_ivsin.minBox();
00076   BaseFab<int> faceLocFab;
00077   BaseFab<bool> doneThat;
00078   if(!a_ivsin.isEmpty())
00079     {
00080       m_fab.resize(     surroundingNodes(minbox,a_direction), 1);
00081       faceLocFab.resize(surroundingNodes(minbox,a_direction), 1);
00082       doneThat.resize(  surroundingNodes(minbox,a_direction), 1);
00083       faceLocFab.setVal(-1);
00084       m_fab.setVal(NULL);
00085       doneThat.setVal(false);
00086     }
00087   //get the data size and save offsets
00088   //need to do this to avoid double-counting faces
00089   m_nFaces = 0;
00090   int nVoFs = 0;
00091   for(IVSIterator ivsit(m_ivs); ivsit.ok(); ++ivsit)
00092     {
00093       const IntVect& iv = ivsit();
00094       IntVect ivLo = iv-BASISV(m_direction);
00095       IntVect ivHi = iv+BASISV(m_direction);
00096       //check to see if on domain boundary
00097       int numVoFsHere = m_ebgraph.numVoFs(iv);
00098      
00099       nVoFs += numVoFsHere;
00100       //we hold the mapping in the high vof of the faces
00101       if(!doneThat(iv, 0))
00102         {
00103           doneThat(iv, 0) = true;
00104 
00105           bool onLoBoundary = false;
00106           if(!domain.isPeriodic(m_direction))
00107             {
00108               onLoBoundary = iv[m_direction] == domain.domainBox().smallEnd(m_direction);
00109             }
00110           int numFacesLoSide = numVoFsHere;
00111           if(!onLoBoundary)
00112             {
00113               int numVoFsLo = m_ebgraph.numVoFs(ivLo);
00114               numFacesLoSide = numVoFsHere*numVoFsLo;
00115             }
00116 
00117           faceLocFab(iv, 0) = m_nFaces;
00118           m_nFaces += numFacesLoSide;
00119         }
00120       if(!doneThat(ivHi, 0))
00121         {
00122           doneThat(ivHi, 0) = true;
00123 
00124           bool onHiBoundary = false;
00125           if(!domain.isPeriodic(m_direction))
00126             {
00127               onHiBoundary = iv[m_direction] == domain.domainBox().bigEnd(m_direction);
00128             }
00129           int numFacesHiSide = numVoFsHere;
00130           if(!onHiBoundary)
00131             {
00132               int numVoFsHi = m_ebgraph.numVoFs(ivHi);
00133               numFacesHiSide = numVoFsHere*numVoFsHi;
00134             }
00135 
00136           faceLocFab(ivHi, 0) = m_nFaces;
00137           m_nFaces += numFacesHiSide;
00138         }
00139     }
00140   // allocate the memory
00141   if (m_nFaces > 0)
00142     {
00143       // const IntVectSet& multi = m_ebgraph.getMultiCells( minbox);
00144       // pout()<< minbox << "nVoFs:"<<nVoFs<<" nFaces:"<<m_nFaces<<"\n";
00145 #ifdef CH_USE_MEMORY_TRACKING
00146       // if (s_Arena == NULL)
00147       //   {
00148       //     s_Arena = new BArena( (typeid(T)).name());
00149       //   }
00150 #else
00151       // if (s_Arena == NULL)
00152       //   {
00153       //     s_Arena = new BArena("");
00154       //   }
00155 #endif
00156       
00157       // if (s_Arena == NULL)
00158       //   {
00159       //     MayDay::Error("malloc in basefab failed");
00160       //   }
00161 
00162       m_truesize = m_nFaces*m_nComp;
00163 
00164       // Note: clear() was called above so this isn't a memory leak
00165       m_data = new T[m_truesize];
00166       // m_data = static_cast<T*>(s_Arena->alloc(m_truesize * sizeof(T)));
00167       
00168 #ifdef CH_USE_MEMORY_TRACKING
00169       // s_Arena->bytes += m_truesize * sizeof(T) + sizeof(BaseIFFAB<T>);
00170       // if (s_Arena->bytes > s_Arena->peak)
00171       //   {
00172       //     s_Arena->peak = s_Arena->bytes;
00173       //   }
00174 #endif
00175       doneThat.setVal(false);
00176     }
00177   //cache pointers to the first component of the first face into m_fab for fast indexing
00178   T* startLoc = &m_data[0];
00179   for(IVSIterator ivsit(m_ivs); ivsit.ok(); ++ivsit)
00180     {
00181       const IntVect& iv = ivsit();
00182       if(!doneThat(iv, 0))
00183         {
00184           doneThat(iv, 0) = true;
00185           int iface = faceLocFab(iv, 0);
00186           if(iface >= 0)
00187             {
00188               m_fab(iv, 0) = startLoc + iface;
00189             }
00190         }
00191       IntVect ivHi = iv+BASISV(m_direction);
00192       if(!doneThat(ivHi, 0))
00193         {
00194           doneThat(ivHi, 0) = true;
00195           int iface = faceLocFab(ivHi, 0);
00196           if(iface >= 0)
00197             {
00198               m_fab(ivHi, 0) = startLoc + iface;
00199             }
00200         }
00201     }
00202 }
00203 /******************/
00204 template <class T> inline
00205 const IntVectSet&
00206 BaseIFFAB<T>::getIVS() const
00207 {
00208   return m_ivs;
00209 }
00210 /******************/
00211 template <class T> inline
00212 int
00213 BaseIFFAB<T>::direction() const
00214 {
00215   return m_direction;
00216 }
00217 /******************/
00218 template <class T> inline
00219 void
00220 BaseIFFAB<T>::setVal(const T& a_value)
00221 {
00222   CH_assert(isDefined());
00223   for(int ivec = 0; ivec < m_nFaces*m_nComp; ivec++)
00224     {
00225       m_data[ivec] = a_value;
00226     }
00227 }
00228 /******************/
00229 template <class T> inline
00230 void
00231 BaseIFFAB<T>::setVal(int ivar, const T& a_value)
00232 {
00233   CH_assert(isDefined());
00234   CH_assert(ivar >= 0);
00235   CH_assert(ivar < m_nComp);
00236 
00237   for(int ivec = ivar; ivec < m_nFaces*m_nComp; ivec += m_nComp)
00238     {
00239       m_data[ivec] = a_value;
00240     }
00241 }
00242 /******************/
00243 template <class T> inline
00244 void
00245 BaseIFFAB<T>::copy(const Box& a_fromBox,
00246                    const Interval& a_dstInterval,
00247                    const Box& a_toBox,
00248                    const BaseIFFAB<T>& a_src,
00249                    const Interval& a_srcInterval)
00250 {
00251   CH_assert(isDefined());
00252   CH_assert(a_src.isDefined());
00253   CH_assert(a_srcInterval.size() == a_dstInterval.size());
00254   CH_assert(a_dstInterval.begin() >= 0);
00255   CH_assert(a_srcInterval.begin() >= 0);
00256   CH_assert(a_dstInterval.end()   < m_nComp);
00257   CH_assert(a_srcInterval.end()   < a_src.m_nComp);
00258 
00259   if((!m_ivs.isEmpty()) && (!a_src.m_ivs.isEmpty()))
00260     {
00261       CH_assert(a_fromBox == a_toBox);
00262       Box intBox = a_toBox;
00263       IntVectSet ivsIntersect = m_ivs;
00264       ivsIntersect &=  a_src.m_ivs;
00265       ivsIntersect &=  intBox;
00266       int compSize = a_srcInterval.size();
00267 
00268       FaceStop::WhichFaces stopCrit = FaceStop::SurroundingWithBoundary;
00269       FaceIterator faceit(ivsIntersect, m_ebgraph, m_direction, stopCrit);
00270       for(faceit.reset(); faceit.ok(); ++faceit)
00271         {
00272           const FaceIndex& face = faceit();
00273           for(int icomp = 0; icomp < compSize; icomp++)
00274             {
00275               int isrccomp = a_srcInterval.begin() + icomp;
00276               int idstcomp = a_dstInterval.begin() + icomp;
00277               (*this)(face, idstcomp) = a_src(face, isrccomp);
00278             }
00279         }
00280     }
00281 }
00282 /*********/
00283 template <class T> inline
00284 T*
00285 BaseIFFAB<T>::getIndex(const FaceIndex& a_face, const int& a_comp) const
00286 {
00287   //which face is this in the local vector (lexi on vof cell indices)
00288 
00289   T* dataPtr =   (T*)(m_fab(a_face.gridIndex(Side::Hi), 0));
00290   int faceLoc = getLocalVecIndex(a_face);
00291   dataPtr += faceLoc;
00292   dataPtr += m_nFaces*a_comp;
00293   return dataPtr;
00294 }
00295 /*********/
00296 template <class T> inline
00297 int
00298 BaseIFFAB<T>::getLocalVecIndex(const FaceIndex& a_face) const
00299 {
00300   int retval;
00301   if(!a_face.isBoundary())
00302     {
00303       //this checks that both vofs are defined
00304       CH_assert(a_face.cellIndex(Side::Lo) >= 0);
00305       CH_assert(a_face.cellIndex(Side::Hi) >= 0);
00306       int xlen = m_ebgraph.numVoFs(a_face.gridIndex(Side::Lo));
00307       int loCell = a_face.cellIndex(Side::Lo);
00308       int hiCell = a_face.cellIndex(Side::Hi);
00309       retval =  loCell + xlen*hiCell;
00310     }
00311   else
00312     {
00313       int loCell = a_face.cellIndex(Side::Lo);
00314       int hiCell = a_face.cellIndex(Side::Hi);
00315       //one should be -1, the other should be larger
00316       CH_assert(((loCell == -1)&&(hiCell > -1))||
00317                 ((hiCell == -1)&&(loCell > -1)));
00318       //return the one that is not -1
00319       retval = Max(loCell, hiCell);
00320     }
00321   return retval;
00322 }
00323 /********************/
00324 template <class T> inline
00325 void
00326 BaseIFFAB<T>::clear()
00327 {
00328   m_nComp = 0;
00329   m_nFaces = 0;
00330   m_direction = -1;
00331   m_ivs.makeEmpty();
00332   m_fab.clear();
00333   if (m_data != NULL)
00334     {
00335       delete[] m_data;
00336 #undef free
00337       // s_Arena->free(m_data);
00338 #ifdef CH_USE_MEMORY_TRACKING
00339       // s_Arena->bytes -= m_truesize * sizeof(T) + sizeof(BaseIFFAB<T>);
00340 #endif
00341       m_data = NULL;
00342     }
00343   m_isDefined = false;
00344 }
00345 /*************************/
00346 template <class T> inline
00347 bool
00348 BaseIFFAB<T>::isDefined() const
00349 {
00350   return (m_isDefined);
00351 }
00352 /*************************/
00353 template <class T> inline
00354 int
00355 BaseIFFAB<T>::numFaces() const
00356 {
00357   return m_nFaces;
00358 }
00359 /*************************/
00360 template <class T> inline
00361 int
00362 BaseIFFAB<T>::nComp() const
00363 {
00364   return m_nComp;
00365 }
00366 /*************************/
00367 template <class T> inline
00368 T&
00369 BaseIFFAB<T>::operator() (const FaceIndex& a_ndin,
00370                           const int& a_comp)
00371 {
00372   CH_assert(isDefined());
00373   CH_assert(a_comp >= 0);
00374   CH_assert(a_comp < m_nComp);
00375   CH_assert((m_ivs.contains(a_ndin.gridIndex(Side::Lo)) ||
00376              m_ivs.contains(a_ndin.gridIndex(Side::Hi))));
00377   CH_assert(a_ndin.direction() == m_direction);
00378 
00379   T* dataPtr =   getIndex(a_ndin, a_comp);
00380   return(*dataPtr);
00381 }
00382 /**************************/
00383 template <class T> inline
00384 const T&
00385 BaseIFFAB<T>::operator() (const FaceIndex& a_ndin,
00386                           const int& a_comp) const
00387 {
00388   CH_assert(isDefined());
00389   CH_assert(a_comp >= 0);
00390   CH_assert(a_comp < m_nComp);
00391   CH_assert((m_ivs.contains(a_ndin.gridIndex(Side::Lo)) ||
00392              m_ivs.contains(a_ndin.gridIndex(Side::Hi))));
00393   CH_assert(a_ndin.direction() == m_direction);
00394 
00395   T* dataPtr =   getIndex(a_ndin, a_comp);
00396   return(*dataPtr);
00397 }
00398 /******************/
00399 template <class T> inline
00400 T*
00401 BaseIFFAB<T>::dataPtr(const int& a_comp)
00402 {
00403   CH_assert(a_comp >= 0);
00404   CH_assert(a_comp <= m_nComp);
00405   return m_data + a_comp*m_nFaces;
00406 }
00407 /******************/
00408 template <class T> inline
00409 const T*
00410 BaseIFFAB<T>::dataPtr(const int& a_comp) const
00411 {
00412   CH_assert(a_comp >= 0);
00413   CH_assert(a_comp <= m_nComp);
00414   return m_data + a_comp*m_nFaces;
00415 }
00416 /******************/
00417 template <class T> inline
00418 void
00419 BaseIFFAB<T>::setDefaultValues()
00420 {
00421   m_isDefined = false;
00422   m_data = NULL;
00423   m_nFaces = 0;
00424   m_nComp = 0;
00425   m_direction = -1;
00426 }
00427 /******************/
00428 template <class T> inline
00429 int BaseIFFAB<T>::size(const Box& a_region,
00430                        const Interval& a_comps) const
00431 {
00432   CH_assert(isDefined());
00433   //create set of cells in fab that are also in the input region
00434   IntVectSet subIVS = m_ivs;
00435   subIVS &= a_region;
00436 
00437   FaceStop::WhichFaces stopCrit = FaceStop::SurroundingWithBoundary;
00438   FaceIterator faceit(subIVS, m_ebgraph, m_direction, stopCrit);
00439   Vector<FaceIndex> faces = faceit.getVector();
00440 
00441   //account for list of faces
00442   int facesize = linearListSize(faces);
00443 
00444   //add for each data pt
00445   int datasize = 0;
00446   for(int iface = 0; iface < faces.size(); iface++)
00447     {
00448       for(int icomp = a_comps.begin(); icomp <= a_comps.end(); icomp++)
00449         {
00450           const T& dataPt = (*this)(faces[iface], icomp);
00451           int pointsize = CH_XD::linearSize(dataPt);
00452           datasize += pointsize;
00453         }
00454     }
00455 
00456   int retval = facesize + datasize;
00457 
00458   if(s_verbose)
00459     {
00460       pout() << "###############" << std::endl;
00461       pout() << "BaseIFFAB size() for region " << m_ebgraph.getRegion() << std::endl;
00462       pout() << " a_comps  = ("   << a_comps.begin() << "," << a_comps.end()  << "),"
00463              << " a_region = "   << a_region << " subIVS = ";
00464       // dumpIVS(&subIVS);
00465       pout() << "m_ivs = ";
00466       // dumpIVS(&m_ivs);
00467       pout() << "size = " << retval << std::endl;
00468       pout() << "###############" << std::endl;
00469     }
00470 
00471   return retval;
00472  }
00473 /********************/
00474 template <class T> inline
00475 void BaseIFFAB<T>::linearOut(void* a_buf,
00476                              const Box& a_region,
00477                              const Interval& a_comps) const
00478 {
00479   CH_assert(isDefined());
00480   //create set of cells in fab that are also in the input region
00481   IntVectSet subIVS = m_ivs;
00482   subIVS &= a_region;
00483 
00484   FaceStop::WhichFaces stopCrit = FaceStop::SurroundingWithBoundary;
00485   FaceIterator faceit(subIVS, m_ebgraph, m_direction, stopCrit);
00486   Vector<FaceIndex> faces = faceit.getVector();
00487 
00488   //output the faces.
00489   unsigned char* charbuffer = (unsigned char *) a_buf;
00490   linearListOut(charbuffer, faces);
00491   charbuffer += linearListSize(faces);
00492 
00493   //output the data
00494   const BaseIFFAB<T>& thisFAB = *this;
00495   for(int iface = 0; iface < faces.size(); iface++)
00496     {
00497       const FaceIndex& face = faces[iface];
00498       for(int icomp = a_comps.begin(); icomp <= a_comps.end(); icomp++)
00499         {
00500           //output the data into the buffer
00501           const T& dataPt =  thisFAB(face, icomp);
00502           CH_XD::linearOut(charbuffer, dataPt);
00503           //increment the buffer offset
00504           charbuffer += CH_XD::linearSize(dataPt);
00505         }
00506     }
00507 }
00508 /********************/
00509 template <class T> inline
00510 void BaseIFFAB<T>::linearIn(void* a_buf, const Box& a_region, const Interval& a_comps)
00511 {
00512   CH_assert(isDefined());
00513   //input the faces
00514   //input the vofs
00515   IntVectSet subIVS = m_ivs;
00516   subIVS &= a_region;
00517 
00518   FaceStop::WhichFaces stopCrit = FaceStop::SurroundingWithBoundary;
00519   FaceIterator faceit(subIVS, m_ebgraph, m_direction, stopCrit);
00520   Vector<FaceIndex> faces = faceit.getVector();
00521 
00522   unsigned char* charbuffer = (unsigned char *) a_buf;
00523   linearListIn(faces, charbuffer);
00524   charbuffer += linearListSize(faces);
00525 
00526   for(int iface = 0; iface < faces.size(); iface++)
00527     {
00528       const FaceIndex& face = faces[iface];
00529       for(int icomp = a_comps.begin(); icomp <= a_comps.end(); icomp++)
00530         {
00531           //input the data
00532           T& dataPt = (*this)(face, icomp);
00533           CH_XD::linearIn(dataPt, charbuffer) ;
00534           //increment the buffer offset
00535           charbuffer += CH_XD::linearSize(dataPt);
00536         }
00537     }
00538 }
00539 /********************/
00540 /******************/
00541 
00542 #include "NamespaceFooter.H"
00543 #endif