amroc/amr/SpecialOutput/AMRAccumulation.h

Go to the documentation of this file.

// -*- C++ -*-

#ifndef AMROC_AMR_ACCUMULATION_H
#define AMROC_AMR_ACCUMULATION_H

#include <string>
#include <cstdio>

template <class DataType, int dim>
class AMRAccumulation : public controlable {
public:
  typedef GridFunction<DataType,dim> grid_fct_type;   
  typedef GridData<DataType,dim> grid_data_type;  

  AMRAccumulation() : _Ghosts(0), _Dim(dim), _AccLevel(0) {
    _Hierarchy = (GridHierarchy*) 0;
    char Name[16];
    std::sprintf(Name,"acc");
    _AccName = Name;
    _Accfunc = (grid_fct_type *) 0;
  }  

  virtual ~AMRAccumulation() {
    if (_Accfunc) delete _Accfunc;
  }

  virtual inline void GridInitializationOp(grid_data_type& accgd) 
  { accgd.equals(0.0); }
  virtual inline void GridAccumulationOp(grid_data_type& accgd, grid_data_type& gd) 
  { accgd.maximum(gd); }
  virtual inline void PointAccumulationOp(DataType& AccPoint, const DataType& value)
  { if (value > AccPoint) AccPoint = value; }


  virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
    LocCtrl = Ctrl.getSubDevice(prefix+"Accumulation");
    RegisterAt(LocCtrl,"GridLevel",_AccLevel);
    RegisterAt(LocCtrl,"OutputName",_AccName);
  }
  virtual void register_at(ControlDevice& Ctrl) { register_at(Ctrl, ""); }

  virtual void finish() { 
    if (_Accfunc) delete _Accfunc;
    _Accfunc = (grid_fct_type *) 0;    
  }

  virtual void SetupData(GridHierarchy* gh, const int& ghosts) {
    _Hierarchy = gh;
    _Ghosts = ghosts;

    int t_sten = 0;   
    int s_sten = 0;      
    if (_AccLevel > TotalLevels(GH())-1) _AccLevel = TotalLevels(GH())-1;
    if (_AccLevel < 0) _AccLevel = 0;
    _Accfunc = new grid_fct_type(AccName(), t_sten, s_sten, GH(),
                                 DAGHCellCentered, -_AccLevel, 1, DAGH_All, 
                                 DAGHNoComm, DAGHNoBoundary, DAGHNoAdaptBoundary,
                                 DAGHNoExternalGhost);
    SetProlongFlag(AccFunc(), DAGHFalse);   
    SetCheckpointFlag(AccFunc(), DAGHTrue);
  }

  virtual void Initialize() {   
    forall (AccFunc(),0,0,c)      
      GridInitializationOp(AccFunc()(0,0,c));
    end_forall
  }
 
  virtual void Accumulate(grid_fct_type& u, const int& Time, 
                          const int& Level, const double& t) {

    int ATime = CurrentTime(GH(),0);
    Coords offset(dim,0);

    forall (u,Time,Level,c)
      forall (AccFunc(),ATime,0,mc)
        AccumulateGrid(u(Time,Level,c),u.interiorbbox(Time,Level,c),
                       AccFunc()(ATime,0,mc),AccFunc().interiorbbox(ATime,0,mc),offset);
      end_forall
    end_forall
  }    

  void AccumulateGrid(grid_data_type& gd, const BBox frombb, 
                      grid_data_type& accgd, BBox tobb, 
                      Coords offset) {
    Coords toss(tobb.stepsize());
    bool sum;
    if (tobb.stepsize(0)<=frombb.stepsize(0)) {
      sum = false;
      tobb *= refine(frombb,frombb.stepsize()/toss); 
    }
    else {
      sum = true;
      tobb *= coarsen(frombb,toss/frombb.stepsize());     
    }
    if (tobb.empty()) return;

    if (dim==1) {
      BeginFastIndex1(u, gd.bbox(), gd.data(), DataType);
      BeginFastIndex1(acc, accgd.bbox(), accgd.data(), DataType);
      for_1 (i, tobb, toss)
        DataType value = 0.0;
        if (!sum)
          value = FastIndex1(u,i-i%frombb.stepsize(0));
        else {
          BBox AddOver(1,i,i,tobb.stepsize(0));
          AddOver.refine(toss/frombb.stepsize());
          DataType count = 0.0;
          for_1 (l, AddOver, AddOver.stepsize())
            if (frombb.inside(l)) {
              value += FastIndex1(u,l);
              count += 1.0;
            }
          end_for
          if (count > 0.0) value /= count;
        }
        if (accgd.bbox().inside(i+offset(0))) 
          PointAccumulationOp(FastIndex1(acc,i+offset(0)),value);
      end_for
      EndFastIndex1(acc);      
      EndFastIndex1(u);
    }
    else if (dim==2) {
      BeginFastIndex2(u, gd.bbox(), gd.data(), DataType);
      BeginFastIndex2(acc, accgd.bbox(), accgd.data(), DataType);
      for_2 (i, j, tobb, toss)
        DataType value = 0.0;
        if (!sum)
          value = FastIndex2(u,i-i%frombb.stepsize(0),j-j%frombb.stepsize(1));
        else {
          BBox AddOver(2,i,j,i,j,tobb.stepsize(0),tobb.stepsize(1));
          AddOver.refine(toss/frombb.stepsize());
          DataType count = 0.0;
          for_2 (l, m, AddOver, AddOver.stepsize())
            if (frombb.inside(l,m)) {
              value += FastIndex2(u,l,m);
              count += 1.0;
            }
          end_for
          if (count > 0.0) value /= count;
        }
        if (accgd.bbox().inside(i+offset(0),j+offset(1))) 
          PointAccumulationOp(FastIndex2(acc,i+offset(0),j+offset(1)),value);
      end_for
      EndFastIndex2(acc);      
      EndFastIndex2(u);
    }
    else if (dim==3) {
      BeginFastIndex3(u, gd.bbox(), gd.data(), DataType);
      BeginFastIndex3(acc, accgd.bbox(), accgd.data(), DataType);
      for_3 (i, j, k, tobb, toss)
        DataType value = 0.0;
        if (!sum)
          value = FastIndex3(u,i-i%frombb.stepsize(0),j-j%frombb.stepsize(1),
                             k-k%frombb.stepsize(2));
        else {
          BBox AddOver(3,i,j,k,i,j,k,tobb.stepsize(0),tobb.stepsize(1),
                       tobb.stepsize(2));
          AddOver.refine(toss/frombb.stepsize());
          DataType count = 0.0;
          for_3 (l, m, n, AddOver, AddOver.stepsize())
            if (frombb.inside(l,m,n)) {
              value += FastIndex3(u,l,m,n);
              count += 1.0;
            }
          end_for
          if (count > 0.0) value /= count;
        }
        if (accgd.bbox().inside(i+offset(0),j+offset(1),k+offset(2))) 
          PointAccumulationOp(FastIndex3(acc,i+offset(0),j+offset(1),k+offset(2)),value);
      end_for
      EndFastIndex3(acc);      
      EndFastIndex3(u);
    }
  }    

