finiteDifference3D.h

00001 /* This file is part of the Palabos library.
00002  *
00003  * Copyright (C) 2011 FlowKit Sarl
00004  * Avenue de Chailly 23
00005  * 1012 Lausanne, Switzerland
00006  * E-mail contact: contact@flowkit.com
00007  *
00008  * The most recent release of Palabos can be downloaded at 
00009  * <http://www.palabos.org/>
00010  *
00011  * The library Palabos is free software: you can redistribute it and/or
00012  * modify it under the terms of the GNU Affero General Public License as
00013  * published by the Free Software Foundation, either version 3 of the
00014  * License, or (at your option) any later version.
00015  *
00016  * The library is distributed in the hope that it will be useful,
00017  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00018  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00019  * GNU Affero General Public License for more details.
00020  *
00021  * You should have received a copy of the GNU Affero General Public License
00022  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
00023 */
00024 
00025 #ifndef FINITE_DIFFERENCE_3D_H
00026 #define FINITE_DIFFERENCE_3D_H
00027 
00028 #include "core/globalDefs.h"
00029 #include "finiteDifference/fdStencils1D.h"
00030 
00031 namespace plb {
00032 
00033 namespace fd {
00034 
00035     template<typename T, template<typename U> class Descriptor,
00036              int direction, int orientation, int deriveDirection,
00037              bool orthogonal>
00038     struct DirectedGradients3D {
00041         static void o1_velocityDerivative( Array<T,Descriptor<T>::d>& velDeriv,
00042                                            BlockLattice3D<T,Descriptor> const& blockLattice,
00043                                            plint iX, plint iY, plint iZ );
00046         static void o1_densityDerivative( T& rhoDeriv,
00047                                           BlockLattice3D<T,Descriptor> const& blockLattice,
00048                                           plint iX, plint iY, plint iZ );
00051         static void o2_velocityDerivative( Array<T,Descriptor<T>::d>& velDeriv,
00052                                            BlockLattice3D<T,Descriptor> const& blockLattice,
00053                                            plint iX, plint iY, plint iZ );
00056         static void o2_densityDerivative( T& rhoDeriv,
00057                                           BlockLattice3D<T,Descriptor> const& blockLattice,
00058                                           plint iX, plint iY, plint iZ );
00059     };
00060 
00061     // Implementation for orthogonal==true; i.e. the derivative is along the boundary normal.
00062     template<typename T, template<typename U> class Descriptor,
00063              int direction, int orientation, int deriveDirection>
00064     struct DirectedGradients3D<T, Descriptor, direction, orientation,
00065                                deriveDirection, true>
00066     {
00067         static void o1_velocityDerivative(Array<T,Descriptor<T>::d>& velDeriv,
00068                                           BlockLattice3D<T,Descriptor> const& blockLattice,
00069                                           plint iX, plint iY, plint iZ)
00070         {
00071             Array<T,Descriptor<T>::d> u0, u1;
00072             
00073             blockLattice.get(iX,iY,iZ).computeVelocity(u0);
00074             blockLattice.get (
00075                 iX+(direction==0 ? (-orientation):0),
00076                 iY+(direction==1 ? (-orientation):0),
00077                 iZ+(direction==2 ? (-orientation):0)  ).computeVelocity(u1);
00078 
00079             for (int iD=0; iD<Descriptor<T>::d; ++iD) {
00080                 velDeriv[iD] = -orientation * fd::o1_fwd_diff(u0[iD], u1[iD]);
00081             }
00082         }
00083 
00084         static void o1_densityDerivative(T& rhoDeriv,
00085                                          BlockLattice3D<T,Descriptor> const& blockLattice,
00086                                          plint iX, plint iY, plint iZ)
00087         {
00088             // note that the derivative runs along direction.
00089             T rho0 = blockLattice.get(iX,iY,iZ).computeRho();
00090             T rho1 = blockLattice.get (
00091                         iX+(direction==0 ? (-orientation):0),
00092                         iY+(direction==1 ? (-orientation):0),
00093                         iZ+(direction==2 ? (-orientation):0)  ).computeRho();
00094             rhoDeriv = -orientation * fd::o1_fwd_diff(rho0, rho1);
00095         }
00096 
00097         static void o2_velocityDerivative(Array<T,Descriptor<T>::d>& velDeriv,
00098                                           BlockLattice3D<T,Descriptor> const& blockLattice,
00099                                           plint iX, plint iY, plint iZ)
00100         {
00101             Array<T,Descriptor<T>::d> u0, u1, u2;
00102             
00103             blockLattice.get(iX,iY,iZ).computeVelocity(u0);
00104             blockLattice.get (
00105                 iX+(direction==0 ? (-orientation):0),
00106                 iY+(direction==1 ? (-orientation):0),
00107                 iZ+(direction==2 ? (-orientation):0)  ).computeVelocity(u1);
00108             blockLattice.get (
00109                 iX+(direction==0 ? (-2*orientation):0),
00110                 iY+(direction==1 ? (-2*orientation):0),
00111                 iZ+(direction==2 ? (-2*orientation):0) ).computeVelocity(u2);
00112 
00113             for (int iD=0; iD<Descriptor<T>::d; ++iD) {
00114                 velDeriv[iD] = -orientation * fd::fwd_diff(u0[iD], u1[iD], u2[iD]);
00115             }
00116         }
00117 
00118         static void o2_densityDerivative(T& rhoDeriv,
00119                                          BlockLattice3D<T,Descriptor> const& blockLattice,
00120                                          plint iX, plint iY, plint iZ)
00121         {
00122             // note that the derivative runs along direction.
00123             T rho0 = blockLattice.get(iX,iY,iZ).computeRho();
00124             T rho1 = blockLattice.get (
00125                         iX+(direction==0 ? (-orientation):0),
00126                         iY+(direction==1 ? (-orientation):0),
00127                         iZ+(direction==2 ? (-orientation):0)  ).computeRho();
00128             T rho2 = blockLattice.get (
00129                         iX+(direction==0 ? (-2*orientation):0),
00130                         iY+(direction==1 ? (-2*orientation):0),
00131                         iZ+(direction==2 ? (-2*orientation):0) ).computeRho();
00132 
00133             rhoDeriv = -orientation * fd::fwd_diff(rho0, rho1, rho2);
00134         }
00135     };
00136 
00137     // Implementation for orthogonal==false; i.e. the derivative is aligned with the boundary.
00138     template<typename T, template<typename U> class Descriptor,
00139              int direction, int orientation, int deriveDirection>
00140     struct DirectedGradients3D<T, Descriptor, direction, orientation,
00141                                deriveDirection, false>
00142     {
00143         static void  o1_velocityDerivative(Array<T,Descriptor<T>::d>& velDeriv,
00144                                            BlockLattice3D<T,Descriptor> const& blockLattice,
00145                                            plint iX, plint iY, plint iZ)
00146         {
00147             // Along the boundary, second-order accuracy is achieved with a nearest-
00148             //   neighbor scheme.
00149             o2_velocityDerivative(velDeriv, blockLattice, iX,iY,iZ);
00150         }
00151         static void  o1_densityDerivative(T& rhoDeriv,
00152                                           BlockLattice3D<T,Descriptor> const& blockLattice,
00153                                           plint iX, plint iY, plint iZ)
00154         {
00155             // Along the boundary, second-order accuracy is achieved with a nearest-
00156             //   neighbor scheme.
00157             o2_densityDerivative(rhoDeriv, blockLattice, iX,iY,iZ);
00158         }
00159         static void  o2_velocityDerivative(Array<T,Descriptor<T>::d>& velDeriv,
00160                                            BlockLattice3D<T,Descriptor> const& blockLattice,
00161                                            plint iX, plint iY, plint iZ)
00162         {
00163             Array<T,Descriptor<T>::d> u_p1, u_m1;
00164             
00165             blockLattice.get (
00166                 iX+(deriveDirection==0 ? 1:0),
00167                 iY+(deriveDirection==1 ? 1:0),
00168                 iZ+(deriveDirection==2 ? 1:0) ).computeVelocity(u_p1);
00169 
00170             blockLattice.get (
00171                 iX+(deriveDirection==0 ? (-1):0),
00172                 iY+(deriveDirection==1 ? (-1):0),
00173                 iZ+(deriveDirection==2 ? (-1):0) ).computeVelocity(u_m1);
00174 
00175             for (int iD=0; iD<Descriptor<T>::d; ++iD) {
00176                 velDeriv[iD] = fd::ctl_diff(u_p1[iD],u_m1[iD]);
00177             }
00178         }
00179 
00180         static void  o2_densityDerivative(T& rhoDeriv,
00181                                           BlockLattice3D<T,Descriptor> const& blockLattice,
00182                                           plint iX, plint iY, plint iZ)
00183         {
00184             typedef Descriptor<T> L;
00185 
00186             T rho_p1 = blockLattice.get (
00187                         iX+(deriveDirection==0 ? 1:0),
00188                         iY+(deriveDirection==1 ? 1:0),
00189                         iZ+(deriveDirection==2 ? 1:0) ).computeRho();
00190 
00191             T rho_m1 = blockLattice.get (
00192                         iX+(deriveDirection==0 ? (-1):0),
00193                         iY+(deriveDirection==1 ? (-1):0),
00194                         iZ+(deriveDirection==2 ? (-1):0) ).computeRho();
00195 
00196             rhoDeriv = fd::ctl_diff(rho_p1, rho_m1);
00197         }
00198     };
00199 
00200 }  // namespace fd
00201 
00202 }  // namespace plb
00203 
00204 
00205 #endif  // FINITE_DIFFERENCE_3D_H