/home/dorival/mechsys/lib/lbm/Cell.h

Go to the documentation of this file.

/************************************************************************
 * MechSys - Open Library for Mechanical Systems                        *
 * Copyright (C) 2005 Dorival M. Pedroso, Raul Durand                   *
 * Copyright (C) 2009 Sergio Galindo                                    *
 *                                                                      *
 * This program is free software: you can redistribute it and/or modify *
 * it under the terms of the GNU General Public License as published by *
 * the Free Software Foundation, either version 3 of the License, or    *
 * any later version.                                                   *
 *                                                                      *
 * This program 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 General Public License for more details.                         *
 *                                                                      *
 * You should have received a copy of the GNU General Public License    *
 * along with this program. If not, see <http://www.gnu.org/licenses/>  *
 ************************************************************************/

#ifndef MECHSYS_LBM_CELL_H
#define MECHSYS_LBM_CELL_H

// Std lib
#ifdef USE_THREAD
    #include <pthread.h>
#endif

// MechSys
#include <mechsys/util/array.h>
#include <mechsys/util/fatal.h>
#include <mechsys/linalg/matvec.h>

enum LBMethod
{
    D2Q9,     
    D3Q15,    
    D3Q27     
};

class Cell
{
public:
    static const double   WEIGHTSD2Q9   [ 9]; 
    static const double   WEIGHTSD3Q15  [15]; 
    //static const double   WEIGHTSD3Q27  [27]; ///< Weights for the equilibrium distribution functions (D3Q27)
    static Vec3_t         LVELOCD2Q9    [ 9]; 
    static Vec3_t         LVELOCD3Q15   [15]; 
    //static const Vec3_t   LVELOCD3Q27   [27]; ///< Local velocities (D3Q27)
    static const size_t   OPPOSITED2Q9  [ 9]; 
    static const size_t   OPPOSITED3Q15 [15]; 
    //static const size_t   OPPOSITED3Q27 [27]; ///< Opposite directions (D3Q27)
   
    //Constructor
    Cell (size_t ID, LBMethod Method, iVec3_t Indexes, iVec3_t Ndim, double Cs, double Tau); 
    
    // Methods
    double       Density();                                        
    void         Velocity(Vec3_t & V);                             
    double       VelDen(Vec3_t & V);                               
    double       Feq(size_t k, Vec3_t const & V, double Rho);      
    void         Initialize(double Rho, Vec3_t const & V);         
    void         BounceBack();                                     

#ifdef USE_THREAD
    pthread_mutex_t lck;
    //std::mutex mtex;       ///< to protect variables in multithreading
#endif

    // Data
    LBMethod     Method;   
    bool         IsSolid;  
    double       Tau;      
    double       Gamma;    
    double       Gs;       
    size_t       Nneigh;   
    size_t       ID;       
    double       Cs;       
    double       Rho;      
    iVec3_t      Index;    
    Vec3_t       Vel;      
    Vec3_t       VelBC;    
    Vec3_t       BForce;   
    Vec3_t       BForcef;  
    double       RhoBC;    

    size_t const  * Op;    
    double const  * W;     
    Vec3_t const  * C;     
    Array<double>   F;     
    Array<double>   Ftemp; 
    Array<double>   Omeis; 
    Array<size_t>   Neighs;
};

inline Cell::Cell(size_t TheID, LBMethod TheMethod, iVec3_t TheIndexes, iVec3_t TheNdim, double TheCs, double TheTau)
{
    IsSolid= false;
    ID     = TheID;
    Method = TheMethod;
    Index  = TheIndexes;
    Cs     = TheCs;
    Gamma  = 0.0;
    Gs     = 1.0;
    Tau    = TheTau;
    if (Method==D2Q9)
    {
        Nneigh = 9;
        W      = WEIGHTSD2Q9;
        C      = LVELOCD2Q9;
        Op     = OPPOSITED2Q9;
    }
    if (Method==D3Q15)
    {
        Nneigh = 15;
        W      = WEIGHTSD3Q15;
        C      = LVELOCD3Q15;
        Op     = OPPOSITED3Q15;
    }
    F.     Resize(Nneigh);
    Ftemp. Resize(Nneigh);
    Neighs.Resize(Nneigh);
    Omeis .Resize(Nneigh);
    BForcef = 0.0,0.0,0.0;

    //Set neighbors
    for (size_t k=0;k<Nneigh;k++)
    {
        //iVec3_t nindex = Index + C[k];
        blitz::TinyVector<int,3>   nindex = Index + C[k];
        if (nindex[0]==        -1) nindex[0] = TheNdim[0]-1;
        if (nindex[0]==TheNdim[0]) nindex[0] = 0;
        if (nindex[1]==        -1) nindex[1] = TheNdim[1]-1;
        if (nindex[1]==TheNdim[1]) nindex[1] = 0;
        if (nindex[2]==        -1) nindex[2] = TheNdim[2]-1;
        if (nindex[2]==TheNdim[2]) nindex[2] = 0;

        Neighs[k] =  nindex[0] + nindex[1]*TheNdim[0] + nindex[2]*TheNdim[0]*TheNdim[1];
    }
#ifdef USE_THREAD
    pthread_mutex_init(&lck,NULL);
#endif
}

