00001
00002
00003
00004
00005
00006
00007
00008
00009 #ifndef AMROC_LIMITERCRITERION_H
00010 #define AMROC_LIMITERCRITERION_H
00011
00018 #include "StdCriterion.h"
00019
00020 #define UROUND 1.0e-15
00021
00022 template <class VectorType, class FlagType, int dim> class LimiterCriterion;
00023
00030 template <class VectorType, class FlagType>
00031 class LimiterCriterion<VectorType,FlagType,1> :
00032 public StdCriterion<VectorType,FlagType,1> {
00033 typedef typename VectorType::InternalDataType DataType;
00034 typedef StdCriterion<VectorType,FlagType,1> base;
00035 public:
00036 typedef typename base::grid_fct_type grid_fct_type;
00037 typedef typename base::flag_fct_type flag_fct_type;
00038 typedef typename base::grid_data_type grid_data_type;
00039
00040 LimiterCriterion() : base() {}
00041
00042 virtual ~LimiterCriterion() {}
00043
00044
00045
00046
00047 virtual DataType phi(const DataType theta) = 0;
00048
00049
00050 virtual void Limiter(grid_fct_type& work, const int& Time, const int& Level,
00051 int* Offset, const DataType& Sc) {
00052
00053 int TStep = TimeStep(work,Level);
00054 forall(work, Time, Level, c)
00055 BBox OpBox = work.interiorbbox(Time,Level,c);
00056 Coords& OpBox_stepsize = OpBox.stepsize();
00057 BeginFastIndex1(work, work(Time,Level,c).bbox(),
00058 work(Time,Level,c).data(), DataType);
00059 BeginFastIndex1(workout, work(Time+TStep,Level,c).bbox(),
00060 work(Time+TStep,Level,c).data(), DataType);
00061
00062 Coords os = Coords(1,Offset)*OpBox_stepsize;
00063 for_1 (n, OpBox, OpBox_stepsize)
00064 DataType ui = FastIndex1(work,n);
00065 DataType uim = FastIndex1(work,n-os(0));
00066 DataType uip = FastIndex1(work,n+os(0));
00067 DataType uc = static_cast<DataType>(1.0);
00068 if (std::fabs(ui-uip)>Sc || std::fabs(ui-uim)>Sc)
00069 if (std::fabs(ui-uip)>UROUND)
00070 uc = phi(std::fabs((ui-uim)/(ui-uip)) );
00071 else uc = static_cast<DataType>(0.0);
00072 if (uc < FastIndex1(workout,n))
00073 FastIndex1(workout,n) = uc;
00074 end_for
00075
00076 EndFastIndex1(work);
00077 EndFastIndex1(workout);
00078 end_forall
00079 }
00080 };
00081
00082
00089 template <class VectorType, class FlagType>
00090 class LimiterCriterion<VectorType,FlagType,2> :
00091 public StdCriterion<VectorType,FlagType,2> {
00092 typedef typename VectorType::InternalDataType DataType;
00093 typedef StdCriterion<VectorType,FlagType,2> base;
00094 public:
00095 typedef typename base::grid_fct_type grid_fct_type;
00096 typedef typename base::flag_fct_type flag_fct_type;
00097 typedef typename base::grid_data_type grid_data_type;
00098
00099 LimiterCriterion() : base() {}
00100
00101 virtual ~LimiterCriterion() {}
00102
00103
00104
00105
00106 virtual DataType phi(const DataType theta) = 0;
00107
00108
00109 virtual void Limiter(grid_fct_type& work, const int& Time, const int& Level,
00110 int* Offset, const DataType& Sc) {
00111
00112 int TStep = TimeStep(work,Level);
00113 forall(work, Time, Level, c)
00114 BBox OpBox = work.interiorbbox(Time,Level,c);
00115 Coords& OpBox_stepsize = OpBox.stepsize();
00116 BeginFastIndex2(work, work(Time,Level,c).bbox(),
00117 work(Time,Level,c).data(), DataType);
00118 BeginFastIndex2(workout, work(Time+TStep,Level,c).bbox(),
00119 work(Time+TStep,Level,c).data(), DataType);
00120
00121 Coords os = Coords(2,Offset)*OpBox_stepsize;
00122 for_2 (n, m, OpBox, OpBox_stepsize)
00123 DataType ui = FastIndex2(work,n,m);
00124 DataType uim = FastIndex2(work,n-os(0),m-os(1));
00125 DataType uip = FastIndex2(work,n+os(0),m+os(1));
00126 DataType uc = static_cast<DataType>(1.0);
00127 if (std::fabs(ui-uip)>Sc || std::fabs(ui-uim)>Sc)
00128 if (std::fabs(ui-uip)>UROUND)
00129 uc = phi(std::fabs((ui-uim)/(ui-uip)) );
00130 else uc = static_cast<DataType>(0.0);
00131 if (uc < FastIndex2(workout,n,m))
00132 FastIndex2(workout,n,m) = uc;
00133 end_for
00134
00135 EndFastIndex2(work);
00136 EndFastIndex2(workout);
00137 end_forall
00138 }
00139 };
00140
00141
00148 template <class VectorType, class FlagType>
00149 class LimiterCriterion<VectorType,FlagType,3> :
00150 public StdCriterion<VectorType,FlagType,3> {
00151 typedef typename VectorType::InternalDataType DataType;
00152 typedef StdCriterion<VectorType,FlagType,3> base;
00153 public:
00154 typedef typename base::grid_fct_type grid_fct_type;
00155 typedef typename base::flag_fct_type flag_fct_type;
00156 typedef typename base::grid_data_type grid_data_type;
00157
00158 LimiterCriterion() : base() {}
00159
00160 virtual ~LimiterCriterion() {}
00161
00162
00163
00164
00165 virtual DataType phi(const DataType theta) = 0;
00166
00167
00168 virtual void Limiter(grid_fct_type& work, const int& Time, const int& Level,
00169 int* Offset, const DataType& Sc) {
00170
00171 int TStep = TimeStep(work,Level);
00172 forall(work, Time, Level, c)
00173 BBox OpBox = work.interiorbbox(Time,Level,c);
00174 Coords& OpBox_stepsize = OpBox.stepsize();
00175 BeginFastIndex3(work, work(Time,Level,c).bbox(),
00176 work(Time,Level,c).data(), DataType);
00177 BeginFastIndex3(workout, work(Time+TStep,Level,c).bbox(),
00178 work(Time+TStep,Level,c).data(), DataType);
00179
00180 Coords os = Coords(3,Offset)*OpBox_stepsize;
00181 for_3 (n, m, l, OpBox, OpBox_stepsize)
00182 DataType ui = FastIndex3(work,n,m,l);
00183 DataType uim = FastIndex3(work,n-os(0),m-os(1),l-os(2));
00184 DataType uip = FastIndex3(work,n+os(0),m+os(1),l+os(2));
00185 DataType uc = static_cast<DataType>(1.0);
00186 if (std::fabs(ui-uip)>Sc || std::fabs(ui-uim)>Sc)
00187 if (std::fabs(ui-uip)>UROUND)
00188 uc = phi(std::fabs((ui-uim)/(ui-uip)) );
00189 else uc = static_cast<DataType>(0.0);
00190 if (uc < FastIndex3(workout,n,m,l))
00191 FastIndex3(workout,n,m,l) = uc;
00192 end_for
00193
00194 EndFastIndex3(work);
00195 EndFastIndex3(workout);
00196 end_forall
00197 }
00198 };
00199
00200
00201 #endif