00001
00002
00003
00004
00005
00006 #ifndef AMROC_AMRFLAGGING_H
00007 #define AMROC_AMRFLAGGING_H
00008
00016 #include "Flagging.h"
00017
00018 template <class VectorType, class FixupType, class FlagType, int dim> class AMRSolverBase;
00019
00020 #define GoodPoint (0)
00021 #define BadPoint (1)
00022
00030 template <class VectorType, class FixupType, class FlagType, int dim>
00031 class AMRFlagging : public Flagging<VectorType,FlagType,dim> {
00032 typedef Flagging<VectorType,FlagType,dim> base;
00033 typedef typename VectorType::InternalDataType DataType;
00034 public:
00035 typedef typename base::vec_grid_fct_type vec_grid_fct_type;
00036 typedef typename base::vec_grid_data_type vec_grid_data_type;
00037 typedef typename base::flag_fct_type flag_fct_type;
00038 typedef typename base::criterion_type criterion_type;
00039 typedef GridFunction<DataType,dim> grid_fct_type;
00040 typedef GridData<DataType,dim> grid_data_type;
00041 typedef GridData<FlagType,dim> flag_data_type;
00042 typedef AMRSolverBase<VectorType,FixupType,FlagType,dim> solver_type;
00043
00044 AMRFlagging(solver_type& solver) : base(), _Solver(solver), _EstimateError(false) {}
00045
00046 virtual ~AMRFlagging() {}
00047
00048 virtual void update() {
00049 base::update();
00050 int nc;
00051 for (nc=0; nc<base::NCrit(); nc++)
00052 if (base::Criterion_(nc).ShadowCriterion() && base::Criterion_(nc).IsUsed() &&
00053 base::Criterion_(nc).UpdateShadow()) {
00054 _EstimateError = true;
00055 break;
00056 }
00057 if (!ErrorEstimation())
00058 for (nc=0; nc<base::NCrit(); nc++)
00059 if (base::Criterion_(nc).ShadowCriterion())
00060 base::Criterion_(nc).SetIsUsed(false);
00061 }
00062
00063 virtual void SetFlags(const int Time, const int Level, double t, double dt) {
00064 if (ErrorEstimation() && Time>0) {
00065
00066
00067
00068
00069 END_WATCH_FLAGGING
00070 Solver_().IntegrateLevel(U(), Time, Level, t, dt, false, 0.0, 1);
00071 START_WATCH
00072 int TStep = TimeStep(U(),Level);
00073 int ShTStep = TimeStep(Ush(),Level);
00074 forall (U(),Time,Level,c)
00075
00076
00077
00078 Ush()(Time+ShTStep,Level,c).equals(Ush()(Time,Level,c));
00079 Solver_().LevelTransfer_().Restrict(U()(Time+TStep,Level,c),Level,
00080 Ush()(Time+ShTStep,Level,c),Level,
00081 U().interiorbbox(Time+TStep,Level,c));
00082 end_forall
00083 }
00084
00085 forall (base::Flags(),Time,Level,c)
00086 base::Flags()(Time,Level,c) = GoodPoint;
00087 end_forall
00088
00089 FlagType FlagValue = BadPoint;
00090 for (int nc=0; nc<base::NCrit(); nc++) {
00091 if (!base::Criterion_(nc).IsUsed()) continue;
00092 vec_grid_fct_type* uh = (!base::Criterion_(nc).ShadowCriterion() ? _u : _u_sh);
00093 grid_fct_type* wh = (!base::Criterion_(nc).ShadowCriterion() ? _work : _work_sh);
00094 if ((uh==(vec_grid_fct_type*) 0) || (wh==(grid_fct_type*)0 )) continue;
00095 for (int cnt=0; cnt<base::Criterion_(nc).Ncnt(); cnt++) {
00096 bool Used = base::Criterion_(nc).SetFlags(*uh, *wh, base::Flags(), cnt,
00097 Time, Level, t, FlagValue);
00098 if (Used && base::Criterion_(nc).OutputFlags() && Solver_().LastOutputTime()==Time) {
00099 int TStep = TimeStep(Work(),Level);
00100 if (base::Criterion_(nc).ShadowCriterion()) {
00101 int ShTStep = TimeStep(Worksh(),Level);
00102 forall(Work(), Time, Level, c)
00103 equals_fromsh(Work()(Time+TStep,Level,c),Worksh()(Time+ShTStep,Level,c));
00104 end_forall
00105 }
00106 SwapTimeLevels(Work(),Level,Time,Time+TStep);
00107 char name[DAGHBktGFNameWidth+16];
00108 base::Criterion_(nc).OutputName(name,cnt);
00109 Solver_().FileOutput_().WritePlain(Work(), name, Time, Level, t);
00110 SwapTimeLevels(Work(),Level,Time,Time+TStep);
00111 }
00112 if (Used) FlagValue++;
00113 }
00114 }
00115 }
00116
00117 void equals_fromsh(grid_data_type& w, grid_data_type& wsh) {
00118 BBox OpBox = w.bbox();
00119 BBox OpBoxsh = wsh.bbox();
00120 Coords& Os = OpBox.stepsize();
00121 Coords& Osh = OpBoxsh.stepsize();
00122
00123 if (dim == 1) {
00124 BeginFastIndex1(w, w.bbox(), w.data(), DataType);
00125 BeginFastIndex1(wsh, wsh.bbox(), wsh.data(), DataType);
00126 for_1 (n, OpBox, Os)
00127 FastIndex1(w,n) = FastIndex1(wsh,n-n%Osh(0));
00128 end_for
00129 EndFastIndex1(w);
00130 EndFastIndex1(wsh);
00131 }
00132 else if (dim == 2) {
00133 BeginFastIndex2(w, w.bbox(), w.data(), DataType);
00134 BeginFastIndex2(wsh, wsh.bbox(), wsh.data(), DataType);
00135 for_2 (n, m, OpBox, Os)
00136 FastIndex2(w,n,m) = FastIndex2(wsh,n-n%Osh(0), m-m%Osh(1));
00137 end_for
00138 EndFastIndex2(w);
00139 EndFastIndex2(wsh);
00140 }
00141 else if (dim == 3) {
00142 BeginFastIndex3(w, w.bbox(), w.data(), DataType);
00143 BeginFastIndex3(wsh, wsh.bbox(), wsh.data(), DataType);
00144 for_3 (n, m, l, OpBox, Os)
00145 FastIndex3(w,n,m,l) = FastIndex3(wsh,n-n%Osh(0), m-m%Osh(1), l-l%Osh(2));
00146 end_for
00147 EndFastIndex3(w);
00148 EndFastIndex3(wsh);
00149 }
00150 }
00151
00152 inline void SetGridFunctions(vec_grid_fct_type* u, vec_grid_fct_type* ush,
00153 grid_fct_type* work, grid_fct_type* worksh) {
00154 _u = u; _u_sh = ush;
00155 _work = work; _work_sh = worksh;
00156 }
00157 inline vec_grid_fct_type& U() { return *_u; }
00158 inline const vec_grid_fct_type& U() const { return *_u; }
00159 inline vec_grid_fct_type& Ush() { return *_u_sh; }
00160 inline const vec_grid_fct_type& Ush() const { return *_u_sh; }
00161 inline grid_fct_type& Work() { return *_work; }
00162 inline const grid_fct_type& Work() const { return *_work; }
00163 inline grid_fct_type& Worksh() { return *_work_sh; }
00164 inline const grid_fct_type& Worksh() const { return *_work_sh; }
00165 inline const bool& ErrorEstimation() const { return _EstimateError;}
00166 void SetErrorEstimation(const bool est) { _EstimateError = est; }
00167 inline solver_type& Solver_() { return _Solver; }
00168 inline const solver_type& Solver_() const { return _Solver; }
00169
00170 protected:
00171 solver_type& _Solver;
00172 vec_grid_fct_type *_u, *_u_sh;
00173 grid_fct_type *_work, *_work_sh;
00174 bool _EstimateError;
00175 };
00176
00177
00178 #endif