LevelDataOps.H

Go to the documentation of this file.

#ifdef CH_LANG_CC
/*
*      _______              __
*     / ___/ /  ___  __ _  / /  ___
*    / /__/ _ \/ _ \/  V \/ _ \/ _ \
*    \___/_//_/\___/_/_/_/_.__/\___/
*    Please refer to Copyright.txt, in Chombo's root directory.
*/
#endif

// BVS, June 26, 2003

// We can assume that template class T has null construction.

#ifndef _LEVELDATAOPS_H_
#define _LEVELDATAOPS_H_

#ifdef CH_MPI
#include <string>
#include <sstream>
#endif
#include "LevelData.H"
#include "RefCountedPtr.H"
#include "SPMD.H"
#include "Copier.H"
#include "NamespaceHeader.H"


// default copy constructor and assign are fine.

template <class T>
class LevelDataOps
{
public:
  LevelDataOps()
    :m_levelFactory( new DefaultDataFactory<T>() )
  {
  }

  LevelDataOps(RefCountedPtr<DataFactory<T> > a_factoryPtr)
    : m_levelFactory(a_factoryPtr){;}
  virtual ~LevelDataOps(){;}
  virtual void define(const RefCountedPtr<DataFactory<T> >& a_factoryPtr)
    {m_levelFactory = a_factoryPtr;}
  virtual void define(DataFactory<T>* a_rawPointer)
    {m_levelFactory = RefCountedPtr<DataFactory<T> >(a_rawPointer);}
  virtual void create(    LevelData<T>& a_lhs, const LevelData<T>& a_rhs);
  virtual void assign(    LevelData<T>& a_lhs, const LevelData<T>& a_rhs) ;
  virtual Real dotProduct(const LevelData<T>& a_1, const LevelData<T>& a_2) ;
  virtual void incr( LevelData<T>& a_lhs, const LevelData<T>& a_x, Real a_scale) ;
  virtual void mult( LevelData<T>& a_lhs, const LevelData<T>& a_x);
  virtual void axby( LevelData<T>& a_lhs, const LevelData<T>& a_x,
                     const LevelData<T>& a_y, Real a_a, Real a_b) ;
  virtual void scale(LevelData<T>& a_lhs, const Real& a_scale) ;
  virtual void plus(LevelData<T>& a_lhs, const Real& a_inc) ;
  virtual void setToZero(LevelData<T>& a_lhs);
  virtual void copyToZero(LevelData<T>& a_lhs, const Copier& a_copier);


protected:
  RefCountedPtr<DataFactory<T> > m_levelFactory;
};

//*******************************************************
// LevelDataOps Implementation
//*******************************************************
// BVS, June 26, 2003

#define ITER(a) for(DataIterator dit = a.dataIterator(); dit.ok(); ++dit) { \
                   DataIndex d = dit();
#define ENDFOR(a) }

template <class T>
void LevelDataOps<T>:: create(LevelData<T>& a_lhs, const LevelData<T>& a_rhs)
{
  a_lhs.define(a_rhs, *m_levelFactory);
}

template <class T>
void LevelDataOps<T>:: assign(LevelData<T>& a_lhs, const LevelData<T>& a_rhs)
{
  Interval interv(0, a_rhs.nComp()-1);
  a_rhs.copyTo(interv, a_lhs, interv);
}

template <class T>
Real  LevelDataOps<T>::dotProduct(const LevelData<T>& a_1, const LevelData<T>& a_2)
{
  const DisjointBoxLayout& dbl = a_1.disjointBoxLayout();
  Real val = 0.0;
  ITER(a_1)
    val += a_1[d].dotProduct(a_2[d], dbl.get(d));
  ENDFOR(a_1);

#ifdef CH_MPI
  Real recv;
  int result = MPI_Allreduce(&val, &recv, 1, MPI_CH_REAL,
                             MPI_SUM, Chombo_MPI::comm);
  if( result != 0 ){
    std::ostringstream msg;
    msg << "LevelDataOps::dotProduct() called MPI_Allreduce() which returned error code " << result ;
    MayDay::Warning( msg.str().c_str() );
  }
  val = recv;
#endif
  return val;

}

template <class T>
void LevelDataOps<T>:: incr( LevelData<T>& a_lhs, const LevelData<T>& a_rhs, Real a_scale)
{
  //  CH_assert(a_lhs.disjointBoxLayout() == a_rhs.disjointBoxLayout());
  int numcomp = a_lhs.nComp();
  int  startcomp = 0;
  ITER(a_lhs)
    Box subbox = a_lhs.disjointBoxLayout()[d];
    a_lhs[d].plus(a_rhs[d],  subbox, subbox, a_scale, startcomp, startcomp, numcomp);
  ENDFOR(a_lhs);
}

template <class T>
void LevelDataOps<T>:: mult( LevelData<T>& a_lhs, const LevelData<T>& a_rhs)
{
  //  CH_assert(a_lhs.disjointBoxLayout() == a_rhs.disjointBoxLayout());
  ITER(a_lhs)
    a_lhs[d] *= a_rhs[d];
  ENDFOR(a_lhs);
}

template <class T>
void LevelDataOps<T>:: setToZero( LevelData<T>& a_lhs)
{
  ITER(a_lhs)
    a_lhs[d].setVal(0.0);
  ENDFOR(a_lhs);
}

template <class T>
void LevelDataOps<T>:: copyToZero( LevelData<T>& a_lhs, const Copier& a_copier)
{
  int nComp = a_lhs.nComp();
  for(CopyIterator it(a_copier, CopyIterator::LOCAL); it.ok(); ++it)
    {
      const MotionItem& item = it();
      a_lhs[item.toIndex].setVal(0, item.toRegion, 0, nComp);

    }
  for(CopyIterator it(a_copier, CopyIterator::TO); it.ok(); ++it)
    {
      const MotionItem& item = it();
      a_lhs[item.toIndex].setVal(0, item.toRegion, 0, nComp);
    }
}

template <class T>
void LevelDataOps<T>:: axby( LevelData<T>& a_lhs, const LevelData<T>& a_x,
                             const LevelData<T>& a_y, Real a, Real b)
{
  //  CH_assert(a_lhs.disjointBoxLayout() == a_x.disjointBoxLayout());
  ITER(a_lhs)
    T& data = a_lhs[d];
  data.copy(a_x[d]);
  data.mult(a);
  data.plus(a_y[d], b);
  ENDFOR(a_lhs);
}

template <class T>
void LevelDataOps<T>:: scale(LevelData<T>& a_lhs, const Real& a_scale)
{
  ITER(a_lhs)
    T& data = a_lhs[d];
  data.mult(a_scale);
  ENDFOR(a_lhs);
}

template <class T>
void LevelDataOps<T>:: plus(LevelData<T>& a_lhs, const Real& a_inc)
{
  ITER(a_lhs)
    T& data = a_lhs[d];
  data += a_inc;
  ENDFOR(a_lhs);
}

#undef ITER
#undef ENDFOR

#include "NamespaceFooter.H"
#endif

---