amroc/amr/Criteria/ResolvePhi.h
Go to the documentation of this file.00001 // -*- C++ -*- 00002 00003 // Copyright (C) 2003-2007 California Institute of Technology 00004 // Ralf Deiterding, ralf@cacr.caltech.edu 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,"res_phi_%d",cnt+1); 00307 else 00308 std::sprintf(name,"res_phish_%d",cnt+1); 00309 } 00310 00311 protected: 00312 gfm_solver_type _solver; 00313 int Use[MAXCOMPONENTS]; 00314 }; 00315 00316 #endif
Generated on Thu Jun 30 02:19:16 2016 for AMROC Fluid-solver Framework by
1.6.3