multiPhaseTemplates3D.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 
00031 #ifndef MULTI_PHASE_TEMPLATES_3D_H
00032 #define MULTI_PHASE_TEMPLATES_3D_H
00033 
00034 #include "core/globalDefs.h"
00035 #include "core/cell.h"
00036 #include "multiPhysics/shanChenLattices3D.h"
00037 
00038 namespace plb {
00039 
00041 template<typename T, template<typename U> class Descriptor>
00042 struct multiPhaseTemplates3D {
00043 
00044 static void shanChenInteraction( BlockLattice3D<T,Descriptor>& lattice,
00045                                  Array<T,Descriptor<T>::d>& rhoContribution,
00046                                  plint iX, plint iY, plint iZ )
00047 {
00048     enum {
00049         densityOffset  = Descriptor<T>::ExternalField::densityBeginsAt
00050     };
00051                      
00052     rhoContribution.resetToZero();
00053     for (plint iPop = 0; iPop < Descriptor<T>::q; ++iPop) {
00054         plint nextX = iX + Descriptor<T>::c[iPop][0];
00055         plint nextY = iY + Descriptor<T>::c[iPop][1];
00056         plint nextZ = iZ + Descriptor<T>::c[iPop][2];
00057         Cell<T,Descriptor> const& cell = lattice.get(nextX,nextY,nextZ);
00058         T rho = *cell.getExternal(densityOffset);
00059         for (int iD = 0; iD < Descriptor<T>::d; ++iD) {
00060            rhoContribution[iD] += Descriptor<T>::t[iPop] * rho * Descriptor<T>::c[iPop][iD];
00061         }
00062     }
00063 }
00064 
00065 };
00066 
00067 template<typename T>
00068 struct multiPhaseTemplates3D<T, descriptors::ForcedShanChenD3Q19Descriptor> {
00069     typedef descriptors::ForcedShanChenD3Q19Descriptor<T> D;
00070 
00071 static void shanChenInteraction (
00072         BlockLattice3D<T,descriptors::ForcedShanChenD3Q19Descriptor>& lattice,
00073         Array<T,D::d>& rhoContribution,
00074         plint iX, plint iY, plint iZ )
00075 {
00076 
00077     enum {
00078         densityOffset  = D::ExternalField::densityBeginsAt
00079     };
00080 
00081     T rho;
00082     rho = *lattice.get(iX-1,iY  ,iZ  ).getExternal(densityOffset);
00083     rhoContribution[0] = -D::t[1] * rho;
00084     rho = *lattice.get(iX  ,iY-1,iZ  ).getExternal(densityOffset);
00085     rhoContribution[1] = -D::t[2] * rho;
00086     rho = *lattice.get(iX  ,iY  ,iZ-1).getExternal(densityOffset);
00087     rhoContribution[2] = -D::t[3] * rho;
00088     rho = *lattice.get(iX-1,iY-1,iZ  ).getExternal(densityOffset);
00089     rhoContribution[0] -= D::t[4] * rho;
00090     rhoContribution[1] -= D::t[4] * rho;
00091     rho = *lattice.get(iX-1,iY+1,iZ  ).getExternal(densityOffset);
00092     rhoContribution[0] -= D::t[5] * rho;
00093     rhoContribution[1] += D::t[5] * rho;
00094     rho = *lattice.get(iX-1,iY  ,iZ-1).getExternal(densityOffset);
00095     rhoContribution[0] -= D::t[6] * rho;
00096     rhoContribution[2] -= D::t[6] * rho;
00097     rho = *lattice.get(iX-1,iY  ,iZ+1).getExternal(densityOffset);
00098     rhoContribution[0] -= D::t[7] * rho;
00099     rhoContribution[2] += D::t[7] * rho;
00100     rho = *lattice.get(iX  ,iY-1,iZ-1).getExternal(densityOffset);
00101     rhoContribution[1] -= D::t[8] * rho;
00102     rhoContribution[2] -= D::t[8] * rho;
00103     rho = *lattice.get(iX  ,iY-1,iZ+1).getExternal(densityOffset);
00104     rhoContribution[1] -= D::t[9] * rho;
00105     rhoContribution[2] += D::t[9] * rho;
00106 
00107     rho = *lattice.get(iX+1,iY  ,iZ  ).getExternal(densityOffset);
00108     rhoContribution[0] += D::t[10] * rho;
00109     rho = *lattice.get(iX  ,iY+1,iZ  ).getExternal(densityOffset);
00110     rhoContribution[1] += D::t[11] * rho;
00111     rho = *lattice.get(iX  ,iY  ,iZ+1).getExternal(densityOffset);
00112     rhoContribution[2] += D::t[12] * rho;
00113     rho = *lattice.get(iX+1,iY+1,iZ  ).getExternal(densityOffset);
00114     rhoContribution[0] += D::t[13] * rho;
00115     rhoContribution[1] += D::t[13] * rho;
00116     rho = *lattice.get(iX+1,iY-1,iZ  ).getExternal(densityOffset);
00117     rhoContribution[0] += D::t[14] * rho;
00118     rhoContribution[1] -= D::t[14] * rho;
00119     rho = *lattice.get(iX+1,iY  ,iZ+1).getExternal(densityOffset);
00120     rhoContribution[0] += D::t[15] * rho;
00121     rhoContribution[2] += D::t[15] * rho;
00122     rho = *lattice.get(iX+1,iY  ,iZ-1).getExternal(densityOffset);
00123     rhoContribution[0] += D::t[16] * rho;
00124     rhoContribution[2] -= D::t[16] * rho;
00125     rho = *lattice.get(iX  ,iY+1,iZ+1).getExternal(densityOffset);
00126     rhoContribution[1] += D::t[17] * rho;
00127     rhoContribution[2] += D::t[17] * rho;
00128     rho = *lattice.get(iX  ,iY+1,iZ-1).getExternal(densityOffset);
00129     rhoContribution[1] += D::t[18] * rho;
00130     rhoContribution[2] -= D::t[18] * rho;
00131 }
00132 
00133 };
00134 
00135 }  // namespace plb
00136 
00137 #endif  // MULTI_PHASE_TEMPLATES_3D_H