mrtTemplates3D.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 /* Orestis Malaspinas contributed this code.
00026  */
00027 
00033 #ifndef MRT_TEMPLATES_3D_H
00034 #define MRT_TEMPLATES_3D_H
00035 
00036 #include "core/globalDefs.h"
00037 
00038 namespace plb {
00039 
00040 // Efficient specialization for D3Q19 lattice
00041 template<typename T> struct mrtTemplatesImpl<T, descriptors::MRTD3Q19DescriptorBase<T> > {
00042     
00043     typedef descriptors::D3Q19DescriptorBase<T> Descriptor; 
00044 
00046     static void computeEquilibrium( Array<T,Descriptor::q>& momentsEq,
00047                                     T rhoBar, Array<T,3> const& j, T jSqr )
00048     {
00049         T invRho = Descriptor::invRho(rhoBar);
00050         
00051         momentsEq[0] = rhoBar;
00052         momentsEq[1] = (T)19*jSqr*invRho-(T)11*rhoBar;
00053         momentsEq[2] = -(T)5.5*jSqr*invRho+(T)3*rhoBar;
00054         momentsEq[3] = j[0];
00055         momentsEq[4] = -((T)2/(T)3)*j[0];
00056         momentsEq[5] = j[1];
00057         momentsEq[6] = -((T)2/3)*j[1];
00058         momentsEq[7] = j[2];
00059         momentsEq[8] = -((T)2/(T)3)*j[2];
00060         momentsEq[9] = ((T)2*j[0]*j[0]-j[1]*j[1]-j[2]*j[2])*invRho;
00061         momentsEq[10] = (-j[0]*j[0]+(T)0.5*j[1]*j[1]+(T)0.5*j[2]*j[2])*invRho;
00062         momentsEq[11] = (j[1]*j[1]-j[2]*j[2])*invRho;
00063         momentsEq[12] = (-(T)0.5*j[1]*j[1]+(T)0.5*j[2]*j[2])*invRho;
00064         momentsEq[13] = j[1]*j[0]*invRho;
00065         momentsEq[14] = j[2]*j[1]*invRho;
00066         momentsEq[15] = j[2]*j[0]*invRho;
00067         momentsEq[16] = T();
00068         momentsEq[17] = T();
00069         momentsEq[18] = T();
00070         
00071     }
00072     
00074     static void computeMoments(Array<T,Descriptor::q>& moments, const Array<T,Descriptor::q>& f)
00075     {
00076         moments[0] = f[0]+f[1]+f[2]+f[3]+f[4]+f[5]+f[6]+f[7]+f[8]+f[9]+f[10]+f[11]+f[12]+f[13]+f[14]+f[15]+f[16]+f[17]+f[18];
00077         moments[1] =
00078             -(T)30*f[0]-(T)11*(f[1]+f[2]+f[3]+f[10]+f[11]+f[12])
00079             +(T)8*(f[4]+f[5]+f[6]+f[7]+f[8]+f[9]+f[13]+f[14]+f[15]+f[16]+f[17]+f[18]);
00080         moments[2] = 12*f[0]-4*f[1]-4*f[2]-4*f[3]+f[4]+f[5]+f[6]+f[7]+f[8]+f[9]-4*f[10]-4*f[11]-4*f[12]+f[13]+f[14]+f[15]+f[16]+f[17]+f[18];
00081         moments[3] = -f[1]-f[4]-f[5]-f[6]-f[7]+f[10]+f[13]+f[14]+f[15]+f[16];
00082         moments[4] = 4*f[1]-f[4]-f[5]-f[6]-f[7]-4*f[10]+f[13]+f[14]+f[15]+f[16];
00083         moments[5] = -f[2]-f[4]+f[5]-f[8]-f[9]+f[11]+f[13]-f[14]+f[17]+f[18];
00084         moments[6] = 4*f[2]-f[4]+f[5]-f[8]-f[9]-4*f[11]+f[13]-f[14]+f[17]+f[18];
00085         moments[7] = -f[3]-f[6]+f[7]-f[8]+f[9]+f[12]+f[15]-f[16]+f[17]-f[18];
00086         moments[8] = 4*f[3]-f[6]+f[7]-f[8]+f[9]-4*f[12]+f[15]-f[16]+f[17]-f[18];
00087         moments[9] = 2*f[1]-f[2]-f[3]+f[4]+f[5]+f[6]+f[7]-2*f[8]-2*f[9]+2*f[10]-f[11]-f[12]+f[13]+f[14]+f[15]+f[16]-2*f[17]-2*f[18];
00088         moments[10] = -4*f[1]+2*f[2]+2*f[3]+f[4]+f[5]+f[6]+f[7]-2*f[8]-2*f[9]-4*f[10]+2*f[11]+2*f[12]+f[13]+f[14]+f[15]+f[16]-2*f[17]-2*f[18];
00089         moments[11] = f[2]-f[3]+f[4]+f[5]-f[6]-f[7]+f[11]-f[12]+f[13]+f[14]-f[15]-f[16];
00090         moments[12] = -2*f[2]+2*f[3]+f[4]+f[5]-f[6]-f[7]-2*f[11]+2*f[12]+f[13]+f[14]-f[15]-f[16];
00091         moments[13] = f[4]-f[5]+f[13]-f[14];
00092         moments[14] = f[8]-f[9]+f[17]-f[18];
00093         moments[15] = f[6]-f[7]+f[15]-f[16];
00094         moments[16] = -f[4]-f[5]+f[6]+f[7]+f[13]+f[14]-f[15]-f[16];
00095         moments[17] = f[4]-f[5]-f[8]-f[9]-f[13]+f[14]+f[17]+f[18];
00096         moments[18] = -f[6]+f[7]+f[8]-f[9]+f[15]-f[16]-f[17]+f[18];
00097         
00098     }
00099     
00101     static T mrtCollision( Array<T,Descriptor::q>& f,
00102                            const T &rhoBar, const Array<T,3> & j,
00103                            const T invM_S[19][19] )
00104     {
00105         
00106         Array<T,19> moments, momentsEq;
00107 
00108         computeMoments(moments,f);
00109         
00110         T jSqr = VectorTemplateImpl<T,3>::normSqr(j);
00111         
00112         computeEquilibrium(momentsEq,rhoBar,j,jSqr);
00113         
00114         T mom1 = moments[1] - momentsEq[1];
00115         T mom2 = moments[2] - momentsEq[2];
00116         T mom3 = moments[3] - momentsEq[3];
00117         T mom4 = moments[4] - momentsEq[4];
00118         T mom5 = moments[5] - momentsEq[5];
00119         T mom6 = moments[6] - momentsEq[6];
00120         T mom7 = moments[7] - momentsEq[7];
00121         T mom8 = moments[8] - momentsEq[8];
00122         T mom9 = moments[9] - momentsEq[9];
00123         T mom10 = moments[10] - momentsEq[10];
00124         T mom11 = moments[11] - momentsEq[11];
00125         T mom12 = moments[12] - momentsEq[12];
00126         T mom13 = moments[13] - momentsEq[13];
00127         T mom14 = moments[14] - momentsEq[14];
00128         T mom15 = moments[15] - momentsEq[15];
00129         T mom16 = moments[16];
00130         T mom17 = moments[17];
00131         T mom18 = moments[18];
00132 
00133         
00134         f[0] -= invM_S[0][1]*mom1
00135                 +invM_S[0][2]*mom2;
00136         
00137         f[1] -= invM_S[1][1]*mom1
00138                 +invM_S[1][2]*mom2
00139                 +invM_S[1][3]*mom3
00140                 +invM_S[1][4]*mom4
00141                 +invM_S[1][9]*mom9
00142                 +invM_S[1][10]*mom10;
00143         
00144         f[2] -= invM_S[2][1]*mom1
00145                 +invM_S[2][2]*mom2
00146                 +invM_S[2][5]*mom5
00147                 +invM_S[2][6]*mom6
00148                 +invM_S[2][9]*mom9
00149                 +invM_S[2][10]*mom10
00150                 +invM_S[2][11]*mom11
00151                 +invM_S[2][12]*mom12;
00152         
00153         f[3] -= invM_S[3][1]*mom1
00154                 +invM_S[3][2]*mom2
00155                 +invM_S[3][7]*mom7
00156                 +invM_S[3][8]*mom8
00157                 +invM_S[3][9]*mom9
00158                 +invM_S[3][10]*mom10
00159                 +invM_S[3][11]*mom11
00160                 +invM_S[3][12]*mom12;
00161         
00162         f[4] -= invM_S[4][1]*mom1
00163                 +invM_S[4][2]*mom2
00164                 +invM_S[4][3]*mom3
00165                 +invM_S[4][4]*mom4
00166                 +invM_S[4][5]*mom5
00167                 +invM_S[4][6]*mom6
00168                 +invM_S[4][9]*mom9
00169                 +invM_S[4][10]*mom10
00170                 +invM_S[4][11]*mom11
00171                 +invM_S[4][12]*mom12
00172                 +invM_S[4][13]*mom13
00173                 +invM_S[4][16]*mom16
00174                 +invM_S[4][17]*mom17;
00175         
00176         f[5] -= invM_S[5][1]*mom1
00177                 +invM_S[5][2]*mom2
00178                 +invM_S[5][3]*mom3
00179                 +invM_S[5][4]*mom4
00180                 +invM_S[5][5]*mom5
00181                 +invM_S[5][6]*mom6
00182                 +invM_S[5][9]*mom9
00183                 +invM_S[5][10]*mom10
00184                 +invM_S[5][11]*mom11
00185                 +invM_S[5][12]*mom12
00186                 +invM_S[5][13]*mom13
00187                 +invM_S[5][16]*mom16
00188                 +invM_S[5][17]*mom17;
00189         
00190         f[6] -= invM_S[6][1]*mom1
00191                 +invM_S[6][2]*mom2
00192                 +invM_S[6][3]*mom3
00193                 +invM_S[6][4]*mom4
00194                 +invM_S[6][7]*mom7
00195                 +invM_S[6][8]*mom8
00196                 +invM_S[6][9]*mom9
00197                 +invM_S[6][10]*mom10
00198                 +invM_S[6][11]*mom11
00199                 +invM_S[6][12]*mom12
00200                 +invM_S[6][15]*mom15
00201                 +invM_S[6][16]*mom16
00202                 +invM_S[6][18]*mom18;
00203         
00204         f[7] -= invM_S[7][1]*mom1
00205                 +invM_S[7][2]*mom2
00206                 +invM_S[7][3]*mom3
00207                 +invM_S[7][4]*mom4
00208                 +invM_S[7][7]*mom7
00209                 +invM_S[7][8]*mom8
00210                 +invM_S[7][9]*mom9
00211                 +invM_S[7][10]*mom10
00212                 +invM_S[7][11]*mom11
00213                 +invM_S[7][12]*mom12
00214                 +invM_S[7][15]*mom15
00215                 +invM_S[7][16]*mom16
00216                 +invM_S[7][18]*mom18;
00217         
00218         f[8] -= invM_S[8][1]*mom1
00219                 +invM_S[8][2]*mom2
00220                 +invM_S[8][5]*mom5
00221                 +invM_S[8][6]*mom6
00222                 +invM_S[8][7]*mom7
00223                 +invM_S[8][8]*mom8
00224                 +invM_S[8][9]*mom9
00225                 +invM_S[8][10]*mom10
00226                 +invM_S[8][14]*mom14
00227                 +invM_S[8][17]*mom17
00228                 +invM_S[8][18]*mom18;
00229         
00230         f[9] -= invM_S[9][1]*mom1
00231                 +invM_S[9][2]*mom2
00232                 +invM_S[9][5]*mom5
00233                 +invM_S[9][6]*mom6
00234                 +invM_S[9][7]*mom7
00235                 +invM_S[9][8]*mom8
00236                 +invM_S[9][9]*mom9
00237                 +invM_S[9][10]*mom10
00238                 +invM_S[9][14]*mom14
00239                 +invM_S[9][17]*mom17
00240                 +invM_S[9][18]*mom18;
00241         
00242         f[10] -= invM_S[10][1]*mom1
00243                 +invM_S[10][2]*mom2
00244                 +invM_S[10][3]*mom3
00245                 +invM_S[10][4]*mom4
00246                 +invM_S[10][9]*mom9
00247                 +invM_S[10][10]*mom10;
00248         
00249         f[11] -= invM_S[11][1]*mom1
00250                 +invM_S[11][2]*mom2
00251                 +invM_S[11][5]*mom5
00252                 +invM_S[11][6]*mom6
00253                 +invM_S[11][9]*mom9
00254                 +invM_S[11][10]*mom10
00255                 +invM_S[11][11]*mom11
00256                 +invM_S[11][12]*mom12;
00257         
00258         f[12] -= invM_S[12][1]*mom1
00259                 +invM_S[12][2]*mom2
00260                 +invM_S[12][7]*mom7
00261                 +invM_S[12][8]*mom8
00262                 +invM_S[12][9]*mom9
00263                 +invM_S[12][10]*mom10
00264                 +invM_S[12][11]*mom11
00265                 +invM_S[12][12]*mom12;
00266         
00267         f[13] -= invM_S[13][1]*mom1
00268                 +invM_S[13][2]*mom2
00269                 +invM_S[13][3]*mom3
00270                 +invM_S[13][4]*mom4
00271                 +invM_S[13][5]*mom5
00272                 +invM_S[13][6]*mom6
00273                 +invM_S[13][9]*mom9
00274                 +invM_S[13][10]*mom10
00275                 +invM_S[13][11]*mom11
00276                 +invM_S[13][12]*mom12
00277                 +invM_S[13][13]*mom13
00278                 +invM_S[13][16]*mom16
00279                 +invM_S[13][17]*mom17;
00280         
00281         f[14] -= invM_S[14][1]*mom1
00282                 +invM_S[14][2]*mom2
00283                 +invM_S[14][3]*mom3
00284                 +invM_S[14][4]*mom4
00285                 +invM_S[14][5]*mom5
00286                 +invM_S[14][6]*mom6
00287                 +invM_S[14][9]*mom9
00288                 +invM_S[14][10]*mom10
00289                 +invM_S[14][11]*mom11
00290                 +invM_S[14][12]*mom12
00291                 +invM_S[14][13]*mom13
00292                 +invM_S[14][16]*mom16
00293                 +invM_S[14][17]*mom17;
00294         
00295         f[15] -= invM_S[15][1]*mom1
00296                 +invM_S[15][2]*mom2
00297                 +invM_S[15][3]*mom3
00298                 +invM_S[15][4]*mom4
00299                 +invM_S[15][7]*mom7
00300                 +invM_S[15][8]*mom8
00301                 +invM_S[15][9]*mom9
00302                 +invM_S[15][10]*mom10
00303                 +invM_S[15][11]*mom11
00304                 +invM_S[15][12]*mom12
00305                 +invM_S[15][15]*mom15
00306                 +invM_S[15][16]*mom16
00307                 +invM_S[15][18]*mom18;
00308         
00309         f[16] -= invM_S[16][1]*mom1
00310                 +invM_S[16][2]*mom2
00311                 +invM_S[16][3]*mom3
00312                 +invM_S[16][4]*mom4
00313                 +invM_S[16][7]*mom7
00314                 +invM_S[16][8]*mom8
00315                 +invM_S[16][9]*mom9
00316                 +invM_S[16][10]*mom10
00317                 +invM_S[16][11]*mom11
00318                 +invM_S[16][12]*mom12
00319                 +invM_S[16][15]*mom15
00320                 +invM_S[16][16]*mom16
00321                 +invM_S[16][18]*mom18;
00322         
00323         f[17] -= invM_S[17][1]*mom1
00324                 +invM_S[17][2]*mom2
00325                 +invM_S[17][5]*mom5
00326                 +invM_S[17][6]*mom6
00327                 +invM_S[17][7]*mom7
00328                 +invM_S[17][8]*mom8
00329                 +invM_S[17][9]*mom9
00330                 +invM_S[17][10]*mom10
00331                 +invM_S[17][14]*mom14
00332                 +invM_S[17][17]*mom17
00333                 +invM_S[17][18]*mom18;
00334         
00335         f[18] -= invM_S[18][1]*mom1
00336                 +invM_S[18][2]*mom2
00337                 +invM_S[18][5]*mom5
00338                 +invM_S[18][6]*mom6
00339                 +invM_S[18][7]*mom7
00340                 +invM_S[18][8]*mom8
00341                 +invM_S[18][9]*mom9
00342                 +invM_S[18][10]*mom10
00343                 +invM_S[18][14]*mom14
00344                 +invM_S[18][17]*mom17
00345                 +invM_S[18][18]*mom18;
00346         
00347         return jSqr;
00348     }
00349     
00350     static void addGuoForce( Array<T,Descriptor::q>& f, const Array<T,Descriptor::d>& force,
00351                              Array<T,Descriptor::d> const& u,
00352                              T invM_S[Descriptor::q][Descriptor::q], T amplitude )
00353     {
00354         Array<T,Descriptor::q> forcing;
00355         T g_u   = force[0] * u[0] + force[1] * u[1] + force[2] * u[2];
00356         
00357         T mom1 = -(T)19*g_u;
00358         T mom2 = (T)5.5*g_u;
00359         T mom3 = -(T)0.5*force[0];
00360         T mom4 = force[0] / (T)3;
00361         T mom5 = -(T)0.5*force[1];
00362         T mom6 = force[1] / (T)3;
00363         T mom7 = -(T)0.5*force[2];
00364         T mom8 = force[2] / (T)3;
00365         T mom9 = -(2*force[0]*u[0]-force[1]*u[1]-force[2]*u[2]);
00366         T mom10 = (T)0.5*(2*force[0]*u[0]-force[1]*u[1]-force[2]*u[2]);
00367         T mom11 = -(force[1]*u[1]-force[2]*u[2]);
00368         T mom12 = (T)0.5*(force[1]*u[1]-force[2]*u[2]);
00369         T mom13 = -(T)0.5*(force[0]*u[1]+force[1]*u[0]);
00370         T mom14 = -(T)0.5*(force[2]*u[1]+force[1]*u[2]);
00371         T mom15 = -(T)0.5*(force[0]*u[2]+force[2]*u[0]);
00372         
00373         
00374         f[0] += invM_S[0][1]*mom1
00375                 +invM_S[0][2]*mom2;
00376         
00377         f[1] += invM_S[1][1]*mom1
00378                 +invM_S[1][2]*mom2
00379                 +invM_S[1][3]*mom3
00380                 +invM_S[1][4]*mom4
00381                 +invM_S[1][9]*mom9
00382                 +invM_S[1][10]*mom10;
00383         
00384         f[2] += invM_S[2][1]*mom1
00385                 +invM_S[2][2]*mom2
00386                 +invM_S[2][5]*mom5
00387                 +invM_S[2][6]*mom6
00388                 +invM_S[2][9]*mom9
00389                 +invM_S[2][10]*mom10
00390                 +invM_S[2][11]*mom11
00391                 +invM_S[2][12]*mom12;
00392         
00393         f[3] += invM_S[3][1]*mom1
00394                 +invM_S[3][2]*mom2
00395                 +invM_S[3][7]*mom7
00396                 +invM_S[3][8]*mom8
00397                 +invM_S[3][9]*mom9
00398                 +invM_S[3][10]*mom10
00399                 +invM_S[3][11]*mom11
00400                 +invM_S[3][12]*mom12;
00401         
00402         f[4] += invM_S[4][1]*mom1
00403                 +invM_S[4][2]*mom2
00404                 +invM_S[4][3]*mom3
00405                 +invM_S[4][4]*mom4
00406                 +invM_S[4][5]*mom5
00407                 +invM_S[4][6]*mom6
00408                 +invM_S[4][9]*mom9
00409                 +invM_S[4][10]*mom10
00410                 +invM_S[4][11]*mom11
00411                 +invM_S[4][12]*mom12
00412                 +invM_S[4][13]*mom13;
00413         
00414         f[5] += invM_S[5][1]*mom1
00415                 +invM_S[5][2]*mom2
00416                 +invM_S[5][3]*mom3
00417                 +invM_S[5][4]*mom4
00418                 +invM_S[5][5]*mom5
00419                 +invM_S[5][6]*mom6
00420                 +invM_S[5][9]*mom9
00421                 +invM_S[5][10]*mom10
00422                 +invM_S[5][11]*mom11
00423                 +invM_S[5][12]*mom12
00424                 +invM_S[5][13]*mom13;
00425         
00426         f[6] += invM_S[6][1]*mom1
00427                 +invM_S[6][2]*mom2
00428                 +invM_S[6][3]*mom3
00429                 +invM_S[6][4]*mom4
00430                 +invM_S[6][7]*mom7
00431                 +invM_S[6][8]*mom8
00432                 +invM_S[6][9]*mom9
00433                 +invM_S[6][10]*mom10
00434                 +invM_S[6][11]*mom11
00435                 +invM_S[6][12]*mom12
00436                 +invM_S[6][15]*mom15;
00437         
00438         f[7] += invM_S[7][1]*mom1
00439                 +invM_S[7][2]*mom2
00440                 +invM_S[7][3]*mom3
00441                 +invM_S[7][4]*mom4
00442                 +invM_S[7][7]*mom7
00443                 +invM_S[7][8]*mom8
00444                 +invM_S[7][9]*mom9
00445                 +invM_S[7][10]*mom10
00446                 +invM_S[7][11]*mom11
00447                 +invM_S[7][12]*mom12
00448                 +invM_S[7][15]*mom15;
00449         
00450         f[8] += invM_S[8][1]*mom1
00451                 +invM_S[8][2]*mom2
00452                 +invM_S[8][5]*mom5
00453                 +invM_S[8][6]*mom6
00454                 +invM_S[8][7]*mom7
00455                 +invM_S[8][8]*mom8
00456                 +invM_S[8][9]*mom9
00457                 +invM_S[8][10]*mom10
00458                 +invM_S[8][14]*mom14;
00459         
00460         f[9] += invM_S[9][1]*mom1
00461                 +invM_S[9][2]*mom2
00462                 +invM_S[9][5]*mom5
00463                 +invM_S[9][6]*mom6
00464                 +invM_S[9][7]*mom7
00465                 +invM_S[9][8]*mom8
00466                 +invM_S[9][9]*mom9
00467                 +invM_S[9][10]*mom10
00468                 +invM_S[9][14]*mom14;
00469         
00470         f[10] += invM_S[10][1]*mom1
00471                 +invM_S[10][2]*mom2
00472                 +invM_S[10][3]*mom3
00473                 +invM_S[10][4]*mom4
00474                 +invM_S[10][9]*mom9
00475                 +invM_S[10][10]*mom10;
00476         
00477         f[11] += invM_S[11][1]*mom1
00478                 +invM_S[11][2]*mom2
00479                 +invM_S[11][5]*mom5
00480                 +invM_S[11][6]*mom6
00481                 +invM_S[11][9]*mom9
00482                 +invM_S[11][10]*mom10
00483                 +invM_S[11][11]*mom11
00484                 +invM_S[11][12]*mom12;
00485         
00486         f[12] += invM_S[12][1]*mom1
00487                 +invM_S[12][2]*mom2
00488                 +invM_S[12][7]*mom7
00489                 +invM_S[12][8]*mom8
00490                 +invM_S[12][9]*mom9
00491                 +invM_S[12][10]*mom10
00492                 +invM_S[12][11]*mom11
00493                 +invM_S[12][12]*mom12;
00494         
00495         f[13] += invM_S[13][1]*mom1
00496                 +invM_S[13][2]*mom2
00497                 +invM_S[13][3]*mom3
00498                 +invM_S[13][4]*mom4
00499                 +invM_S[13][5]*mom5
00500                 +invM_S[13][6]*mom6
00501                 +invM_S[13][9]*mom9
00502                 +invM_S[13][10]*mom10
00503                 +invM_S[13][11]*mom11
00504                 +invM_S[13][12]*mom12
00505                 +invM_S[13][13]*mom13;
00506         
00507         f[14] += invM_S[14][1]*mom1
00508                 +invM_S[14][2]*mom2
00509                 +invM_S[14][3]*mom3
00510                 +invM_S[14][4]*mom4
00511                 +invM_S[14][5]*mom5
00512                 +invM_S[14][6]*mom6
00513                 +invM_S[14][9]*mom9
00514                 +invM_S[14][10]*mom10
00515                 +invM_S[14][11]*mom11
00516                 +invM_S[14][12]*mom12
00517                 +invM_S[14][13]*mom13;
00518         
00519         f[15] += invM_S[15][1]*mom1
00520                 +invM_S[15][2]*mom2
00521                 +invM_S[15][3]*mom3
00522                 +invM_S[15][4]*mom4
00523                 +invM_S[15][7]*mom7
00524                 +invM_S[15][8]*mom8
00525                 +invM_S[15][9]*mom9
00526                 +invM_S[15][10]*mom10
00527                 +invM_S[15][11]*mom11
00528                 +invM_S[15][12]*mom12
00529                 +invM_S[15][15]*mom15;
00530         
00531         f[16] += invM_S[16][1]*mom1
00532                 +invM_S[16][2]*mom2
00533                 +invM_S[16][3]*mom3
00534                 +invM_S[16][4]*mom4
00535                 +invM_S[16][7]*mom7
00536                 +invM_S[16][8]*mom8
00537                 +invM_S[16][9]*mom9
00538                 +invM_S[16][10]*mom10
00539                 +invM_S[16][11]*mom11
00540                 +invM_S[16][12]*mom12
00541                 +invM_S[16][15]*mom15;
00542         
00543         f[17] += invM_S[17][1]*mom1
00544                 +invM_S[17][2]*mom2
00545                 +invM_S[17][5]*mom5
00546                 +invM_S[17][6]*mom6
00547                 +invM_S[17][7]*mom7
00548                 +invM_S[17][8]*mom8
00549                 +invM_S[17][9]*mom9
00550                 +invM_S[17][10]*mom10
00551                 +invM_S[17][14]*mom14;
00552         
00553         f[18] += invM_S[18][1]*mom1
00554                 +invM_S[18][2]*mom2
00555                 +invM_S[18][5]*mom5
00556                 +invM_S[18][6]*mom6
00557                 +invM_S[18][7]*mom7
00558                 +invM_S[18][8]*mom8
00559                 +invM_S[18][9]*mom9
00560                 +invM_S[18][10]*mom10
00561                 +invM_S[18][14]*mom14;
00562                 
00563                 
00564         static const T oneOver6 = (T)1/(T)6;
00565         static const T oneOver12 = (T)1/(T)12;
00566         
00567         f[0] += -g_u;
00568         
00569         f[1] += oneOver6*(force[0]*(-(T)1+2*u[0])-force[1]*u[1]-force[2]*u[2]);
00570         
00571         f[2] += -oneOver6*(force[0]*u[0]+force[1]*((T)1-2*u[1])+force[2]*u[2]);
00572         
00573         f[3] += -oneOver6*(force[0]*u[0]
00574                             + force[1]*u[1]
00575                             + force[2]*((T)1-2*u[2]));
00576         
00577         f[4] += oneOver12*(force[0]*(-(T)1+2*u[0]+3*u[1])
00578                             + force[1]*(-(T)1+2*u[1]+3*u[0])
00579                             - force[2]*u[2]);
00580         
00581         f[5] += oneOver12*( force[0]*(-(T)1+2*u[0]-3*u[1])
00582                                 + force[1]*((T)1+2*u[1]-3*u[0])
00583                                 - force[2]*u[2]);
00584         
00585         f[6] += oneOver12*(force[0]*(-(T)1+2*u[0]+3*u[2])
00586                             - force[1]*u[1]
00587                             + force[2]*(-(T)1+2*u[2]+3*u[0]));
00588         
00589         f[7] += oneOver12*(force[0]*(-(T)1+2*u[0]-3*u[2])
00590                             - force[1]*u[1]
00591                             + force[2]*((T)1+2*u[2]-3*u[0]));
00592         
00593         f[8] += -oneOver12*(force[0]*u[0]
00594                                 + force[1]*((T)1-2*u[1]-3*u[2])
00595                                 + force[2]*((T)1-2*u[2]-3*u[1]));
00596         
00597         f[9] += -oneOver12*(force[0]*u[0]
00598                                 + force[1]*((T)1-2*u[1]+3*u[2])
00599                                 + force[2]*(-(T)1-2*u[2]+3*u[1]));
00600         
00601         f[10] += oneOver6*(force[0]*((T)1+2*u[0])
00602                                 -force[1]*u[1]
00603                                 -force[2]*u[2]);
00604         
00605         f[11] += -oneOver6*(force[0]*u[0]
00606                                 +force[1]*(-(T)1-2*u[1])
00607                                 +force[2]*u[2]);
00608         
00609         f[12] += -oneOver6*(force[0]*u[0]
00610                                 +force[1]*u[1]
00611                                 +force[2]*(-(T)1-2*u[2]));
00612         
00613         f[13] += oneOver12*(force[0]*((T)1+2*u[0]+3*u[1])
00614                                 +force[1]*((T)1+2*u[1]+3*u[0])
00615                                 -force[2]*u[2]);
00616         
00617         f[14] += oneOver12*(force[0]*((T)1+2*u[0]-3*u[1])
00618                                 +force[1]*(-(T)1+2*u[1]-3*u[0])
00619                                 -force[2]*u[2]);
00620         
00621         f[15] += oneOver12*(force[0]*((T)1+2*u[0]+3*u[2])
00622                                 -force[1]*u[1]
00623                                 +force[2]*((T)1+2*u[2]+3*u[0]));
00624         
00625         f[16] += oneOver12*(force[0]*((T)1+2*u[0]-3*u[2])
00626                                 -force[1]*u[1]
00627                                 +force[2]*(-(T)1+2*u[2]-3*u[0]));
00628         
00629         f[17] += -oneOver12*(force[0]*u[0]
00630                                 +force[1]*(-(T)1-2*u[1]-3*u[2])
00631                                 +force[2]*(-(T)1-2*u[2]-3*u[1]));
00632         
00633         f[18] += -oneOver12*(force[0]*u[0]
00634                                 +force[1]*(-(T)1-2*u[1]+3*u[2])
00635                                 +force[2]*((T)1-2*u[2]+3*u[1]));
00636     }
00637     
00639     static T mrtCollisionWithForce( Array<T,Descriptor::q>& f,
00640                                     const T &rhoBar, const Array<T,Descriptor::d> & u,
00641                                     T invM_S[Descriptor::q][Descriptor::q], 
00642                                     const Array<T,Descriptor::d> &force, T amplitude) 
00643     {
00644         Array<T,Descriptor::d> j = Descriptor::fullRho(rhoBar)*u;
00645         T jSqr = mrtCollision( f, rhoBar, j, invM_S );
00646         addGuoForce( f, force, u, invM_S, amplitude );
00647         
00648         return jSqr;
00649     }
00650     
00652     static void computeQuasiIncEquilibrium( Array<T,Descriptor::q>& momentsEq,
00653                                             T rhoBar, Array<T,3> const& j, T jSqr )
00654     {
00655         momentsEq[0] = rhoBar;
00656         momentsEq[1] = (T)19*jSqr-(T)11*rhoBar;
00657         momentsEq[2] = -(T)5.5*jSqr+(T)3*rhoBar;
00658         momentsEq[3] = j[0];
00659         momentsEq[4] = -((T)2/(T)3)*j[0];
00660         momentsEq[5] = j[1];
00661         momentsEq[6] = -((T)2/3)*j[1];
00662         momentsEq[7] = j[2];
00663         momentsEq[8] = -((T)2/(T)3)*j[2];
00664         momentsEq[9] = ((T)2*j[0]*j[0]-j[1]*j[1]-j[2]*j[2]);
00665         momentsEq[10] = (-j[0]*j[0]+(T)0.5*j[1]*j[1]+(T)0.5*j[2]*j[2]);
00666         momentsEq[11] = (j[1]*j[1]-j[2]*j[2]);
00667         momentsEq[12] = (-(T)0.5*j[1]*j[1]+(T)0.5*j[2]*j[2]);
00668         momentsEq[13] = j[1]*j[0];
00669         momentsEq[14] = j[2]*j[1];
00670         momentsEq[15] = j[2]*j[0];
00671         momentsEq[16] = T();
00672         momentsEq[17] = T();
00673         momentsEq[18] = T();
00674     }
00675     
00677     static T quasiIncMrtCollision( Array<T,Descriptor::q>& f,
00678                                    const T &rhoBar, const Array<T,3> & j,
00679                                    const T invM_S[19][19] )
00680     {
00681         
00682         Array<T,19> moments, momentsEq;
00683 
00684         computeMoments(moments,f);
00685         moments[0] = rhoBar;
00686         moments[3] = j[0];
00687         moments[5] = j[1];
00688         moments[7] = j[2];
00689         
00690         T jSqr = VectorTemplateImpl<T,3>::normSqr(j);
00691         
00692         computeQuasiIncEquilibrium(momentsEq,rhoBar,j,jSqr);
00693         
00694         T mom1 = moments[1] - momentsEq[1];
00695         T mom2 = moments[2] - momentsEq[2];
00696         T mom4 = moments[4] - momentsEq[4];
00697         T mom6 = moments[6] - momentsEq[6];
00698         T mom8 = moments[8] - momentsEq[8];
00699         T mom9 = moments[9] - momentsEq[9];
00700         T mom10 = moments[10] - momentsEq[10];
00701         T mom11 = moments[11] - momentsEq[11];
00702         T mom12 = moments[12] - momentsEq[12];
00703         T mom13 = moments[13] - momentsEq[13];
00704         T mom14 = moments[14] - momentsEq[14];
00705         T mom15 = moments[15] - momentsEq[15];
00706         T mom16 = moments[16];
00707         T mom17 = moments[17];
00708         T mom18 = moments[18];
00709 
00710         
00711         f[0] -= invM_S[0][1]*mom1
00712                 +invM_S[0][2]*mom2;
00713         
00714         f[1] -= invM_S[1][1]*mom1
00715                 +invM_S[1][2]*mom2
00716                 +invM_S[1][4]*mom4
00717                 +invM_S[1][9]*mom9
00718                 +invM_S[1][10]*mom10;
00719         
00720         f[2] -= invM_S[2][1]*mom1
00721                 +invM_S[2][2]*mom2
00722                 +invM_S[2][6]*mom6
00723                 +invM_S[2][9]*mom9
00724                 +invM_S[2][10]*mom10
00725                 +invM_S[2][11]*mom11
00726                 +invM_S[2][12]*mom12;
00727         
00728         f[3] -= invM_S[3][1]*mom1
00729                 +invM_S[3][2]*mom2
00730                 +invM_S[3][8]*mom8
00731                 +invM_S[3][9]*mom9
00732                 +invM_S[3][10]*mom10
00733                 +invM_S[3][11]*mom11
00734                 +invM_S[3][12]*mom12;
00735         
00736         f[4] -= invM_S[4][1]*mom1
00737                 +invM_S[4][2]*mom2
00738                 +invM_S[4][4]*mom4
00739                 +invM_S[4][6]*mom6
00740                 +invM_S[4][9]*mom9
00741                 +invM_S[4][10]*mom10
00742                 +invM_S[4][11]*mom11
00743                 +invM_S[4][12]*mom12
00744                 +invM_S[4][13]*mom13
00745                 +invM_S[4][16]*mom16
00746                 +invM_S[4][17]*mom17;
00747         
00748         f[5] -= invM_S[5][1]*mom1
00749                 +invM_S[5][2]*mom2
00750                 +invM_S[5][4]*mom4
00751                 +invM_S[5][6]*mom6
00752                 +invM_S[5][9]*mom9
00753                 +invM_S[5][10]*mom10
00754                 +invM_S[5][11]*mom11
00755                 +invM_S[5][12]*mom12
00756                 +invM_S[5][13]*mom13
00757                 +invM_S[5][16]*mom16
00758                 +invM_S[5][17]*mom17;
00759         
00760         f[6] -= invM_S[6][1]*mom1
00761                 +invM_S[6][2]*mom2
00762                 +invM_S[6][4]*mom4
00763                 +invM_S[6][8]*mom8
00764                 +invM_S[6][9]*mom9
00765                 +invM_S[6][10]*mom10
00766                 +invM_S[6][11]*mom11
00767                 +invM_S[6][12]*mom12
00768                 +invM_S[6][15]*mom15
00769                 +invM_S[6][16]*mom16
00770                 +invM_S[6][18]*mom18;
00771         
00772         f[7] -= invM_S[7][1]*mom1
00773                 +invM_S[7][2]*mom2
00774                 +invM_S[7][4]*mom4
00775                 +invM_S[7][8]*mom8
00776                 +invM_S[7][9]*mom9
00777                 +invM_S[7][10]*mom10
00778                 +invM_S[7][11]*mom11
00779                 +invM_S[7][12]*mom12
00780                 +invM_S[7][15]*mom15
00781                 +invM_S[7][16]*mom16
00782                 +invM_S[7][18]*mom18;
00783         
00784         f[8] -= invM_S[8][1]*mom1
00785                 +invM_S[8][2]*mom2
00786                 +invM_S[8][6]*mom6
00787                 +invM_S[8][8]*mom8
00788                 +invM_S[8][9]*mom9
00789                 +invM_S[8][10]*mom10
00790                 +invM_S[8][14]*mom14
00791                 +invM_S[8][17]*mom17
00792                 +invM_S[8][18]*mom18;
00793         
00794         f[9] -= invM_S[9][1]*mom1
00795                 +invM_S[9][2]*mom2
00796                 +invM_S[9][6]*mom6
00797                 +invM_S[9][8]*mom8
00798                 +invM_S[9][9]*mom9
00799                 +invM_S[9][10]*mom10
00800                 +invM_S[9][14]*mom14
00801                 +invM_S[9][17]*mom17
00802                 +invM_S[9][18]*mom18;
00803         
00804         f[10] -= invM_S[10][1]*mom1
00805                 +invM_S[10][2]*mom2
00806                 +invM_S[10][4]*mom4
00807                 +invM_S[10][9]*mom9
00808                 +invM_S[10][10]*mom10;
00809         
00810         f[11] -= invM_S[11][1]*mom1
00811                 +invM_S[11][2]*mom2
00812                 +invM_S[11][6]*mom6
00813                 +invM_S[11][9]*mom9
00814                 +invM_S[11][10]*mom10
00815                 +invM_S[11][11]*mom11
00816                 +invM_S[11][12]*mom12;
00817         
00818         f[12] -= invM_S[12][1]*mom1
00819                 +invM_S[12][2]*mom2
00820                 +invM_S[12][8]*mom8
00821                 +invM_S[12][9]*mom9
00822                 +invM_S[12][10]*mom10
00823                 +invM_S[12][11]*mom11
00824                 +invM_S[12][12]*mom12;
00825         
00826         f[13] -= invM_S[13][1]*mom1
00827                 +invM_S[13][2]*mom2
00828                 +invM_S[13][4]*mom4
00829                 +invM_S[13][6]*mom6
00830                 +invM_S[13][9]*mom9
00831                 +invM_S[13][10]*mom10
00832                 +invM_S[13][11]*mom11
00833                 +invM_S[13][12]*mom12
00834                 +invM_S[13][13]*mom13
00835                 +invM_S[13][16]*mom16
00836                 +invM_S[13][17]*mom17;
00837         
00838         f[14] -= invM_S[14][1]*mom1
00839                 +invM_S[14][2]*mom2
00840                 +invM_S[14][4]*mom4
00841                 +invM_S[14][6]*mom6
00842                 +invM_S[14][9]*mom9
00843                 +invM_S[14][10]*mom10
00844                 +invM_S[14][11]*mom11
00845                 +invM_S[14][12]*mom12
00846                 +invM_S[14][13]*mom13
00847                 +invM_S[14][16]*mom16
00848                 +invM_S[14][17]*mom17;
00849         
00850         f[15] -= invM_S[15][1]*mom1
00851                 +invM_S[15][2]*mom2
00852                 +invM_S[15][4]*mom4
00853                 +invM_S[15][8]*mom8
00854                 +invM_S[15][9]*mom9
00855                 +invM_S[15][10]*mom10
00856                 +invM_S[15][11]*mom11
00857                 +invM_S[15][12]*mom12
00858                 +invM_S[15][15]*mom15
00859                 +invM_S[15][16]*mom16
00860                 +invM_S[15][18]*mom18;
00861         
00862         f[16] -= invM_S[16][1]*mom1
00863                 +invM_S[16][2]*mom2
00864                 +invM_S[16][4]*mom4
00865                 +invM_S[16][8]*mom8
00866                 +invM_S[16][9]*mom9
00867                 +invM_S[16][10]*mom10
00868                 +invM_S[16][11]*mom11
00869                 +invM_S[16][12]*mom12
00870                 +invM_S[16][15]*mom15
00871                 +invM_S[16][16]*mom16
00872                 +invM_S[16][18]*mom18;
00873         
00874         f[17] -= invM_S[17][1]*mom1
00875                 +invM_S[17][2]*mom2
00876                 +invM_S[17][6]*mom6
00877                 +invM_S[17][8]*mom8
00878                 +invM_S[17][9]*mom9
00879                 +invM_S[17][10]*mom10
00880                 +invM_S[17][14]*mom14
00881                 +invM_S[17][17]*mom17
00882                 +invM_S[17][18]*mom18;
00883         
00884         f[18] -= invM_S[18][1]*mom1
00885                 +invM_S[18][2]*mom2
00886                 +invM_S[18][6]*mom6
00887                 +invM_S[18][8]*mom8
00888                 +invM_S[18][9]*mom9
00889                 +invM_S[18][10]*mom10
00890                 +invM_S[18][14]*mom14
00891                 +invM_S[18][17]*mom17
00892                 +invM_S[18][18]*mom18;
00893         
00894         return jSqr;
00895     }
00896 
00897 };
00898 
00899 
00900 }  // namespace plb
00901 
00902 #endif