amroc/amr/AMRFlagging.h

Go to the documentation of this file.

// -*- C++ -*-

// Copyright (C) 2002 Ralf Deiterding
// Brandenburgische Universitaet Cottbus

#ifndef AMROC_AMRFLAGGING_H
#define AMROC_AMRFLAGGING_H

#include "Flagging.h"

template <class VectorType, class FixupType, class FlagType, int dim> class AMRSolverBase;

#define GoodPoint (0)
#define BadPoint  (1)

template <class VectorType, class FixupType, class FlagType, int dim>
class AMRFlagging : public Flagging<VectorType,FlagType,dim> {
  typedef Flagging<VectorType,FlagType,dim> base;
  typedef typename VectorType::InternalDataType DataType;
public:
  typedef typename base::vec_grid_fct_type vec_grid_fct_type;  
  typedef typename base::vec_grid_data_type vec_grid_data_type;
  typedef typename base::flag_fct_type flag_fct_type;
  typedef typename base::criterion_type criterion_type;
  typedef GridFunction<DataType,dim> grid_fct_type;
  typedef GridData<DataType,dim> grid_data_type;
  typedef GridData<FlagType,dim> flag_data_type;
  typedef AMRSolverBase<VectorType,FixupType,FlagType,dim> solver_type;

  AMRFlagging(solver_type& solver) : base(), _Solver(solver), _EstimateError(false) {}

  virtual ~AMRFlagging() {}

  virtual void update() {
    base::update();
    int nc;
    for (nc=0; nc<base::NCrit(); nc++) 
      if (base::Criterion_(nc).ShadowCriterion() && base::Criterion_(nc).IsUsed() &&
          base::Criterion_(nc).UpdateShadow()) {
        _EstimateError = true;
        break;
      }
    if (!ErrorEstimation())
      for (nc=0; nc<base::NCrit(); nc++) 
        if (base::Criterion_(nc).ShadowCriterion())
          base::Criterion_(nc).SetIsUsed(false);
  }

  virtual void SetFlags(const int Time, const int Level, double t, double dt) {
    if (ErrorEstimation() && Time>0) {
      //---------------------------------------------------------------
      // Take a step on the U(). U()(Time,Level) has to remain 
      // unchanged in a fractional step method!
      //---------------------------------------------------------------
      END_WATCH_FLAGGING
        Solver_().IntegrateLevel(U(), Time, Level, t, dt, false, 0.0, 1); 
      START_WATCH
      int TStep = TimeStep(U(),Level);
      int ShTStep = TimeStep(Ush(),Level);
      forall (U(),Time,Level,c)      
        //------------------------------------------------------------------------
        // First update Ush() from U()
        //------------------------------------------------------------------------
        Ush()(Time+ShTStep,Level,c).equals(Ush()(Time,Level,c));
        Solver_().LevelTransfer_().Restrict(U()(Time+TStep,Level,c),Level,
                                            Ush()(Time+ShTStep,Level,c),Level,
                                            U().interiorbbox(Time+TStep,Level,c)); 
      end_forall
    }

    forall (base::Flags(),Time,Level,c)  
      base::Flags()(Time,Level,c) = GoodPoint;
    end_forall
    
    FlagType FlagValue = BadPoint;
    for (int nc=0; nc<base::NCrit(); nc++) {
      if (!base::Criterion_(nc).IsUsed()) continue;
      vec_grid_fct_type* uh = (!base::Criterion_(nc).ShadowCriterion() ? _u : _u_sh);
      grid_fct_type* wh = (!base::Criterion_(nc).ShadowCriterion() ? _work : _work_sh);
      if ((uh==(vec_grid_fct_type*) 0) || (wh==(grid_fct_type*)0 )) continue;
      for (int cnt=0; cnt<base::Criterion_(nc).Ncnt(); cnt++) {
        bool Used = base::Criterion_(nc).SetFlags(*uh, *wh, base::Flags(), cnt,
                                                  Time, Level, t, FlagValue);
        if (Used && base::Criterion_(nc).OutputFlags() && Solver_().LastOutputTime()==Time) {
          int TStep = TimeStep(Work(),Level);
          if (base::Criterion_(nc).ShadowCriterion()) {
            int ShTStep = TimeStep(Worksh(),Level);
            forall(Work(), Time, Level, c)
              equals_fromsh(Work()(Time+TStep,Level,c),Worksh()(Time+ShTStep,Level,c));
            end_forall
          }
          SwapTimeLevels(Work(),Level,Time,Time+TStep);
          char name[DAGHBktGFNameWidth+16];
          base::Criterion_(nc).OutputName(name,cnt);
          bool flushfile=(Level==0 ? true : false);
          Solver_().FileOutput_().WritePlain(Work(), name, Time, Level, t, BBox::_empty_bbox, false, flushfile);
          SwapTimeLevels(Work(),Level,Time,Time+TStep);
        }
        if (Used) FlagValue++;
      }
    }
  }