//   void Grad(grid_fct_type& u, const int Time, const int Level) {
//     int TStep = TimeStep(u,Level);
//     forall (u,Time,Level,c)
//       BBox bb(u(Time,Level,c).bbox());
//       BBox bbi(u.interiorbbox(Time,Level,c));
//       Coords ss(u(Time,Level,c).bbox().stepsize());
//       DCoords dx = GH().worldStep(ss);
//       Coords lb(bb.lower());
      
//       if (dim == 1) {
//      BeginFastIndex1(U, bb, u(Time,Level,c).data(), DataType)
//      BeginFastIndex1(Un, bb, u(Time+TStep,Level,c).data(), DataType)
//      for_1 (i, bbi, ss)
//        Fastindex1(Un,i) = std::fabs(FastIndex1(U,i+ss(0))-FastIndex1(U,i-ss(0)))/(2.0*dx[0]);
//      end_for
//      EndFastIndex1(U)   
//      EndFastIndex1(Un)    
//       }
//       else if if (dim == 2) {
//      BeginFastIndex2(U, bb, u(Time,Level,c).data(), DataType);
//      BeginFastIndex2(Un, bb, u(Time+TStep,Level,c).data(), DataType);
//      double xgrad, ygrad;
//      for_2 (i, j, bbi, ss)
//        xgrad = (FastIndex2(U,i+ss(0),j)-FastIndex2(U,i-ss(0),j))/(2.0*dx[0]);
//        ygrad = (FastIndex2(U,i,j+ss(1))-FastIndex2(U,i,j-ss(1)))/(2.0*dx[1]);
//        FastIndex2(Un,i,j) = std::sqrt(xgrad*xgrad+ygrad*ygrad);
//      end_for
//      EndFastIndex2(U)     
//         EndFastIndex2(Un)   
//       }  
//       else if if (dim == 3) {
//      BeginFastIndex3(U, bb, u(Time,Level,c).data(), DataType)
//      BeginFastIndex3(Un, bb, u(Time+TStep,Level,c).data(), DataType)
//      double xgrad, ygrad, zgrad;
//      for_3 (i, j, k, bbi, ss)
//        xgrad = (FastIndex3(U,i+ss(0),j,k)-FastIndex3(U,i-ss(0),j,k))/(2.0*dx[0]);
//           ygrad = (FastIndex3(U,i,j+ss(1),k)-FastIndex3(U,i,j-ss(1),k))/(2.0*dx[1]);
//        zgrad = (FastIndex3(U,i,j,k+ss(2))-FastIndex3(U,i,j,k-ss(2)))/(2.0*dx[2]);
//        FastIndex3(Un,i,j,k) = std::sqrt(xgrad*xgrad+ygrad*ygrad+zgrad*zgrad);
//      end_for
//      EndFastIndex3(U)      
//      EndFastIndex3(Un)  
//       } 

//     end_forall
//   }
    
  GridHierarchy& GH() { return *_Hierarchy; }
  const GridHierarchy& GH() const { return *_Hierarchy; }

  const int& NGhosts() const { return _Ghosts; }
  const int& Dim() const { return _Dim; }

  inline grid_fct_type& AccFunc() { return *_Accfunc;}
  inline const char* AccName() { return _AccName.c_str(); }
  inline const int& AccLevel() const { return _AccLevel; }
  inline int AccLevel() { return _AccLevel; }

protected:
  int _Ghosts;
  int _Dim;
  int _AccLevel;
  std::string _AccName;
  grid_fct_type *_Accfunc;
  GridHierarchy* _Hierarchy;
  ControlDevice LocCtrl;
}; 
 

#endif

---