amroc/amr/Criteria/MRVectorPrediction.h

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 
00003 // Copyright (C) 2015 Ralf Deiterding
00004 // University of Southampton
00005 
00006 #ifndef AMROC_MRPREDICTION_H
00007 #define AMROC_MRPREDICTION_H
00008 
00015 #include "StdCriterion.h"
00016 
00023 template <class VectorType, class FlagType, int dim>
00024 class MRPrediction : 
00025   public StdCriterion<VectorType,FlagType,dim> {
00026   typedef typename VectorType::InternalDataType DataType;
00027   typedef StdCriterion<VectorType,FlagType,dim> base;
00028 public:
00029   typedef typename base::vec_grid_fct_type vec_grid_fct_type;  
00030   typedef typename base::grid_fct_type grid_fct_type;
00031   typedef typename base::flag_fct_type flag_fct_type;
00032   typedef typename base::max_vector_type max_vector_type;
00033   typedef typename base::grid_data_type grid_data_type;
00034 
00035   MRPrediction() : base() {
00036     base::SetShadowCriterion(false);
00037     Tol = static_cast<DataType>(0.0);
00038   }
00039 
00040   virtual ~MRPrediction() {}
00041   
00042   virtual void register_at(ControlDevice& Ctrl) {}
00043   virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00044     base::LocCtrl = Ctrl.getSubDevice(prefix+"MRPrediction");
00045     char Name[16];
00046     for (int d=0; d<max_vector_type::Length(); d++) {
00047       std::sprintf(Name,"Tol(%d)",d+1);
00048       RegisterAt(base::LocCtrl,Name,Tol(d));
00049       std::sprintf(Name,"MaxLev(%d)",d+1);
00050       RegisterAt(base::LocCtrl,Name,base::_MaxLevel[d]);
00051     }    
00052     RegisterAt(base::LocCtrl,"Output",base::_Output);
00053   }
00054 
00055   virtual void update() { 
00056     int d=max_vector_type::Length();
00057     for (; d>0; d--) 
00058       if (Tol(d-1) > static_cast<DataType>(0.0)) 
00059         break;
00060     if (d>0) {
00061       base::SetNcnt(d);
00062       base::SetIsUsed(true);
00063     }
00064     else 
00065       base::SetIsUsed(false);
00066   }
00067 
00068   virtual bool SetFlags(vec_grid_fct_type& u, grid_fct_type& work, flag_fct_type& flags, 
00069                         const int& cnt, const int& Time, const int& Level, const double& t, 
00070                         const FlagType& FlagValue) {
00071     if (Tol(cnt) <= 0.0) return false;
00072     int TStep = TimeStep(work,Level);
00073     if (base::MaxLevel(cnt)<Level && base::MaxLevel(cnt)>=0) {
00074       if (base::OutputFlags()) {
00075         forall(work,Time,Level,c)
00076           work(Time+TStep,Level,c) = static_cast<DataType>(0.0);
00077         end_forall 
00078       }
00079       return true;
00080     }
00081 
00082     forall(work,Time,Level,c)
00083       equals_from(work(Time,Level,c), u(Time,Level,c), cnt);
00084     end_forall 
00085     return FlagByMRPrediction(work, flags, Time, Level, Tol(cnt), FlagValue);
00086   }
00087 
00088   virtual bool FlagByMRPrediction(grid_fct_type& work, flag_fct_type& flags,
00089                                   const int& Time, const int& Level, 
00090                                   const DataType& TolVal, const FlagType& FlagValue) { 
00091     if (TolVal <= 0.0) return false;
00092 
00093     int TStep = TimeStep(work,Level);
00094     forall(work,Time,Level,c)
00095       BBox bb = work(Time,Level,c).bbox();
00096       BBox bbi = work.interiorbbox(Time,Level,c);
00097       BBox avbb = coarsen(bb,2);
00098       Coords avss = avbb.stepsize();
00099       grid_data_type average(avbb);
00100 
00101       if (dim == 1) {
00102         BeginFastIndex1(work, bb, work(Time,Level,c).data(), DataType);
00103         BeginFastIndex1(workn, bb, work(Time+TStep,Level,c).data(), DataType);
00104         BeginFastIndex1(average, avbb, average.data(), DataType);
00105         int dn=avss(0);  
00106         for_1 (n,avbb,avss)
00107           FastIndex1(average,n)=0.;
00108           BBox bbadd(1,n,n,dn);
00109           bbadd.refine(2); bbadd*=bb;
00110           for_1 (nn,bbadd,bbadd.stepsize())
00111             FastIndex1(average,n)+=FastIndex1(work,nn);
00112           end_for
00113           FastIndex1(average,n)/=bbadd.size();
00114         end_for
00115 
00116         for_1 (nn,bbi,bbi.stepsize())
00117           int n=(nn/dn)*dn;     
00118           FastIndex1(workn,nn) = FastIndex1(average,n);
00119 
00120           DataType pn = (n==nn) ? static_cast<DataType>(1.): static_cast<DataType>(-1.);
00121           FastIndex1(workn,nn) += (pn*-static_cast<DataType>(.125)) * FastIndex1(average,n+dn);
00122           FastIndex1(workn,nn) -= (pn*-static_cast<DataType>(.125)) * FastIndex1(average,n-dn);
00123 
00124           FastIndex1(workn,nn) = std::fabs(FastIndex1(workn,nn)-FastIndex1(work,nn));
00125         end_for
00126         EndFastIndex1(work); 
00127         EndFastIndex1(workn); 
00128         EndFastIndex1(average); 
00129       }
00130       
00131       else if (dim == 2) {
00132         BeginFastIndex2(work, bb, work(Time,Level,c).data(), DataType);
00133         BeginFastIndex2(workn, bb, work(Time+TStep,Level,c).data(), DataType);
00134         BeginFastIndex2(average, avbb, average.data(), DataType);
00135         int dn=avss(0), dm=avss(1);  
00136         for_2 (n,m,avbb,avss)
00137           FastIndex2(average,n,m)=0.;
00138           BBox bbadd(2,n,m,n,m,dn,dm);
00139           bbadd.refine(2); bbadd*=bb;
00140           for_2 (nn,mm,bbadd,bbadd.stepsize())
00141             FastIndex2(average,n,m)+=FastIndex2(work,nn,mm);
00142           end_for
00143           FastIndex2(average,n,m)/=bbadd.size();
00144         end_for
00145 
00146         for_2 (nn,mm,bbi,bbi.stepsize())
00147           int n=(nn/dn)*dn;
00148           int m=(mm/dm)*dm;
00149         
00150           FastIndex2(workn,nn,mm) = FastIndex2(average,n,m);
00151           
00152           DataType pn = (n==nn) ? static_cast<DataType>(1.): static_cast<DataType>(-1.);
00153           FastIndex2(workn,nn,mm) += (pn*-static_cast<DataType>(.125)) * FastIndex2(average,n+dn,m);
00154           FastIndex2(workn,nn,mm) -= (pn*-static_cast<DataType>(.125)) * FastIndex2(average,n-dn,m);
00155           
00156           DataType pm = (m==mm) ? static_cast<DataType>(1.): static_cast<DataType>(-1.);        
00157           FastIndex2(workn,nn,mm) += (pm*-static_cast<DataType>(.125)) * FastIndex2(average,n,m+dm); 
00158           FastIndex2(workn,nn,mm) -= (pm*-static_cast<DataType>(.125)) * FastIndex2(average,n,m-dm);
00159 
00160           FastIndex2(workn,nn,mm) += (pn*pm*static_cast<DataType>(.015625)) * FastIndex2(average,n+dn,m+dm);
00161           FastIndex2(workn,nn,mm) -= (pn*pm*static_cast<DataType>(.015625)) * FastIndex2(average,n+dn,m-dm);
00162           FastIndex2(workn,nn,mm) -= (pn*pm*static_cast<DataType>(.015625)) * FastIndex2(average,n-dn,m+dm);
00163           FastIndex2(workn,nn,mm) += (pn*pm*static_cast<DataType>(.015625)) * FastIndex2(average,n-dn,m-dm);
00164 
00165           FastIndex2(workn,nn,mm) = std::fabs(FastIndex2(workn,nn,mm)-FastIndex2(work,nn,mm));
00166         end_for
00167         EndFastIndex2(work); 
00168         EndFastIndex2(workn); 
00169         EndFastIndex2(average); 
00170       }
00171       
00172       else if (dim == 3) {
00173         BeginFastIndex3(work, bb, work(Time,Level,c).data(), DataType);
00174         BeginFastIndex3(workn, bb, work(Time+TStep,Level,c).data(), DataType);
00175         BeginFastIndex3(average, avbb, average.data(), DataType);
00176         int dn=avss(0), dm=avss(1), dl=avss(2);  
00177         for_3 (n,m,l,avbb,avss)
00178           FastIndex3(average,n,m,l)=0.;
00179           BBox bbadd(3,n,m,l,n,m,l,dn,dm,dl);
00180           bbadd.refine(2); bbadd*=bb;
00181           for_3 (nn,mm,ll,bbadd,bbadd.stepsize())
00182             FastIndex3(average,n,m,l)+=FastIndex3(work,nn,mm,ll);
00183           end_for
00184           FastIndex3(average,n,m,l)/=bbadd.size();
00185         end_for
00186 
00187         for_3 (nn,mm,ll,bbi,bbi.stepsize())
00188           int n=(nn/avss(0))*dn;
00189           int m=(mm/avss(1))*dm;
00190           int l=(mm/avss(2))*dl;
00191           
00192           FastIndex3(workn,nn,mm,ll) = FastIndex3(average,n,m,l);
00193           
00194           DataType pn = (n==nn) ? static_cast<DataType>(1.): static_cast<DataType>(-1.);
00195           FastIndex3(workn,nn,mm,ll) += (pn*-static_cast<DataType>(.125)) * FastIndex3(average,n+dn,m,l);
00196           FastIndex3(workn,nn,mm,ll) -= (pn*-static_cast<DataType>(.125)) * FastIndex3(average,n-dn,m,l);
00197           
00198           DataType pm = (m==mm) ? static_cast<DataType>(1.): static_cast<DataType>(-1.);        
00199           FastIndex3(workn,nn,mm,ll) += (pm*-static_cast<DataType>(.125)) * FastIndex3(average,n,m+dm,l); 
00200           FastIndex3(workn,nn,mm,ll) -= (pm*-static_cast<DataType>(.125)) * FastIndex3(average,n,m-dm,l);
00201           
00202           DataType pl = (l==ll) ? static_cast<DataType>(1.): static_cast<DataType>(-1.);        
00203           FastIndex3(workn,nn,mm,ll) += (pl*-static_cast<DataType>(.125)) * FastIndex3(average,n,m,l+dl); 
00204           FastIndex3(workn,nn,mm,ll) -= (pl*-static_cast<DataType>(.125)) * FastIndex3(average,n,m,l-dl);
00205           
00206           FastIndex3(workn,nn,mm,ll) += (pn*pm*static_cast<DataType>(.015625)) * FastIndex3(average,n+dn,m+dm,l);
00207           FastIndex3(workn,nn,mm,ll) -= (pn*pm*static_cast<DataType>(.015625)) * FastIndex3(average,n+dn,m-dm,l);
00208           FastIndex3(workn,nn,mm,ll) -= (pn*pm*static_cast<DataType>(.015625)) * FastIndex3(average,n-dn,m+dm,l);
00209           FastIndex3(workn,nn,mm,ll) += (pn*pm*static_cast<DataType>(.015625)) * FastIndex3(average,n-dn,m-dm,l);
00210           
00211           FastIndex3(workn,nn,mm,ll) += (pn*pl*static_cast<DataType>(.015625)) * FastIndex3(average,n+dn,m,l+dl);
00212           FastIndex3(workn,nn,mm,ll) -= (pn*pl*static_cast<DataType>(.015625)) * FastIndex3(average,n+dn,m,l-dl);
00213           FastIndex3(workn,nn,mm,ll) -= (pn*pl*static_cast<DataType>(.015625)) * FastIndex3(average,n-dn,m,l+dl);
00214           FastIndex3(workn,nn,mm,ll) += (pn*pl*static_cast<DataType>(.015625)) * FastIndex3(average,n-dn,m,l-dl);
00215 
00216           FastIndex3(workn,nn,mm,ll) += (pm*pl*static_cast<DataType>(.015625)) * FastIndex3(average,n,m+dm,l+dl);
00217           FastIndex3(workn,nn,mm,ll) -= (pm*pl*static_cast<DataType>(.015625)) * FastIndex3(average,n,m+dm,l-dl);
00218           FastIndex3(workn,nn,mm,ll) -= (pm*pl*static_cast<DataType>(.015625)) * FastIndex3(average,n,m-dm,l+dl);
00219           FastIndex3(workn,nn,mm,ll) += (pm*pl*static_cast<DataType>(.015625)) * FastIndex3(average,n,m-dm,l-dl);
00220 
00221           FastIndex3(workn,nn,mm,ll) += (pn*pm*pl-static_cast<DataType>(.001953125)) * FastIndex3(average,n+dn,m+dm,l+dl); 
00222           FastIndex3(workn,nn,mm,ll) -= (pn*pm*pl-static_cast<DataType>(.001953125)) * FastIndex3(average,n+dn,m+dm,l-dl);
00223           FastIndex3(workn,nn,mm,ll) -= (pn*pm*pl-static_cast<DataType>(.001953125)) * FastIndex3(average,n+dn,m-dm,l+dl);
00224           FastIndex3(workn,nn,mm,ll) += (pn*pm*pl-static_cast<DataType>(.001953125)) * FastIndex3(average,n+dn,m-dm,l-dl);
00225           FastIndex3(workn,nn,mm,ll) -= (pn*pm*pl-static_cast<DataType>(.001953125)) * FastIndex3(average,n-dn,m+dm,l+dl);
00226           FastIndex3(workn,nn,mm,ll) += (pn*pm*pl-static_cast<DataType>(.001953125)) * FastIndex3(average,n-dn,m+dm,l-dl);
00227           FastIndex3(workn,nn,mm,ll) += (pn*pm*pl-static_cast<DataType>(.001953125)) * FastIndex3(average,n-dn,m-dm,l+dl);
00228           FastIndex3(workn,nn,mm,ll) -= (pn*pm*pl-static_cast<DataType>(.001953125)) * FastIndex3(average,n-dn,m-dm,l-dl);
00229 
00230           FastIndex3(workn,nn,mm,ll) = std::fabs(FastIndex3(workn,nn,mm,ll)-FastIndex3(work,nn,mm,ll));
00231         end_for
00232         EndFastIndex3(work); 
00233         EndFastIndex3(workn); 
00234         EndFastIndex3(average); 
00235       }
00236     end_forall 
00237       
00238     base::FlagByValue(work, flags, Time, Level, TolVal, FlagValue);
00239     return true;
00240   }
00241 
00242   virtual void OutputName(char* name, int cnt) { std::sprintf(name,"mrpredict_%d",cnt+1); }
00243 
00244   max_vector_type Tol;
00245 };
00246 
00247 #endif
Generated on Thu Jun 30 02:19:15 2016 for AMROC Fluid-solver Framework by 1.6.3