amroc/amr/Criteria/MRPrediction.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(), GlobalVolume(1.) {
00036     base::SetShadowCriterion(false);
00037     Tol = static_cast<DataType>(0.0);
00038     for (register int i=0; i<MAXCOMPONENTS; i++) {
00039       ThresholdNorm[i] = 0; Relative[i] = 0;
00040     }
00041   }
00042 
00043   virtual ~MRPrediction() {}
00044   
00045   virtual void register_at(ControlDevice& Ctrl) {}
00046   virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00047     base::LocCtrl = Ctrl.getSubDevice(prefix+"MRPrediction");
00048     char Name[32];
00049     for (int d=0; d<max_vector_type::Length(); d++) {
00050       std::sprintf(Name,"Tol(%d)",d+1);
00051       RegisterAt(base::LocCtrl,Name,Tol(d));
00052       std::sprintf(Name,"MaxLev(%d)",d+1);
00053       RegisterAt(base::LocCtrl,Name,base::_MaxLevel[d]);
00054       std::sprintf(Name,"ThresholdNorm(%d)",d+1);
00055       RegisterAt(base::LocCtrl,Name,ThresholdNorm[d]);
00056       std::sprintf(Name,"Relative(%d)",d+1);
00057       RegisterAt(base::LocCtrl,Name,Relative[d]);
00058     }    
00059     RegisterAt(base::LocCtrl,"Output",base::_Output);
00060   }
00061 
00062   virtual void update() { 
00063     int d=max_vector_type::Length();
00064     for (; d>0; d--) 
00065       if (Tol(d-1) > static_cast<DataType>(0.0)) 
00066         break;
00067     if (d>0) {
00068       base::SetNcnt(d);
00069       base::SetIsUsed(true);
00070     }
00071     else 
00072       base::SetIsUsed(false);
00073   }
00074 
00075   virtual void SetupData(GridHierarchy* gh, const int& ghosts) {
00076     base::SetupData(gh,ghosts);
00077     if (base::IsUsed()) {
00078       const double* bndry = base::GH().wholebndry();
00079       for (int d=0; d<dim; d++)      
00080         GlobalVolume *= bndry[1+dim*d]-bndry[0+dim*d];
00081     }
00082   }
00083   
00084   virtual bool SetFlags(vec_grid_fct_type& u, grid_fct_type& work, flag_fct_type& flags, 
00085                         const int& cnt, const int& Time, const int& Level, const double& t, 
00086                         const FlagType& FlagValue) {
00087     if (Tol(cnt) <= 0.0) return false;
00088     int TStep = TimeStep(work,Level);
00089     if (base::MaxLevel(cnt)<Level && base::MaxLevel(cnt)>=0) {
00090       if (base::OutputFlags()) {
00091         forall(work,Time,Level,c)
00092           work(Time+TStep,Level,c) = static_cast<DataType>(0.0);
00093         end_forall 
00094       }
00095       return true;
00096     }
00097 
00098     forall(work,Time,Level,c)
00099       equals_from(work(Time,Level,c), u(Time,Level,c), cnt);
00100     end_forall 
00101     return FlagByMRPrediction(work, flags, Time, Level, Tol[cnt], ThresholdNorm[cnt], Relative[cnt], FlagValue);
00102   }
00103 
00104   virtual bool FlagByMRPrediction(grid_fct_type& work, flag_fct_type& flags,
00105                                   const int& Time, const int& Level, 
00106                                   const DataType& TolVal, const int& ThresholdVal,
00107                                   const int& RelVal, const FlagType& FlagValue) { 
00108     if (TolVal <= 0.0) return false;
00109 
00110     DataType Tol = TolVal;
00111     if (ThresholdVal)
00112       Tol *= std::exp(dim*(Level-MaxLevel(base::GH())+1)*std::log(static_cast<DataType>(2.)))/GlobalVolume;
00113     DataType sc(1.);
00114     if (RelVal) {
00115       DataType max=MaxVal(work,Time,Level);
00116       DataType min=MinVal(work,Time,Level);
00117       sc = std::max(std::fabs(max),std::fabs(min));
00118     }
00119 
00120     int TStep = TimeStep(work,Level);
00121     forall(work,Time,Level,c)
00122       BBox bb = work(Time,Level,c).bbox();
00123       BBox bbi = work.interiorbbox(Time,Level,c);
00124       BBox avbb = coarsen(bb,2);
00125       Coords avss = avbb.stepsize();
00126       grid_data_type average(avbb);
00127 
00128       if (dim == 1) {
00129         BeginFastIndex1(work, bb, work(Time,Level,c).data(), DataType);
00130         BeginFastIndex1(workn, bb, work(Time+TStep,Level,c).data(), DataType);
00131         BeginFastIndex1(average, avbb, average.data(), DataType);
00132         int dn=avss(0);  
00133         for_1 (n,avbb,avss)
00134           FastIndex1(average,n)=0.;
00135           BBox bbadd(1,n,n,dn);
00136           bbadd.refine(2); bbadd*=bb;
00137           for_1 (nn,bbadd,bbadd.stepsize())
00138             FastIndex1(average,n)+=FastIndex1(work,nn);
00139           end_for
00140           FastIndex1(average,n)/=bbadd.size();
00141         end_for
00142 
00143         for_1 (nn,bbi,bbi.stepsize())
00144           int n=(nn/dn)*dn;
00145           DataType pred = FastIndex1(average,n);
00146 
00147           DataType pn = (n==nn) ? static_cast<DataType>(1.): static_cast<DataType>(-1.);
00148           pred += (pn*-static_cast<DataType>(.125)) * FastIndex1(average,n+dn);
00149           pred -= (pn*-static_cast<DataType>(.125)) * FastIndex1(average,n-dn);
00150 
00151           FastIndex1(workn,nn) = std::fabs(pred-FastIndex1(work,nn))/sc;
00152         end_for
00153         EndFastIndex1(work); 
00154         EndFastIndex1(workn); 
00155         EndFastIndex1(average); 
00156       }
00157       
00158       else if (dim == 2) {
00159         BeginFastIndex2(work, bb, work(Time,Level,c).data(), DataType);
00160         BeginFastIndex2(workn, bb, work(Time+TStep,Level,c).data(), DataType);
00161         BeginFastIndex2(average, avbb, average.data(), DataType);
00162         int dn=avss(0), dm=avss(1);  
00163         for_2 (n,m,avbb,avss)
00164           FastIndex2(average,n,m)=0.;
00165           BBox bbadd(2,n,m,n,m,dn,dm);
00166           bbadd.refine(2); bbadd*=bb;
00167           for_2 (nn,mm,bbadd,bbadd.stepsize())
00168             FastIndex2(average,n,m)+=FastIndex2(work,nn,mm);
00169           end_for
00170           FastIndex2(average,n,m)/=bbadd.size();
00171         end_for
00172 
00173         for_2 (nn,mm,bbi,bbi.stepsize())
00174           int n=(nn/dn)*dn;
00175           int m=(mm/dm)*dm;
00176         
00177           DataType pred = FastIndex2(average,n,m);
00178           
00179           DataType pn = (n==nn) ? static_cast<DataType>(1.): static_cast<DataType>(-1.);
00180           pred += (pn*-static_cast<DataType>(.125)) * FastIndex2(average,n+dn,m);
00181           pred -= (pn*-static_cast<DataType>(.125)) * FastIndex2(average,n-dn,m);
00182           
00183           DataType pm = (m==mm) ? static_cast<DataType>(1.): static_cast<DataType>(-1.);        
00184           pred += (pm*-static_cast<DataType>(.125)) * FastIndex2(average,n,m+dm); 
00185           pred -= (pm*-static_cast<DataType>(.125)) * FastIndex2(average,n,m-dm);
00186 
00187           pred += (pn*pm*static_cast<DataType>(.015625)) * FastIndex2(average,n+dn,m+dm);
00188           pred -= (pn*pm*static_cast<DataType>(.015625)) * FastIndex2(average,n+dn,m-dm);
00189           pred -= (pn*pm*static_cast<DataType>(.015625)) * FastIndex2(average,n-dn,m+dm);
00190           pred += (pn*pm*static_cast<DataType>(.015625)) * FastIndex2(average,n-dn,m-dm);
00191 
00192           FastIndex2(workn,nn,mm) = std::fabs(pred-FastIndex2(work,nn,mm))/sc;
00193         end_for
00194         EndFastIndex2(work); 
00195         EndFastIndex2(workn); 
00196         EndFastIndex2(average); 
00197       }
00198       
00199       else if (dim == 3) {
00200         BeginFastIndex3(work, bb, work(Time,Level,c).data(), DataType);
00201         BeginFastIndex3(workn, bb, work(Time+TStep,Level,c).data(), DataType);
00202         BeginFastIndex3(average, avbb, average.data(), DataType);
00203         int dn=avss(0), dm=avss(1), dl=avss(2);  
00204         for_3 (n,m,l,avbb,avss)
00205           FastIndex3(average,n,m,l)=0.;
00206           BBox bbadd(3,n,m,l,n,m,l,dn,dm,dl);
00207           bbadd.refine(2); bbadd*=bb;
00208           for_3 (nn,mm,ll,bbadd,bbadd.stepsize())
00209             FastIndex3(average,n,m,l)+=FastIndex3(work,nn,mm,ll);
00210           end_for
00211           FastIndex3(average,n,m,l)/=bbadd.size();
00212         end_for
00213 
00214         for_3 (nn,mm,ll,bbi,bbi.stepsize())
00215           int n=(nn/avss(0))*dn;
00216           int m=(mm/avss(1))*dm;
00217           int l=(mm/avss(2))*dl;
00218           
00219           DataType pred = FastIndex3(average,n,m,l);
00220           
00221           DataType pn = (n==nn) ? static_cast<DataType>(1.): static_cast<DataType>(-1.);
00222           pred += (pn*-static_cast<DataType>(.125)) * FastIndex3(average,n+dn,m,l);
00223           pred -= (pn*-static_cast<DataType>(.125)) * FastIndex3(average,n-dn,m,l);
00224           
00225           DataType pm = (m==mm) ? static_cast<DataType>(1.): static_cast<DataType>(-1.);        
00226           pred += (pm*-static_cast<DataType>(.125)) * FastIndex3(average,n,m+dm,l); 
00227           pred -= (pm*-static_cast<DataType>(.125)) * FastIndex3(average,n,m-dm,l);
00228           
00229           DataType pl = (l==ll) ? static_cast<DataType>(1.): static_cast<DataType>(-1.);        
00230           pred += (pl*-static_cast<DataType>(.125)) * FastIndex3(average,n,m,l+dl); 
00231           pred -= (pl*-static_cast<DataType>(.125)) * FastIndex3(average,n,m,l-dl);
00232           
00233           pred += (pn*pm*static_cast<DataType>(.015625)) * FastIndex3(average,n+dn,m+dm,l);
00234           pred -= (pn*pm*static_cast<DataType>(.015625)) * FastIndex3(average,n+dn,m-dm,l);
00235           pred -= (pn*pm*static_cast<DataType>(.015625)) * FastIndex3(average,n-dn,m+dm,l);
00236           pred += (pn*pm*static_cast<DataType>(.015625)) * FastIndex3(average,n-dn,m-dm,l);
00237           
00238           pred += (pn*pl*static_cast<DataType>(.015625)) * FastIndex3(average,n+dn,m,l+dl);
00239           pred -= (pn*pl*static_cast<DataType>(.015625)) * FastIndex3(average,n+dn,m,l-dl);
00240           pred -= (pn*pl*static_cast<DataType>(.015625)) * FastIndex3(average,n-dn,m,l+dl);
00241           pred += (pn*pl*static_cast<DataType>(.015625)) * FastIndex3(average,n-dn,m,l-dl);
00242 
00243           pred += (pm*pl*static_cast<DataType>(.015625)) * FastIndex3(average,n,m+dm,l+dl);
00244           pred -= (pm*pl*static_cast<DataType>(.015625)) * FastIndex3(average,n,m+dm,l-dl);
00245           pred -= (pm*pl*static_cast<DataType>(.015625)) * FastIndex3(average,n,m-dm,l+dl);
00246           pred += (pm*pl*static_cast<DataType>(.015625)) * FastIndex3(average,n,m-dm,l-dl);
00247 
00248           pred += (pn*pm*pl-static_cast<DataType>(.001953125)) * FastIndex3(average,n+dn,m+dm,l+dl); 
00249           pred -= (pn*pm*pl-static_cast<DataType>(.001953125)) * FastIndex3(average,n+dn,m+dm,l-dl);
00250           pred -= (pn*pm*pl-static_cast<DataType>(.001953125)) * FastIndex3(average,n+dn,m-dm,l+dl);
00251           pred += (pn*pm*pl-static_cast<DataType>(.001953125)) * FastIndex3(average,n+dn,m-dm,l-dl);
00252           pred -= (pn*pm*pl-static_cast<DataType>(.001953125)) * FastIndex3(average,n-dn,m+dm,l+dl);
00253           pred += (pn*pm*pl-static_cast<DataType>(.001953125)) * FastIndex3(average,n-dn,m+dm,l-dl);
00254           pred += (pn*pm*pl-static_cast<DataType>(.001953125)) * FastIndex3(average,n-dn,m-dm,l+dl);
00255           pred -= (pn*pm*pl-static_cast<DataType>(.001953125)) * FastIndex3(average,n-dn,m-dm,l-dl);
00256 
00257           FastIndex3(workn,nn,mm,ll) = std::fabs(pred-FastIndex3(work,nn,mm,ll))/sc;
00258         end_for
00259         EndFastIndex3(work); 
00260         EndFastIndex3(workn); 
00261         EndFastIndex3(average); 
00262       }
00263     end_forall 
00264 
00265     base::FlagByValue(work, flags, Time, Level, Tol, FlagValue);
00266     return true;
00267   }
00268 
00269   virtual void OutputName(char* name, int cnt) { std::sprintf(name,"mrpr_%d",cnt+1); }
00270 
00271   max_vector_type Tol;
00272 
00273   int ThresholdNorm[MAXCOMPONENTS], Relative[MAXCOMPONENTS];
00274 
00275 protected:
00276   double GlobalVolume;
00277 };
00278 
00279 
00286 template <class VectorType, class FlagType, int dim>
00287 class MRVectorPrediction : 
00288   public StdCriterion<VectorType,FlagType,dim> {
00289   typedef typename VectorType::InternalDataType DataType;
00290   typedef StdCriterion<VectorType,FlagType,dim> base;
00291 public:
00292   typedef typename base::vec_grid_fct_type vec_grid_fct_type;  
00293   typedef typename base::vec_grid_data_type vec_grid_data_type;  
00294   typedef typename base::grid_fct_type grid_fct_type;
00295   typedef typename base::flag_fct_type flag_fct_type;
00296   typedef typename base::max_vector_type max_vector_type;
00297   typedef typename base::grid_data_type grid_data_type;
00298 
00299   MRVectorPrediction() : base(), MaxMomentum(1), GlobalVolume(1.) {
00300     base::SetShadowCriterion(false);
00301     Tol = static_cast<DataType>(0.0);
00302     for (register int i=0; i<MAXCOMPONENTS; i++) {
00303       ThresholdNorm[i] = 0; Relative[i] = 0;
00304     }
00305   }
00306 
00307   virtual ~MRVectorPrediction() {}
00308    
00309   virtual void register_at(ControlDevice& Ctrl) {}
00310   virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00311     base::LocCtrl = Ctrl.getSubDevice(prefix+"MRVectorPrediction");
00312     char Name[32];
00313     for (int d=0; d<max_vector_type::Length(); d++) {
00314       std::sprintf(Name,"Tol(%d)",d+1);
00315       RegisterAt(base::LocCtrl,Name,Tol(d));
00316       std::sprintf(Name,"MaxLev(%d)",d+1);
00317       RegisterAt(base::LocCtrl,Name,base::_MaxLevel[d]);
00318       std::sprintf(Name,"ThresholdNorm(%d)",d+1);
00319       RegisterAt(base::LocCtrl,Name,ThresholdNorm[d]);
00320       std::sprintf(Name,"Relative(%d)",d+1);
00321       RegisterAt(base::LocCtrl,Name,Relative[d]);
00322     }    
00323     RegisterAt(base::LocCtrl,"Output",base::_Output);
00324     RegisterAt(base::LocCtrl,"MaxMomentum",MaxMomentum);
00325   }
00326 
00327   virtual void update() { 
00328     int d=max_vector_type::Length();
00329     for (; d>0; d--) 
00330       if (Tol(d-1) > static_cast<DataType>(0.0)) 
00331         break;
00332     if (d>0) {
00333       base::SetNcnt(d);
00334       base::SetIsUsed(true);
00335     }
00336     else 
00337       base::SetIsUsed(false);
00338   }
00339 
00340   virtual void SetupData(GridHierarchy* gh, const int& ghosts) {
00341     base::SetupData(gh,ghosts);
00342     if (base::IsUsed()) {
00343       const double* bndry = base::GH().wholebndry();
00344       for (int d=0; d<dim; d++)      
00345         GlobalVolume *= bndry[1+dim*d]-bndry[0+dim*d];
00346     }
00347   }
00348   
00349   virtual bool SetFlags(vec_grid_fct_type& u, grid_fct_type& work, flag_fct_type& flags, 
00350                         const int& cnt, const int& Time, const int& Level, const double& t, 
00351                         const FlagType& FlagValue) {
00352     if (Tol(cnt)<=0.0 || cnt!=0) return false;
00353     int TStep = TimeStep(work,Level);
00354     if (base::MaxLevel(cnt)<Level && base::MaxLevel(cnt)>=0) {
00355       if (base::OutputFlags()) {
00356         forall(work,Time,Level,c)
00357           work(Time+TStep,Level,c) = static_cast<DataType>(0.0);
00358         end_forall 
00359       }
00360       return true;
00361     }
00362 
00363     return FlagByMRVectorPrediction(u, work, flags, Time, Level, Tol[cnt], ThresholdNorm[cnt], Relative[cnt], FlagValue);
00364   }
00365 
00366   virtual bool FlagByMRVectorPrediction(vec_grid_fct_type& u, grid_fct_type& work,
00367                                         flag_fct_type& flags, const int& Time, const int& Level, 
00368                                         const DataType& TolVal, const int& ThresholdVal,
00369                                         const int& RelVal, const FlagType& FlagValue) { 
00370     if (TolVal <= 0.0) return false;
00371 
00372     DataType Tol = TolVal;
00373     if (ThresholdVal)
00374       Tol *= std::exp(dim*(Level-MaxLevel(base::GH())+1)*std::log(static_cast<DataType>(2.)))/GlobalVolume; 
00375     VectorType sc(1.);
00376     if (RelVal) {
00377       VectorType max=MaxVal(u,Time,Level);
00378       VectorType min=MinVal(u,Time,Level);
00379       sc = Max(fabs(max),fabs(min));
00380     }
00381     
00382     int TStep = TimeStep(work,Level);
00383     forall(u,Time,Level,c)
00384       BBox bb = u(Time,Level,c).bbox();
00385       BBox bbi = u.interiorbbox(Time,Level,c);
00386       BBox avbb = coarsen(bb,2);
00387       Coords avss = avbb.stepsize();
00388       vec_grid_data_type average(avbb);
00389 
00390       if (dim == 1) {
00391         BeginFastIndex1(u, bb, u(Time,Level,c).data(), VectorType);
00392         BeginFastIndex1(workn, bb, work(Time+TStep,Level,c).data(), DataType);
00393         BeginFastIndex1(average, avbb, average.data(), VectorType);
00394         int dn=avss(0);  
00395         for_1 (n,avbb,avss)
00396           FastIndex1(average,n)=0.;
00397           BBox bbadd(1,n,n,dn);
00398           bbadd.refine(2); bbadd*=bb;
00399           for_1 (nn,bbadd,bbadd.stepsize())
00400             FastIndex1(average,n)+=FastIndex1(u,nn);
00401           end_for
00402           FastIndex1(average,n)/=bbadd.size();
00403         end_for
00404 
00405         for_1 (nn,bbi,bbi.stepsize())
00406           int n=(nn/dn)*dn;     
00407           VectorType pred = FastIndex1(average,n);
00408 
00409           DataType pn = (n==nn) ? static_cast<DataType>(1.): static_cast<DataType>(-1.);
00410           pred += (pn*-static_cast<DataType>(.125)) * FastIndex1(average,n+dn);
00411           pred -= (pn*-static_cast<DataType>(.125)) * FastIndex1(average,n-dn);
00412 
00413           FastIndex1(workn,nn) = maxcomp(fabs(pred-FastIndex1(u,nn))/sc);
00414         end_for
00415         EndFastIndex1(u); 
00416         EndFastIndex1(workn); 
00417         EndFastIndex1(average); 
00418       }
00419       
00420       else if (dim == 2) {
00421         BeginFastIndex2(u, bb, u(Time,Level,c).data(), VectorType);
00422         BeginFastIndex2(workn, bb, work(Time+TStep,Level,c).data(), DataType);
00423         BeginFastIndex2(average, avbb, average.data(), VectorType);
00424         int dn=avss(0), dm=avss(1);  
00425         for_2 (n,m,avbb,avss)
00426           FastIndex2(average,n,m)=0.;
00427           BBox bbadd(2,n,m,n,m,dn,dm);
00428           bbadd.refine(2); bbadd*=bb;
00429           for_2 (nn,mm,bbadd,bbadd.stepsize())
00430             FastIndex2(average,n,m)+=FastIndex2(u,nn,mm);
00431           end_for
00432           FastIndex2(average,n,m)/=bbadd.size();
00433         end_for
00434 
00435         if (RelVal && MaxMomentum) { sc(1) = std::max(sc(1),sc(2)); sc(2) = sc(1); }
00436         for_2 (nn,mm,bbi,bbi.stepsize())
00437           int n=(nn/dn)*dn;
00438           int m=(mm/dm)*dm;
00439         
00440           VectorType pred = FastIndex2(average,n,m);
00441           
00442           DataType pn = (n==nn) ? static_cast<DataType>(1.): static_cast<DataType>(-1.);
00443           pred += (pn*-static_cast<DataType>(.125)) * FastIndex2(average,n+dn,m);
00444           pred -= (pn*-static_cast<DataType>(.125)) * FastIndex2(average,n-dn,m);
00445           
00446           DataType pm = (m==mm) ? static_cast<DataType>(1.): static_cast<DataType>(-1.);        
00447           pred += (pm*-static_cast<DataType>(.125)) * FastIndex2(average,n,m+dm); 
00448           pred -= (pm*-static_cast<DataType>(.125)) * FastIndex2(average,n,m-dm);
00449 
00450           pred += (pn*pm*static_cast<DataType>(.015625)) * FastIndex2(average,n+dn,m+dm);
00451           pred -= (pn*pm*static_cast<DataType>(.015625)) * FastIndex2(average,n+dn,m-dm);
00452           pred -= (pn*pm*static_cast<DataType>(.015625)) * FastIndex2(average,n-dn,m+dm);
00453           pred += (pn*pm*static_cast<DataType>(.015625)) * FastIndex2(average,n-dn,m-dm);
00454 
00455           pred = fabs(pred-FastIndex2(u,nn,mm));
00456           if (RelVal && MaxMomentum)
00457             { pred(1) = sqrt(pred(1)*pred(1)+pred(2)*pred(2)); pred(2) = pred(1); }       
00458           FastIndex2(workn,nn,mm) = maxcomp(pred/sc);
00459         end_for
00460         EndFastIndex2(u); 
00461         EndFastIndex2(workn); 
00462         EndFastIndex2(average); 
00463       }
00464       
00465       else if (dim == 3) {
00466         BeginFastIndex3(u, bb, u(Time,Level,c).data(), VectorType);
00467         BeginFastIndex3(workn, bb, work(Time+TStep,Level,c).data(), DataType);
00468         BeginFastIndex3(average, avbb, average.data(), VectorType);
00469         int dn=avss(0), dm=avss(1), dl=avss(2);  
00470         for_3 (n,m,l,avbb,avss)
00471           FastIndex3(average,n,m,l)=0.;
00472           BBox bbadd(3,n,m,l,n,m,l,dn,dm,dl);
00473           bbadd.refine(2); bbadd*=bb;
00474           for_3 (nn,mm,ll,bbadd,bbadd.stepsize())
00475             FastIndex3(average,n,m,l)+=FastIndex3(u,nn,mm,ll);
00476           end_for
00477           FastIndex3(average,n,m,l)/=bbadd.size();
00478         end_for
00479 
00480         if (RelVal && MaxMomentum)
00481           { sc(1) = std::max(sc(1),sc(2)); sc(1) = std::max(sc(1),sc(3)); sc(2) = sc(1); sc(3) = sc(1); }
00482         for_3 (nn,mm,ll,bbi,bbi.stepsize())
00483           int n=(nn/avss(0))*dn;
00484           int m=(mm/avss(1))*dm;
00485           int l=(mm/avss(2))*dl;
00486           
00487           VectorType pred = FastIndex3(average,n,m,l);
00488           
00489           DataType pn = (n==nn) ? static_cast<DataType>(1.): static_cast<DataType>(-1.);
00490           pred += (pn*-static_cast<DataType>(.125)) * FastIndex3(average,n+dn,m,l);
00491           pred -= (pn*-static_cast<DataType>(.125)) * FastIndex3(average,n-dn,m,l);
00492           
00493           DataType pm = (m==mm) ? static_cast<DataType>(1.): static_cast<DataType>(-1.);        
00494           pred += (pm*-static_cast<DataType>(.125)) * FastIndex3(average,n,m+dm,l); 
00495           pred -= (pm*-static_cast<DataType>(.125)) * FastIndex3(average,n,m-dm,l);
00496           
00497           DataType pl = (l==ll) ? static_cast<DataType>(1.): static_cast<DataType>(-1.);        
00498           pred += (pl*-static_cast<DataType>(.125)) * FastIndex3(average,n,m,l+dl); 
00499           pred -= (pl*-static_cast<DataType>(.125)) * FastIndex3(average,n,m,l-dl);
00500           
00501           pred += (pn*pm*static_cast<DataType>(.015625)) * FastIndex3(average,n+dn,m+dm,l);
00502           pred -= (pn*pm*static_cast<DataType>(.015625)) * FastIndex3(average,n+dn,m-dm,l);
00503           pred -= (pn*pm*static_cast<DataType>(.015625)) * FastIndex3(average,n-dn,m+dm,l);
00504           pred += (pn*pm*static_cast<DataType>(.015625)) * FastIndex3(average,n-dn,m-dm,l);
00505           
00506           pred += (pn*pl*static_cast<DataType>(.015625)) * FastIndex3(average,n+dn,m,l+dl);
00507           pred -= (pn*pl*static_cast<DataType>(.015625)) * FastIndex3(average,n+dn,m,l-dl);
00508           pred -= (pn*pl*static_cast<DataType>(.015625)) * FastIndex3(average,n-dn,m,l+dl);
00509           pred += (pn*pl*static_cast<DataType>(.015625)) * FastIndex3(average,n-dn,m,l-dl);
00510 
00511           pred += (pm*pl*static_cast<DataType>(.015625)) * FastIndex3(average,n,m+dm,l+dl);
00512           pred -= (pm*pl*static_cast<DataType>(.015625)) * FastIndex3(average,n,m+dm,l-dl);
00513           pred -= (pm*pl*static_cast<DataType>(.015625)) * FastIndex3(average,n,m-dm,l+dl);
00514           pred += (pm*pl*static_cast<DataType>(.015625)) * FastIndex3(average,n,m-dm,l-dl);
00515 
00516           pred += (pn*pm*pl-static_cast<DataType>(.001953125)) * FastIndex3(average,n+dn,m+dm,l+dl); 
00517           pred -= (pn*pm*pl-static_cast<DataType>(.001953125)) * FastIndex3(average,n+dn,m+dm,l-dl);
00518           pred -= (pn*pm*pl-static_cast<DataType>(.001953125)) * FastIndex3(average,n+dn,m-dm,l+dl);
00519           pred += (pn*pm*pl-static_cast<DataType>(.001953125)) * FastIndex3(average,n+dn,m-dm,l-dl);
00520           pred -= (pn*pm*pl-static_cast<DataType>(.001953125)) * FastIndex3(average,n-dn,m+dm,l+dl);
00521           pred += (pn*pm*pl-static_cast<DataType>(.001953125)) * FastIndex3(average,n-dn,m+dm,l-dl);
00522           pred += (pn*pm*pl-static_cast<DataType>(.001953125)) * FastIndex3(average,n-dn,m-dm,l+dl);
00523           pred -= (pn*pm*pl-static_cast<DataType>(.001953125)) * FastIndex3(average,n-dn,m-dm,l-dl);
00524 
00525           pred = fabs(pred-FastIndex2(u,nn,mm));
00526           if (RelVal && MaxMomentum)
00527             { pred(1) = sqrt(pred(1)*pred(1)+pred(2)*pred(2)+pred(3)*pred(3)); pred(2) = pred(1); pred(3) = pred(1); }    
00528           FastIndex3(workn,nn,mm,ll) = maxcomp(pred/sc);
00529         end_for
00530         EndFastIndex3(u); 
00531         EndFastIndex3(workn); 
00532         EndFastIndex3(average); 
00533       }
00534     end_forall 
00535       
00536     base::FlagByValue(work, flags, Time, Level, Tol, FlagValue);
00537     return true;
00538   }
00539 
00540   virtual void OutputName(char* name, int cnt) { std::sprintf(name,"mrvecpr_%d",cnt+1); }
00541 
00542   max_vector_type Tol;
00543 
00544   int ThresholdNorm[MAXCOMPONENTS], Relative[MAXCOMPONENTS];
00545 
00546 protected:
00547   int MaxMomentum;
00548   double GlobalVolume;
00549 };
00550 
00551 #endif
Generated on Thu Jun 30 02:19:15 2016 for AMROC Fluid-solver Framework by 1.6.3