00001
00002
00003
00004
00005
00006 #ifndef AMROC_RESOLVEPHI_H
00007 #define AMROC_RESOLVEPHI_H
00008
00015 #include "StdCriterion.h"
00016
00017 template <class VectorType, class FlagType, int dim> class PhiCriterion;
00018 template <class VectorType, class FixupType, class FlagType, int dim> class AMRGFMSolver;
00019
00026 template <class VectorType, class FlagType>
00027 class PhiCriterion<VectorType,FlagType,1> :
00028 public StdCriterion<VectorType,FlagType,1> {
00029 typedef typename VectorType::InternalDataType DataType;
00030 typedef StdCriterion<VectorType,FlagType,1> base;
00031 public:
00032 typedef typename base::grid_fct_type grid_fct_type;
00033 typedef typename base::grid_data_type grid_data_type;
00034
00035 PhiCriterion() : base() {}
00036
00037 virtual ~PhiCriterion() {}
00038
00039 virtual void IsolatedValue(grid_fct_type& work, const int& Time, const int& Level,
00040 int* Offset1, int* Offset2) {
00041
00042 int TStep = TimeStep(work,Level);
00043 forall(work, Time, Level, c)
00044 BBox OpBox = work.interiorbbox(Time,Level,c);
00045 Coords& OpBox_stepsize = OpBox.stepsize();
00046 BeginFastIndex1(work, work(Time,Level,c).bbox(),
00047 work(Time,Level,c).data(), DataType);
00048 BeginFastIndex1(workout, work(Time+TStep,Level,c).bbox(),
00049 work(Time+TStep,Level,c).data(), DataType);
00050
00051 Coords o1 = Coords(1,Offset1)*OpBox_stepsize;
00052 Coords o2 = Coords(1,Offset2)*OpBox_stepsize;
00053 for_1 (n, OpBox, OpBox_stepsize)
00054 if (FastIndex1(work,n+o1(0))!=FastIndex1(work,n) &&
00055 FastIndex1(work,n+o2(0))!=FastIndex1(work,n))
00056 FastIndex1(workout,n) = static_cast<DataType>(1.0);
00057 end_for
00058
00059 EndFastIndex1(work);
00060 EndFastIndex1(workout);
00061 end_forall
00062 }
00063
00064 virtual void Sgn(grid_data_type& work, const DataType& c) {
00065 BBox OpBox = work.bbox();
00066 Coords& OpBox_stepsize = OpBox.stepsize();
00067 BeginFastIndex1(work, OpBox, work.data(), DataType);
00068 for_1 (n, OpBox, OpBox_stepsize)
00069 FastIndex1(work,n) = (FastIndex1(work,n)>=c ? static_cast<DataType>(1.0) :
00070 static_cast<DataType>(-1.0));
00071 end_for
00072 EndFastIndex1(work);
00073 }
00074 };
00075
00082 template <class VectorType, class FlagType>
00083 class PhiCriterion<VectorType,FlagType,2> :
00084 public StdCriterion<VectorType,FlagType,2> {
00085 typedef typename VectorType::InternalDataType DataType;
00086 typedef StdCriterion<VectorType,FlagType,2> base;
00087 public:
00088 typedef typename base::grid_fct_type grid_fct_type;
00089 typedef typename base::grid_data_type grid_data_type;
00090
00091 PhiCriterion() : base() {}
00092
00093 virtual ~PhiCriterion() {}
00094
00095 virtual void IsolatedValue(grid_fct_type& work, const int& Time, const int& Level,
00096 int* Offset1, int* Offset2) {
00097
00098 int TStep = TimeStep(work,Level);
00099 forall(work, Time, Level, c)
00100 BBox OpBox = work.interiorbbox(Time,Level,c);
00101 Coords& OpBox_stepsize = OpBox.stepsize();
00102 BeginFastIndex2(work, work(Time,Level,c).bbox(),
00103 work(Time,Level,c).data(), DataType);
00104 BeginFastIndex2(workout, work(Time+TStep,Level,c).bbox(),
00105 work(Time+TStep,Level,c).data(), DataType);
00106
00107 Coords o1 = Coords(2,Offset1)*OpBox_stepsize;
00108 Coords o2 = Coords(2,Offset2)*OpBox_stepsize;
00109 for_2 (n, m, OpBox, OpBox_stepsize)
00110 if (FastIndex2(work,n+o1(0),m+o1(1))!=FastIndex2(work,n,m) &&
00111 FastIndex2(work,n+o2(0),m+o2(1))!=FastIndex2(work,n,m))
00112 FastIndex2(workout,n,m) = static_cast<DataType>(1.0);
00113 end_for
00114
00115 EndFastIndex2(work);
00116 EndFastIndex2(workout);
00117 end_forall
00118 }
00119
00120 virtual void Sgn(grid_data_type& work, const DataType& c) {
00121 BBox OpBox = work.bbox();
00122 Coords& OpBox_stepsize = OpBox.stepsize();
00123 BeginFastIndex2(work, OpBox, work.data(), DataType);
00124 for_2 (n, m, OpBox, OpBox_stepsize)
00125 FastIndex2(work,n,m) = (FastIndex2(work,n,m)>=c ? static_cast<DataType>(1.0) :
00126 static_cast<DataType>(-1.0));
00127 end_for
00128 EndFastIndex2(work);
00129 }
00130 };
00131
00138 template <class VectorType, class FlagType>
00139 class PhiCriterion<VectorType,FlagType,3> :
00140 public StdCriterion<VectorType,FlagType,3> {
00141 typedef typename VectorType::InternalDataType DataType;
00142 typedef StdCriterion<VectorType,FlagType,3> base;
00143 public:
00144 typedef typename base::grid_fct_type grid_fct_type;
00145 typedef typename base::grid_data_type grid_data_type;
00146
00147 PhiCriterion() : base() {}
00148
00149 virtual ~PhiCriterion() {}
00150
00151 virtual void IsolatedValue(grid_fct_type& work, const int& Time, const int& Level,
00152 int* Offset1, int* Offset2) {
00153
00154 int TStep = TimeStep(work,Level);
00155 forall(work, Time, Level, c)
00156 BBox OpBox = work.interiorbbox(Time,Level,c);
00157 Coords& OpBox_stepsize = OpBox.stepsize();
00158 BeginFastIndex3(work, work(Time,Level,c).bbox(),
00159 work(Time,Level,c).data(), DataType);
00160 BeginFastIndex3(workout, work(Time+TStep,Level,c).bbox(),
00161 work(Time+TStep,Level,c).data(), DataType);
00162
00163 Coords o1 = Coords(3,Offset1)*OpBox_stepsize;
00164 Coords o2 = Coords(3,Offset2)*OpBox_stepsize;
00165 for_3 (n, m, l, OpBox, OpBox_stepsize)
00166 if (FastIndex3(work,n+o1(0),m+o1(1),l+o1(2))!=FastIndex3(work,n,m,l) &&
00167 FastIndex3(work,n+o2(0),m+o2(1),l+o2(2))!=FastIndex3(work,n,m,l))
00168 FastIndex3(workout,n,m,l) = static_cast<DataType>(1.0);
00169 end_for
00170
00171 EndFastIndex3(work);
00172 EndFastIndex3(workout);
00173 end_forall
00174 }
00175
00176 virtual void Sgn(grid_data_type& work, const DataType& c) {
00177 BBox OpBox = work.bbox();
00178 Coords& OpBox_stepsize = OpBox.stepsize();
00179 BeginFastIndex3(work, OpBox, work.data(), DataType);
00180 for_3 (n, m, l, OpBox, OpBox_stepsize)
00181 FastIndex3(work,n,m,l) = (FastIndex3(work,n,m,l)>=c ? static_cast<DataType>(1.0) :
00182 static_cast<DataType>(-1.0));
00183 end_for
00184 EndFastIndex3(work);
00185 }
00186 };
00187
00188
00195 template <class VectorType, class Fixup, class FlagType, int dim>
00196 class ResolvePhi :
00197 public PhiCriterion<VectorType,FlagType,dim> {
00198 typedef typename VectorType::InternalDataType DataType;
00199 typedef PhiCriterion<VectorType,FlagType,dim> base;
00200 public:
00201 typedef AMRGFMSolver<VectorType,FixupType,FlagType,dim> gfm_solver_type;
00202 typedef GhostFluidMethod<VectorType,dim> gfm_type;
00203 typedef typename base::vec_grid_fct_type vec_grid_fct_type;
00204 typedef typename base::grid_fct_type grid_fct_type;
00205 typedef typename base::flag_fct_type flag_fct_type;
00206 typedef typename base::max_vector_type max_vector_type;
00207
00208 ResolvePhi(gfm_solver_type& solver, bool shadow) : base(), _solver(solver) {
00209 base::SetShadowCriterion(shadow);
00210 for (register int i=0; i<MAXCOMPONENTS; i++) {
00211 if (i<solver.NGFM()) Use[i] = 1;
00212 else Use[i] = 0;
00213 }
00214 }
00215
00216 virtual ~ResolvePhi() {}
00217
00218 virtual void register_at(ControlDevice& Ctrl) {}
00219 virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00220 if (!base::ShadowCriterion())
00221 base::LocCtrl = Ctrl.getSubDevice(prefix+"ResolvePhi");
00222 else
00223 base::LocCtrl = Ctrl.getSubDevice(prefix+"ResolvePhish");
00224 char Name[16];
00225 for (int d=0; d<max_vector_type::Length(); d++) {
00226 std::sprintf(Name,"Use(%d)",d+1);
00227 RegisterAt(base::LocCtrl,Name,Use[d]);
00228 std::sprintf(Name,"MaxLev(%d)",d+1);
00229 RegisterAt(base::LocCtrl,Name,base::_MaxLevel[d]);
00230 }
00231 RegisterAt(base::LocCtrl,"Output",base::_Output);
00232 }
00233
00234 virtual void update() {
00235 int d=max_vector_type::Length();
00236 for (; d>0; d--)
00237 if (Use[d-1] > 0)
00238 break;
00239 if (d>0) {
00240 base::SetNcnt(d);
00241 base::SetIsUsed(true);
00242 }
00243 else
00244 base::SetIsUsed(false);
00245 }
00246
00247 virtual bool SetFlags(vec_grid_fct_type& u, grid_fct_type& work, flag_fct_type& flags,
00248 const int& cnt, const int& Time, const int& Level, const double& t,
00249 const FlagType& FlagValue) {
00250
00251 if (Use[cnt] <= 0 || cnt>=NGFM() || !GFM(cnt).IsUsed()) return false;
00252 int TStep = TimeStep(work,Level);
00253 if (base::MaxLevel(cnt)<Level && base::MaxLevel(cnt)>=0) {
00254 if (base::OutputFlags()) {
00255 forall(work,Time,Level,c)
00256 work(Time+TStep,Level,c) = static_cast<DataType>(0.0);
00257 end_forall
00258 }
00259 return true;
00260 }
00261
00262 DataType md = -GFM(cnt).Boundary().Cutoff();
00263 forall(work,Time,Level,c)
00264 work(Time,Level,c).equals((GFM(cnt).Phi())(Time,Level,c));
00265 DCoords dx = base::GH().worldStep((GFM(cnt).Phi())(Time,Level,c).stepsize());
00266 if (!base::ShadowCriterion())
00267 work(Time,Level,c).equals((GFM(cnt).Phi())(Time,Level,c));
00268 else
00269 work(Time,Level,c).equals((GFM(cnt).Phish())(Time,Level,c));
00270 base::Sgn(work(Time,Level,c),md);
00271 end_forall
00272
00273 forall(work,Time,Level,c)
00274 work(Time+TStep,Level,c) = static_cast<DataType>(0.0);
00275 end_forall
00276
00277 if (dim == 1) {
00278 int East[1] = { 1 }; int West[1] = { -1 };
00279 base::IsolatedValue(work, Time, Level, East, West);
00280 }
00281 else if (dim == 2) {
00282 int East[2] = { 1, 0 }; int West[2] = { -1, 0 };
00283 int North[2] = { 0, 1 }; int South[2] = { 0, -1 };
00284 base::IsolatedValue(work, Time, Level, East, West);
00285 base::IsolatedValue(work, Time, Level, North, South);
00286 }
00287 else if (dim == 3) {
00288 int East[3] = { 1, 0, 0 }; int West[3] = {-1, 0, 0 };
00289 int North[3] = { 0, 1, 0 }; int South[3] = { 0, -1, 0 };
00290 int Top[3] = { 0, 0, 1 }; int Bottom[3] = { 0, 0, -1 };
00291 base::IsolatedValue(work, Time, Level, East, West);
00292 base::IsolatedValue(work, Time, Level, North, South);
00293 base::IsolatedValue(work, Time, Level, Top, Bottom);
00294 }
00295 int growby = (_solver.Bufferwidth()>0 ? 0 : 1);
00296 base::FlagByValue(work, flags, Time, Level, static_cast<DataType>(0.5),
00297 FlagValue, 0, growby);
00298 return true;
00299 }
00300
00301 inline gfm_type& GFM(const int& n) { return _solver.GFM(n); }
00302 inline const int& NGFM() { return _solver.NGFM(); }
00303
00304 virtual void OutputName(char* name, int cnt) {
00305 if (!base::ShadowCriterion())
00306 std::sprintf(name,"resolve_phi_%d",cnt+1);
00307 else
00308 std::sprintf(name,"resolve_phish_%d",cnt+1);
00309 }
00310
00311 protected:
00312 gfm_solver_type _solver;
00313 int Use[MAXCOMPONENTS];
00314 };
00315
00316 #endif
