cell.h

Go to the documentation of this file.

/* This file is part of the Palabos library.
 *
 * Copyright (C) 2011 FlowKit Sarl
 * Avenue de Chailly 23
 * 1012 Lausanne, Switzerland
 * E-mail contact: contact@flowkit.com
 *
 * The most recent release of Palabos can be downloaded at 
 * <http://www.palabos.org/>
 *
 * The library Palabos is free software: you can redistribute it and/or
 * modify it under the terms of the GNU Affero General Public License as
 * published by the Free Software Foundation, either version 3 of the
 * License, or (at your option) any later version.
 *
 * The library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU Affero General Public License for more details.
 *
 * You should have received a copy of the GNU Affero General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
*/

#ifndef CELL_H
#define CELL_H

#include "core/globalDefs.h"
#include "core/plbDebug.h"
#include "latticeBoltzmann/externalFields.h"
#include "core/dynamics.h"
#include "core/array.h"
#include "core/hierarchicSerializer.h"

namespace plb {

template<typename T, typename ExternalField>
class ExternalFieldArray {
public:
    T* get(plint index) {
        PLB_PRECONDITION( index < ExternalField::numScalars );
        return data+index;
    }
    T const* get(plint index) const {
        PLB_PRECONDITION( index < ExternalField::numScalars );
        return data+index;
    }
private:
    T data[ExternalField::numScalars];
};

template<typename T>
class ExternalFieldArray<T,descriptors::NoExternalField> {
public:
    T* get(pluint index) {
        PLB_PRECONDITION( false );
        static T data = T();
        return &data;
    }
    T const* get(pluint index) const {
        PLB_PRECONDITION( false );
        static T data = T();
        return &data;
    }
};

template<typename T, template<typename U> class Descriptor>
struct CellInfo {
    static const int n = Descriptor<T>::numPop + Descriptor<T>::ExternalField::numScalars;
};


template<typename T, template<typename U> class Descriptor>
class Cell {
public:
    typedef ExternalFieldArray<T, typename Descriptor<T>::ExternalField> External;
public:
    Cell();
    Cell(Dynamics<T,Descriptor>* dynamics_);
public:

    T& operator[](plint iPop) {
        PLB_PRECONDITION( iPop < Descriptor<T>::numPop );
        return f[iPop];
    }

    T const& operator[](plint iPop) const {
        PLB_PRECONDITION( iPop < Descriptor<T>::numPop );
        return f[iPop];
    }
    Array<T,Descriptor<T>::numPop>& getRawPopulations() {
        return f;
    }
    Array<T,Descriptor<T>::numPop> const& getRawPopulations() const {
        return f;
    }

    Cell<T,Descriptor>& attributeF(Cell<T,Descriptor> const& rhs) {
        f = rhs.getRawPopulations();
        return *this;
    };

    Cell<T,Descriptor>& attributeValues(Cell<T,Descriptor> const& rhs) {
        attributeF(rhs);
        for (plint iExt=0; iExt < Descriptor<T>::ExternalField::numScalars; ++iExt) {
            *external.get(iExt) = *rhs.external.get(iExt);
        }
        return *this;
    };
    T* getExternal(plint offset) {
        PLB_PRECONDITION( offset < Descriptor<T>::ExternalField::numScalars );
        return external.get(offset);
    }
    T const* getExternal(plint offset) const {
        PLB_PRECONDITION( offset < Descriptor<T>::ExternalField::numScalars );
        return external.get(offset);
    }
    Dynamics<T,Descriptor> const& getDynamics() const;
    Dynamics<T,Descriptor>& getDynamics();
    bool takesStatistics() const {
        return takesStat;
    }
    void specifyStatisticsStatus(bool status) {
        takesStat = status;
    }
private:

    void attributeDynamics(Dynamics<T,Descriptor>* dynamics_);
    // Declare the BlockLatticeXD as a friend, to enable access to attributeDynamics.
    template<typename T_, template<typename U_> class Descriptor_> friend class BlockLattice2D;
    template<typename T_, template<typename U_> class Descriptor_> friend class BlockLattice3D;
#ifdef PLB_MPI_PARALLEL
    template<typename T_, template<typename U_> class Descriptor_> friend class ParallelCellAccess2D;
    template<typename T_, template<typename U_> class Descriptor_> friend class ParallelCellAccess3D;
#endif

// The following helper functions forward the function call
// to the Dynamics object
public:
    void collide(BlockStatistics& statistics) {
        PLB_PRECONDITION( dynamics );
        dynamics->collide(*this, statistics);
    }

