basefabimplem.h

自适应网格划分通用程序包

// This software is copyright (C) by the Lawrence Berkeley
// National Laboratory.  Permission is granted to reproduce
// this software for non-commercial purposes provided that
// this notice is left intact.
// 
// It is acknowledged that the U.S. Government has rights to
// this software under Contract DE-AC03-765F00098 between
// the U.S.  Department of Energy and the University of
// California.
//
// This software is provided as a professional and academic
// contribution for joint exchange. Thus it is experimental,
// is provided ``as is'', with no warranties of any kind
// whatsoever, no support, no promise of updates, or printed
// documentation. By using this software, you acknowledge
// that the Lawrence Berkeley National Laboratory and
// Regents of the University of California shall have no
// liability with respect to the infringement of other
// copyrights by any part of this software.
//

#ifndef CH_BASEFABIMPLEM_H
#define CH_BASEFABIMPLEM_H
#include "BaseFabMacros.H"

#include "MayDay.H"
//
// BaseFab<Real> spaecializations
//

template < >
void
BaseFab<Real>::define();

template < >
void 
BaseFab<Real>::undefine();

template < >
void
BaseFab<Real>::setVal (Real val);


template < >
void
BaseFab<int>::define();

template < >
void 
BaseFab<int>::undefine();

//
// Implementation.=====================================
//

template <class T>
Arena* BaseFab<T>::s_Arena = NULL;

template <class T>
inline
int
BaseFab<T>::nComp () const
{
  return nvar;
}

template <class T>
inline
const Box&
BaseFab<T>::box () const
{
  return domain;
}

template <class T>
inline
const int*
BaseFab<T>::size () const
{
  return domain.size().getVect();
}

template <class T>
inline
const IntVect&
BaseFab<T>::smallEnd () const
{
  return domain.smallEnd();
}

template <class T>
inline
const IntVect&
BaseFab<T>::bigEnd () const
{
  return domain.bigEnd();
}

template <class T>
inline
const int*
BaseFab<T>::loVect () const
{
  return domain.loVect();
}

template <class T>
inline
const int*
BaseFab<T>::hiVect () const
{
  return domain.hiVect();
}

template <class T>
inline
bool
BaseFab<T>::contains (const BaseFab<T>& fab) const
{
  return box().contains(fab.box()) && nvar <= fab.nvar;
}

template <class T>
inline
bool
BaseFab<T>::contains (const Box& bx) const
{
  return box().contains(bx);
}

template <class T>
inline
T*
BaseFab<T>::dataPtr (int n)
{
  assert(!(dptr == 0));
  assert((n >= 0) && (n < nvar));
  return &dptr[n*numpts];
}

template <class T>
inline
const int*
BaseFab<T>::nCompPtr () const
{
  assert(!(dptr == 0));
  return &nvar;
}

template <class T>
inline
const T*
BaseFab<T>::dataPtr (int n) const
{
  assert(!(dptr == 0));
  assert((n >= 0) && (n < nvar));
  return &dptr[n*numpts];
}

template <class T>
inline
T&
BaseFab<T>::operator() (const IntVect& p,
                        int            n)
{
  assert(n >= 0);
  assert(n < nvar);
  assert(!(dptr == 0));
  assert(domain.contains(p));

  return dptr[domain.index(p)+n*numpts];
}

template <class T>
inline
T&
BaseFab<T>::operator() (const IntVect& p)
{
  assert(!(dptr == 0));
  assert(domain.contains(p));

  return dptr[domain.index(p)];
}

template <class T>
inline
const T&
BaseFab<T>::operator() (const IntVect& p,
                        int            n) const
{
  assert(n >= 0);
  assert(n < nvar);
  assert(!(dptr == 0));
  assert(domain.contains(p));

  return dptr[domain.index(p)+n*numpts];
}

template <class T>
inline
const T&
BaseFab<T>::operator() (const IntVect& p) const
{
  assert(!(dptr == 0));
  assert(domain.contains(p));

  return dptr[domain.index(p)];
}

template <class T>
inline
void
BaseFab<T>::getVal  (T*             data,
                     const IntVect& pos,
                     int            n,
                     int            numcomp) const
{
  const int loc      = domain.index(pos);
  const long size    = domain.numPts();

  assert(!(dptr == 0));
  assert(n >= 0 && n + numcomp <= nvar);

  for (int k = 0; k < numcomp; k++)
    data[k] = dptr[loc+(n+k)*size];
}

template <class T>
inline
void
BaseFab<T>::getVal (T*             data,
                    const IntVect& pos) const
{
  getVal(data,pos,0,nvar);
}

template <class T>
inline
BaseFab<T>&
BaseFab<T>::shift (const IntVect& v)
{
  domain += v;
  return *this;
}

template <class T>
inline
BaseFab<T>&
BaseFab<T>::shift (int idir,
                   int n_cell)
{
  domain.shift(idir,n_cell);
  return *this;
}