inline double Cell::Density()
{
    if (IsSolid) return 0.0;
    double rho = 0.0;
    for (size_t k=0;k<Nneigh;k++) rho += F[k];
    return rho;
}

inline void Cell::Velocity(Vec3_t & V)
{
    V = 0.0, 0.0, 0.0;
    if (IsSolid) return;
    double rho = Density();
    if (rho<1.0e-12) return;
    for (size_t k=0;k<Nneigh;k++) V += F[k]*C[k];
    V *= Cs/rho;
}

inline double Cell::VelDen(Vec3_t & V)
{
    V = 0.0, 0.0, 0.0;
    if (IsSolid) return 0.0;
    double rho = Density();
    if (rho<1.0e-12) return 0.0;
    for (size_t k=1;k<Nneigh;k++) V += F[k]*C[k];
    V *= Cs/rho;
    return rho;
}

inline double Cell::Feq(size_t k, Vec3_t const & V, double TheRho)
{
    double VdotC = dot(V,C[k]);
    double VdotV = dot(V,V);
    return W[k]*TheRho*(1.0 + 3.0*VdotC/Cs + 4.5*VdotC*VdotC/(Cs*Cs) - 1.5*VdotV/(Cs*Cs));
}

inline void Cell::Initialize(double TheRho, Vec3_t const & V)
{
    for (size_t k=0;k<Nneigh;k++) F[k] = Feq(k,V,TheRho);
    Rho = VelDen(Vel);
}

inline void Cell::BounceBack()
{
    for (size_t j = 0;j<Nneigh;j++) Ftemp[j] = F[j];
    for (size_t j = 0;j<Nneigh;j++) F[j]     = Ftemp[Op[j]];
}


const double Cell::WEIGHTSD2Q9   [ 9] = { 4./9., 1./9., 1./9., 1./9., 1./9., 1./36., 1./36., 1./36., 1./36. };
const double Cell::WEIGHTSD3Q15  [15] = { 2./9., 1./9., 1./9., 1./9., 1./9.,  1./9.,  1./9., 1./72., 1./72. , 1./72., 1./72., 1./72., 1./72., 1./72., 1./72.};
const size_t Cell::OPPOSITED2Q9  [ 9] = { 0, 3, 4, 1, 2, 7, 8, 5, 6 };                       
const size_t Cell::OPPOSITED3Q15 [15] = { 0, 2, 1, 4, 3, 6, 5, 8, 7, 10, 9, 12, 11, 14, 13}; 
Vec3_t Cell::LVELOCD2Q9  [ 9]; 
Vec3_t Cell::LVELOCD3Q15 [15]; 


/*TODO: Go back to the previous version since this is not nice
 *      Changes made for g++ 4.3.2 */
int __allocate_constants ()
{
    Cell::LVELOCD2Q9  [0] =  0 ,  0 , 0;
    Cell::LVELOCD2Q9  [1] =  1 ,  0 , 0;
    Cell::LVELOCD2Q9  [2] =  0 ,  1 , 0;
    Cell::LVELOCD2Q9  [3] = -1 ,  0 , 0;
    Cell::LVELOCD2Q9  [4] =  0 , -1 , 0;
    Cell::LVELOCD2Q9  [5] =  1 ,  1 , 0;
    Cell::LVELOCD2Q9  [6] = -1 ,  1 , 0;
    Cell::LVELOCD2Q9  [7] = -1 , -1 , 0;
    Cell::LVELOCD2Q9  [8] =  1 , -1 , 0;

    Cell::LVELOCD3Q15[ 0]= 0, 0, 0;  Cell::LVELOCD3Q15[ 1]= 1, 0, 0;  Cell::LVELOCD3Q15[ 2]=-1, 0, 0;  Cell::LVELOCD3Q15[ 3]= 0, 1, 0;  Cell::LVELOCD3Q15[ 4]= 0,-1, 0;
    Cell::LVELOCD3Q15[ 5]= 0, 0, 1;  Cell::LVELOCD3Q15[ 6]= 0, 0,-1;  Cell::LVELOCD3Q15[ 7]= 1, 1, 1;  Cell::LVELOCD3Q15[ 8]=-1,-1,-1;  Cell::LVELOCD3Q15[ 9]= 1, 1,-1; 
    Cell::LVELOCD3Q15[10]=-1,-1, 1;  Cell::LVELOCD3Q15[11]= 1,-1, 1;  Cell::LVELOCD3Q15[12]=-1, 1,-1;  Cell::LVELOCD3Q15[13]= 1,-1,-1;  Cell::LVELOCD3Q15[14]=-1, 1, 1;

    return 0;
}
int __dummy__allocate_constants = __allocate_constants();

#endif // MECHSYS_LBM_CELL_H