    T computeEquilibrium(plint iPop, T rhoBar, Array<T,Descriptor<T>::d> const& j,
                         T jSqr, T thetaBar=T()) const
    {
        PLB_PRECONDITION( dynamics );
        return dynamics->computeEquilibrium(iPop, rhoBar, j, jSqr, thetaBar);
    }

    void regularize(T rhoBar, Array<T,Descriptor<T>::d> const& j, T jSqr,
                    Array<T,SymmetricTensor<T,Descriptor>::n> const& PiNeq, T thetaBar=T() )
    {
        PLB_PRECONDITION( dynamics );
        dynamics->regularize(*this, rhoBar, j, jSqr, PiNeq, thetaBar);
    }


    T computeDensity() const {
        PLB_PRECONDITION( dynamics );
        return dynamics->computeDensity(*this);
    }


    T computePressure() const {
        PLB_PRECONDITION( dynamics );
        return dynamics->computePressure(*this);
    }


    void computeVelocity(Array<T,Descriptor<T>::d>& u) const {
        PLB_PRECONDITION( dynamics );
        dynamics->computeVelocity(*this, u);
    }


    T computeTemperature() const {
        PLB_PRECONDITION( dynamics );
        return dynamics->computeTemperature(*this);
    }


    void computeDeviatoricStress (
            Array<T,SymmetricTensor<T,Descriptor>::n>& PiNeq ) const
    {
        PLB_PRECONDITION( dynamics );
        dynamics->computeDeviatoricStress(*this, PiNeq);
    }


    void computeHeatFlux(Array<T,Descriptor<T>::d>& q) const {
        PLB_PRECONDITION( dynamics );
        dynamics->computeHeatFlux(*this, q);
    }


    void computeMoment(plint momentId, T* moment) {
        PLB_PRECONDITION( dynamics );
        dynamics->computeMoment(*this, momentId, moment);
    }


    void getPopulations(Array<T,Descriptor<T>::numPop>& f) const {
        PLB_PRECONDITION( dynamics );
        dynamics->getPopulations(*this, f);
    }


    void getExternalField(plint pos, plint size, T* ext) const {
        PLB_PRECONDITION( dynamics );
        dynamics->getExternalField(*this, pos, size, ext);
    }   


    void defineDensity(T rho) {
        PLB_PRECONDITION( dynamics );
        dynamics->defineDensity(*this, rho);
    }


    void defineVelocity(Array<T,Descriptor<T>::d> const& u) {
        PLB_PRECONDITION( dynamics );
        dynamics->defineVelocity(*this, u);
    }


    void defineTemperature(T temperature) {
        PLB_PRECONDITION( dynamics );
        dynamics->defineTemperature(*this, temperature);
    }


    void defineHeatFlux(Array<T,Descriptor<T>::d> const& q) {
        PLB_PRECONDITION( dynamics );
        dynamics->defineHeatFlux(*this, q);
    }


    void defineDeviatoricStress (
            Array<T,SymmetricTensor<T,Descriptor>::n> const& PiNeq )
    {
        PLB_PRECONDITION( dynamics );
        dynamics->defineDeviatoricStress(*this, PiNeq);
    }


    void defineMoment(plint momentId, T const* value)
    {
        PLB_PRECONDITION( dynamics );
        dynamics->defineMoment(*this, momentId, value);
    }


    void setPopulations(Array<T,Descriptor<T>::numPop> const& f) {
        PLB_PRECONDITION( dynamics );
        dynamics->setPopulations(*this, f);
    }

    void setExternalField(plint pos, plint size, const T* ext) {
        PLB_PRECONDITION( dynamics );
        dynamics->setExternalField(*this, pos, size, ext);
    }
    void revert();
    void serialize(char* data) const;
    void unSerialize(char const* data);
private:
    void iniPop();
    void iniExternal();
private:
    Array<T,Descriptor<T>::numPop> f;         
    External                       external;  
    bool                           takesStat; 
    Dynamics<T,Descriptor>*        dynamics;  
};

template<typename T, template<typename U> class Descriptor>
void iniCellAtEquilibrium(Cell<T,Descriptor>& cell, T density, Array<T,Descriptor<T>::d> const& velocity);

template<typename T, template<typename U> class Descriptor>
void iniCellAtEquilibrium(Cell<T,Descriptor>& cell, T density, Array<T,Descriptor<T>::d> const& velocity, T temperature);

}  // namespace plb

#endif  // CELL_H