generalizedBoundaryDynamics.hh

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 GENERALIZED_BOUNDARY_DYNAMICS_HH
00030 #define GENERALIZED_BOUNDARY_DYNAMICS_HH
00031 
00032 #include "generalizedBoundaryDynamics.h"
00033 #include "generalizedBoundaryDynamicsSolvers.h"
00034 #include "generalizedBoundaryDynamicsSolvers.hh"
00035 #include "generalizedIncompressibleBoundaryTemplates.h"
00036 #include "generalizedCompressibleBoundaryTemplates.h"
00037 #include "core/cell.h"
00038 #include "core/dynamicsIdentifiers.h"
00039 #include "latticeBoltzmann/indexTemplates.h"
00040 #include <Eigen3/Core>
00041 #include <Eigen3/LU>
00042 #include <Eigen3/Cholesky>
00043 
00044 namespace plb {
00045 
00046 /* *************** Class GeneralizedVelocityBoundaryDynamics ************* */
00047 
00048 template<typename T, template<typename U> class Descriptor>
00049 int GeneralizedVelocityBoundaryDynamics<T,Descriptor>::id =
00050     meta::registerGeneralDynamics<T,Descriptor,GeneralizedVelocityBoundaryDynamics<T,Descriptor> >
00051                                      ("Boundary_GeneralizedVelocity");
00052 
00053 template<typename T, template<typename U> class Descriptor>
00054 GeneralizedVelocityBoundaryDynamics<T,Descriptor>::
00055     GeneralizedVelocityBoundaryDynamics(Dynamics<T,Descriptor>* baseDynamics_,
00056                                         std::vector<plint> missingIndices_,
00057                                         bool automaticPrepareCollision)
00058         : StoreVelocityDynamics<T,Descriptor>(baseDynamics_,automaticPrepareCollision),
00059           missingIndices(missingIndices_)
00060 { 
00061     knownIndices = indexTemplates::remainingIndexes<Descriptor<T> >(missingIndices);
00062 }
00063 
00064 template<typename T, template<typename U> class Descriptor>
00065 GeneralizedVelocityBoundaryDynamics<T,Descriptor>::
00066     GeneralizedVelocityBoundaryDynamics(HierarchicUnserializer& unserializer)
00067         : StoreVelocityDynamics<T,Descriptor>(0, false)
00068 {
00069     unserialize(unserializer);
00070 }
00071 
00072 template<typename T, template<typename U> class Descriptor>
00073 GeneralizedVelocityBoundaryDynamics<T,Descriptor>*
00074     GeneralizedVelocityBoundaryDynamics<T,Descriptor>::clone() const
00075 {
00076     return new GeneralizedVelocityBoundaryDynamics<T,Descriptor>(*this);
00077 }
00078  
00079 template<typename T, template<typename U> class Descriptor>
00080 int GeneralizedVelocityBoundaryDynamics<T,Descriptor>::getId() const {
00081     return id;
00082 }
00083 
00084 template<typename T, template<typename U> class Descriptor>
00085 void GeneralizedVelocityBoundaryDynamics<T,Descriptor>::serialize(HierarchicSerializer& serializer) const
00086 {
00087     serializer.addValue((int)(missingIndices.size()));
00088     serializer.addValues(missingIndices);
00089     serializer.addValue((int)(knownIndices.size()));
00090     serializer.addValues(knownIndices);
00091     StoreVelocityDynamics<T,Descriptor>::serialize(serializer);
00092 }
00093 
00094 template<typename T, template<typename U> class Descriptor>
00095 void GeneralizedVelocityBoundaryDynamics<T,Descriptor>::unserialize(HierarchicUnserializer& unserializer)
00096 {
00097     PLB_PRECONDITION( unserializer.getId() == this->getId() );
00098     missingIndices.resize(unserializer.readValue<int>());
00099     unserializer.readValues(missingIndices);
00100     knownIndices.resize(unserializer.readValue<int>());
00101     unserializer.readValues(knownIndices);
00102     StoreVelocityDynamics<T,Descriptor>::unserialize(unserializer);
00103 }
00104 
00105 template<typename T, template<typename U> class Descriptor>
00106 void GeneralizedVelocityBoundaryDynamics<T,Descriptor>::
00107     computeUlb(Cell<T,Descriptor> const& cell, Array<T,Descriptor<T>::d> &uLb) const {
00108         
00109     this -> computeVelocity(cell, uLb);
00110     for (plint iDim = 0; iDim < Descriptor<T>::d; ++iDim) {
00111         uLb[iDim] -= 0.5*getExternalForceComponent(cell,iDim);
00112     }
00113 }
00114 
00115 template<typename T, template<typename U> class Descriptor>
00116 void GeneralizedVelocityBoundaryDynamics<T,Descriptor>::
00117     completePopulations(Cell<T,Descriptor>& cell) const
00118 {
00119     Array<T,Descriptor<T>::d> uLb;
00120     computeUlb(cell, uLb);
00121     
00122     DirichletVelocityBoundarySolver<T,Descriptor> bc(missingIndices, knownIndices, uLb);
00123     bc.apply(cell,this->getBaseDynamics());
00124 }
00125 
00126 /* *************** Class GeneralizedMassConservingVelocityBoundaryDynamics ************* */
00127 
00128 template<typename T, template<typename U> class Descriptor>
00129 int GeneralizedMassConservingVelocityBoundaryDynamics<T,Descriptor>::id =
00130     meta::registerGeneralDynamics<T,Descriptor,GeneralizedMassConservingVelocityBoundaryDynamics<T,Descriptor> >
00131                                      ("Boundary_MassConservingGeneralizedVelocity");
00132 
00133 template<typename T, template<typename U> class Descriptor>
00134 GeneralizedMassConservingVelocityBoundaryDynamics<T,Descriptor>::
00135     GeneralizedMassConservingVelocityBoundaryDynamics(Dynamics<T,Descriptor>* baseDynamics_,
00136                                         std::vector<plint> missingIndices_,
00137                                         std::vector<plint> knownIndices_,
00138                                         std::vector<plint> inGoingIndices_,
00139                                         bool automaticPrepareCollision)
00140         : StoreVelocityDynamics<T,Descriptor>(baseDynamics_,automaticPrepareCollision),
00141             missingIndices(missingIndices_), knownIndices(knownIndices_), inGoingIndices(inGoingIndices_)
00142 {  }
00143 
00144 template<typename T, template<typename U> class Descriptor>
00145 GeneralizedMassConservingVelocityBoundaryDynamics<T,Descriptor>::
00146     GeneralizedMassConservingVelocityBoundaryDynamics(HierarchicUnserializer& unserializer)
00147         : StoreVelocityDynamics<T,Descriptor>(0, false)
00148 {
00149     unserialize(unserializer);
00150 }
00151 
00152 template<typename T, template<typename U> class Descriptor>
00153 GeneralizedMassConservingVelocityBoundaryDynamics<T,Descriptor>*
00154     GeneralizedMassConservingVelocityBoundaryDynamics<T,Descriptor>::clone() const
00155 {
00156     return new GeneralizedMassConservingVelocityBoundaryDynamics<T,Descriptor>(*this);
00157 }
00158  
00159 template<typename T, template<typename U> class Descriptor>
00160 int GeneralizedMassConservingVelocityBoundaryDynamics<T,Descriptor>::getId() const {
00161     return id;
00162 }
00163 
00164 template<typename T, template<typename U> class Descriptor>
00165 void GeneralizedMassConservingVelocityBoundaryDynamics<T,Descriptor>::serialize(HierarchicSerializer& serializer) const
00166 {
00167     serializer.addValue((int)(missingIndices.size()));
00168     serializer.addValues(missingIndices);
00169     serializer.addValue((int)(knownIndices.size()));
00170     serializer.addValues(knownIndices);
00171     serializer.addValue((int)(inGoingIndices.size()));
00172     serializer.addValues(inGoingIndices);
00173     StoreVelocityDynamics<T,Descriptor>::serialize(serializer);
00174 }
00175 
00176 template<typename T, template<typename U> class Descriptor>
00177 void GeneralizedMassConservingVelocityBoundaryDynamics<T,Descriptor>::unserialize(HierarchicUnserializer& unserializer)
00178 {
00179     PLB_PRECONDITION( unserializer.getId() == this->getId() );
00180     missingIndices.resize(unserializer.readValue<int>());
00181     unserializer.readValues(missingIndices);
00182     knownIndices.resize(unserializer.readValue<int>());
00183     unserializer.readValues(knownIndices);
00184     inGoingIndices.resize(unserializer.readValue<int>());
00185     unserializer.readValues(inGoingIndices);
00186     StoreVelocityDynamics<T,Descriptor>::unserialize(unserializer);
00187 }
00188 
00189 template<typename T, template<typename U> class Descriptor>
00190 void GeneralizedMassConservingVelocityBoundaryDynamics<T,Descriptor>::
00191     computeUlb(Cell<T,Descriptor> const& cell, Array<T,Descriptor<T>::d> &uLb) const {
00192         
00193     this -> computeVelocity(cell, uLb);
00194     for (plint iDim = 0; iDim < Descriptor<T>::d; ++iDim) {
00195         uLb[iDim] -= 0.5*getExternalForceComponent(cell,iDim);
00196     }
00197 }
00198 
00199 template<typename T, template<typename U> class Descriptor>
00200 void GeneralizedMassConservingVelocityBoundaryDynamics<T,Descriptor>::
00201     completePopulations(Cell<T,Descriptor>& cell) const
00202 {
00203     Array<T,Descriptor<T>::d> uLb;
00204     computeUlb(cell, uLb);
00205         
00206     DirichletMassConservingVelocityBoundarySolver<T,Descriptor> bc(missingIndices, knownIndices, 
00207                                                                    inGoingIndices, uLb);
00208     bc.apply(cell,this->getBaseDynamics());
00209 }
00210 
00211 /* *************** Class GeneralizedDensityBoundaryDynamics ************* */
00212 
00213 template<typename T, template<typename U> class Descriptor,
00214          int direction, int orientation>
00215 int GeneralizedDensityBoundaryDynamics<T,Descriptor,direction,orientation>::id =
00216     meta::registerGeneralDynamics<T,Descriptor, GeneralizedDensityBoundaryDynamics<T,Descriptor,direction,orientation> >
00217             ( std::string("Boundary_GeneralizedDensity_")+util::val2str(direction) +
00218               std::string("_")+util::val2str(orientation) );
00219 
00220 template<typename T, template<typename U> class Descriptor,
00221          int direction, int orientation>
00222 GeneralizedDensityBoundaryDynamics<T,Descriptor,direction,orientation>::
00223     GeneralizedDensityBoundaryDynamics(Dynamics<T,Descriptor>* baseDynamics_, 
00224                                        std::vector<plint> missingIndices_, 
00225                                        bool automaticPrepareCollision_)
00226         : DensityDirichletBoundaryDynamics<T,Descriptor,direction,orientation>(baseDynamics_, automaticPrepareCollision_),
00227           missingIndices(missingIndices_)
00228 { 
00229     knownIndices = indexTemplates::remainingIndexes<Descriptor<T> >(missingIndices);
00230 }
00231 
00232 template<typename T, template<typename U> class Descriptor,
00233             int direction, int orientation>
00234 GeneralizedDensityBoundaryDynamics<T,Descriptor,direction,orientation>::
00235     GeneralizedDensityBoundaryDynamics(HierarchicUnserializer& unserializer)
00236         : DensityDirichletBoundaryDynamics<T,Descriptor,direction,orientation>(0, false)
00237 {
00238     unserialize(unserializer);
00239 }
00240 
00241 template<typename T, template<typename U> class Descriptor,
00242          int direction, int orientation>
00243 GeneralizedDensityBoundaryDynamics<T,Descriptor,direction,orientation>*
00244     GeneralizedDensityBoundaryDynamics<T,Descriptor,direction,orientation>::clone() const
00245 {
00246     return new GeneralizedDensityBoundaryDynamics<T,Descriptor,direction,orientation>(*this);
00247 }
00248  
00249 template<typename T, template<typename U> class Descriptor,
00250          int direction, int orientation>
00251 int GeneralizedDensityBoundaryDynamics<T,Descriptor,direction,orientation>::getId() const {
00252     return id;
00253 }
00254  
00255 template<typename T, template<typename U> class Descriptor,
00256          int direction, int orientation>
00257 void GeneralizedDensityBoundaryDynamics<T,Descriptor,direction,orientation>::serialize(HierarchicSerializer& serializer) const
00258 {
00259     serializer.addValue((int)(missingIndices.size()));
00260     serializer.addValues(missingIndices);
00261     serializer.addValue((int)(knownIndices.size()));
00262     serializer.addValues(knownIndices);
00263 //     serializer.addValues(u);
00264 //     serializer.addValues(PiNeq);
00265     DensityDirichletBoundaryDynamics<T,Descriptor,direction,orientation>::serialize(serializer);
00266 }
00267 
00268 template<typename T, template<typename U> class Descriptor,
00269          int direction, int orientation>
00270 void GeneralizedDensityBoundaryDynamics<T,Descriptor,direction,orientation>::unserialize(HierarchicUnserializer& unserializer)
00271 {
00272     PLB_PRECONDITION( unserializer.getId() == this->getId() );
00273     missingIndices.resize(unserializer.readValue<int>());
00274     unserializer.readValues(missingIndices);
00275     knownIndices.resize(unserializer.readValue<int>());
00276     unserializer.readValues(knownIndices);
00277 //     unserializer.readValues(u);
00278 //     unserializer.readValues(PiNeq);
00279     DensityDirichletBoundaryDynamics<T,Descriptor,direction,orientation>::unserialize(unserializer);
00280 }
00281 
00282 template<typename T, template<typename U> class Descriptor,
00283          int direction, int orientation>
00284 void GeneralizedDensityBoundaryDynamics<T,Descriptor,direction,orientation>::
00285     completePopulations(Cell<T,Descriptor>& cell) const
00286 {
00287     T rhoBar;
00288     Array<T,Descriptor<T>::d> u,j;
00289     this -> computeRhoBarJ(cell, rhoBar, j);
00290     T rho = Descriptor<T>::fullRho(rhoBar);
00291     T invRho = Descriptor<T>::invRho(rhoBar);
00292     u = invRho * j;
00293     
00294     Array<T,SymmetricTensor<T,Descriptor>::n> PiNeq;
00295     PiNeq.resetToZero();
00296     
00297     DirichletDensityBoundarySolver<T,Descriptor,direction> bc(missingIndices, knownIndices, rho, u, PiNeq, 1.0e-7);
00298     bc.apply(cell,this->getBaseDynamics());
00299 
00300 }
00301 
00302 /* ********** Class GeneralizedVelocityTemperatureBoundaryDynamics ********* */
00303 
00304 template<typename T, template<typename U> class Descriptor>
00305 int GeneralizedVelocityTemperatureBoundaryDynamics<T,Descriptor>::id =
00306     meta::registerGeneralDynamics<T,Descriptor,GeneralizedVelocityTemperatureBoundaryDynamics<T,Descriptor> >
00307                                      ("Boundary_GeneralizedVelocityTemperature");
00308 
00309 template<typename T, template<typename U> class Descriptor>
00310 GeneralizedVelocityTemperatureBoundaryDynamics<T,Descriptor>::
00311     GeneralizedVelocityTemperatureBoundaryDynamics(Dynamics<T,Descriptor>* baseDynamics_,
00312                                         std::vector<plint> missingIndices_, bool massConserving_,
00313                                         bool automaticPrepareCollision)
00314         : StoreTemperatureAndVelocityDynamics<T,Descriptor>(baseDynamics_,automaticPrepareCollision),
00315           missingIndices(missingIndices_), massConserving(massConserving_)
00316 { 
00317     knownIndices = indexTemplates::remainingIndexes<Descriptor<T> >(missingIndices);
00318 }
00319 
00320 template<typename T, template<typename U> class Descriptor>
00321 GeneralizedVelocityTemperatureBoundaryDynamics<T,Descriptor>::
00322     GeneralizedVelocityTemperatureBoundaryDynamics(Dynamics<T,Descriptor>* baseDynamics_,
00323                                         std::vector<plint> missingIndices_, std::vector<plint> knownIndices_, 
00324                                         bool massConserving_,
00325                                         bool automaticPrepareCollision)
00326         : StoreTemperatureAndVelocityDynamics<T,Descriptor>(baseDynamics_,automaticPrepareCollision),
00327             missingIndices(missingIndices_), knownIndices(knownIndices_), massConserving(massConserving_)
00328 {  }
00329 
00330 template<typename T, template<typename U> class Descriptor>
00331 GeneralizedVelocityTemperatureBoundaryDynamics<T,Descriptor>::
00332     GeneralizedVelocityTemperatureBoundaryDynamics(HierarchicUnserializer& unserializer)
00333         : StoreTemperatureAndVelocityDynamics<T,Descriptor>(0, false)
00334 {
00335     this->unserialize(unserializer);
00336 }
00337 
00338 template<typename T, template<typename U> class Descriptor>
00339 GeneralizedVelocityTemperatureBoundaryDynamics<T,Descriptor>*
00340     GeneralizedVelocityTemperatureBoundaryDynamics<T,Descriptor>::clone() const
00341 {
00342     return new GeneralizedVelocityTemperatureBoundaryDynamics<T,Descriptor>(*this);
00343 }
00344  
00345 template<typename T, template<typename U> class Descriptor>
00346 int GeneralizedVelocityTemperatureBoundaryDynamics<T,Descriptor>::getId() const {
00347     return id;
00348 }
00349 
00350 template<typename T, template<typename U> class Descriptor>
00351 void GeneralizedVelocityTemperatureBoundaryDynamics<T,Descriptor>::serialize(HierarchicSerializer& serializer) const
00352 {
00353     serializer.addValue((int)(missingIndices.size()));
00354     serializer.addValues(missingIndices);
00355     serializer.addValue((int)(knownIndices.size()));
00356     serializer.addValues(knownIndices);
00357     serializer.addValue(massConserving);
00358     StoreTemperatureAndVelocityDynamics<T,Descriptor>::serialize(serializer);
00359 }
00360 
00361 template<typename T, template<typename U> class Descriptor>
00362 void GeneralizedVelocityTemperatureBoundaryDynamics<T,Descriptor>::unserialize(HierarchicUnserializer& unserializer)
00363 {
00364     PLB_PRECONDITION( unserializer.getId() == this->getId() );
00365     missingIndices.resize(unserializer.readValue<int>());
00366     unserializer.readValues(missingIndices);
00367     knownIndices.resize(unserializer.readValue<int>());
00368     unserializer.readValues(knownIndices);
00369     massConserving = unserializer.readValue<bool>();
00370     StoreTemperatureAndVelocityDynamics<T,Descriptor>::unserialize(unserializer);
00371 }
00372 
00373 template<typename T, template<typename U> class Descriptor>
00374 void GeneralizedVelocityTemperatureBoundaryDynamics<T,Descriptor>::
00375     computeUlb(Cell<T,Descriptor> const& cell, Array<T,Descriptor<T>::d> &uLb) const {
00376         
00377     this -> computeVelocity(cell, uLb);
00378     for (plint iDim = 0; iDim < Descriptor<T>::d; ++iDim) {
00379         uLb[iDim] -= 0.5*getExternalForceComponent(cell,iDim);
00380     }
00381 }
00382 
00383 template<typename T, template<typename U> class Descriptor>
00384 void GeneralizedVelocityTemperatureBoundaryDynamics<T,Descriptor>::
00385     completePopulations(Cell<T,Descriptor>& cell) const
00386 {
00387     Array<T,Descriptor<T>::d> uLb;
00388     computeUlb(cell, uLb);
00389     T thetaBar = this->computeTemperature(cell)*Descriptor<T>::invCs2-(T)1;
00390     
00391     T rho;
00392     Array<T,SymmetricTensor<T,Descriptor>::n> PiNeq;
00393     Array<T,SymmetricRankThreeTensor<T,Descriptor>::n> qNeq;
00394 //     generalizedComprTempBoundaryTemplates<T,Descriptor>::
00395 //         solveLinearSystemTrLessPiNeq(cell, uLb, thetaBar,missingIndices, knownIndices, rho, PiNeq, qNeq, massConserving);
00396         
00397     generalizedComprTempBoundaryTemplates<T,Descriptor>::
00398         solveLinearSystem(cell, uLb, thetaBar,missingIndices, knownIndices, rho, PiNeq, qNeq, massConserving);
00399         
00400     T rhoBar = Descriptor<T>::rhoBar(rho);
00401     
00402     Array<T,Descriptor<T>::d> j;
00403     for (plint iDim = 0; iDim < Descriptor<T>::d; ++iDim) {
00404         j[iDim] = rho * uLb[iDim];
00405     }
00406     T jSqr = VectorTemplate<T,Descriptor>::normSqr(j);
00407     
00408     // TODO not elegant. Regularize functions should be made more generic
00409     typedef Descriptor<T> L;
00410     Array<T,Descriptor<T>::q> f, feq;
00411     for (plint iPop=0; iPop < L::q; ++iPop) {
00412         cell[iPop] = this->computeEquilibrium(iPop, rhoBar, j, jSqr, thetaBar);
00413         cell[iPop] += offEquilibriumTemplates<T,Descriptor>::fromPiAndQtoFneq(iPop, PiNeq, qNeq);
00414     }
00415 }
00416 
00417 /* *************** Class GeneralizedNextToWallBoundaryDynamics ************* */
00418 
00419 template<typename T, template<typename U> class Descriptor>
00420 int GeneralizedNextToBoundaryDynamics<T,Descriptor>::id =
00421     meta::registerGeneralDynamics<T,Descriptor,GeneralizedNextToBoundaryDynamics<T,Descriptor> >
00422         ("Boundary_GeneralizedNextToBoundaryDynamics");
00423 
00424 template<typename T, template<typename U> class Descriptor>
00425 GeneralizedNextToBoundaryDynamics<T,Descriptor>::
00426     GeneralizedNextToBoundaryDynamics(Dynamics<T,Descriptor>* baseDynamics_,
00427                                         std::vector<plint> missingIndices_,
00428                                         bool automaticPrepareCollision)
00429         : BoundaryCompositeDynamics<T,Descriptor>(baseDynamics_,automaticPrepareCollision),
00430             missingIndices(missingIndices_)
00431 { 
00432     knownIndices = indexTemplates::remainingIndexes<Descriptor<T> >(missingIndices);
00433 
00434     rho = (T)1;
00435     thetaBar = T();
00436     u.resetToZero();
00437     PiNeq.resetToZero();
00438     Qneq.resetToZero();
00439 }
00440 
00441 template<typename T, template<typename U> class Descriptor>
00442 GeneralizedNextToBoundaryDynamics<T,Descriptor>::
00443     GeneralizedNextToBoundaryDynamics(Dynamics<T,Descriptor>* baseDynamics_,
00444                                       std::vector<plint> missingIndices_,
00445                                       std::vector<plint> knownIndices_,
00446                                       bool automaticPrepareCollision)
00447         : BoundaryCompositeDynamics<T,Descriptor>(baseDynamics_,automaticPrepareCollision),
00448             missingIndices(missingIndices_), knownIndices(knownIndices_)
00449 {  }
00450 
00451 template<typename T, template<typename U> class Descriptor>
00452 GeneralizedNextToBoundaryDynamics<T,Descriptor>::
00453     GeneralizedNextToBoundaryDynamics(HierarchicUnserializer& unserializer)
00454         : BoundaryCompositeDynamics<T,Descriptor>(0, false)
00455 {
00456     this->unserialize(unserializer);
00457 }
00458 
00459 template<typename T, template<typename U> class Descriptor>
00460 GeneralizedNextToBoundaryDynamics<T,Descriptor>*
00461     GeneralizedNextToBoundaryDynamics<T,Descriptor>::clone() const
00462 {
00463     return new GeneralizedNextToBoundaryDynamics<T,Descriptor>(*this);
00464 }
00465  
00466 template<typename T, template<typename U> class Descriptor>
00467 int GeneralizedNextToBoundaryDynamics<T,Descriptor>::getId() const {
00468     return id;
00469 }
00470 
00471 template<typename T, template<typename U> class Descriptor>
00472 void GeneralizedNextToBoundaryDynamics<T,Descriptor>::serialize(HierarchicSerializer& serializer) const
00473 {
00474     serializer.addValue((int)(missingIndices.size()));
00475     serializer.addValues(missingIndices);
00476     serializer.addValue((int)(knownIndices.size()));
00477     serializer.addValues(knownIndices);
00478     BoundaryCompositeDynamics<T,Descriptor>::serialize(serializer);
00479 }
00480 
00481 template<typename T, template<typename U> class Descriptor>
00482 void GeneralizedNextToBoundaryDynamics<T,Descriptor>::unserialize(HierarchicUnserializer& unserializer)
00483 {
00484     PLB_PRECONDITION( unserializer.getId() == this->getId() );
00485     missingIndices.resize(unserializer.readValue<int>());
00486     unserializer.readValues(missingIndices);
00487     knownIndices.resize(unserializer.readValue<int>());
00488     unserializer.readValues(knownIndices);
00489     BoundaryCompositeDynamics<T,Descriptor>::unserialize(unserializer);
00490 }
00491 
00492 template<typename T, template<typename U> class Descriptor>
00493 void GeneralizedNextToBoundaryDynamics<T,Descriptor>::
00494     completePopulations(Cell<T,Descriptor>& cell) 
00495 {
00496     
00497     T epsilon = 1.0e-10;
00498     T resSum;
00499     bool converged = generalizedComprTempBoundaryTemplates<T,Descriptor>::
00500         iterativelySolveSystem(cell, rho, u, thetaBar, PiNeq, Qneq, knownIndices, missingIndices, epsilon, resSum);
00501 //     bool converged = generalizedComprTempBoundaryTemplates<T,Descriptor>::
00502 //         iterativelySolveSystemTrLessPiNeq(cell, rho, u, thetaBar, PiNeq, Qneq, knownIndices, missingIndices, epsilon, resSum);
00503     if (!converged) {
00504         pcout << "Never converged... Exiting." << std::endl;
00505         exit(-1);
00506     }
00507     Array<T,Descriptor<T>::d> j;
00508     j = rho*u;
00509     T rhoBar = Descriptor<T>::rhoBar(rho);
00510     T jSqr = VectorTemplate<T,Descriptor>::normSqr(j);
00511     for (plint iPop=0; iPop < Descriptor<T>::q; ++iPop) {
00512         cell[iPop] = this->computeEquilibrium(iPop, rhoBar, j, jSqr, thetaBar);
00513         cell[iPop] += offEquilibriumTemplates<T,Descriptor>::fromPiAndQtoFneq(iPop, PiNeq, Qneq);
00514     }
00515 }
00516 
00517 }  // namespace plb
00518 
00519 #endif  // GENERALIZED_BOUNDARY_DYNAMICS_HH