00001
00002
00003
00004
00005
00006 #ifndef AMROC_PROBLEM_H
00007 #define AMROC_PROBLEM_H
00008
00009 #include "euler3.h"
00010 #include "ClpProblem.h"
00011
00012 class F77ScaledGradientSpecific : public F77ScaledGradient<VectorType,FlagType,DIM> {
00013 typedef F77ScaledGradient<VectorType,FlagType,DIM> base;
00014 typedef VectorType::InternalDataType DataType;
00015 public:
00016 typedef base::grid_fct_type grid_fct_type;
00017 typedef base::flag_fct_type flag_fct_type;
00018
00019 F77ScaledGradientSpecific(base::generic_func_type out) : base(out) {}
00020
00021 virtual bool FlagByScaledGradient(grid_fct_type& work, flag_fct_type& flags,
00022 const int& Time, const int& Level,
00023 const DataType& TolVal, const FlagType& FlagValue) {
00024 if (TolVal <= 0.0) return false;
00025
00026 int TStep = TimeStep(work,Level);
00027 forall(work,Time,Level,c)
00028 work(Time+TStep,Level,c) = static_cast<DataType>(0.0);
00029 end_forall
00030
00031 int East[3] = { 1, 0, 0 }; int West[3] = {-1, 0, 0 };
00032 int North[3] = { 0, 1, 0 }; int South[3] = { 0, -1, 0 };
00033 int Top[3] = { 0, 0, 1 }; int Bottom[3] = { 0, 0, -1 };
00034 base::Difference(work, Time, Level, East, West);
00035 base::Difference(work, Time, Level, North, South);
00036 base::Difference(work, Time, Level, Top, Bottom);
00037
00038 base::FlagByValue(work, flags, Time, Level, TolVal, FlagValue);
00039 return true;
00040 }
00041
00042 virtual void Difference(grid_fct_type& work, const int& Time, const int& Level,
00043 int* Offset1, int* Offset2) {
00044
00045 int TStep = TimeStep(work,Level);
00046 forall(work, Time, Level, c)
00047 BBox OpBox = work.interiorbbox(Time,Level,c);
00048 Coords& OpBox_stepsize = OpBox.stepsize();
00049 BeginFastIndex3(work, work(Time,Level,c).bbox(),
00050 work(Time,Level,c).data(), DataType);
00051 BeginFastIndex3(workout, work(Time+TStep,Level,c).bbox(),
00052 work(Time+TStep,Level,c).data(), DataType);
00053
00054 Coords o1 = Coords(3,Offset1)*OpBox_stepsize;
00055 Coords o2 = Coords(3,Offset2)*OpBox_stepsize;
00056 DataType value;
00057 for_3 (n, m, l, OpBox, OpBox_stepsize)
00058 value = std::max(FastIndex3(work,n+o1(0),m+o1(1),l+o1(2)),
00059 FastIndex3(work,n+o2(0),m+o2(1),l+o2(2))) /
00060 std::min(FastIndex3(work,n+o1(0),m+o1(1),l+o1(2)),
00061 FastIndex3(work,n+o2(0),m+o2(1),l+o2(2))) - 1.0;
00062 FastIndex3(workout,n,m,l) = std::max(FastIndex3(workout,n,m,l),value);
00063 end_for
00064
00065 EndFastIndex3(work);
00066 EndFastIndex3(workout);
00067 end_forall
00068 }
00069 };
00070
00071 class FlaggingSpecific :
00072 public AMRFlagging<VectorType,FixupType,FlagType,DIM> {
00073 typedef AMRFlagging<VectorType,FixupType,FlagType,DIM> base;
00074 public:
00075 FlaggingSpecific(base::solver_type& solver) : base(solver) {
00076 base::AddCriterion(new F77ScaledGradientSpecific(f_flgout));
00077 }
00078
00079 ~FlaggingSpecific() { DeleteAllCriterions(); }
00080
00081 virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00082 base::register_at(Ctrl, prefix);
00083 char VariableName[32];
00084 for (int d=0; d<DIM; d++) {
00085 sprintf(VariableName,"DMin(%d)",d+1);
00086 RegisterAt(base::LocCtrl,VariableName,dcmin[d]);
00087 sprintf(VariableName,"DMax(%d)",d+1);
00088 RegisterAt(base::LocCtrl,VariableName,dcmax[d]);
00089 }
00090 }
00091 virtual void register_at(ControlDevice& Ctrl) {
00092 register_at(Ctrl, "");
00093 }
00094
00095 virtual void SetFlags(const int Time, const int Level, double t, double dt) {
00096 base::SetFlags(Time, Level, t, dt);
00097 Coords lb = base::GH().localCoords(dcmin);
00098 Coords ub = base::GH().localCoords(dcmax);
00099 Coords ss(base::Dim(),StepSize(base::GH(),Level));
00100 lb /= ss; lb *= ss;
00101 ub /= ss; ub *= ss;
00102 BBox bbox(lb,ub,ss);
00103 forall (base::Flags(),Time,Level,c)
00104 base::Flags()(Time,Level,c).equals(BadPoint, bbox);
00105 end_forall
00106 }
00107
00108 private:
00109 double dcmin[DIM], dcmax[DIM];
00110 };
00111
00112 #define OWN_FLAGGING