amroc/amr/Interpolation.h

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 
00003 // Copyright (C) 2002 Ralf Deiterding
00004 // Brandenburgische Universitaet Cottbus
00005 //
00006 // Copyright (C) 2003-2007 California Institute of Technology
00007 // Ralf Deiterding, ralf@cacr.caltech.edu
00008 
00009 #ifndef AMROC_INTERPOLATION_H
00010 #define AMROC_INTERPOLATION_H
00011 
00019 #include "numerical/grid_interp_extrap.h"
00020 
00027 template <class VectorType, int dim>
00028 class Interpolation {
00029   typedef typename VectorType::InternalDataType DataType;
00030 public:
00031   typedef GridData<VectorType,dim> vec_grid_data_type;  
00032   typedef GridData<DataType,dim> grid_data_type;  
00033   typedef Vector<DataType,dim> point_type;
00034   typedef GridData<bool,dim> bool_grid_data_type;
00035 
00036   Interpolation() : _Equations(VectorType::Length()) {}
00037 
00038   virtual ~Interpolation() {}
00039 
00040   virtual void Interpolate(GridHierarchy& GH, vec_grid_data_type& gdu, const int& nc, 
00041                            const point_type* xc, VectorType* uv, const DataType& ErrorValue) {
00042 
00043     grid_data_type phihelp(grow(gdu.bbox(),2));
00044     phihelp.equals(static_cast<DataType>(-1.0));
00045     phihelp.equals(static_cast<DataType>(1.0),gdu.bbox());
00046 
00047     Coords ex = phihelp.extents();
00048     DCoords dxh2(GH.worldStep(phihelp.stepsize())); dxh2 *= static_cast<DataType>(0.5);
00049     DCoords lbcorner = GH.worldCoords(phihelp.lower(), phihelp.stepsize()) + dxh2;
00050     DCoords ubcorner = GH.worldCoords(phihelp.upper(), phihelp.stepsize()) + dxh2;
00051       
00052     DataType domain[2*dim];
00053     for (register int n=0; n<dim; n++) {
00054       domain[n] = lbcorner(n); domain[dim+n] = ubcorner(n);
00055     }
00056 
00057     register int cnt;
00058     VectorType default_values = ErrorValue;
00059     vec_grid_data_type Field(grow(gdu.bbox(),2));
00060     VectorType *values = new VectorType[nc];
00061     for (cnt=0; cnt<nc; cnt++) values[cnt] = ErrorValue;
00062     Field.copy(gdu,gdu.bbox());
00063      
00064     numerical::grid_interp_extrap<dim,VectorType::Size,DataType>( nc, reinterpret_cast<DataType *>(values), xc[0].data(), 
00065                                                                   reinterpret_cast<DataType *>(&default_values), ex(), domain, 
00066                                                                   phihelp.data(), reinterpret_cast<const DataType *>(Field.data()) );   
00067     for (cnt = 0; cnt<nc; cnt++) 
00068       if (values[cnt](0) != ErrorValue)   
00069         uv[cnt] = values[cnt];
00070 
00071     delete [] values;
00072   }
00073 
00074   virtual void Interpolate(GridHierarchy& GH, vec_grid_data_type& gdu, 
00075                            grid_data_type& gdphi, const int& nc, const point_type* xc, 
00076                            VectorType* uv, const DataType& ErrorValue) {
00077 
00078     grid_data_type phihelp(grow(gdphi.bbox(),2));
00079     phihelp.equals(static_cast<DataType>(-1.0));
00080     phihelp.equals(gdphi,gdphi.bbox());
00081 
00082     Coords ex = phihelp.extents();
00083     DCoords dxh2(GH.worldStep(phihelp.stepsize())); dxh2 *= static_cast<DataType>(0.5);
00084     DCoords lbcorner = GH.worldCoords(phihelp.lower(), phihelp.stepsize()) + dxh2;
00085     DCoords ubcorner = GH.worldCoords(phihelp.upper(), phihelp.stepsize()) + dxh2;
00086       
00087     DataType domain[2*dim];
00088     for (register int n=0; n<dim; n++) {
00089       domain[n] = lbcorner(n); domain[dim+n] = ubcorner(n);
00090     }
00091 
00092     register int cnt;
00093     VectorType default_values = ErrorValue;
00094     vec_grid_data_type Field(grow(gdu.bbox(),2));
00095     VectorType *values = new VectorType[nc];
00096     for (cnt=0; cnt<nc; cnt++) values[cnt] = ErrorValue;
00097     Field.copy(gdu,gdu.bbox());
00098      
00099     numerical::grid_interp_extrap<dim,VectorType::Size,DataType>( nc, reinterpret_cast<DataType *>(values), xc[0].data(), 
00100                                                                   reinterpret_cast<DataType *>(&default_values), ex(), domain, 
00101                                                                   phihelp.data(), reinterpret_cast<const DataType *>(Field.data()) );   
00102     for (cnt = 0; cnt<nc; cnt++) 
00103       if (values[cnt](0) != ErrorValue)   
00104         uv[cnt] = values[cnt];
00105 
00106     delete [] values;
00107   }
00108 
00109   virtual void Interpolate(GridHierarchy& GH, vec_grid_data_type& gdu, 
00110                            bool_grid_data_type& gdflg, const int& nc, const point_type* xc, 
00111                            VectorType* uv, const DataType& ErrorValue) {
00112 
00113     grid_data_type phihelp(grow(gdflg.bbox(),2));
00114     SetPhiFromBool(phihelp,gdflg);
00115 
00116     Coords ex = phihelp.extents();
00117     DCoords dxh2(GH.worldStep(phihelp.stepsize())); dxh2 *= static_cast<DataType>(0.5);
00118     DCoords lbcorner = GH.worldCoords(phihelp.lower(), phihelp.stepsize()) + dxh2;
00119     DCoords ubcorner = GH.worldCoords(phihelp.upper(), phihelp.stepsize()) + dxh2;
00120       
00121     DataType domain[2*dim];
00122     for (register int n=0; n<dim; n++) {
00123       domain[n] = lbcorner(n); domain[dim+n] = ubcorner(n);
00124     }
00125 
00126     register int cnt;
00127     VectorType default_values = ErrorValue;
00128     vec_grid_data_type Field(grow(gdu.bbox(),2));
00129     VectorType *values = new VectorType[nc];
00130     for (cnt=0; cnt<nc; cnt++) values[cnt] = ErrorValue;
00131     Field.copy(gdu,gdu.bbox());
00132      
00133     numerical::grid_interp_extrap<dim,VectorType::Size,DataType>( nc, reinterpret_cast<DataType *>(values), xc[0].data(), 
00134                                                                   reinterpret_cast<DataType *>(&default_values), ex(), domain, 
00135                                                                   phihelp.data(), reinterpret_cast<const DataType *>(Field.data()) );   
00136     for (cnt = 0; cnt<nc; cnt++) 
00137       if (values[cnt](0) != ErrorValue)   
00138         uv[cnt] = values[cnt];
00139 
00140     delete [] values;
00141   }
00142 
00143   virtual void Interpolate(GridHierarchy& GH, grid_data_type& gd, const int& nc, 
00144                            const point_type* xc, DataType* u, const DataType& ErrorValue) {
00145 
00146     grid_data_type phihelp(grow(gd.bbox(),2));
00147     phihelp.equals(static_cast<DataType>(-1.0));
00148     phihelp.equals(static_cast<DataType>(1.0),gd.bbox());
00149 
00150     Coords ex = phihelp.extents();
00151     DCoords dxh2(GH.worldStep(phihelp.stepsize())); dxh2 *= static_cast<DataType>(0.5);
00152     DCoords lbcorner = GH.worldCoords(phihelp.lower(), phihelp.stepsize()) + dxh2;
00153     DCoords ubcorner = GH.worldCoords(phihelp.upper(), phihelp.stepsize()) + dxh2;
00154       
00155     DataType domain[2*dim];
00156     for (register int n=0; n<dim; n++) {
00157       domain[n] = lbcorner(n); domain[dim+n] = ubcorner(n);
00158     }
00159 
00160     const DataType default_values[1] = { ErrorValue };
00161     grid_data_type field(grow(gd.bbox(),2));
00162     field.equals(gd,gd.bbox());
00163     const DataType* fields[1] = { field.data() };
00164 
00165     numerical::grid_interp_extrap<dim,1,DataType>( nc, u, xc[0].data(), 
00166                                                    default_values, ex(), domain, 
00167                                                    phihelp.data(), fields );
00168   }
00169 
00170   virtual void Interpolate(GridHierarchy& GH, grid_data_type& gd, 
00171                            grid_data_type& gdphi, const int& nc, const point_type* xc, 
00172                            DataType* u, const DataType& ErrorValue) {
00173 
00174     grid_data_type phihelp(grow(gdphi.bbox(),2));
00175     phihelp.equals(static_cast<DataType>(-1.0));
00176     phihelp.equals(gdphi,gdphi.bbox());
00177 
00178     Coords ex = phihelp.extents();
00179     DCoords dxh2(GH.worldStep(phihelp.stepsize())); dxh2 *= static_cast<DataType>(0.5);
00180     DCoords lbcorner = GH.worldCoords(phihelp.lower(), phihelp.stepsize()) + dxh2;
00181     DCoords ubcorner = GH.worldCoords(phihelp.upper(), phihelp.stepsize()) + dxh2;
00182       
00183     DataType domain[2*dim];
00184     for (register int n=0; n<dim; n++) {
00185       domain[n] = lbcorner(n); domain[dim+n] = ubcorner(n);
00186     }
00187 
00188     const DataType default_values[1] = { ErrorValue };
00189     grid_data_type field(grow(gd.bbox(),2));
00190     field.equals(gd,gd.bbox());
00191     const DataType* fields[1] = { field.data() };
00192 
00193     numerical::grid_interp_extrap<dim,1,DataType>( nc, u, xc[0].data(), 
00194                                                    default_values, ex(), domain, 
00195                                                    phihelp.data(), fields );
00196   }
00197 
00198   virtual void Interpolate(GridHierarchy& GH, grid_data_type& gd, 
00199                            bool_grid_data_type& gdflg, const int& nc, const point_type* xc, 
00200                            DataType* u, const DataType& ErrorValue) {
00201 
00202     grid_data_type phihelp(grow(gdflg.bbox(),2));
00203     SetPhiFromBool(phihelp,gdflg);
00204 
00205     Coords ex = phihelp.extents();
00206     DCoords dxh2(GH.worldStep(phihelp.stepsize())); dxh2 *= static_cast<DataType>(0.5);
00207     DCoords lbcorner = GH.worldCoords(phihelp.lower(), phihelp.stepsize()) + dxh2;
00208     DCoords ubcorner = GH.worldCoords(phihelp.upper(), phihelp.stepsize()) + dxh2;
00209       
00210     DataType domain[2*dim];
00211     for (register int n=0; n<dim; n++) {
00212       domain[n] = lbcorner(n); domain[dim+n] = ubcorner(n);
00213     }
00214 
00215     const DataType default_values[1] = { ErrorValue };
00216     grid_data_type field(grow(gd.bbox(),2));
00217     field.equals(gd,gd.bbox());
00218     const DataType* fields[1] = { field.data() };
00219 
00220     numerical::grid_interp_extrap<dim,1,DataType>( nc, u, xc[0].data(), 
00221                                                    default_values, ex(), domain, 
00222                                                    phihelp.data(), fields );
00223   }
00224 
00225   const int& NEquations() const { return _Equations; }
00226 
00227 protected:
00228   void SetPhiFromBool(grid_data_type& gdphi, bool_grid_data_type& gdflg) {
00229     gdphi.equals(static_cast<DataType>(-1.0));
00230     if (dim == 1) {
00231       BeginFastIndex1(phi, gdphi.bbox(), gdphi.data(), DataType);
00232       BeginFastIndex1(flg, gdflg.bbox(), gdflg.data(), bool);
00233       for_1 (n, gdflg.bbox(), gdflg.bbox().stepsize())
00234         if (!FastIndex1(flg,n)) 
00235           FastIndex1(phi,n) = static_cast<DataType>(1.0);
00236       end_for
00237       EndFastIndex1(phi);
00238       EndFastIndex1(flg);
00239     }      
00240     else if (dim == 2) {
00241       BeginFastIndex2(phi, gdphi.bbox(), gdphi.data(), DataType);
00242       BeginFastIndex2(flg, gdflg.bbox(), gdflg.data(), bool);
00243       for_2 (n, m, gdflg.bbox(), gdflg.bbox().stepsize())
00244         if (!FastIndex2(flg,n,m)) 
00245           FastIndex2(phi,n,m) = static_cast<DataType>(1.0);
00246       end_for
00247       EndFastIndex2(phi);
00248       EndFastIndex2(flg);
00249     }
00250     else if (dim == 3) {
00251       BeginFastIndex3(phi, gdphi.bbox(), gdphi.data(), DataType);
00252       BeginFastIndex3(flg, gdflg.bbox(), gdflg.data(), bool);
00253       for_3 (n, m, l, gdflg.bbox(), gdflg.bbox().stepsize())
00254         if (!FastIndex3(flg,n,m,l)) 
00255           FastIndex3(phi,n,m,l) = static_cast<DataType>(1.0);
00256       end_for
00257       EndFastIndex3(phi);
00258       EndFastIndex3(flg);
00259     }
00260   }
00261     
00262 protected:
00263   int _Equations;
00264 };
00265 
00266 #endif
Generated on Thu Jun 30 02:19:27 2016 for AMROC Fluid-solver Framework by 1.6.3