hermitePolynomialsTemplates.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 
00030 #ifndef HERMITE_POLYNOMIALS_TEMPLATES_H
00031 #define HERMITE_POLYNOMIALS_TEMPLATES_H
00032 
00033 #include "core/globalDefs.h"
00034 #include "core/array.h"
00035 #include "geometricOperationTemplates.h"
00036 #include "momentTemplates.h"
00037 #include <algorithm>
00038 #include <vector>
00039 
00040 namespace plb {
00041 
00042 template <typename T, class Descriptor, int d> struct HermiteTemplateImpl;
00043 
00044 template <typename T, template<typename U> class Descriptor>
00045 struct HermiteTemplate {
00047     static const int d = Descriptor<T>::d;
00049     static T order0(plint iPop) {
00050         PLB_ASSERT(iPop < Descriptor<T>::q && iPop >= 0);
00051         return (T)1;
00052     }
00054     static Array<T,d> order1(plint iPop) {
00055         PLB_ASSERT(iPop < Descriptor<T>::q && iPop >= 0);
00056         return HermiteTemplateImpl<T,typename Descriptor<T>::BaseDescriptor,Descriptor<T>::d>::
00057             order1(iPop);
00058     }
00059     
00061     static Array<T,SymmetricTensorImpl<T,d>::n > order2(plint iPop) {
00062         PLB_ASSERT(iPop < Descriptor<T>::q && iPop >= 0);
00063         return HermiteTemplateImpl<T,typename Descriptor<T>::BaseDescriptor,Descriptor<T>::d>::
00064             order2(iPop);
00065     }
00066     
00069     static Array<T,SymmetricTensorImpl<T,d>::n > contractedOrder2(plint iPop) {
00070         PLB_ASSERT(iPop < Descriptor<T>::q && iPop >= 0);
00071         return HermiteTemplateImpl<T,typename Descriptor<T>::BaseDescriptor,Descriptor<T>::d>::
00072             contractedOrder2(iPop);
00073     }
00074     
00076     static Array<T,SymmetricRankThreeTensorImpl<T,d>::n > order3(plint iPop) {
00077         PLB_ASSERT(iPop < Descriptor<T>::q && iPop >= 0);
00078         return HermiteTemplateImpl<T,typename Descriptor<T>::BaseDescriptor,Descriptor<T>::d>::
00079             order3(iPop);
00080     }
00081     
00084     static Array<T,SymmetricRankThreeTensorImpl<T,d>::n > contractedOrder3(plint iPop) {
00085         PLB_ASSERT(iPop < Descriptor<T>::q && iPop >= 0);
00086         return HermiteTemplateImpl<T,typename Descriptor<T>::BaseDescriptor,Descriptor<T>::d>::
00087             contractedOrder3(iPop);
00088     }
00089     
00091     static Array<T,SymmetricRankFourTensorImpl<T,d>::n > order4(plint iPop) {
00092     PLB_ASSERT(iPop < Descriptor<T>::q && iPop >= 0);
00093         return HermiteTemplateImpl<T,typename Descriptor<T>::BaseDescriptor,Descriptor<T>::d>::
00094             order4(iPop);
00095     }
00096     
00099     static Array<T,SymmetricRankFourTensorImpl<T,d>::n > contractedOrder4(plint iPop) {
00100         PLB_ASSERT(iPop < Descriptor<T>::q && iPop >= 0);
00101         return HermiteTemplateImpl<T,typename Descriptor<T>::BaseDescriptor,Descriptor<T>::d>::
00102             contractedOrder4(iPop);
00103     }
00104     
00105     static plint numDecomposedVariables(plint order) {
00106         switch (order) {
00107             case 0 : {
00108                 return 1;
00109                 break;
00110             }
00111             case 1: {
00112                 return 1+Descriptor<T>::d;
00113                 break;
00114             }
00115             case 2: {
00116                 return 1+Descriptor<T>::d+SymmetricTensor<T,Descriptor>::n;
00117                 break;
00118             }
00119             case 3: {
00120                 return 1+Descriptor<T>::d+SymmetricTensor<T,Descriptor>::n+SymmetricRankThreeTensor<T,Descriptor>::n;
00121                 break;
00122             }
00123             case 4: {
00124                 return 1+Descriptor<T>::d+SymmetricTensor<T,Descriptor>::n
00125                                          +SymmetricRankThreeTensor<T,Descriptor>::n
00126                                          +SymmetricRankFourTensor<T,Descriptor>::n;
00127                 break;
00128             }
00129             default :
00130                 return 0;
00131         }
00132     }
00133     
00134     // specialization in 2D of the decomposition in hermite polynomials
00135     static std::vector<T> decompose(const Cell<T,Descriptor> &cell, plint order) {
00136         PLB_ASSERT(order >= 0 && order <= 4 );
00137         std::vector<T> coeffs;
00138         coeffs.resize(numDecomposedVariables(order));
00139         
00140         if (order >= 0) {
00141             coeffs[0] = momentTemplates<T,Descriptor>::get_rhoBar(cell);
00142             if (order >= 1) {
00143                 Array<T,d> a1;
00144                 momentTemplates<T,Descriptor>::get_j(cell,a1);
00145                 a1.to_cArray(&coeffs[1]);
00146             }
00147             if (order >= 2) {
00148                 Array<T,SymmetricTensorImpl<T,d>::n> a2;
00149                 a2.resetToZero();
00150                 for (plint iPop = 0; iPop < Descriptor<T>::q; ++iPop) {
00151                     Array<T,SymmetricTensorImpl<T,d>::n> H2 = order2(iPop);
00152                     for (plint iPi = 0; iPi < SymmetricTensorImpl<T,d>::n; ++iPi) {
00153                         a2[iPi] += H2[iPi]*cell[iPop];
00154                     }
00155                 }
00156                 a2.to_cArray(&coeffs[1+Descriptor<T>::d]);
00157                 if (order >= 3) {
00158                     Array<T,SymmetricRankThreeTensorImpl<T,d>::n> a3;
00159                     a3.resetToZero();
00160                     for (plint iPop = 0; iPop < Descriptor<T>::q; ++iPop) {
00161                         Array<T,SymmetricRankThreeTensorImpl<T,d>::n> H3 = order3(iPop);
00162                         for (plint iPi = 0; iPi < SymmetricRankThreeTensorImpl<T,d>::n; ++iPi) {
00163                             a3[iPi] += H3[iPi]*cell[iPop];
00164                         }
00165                     }
00166                     a3.to_cArray(&coeffs[1+Descriptor<T>::d]+SymmetricTensorImpl<T,d>::n);
00167                     if (order >= 4) {
00168                         Array<T,SymmetricRankFourTensorImpl<T,d>::n> a4;
00169                         a4.resetToZero();
00170                         for (plint iPop = 0; iPop < Descriptor<T>::q; ++iPop) {
00171                             Array<T,SymmetricRankFourTensorImpl<T,d>::n> H4 = order4(iPop);
00172                             for (plint iPi = 0; iPi < SymmetricRankFourTensorImpl<T,d>::n; ++iPi) {
00173                                 a4[iPi] += H4[iPi]*cell[iPop];
00174                             }
00175                         }
00176                         a4.to_cArray(&coeffs[1+Descriptor<T>::d]+SymmetricTensorImpl<T,d>::n+SymmetricRankThreeTensorImpl<T,d>::n);
00177                     }
00178                 }
00179             }
00180         }
00181         
00182         return coeffs;
00183     }
00184     
00185     // specialization in 2D of the decomposition in hermite polynomials
00186     static void recompose(const std::vector<T> &coeffs, plint order, Cell<T,Descriptor> &cell) {
00187         PLB_ASSERT(order >= 0 && order <= 4 );
00188         
00189         static T invCs4 = Descriptor<T>::invCs2* Descriptor<T>::invCs2 / (T)2;
00190         static T invCs6 = invCs4 * Descriptor<T>::invCs2 / (T)3;
00191         static T invCs8 = invCs6 * Descriptor<T>::invCs2 / (T)4;
00192         
00193         for (plint iPop = 0; iPop < Descriptor<T>::q; ++iPop) {
00194             if (order >= 0) {
00195                 cell[iPop] = coeffs[0];
00196                 if (order >= 1) {
00197                     Array<T,d> a1;
00198                     a1.from_cArray(&coeffs[1]);
00199                     T c_u = Descriptor<T>::c[iPop][0] * a1[0];
00200                     for (plint iD = 1; iD < d; ++iD) c_u += Descriptor<T>::c[iPop][iD] * a1[iD];
00201                     cell[iPop] += Descriptor<T>::invCs2 * c_u;
00202                 }
00203                 if (order >= 2) {
00204                     Array<T,SymmetricTensor<T,Descriptor>::n> a2;
00205                     a2.from_cArray(&coeffs[1+d]);
00206                     Array<T,SymmetricTensor<T,Descriptor>::n> H2 = order2(iPop);
00207                     cell[iPop] += invCs4 * SymmetricTensor<T,Descriptor>::contractIndexes(H2,a2);
00208                     if (order >= 3) {
00209                         Array<T,SymmetricRankThreeTensor<T,Descriptor>::n> a3;
00210                         Array<T,SymmetricRankThreeTensor<T,Descriptor>::n> H3 = order3(iPop);
00211                         a3.from_cArray(&coeffs[1+d+SymmetricTensor<T,Descriptor>::n]);
00212                         cell[iPop] += invCs6 * SymmetricRankThreeTensor<T,Descriptor>::contractIndexes(H3,a3);
00213                         if (order >= 4) {
00214                             Array<T,SymmetricRankFourTensor<T,Descriptor>::n> a4;
00215                             Array<T,SymmetricRankFourTensor<T,Descriptor>::n> H4 = order4(iPop);
00216                             a4.from_cArray(&coeffs[1+d+SymmetricTensor<T,Descriptor>::n
00217                                                       +SymmetricRankThreeTensor<T,Descriptor>::n]);
00218                             cell[iPop] += invCs8 * SymmetricRankFourTensor<T,Descriptor>::contractIndexes(H4,a4);
00219                         }
00220                     }
00221                 }
00222             }
00223             cell[iPop] *= Descriptor<T>::t[iPop];
00224         }
00225     }
00226 };
00227 
00228 template <typename T, class Descriptor>
00229 struct HermiteTemplateImpl<T,Descriptor,2> {
00230     
00231     typedef Descriptor D;
00232     
00233     static T order0(plint iPop) {
00234         return (T)1;
00235     }
00236     
00237     static Array<T,2> order1(plint iPop) {
00238         Array<T,2> H1(D::c[iPop][0],D::c[iPop][1]);
00239         return H1;
00240     }
00241     
00243     static Array<T,SymmetricTensorImpl<T,2>::n > order2(plint iPop) {
00244         Array<T,SymmetricTensorImpl<T,2>::n> H2;
00245         H2[SymmetricTensorImpl<T,2>::xx] = D::c[iPop][0]*D::c[iPop][0]-D::cs2;
00246         H2[SymmetricTensorImpl<T,2>::xy] = D::c[iPop][0]*D::c[iPop][1];
00247         H2[SymmetricTensorImpl<T,2>::yy] = D::c[iPop][1]*D::c[iPop][1]-D::cs2;
00248         
00249         return H2;
00250     }
00251     
00254     static Array<T,SymmetricTensorImpl<T,2>::n > contractedOrder2(plint iPop) {
00255         Array<T,SymmetricTensorImpl<T,2>::n> H2;
00256         H2[SymmetricTensorImpl<T,2>::xx] = D::c[iPop][0]*D::c[iPop][0]-D::cs2;
00257         H2[SymmetricTensorImpl<T,2>::xy] = (T)2*D::c[iPop][0]*D::c[iPop][1];
00258         H2[SymmetricTensorImpl<T,2>::yy] = D::c[iPop][1]*D::c[iPop][1]-D::cs2;
00259         
00260         return H2;
00261     }
00262     
00264     static Array<T,SymmetricRankThreeTensorImpl<T,2>::n > order3(plint iPop) {
00265         plint cx2 = D::c[iPop][0]*D::c[iPop][0];
00266         plint cy2 = D::c[iPop][1]*D::c[iPop][1];
00267         
00268         Array<T,SymmetricRankThreeTensorImpl<T,2>::n> H3;
00269         H3[SymmetricRankThreeTensorImpl<T,2>::xxx] = D::c[iPop][0]*(cx2-(T)3*D::cs2);
00270         H3[SymmetricRankThreeTensorImpl<T,2>::xxy] = D::c[iPop][1]*(cx2-D::cs2);
00271         H3[SymmetricRankThreeTensorImpl<T,2>::xyy] = D::c[iPop][0]*(cy2-D::cs2);
00272         H3[SymmetricRankThreeTensorImpl<T,2>::yyy] = D::c[iPop][1]*(cy2-(T)3*D::cs2);
00273         
00274         return H3;
00275     }
00276     
00279     static Array<T,SymmetricRankThreeTensorImpl<T,2>::n > contractedOrder3(plint iPop) {
00280         plint cx2 = D::c[iPop][0]*D::c[iPop][0];
00281         plint cy2 = D::c[iPop][1]*D::c[iPop][1];
00282         
00283         Array<T,SymmetricRankThreeTensorImpl<T,2>::n> H3;
00284         H3[SymmetricRankThreeTensorImpl<T,2>::xxx] = D::c[iPop][0]*(cx2-(T)3*D::cs2);
00285         H3[SymmetricRankThreeTensorImpl<T,2>::xxy] = (T)3*D::c[iPop][1]*(cx2-D::cs2);
00286         H3[SymmetricRankThreeTensorImpl<T,2>::xyy] = (T)3*D::c[iPop][0]*(cy2-D::cs2);
00287         H3[SymmetricRankThreeTensorImpl<T,2>::yyy] = D::c[iPop][1]*(cy2-(T)3*D::cs2);
00288         
00289         return H3;
00290     }
00291     
00293     static Array<T,SymmetricRankFourTensorImpl<T,2>::n > order4(plint iPop) {
00294         plint cx2 = D::c[iPop][0]*D::c[iPop][0];
00295         plint cy2 = D::c[iPop][1]*D::c[iPop][1];
00296         
00297         Array<T,SymmetricRankFourTensorImpl<T,2>::n> H4;
00298         H4[SymmetricRankFourTensorImpl<T,2>::xxxx] = 
00299             cx2 * (cx2 - (T)6*D::cs2) + (T)3*D::cs2*D::cs2;
00300         H4[SymmetricRankFourTensorImpl<T,2>::xxxy] = 
00301             D::c[iPop][0] * D::c[iPop][1] * (cx2 - (T)3*D::cs2);
00302         H4[SymmetricRankFourTensorImpl<T,2>::xxyy] = 
00303             cx2*cy2 - D::cs2*cx2 - D::cs2*cy2 + D::cs2*D::cs2;;
00304         H4[SymmetricRankFourTensorImpl<T,2>::xyyy] = 
00305             D::c[iPop][0] * D::c[iPop][1] * (cy2 - (T)3*D::cs2);
00306         H4[SymmetricRankFourTensorImpl<T,2>::yyyy] = 
00307             cy2 * (cy2 - (T)6*D::cs2) + (T)3*D::cs2*D::cs2;
00308         
00309         return H4;
00310     }
00311     
00314     static Array<T,SymmetricRankFourTensorImpl<T,2>::n > contractedOrder4(plint iPop) {
00315         plint cx2 = D::c[iPop][0]*D::c[iPop][0];
00316         plint cy2 = D::c[iPop][1]*D::c[iPop][1];
00317         
00318         Array<T,SymmetricRankFourTensorImpl<T,2>::n> H4;
00319         H4[SymmetricRankFourTensorImpl<T,2>::xxxx] = 
00320             (T)cx2 * (cx2 - (T)6*D::cs2) + (T)3*D::cs2*D::cs2;
00321         H4[SymmetricRankFourTensorImpl<T,2>::xxxy] = 
00322             (T)4*D::c[iPop][0] * D::c[iPop][1] * (cx2 - (T)3*D::cs2);
00323         H4[SymmetricRankFourTensorImpl<T,2>::xxyy] = 
00324             (T)6*(cx2*cy2 - D::cs2*cx2 - D::cs2*cy2 + D::cs2*D::cs2);
00325         H4[SymmetricRankFourTensorImpl<T,2>::xyyy] = 
00326             (T)4*D::c[iPop][0] * D::c[iPop][1] * (cy2 - (T)3*D::cs2);
00327         H4[SymmetricRankFourTensorImpl<T,2>::yyyy] = 
00328             (T)cy2 * (cy2 - (T)6*D::cs2) + (T)3*D::cs2*D::cs2;
00329         
00330         return H4;
00331     }
00332 };
00333 
00334 
00335 template <typename T, class Descriptor>
00336 struct HermiteTemplateImpl<T,Descriptor,3> {
00337     
00338     typedef Descriptor D;
00339     
00340     static T order0(plint iPop) {
00341         return (T)1;
00342     }
00343     
00344     static Array<T,3> order1(plint iPop) {
00345         Array<T,3> H1(D::c[iPop][0],D::c[iPop][1],D::c[iPop][2]);
00346         return H1;
00347     }
00348     
00350     static Array<T,SymmetricTensorImpl<T,3>::n > order2(plint iPop) {
00351         Array<T,SymmetricTensorImpl<T,3>::n> H2;
00352         H2[SymmetricTensorImpl<T,3>::xx] = D::c[iPop][0]*D::c[iPop][0]-D::cs2;
00353         H2[SymmetricTensorImpl<T,3>::xy] = D::c[iPop][0]*D::c[iPop][1];
00354         H2[SymmetricTensorImpl<T,3>::xz] = D::c[iPop][0]*D::c[iPop][2];
00355         H2[SymmetricTensorImpl<T,3>::yy] = D::c[iPop][1]*D::c[iPop][1]-D::cs2;
00356         H2[SymmetricTensorImpl<T,3>::yz] = D::c[iPop][1]*D::c[iPop][2];
00357         H2[SymmetricTensorImpl<T,3>::zz] = D::c[iPop][2]*D::c[iPop][2]-D::cs2;
00358         
00359         return H2;
00360     }
00361     
00364     static Array<T,SymmetricTensorImpl<T,3>::n > contractedOrder2(plint iPop) {
00365         Array<T,SymmetricTensorImpl<T,3>::n> H2;
00366         H2[SymmetricTensorImpl<T,3>::xx] =      D::c[iPop][0]*D::c[iPop][0]-D::cs2;
00367         H2[SymmetricTensorImpl<T,3>::xy] = (T)2*D::c[iPop][0]*D::c[iPop][1];
00368         H2[SymmetricTensorImpl<T,3>::xz] = (T)2*D::c[iPop][0]*D::c[iPop][2];
00369         H2[SymmetricTensorImpl<T,3>::yy] =      D::c[iPop][1]*D::c[iPop][1]-D::cs2;
00370         H2[SymmetricTensorImpl<T,3>::yz] = (T)2*D::c[iPop][1]*D::c[iPop][2];
00371         H2[SymmetricTensorImpl<T,3>::zz] =      D::c[iPop][2]*D::c[iPop][2]-D::cs2;
00372         
00373         return H2;
00374     }
00375     
00377     static Array<T,SymmetricRankThreeTensorImpl<T,3>::n > order3(plint iPop) {
00378         plint cx2 = D::c[iPop][0]*D::c[iPop][0];
00379         plint cy2 = D::c[iPop][1]*D::c[iPop][1];
00380         plint cz2 = D::c[iPop][2]*D::c[iPop][2];
00381         
00382         Array<T,SymmetricRankThreeTensorImpl<T,3>::n> H3;
00383         H3[SymmetricRankThreeTensorImpl<T,3>::xxx] = D::c[iPop][0]*(cx2-(T)3*D::cs2);
00384         H3[SymmetricRankThreeTensorImpl<T,3>::xxy] = D::c[iPop][1]*(cx2-D::cs2);
00385         H3[SymmetricRankThreeTensorImpl<T,3>::xxz] = D::c[iPop][2]*(cx2-D::cs2);
00386         
00387         H3[SymmetricRankThreeTensorImpl<T,3>::xyy] = D::c[iPop][0]*(cy2-D::cs2);
00388         H3[SymmetricRankThreeTensorImpl<T,3>::xzz] = D::c[iPop][0]*(cz2-D::cs2);
00389         
00390         H3[SymmetricRankThreeTensorImpl<T,3>::xyz] = D::c[iPop][0]*D::c[iPop][1]*D::c[iPop][2];
00391         
00392         H3[SymmetricRankThreeTensorImpl<T,3>::yyz] = D::c[iPop][2]*(cy2-D::cs2);
00393         H3[SymmetricRankThreeTensorImpl<T,3>::yzz] = D::c[iPop][1]*(cz2-D::cs2);
00394 
00395         H3[SymmetricRankThreeTensorImpl<T,3>::yyy] = D::c[iPop][1]*(cy2-(T)3*D::cs2);
00396         
00397         H3[SymmetricRankThreeTensorImpl<T,3>::zzz] = D::c[iPop][2]*(cz2-(T)3*D::cs2);
00398         
00399         return H3;
00400     }
00401     
00404     static Array<T,SymmetricRankThreeTensorImpl<T,3>::n > contractedOrder3(plint iPop) {
00405         plint cx2 = D::c[iPop][0]*D::c[iPop][0];
00406         plint cy2 = D::c[iPop][1]*D::c[iPop][1];
00407         plint cz2 = D::c[iPop][2]*D::c[iPop][2];
00408         
00409         Array<T,SymmetricRankThreeTensorImpl<T,3>::n> H3;
00410         H3[SymmetricRankThreeTensorImpl<T,3>::xxx] =      D::c[iPop][0]*(cx2-(T)3*D::cs2);
00411         H3[SymmetricRankThreeTensorImpl<T,3>::xxy] = (T)3*D::c[iPop][1]*(cx2-D::cs2);
00412         H3[SymmetricRankThreeTensorImpl<T,3>::xxz] = (T)3*D::c[iPop][2]*(cx2-D::cs2);
00413         
00414         H3[SymmetricRankThreeTensorImpl<T,3>::xyy] = (T)3*D::c[iPop][0]*(cy2-D::cs2);
00415         H3[SymmetricRankThreeTensorImpl<T,3>::xzz] = (T)3*D::c[iPop][0]*(cz2-D::cs2);
00416         
00417         H3[SymmetricRankThreeTensorImpl<T,3>::xyz] = (T)6*D::c[iPop][0]*D::c[iPop][1]*D::c[iPop][2];
00418         
00419         H3[SymmetricRankThreeTensorImpl<T,3>::yyz] = (T)3*D::c[iPop][2]*(cy2-D::cs2);
00420         H3[SymmetricRankThreeTensorImpl<T,3>::yzz] = (T)3*D::c[iPop][1]*(cz2-D::cs2);
00421         
00422         H3[SymmetricRankThreeTensorImpl<T,3>::yyy] =      D::c[iPop][1]*(cy2-(T)3*D::cs2);
00423         
00424         H3[SymmetricRankThreeTensorImpl<T,3>::zzz] =      D::c[iPop][2]*(cz2-(T)3*D::cs2);
00425         
00426         return H3;
00427     }
00428     
00430     static Array<T,SymmetricRankFourTensorImpl<T,3>::n > order4(plint iPop) {
00431         plint cx2 = D::c[iPop][0]*D::c[iPop][0];
00432         plint cy2 = D::c[iPop][1]*D::c[iPop][1];
00433         plint cz2 = D::c[iPop][2]*D::c[iPop][2];
00434         
00435         plint cxy = D::c[iPop][0]*D::c[iPop][1];
00436         plint cxz = D::c[iPop][0]*D::c[iPop][2];
00437         plint cyz = D::c[iPop][1]*D::c[iPop][2];
00438         
00439         Array<T,SymmetricRankFourTensorImpl<T,3>::n> H4;
00440         H4[SymmetricRankFourTensorImpl<T,3>::xxxx] = cx2 * (cx2 - (T)6*D::cs2) + (T)3*D::cs2*D::cs2;
00441         H4[SymmetricRankFourTensorImpl<T,3>::yyyy] = cy2 * (cy2 - (T)6*D::cs2) + (T)3*D::cs2*D::cs2;
00442         H4[SymmetricRankFourTensorImpl<T,3>::zzzz] = cz2 * (cz2 - (T)6*D::cs2) + (T)3*D::cs2*D::cs2;
00443         
00444         H4[SymmetricRankFourTensorImpl<T,3>::xxxy] = D::c[iPop][0] * D::c[iPop][1] * (cx2 - (T)3*D::cs2);
00445         H4[SymmetricRankFourTensorImpl<T,3>::xxxz] = D::c[iPop][0] * D::c[iPop][2] * (cx2 - (T)3*D::cs2);
00446         H4[SymmetricRankFourTensorImpl<T,3>::xyyy] = D::c[iPop][0] * D::c[iPop][1] * (cy2 - (T)3*D::cs2);
00447         H4[SymmetricRankFourTensorImpl<T,3>::xzzz] = D::c[iPop][0] * D::c[iPop][2] * (cz2 - (T)3*D::cs2);
00448         H4[SymmetricRankFourTensorImpl<T,3>::yyyz] = D::c[iPop][2] * D::c[iPop][1] * (cy2 - (T)3*D::cs2);
00449         H4[SymmetricRankFourTensorImpl<T,3>::yzzz] = D::c[iPop][2] * D::c[iPop][1] * (cz2 - (T)3*D::cs2);
00450         
00451         H4[SymmetricRankFourTensorImpl<T,3>::xxyy] = cx2*cy2 - D::cs2*(cx2 - cy2) + D::cs2*D::cs2;
00452         H4[SymmetricRankFourTensorImpl<T,3>::xxzz] = cx2*cz2 - D::cs2*(cx2 - cz2) + D::cs2*D::cs2;
00453         H4[SymmetricRankFourTensorImpl<T,3>::yyzz] = cy2*cz2 - D::cs2*(cy2 - cz2) + D::cs2*D::cs2;
00454         
00455         H4[SymmetricRankFourTensorImpl<T,3>::xxyz] = D::c[iPop][1] * D::c[iPop][2] * (cx2 - D::cs2);
00456         H4[SymmetricRankFourTensorImpl<T,3>::xyyz] = D::c[iPop][0] * D::c[iPop][2] * (cy2 - D::cs2);
00457         H4[SymmetricRankFourTensorImpl<T,3>::xyzz] = D::c[iPop][0] * D::c[iPop][1] * (cz2 - D::cs2);
00458         
00459         return H4;
00460     }
00461     
00464     static Array<T,SymmetricRankFourTensorImpl<T,3>::n > contractedOrder4(plint iPop) {
00465         plint cx2 = D::c[iPop][0]*D::c[iPop][0];
00466         plint cy2 = D::c[iPop][1]*D::c[iPop][1];
00467         plint cz2 = D::c[iPop][2]*D::c[iPop][2];
00468         
00469         plint cxy = D::c[iPop][0]*D::c[iPop][1];
00470         plint cxz = D::c[iPop][0]*D::c[iPop][2];
00471         plint cyz = D::c[iPop][1]*D::c[iPop][2];
00472         
00473         Array<T,SymmetricRankFourTensorImpl<T,3>::n> H4;
00474         H4[SymmetricRankFourTensorImpl<T,3>::xxxx] = cx2 * (cx2 - (T)6*D::cs2) + (T)3*D::cs2*D::cs2;
00475         H4[SymmetricRankFourTensorImpl<T,3>::yyyy] = cy2 * (cy2 - (T)6*D::cs2) + (T)3*D::cs2*D::cs2;
00476         H4[SymmetricRankFourTensorImpl<T,3>::zzzz] = cz2 * (cz2 - (T)6*D::cs2) + (T)3*D::cs2*D::cs2;
00477         
00478         H4[SymmetricRankFourTensorImpl<T,3>::xxxy] = (T)4*D::c[iPop][0] * D::c[iPop][1] * (cx2 - (T)3*D::cs2);
00479         H4[SymmetricRankFourTensorImpl<T,3>::xxxz] = (T)4*D::c[iPop][0] * D::c[iPop][2] * (cx2 - (T)3*D::cs2);
00480         H4[SymmetricRankFourTensorImpl<T,3>::xyyy] = (T)4*D::c[iPop][0] * D::c[iPop][1] * (cy2 - (T)3*D::cs2);
00481         H4[SymmetricRankFourTensorImpl<T,3>::xzzz] = (T)4*D::c[iPop][0] * D::c[iPop][2] * (cz2 - (T)3*D::cs2);
00482         H4[SymmetricRankFourTensorImpl<T,3>::yyyz] = (T)4*D::c[iPop][2] * D::c[iPop][1] * (cy2 - (T)3*D::cs2);
00483         H4[SymmetricRankFourTensorImpl<T,3>::yzzz] = (T)4*D::c[iPop][2] * D::c[iPop][1] * (cz2 - (T)3*D::cs2);
00484         
00485         H4[SymmetricRankFourTensorImpl<T,3>::xxyy] = (T)6*cx2*cy2 - D::cs2*(cx2 - cy2) + D::cs2*D::cs2;
00486         H4[SymmetricRankFourTensorImpl<T,3>::xxzz] = (T)6*cx2*cz2 - D::cs2*(cx2 - cz2) + D::cs2*D::cs2;
00487         H4[SymmetricRankFourTensorImpl<T,3>::yyzz] = (T)6*cy2*cz2 - D::cs2*(cy2 - cz2) + D::cs2*D::cs2;
00488         
00489         H4[SymmetricRankFourTensorImpl<T,3>::xxyz] = (T)12*D::c[iPop][1] * D::c[iPop][2] * (cx2 - D::cs2);
00490         H4[SymmetricRankFourTensorImpl<T,3>::xyyz] = (T)12*D::c[iPop][0] * D::c[iPop][2] * (cy2 - D::cs2);
00491         H4[SymmetricRankFourTensorImpl<T,3>::xyzz] = (T)12*D::c[iPop][0] * D::c[iPop][1] * (cz2 - D::cs2);
00492         
00493         return H4;
00494     }
00495 };
00496 
00497 }  // namespace plb
00498 
00499 #endif