AMRParticlePaths: BaseFab.H Source File

00001 /*
00002 ** (c) 1996-2000 The Regents of the University of California (through
00003 ** E.O. Lawrence Berkeley National Laboratory), subject to approval by
00004 ** the U.S. Department of Energy.  Your use of this software is under
00005 ** license -- the license agreement is attached and included in the
00006 ** directory as license.txt or you may contact Berkeley Lab's Technology
00007 ** Transfer Department at TTD@lbl.gov.  NOTICE OF U.S. GOVERNMENT RIGHTS.
00008 ** The Software was developed under funding from the U.S. Government
00009 ** which consequently retains certain rights as follows: the
00010 ** U.S. Government has been granted for itself and others acting on its
00011 ** behalf a paid-up, nonexclusive, irrevocable, worldwide license in the
00012 ** Software to reproduce, prepare derivative works, and perform publicly
00013 ** and display publicly.  Beginning five (5) years after the date
00014 ** permission to assert copyright is obtained from the U.S. Department of
00015 ** Energy, and subject to any subsequent five (5) year renewals, the
00016 ** U.S. Government is granted for itself and others acting on its behalf
00017 ** a paid-up, nonexclusive, irrevocable, worldwide license in the
00018 ** Software to reproduce, prepare derivative works, distribute copies to
00019 ** the public, perform publicly and display publicly, and to permit
00020 ** others to do so.
00021 */
00022 
00023 #ifndef BL_BASEFAB_H
00024 #define BL_BASEFAB_H
00025 //
00026 // $Id: BaseFab.H,v 1.48 2002/03/26 20:04:06 lijewski Exp $
00027 //
00028 #include <winstd.H>
00029 
00030 #if defined(BL_OLD_STL)
00031 #include <math.h>
00032 #include <stdlib.h>
00033 #else
00034 #include <cmath>
00035 #include <cstdlib>
00036 #endif
00037 
00038 #include <algorithm>
00039 
00040 #include <BLassert.H>
00041 #include <Box.H>
00042 #include <BoxList.H>
00043 #include <CArena.H>
00044 #include <Looping.H>
00045 #include <REAL.H>
00046 //
00047 // :
00049 
00095 template <class T>
00096 class BaseFab
00097 {
00098 public:
00099     typedef T value_type;
00100 
00106     BaseFab ();
00107 
00108     BaseFab (const BaseFab&);
00109     BaseFab& operator= (const BaseFab&);
00110     BaseFab& operator= (const T&);
00111     //
00113     //
00114     explicit BaseFab (const Box& bx,
00115                       int        n = 1);
00116     //
00118     //
00119     ~BaseFab ();
00120 
00130     void resize (const Box& b,
00131                  int        N = 1);
00132 
00136     void clear ();
00137     //
00139     //
00140     int nComp () const;
00141     //
00143     //
00144     const Box& box () const;
00145 
00149     const int* length () const;
00150 
00154     const IntVect& smallEnd () const;
00155  
00159     const IntVect& bigEnd () const;
00160 
00167     const int* loVect () const;
00168  
00175     const int* hiVect () const;
00176  
00180     bool contains (const BaseFab<T>& fab) const;
00181  
00185     bool contains (const Box& bx) const;
00186  
00195     T* dataPtr (int N = 0);
00196     //
00198     //
00199     const T* dataPtr (int N = 0) const;
00200     //
00202     //
00203     bool isAllocated () const;
00204 
00210     T& operator() (const IntVect& p,
00211                    int            N);
00212     //
00214     //
00215     T& operator() (const IntVect& p);
00216     //
00218     //
00219      const T& operator() (const IntVect& p,
00220                           int            N) const;
00221     //
00223     //
00224     const T& operator() (const IntVect& p) const;
00225 
00226 
00231     void getVal (T*             data,
00232                  const IntVect& pos,
00233                  int            N,
00234                  int            numcomp) const;
00235 
00239     void getVal (T*             data,
00240                  const IntVect& pos) const;
00241 
00247     void setVal (T          x,
00248                  const Box& bx,
00249                  int        nstart,
00250                  int        ncomp);
00251 
00255     void setVal (T          x,
00256                  const Box& bx,
00257                  int        N);
00258     //
00260     //
00261     void setVal (T   x,
00262                  int N);
00263     //
00265     //
00266     void setVal (T x);
00267 
00272     void setComplement (T          x,
00273                         const Box& b,
00274                         int        ns,
00275                         int        num);
00276 
00294     BaseFab<T>& copy (const BaseFab<T>& src,
00295                       const Box&        srcbox,
00296                       int               srccomp,
00297                       const Box&        destbox,
00298                       int               destcomp,
00299                       int               numcomp);
00300 
00307     BaseFab<T>& copy (const BaseFab<T>& src,
00308                       int               srccomp,
00309                       int               destcomp,
00310                       int               numcomp = 1);
00311 
00318     BaseFab<T>& copy (const BaseFab<T>& src,
00319                       const Box&        destbox);
00320 
00326     BaseFab<T>& copy (const BaseFab<T>& src);
00327 
00332     BaseFab<T>& shift (const IntVect& v);
00333 
00338     BaseFab<T>& shift (int idir,
00339                        int n_cell);
00340 
00345     BaseFab<T>& shiftHalf (int dir,
00346                            int num_halfs);
00347 
00352     BaseFab<T>& shiftHalf (const IntVect& num_halfs);
00353 
00360     Real norm (int p,
00361                int scomp = 0,
00362                int ncomp = 1) const;
00363     //
00365     //
00366     Real norm (const Box& subbox,
00367                int        p,
00368                int        scomp = 0,
00369                int        ncomp = 1) const;
00370     //
00372     //
00373     void abs ();
00374     //
00376     //
00377     void abs (int comp,
00378               int numcomp=1);
00379     //
00381     //
00382     void abs (const Box& subbox,
00383               int        comp = 0,
00384               int        numcomp=1);
00385     //
00387     //
00388     T min (int comp = 0) const;
00389     //
00391     //
00392     T min (const Box& subbox,
00393            int        comp = 0) const;
00394     //
00396     //
00397     T max (int comp = 0) const;
00398     //
00400     //
00401     T max (const Box& subbox,
00402            int        comp = 0) const;
00403     //
00405     //
00406     IntVect minIndex (int comp = 0) const;
00407 
00411     IntVect minIndex (const Box& subbox,
00412                       int        comp = 0) const;
00413     //
00415     //
00416     IntVect maxIndex (int comp = 0) const;
00417  
00421     IntVect maxIndex (const Box& subbox,
00422                       int        comp = 0) const;
00423 
00429     int maskLT (BaseFab<int>& mask,
00430                 T             val,
00431                 int           comp = 0) const;
00432  
00436     int maskLE (BaseFab<int>& mask,
00437                 T             val,
00438                 int           comp = 0) const;
00439     //
00441     //
00442     int maskEQ (BaseFab<int>& mask,
00443                 T             val,
00444                 int           comp = 0) const;
00445     //
00447     //
00448     int maskGT (BaseFab<int>& mask,
00449                 T             val,
00450                 int           comp = 0) const;
00451  
00455     int maskGE (BaseFab<int>& mask,
00456                 T             val,
00457                 int           comp = 0) const;
00458 
00459     //
00461     //
00462     void patternFill (int mark = 0);
00463 
00468     void copyRev (const Box&        destbox,
00469                   const BaseFab<T>& src,
00470                   const Box&        srcbox,
00471                   int               reversal_index,
00472                   T*                multiplier);
00473     //
00475     //
00476     T sum (int comp,
00477            int numcomp = 1) const;
00478 
00482     T sum (const Box& subbox,
00483            int        comp,
00484            int        numcomp = 1) const;
00485     //
00487     //
00488     BaseFab<T>& invert (T          v,
00489                         const Box& subbox,
00490                         int        comp=0,
00491                         int        numcomp=1);
00492     //
00494     //
00495     BaseFab<T>& invert (T   v,
00496                         int comp,
00497                         int numcomp=1);
00498     //
00500     //
00501     BaseFab<T>& invert (T v);
00502     //
00504     //
00505     BaseFab<T>& negate (const Box& subbox,
00506                         int        comp=0,
00507                         int        numcomp=1);
00508     //
00510     //
00511     BaseFab<T>& negate (int comp,
00512                         int numcomp=1);
00513     //
00515     //
00516     BaseFab<T>& negate ();
00517     //
00519     //
00520     BaseFab<T>& plus (T          r,
00521                       const Box& b,
00522                       int        comp=0,
00523                       int        numcomp=1);
00524     //
00526     //
00527     BaseFab<T>& plus (T   r,
00528                       int comp,
00529                       int numcomp=1);
00530     //
00532     //
00533     BaseFab<T>& plus (T r);
00534 
00538     BaseFab<T>& operator+= (T r);
00539     //
00541     //
00542     BaseFab<T>& operator+= (const BaseFab<T>& f);
00543     //
00545     //
00546     BaseFab<T>& plus (const BaseFab<T>& src);
00547 
00552     BaseFab<T>& plus (const BaseFab<T>& src,
00553                       int               srccomp,
00554                       int               destcomp,
00555                       int               numcomp=1);
00556 
00561     BaseFab<T>& plus (const BaseFab<T>& src,
00562                       const Box&        subbox,
00563                       int               srccomp,
00564                       int               destcomp,
00565                       int               numcomp=1);
00566 
00570     BaseFab<T>& plus (const BaseFab<T>& src,
00571                       const Box&        srcbox,
00572                       const Box&        destbox,
00573                       int               srccomp,
00574                       int               destcomp,
00575                       int               numcomp=1);
00576 
00580     BaseFab<T>& operator-= (T r);
00581     //
00583     //
00584     BaseFab<T>& operator-= (const BaseFab<T>& f);
00585     //
00587     //
00588     BaseFab<T>& minus (const BaseFab<T>& src);
00589 
00594     BaseFab<T>& minus (const BaseFab<T>& src,
00595                        int               srccomp,
00596                        int               destcomp,
00597                        int               numcomp=1);
00598 
00603     BaseFab<T>& minus (const BaseFab<T>& src,
00604                        const Box&        subbox,
00605                        int               srccomp,
00606                        int               destcomp,
00607                        int               numcomp=1);
00608 
00612     BaseFab<T>& minus (const BaseFab<T>& src,
00613                        const Box&        srcbox,
00614                        const Box&        destbox,
00615                        int               srccomp,
00616                        int               destcomp,
00617                        int               numcomp=1);
00618     //
00620     //
00621     BaseFab<T>& operator*= (T r);
00622     //
00624     //
00625     BaseFab<T>& mult (T r);
00626 
00630     BaseFab<T>& mult (T   r,
00631                       int comp,
00632                       int numcomp=1);
00633     //
00635     //
00636     BaseFab<T>& mult (T          r,
00637                       const Box& b,
00638                       int        comp=0,
00639                       int        numcomp=1);
00640     //
00642     //
00643     BaseFab<T>& operator*= (const BaseFab<T>& f);
00644     //
00646     //
00647     BaseFab<T>& mult (const BaseFab<T>& src);
00648 
00653     BaseFab<T>& mult (const BaseFab<T>& src,
00654                       int               srccomp,
00655                       int               destcomp,
00656                       int               numcomp=1);
00657 
00662     BaseFab<T>& mult (const BaseFab<T>& src,
00663                       const Box&        subbox,
00664                       int               srccomp,
00665                       int               destcomp,
00666                       int               numcomp=1);
00667 
00671     BaseFab<T>& mult (const BaseFab<T>& src,
00672                       const Box&        srcbox,
00673                       const Box&        destbox,
00674                       int               srccomp,
00675                       int               destcomp,
00676                       int               numcomp=1);
00677     //
00679     //
00680     BaseFab<T>& operator/= (T r);
00681     //
00683     //
00684     BaseFab<T>& divide (T r);
00685     //
00687     //
00688     BaseFab<T>& divide (T   r,
00689                         int comp,
00690                         int numcomp=1);
00691     //
00693     //
00694     BaseFab<T>& divide (T          r,
00695                         const Box& b,
00696                         int        comp=0,
00697                         int        numcomp=1);
00698     //
00700     //
00701     BaseFab<T>& operator/= (const BaseFab<T>& src);
00702     //
00704     //
00705     BaseFab<T>& divide (const BaseFab<T>& src);
00706 
00712     BaseFab<T>& divide (const BaseFab<T>& src,
00713                         int               srccomp,
00714                         int               destcomp,
00715                         int               numcomp=1);
00716 
00721     BaseFab<T>& divide (const BaseFab<T>& src,
00722                         const Box&        subbox,
00723                         int               srccomp,
00724                         int               destcomp,
00725                         int               numcomp=1);
00726  
00730     BaseFab<T>& divide (const BaseFab<T>& src,
00731                         const Box&        srcbox,
00732                         const Box&        destbox,
00733                         int               srccomp,
00734                         int               destcomp,
00735                         int               numcomp=1);
00736 
00746     BaseFab<T>& linInterp (const BaseFab<T>& f1,
00747                            const Box&        b1,
00748                            int               comp1,
00749                            const BaseFab<T>& f2,
00750                            const Box&        b2,
00751                            int               comp2,
00752                            Real              t1,
00753                            Real              t2,
00754                            Real              t,
00755                            const Box&        b,
00756                            int               comp,
00757                            int               numcomp = 1);
00758 
00767     BaseFab<T>& linComb (const BaseFab<T>& f1,
00768                          const Box&        b1,
00769                          int               comp1,
00770                          const BaseFab<T>& f2,
00771                          const Box&        b2,
00772                          int               comp2,
00773                          Real              alpha,
00774                          Real              beta,
00775                          const Box&        b,
00776                          int               comp,
00777                          int               numcomp = 1);
00778 
00779 protected:
00780     //
00781     // Allocates memory for the `BaseFab<T>'.
00782     //
00783     void define ();
00784     //
00785     // Deallocates memory for the `BaseFab<T>'.
00786     //
00787     void undefine ();
00788     //
00789     // The function called by BaseFab copy operations.
00790     //
00791     void performCopy (const BaseFab<T>& src,
00792                       const Box&        srcbox,
00793                       int               srccomp,
00794                       const Box&        destbox,
00795                       int               destcomp,
00796                       int               numcomp);
00797     //
00798     // This function is called by the `BaseFab' setVal operations.
00799     //
00800     void performSetVal (T          x,
00801                         const Box& bx,
00802                         int        nstart,
00803                         int        numcomp);
00804 protected:
00805 
00806     Box  domain;   // My index space.
00807     int  nvar;     // Number components.
00808     long numpts;   // Cached number of points in FAB.
00809     long truesize; // nvar*numpts that was allocated on heap.
00810     T*   dptr;     // The data pointer.
00811 };
00812 
00813 template <class T>
00814 inline
00815 int
00816 BaseFab<T>::nComp () const
00817 {
00818     return nvar;
00819 }
00820 
00821 template <class T>
00822 inline
00823 const Box&
00824 BaseFab<T>::box () const
00825 {
00826     return domain;
00827 }
00828 
00829 template <class T>
00830 inline
00831 const int*
00832 BaseFab<T>::length () const
00833 {
00834     return domain.length().getVect();
00835 }
00836 
00837 template <class T>
00838 inline
00839 const IntVect&
00840 BaseFab<T>::smallEnd () const
00841 {
00842     return domain.smallEnd();
00843 }
00844 
00845 template <class T>
00846 inline
00847 const IntVect&
00848 BaseFab<T>::bigEnd () const
00849 {
00850     return domain.bigEnd();
00851 }
00852 
00853 template <class T>
00854 inline
00855 const int*
00856 BaseFab<T>::loVect () const
00857 {
00858     return domain.loVect();
00859 }
00860 
00861 template <class T>
00862 inline
00863 const int*
00864 BaseFab<T>::hiVect () const
00865 {
00866     return domain.hiVect();
00867 }
00868 
00869 template <class T>
00870 bool
00871 BaseFab<T>::contains (const BaseFab<T>& fab) const
00872 {
00873     return box().contains(fab.box()) && nvar <= fab.nvar;
00874 }
00875 
00876 template <class T>
00877 bool
00878 BaseFab<T>::contains (const Box& bx) const
00879 {
00880     return box().contains(bx);
00881 }
00882 
00883 template <class T>
00884 inline
00885 T*
00886 BaseFab<T>::dataPtr (int n)
00887 {
00888     BL_ASSERT(!(dptr == 0));
00889     return &dptr[n*numpts];
00890 }
00891 
00892 template <class T>
00893 inline
00894 const T*
00895 BaseFab<T>::dataPtr (int n) const
00896 {
00897     BL_ASSERT(!(dptr == 0));
00898     return &dptr[n*numpts];
00899 }
00900 
00901 template <class T>
00902 inline
00903 bool
00904 BaseFab<T>::isAllocated () const
00905 {
00906     return dptr != 0;
00907 }
00908 
00909 template <class T>
00910 inline
00911 T&
00912 BaseFab<T>::operator() (const IntVect& p,
00913                         int            n)
00914 {
00915     BL_ASSERT(n >= 0);
00916     BL_ASSERT(n < nvar);
00917     BL_ASSERT(!(dptr == 0));
00918     BL_ASSERT(domain.contains(p));
00919     return dptr[domain.index(p)+n*numpts];
00920 }
00921 
00922 template <class T>
00923 inline
00924 T&
00925 BaseFab<T>::operator() (const IntVect& p)
00926 {
00927     BL_ASSERT(!(dptr == 0));
00928     BL_ASSERT(domain.contains(p));
00929 
00930     return dptr[domain.index(p)];
00931 }
00932 
00933 template <class T>
00934 inline
00935 const T&
00936 BaseFab<T>::operator() (const IntVect& p,
00937                         int            n) const
00938 {
00939     BL_ASSERT(n >= 0);
00940     BL_ASSERT(n < nvar);
00941     BL_ASSERT(!(dptr == 0));
00942     BL_ASSERT(domain.contains(p));
00943 
00944     return dptr[domain.index(p)+n*numpts];
00945 }
00946 
00947 template <class T>
00948 inline
00949 const T&
00950 BaseFab<T>::operator() (const IntVect& p) const
00951 {
00952     BL_ASSERT(!(dptr == 0));
00953     BL_ASSERT(domain.contains(p));
00954 
00955     return dptr[domain.index(p)];
00956 }
00957 
00958 template <class T>
00959 void
00960 BaseFab<T>::getVal  (T*             data,
00961                      const IntVect& pos,
00962                      int            n,
00963                      int            numcomp) const
00964 {
00965     const int loc      = domain.index(pos);
00966     const long size    = domain.numPts();
00967 
00968     BL_ASSERT(!(dptr == 0));
00969     BL_ASSERT(n >= 0 && n + numcomp <= nvar);
00970 
00971     for (int k = 0; k < numcomp; k++)
00972         data[k] = dptr[loc+(n+k)*size];
00973 }
00974 
00975 template <class T>
00976 void
00977 BaseFab<T>::getVal (T*             data,
00978                     const IntVect& pos) const
00979 {
00980     getVal(data,pos,0,nvar);
00981 }
00982 
00983 template <class T>
00984 BaseFab<T>&
00985 BaseFab<T>::shift (const IntVect& v)
00986 {
00987     domain += v;
00988     return *this;
00989 }
00990 
00991 template <class T>
00992 BaseFab<T>&
00993 BaseFab<T>::shift (int idir,
00994                    int n_cell)
00995 {
00996     domain.shift(idir,n_cell);
00997     return *this;
00998 }
00999 
01000 template <class T>
01001 BaseFab<T> &
01002 BaseFab<T>::shiftHalf (const IntVect& v)
01003 {
01004     domain.shiftHalf(v);
01005     return *this;
01006 }
01007 
01008 template <class T>
01009 BaseFab<T> &
01010 BaseFab<T>::shiftHalf (int idir,
01011                        int n_cell)
01012 {
01013     domain.shiftHalf(idir,n_cell);
01014     return *this;
01015 }
01016 
01017 template <class T>
01018 void
01019 BaseFab<T>::setVal (T val)
01020 {
01021     performSetVal(val,box(), 0, nvar);
01022 }
01023 
01024 template <class T>
01025 void
01026 BaseFab<T>::setVal (T          x,
01027                     const Box& bx,
01028                     int        n)
01029 {
01030     performSetVal(x,bx,n,1);
01031 }
01032 
01033 template <class T>
01034 void
01035 BaseFab<T>::setVal (T   x,
01036                     int n)
01037 {
01038     performSetVal(x,domain,n,1);
01039 }
01040 
01041 template <class T>
01042 void
01043 BaseFab<T>::setVal (T          x,
01044                     const Box& b,
01045                     int        ns,
01046                     int        num)
01047 {
01048     performSetVal(x,b,ns,num);
01049 }
01050 
01051 template <class T>
01052 BaseFab<T>&
01053 BaseFab<T>::copy (const BaseFab<T>& src,
01054                   const Box&        srcbox,
01055                   int               srccomp,
01056                   const Box&        destbox,
01057                   int               destcomp,
01058                   int               numcomp)
01059 {
01060     if (!srcbox.sameSize(destbox))
01061      BoxLib::Error("srcbox<>destbox in fab copy");
01062     if (!src.box().contains(srcbox))
01063      BoxLib::Error("srcbox bigger then FAB src box");
01064     if (!domain.contains(destbox))
01065      BoxLib::Error("destbox bigger then FAB dest box");
01066     if ((srccomp<0)||(srccomp+numcomp>src.nComp()))
01067      BoxLib::Error("srccomp+numcomp out of bounds");
01068     if ((destcomp<0)||(destcomp+numcomp>nvar))
01069      BoxLib::Error("destcomp+numcomp out of bounds");
01070 
01071     BL_ASSERT(srcbox.sameSize(destbox));
01072     BL_ASSERT(src.box().contains(srcbox));
01073     BL_ASSERT(domain.contains(destbox));
01074     BL_ASSERT(srccomp >= 0 && srccomp+numcomp <= src.nComp());
01075     BL_ASSERT(destcomp >= 0 && destcomp+numcomp <= nvar);
01076     performCopy(src,srcbox,srccomp,destbox,destcomp,numcomp);
01077     return *this;
01078 }
01079 
01080 template <class T>
01081 BaseFab<T>&
01082 BaseFab<T>::copy (const BaseFab<T>& src)
01083 {
01084     BL_ASSERT(nvar <= src.nvar);
01085     BL_ASSERT(domain.sameType(src.domain));
01086     Box overlap(domain);
01087     overlap &= src.domain;
01088     if (overlap.ok())
01089         performCopy(src,overlap,0,overlap,0,nvar);
01090     return *this;
01091 }
01092 
01093 template <class T>
01094 BaseFab<T>&
01095 BaseFab<T>::copy (const BaseFab<T>& src,
01096                   const Box&        destbox)
01097 {
01098     BL_ASSERT(nvar <= src.nvar);
01099     BL_ASSERT(domain.contains(destbox));
01100     Box overlap(destbox);
01101     overlap &= src.domain;
01102     if (overlap.ok())
01103         performCopy(src,overlap,0,overlap,0,nvar);
01104     return *this;
01105 }
01106 
01107 template <class T>
01108 BaseFab<T>&
01109 BaseFab<T>::copy (const BaseFab<T>& src,
01110                   int               srccomp,
01111                   int               destcomp,
01112                   int               numcomp)
01113 {
01114     BL_ASSERT(srccomp >= 0 && srccomp + numcomp <= src.nvar);
01115     BL_ASSERT(destcomp >= 0 && destcomp + numcomp <= nvar);
01116     Box overlap(domain);
01117     overlap &= src.domain;
01118     if (overlap.ok())
01119         performCopy(src,overlap,srccomp,overlap,destcomp,numcomp);
01120     return *this;
01121 }
01122 
01123 template <class T>
01124 void
01125 BaseFab<T>::define ()
01126 {
01127     BL_ASSERT(nvar > 0);
01128     BL_ASSERT(dptr == 0);
01129     BL_ASSERT(numpts > 0);
01130 
01131     truesize = nvar*numpts;
01132     dptr     = static_cast<T*>(BoxLib::The_Arena()->alloc(truesize*sizeof(T)));
01133     //
01134     // Now call T::T() on the raw memory so we have valid Ts.
01135     //
01136     T* ptr = dptr;
01137 
01138     for (int i = 0; i < truesize; i++, ptr++)
01139         new (ptr) T;
01140 }
01141 
01142 template <class T>
01143 void
01144 BaseFab<T>::undefine ()
01145 {
01146     //
01147     // Call T::~T() on the to-be-destroyed memory.
01148     //
01149     T* ptr = dptr;
01150 
01151     for (int i = 0; i < truesize; i++, ptr++)
01152     {
01153         ptr->~T();
01154     }
01155     BoxLib::The_Arena()->free(dptr);
01156 
01157     dptr = 0;
01158 }
01159 
01160 template <class T>
01161 BaseFab<T>::BaseFab ()
01162     :
01163     domain(Box()),
01164     nvar(0),
01165     numpts(0),
01166     truesize(0),
01167     dptr(0)
01168 {}
01169 
01170 template <class T>
01171 BaseFab<T>::BaseFab (const BaseFab& fab)
01172     : domain(fab.box()),
01173       nvar(fab.nComp())
01174 {
01175     numpts = domain.numPts();
01176     dptr = 0;
01177     define();
01178     for ( int i = 0; i < numpts*nvar; ++i )
01179     {
01180     dptr[i] = fab.dptr[i];
01181     }
01182 }
01183 
01184 template <class T>
01185 BaseFab<T>&
01186 BaseFab<T>::operator= (const BaseFab& fab)
01187 {
01188     if ( this == &fab ) return *this;
01189     domain = fab.box();
01190     nvar   = fab.nComp();
01191     numpts = domain.numPts();
01192     dptr = 0;
01193     define();
01194     for ( int i = 0; i < numpts*nvar; ++i )
01195     {
01196     dptr[i] = fab.dptr[i];
01197     }
01198     return *this;
01199 }
01200 
01201 template <class T>
01202 BaseFab<T>&
01203 BaseFab<T>::operator= (const T& t)
01204 {
01205     setVal(t);
01206     return *this;
01207 }
01208 
01209 template <class T>
01210 BaseFab<T>::BaseFab (const Box& bx,
01211                      int        n)
01212     :
01213     domain(bx),
01214     nvar(n),
01215     numpts(bx.numPts()),
01216     dptr(0)
01217 {
01218     define();
01219 }
01220 
01221 template <class T>
01222 void
01223 BaseFab<T>::resize (const Box& b,
01224                     int        n)
01225 {
01226     nvar   = n;
01227     domain = b;
01228     numpts = domain.numPts();
01229 
01230     if (dptr == 0)
01231     {
01232         define();
01233     }
01234     else if (nvar*numpts > truesize)
01235     {
01236         undefine();
01237 
01238         define();
01239     }
01240 }
01241 
01242 template <class T>
01243 BaseFab<T>::~BaseFab ()
01244 {
01245     undefine();
01246 }
01247 
01248 template <class T>
01249 void
01250 BaseFab<T>::clear ()
01251 {
01252     undefine();
01253 
01254     dptr   = 0;
01255     domain = Box();
01256     nvar   = 0;
01257     numpts = 0;
01258 }
01259 
01260 //
01261 // performCopy() has been rewritten here so we can insert pragma's which
01262 // will enhance vectorization on the Cray's.  The downside is that the
01263 // code is greatly expanded and rather incomprehensible.
01264 //
01265 
01266 template <class T>
01267 void
01268 BaseFab<T>::performCopy (const BaseFab<T>& src,
01269                          const Box&        srcbox,
01270                          int               srccomp,
01271                          const Box&        destbox,
01272                          int               destcomp,
01273                          int               numcomp)
01274 {
01275     if (!srcbox.sameSize(destbox))
01276      BoxLib::Error("srcbox<>destbox in performCopy");
01277     if (!src.box().contains(srcbox))
01278      BoxLib::Error("srcbox bigger then FAB src box performCopy");
01279     if (!box().contains(destbox))
01280      BoxLib::Error("destbox bigger then FAB dest box performCopy");
01281     if ((srccomp<0)||(srccomp+numcomp>src.nComp()))
01282      BoxLib::Error("srccomp+numcomp out of bounds performCopy");
01283     if ((destcomp<0)||(destcomp+numcomp>nComp()))
01284      BoxLib::Error("destcomp+numcomp out of bounds performCopy");
01285 
01286 
01287 #if (BL_SPACEDIM == 1)
01288 {                                                                       
01289     BL_ASSERT((destcomp) >= 0 && (destcomp) + (numcomp) <= nComp());
01290     BL_ASSERT((srccomp) >= 0 && (srccomp) + (numcomp) <= (src).nComp());
01291     Box _subbox_ = box(); 
01292     _subbox_ &= destbox; 
01293     BL_ASSERT(srcbox.sameSize(_subbox_)); 
01294     if (_subbox_.ok()) 
01295     { 
01296         const int *_th_plo = loVect(); 
01297         const int *_th_plen = length(); 
01298         const int *_x_plo = (src).loVect(); 
01299         const int *_x_plen = (src).length(); 
01300         const int *_subbox_lo = _subbox_.loVect(); 
01301         const int *_subbox_len = _subbox_.length().getVect(); 
01302         const int *_bx_lo = (srcbox).loVect(); 
01303         // const int *_bx_len = (srcbox).length().getVect(); 
01304         T* _th_p = dataPtr(destcomp); 
01305         const T* _x_p  = (src).dataPtr(srccomp); 
01306         for (int _n = 0; _n < (numcomp); ++_n)
01307         { 
01308             T *_th_pp = _th_p + ((_subbox_lo[0]-_th_plo[0])+_n*_th_plen[0]); 
01309             const T *_x_pp = _x_p + ((_bx_lo[0]-_x_plo[0])+_n*_x_plen[0]);
01310 #ifdef BL_ARCH_CRAY
01311 #pragma _CRI ivdep
01312 #endif
01313             for (int _i = 0; _i < _subbox_len[0]; ++_i, ++_th_pp)
01314             { 
01315                 int iR = _i + _subbox_lo[0]; iR += 0; 
01316                 int ixR = _i + _bx_lo[0]; ixR += 0; 
01317                 T &thisR = * _th_pp; const T & srcR = _x_pp[_i];
01318 #elif (BL_SPACEDIM == 2)
01319 {                                                                       
01320     BL_ASSERT((destcomp) >= 0 && (destcomp) + (numcomp) <= nComp()); 
01321     BL_ASSERT((srccomp) >= 0 && (srccomp) + (numcomp) <= (src).nComp()); 
01322     Box _subbox_ = box(); 
01323     _subbox_ &= destbox; 
01324     BL_ASSERT(srcbox.sameSize(_subbox_)); 
01325     if (_subbox_.ok())
01326     { 
01327         const int *_th_plo = loVect(); 
01328         const int *_th_plen = length(); 
01329         const int *_x_plo = (src).loVect(); 
01330         const int *_x_plen = (src).length(); 
01331         const int *_subbox_lo = _subbox_.loVect(); 
01332         const int *_subbox_len = _subbox_.length().getVect(); 
01333         const int *_bx_lo = (srcbox).loVect(); 
01334         // const int *_bx_len = (srcbox).length().getVect(); 
01335         T* _th_p = dataPtr(destcomp); 
01336         const T* _x_p  = (src).dataPtr(srccomp); 
01337         for (int _n = 0; _n < (numcomp); ++_n)
01338         { 
01339             int nR = _n + destcomp; nR += 0; 
01340             int nxR = _n + srccomp; nxR += 0; 
01341             for(int _j = 0; _j < _subbox_len[1]; ++_j)
01342             { 
01343                 const int jR = _j + _subbox_lo[1]; 
01344                 const int jxR = _j + _bx_lo[1]; 
01345                 T *_th_pp = _th_p + ((_subbox_lo[0] - _th_plo[0]) 
01346                                      + _th_plen[0]*((jR - _th_plo[1]) 
01347                                                     + _n * _th_plen[1])); 
01348                 const T *_x_pp = _x_p + ((_bx_lo[0] - _x_plo[0]) 
01349                                          + _x_plen[0]*((jxR - _x_plo[1]) 
01350                                                        + _n * _x_plen[1])); 
01351 #ifdef BL_ARCH_CRAY
01352 #     pragma _CRI ivdep
01353 #endif
01354                 for (int _i = 0; _i < _subbox_len[0]; ++_i, ++_th_pp)
01355                 {  
01356                     T &thisR = * _th_pp; const T & srcR = _x_pp[_i]; 
01357 #elif (BL_SPACEDIM == 3)
01358 { 
01359     BL_ASSERT((destcomp) >= 0 && (destcomp) + (numcomp) <= nComp()); 
01360     BL_ASSERT((srccomp) >= 0 && (srccomp) + (numcomp) <= (src).nComp()); 
01361     Box _subbox_(box()); 
01362     _subbox_ &= destbox; 
01363     BL_ASSERT((srcbox).sameSize(_subbox_)); 
01364     if (_subbox_.ok())
01365     { 
01366         const int *_th_plo = loVect(); 
01367         const int *_th_plen = length(); 
01368         const int *_x_plo = (src).loVect(); 
01369         const int *_x_plen = (src).length(); 
01370         const int *_subbox_lo = _subbox_.loVect(); 
01371         const int *_subbox_len = _subbox_.length().getVect(); 
01372         const int *_bx_lo = (srcbox).loVect(); 
01373         // const int *_bx_len = (srcbox).length().getVect(); 
01374         T* _th_p = dataPtr(destcomp); 
01375         const T* _x_p  = (src).dataPtr(srccomp); 
01376         for (int _n = 0; _n < (numcomp); ++_n)
01377         { 
01378             for (int _k = 0; _k < _subbox_len[2]; ++_k)
01379             { 
01380                 const int kR = _k + _subbox_lo[2]; 
01381                 const int kxR = _k + _bx_lo[2]; 
01382                 for(int _j = 0; _j < _subbox_len[1]; ++_j)
01383                 { 
01384                     const int jR = _j + _subbox_lo[1]; 
01385                     const int jxR = _j + _bx_lo[1]; 
01386                     T *_th_pp = _th_p + ((_subbox_lo[0] - _th_plo[0]) 
01387                                          + _th_plen[0]*((jR - _th_plo[1]) 
01388                                                         + _th_plen[1]*(
01389                                                             (kR - _th_plo[2]) 
01390                                                         + _n * _th_plen[2]))); 
01391                     const T *_x_pp = _x_p + ((_bx_lo[0] - _x_plo[0]) 
01392                                              + _x_plen[0]*((jxR - _x_plo[1]) 
01393                                                            + _x_plen[1]*( 
01394                                                             (kxR - _x_plo[2])
01395                                                          + _n * _x_plen[2])));
01396 #ifdef BL_ARCH_CRAY
01397 #     pragma _CRI ivdep
01398 #endif
01399                     for (int _i = 0; _i < _subbox_len[0]; ++_i, ++_th_pp)
01400                     {
01401                         T &thisR = * _th_pp; const T & srcR = _x_pp[_i]; 
01402 #endif
01403     {
01404         thisR = srcR;
01405     }
01406 #if (BL_SPACEDIM == 1)
01407      }}}}
01408 #elif (BL_SPACEDIM == 2)
01409      }}}}}
01410 #elif (BL_SPACEDIM == 3)
01411      }}}}}}
01412 #endif
01413 }
01414 
01415 template <class T>
01416 void
01417 BaseFab<T>::performSetVal (T         val,
01418                            const Box& bx,
01419                            int        ns,
01420                            int        num)
01421 {
01422     BL_ASSERT(domain.contains(bx));
01423     BL_ASSERT(ns >= 0 && ns + num <= nvar);
01424 
01425     if (bx == domain)
01426     {
01427         T* data = &dptr[ns*numpts];
01428         for (long i = 0, N = num*numpts; i < N; i++)
01429         {
01430             *data++ = val;
01431         }
01432     }
01433     else
01434     {
01435         ForAllThisBNN(T,bx,ns,num)
01436         {
01437             thisR = val;
01438         } EndFor
01439     }
01440 }
01441 
01442 template <class T>
01443 void
01444 BaseFab<T>::setComplement (T          x,
01445                            const Box& b,
01446                            int        ns,
01447                            int        num)
01448 {
01449     BoxList b_lst = BoxLib::boxDiff(domain,b);
01450     for (BoxList::iterator bli = b_lst.begin(); bli != b_lst.end(); ++bli)
01451         performSetVal(x, *bli, ns, num);
01452 }
01453 
01454 template <class T>
01455 void
01456 BaseFab<T>::abs ()
01457 {
01458     ForAllThis(T)
01459     {
01460         thisR = std::abs(thisR);
01461     } EndFor
01462 }
01463 
01464 template <class T>
01465 void
01466 BaseFab<T>::abs (int comp,
01467                  int numcomp)
01468 {
01469     ForAllThisNN(T,comp,numcomp)
01470     {
01471         thisR = std::abs(thisR);
01472     } EndFor
01473 }
01474 
01475 template <class T>
01476 void
01477 BaseFab<T>::abs (const Box& subbox,
01478                  int        comp,
01479                  int        numcomp)
01480 {
01481     ForAllThisBNN(T,subbox,comp,numcomp)
01482     {
01483         thisR = std::abs(thisR);
01484     } EndFor
01485 }
01486 
01487 template <class T>
01488 Real
01489 BaseFab<T>::norm (int p,
01490                   int comp,
01491                   int numcomp) const
01492 {
01493     return norm(domain,p,comp,numcomp);
01494 }
01495 
01496 template <class T>
01497 Real
01498 BaseFab<T>::norm (const Box& subbox,
01499                   int        p,
01500                   int        comp,
01501                   int        numcomp) const
01502 {
01503     BL_ASSERT(comp >= 0 && comp+numcomp <= nComp());
01504     BL_ASSERT(p >= 0);
01505 
01506     Real  nrm    = 0;
01507     Real* tmp    = 0;
01508     int   tmplen = 0;
01509 
01510     if (p == 0)
01511     {
01512         ForAllThisCPencil(T,subbox,comp,numcomp)
01513         {
01514             const T* row = &thisR;
01515             if (tmp == 0)
01516             {
01517                 tmp    = new Real[thisLen];
01518                 tmplen = thisLen;
01519                 for (int i = 0; i < thisLen; i++)
01520                     tmp[i] = std::abs(row[i]);
01521             }
01522             else
01523             {
01524                 for (int i = 0; i < thisLen; i++) {
01525                 Real kluge=std::abs(row[i]);
01526                     tmp[i] = std::max(tmp[i],kluge);
01527         }
01528             }
01529         } EndForPencil
01530         nrm = tmp[0];
01531         for (int i = 1; i < tmplen; i++)
01532             nrm = std::max(nrm, tmp[i]);
01533     }
01534     else if (p == 1)
01535     {
01536         ForAllThisCPencil(T,subbox,comp,numcomp)
01537         {
01538             const T* row = &thisR;
01539             if (tmp == 0)
01540             {
01541                 tmp    = new Real[thisLen];
01542                 tmplen = thisLen;
01543                 for (int i = 0; i < thisLen; i++)
01544                     tmp[i] = std::abs(row[i]);
01545             }
01546             else
01547             {
01548                 for (int i = 0; i < thisLen; i++)
01549                     tmp[i] += std::abs(row[i]);
01550             }
01551         } EndForPencil
01552         nrm = tmp[0];
01553         for (int i = 1; i < tmplen; i++)
01554             nrm += tmp[i];
01555     }
01556     else
01557     {
01558       BoxLib::Error("BaseFab::norm(): only p == 0 or p == 1 are supported");
01559     }
01560 
01561     delete [] tmp;
01562 
01563     return nrm;
01564 }
01565 
01566 
01567 #ifdef BL_ARCH_CRAY
01568 #define RESTRICT restrict
01569 //
01570 // Template specialization for Real.
01571 //
01572 template<>
01573 Real
01574 BaseFab<Real>::norm (const Box& subbox,
01575                      int        p,
01576                      int        comp,
01577                      int        numcomp) const
01578 {
01579     BL_ASSERT(comp >= 0 && comp+numcomp <= nComp());
01580     BL_ASSERT(p >= 0);
01581 
01582     Real* RESTRICT tmp  = 0;
01583     int tmplen = 0;
01584     Real nrm   = 0;
01585     if (p == 0)
01586     {
01587         ForAllThisCPencil(Real,subbox,comp,numcomp)
01588         {
01589             const Real* RESTRICT row = &thisR;
01590             if (tmp == 0)
01591             {
01592                 tmp = new Real[thisLen];
01593                 tmplen = thisLen;
01594 #pragma _CRI ivdep
01595                 for (int i = 0; i < thisLen; i++)
01596                     tmp[i] = fabs(row[i]);
01597             }
01598             else
01599             {
01600 #pragma _CRI ivdep
01601               for (int i = 0; i < thisLen; i++)
01602               {
01603                 Real a = fabs(row[i]);
01604                 tmp[i] = tmp[i] > a ? tmp[i] : a ;
01605               }
01606             }
01607         } EndForPencil
01608         nrm = tmp[0];
01609         for (int i = 1; i < tmplen; i++)
01610         {
01611             Real a = tmp[i];
01612             nrm = nrm > a ? nrm : a ;
01613         }
01614     }
01615     else if (p == 1)
01616     {
01617         ForAllThisCPencil(Real,subbox,comp,numcomp)
01618         {
01619             const Real* row = &thisR;
01620             if (tmp == 0)
01621             {
01622                 tmp = new Real[thisLen];
01623                 tmplen = thisLen;
01624 #pragma _CRI ivdep
01625                 for (int i = 0; i < thisLen; i++)
01626                     tmp[i] = fabs(row[i]);
01627             }
01628             else
01629             {
01630 #pragma _CRI ivdep
01631               for (int i = 0; i < thisLen; i++)
01632               {
01633                     tmp[i] += fabs( row[i] );
01634               }
01635             }
01636         } EndForPencil
01637         nrm = tmp[0];
01638 #pragma _CRI ivdep
01639         for (int i = 1; i < tmplen; i++)
01640             nrm += tmp[i];
01641     }
01642     else
01643       BoxLib::Error("BaseFab<Real>::norm(): only p == 0 or p == 1 are supported");
01644 
01645     delete [] tmp;
01646 
01647     return nrm;
01648 }
01649 #endif
01650 
01651 template <class T>
01652 T
01653 BaseFab<T>::min (int comp) const
01654 {
01655     T*  _min_row = 0;
01656     int _X_len   = 0;
01657     ForAllThisCPencil(T,domain,comp,1)
01658     {
01659         const T* _row = &thisR;
01660         if (_min_row == 0)
01661         {
01662             _min_row = new T[thisLen];
01663             _X_len = thisLen;
01664             for (int i = 0; i < thisLen; i++)
01665                 _min_row[i] = _row[i];
01666         }
01667         else
01668         {
01669             for (int i = 0; i < thisLen; i++)
01670                 _min_row[i] = std::min(_row[i],_min_row[i]);
01671         }
01672     } EndForPencil;
01673 
01674     T _min = _min_row[0];
01675     for (int i = 1; i < _X_len; i++)
01676         _min = std::min(_min,_min_row[i]);
01677 
01678     delete [] _min_row;
01679 
01680     return _min;
01681 }
01682 
01683 template <class T>
01684 T
01685 BaseFab<T>::min (const Box& subbox,
01686                  int        comp) const
01687 {
01688     T *_min_row = 0;
01689     int _X_len = 0;
01690     ForAllThisCPencil(T,subbox,comp,1)
01691     {
01692         const T* _row = &thisR;
01693         if (_min_row == 0)
01694         {
01695             _min_row = new T[thisLen];
01696             _X_len = thisLen;
01697             for (int i = 0; i < thisLen; i++)
01698                 _min_row[i] = _row[i];
01699         }
01700         else
01701         {
01702             for (int i = 0; i < thisLen; i++)
01703                 _min_row[i] = std::min(_row[i],_min_row[i]);
01704         }
01705     } EndForPencil;
01706 
01707     T _min = _min_row[0];
01708     for (int i = 1; i < _X_len; i++)
01709         _min = std::min(_min,_min_row[i]);
01710 
01711     delete [] _min_row;
01712 
01713     return _min;
01714 }
01715 
01716 template <class T>
01717 T
01718 BaseFab<T>::max (int comp) const
01719 {
01720     T* _max_row = 0;
01721     int _X_len  = 0;
01722     ForAllThisCPencil(T,domain,comp,1)
01723     {
01724         const T* _row = &thisR;
01725         if (_max_row== 0)
01726         {
01727             _max_row = new T[thisLen];
01728             _X_len = thisLen;
01729             for (int i = 0; i < thisLen; i++)
01730                 _max_row[i] = _row[i];
01731         }
01732         else
01733         {
01734             for (int i = 0; i < thisLen; i++)
01735                 _max_row[i] = std::max(_row[i],_max_row[i]);
01736         }
01737     } EndForPencil;
01738 
01739     T _max = _max_row[0];
01740     for (int i = 1; i < _X_len; i++)
01741         _max = std::max(_max,_max_row[i]);
01742 
01743     delete [] _max_row;
01744 
01745     return _max;
01746 }
01747 
01748 template <class T>
01749 T
01750 BaseFab<T>::max (const Box& subbox,
01751                  int        comp) const
01752 {
01753     T*  _max_row = 0;
01754     int _X_len   = 0;
01755     ForAllThisCPencil(T,subbox,comp,1)
01756     {
01757         const T* _row = &thisR;
01758         if (_max_row == 0)
01759         {
01760             _max_row = new T[thisLen];
01761             _X_len = thisLen;
01762             for (int i = 0; i < thisLen; i++)
01763                 _max_row[i] = _row[i];
01764         }
01765         else
01766         {
01767             for (int i = 0; i < thisLen; i++)
01768                 _max_row[i] = std::max(_row[i],_max_row[i]);
01769         }
01770     } EndForPencil;
01771 
01772     T _max = _max_row[0];
01773     for (int i = 1; i < _X_len; i++)
01774         _max = std::max(_max,_max_row[i]);
01775 
01776     delete [] _max_row;
01777 
01778     return _max;
01779 }
01780 
01781 template <class T>
01782 IntVect
01783 BaseFab<T>::minIndex (int comp) const
01784 {
01785     IntVect _min_loc(domain.smallEnd());
01786     T _min_val = (*this).operator()(_min_loc,comp);
01787     ForAllThisCBNN(T,domain,comp,1)
01788     {
01789         if (thisR < _min_val)
01790         {
01791             _min_val = thisR;
01792             D_EXPR(_min_loc[0] = iR,
01793                    _min_loc[1] = jR,
01794                    _min_loc[2] = kR);
01795         }
01796     } EndFor;
01797 
01798     return _min_loc;
01799 }
01800 
01801 template <class T>
01802 IntVect
01803 BaseFab<T>::minIndex (const Box& subbox,
01804                       int        comp) const
01805 {
01806     IntVect _min_loc(subbox.smallEnd());
01807     T _min_val = (*this).operator()(_min_loc,comp);
01808     ForAllThisCBNN(T,subbox,comp,1)
01809     {
01810         if (thisR < _min_val)
01811         {
01812             _min_val = thisR;
01813             D_EXPR(_min_loc[0] = iR,
01814                    _min_loc[1] = jR,
01815                    _min_loc[2] = kR);
01816         }
01817     } EndFor;
01818 
01819     return _min_loc;
01820 }
01821 
01822 template <class T>
01823 IntVect
01824 BaseFab<T>::maxIndex (int comp) const
01825 {
01826     IntVect _max_loc(domain.smallEnd());
01827     T _max_val = (*this).operator()(_max_loc,comp);
01828     ForAllThisCBNN(T,domain,comp,1)
01829     {
01830         if (thisR > _max_val)
01831         {
01832             _max_val = thisR;
01833             D_EXPR(_max_loc[0] = iR,
01834                    _max_loc[1] = jR,
01835                    _max_loc[2] = kR);
01836         }
01837     } EndFor;
01838 
01839     return _max_loc;
01840 }
01841 
01842 template <class T>
01843 IntVect
01844 BaseFab<T>::maxIndex (const Box& subbox,
01845                       int        comp) const
01846 {
01847     IntVect _max_loc(subbox.smallEnd());
01848     T _max_val = (*this).operator()(_max_loc,comp);
01849     ForAllThisCBNN(T,subbox,comp,1)
01850     {
01851         if (thisR > _max_val)
01852         {
01853             _max_val = thisR;
01854             D_EXPR(_max_loc[0] = iR,
01855                    _max_loc[1] = jR,
01856                    _max_loc[2] = kR);
01857         }
01858     } EndFor;
01859 
01860     return _max_loc;
01861 }
01862 
01863 template <class T>
01864 int
01865 BaseFab<T>::maskLT (BaseFab<int>& mask,
01866                     T             val,
01867                     int           comp) const
01868 {
01869     mask.resize(domain,1);
01870     mask.setVal(0);
01871 
01872     int* mptr = mask.dataPtr();
01873     int  cnt  = 0;
01874 
01875     ForAllThisCBNN(T,domain,comp,1)
01876     {
01877         int ix = D_TERM(_i, +_j*_b_len[0], +_k*_b_len[0]*_b_len[1]);
01878         if (thisR < val)
01879         {
01880             mptr[ix] = 1;
01881             cnt++;
01882         }
01883     } EndFor;
01884 
01885     return cnt;
01886 }
01887 
01888 template <class T>
01889 int
01890 BaseFab<T>::maskLE (BaseFab<int>& mask,
01891                     T             val,
01892                     int           comp) const
01893 {
01894     mask.resize(domain,1);
01895     mask.setVal(0);
01896 
01897     int* mptr = mask.dataPtr();
01898     int  cnt  = 0;
01899 
01900     ForAllThisCBNN(T,domain,comp,1)
01901     {
01902         int ix = D_TERM(_i, +_j*_b_len[0], +_k*_b_len[0]*_b_len[1]);
01903         if (thisR <= val)
01904         {
01905             mptr[ix] = 1;
01906             cnt++;
01907         }
01908     } EndFor;
01909 
01910     return cnt;
01911 }
01912 
01913 template <class T>
01914 int
01915 BaseFab<T>::maskEQ (BaseFab<int>& mask,
01916                     T             val,
01917                     int           comp) const
01918 {
01919     mask.resize(domain,1);
01920     mask.setVal(0);
01921 
01922     int* mptr = mask.dataPtr();
01923     int  cnt  = 0;
01924 
01925     ForAllThisCBNN(T,domain,comp,1)
01926     {
01927         int ix = D_TERM(_i, +_j*_b_len[0], +_k*_b_len[0]*_b_len[1]);
01928         if (thisR == val)
01929         {
01930             mptr[ix] = 1;
01931             cnt++;
01932         }
01933     } EndFor;
01934 
01935     return cnt;
01936 }
01937 
01938 template <class T>
01939 int
01940 BaseFab<T>::maskGT (BaseFab<int>& mask,
01941                     T             val,
01942                     int           comp) const
01943 {
01944     mask.resize(domain,1);
01945     mask.setVal(0);
01946 
01947     int* mptr = mask.dataPtr();
01948     int  cnt  = 0;
01949 
01950     ForAllThisCBNN(T,domain,comp,1)
01951     {
01952         int ix = D_TERM(_i, +_j*_b_len[0], +_k*_b_len[0]*_b_len[1]);
01953         if (thisR > val)
01954         {
01955             mptr[ix] = 1;
01956             cnt++;
01957         }
01958     } EndFor;
01959 
01960     return cnt;
01961 }
01962 
01963 template <class T>
01964 int
01965 BaseFab<T>::maskGE(BaseFab<int>& mask,
01966                    T             val,
01967                    int           comp) const
01968 {
01969     mask.resize(domain,1);
01970     mask.setVal(0);
01971 
01972     int* mptr = mask.dataPtr();
01973     int  cnt  = 0;
01974 
01975     ForAllThisCBNN(T,domain,comp,1)
01976     {
01977         int ix = D_TERM(_i, +_j*_b_len[0], +_k*_b_len[0]*_b_len[1]);
01978         if (thisR >= val)
01979         {
01980             mptr[ix] = 1;
01981             cnt++;
01982         }
01983     } EndFor;
01984 
01985     return cnt;
01986 }
01987 
01988 template <class T>
01989 BaseFab<T>&
01990 BaseFab<T>::plus (T r)
01991 {
01992     return operator+=(r);
01993 }
01994 
01995 template <class T>
01996 BaseFab<T>&
01997 BaseFab<T>::plus (const BaseFab<T>& x)
01998 {
01999     return operator+=(x);
02000 }
02001 
02002 template <class T>
02003 BaseFab<T>&
02004 BaseFab<T>::operator-= (T r)
02005 {
02006     return operator+=(-r);
02007 }
02008 
02009 template <class T>
02010 BaseFab<T>&
02011 BaseFab<T>::minus (const BaseFab<T>& x)
02012 {
02013     return operator-=(x);
02014 }
02015 
02016 
02017 template <class T>
02018 BaseFab<T>&
02019 BaseFab<T>::mult (T r)
02020 {
02021     return operator*=(r);
02022 }
02023 
02024 template <class T>
02025 BaseFab<T>&
02026 BaseFab<T>::mult (const BaseFab<T>& x)
02027 {
02028     return operator*=(x);
02029 }
02030 
02031 template <class T>
02032 BaseFab<T>&
02033 BaseFab<T>::divide (T r)
02034 {
02035     return operator/=(r);
02036 }
02037 
02038 template <class T>
02039 BaseFab<T>&
02040 BaseFab<T>::divide (const BaseFab<T>& x)
02041 {
02042     return operator/=(x);
02043 }
02044 
02045 
02046 template <class T>
02047 void
02048 BaseFab<T>::patternFill (int mark)
02049 {
02050     ForAllThis(T)
02051     {
02052         thisR = D_TERM(iR*100, +jR*10, + kR) + 1000*nR + 10000*mark;
02053     } EndFor
02054 }
02055 
02056 template <class T>
02057 void
02058 BaseFab<T>::copyRev (const Box&        destbox,
02059                      const BaseFab<T>& src,
02060                      const Box&        srcbox,
02061                      int               reversal_index,
02062                      T*                multiplier)
02063 {
02064     BaseFab<T>& dest = *this;
02065 
02066     ForAllRevXBNYCBNNN(T,dest,destbox,0,src,srcbox,0,nComp(),reversal_index)
02067     {
02068         destR = multiplier[_n]*srcR;
02069     } EndFor
02070 }
02071 
02072 template <class T>
02073 T
02074 BaseFab<T>::sum (int comp,
02075                  int numcomp) const
02076 {
02077     T*  _sum_row = 0;
02078     int _sum_len = 0;
02079     ForAllThisCPencil(T,domain,comp,numcomp)
02080     {
02081         const T* _row = &thisR;
02082         if (_sum_row == 0)
02083         {
02084             _sum_row = new T[thisLen];
02085             _sum_len = thisLen;
02086             for (int i = 0; i < thisLen; i++)
02087                 _sum_row[i] = _row[i];
02088         }
02089         else
02090         {
02091             for (int i = 0; i < thisLen; i++)
02092                 _sum_row[i] += _row[i];
02093         }
02094     } EndForPencil;
02095 
02096     T _sum = _sum_row[0];
02097     for (int i = 1; i < _sum_len; i++)
02098         _sum += _sum_row[i];
02099 
02100     delete [] _sum_row;
02101 
02102     return _sum;
02103 }
02104 
02105 template <class T>
02106 T
02107 BaseFab<T>::sum (const Box& subbox,
02108                  int        comp,
02109                  int        numcomp) const
02110 {
02111     T*  _sum_row = 0;
02112     int _sum_len = 0;
02113     ForAllThisCPencil(T,subbox,comp,numcomp)
02114     {
02115         const T* _row = &thisR;
02116         if (_sum_row == 0)
02117         {
02118             _sum_row = new T[thisLen];
02119             _sum_len = thisLen;
02120             for (int i = 0; i < thisLen; i++)
02121                 _sum_row[i] = _row[i];
02122         }
02123         else
02124         {
02125             for (int i = 0; i < thisLen; i++)
02126             {
02127                 _sum_row[i] += _row[i];
02128             }
02129         }
02130     } EndForPencil;
02131 
02132     T _sum = _sum_row[0];
02133     for (int i = 1; i < _sum_len; i++)
02134         _sum += _sum_row[i];
02135 
02136     delete [] _sum_row;
02137 
02138     return _sum;
02139 }
02140 
02141 template <class T>
02142 BaseFab<T>&
02143 BaseFab<T>::negate ()
02144 {
02145     ForAllThis(T)
02146     {
02147         thisR = - thisR;
02148     } EndFor
02149     return *this;
02150 }
02151 
02152 template <class T>
02153 BaseFab<T>&
02154 BaseFab<T>::negate (int comp,
02155                     int numcomp)
02156 {
02157     ForAllThisNN(T,comp,numcomp)
02158     {
02159         thisR = - thisR;
02160     } EndFor
02161     return *this;
02162 }
02163 
02164 template <class T>
02165 BaseFab<T>&
02166 BaseFab<T>::negate (const Box& b,
02167                     int        comp,
02168                     int        numcomp)
02169 {
02170     ForAllThisBNN(T,b,comp,numcomp)
02171     {
02172         thisR = - thisR;
02173     } EndFor
02174     return *this;
02175 }
02176 
02177 template <class T>
02178 BaseFab<T>&
02179 BaseFab<T>::invert (T r)
02180 {
02181     ForAllThis(T)
02182     {
02183         thisR = r/thisR;
02184     } EndFor
02185     return *this;
02186 }
02187 
02188 template <class T>
02189 BaseFab<T>&
02190 BaseFab<T>::invert (T   r,
02191                     int comp,
02192                     int numcomp)
02193 {
02194     ForAllThisNN(T,comp,numcomp)
02195     {
02196         thisR = r/thisR;
02197     } EndFor
02198     return *this;
02199 }
02200 
02201 template <class T>
02202 BaseFab<T>&
02203 BaseFab<T>::invert (T          r,
02204                     const Box& b,
02205                     int        comp,
02206                     int        numcomp)
02207 {
02208     ForAllThisBNN(T,b,comp,numcomp)
02209     {
02210         thisR = r/thisR;
02211     } EndFor
02212     return *this;
02213 }
02214 
02215 template <class T>
02216 BaseFab<T>&
02217 BaseFab<T>::operator+= (T r)
02218 {
02219     ForAllThis(T)
02220     {
02221         thisR += r;
02222     } EndFor
02223     return *this;
02224 }
02225 
02226 template <class T>
02227 BaseFab<T>&
02228 BaseFab<T>::plus (T   r,
02229                   int comp,
02230                   int numcomp)
02231 {
02232     ForAllThisNN(T,comp,numcomp)
02233     {
02234         thisR += r;
02235     } EndFor
02236     return *this;
02237 }
02238 
02239 template <class T>
02240 BaseFab<T>&
02241 BaseFab<T>::plus (T          r,
02242                   const Box& b,
02243                   int        comp,
02244                   int        numcomp)
02245 {
02246     ForAllThisBNN(T,b,comp,numcomp)
02247     {
02248         thisR += r;
02249     } EndFor
02250     return *this;
02251 }
02252 
02253 template <class T>
02254 BaseFab<T>&
02255 BaseFab<T>::operator+= (const BaseFab<T>& x)
02256 {
02257     ForAllThisXC(T,x)
02258     {
02259         thisR += xR;
02260     } EndForTX
02261     return *this;
02262 }
02263 
02264 template <class T>
02265 BaseFab<T>&
02266 BaseFab<T>::plus (const BaseFab<T>& src,
02267                   int                srccomp,
02268                   int                destcomp,
02269                   int                numcomp)
02270 {
02271     ForAllThisBNNXC(T,domain,destcomp,numcomp,src,srccomp)
02272     {
02273         thisR += srcR;
02274     } EndForTX
02275     return *this;
02276 }
02277 
02278 template <class T>
02279 BaseFab<T>&
02280 BaseFab<T>::plus (const BaseFab<T>& src,
02281                   const Box&         subbox,
02282                   int                srccomp,
02283                   int                destcomp,
02284                   int                numcomp)
02285 {
02286     ForAllThisBNNXC(T,subbox,destcomp,numcomp,src,srccomp)
02287     {
02288         thisR += srcR;
02289     } EndForTX
02290     return *this;
02291 }
02292 
02293 template <class T>
02294 BaseFab<T>&
02295 BaseFab<T>::plus (const BaseFab<T>& src,
02296                   const Box&         srcbox,
02297                   const Box&         destbox,
02298                   int                srccomp,
02299                   int                destcomp,
02300                   int                numcomp)
02301 {
02302     ForAllThisBNNXCBN(T,destbox,destcomp,numcomp,src,srcbox,srccomp)
02303     {
02304         thisR += srcR;
02305     } EndForTX
02306     return *this;
02307 }
02308 
02309 template <class T>
02310 BaseFab<T>&
02311 BaseFab<T>::operator-= (const BaseFab<T>& x)
02312 {
02313     ForAllThisXC(T,x)
02314     {
02315         thisR -= xR;
02316     } EndForTX
02317     return *this;
02318 }
02319 
02320 template <class T>
02321 BaseFab<T>&
02322 BaseFab<T>::minus (const BaseFab<T>& src,
02323                    int                srccomp,
02324                    int                destcomp,
02325                    int                numcomp)
02326 {
02327     ForAllThisBNNXC(T,domain,destcomp,numcomp,src,srccomp)
02328     {
02329         thisR -= srcR;
02330     } EndForTX
02331     return *this;
02332 }
02333 
02334 template <class T>
02335 BaseFab<T>&
02336 BaseFab<T>::minus (const BaseFab<T>& src,
02337                    const Box&         subbox,
02338                    int                srccomp,
02339                    int                destcomp,
02340                    int                numcomp)
02341 {
02342     ForAllThisBNNXC(T,subbox,destcomp,numcomp,src,srccomp)
02343     {
02344         thisR -= srcR;
02345     } EndForTX
02346     return *this;
02347 }
02348 
02349 template <class T>
02350 BaseFab<T>&
02351 BaseFab<T>::minus (const BaseFab<T>& src,
02352                    const Box&         srcbox,
02353                    const Box&         destbox,
02354                    int                srccomp,
02355                    int                destcomp,
02356                    int                numcomp)
02357 {
02358     ForAllThisBNNXCBN(T,destbox,destcomp,numcomp,src,srcbox,srccomp)
02359     {
02360         thisR -= srcR;
02361     } EndForTX
02362     return *this;
02363 }
02364 
02365 template <class T>
02366 BaseFab<T>&
02367 BaseFab<T>::operator*= (T r)
02368 {
02369     ForAllThis(T)
02370     {
02371         thisR *= r;
02372     } EndFor
02373     return *this;
02374 }
02375 
02376 template <class T>
02377 BaseFab<T>&
02378 BaseFab<T>::mult (T   r,
02379                   int comp,
02380                   int numcomp)
02381 {
02382     ForAllThisNN(T,comp,numcomp)
02383     {
02384         thisR *= r;
02385     } EndFor
02386     return *this;
02387 }
02388 
02389 template <class T>
02390 BaseFab<T>&
02391 BaseFab<T>::mult (T          r,
02392                   const Box& b,
02393                   int        comp,
02394                   int        numcomp)
02395 {
02396     ForAllThisBNN(T,b,comp,numcomp)
02397     {
02398         thisR *= r;
02399     } EndFor
02400     return *this;
02401 }
02402 
02403 template <class T>
02404 BaseFab<T>&
02405 BaseFab<T>::operator*= (const BaseFab<T> &x)
02406 {
02407     ForAllThisXC(T,x)
02408     {
02409         thisR *= xR;
02410     } EndForTX
02411     return *this;
02412 }
02413 
02414 template <class T>
02415 BaseFab<T>&
02416 BaseFab<T>::mult (const BaseFab<T>& src,
02417                   int                srccomp,
02418                   int                destcomp,
02419                   int                numcomp)
02420 {
02421     ForAllThisBNNXC(T,domain,destcomp,numcomp,src,srccomp)
02422     {
02423         thisR *= srcR;
02424     } EndForTX
02425     return *this;
02426 }
02427 
02428 template <class T>
02429 BaseFab<T>&
02430 BaseFab<T>::mult (const BaseFab<T>& src,
02431                   const Box&         subbox,
02432                   int                srccomp,
02433                   int                destcomp,
02434                   int                numcomp)
02435 {
02436     ForAllThisBNNXC(T,subbox,destcomp,numcomp,src,srccomp)
02437     {
02438         thisR *= srcR;
02439     } EndForTX
02440     return *this;
02441 }
02442 
02443 template <class T>
02444 BaseFab<T>&
02445 BaseFab<T>::mult (const BaseFab<T>& src,
02446                   const Box&         srcbox,
02447                   const Box&         destbox,
02448                   int                srccomp,
02449                   int                destcomp,
02450                   int                numcomp)
02451 {
02452     ForAllThisBNNXCBN(T,destbox,destcomp,numcomp,src,srcbox,srccomp)
02453     {
02454         thisR *= srcR;
02455     } EndForTX
02456     return *this;
02457 }
02458 
02459 template <class T>
02460 BaseFab<T>&
02461 BaseFab<T>::operator/= (T r)
02462 {
02463     ForAllThis(T)
02464     {
02465         thisR /= r;
02466     } EndFor
02467     return *this;
02468 }
02469 
02470 template <class T>
02471 BaseFab<T>&
02472 BaseFab<T>::divide (T   r,
02473                     int comp,
02474                     int numcomp)
02475 {
02476     ForAllThisNN(T,comp,numcomp)
02477     {
02478         thisR /= r;
02479     } EndFor
02480     return *this;
02481 }
02482 
02483 template <class T>
02484 BaseFab<T>&
02485 BaseFab<T>::divide (T          r,
02486                     const Box& b,
02487                     int        comp,
02488                     int        numcomp)
02489 {
02490     ForAllThisBNN(T,b,comp,numcomp)
02491     {
02492         thisR /= r;
02493     } EndFor
02494     return *this;
02495 }
02496 
02497 template <class T>
02498 BaseFab<T>&
02499 BaseFab<T>::operator/= (const BaseFab<T> &x)
02500 {
02501     ForAllThisXC(T,x)
02502     {
02503         thisR /= xR;
02504     } EndForTX
02505     return *this;
02506 }
02507 
02508 template <class T>
02509 BaseFab<T>&
02510 BaseFab<T>::divide (const BaseFab<T>& src,
02511                     int                srccomp,
02512                     int                destcomp,
02513                     int                numcomp)
02514 {
02515     ForAllThisBNNXC(T,domain,destcomp,numcomp,src,srccomp)
02516     {
02517         thisR /= srcR;
02518     } EndForTX
02519     return *this;
02520 }
02521 
02522 template <class T>
02523 BaseFab<T>&
02524 BaseFab<T>::divide (const BaseFab<T>& src,
02525                     const Box&         subbox,
02526                     int                srccomp,
02527                     int                destcomp,
02528                     int                numcomp)
02529 {
02530     ForAllThisBNNXC(T,subbox,destcomp,numcomp,src,srccomp)
02531     {
02532         thisR /= srcR;
02533     } EndForTX
02534     return *this;
02535 }
02536 
02537 template <class T>
02538 BaseFab<T>&
02539 BaseFab<T>::divide (const BaseFab<T>& src,
02540                     const Box&         srcbox,
02541                     const Box&         destbox,
02542                     int                srccomp,
02543                     int                destcomp,
02544                     int                numcomp)
02545 {
02546     ForAllThisBNNXCBN(T,destbox,destcomp,numcomp,src,srcbox,srccomp)
02547     {
02548         thisR /= srcR;
02549     } EndForTX
02550     return *this;
02551 }
02552 
02553 //
02554 // Linear Interpolation / Extrapolation
02555 // Result is (t2-t)/(t2-t1)*f1 + (t-t1)/(t2-t1)*f2
02556 // Data is taken from b1 region of f1, b2 region of f2
02557 // and stored in b region of this FAB.
02558 // Boxes b, b1 and b2 must be the same size.
02559 // Data is taken from component comp1 of f1, comp2 of f2,
02560 // and stored in component comp of this FAB.
02561 // This fab is returned as a reference for chaining.
02562 //
02563 
02564 template <class T>
02565 BaseFab<T>&
02566 BaseFab<T>::linInterp (const BaseFab<T>& f1,
02567                        const Box&         b1,
02568                        int                comp1,
02569                        const BaseFab<T>& f2,
02570                        const Box&         b2,
02571                        int                comp2,
02572                        Real               t1,
02573                        Real               t2,
02574                        Real               t,
02575                        const Box&         b,
02576                        int                comp,
02577                        int                numcomp)
02578 {
02579     Real alpha = (t2-t)/(t2-t1);
02580     Real beta = (t-t1)/(t2-t1);
02581     ForAllThisBNNXCBNYCBN(T,b,comp,numcomp,f1,b1,comp1,f2,b2,comp2)
02582     {
02583         thisR = (T) (alpha*Real(f1R) + beta*Real(f2R));
02584     } EndForTX
02585     return *this;
02586 }
02587 
02588 //
02589 // Linear combination, Result is alpha*f1 + beta*f2
02590 // Data is taken from b1 region of f1, b2 region of f2
02591 // and stored in b region of this FAB.
02592 // Boxes b, b1 and b2 must be the same size.
02593 // Data is taken from component comp1 of f1, comp2 of f2,
02594 // and stored in component comp of this FAB.
02595 // This fab is returned as a reference for chaining.
02596 //
02597 
02598 template <class T>
02599 BaseFab<T>&
02600 BaseFab<T>::linComb (const BaseFab<T>&  f1,
02601                      const Box&         b1,
02602                      int                comp1,
02603                      const BaseFab<T>&  f2,
02604                      const Box&         b2,
02605                      int                comp2,
02606                      Real               alpha,
02607                      Real               beta,
02608                      const Box&         b,
02609                      int                comp,
02610                      int                numcomp)
02611 {
02612     ForAllThisBNNXCBNYCBN(T,b,comp,numcomp,f1,b1,comp1,f2,b2,comp2)
02613     {
02614         thisR = (T) (alpha*Real(f1R) + beta*Real(f2R));
02615     } EndForTX
02616     return *this;
02617 }
02618 
02619 namespace BoxLib
02620 {
02621     class BF_init
02622     {
02623     public:
02624         BF_init ();
02625         ~BF_init ();
02626     private:
02627         static int m_cnt;
02628     };
02629 }
02630 
02631 static BoxLib::BF_init file_scope_BF_init_object;
02632 
02633 
02634 #endif /*BL_BASEFAB_H*/