template <class T>
inline
BaseFab<T> &
BaseFab<T>::shiftHalf (const IntVect& v)
{
  domain.shiftHalf(v);
  return *this;
}

template <class T>
inline
BaseFab<T> &
BaseFab<T>::shiftHalf (int idir,
                       int n_cell)
{
  domain.shiftHalf(idir,n_cell);
  return *this;
}

template <class T>
inline
void
BaseFab<T>::setVal (T val)
{
  performSetVal(val,box(), 0, nvar);
}

template <class T>
inline
void
BaseFab<T>::setVal (T          x,
                    const Box& bx,
                    int        n)
{
  performSetVal(x,bx,n,1);
}

template <class T>
inline
void
BaseFab<T>::setVal (T   x,
                    int n)
{
  performSetVal(x,domain,n,1);
}

template <class T>
inline
void
BaseFab<T>::setVal (T          x,
                    const Box& b,
                    int        ns,
                    int        num)
{
  performSetVal(x,b,ns,num);
}

template <class T>
inline
BaseFab<T>&
BaseFab<T>::copy (const BaseFab<T>& src,
                  const Box&        srcbox,
                  int               srccomp,
                  const Box&        destbox,
                  int               destcomp,
                  int               numcomp)
{
  assert(srcbox.sameSize(destbox));
  assert(src.box().contains(srcbox));
  assert(domain.contains(destbox));
  assert(srccomp >= 0 && srccomp+numcomp <= src.nComp());
  assert(destcomp >= 0 && destcomp+numcomp <= nvar);
  performCopy(src,srcbox,srccomp,destbox,destcomp,numcomp);
  return *this;
}

template <class T>
inline
BaseFab<T>&
BaseFab<T>::copy (const BaseFab<T>& src)
{
  assert(nvar <= src.nvar);
  assert(domain.sameType(src.domain));
  Box overlap(domain);
  overlap &= src.domain;
  if (!overlap.isEmpty())
    performCopy(src,overlap,0,overlap,0,nvar);
  return *this;
}

template <class T>
inline
BaseFab<T>&
BaseFab<T>::copy (const BaseFab<T>& src,
                  const Box&        destbox)
{
  assert(nvar <= src.nvar);
  assert(domain.contains(destbox));
  Box overlap(destbox);
  overlap &= src.domain;
  if (!overlap.isEmpty())
    performCopy(src,overlap,0,overlap,0,nvar);
  return *this;
}

template <class T>
inline
BaseFab<T>&
BaseFab<T>::copy (const BaseFab<T>& src,
                  int               srccomp,
                  int               destcomp,
                  int               numcomp)
{
  assert(srccomp >= 0 && srccomp + numcomp <= src.nvar);
  assert(destcomp >= 0 && destcomp + numcomp <= nvar);
  Box overlap(domain);
  overlap &= src.domain;
  if (!overlap.isEmpty())
    performCopy(src,overlap,srccomp,overlap,destcomp,numcomp);
  return *this;
}

template <class T>
inline
void
BaseFab<T>::define (const Interval& a_comps, BaseFab<T>& a_original)
{
  undefine();
  domain   = a_original.domain;
  numpts   = a_original.numpts;
  truesize = a_original.numpts*a_comps.size();
  nvar     = a_comps.size();
  dptr     = a_original.dataPtr(a_comps.begin());
  aliased  = true;

  //resize(a_original.domain, a_comps.size(), a_original.dataPtr(a_comps.begin()));
}


template <class T>
inline
void
BaseFab<T>::define ()
{
  assert(nvar > 0);
  assert(dptr == 0);
  assert(numpts > 0);
  assert(!aliased);
  //    assert(!(The_FAB_Arena == 0));// not a sufficient test !!!
#ifdef ENABLE_MEMORY_TRACKING  
  if(s_Arena == NULL) s_Arena = new BArena(name().c_str());
#else
  if(s_Arena == NULL) s_Arena = new BArena("");
#endif
  if(s_Arena == NULL)
    MayDay::Error("malloc in basefab failed");
  truesize = nvar*numpts;
  dptr     = static_cast<T*>(s_Arena->alloc(truesize*sizeof(T)));

#ifdef ENABLE_MEMORY_TRACKING  
  s_Arena->bytes+=truesize*sizeof(T)+sizeof(BaseFab<T>);
  if(s_Arena->bytes > s_Arena->peak) 
    s_Arena->peak = s_Arena->bytes;
#endif

  //
  // Now call T::T() on the raw memo'ry so we have valid Ts.
  //
  T* ptr = dptr;

  for (int i = 0; i < truesize; i++, ptr++)
    new (ptr) T;
}

