indexTemplates.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 
00029 #ifndef INDEX_TEMPLATES_H
00030 #define INDEX_TEMPLATES_H
00031 
00032 #include "core/globalDefs.h"
00033 #include "core/array.h"
00034 #include <algorithm>
00035 #include <vector>
00036 
00037 namespace plb {
00038 
00039 namespace indexTemplates {
00040 
00042 template <typename Descriptor> inline plint opposite(plint iPop) {
00043     if (iPop==0) return 0;
00044     if (iPop<=Descriptor::q/2) return iPop + Descriptor::q/2;
00045     return iPop - Descriptor::q/2;
00046 }
00047 
00049 template <typename Descriptor> inline std::vector<plint> opposite(const std::vector<plint> &pops_) {
00050     std::vector<plint> pops(pops_);
00051     for (pluint iPop = 0; iPop < pops_.size(); ++iPop) {
00052         pops[iPop] = opposite<Descriptor>(pops_[iPop]);
00053     }
00054     return pops;
00055 }
00056 
00057 template <typename Descriptor>
00058 plint findVelocity(Array<int,Descriptor::d> const& v) {
00059     for (plint iPop=0; iPop<Descriptor::q; ++iPop) {
00060         bool fit = true;
00061         for (int iD=0; iD<Descriptor::d; ++iD) {
00062             if (Descriptor::c[iPop][iD] != v[iD]) {
00063                 fit = false;
00064                 break;
00065             }
00066         }
00067         if (fit) return iPop;
00068     }
00069     return Descriptor::q;
00070 }
00071 
00073 template <typename Descriptor, int orientation> inline plint specularReflexion(plint iPop) {
00074     if (iPop==0) return 0;
00075     Array<int,Descriptor::d> v; 
00076     for (plint iD = 0; iD < Descriptor::d; ++iD) v[iD] = Descriptor::c[iPop][iD];
00077     v[orientation] = -v[orientation];
00078     
00079     return findVelocity<Descriptor>(v);
00080 }
00081 
00082 template <typename Descriptor, plint index, plint value>
00083 class SubIndex {
00084 private:
00085     SubIndex() {
00086         for (int iVel=0; iVel<Descriptor::q; ++iVel) {
00087             if (Descriptor::c[iVel][index]==value) {
00088                 indices.push_back(iVel);
00089             }
00090         }
00091     }
00092 
00093     std::vector<plint> indices;
00094 
00095     template <typename Descriptor_, plint index_, plint value_>
00096     friend std::vector<plint> const& subIndex();
00097 };
00098 
00099 template <typename Descriptor, plint index, plint value>
00100 std::vector<plint> const& subIndex() {
00101     static SubIndex<Descriptor, index, value> subIndexSingleton;
00102     return subIndexSingleton.indices;
00103 }
00104 
00110 template <typename Descriptor, int direction, int orientation>
00111 class SubIndexOutgoing {
00112 private:
00113     SubIndexOutgoing() // finds the indexes outgoing from the walls
00114     {
00115         indices = indexTemplates::subIndex<Descriptor,direction,orientation>();
00116 
00117         for (pluint iPop = 0; iPop < indices.size(); ++iPop) {
00118             indices[iPop] = indexTemplates::opposite<Descriptor>(indices[iPop]);
00119         }
00120 
00121     }
00122 
00123     std::vector<plint> indices;
00124 
00125     template <typename Descriptor_, int direction_, int orientation_>
00126     friend std::vector<plint> const& subIndexOutgoing();
00127 };
00128 
00129 template <typename Descriptor, int direction, int orientation>
00130 std::vector<plint> const& subIndexOutgoing() {
00131     static SubIndexOutgoing<Descriptor, direction, orientation> subIndexOutgoingSingleton;
00132     return subIndexOutgoingSingleton.indices;
00133 }
00134 
00136 template <typename Descriptor>
00137 std::vector<plint> remainingIndexes(const std::vector<plint> &indices)
00138 {
00139     std::vector<plint> remaining;
00140     for (plint iPop = 0; iPop < Descriptor::q; ++iPop)
00141     {
00142         bool found = false;
00143         for (pluint jPop = 0; jPop < indices.size(); ++jPop)
00144         {
00145             if (indices[jPop] == iPop)
00146             {
00147                 found = true;
00148             }
00149         }
00150         if (!found)
00151         {
00152             remaining.push_back(iPop);
00153         }
00154     }
00155     return remaining;
00156 }
00157 
00159 template <typename Descriptor, int xNormal, int yNormal>
00160 class SubIndexIngoingCorner2D {
00161 private:
00162     SubIndexIngoingCorner2D()
00163     {
00164         typedef Descriptor L;
00165 
00166         Array<int, L::d> vect ( -xNormal, -yNormal );
00167         indices.push_back(indexTemplates::findVelocity<L>(vect));
00168         vect[0] = -xNormal; vect[1] = 0;
00169         indices.push_back(indexTemplates::findVelocity<L>(vect));
00170         vect[0] = 0; vect[1] = -yNormal;
00171         indices.push_back(indexTemplates::findVelocity<L>(vect));
00172     }
00173 
00174     std::vector<plint> indices;
00175 
00176     template <typename Descriptor_, int xNormal_, int yNormal_>
00177     friend std::vector<plint> const& subIndexIngoingCorner2D();
00178 };
00179 
00180 template <typename Descriptor, int xNormal, int yNormal>
00181 std::vector<plint> const& subIndexIngoingCorner2D() {
00182     static SubIndexIngoingCorner2D<Descriptor, xNormal, yNormal> subIndexIngoingCorner2DSingleton;
00183     return subIndexIngoingCorner2DSingleton.indices;
00184 }
00185 
00187 template <typename Descriptor, int xNormal, int yNormal>
00188 class SubIndexOutgoingCorner2D {
00189 private:
00190     SubIndexOutgoingCorner2D()
00191     {
00192         typedef Descriptor L;
00193 
00194         Array<int, L::d> vect ( xNormal, yNormal );
00195         std::vector<plint> knownIndexes;
00196         knownIndexes.push_back(indexTemplates::findVelocity<L>(vect));
00197         vect[0] = xNormal; vect[1] = 0;
00198         knownIndexes.push_back(indexTemplates::findVelocity<L>(vect));
00199         vect[0] = 0; vect[1] = yNormal;
00200         knownIndexes.push_back(indexTemplates::findVelocity<L>(vect));
00201         vect[0] = 0; vect[1] = 0;
00202         knownIndexes.push_back(indexTemplates::findVelocity<L>(vect));
00203         
00204         indices = indexTemplates::remainingIndexes<L>(knownIndexes);
00205     }
00206 
00207     std::vector<plint> indices;
00208 
00209     template <typename Descriptor_, int xNormal_, int yNormal_>
00210     friend std::vector<plint> const& subIndexOutgoingCorner2D();
00211 };
00212 
00213 template <typename Descriptor, int xNormal, int yNormal>
00214 std::vector<plint> const& subIndexOutgoingCorner2D() {
00215     static SubIndexOutgoingCorner2D<Descriptor, xNormal, yNormal> subIndexOutgoingCorner2DSingleton;
00216     return subIndexOutgoingCorner2DSingleton.indices;
00217 }
00218 
00220 template <typename Descriptor, int xNormal, int yNormal>
00221 class SubIndexOutgoingInternalCorner2D {
00222 private:
00223     SubIndexOutgoingInternalCorner2D()
00224     {
00225         typedef Descriptor L;
00226 
00227         Array<int, L::d> vect ( -xNormal, -yNormal );
00228         indices.push_back(indexTemplates::findVelocity<L>(vect));
00229     }
00230 
00231     std::vector<plint> indices;
00232 
00233     template <typename Descriptor_, int xNormal_, int yNormal_>
00234     friend std::vector<plint> const& subIndexOutgoingInternalCorner2D();
00235 };
00236 
00237 template <typename Descriptor, int xNormal, int yNormal>
00238 std::vector<plint> const& subIndexOutgoingInternalCorner2D() {
00239     static SubIndexOutgoingInternalCorner2D<Descriptor, xNormal, yNormal> subIndexOutgoingInternalCorner2DSingleton;
00240     return subIndexOutgoingInternalCorner2DSingleton.indices;
00241 }
00242 
00244 template <typename Descriptor, int xNormal, int yNormal, int zNormal>
00245 class SubIndexOutgoingExternalCorner3D {
00246 private:
00247     SubIndexOutgoingExternalCorner3D() {
00248         typedef Descriptor L;
00249         
00250         Array<int, L::d> vect ( xNormal, yNormal, zNormal);
00251         std::vector<plint> knownIndexes;
00252         knownIndexes.push_back(indexTemplates::findVelocity<L>(vect));
00253         
00254         vect[0] = xNormal; vect[1] = 0;       vect[2] = 0;
00255         knownIndexes.push_back(indexTemplates::findVelocity<L>(vect));
00256         vect[0] = 0;       vect[1] = yNormal; vect[2] = 0;
00257         knownIndexes.push_back(indexTemplates::findVelocity<L>(vect));
00258         vect[0] = 0;       vect[1] = 0;       vect[2] = zNormal;
00259         knownIndexes.push_back(indexTemplates::findVelocity<L>(vect));
00260         
00261         vect[0] = xNormal; vect[1] = yNormal; vect[2] = 0;
00262         knownIndexes.push_back(indexTemplates::findVelocity<L>(vect));
00263         vect[0] = 0;       vect[1] = yNormal; vect[2] = zNormal;
00264         knownIndexes.push_back(indexTemplates::findVelocity<L>(vect));
00265         vect[0] = xNormal; vect[1] = 0;       vect[2] = zNormal;
00266         knownIndexes.push_back(indexTemplates::findVelocity<L>(vect));
00267         
00268         vect[0] = 0;       vect[1] = 0;       vect[2] = 0;
00269         knownIndexes.push_back(indexTemplates::findVelocity<L>(vect));
00270         
00271         indices = indexTemplates::remainingIndexes<L>(knownIndexes);
00272     }
00273     
00274     std::vector<plint> indices;
00275     
00276     template <typename Descriptor_, int xNormal_, int yNormal_, int zNormal_>
00277     friend std::vector<plint> const& subIndexOutgoingExternalCorner3D();
00278 };
00279 
00280 template <typename Descriptor, int xNormal, int yNormal, int zNormal>
00281 std::vector<plint> const& subIndexOutgoingExternalCorner3D() {
00282     static SubIndexOutgoingExternalCorner3D<Descriptor, xNormal, yNormal,zNormal>
00283         subIndexOutgoingExternalCorner3DSingleton;
00284     return subIndexOutgoingExternalCorner3DSingleton.indices;
00285 }
00286 
00288 template <typename Descriptor, int plane, int normal1, int normal2>
00289 class SubIndexOutgoingExternalEdge3D {
00290 private:
00291     SubIndexOutgoingExternalEdge3D() {
00292         enum { direction1 = (plane+1)%3, direction2 = (plane+2)%3 };
00293         typedef Descriptor L;
00294         
00295         Array<int, L::d> vect(0, 0, 0);
00296         std::vector<plint> knownIndexes;
00297         knownIndexes.push_back(indexTemplates::findVelocity<L>(vect));
00298         
00299         vect[direction1] = 0; vect[direction2] = 0; vect[plane] = -1;
00300         knownIndexes.push_back(indexTemplates::findVelocity<L>(vect));
00301         vect[direction1] = 0; vect[direction2] = 0; vect[plane] = +1;
00302         knownIndexes.push_back(indexTemplates::findVelocity<L>(vect));
00303         
00304         vect[direction1] = normal1; vect[direction2] = 0; vect[plane] = -1;
00305         knownIndexes.push_back(indexTemplates::findVelocity<L>(vect));
00306         vect[direction1] = normal1; vect[direction2] = 0; vect[plane] = 0;
00307         knownIndexes.push_back(indexTemplates::findVelocity<L>(vect));
00308         vect[direction1] = normal1; vect[direction2] = 0; vect[plane] = +1;
00309         knownIndexes.push_back(indexTemplates::findVelocity<L>(vect));
00310         
00311         vect[direction1] = 0; vect[direction2] = normal2; vect[plane] = -1;
00312         knownIndexes.push_back(indexTemplates::findVelocity<L>(vect));
00313         vect[direction1] = 0; vect[direction2] = normal2; vect[plane] = 0;
00314         knownIndexes.push_back(indexTemplates::findVelocity<L>(vect));
00315         vect[direction1] = 0; vect[direction2] = normal2; vect[plane] = +1;
00316         knownIndexes.push_back(indexTemplates::findVelocity<L>(vect));
00317         
00318         vect[direction1] = normal1; vect[direction2] = normal2; vect[plane] = -1;
00319         knownIndexes.push_back(indexTemplates::findVelocity<L>(vect));
00320         vect[direction1] = normal1; vect[direction2] = normal2; vect[plane] = 0;
00321         knownIndexes.push_back(indexTemplates::findVelocity<L>(vect));
00322         vect[direction1] = normal1; vect[direction2] = normal2; vect[plane] = +1;
00323         knownIndexes.push_back(indexTemplates::findVelocity<L>(vect));
00324         
00325         indices = indexTemplates::remainingIndexes<L>(knownIndexes);
00326     }
00327     
00328     std::vector<plint> indices;
00329     
00330     template <typename Descriptor_, int plane_, int normal1_, int normal2_>
00331     friend std::vector<plint> const& subIndexOutgoingExternalEdge3D();
00332 };
00333 
00334 template <typename Descriptor, int plane, int normal1, int normal2>
00335 std::vector<plint> const& subIndexOutgoingExternalEdge3D() {
00336     static SubIndexOutgoingExternalEdge3D<Descriptor, plane, normal1, normal2>
00337         subIndexOutgoingExternalEdge3DSingleton;
00338     return subIndexOutgoingExternalEdge3DSingleton.indices;
00339 }
00340 
00342 template <typename Descriptor, int plane, int normal1, int normal2>
00343 class SubIndexOutgoingInternalEdge3D {
00344 private:
00345     SubIndexOutgoingInternalEdge3D() {
00346         enum { direction1 = (plane+1)%3, direction2 = (plane+2)%3 };
00347         
00348         typedef Descriptor L;
00349         
00350         Array<int, L::d> vect(0,0,0);
00351         vect[direction1] = -normal1;
00352         vect[direction2] = -normal2;
00353         vect[plane] = -1;
00354         indices.push_back(indexTemplates::findVelocity<L>(vect));
00355         vect[plane] = 0;
00356         indices.push_back(indexTemplates::findVelocity<L>(vect));
00357         vect[plane] = +1;
00358         indices.push_back(indexTemplates::findVelocity<L>(vect));
00359     }
00360     
00361     std::vector<plint> indices;
00362     
00363     template <typename Descriptor_, int plane_, int normal1_, int normal2_>
00364     friend std::vector<plint> const& subIndexOutgoingInternalEdge3D();
00365 };
00366 
00367 template <typename Descriptor, int plane, int normal1, int normal2>
00368 std::vector<plint> const& subIndexOutgoingInternalEdge3D() {
00369     static SubIndexOutgoingInternalEdge3D<Descriptor, plane, normal1, normal2>
00370         subIndexOutgoingInternalEdge3DSingleton;
00371     return subIndexOutgoingInternalEdge3DSingleton.indices;
00372 }
00373 
00375 template <typename Descriptor, int xNormal, int yNormal, int zNormal>
00376 class SubIndexOutgoingInternalCorner3D {
00377 private:
00378     SubIndexOutgoingInternalCorner3D()
00379     {
00380         typedef Descriptor L;
00381         
00382         Array<int, L::d> vect ( -xNormal, -yNormal, -zNormal );
00383         indices.push_back(indexTemplates::findVelocity<L>(vect));
00384     }
00385     
00386     std::vector<plint> indices;
00387     
00388     template <typename Descriptor_, int xNormal_, int yNormal_, int zNormal_>
00389     friend std::vector<plint> const& subIndexOutgoingInternalCorner3D();
00390 };
00391 
00392 template <typename Descriptor, int xNormal, int yNormal, int zNormal>
00393 std::vector<plint> const& subIndexOutgoingInternalCorner3D() {
00394     static SubIndexOutgoingInternalCorner3D<Descriptor, xNormal, yNormal, zNormal>
00395         subIndexOutgoingInternalCorner3DSingleton;
00396     return subIndexOutgoingInternalCorner3DSingleton.indices;
00397 }
00398 
00399 }  // namespace indexTemplates
00400 
00401 }  // namespace plb
00402 
00403 #endif  // INDEX_TEMPLATES_H