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