template <class T>
inline
void
BaseFab<T>::undefine ()
{
  //    assert(!(The_FAB_Arena == 0));
  //
  // Call T::~T() on the to-be-destroyed memory.
  //
  if(aliased)
        {
          dptr = 0;
          return;
        }
  if(dptr == 0) return;
  T* ptr = dptr;

  for (int i = 0; i < truesize; i++, ptr++)
    ptr->~T();

  s_Arena->free(dptr);

#ifdef ENABLE_MEMORY_TRACKING  
  s_Arena->bytes-= truesize*sizeof(T)+sizeof(BaseFab<T>);
#endif

  dptr = 0;
}

template <class T>
inline
BaseFab<T>::BaseFab ()
  :
  domain(Box()),
  nvar(0),
  numpts(0),
  truesize(0),
  dptr(0),
  aliased(false)
{}

template <class T>
inline
BaseFab<T>::BaseFab (const Interval& a_comps, BaseFab<T>& a_original)
  :
  domain(a_original.domain),
  nvar(a_comps.size()),
  numpts(a_original.numpts),
  truesize(a_original.numpts*nvar),
  dptr(a_original.dataPtr(a_comps.begin())),
  aliased(true)
{}

template <class T>
inline
BaseFab<T>::BaseFab (const Box& bx,
                     int        n,
                                         T*         alias)
  :
  domain(bx),
  nvar(n),
  numpts(bx.numPts()),
  dptr(0),
  aliased(false)
{
  if(alias != NULL)
        {
          dptr = alias;
          aliased = true;
        }
  else
        define();
}


template <class T>
void
BaseFab<T>::resize (const Box& b,
                    int        n,
                                        T*         alias)
{
  //     nvar   = n;
  //     domain = b;
  //     numpts = domain.numPts();

  //     if (dptr == 0)
  //     {
  //         define();
  //     }
  //     else if (nvar*numpts > truesize)
  //     {
  //         undefine();

//         define();
//     }
  undefine();
  nvar = n;
  domain = b;
  numpts = domain.numPts();
  if(alias != NULL)
        {
          dptr = alias;
          aliased = true;
        }
  else
        {
          aliased = false;
          define();
        }
}

template <class T>
inline
BaseFab<T>::~BaseFab ()
{
  undefine();
}

template <class T>
inline
void
BaseFab<T>::clear ()
{
  undefine();

  domain = Box();
  nvar   = 0;
  numpts = 0;
}

template <class T>
inline
void
BaseFab<T>::performCopy (const BaseFab<T>& src,
                         const Box&        srcbox,
                         int               srccomp,
                         const Box&        destbox,
                         int               destcomp,
                         int               numcomp)
{
  assert(src.box().contains(srcbox));
  assert(box().contains(destbox));
  assert(destbox.sameSize(srcbox));
  assert(srccomp >= 0 && srccomp+numcomp <= src.nComp());
  assert(destcomp >= 0 && destcomp+numcomp <= nComp());

  ForAllThisBNNXCBN(T, destbox, destcomp, numcomp, src, srcbox, srccomp)
    {
      thisR = srcR;
    } EndForTX
  ;
}

template <class T>
inline
void
BaseFab<T>::performSetVal (T         val,
                           const Box& bx,
                           int        ns,
                           int        num)
{
  assert(domain.contains(bx));
  assert(ns >= 0 && ns + num <= nvar);

  if (bx == domain)
    {
      T* data = &dptr[ns*numpts];
      for (long i = 0, N = num*numpts; i < N; i++)
        {
          *data++ = val;
        }
    }
  else
    {
      ForAllThisBNN(T,bx,ns,num)
        {
          thisR = val;
        } EndFor
            }
}


template <class T>
inline
void
BaseFab<T>::copy(const Box& RegionFrom, 
                 const Interval& Cdest, 
                 const Box& RegionTo, 
                 const BaseFab<T>& src, 
                 const Interval& Csrc)
{
  if((this == &src) && (RegionFrom == RegionTo)) return;
  assert(Cdest.size() == Csrc.size()); 
  copy(src, RegionFrom, Csrc.begin(), RegionTo, 
       Cdest.begin(), Cdest.size());
}

template <class T>
inline
void
BaseFab<T>::linearOut(void* buf, const Box& R, const Interval& comps) const
{
  T* buffer = (T*)buf;
  ForAllThisBNN(T,R,comps.begin(), comps.size())
    {
      *buffer = thisR;
      ++buffer;
    } EndFor;
}

template <class T>
inline
void
BaseFab<T>::linearIn(void* buf, const Box& R, const Interval& comps)
{
  T* buffer = (T*)buf;
  ForAllThisBNN(T,R,comps.begin(), comps.size())
    {
      thisR = *buffer;
      ++buffer;
    } EndFor;

}

template <class T>
inline
int
BaseFab<T>::size(const Box& box, const Interval& comps) const
{
  return box.numPts()*sizeof(T)*comps.size();
}

template <class T>
inline
std::string
BaseFab<T>::name()
{
  std::string rtn = (typeid(T)).name();
  return rtn;
}

#endif /*CH_BASEFAB_H*/