  void equals_fromsh(grid_data_type& w, grid_data_type& wsh) {
    BBox OpBox = w.bbox();
    BBox OpBoxsh = wsh.bbox();
    Coords& Os = OpBox.stepsize();
    Coords& Osh = OpBoxsh.stepsize();

    if (dim == 1) {
      BeginFastIndex1(w, w.bbox(), w.data(), DataType);        
      BeginFastIndex1(wsh, wsh.bbox(), wsh.data(), DataType); 
      for_1 (n, OpBox, Os)
        FastIndex1(w,n) = FastIndex1(wsh,n-n%Osh(0));
      end_for
      EndFastIndex1(w);
      EndFastIndex1(wsh);
    }
    else if (dim == 2) {
      BeginFastIndex2(w, w.bbox(), w.data(), DataType);        
      BeginFastIndex2(wsh, wsh.bbox(), wsh.data(), DataType); 
      for_2 (n, m, OpBox, Os)
        FastIndex2(w,n,m) = FastIndex2(wsh,n-n%Osh(0), m-m%Osh(1));
      end_for
      EndFastIndex2(w);
      EndFastIndex2(wsh);
    }
    else if (dim == 3) {
      BeginFastIndex3(w, w.bbox(), w.data(), DataType);        
      BeginFastIndex3(wsh, wsh.bbox(), wsh.data(), DataType); 
      for_3 (n, m, l, OpBox, Os)
        FastIndex3(w,n,m,l) = FastIndex3(wsh,n-n%Osh(0), m-m%Osh(1), l-l%Osh(2));
      end_for
      EndFastIndex3(w);
      EndFastIndex3(wsh);
    }
  }

  inline void SetGridFunctions(vec_grid_fct_type* u, vec_grid_fct_type* ush, 
                               grid_fct_type* work, grid_fct_type* worksh) {
    _u = u; _u_sh = ush;
    _work = work; _work_sh = worksh; 
  }
  inline vec_grid_fct_type& U() { return *_u; }
  inline const vec_grid_fct_type& U() const { return *_u; }
  inline vec_grid_fct_type& Ush() { return *_u_sh; }
  inline const vec_grid_fct_type& Ush() const { return *_u_sh; }
  inline grid_fct_type& Work() { return *_work; }
  inline const grid_fct_type& Work() const { return *_work; }
  inline grid_fct_type& Worksh() { return *_work_sh; }
  inline const grid_fct_type& Worksh() const { return *_work_sh; }
  inline const bool& ErrorEstimation() const { return _EstimateError;}
  void SetErrorEstimation(const bool est) { _EstimateError = est; }
  inline solver_type& Solver_() { return _Solver; }
  inline const solver_type& Solver_() const { return _Solver; }

protected:
  solver_type& _Solver;
  vec_grid_fct_type *_u, *_u_sh;
  grid_fct_type *_work, *_work_sh;
  bool _EstimateError;
};


#endif

---