/home/dorival/mechsys/lib/lbm/Domain.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_DOMAIN_H
#define MECHSYS_LBM_DOMAIN_H

// STD
#include <map>
#include <vector>
#include <utility>
#include <set>

// MechSys
#include <mechsys/lbm/Interacton.h>

using std::set;
using std::map;
using std::pair;
using std::make_pair;

namespace LBM
{

class Domain
{
public:
    //typedefs
    typedef void (*ptDFun_t) (Domain & Dom, void * UserData);

    //Constructors
    Domain (LBMethod      Method, 
    Array<double>         nu,     
    iVec3_t               Ndim,   
    double                dx,     
    double                dt);    

    //Special constructor with only one component, the parameters are the same as above
    Domain (LBMethod      Method, 
    double                nu,     
    iVec3_t               Ndim,   
    double                dx,     
    double                dt);    


    //Methods
#ifdef USE_HDF5
    void WriteXDMF (char const * FileKey);           
#endif
    void ApplyForce (size_t n = 0, size_t Np = 1);   
    void Collide    (size_t n = 0, size_t Np = 1);   
    void ImprintLattice ();                          
    void Solve(double Tf, double dtOut, ptDFun_t ptSetup=NULL, ptDFun_t ptReport=NULL,
               char const * FileKey=NULL, bool RenderVideo=true, size_t Nproc=1);                                                          
    void AddDisk  (int TheTag, Vec3_t const & TheX, Vec3_t const & TheV, Vec3_t const & TheW, double Therho, double TheR, double dt);        
    void AddSphere(int TheTag, Vec3_t const & TheX, Vec3_t const & TheV, Vec3_t const & TheW, double Therho, double TheR, double dt);        
    void ResetContacts();                            
    void ResetDisplacements();                       
    double  MaxDisplacement();                       
    void BoundingBox(Vec3_t & Xmin, Vec3_t & Xmax);  



    //Data
    bool                                    PrtVec;         
    Array<Lattice>                             Lat;         
    Array <Particle *>                   Particles;         
    Array <Interacton *>               Interactons;         
    Array <Interacton *>              CInteractons;         
    Array <iVec3_t>                      CellPairs;         
    set<pair<Particle *, Particle *> > Listofpairs;         
    double                                    Time;         
    double                                      dt;         
    double                                   Alpha;         
    double                                    Gmix;         
    void *                                UserData;         
    size_t                                 idx_out;         
    
};

#ifdef USE_THREAD
struct MtData
{
    size_t   ProcRank; 
    size_t     N_Proc; 
    LBM::Domain * Dom; 
};

void * GlobalSetZeroGamma(void * Data)
{
    MtData & dat = (*static_cast<MtData *>(Data));
    for (size_t i=0;i<dat.Dom->Lat.Size();i++)
    {
        dat.Dom->Lat[i].SetZeroGamma(dat.ProcRank, dat.N_Proc);
    }
}

void * GlobalApplyForce (void * Data)
{
    MtData & dat = (*static_cast<MtData *>(Data));
    dat.Dom->ApplyForce(dat.ProcRank, dat.N_Proc);
}

void * GlobalCollide (void * Data)
{
    MtData & dat = (*static_cast<MtData *>(Data));
    dat.Dom->Collide(dat.ProcRank, dat.N_Proc);
}

void * GlobalBounceBack (void * Data)
{
    MtData & dat = (*static_cast<MtData *>(Data));
    for (size_t i=0;i<dat.Dom->Lat.Size();i++)
    {
        dat.Dom->Lat[i].BounceBack(dat.ProcRank, dat.N_Proc);
    }
}

void * GlobalStream (void * Data)
{
    MtData & dat = (*static_cast<MtData *>(Data));
    for (size_t i=0;i<dat.Dom->Lat.Size();i++)
    {
        dat.Dom->Lat[i].Stream(dat.ProcRank, dat.N_Proc);
    }
}

void * GlobalStream1 (void * Data)
{
    MtData & dat = (*static_cast<MtData *>(Data));
    for (size_t i=0;i<dat.Dom->Lat.Size();i++)
    {
        dat.Dom->Lat[i].Stream1(dat.ProcRank, dat.N_Proc);
    }
}

void * GlobalStream2 (void * Data)
{
    MtData & dat = (*static_cast<MtData *>(Data));
    for (size_t i=0;i<dat.Dom->Lat.Size();i++)
    {
        dat.Dom->Lat[i].Stream2(dat.ProcRank, dat.N_Proc);
    }
}

#endif

inline Domain::Domain(LBMethod Method, Array<double> nu, iVec3_t Ndim, double dx, double Thedt)
{
    if (nu.Size()==0) throw new Fatal("LBM::Domain: Declare at leat one fluid please");
    if (Ndim(2) >1&&Method==D2Q9)  throw new Fatal("LBM::Domain: D2Q9 scheme does not allow for a third dimension, please set Ndim(2)=1 or change to D3Q15");
    if (Ndim(2)==1&&Method==D3Q15) throw new Fatal("LBM::Domain: Ndim(2) is greater than 1. Either change the method to D2Q9 or increse the z-dimension");
    for (size_t i=0;i<nu.Size();i++)
    {
        Lat.Push(Lattice(Method,nu[i],Ndim,dx,Thedt));
    }
    Time   = 0.0;
    dt     = Thedt;
    Alpha  = 10.0;
    PrtVec = false;
}

inline Domain::Domain(LBMethod Method, double nu, iVec3_t Ndim, double dx, double Thedt)
{
    Lat.Push(Lattice(Method,nu,Ndim,dx,Thedt));
    if (Ndim(2) >1&&Method==D2Q9)  throw new Fatal("LBM::Domain: D2Q9 scheme does not allow for a third dimension, please set Ndim(2)=1 or change to D3Q15");
    if (Ndim(2)==1&&Method==D3Q15) throw new Fatal("LBM::Domain: Ndim(2) is greater than 1. Either change the method to D2Q9 or increse the z-dimension");
    Time   = 0.0;
    dt     = Thedt;
    Alpha  = 10.0;
    PrtVec = false;
}

#ifdef USE_HDF5
inline void Domain::WriteXDMF(char const * FileKey)
{
    String fn(FileKey);
    fn.append(".h5");
    hid_t     file_id;
    file_id = H5Fcreate(fn.CStr(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
    for (size_t j=0;j<Lat.Size();j++)
    {
        // Creating data sets
        float * Density   = new float[Lat[0].Ndim[0]*Lat[0].Ndim[1]*Lat[0].Ndim[2]];
        float * Gamma     = new float[Lat[0].Ndim[0]*Lat[0].Ndim[1]*Lat[0].Ndim[2]];
        //float * Velocity  = new float[Lat[0].Ndim[0]*Lat[0].Ndim[1]*Lat[0].Ndim[2]];
        //float * MassFlux  = new float[Lat[0].Ndim[0]*Lat[0].Ndim[1]*Lat[0].Ndim[2]];
        float * Vvec      = new float[3*Lat[0].Ndim[0]*Lat[0].Ndim[1]*Lat[0].Ndim[2]];
        for (size_t i=0;i<Lat[j].Cells.Size();i++)
        {
            double rho = Lat[j].Cells[i]->Rho;
            Vec3_t vel = Lat[j].Cells[i]->Vel;
            Density  [i] = (float) rho;
            Gamma    [i] = (float) Lat[j].Cells[i]->IsSolid? 1.0:Lat[j].Cells[i]->Gamma;
            //Velocity [i] = (float) norm(vel);
            //MassFlux [i] = (float) rho*norm(vel);
            Vvec[3*i  ]  = (float) vel(0);
            Vvec[3*i+1]  = (float) vel(1);
            Vvec[3*i+2]  = (float) vel(2);
        }

        //Write the data
        hsize_t dims[1];
        dims[0] = Lat[0].Ndim(0)*Lat[0].Ndim(1)*Lat[0].Ndim(2);
        String dsname;
        dsname.Printf("Density_%d",j);
        H5LTmake_dataset_float(file_id,dsname.CStr(),1,dims,Density );
        dsname.Printf("Gamma_%d",j);
        H5LTmake_dataset_float(file_id,dsname.CStr(),1,dims,Gamma   );
        //dsname.Printf("Velocity_%d",j);
        //H5LTmake_dataset_float(file_id,dsname.CStr(),1,dims,Velocity);
        //dsname.Printf("MassFlux_%d",j);
        //H5LTmake_dataset_float(file_id,dsname.CStr(),1,dims,MassFlux);
        dims[0] = 3*Lat[0].Ndim(0)*Lat[0].Ndim(1)*Lat[0].Ndim(2);
        dsname.Printf("Velocity_%d",j);
        H5LTmake_dataset_float(file_id,dsname.CStr(),1,dims,Vvec    );

        delete [] Density ;
        delete [] Gamma   ;
        //delete [] Velocity;
        //delete [] MassFlux;
        delete [] Vvec    ;
    }

    //Closing the file
    H5Fclose(file_id);

    // Writing xmf file
    std::ostringstream oss;

    if (Lat[0].Ndim(2)==1)
    {
        oss << "<?xml version=\"1.0\" ?>\n";
        oss << "<!DOCTYPE Xdmf SYSTEM \"Xdmf.dtd\" []>\n";
        oss << "<Xdmf Version=\"2.0\">\n";
        oss << " <Domain>\n";
        oss << "   <Grid Name=\"mesh1\" GridType=\"Uniform\">\n";
        oss << "     <Topology TopologyType=\"2DCoRectMesh\" Dimensions=\"" << Lat[0].Ndim(1) << " " << Lat[0].Ndim(0) << "\"/>\n";
        oss << "     <Geometry GeometryType=\"ORIGIN_DXDY\">\n";
        oss << "       <DataItem Format=\"XML\" NumberType=\"Float\" Dimensions=\"2\"> 0.0 0.0\n";
        oss << "       </DataItem>\n";
        oss << "       <DataItem Format=\"XML\" NumberType=\"Float\" Dimensions=\"2\"> 1.0 1.0\n";
        oss << "       </DataItem>\n";
        oss << "     </Geometry>\n";
        for (size_t j=0;j<Lat.Size();j++)
        {
        oss << "     <Attribute Name=\"Density_" << j << "\" AttributeType=\"Scalar\" Center=\"Node\">\n";
        oss << "       <DataItem Dimensions=\"" << Lat[0].Ndim(0) << " " << Lat[0].Ndim(1) << "\" NumberType=\"Float\" Precision=\"4\" Format=\"HDF\">\n";
        oss << "        " << fn.CStr() <<":/Density_" << j << "\n";
        oss << "       </DataItem>\n";
        oss << "     </Attribute>\n";
        oss << "     <Attribute Name=\"Gamma_" << j << "\" AttributeType=\"Scalar\" Center=\"Node\">\n";
        oss << "       <DataItem Dimensions=\"" << Lat[0].Ndim(0) << " " << Lat[0].Ndim(1) << "\" NumberType=\"Float\" Precision=\"4\" Format=\"HDF\">\n";
        oss << "        " << fn.CStr() <<":/Gamma_" << j << "\n";
        oss << "       </DataItem>\n";
        oss << "     </Attribute>\n";
        //oss << "     <Attribute Name=\"Velocity_" << j << "\" AttributeType=\"Scalar\" Center=\"Node\">\n";
        //oss << "       <DataItem Dimensions=\"" << Lat[0].Ndim(0) << " " << Lat[0].Ndim(1) << "\" NumberType=\"Float\" Precision=\"4\" Format=\"HDF\">\n";
        //oss << "        " << fn.CStr() <<":/Velocity_" << j << "\n";
        //oss << "       </DataItem>\n";
        //oss << "     </Attribute>\n";
        //oss << "     <Attribute Name=\"MassFlux_" << j << "\" AttributeType=\"Scalar\" Center=\"Node\">\n";
        //oss << "       <DataItem Dimensions=\"" << Lat[0].Ndim(0) << " " << Lat[0].Ndim(1) << "\" NumberType=\"Float\" Precision=\"4\" Format=\"HDF\">\n";
        //oss << "        " << fn.CStr() <<":/MassFlux_" << j << "\n";
        //oss << "       </DataItem>\n";
        //oss << "     </Attribute>\n";
        oss << "     <Attribute Name=\"Velocity_" << j << "\" AttributeType=\"Vector\" Center=\"Node\">\n";
        oss << "       <DataItem Dimensions=\"" << Lat[0].Ndim(0) << " " << Lat[0].Ndim(1) << " 3\" NumberType=\"Float\" Precision=\"4\" Format=\"HDF\">\n";
        oss << "        " << fn.CStr() <<":/Velocity_" << j << "\n";
        oss << "       </DataItem>\n";
        oss << "     </Attribute>\n";
        }
        oss << "   </Grid>\n";
        oss << " </Domain>\n";
        oss << "</Xdmf>\n";
    }
    else
    {
        oss << "<?xml version=\"1.0\" ?>\n";
        oss << "<!DOCTYPE Xdmf SYSTEM \"Xdmf.dtd\" []>\n";
        oss << "<Xdmf Version=\"2.0\">\n";
        oss << " <Domain>\n";
        oss << "   <Grid Name=\"mesh1\" GridType=\"Uniform\">\n";
        oss << "     <Topology TopologyType=\"3DCoRectMesh\" Dimensions=\"" << Lat[0].Ndim(2) << " " << Lat[0].Ndim(1) << " " << Lat[0].Ndim(0) << "\"/>\n";
        oss << "     <Geometry GeometryType=\"ORIGIN_DXDYDZ\">\n";
        oss << "       <DataItem Format=\"XML\" NumberType=\"Float\" Dimensions=\"3\"> 0.0 0.0 0.0\n";
        oss << "       </DataItem>\n";
        oss << "       <DataItem Format=\"XML\" NumberType=\"Float\" Dimensions=\"3\"> 1.0 1.0 1.0\n";
        oss << "       </DataItem>\n";
        oss << "     </Geometry>\n";
        for (size_t j=0;j<Lat.Size();j++)
        {
        oss << "     <Attribute Name=\"Density_" << j << "\" AttributeType=\"Scalar\" Center=\"Node\">\n";
        oss << "       <DataItem Dimensions=\"" << Lat[0].Ndim(0) << " " << Lat[0].Ndim(1) << " " << Lat[0].Ndim(2) << "\" NumberType=\"Float\" Precision=\"4\" Format=\"HDF\">\n";
        oss << "        " << fn.CStr() <<":/Density_" << j << "\n";
        oss << "       </DataItem>\n";
        oss << "     </Attribute>\n";
        oss << "     <Attribute Name=\"Gamma_" << j << "\" AttributeType=\"Scalar\" Center=\"Node\">\n";
        oss << "       <DataItem Dimensions=\"" << Lat[0].Ndim(0) << " " << Lat[0].Ndim(1) << " " << Lat[0].Ndim(2) << "\" NumberType=\"Float\" Precision=\"4\" Format=\"HDF\">\n";
        oss << "        " << fn.CStr() <<":/Gamma_" << j << "\n";
        oss << "       </DataItem>\n";
        oss << "     </Attribute>\n";
        //oss << "     <Attribute Name=\"Velocity_" << j << "\" AttributeType=\"Scalar\" Center=\"Node\">\n";
        //oss << "       <DataItem Dimensions=\"" << Lat[0].Ndim(0) << " " << Lat[0].Ndim(1) << " " << Lat[0].Ndim(2) << "\" NumberType=\"Float\" Precision=\"4\" Format=\"HDF\">\n";
        //oss << "        " << fn.CStr() <<":/Velocity_" << j << "\n";
        //oss << "       </DataItem>\n";
        //oss << "     </Attribute>\n";
        //oss << "     <Attribute Name=\"MassFlux_" << j << "\" AttributeType=\"Scalar\" Center=\"Node\">\n";
        //oss << "       <DataItem Dimensions=\"" << Lat[0].Ndim(0) << " " << Lat[0].Ndim(1) << " " << Lat[0].Ndim(2) << "\" NumberType=\"Float\" Precision=\"4\" Format=\"HDF\">\n";
        //oss << "        " << fn.CStr() <<":/MassFlux_" << j << "\n";
        //oss << "       </DataItem>\n";
        //oss << "     </Attribute>\n";
        oss << "     <Attribute Name=\"Velocity_" << j << "\" AttributeType=\"Vector\" Center=\"Node\">\n";
        oss << "       <DataItem Dimensions=\"" << Lat[0].Ndim(0) << " " << Lat[0].Ndim(1) << " " << Lat[0].Ndim(2) << " 3\" NumberType=\"Float\" Precision=\"4\" Format=\"HDF\">\n";
        oss << "        " << fn.CStr() <<":/Velocity_" << j << "\n";
        oss << "       </DataItem>\n";
        oss << "     </Attribute>\n";
        }
        oss << "   </Grid>\n";
        oss << " </Domain>\n";
        oss << "</Xdmf>\n";
    }
    fn = FileKey;
    fn.append(".xmf");
    std::ofstream of(fn.CStr(), std::ios::out);
    of << oss.str();
    of.close();
}
#endif

inline void Domain::ApplyForce(size_t n, size_t Np)
{
    size_t Ni = CellPairs.Size()/Np;
    size_t In = n*Ni;
    size_t Fn;
    n == Np-1 ? Fn = CellPairs.Size() : Fn = (n+1)*Ni;

    for (size_t i=In;i<Fn;i++)
    {
        size_t ind1 = CellPairs[i](0);
        size_t ind2 = CellPairs[i](1);
        size_t vec  = CellPairs[i](2);
        for (size_t j=0;j<Lat.Size();j++)
        {
            bool solid = false;
            Cell * c = Lat[j].Cells[ind1];
            double psi;
            if (c->IsSolid||c->Gamma>0.0)
            {
                psi   = 1.0;
                solid = true;
            }
            else fabs(Lat[j].G)>1.0e-12 ? psi = Lat[j].Psi(c->Rho) : psi = c->Rho;
            for (size_t k=0;k<Lat.Size();k++)
            {
                Cell * nb = Lat[k].Cells[ind2];
                double nb_psi;
                if (nb->IsSolid||nb->Gamma>0.0)
                {
                    nb_psi = 1.0;
                    solid  = true;
                }
                else fabs(Lat[j].G)>1.0e-12 ? nb_psi = Lat[k].Psi(nb->Rho) : nb_psi = nb->Rho;
                double G;
                solid ? G = Lat[j].Gs*2.0*ReducedValue(nb->Gs,c->Gs) : G = Lat[j].G; 
                Vec3_t BF(OrthoSys::O);
                if (j==k)
                {
                    BF += -G*psi*nb_psi*c->W[vec]*c->C[vec];
                }
                else if(!solid)
                {
                    BF += -Gmix*c->Rho*nb->Rho*c->W[vec]*c->C[vec];
                }
#ifdef USE_THREAD
                pthread_mutex_lock(&c ->lck);
                pthread_mutex_lock(&nb->lck);
#endif
                c ->BForce += BF;
                nb->BForce -= BF;
#ifdef USE_THREAD
                pthread_mutex_unlock(&c ->lck);
                pthread_mutex_unlock(&nb->lck);
#endif
            }
        }
    }
}

void Domain::Collide (size_t n, size_t Np)
{
    size_t Ni = Lat[0].Cells.Size()/Np;
    size_t In = n*Ni;
    size_t Fn;
    n == Np-1 ? Fn = Lat[0].Cells.Size() : Fn = (n+1)*Ni;
    for (size_t i=In;i<Fn;i++)
    {
        Vec3_t num(0.0,0.0,0.0);
        double den = 0.0;
        for (size_t j=0;j<Lat.Size();j++)
        {
            Cell * c = Lat[j].Cells[i];
            double tau = Lat[j].Tau;
            num += c->Vel*c->Rho/tau;
            den += c->Rho/tau;
        }
        Vec3_t Vmix = num/den;


        for (size_t j=0;j<Lat.Size();j++)
        {
            Cell * c = Lat[j].Cells[i];
            double rho = c->Rho;
            if (c->IsSolid||rho<1.0e-12) continue;
            if (fabs(c->Gamma-1.0)<1.0e-12&&fabs(Lat[j].G)>1.0e-12) continue;
            double Tau = Lat[j].Tau;
            Vec3_t DV  = Vmix + c->BForce*dt/rho;
            double Bn  = (c->Gamma*(Tau-0.5))/((1.0-c->Gamma)+(Tau-0.5));
            bool valid  = true;
            double alphal = 1.0;
            double alphat = 1.0;
            size_t num  = 0;
            while (valid)
            {
                valid = false;
                alphal  = alphat;
                for (size_t k=0;k<c->Nneigh;k++)
                {
                    double FDeqn = c->Feq(k,DV,rho);
                    c->Ftemp[k] = c->F[k] - alphal*((1 - Bn)*(c->F[k] - FDeqn)/Tau - Bn*c->Omeis[k]);
                    if (c->Ftemp[k]<0.0&&num<1)
                    {
                        double temp = fabs(c->F[k]/((1 - Bn)*(c->F[k] - FDeqn)/Tau - Bn*c->Omeis[k]));
                        if (temp<alphat) alphat = temp;
                        valid = true;
                    }
                }
                num++;
                if (num>2) 
                {
                    throw new Fatal("Domain::Collide: Redefine your time step, the current value ensures unstability");
                }
            }
            for (size_t k=0;k<c->Nneigh;k++)
            {
                if (isnan(c->Ftemp[k]))
                {
                    c->Gamma = 2.0;
                    #ifdef USE_HDF5
                    WriteXDMF("error");
                    #endif
                    std::cout << c->Density() << " " << c->BForce << " " << num << " " << alphat << " " << c->Index << " " << c->IsSolid << " " << j << " " << k << std::endl;
                    throw new Fatal("Domain::Collide: Body force gives nan value, check parameters");
                }
                c->F[k] = fabs(c->Ftemp[k]);
            }
        }
    }   
}

void Domain::ImprintLattice ()
{
    
    // 2D imprint
    if (Lat[0].Ndim(2)==1)
    {
        for (size_t i = 0;i<Particles.Size();i++)
        {
            LBM::Particle * Pa = Particles[i];
            for (size_t n=std::max(0.0,double(Pa->X(0)-Pa->R-Lat[0].dx)/Lat[0].dx);n<=std::min(double(Lat[0].Ndim(0)-1),double(Pa->X(0)+Pa->R+Lat[0].dx)/Lat[0].dx);n++)
            for (size_t m=std::max(0.0,double(Pa->X(1)-Pa->R-Lat[0].dx)/Lat[0].dx);m<=std::min(double(Lat[0].Ndim(1)-1),double(Pa->X(1)+Pa->R+Lat[0].dx)/Lat[0].dx);m++)
            {
                Cell  * cell = Lat[0].GetCell(iVec3_t(n,m,0));
                double x     = Lat[0].dx*cell->Index(0);
                double y     = Lat[0].dx*cell->Index(1);
                double z     = Lat[0].dx*cell->Index(2);
                Vec3_t  C(x,y,z);
                Array<Vec3_t> P(4);

                P[0] = C - 0.5*Lat[0].dx*OrthoSys::e0 - 0.5*Lat[0].dx*OrthoSys::e1;
                P[1] = C + 0.5*Lat[0].dx*OrthoSys::e0 - 0.5*Lat[0].dx*OrthoSys::e1;
                P[2] = C + 0.5*Lat[0].dx*OrthoSys::e0 + 0.5*Lat[0].dx*OrthoSys::e1;
                P[3] = C - 0.5*Lat[0].dx*OrthoSys::e0 + 0.5*Lat[0].dx*OrthoSys::e1;

                double len = 0.0;
                
                for (size_t j=0;j<4;j++)
                {
                    Vec3_t D = P[(j+1)%4] - P[j];
                    double a = dot(D,D);
                    double b = 2*dot(P[j]-Pa->X,D);
                    double c = dot(P[j]-Pa->X,P[j]-Pa->X) - Pa->R*Pa->R;
                    if (b*b-4*a*c>0.0)
                    {
                        double ta = (-b - sqrt(b*b-4*a*c))/(2*a);
                        double tb = (-b + sqrt(b*b-4*a*c))/(2*a);
                        if (ta>1.0&&tb>1.0) continue;
                        if (ta<0.0&&tb<0.0) continue;
                        if (ta<0.0) ta = 0.0;
                        if (tb>1.0) tb = 1.0;
                        len += norm((tb-ta)*D);
                    }
                }

                if (len>0.0)
                {
                    for (size_t j=0;j<Lat.Size();j++)
                    {
                        cell = Lat[j].GetCell(iVec3_t(n,m,0));
                        cell->Gamma   = std::max(len/(4.0*Lat[0].dx),cell->Gamma);
                        if (fabs(cell->Gamma-1.0)<1.0e-12&&fabs(Lat[0].G)>1.0e-12) continue;
                        Vec3_t B      = C - Pa->X;
                        Vec3_t VelP   = Pa->V + cross(Pa->W,B);
                        double rho = cell->Rho;
                        double Bn  = (cell->Gamma*(cell->Tau-0.5))/((1.0-cell->Gamma)+(cell->Tau-0.5));
                        for (size_t k=0;k<cell->Nneigh;k++)
                        {
                            double Fvpp    = cell->Feq(cell->Op[k],VelP,rho);
                            double Fvp     = cell->Feq(k          ,VelP,rho);
                            cell->Omeis[k] = cell->F[cell->Op[k]] - Fvpp - (cell->F[k] - Fvp);
                            Vec3_t Flbm    = -Bn*cell->Omeis[k]*cell->C[k];
                            Pa->F          += Flbm;
                            Pa->T          += cross(B,Flbm);
                            //std::cout << i << " " << Pa->F << " " << Flbm << " " << Bn << std::endl;
                        }
                    }
                }
            }
        }
    }

    //3D imprint
    else
    {
        for (size_t i = 0;i<Particles.Size();i++)
        {
            LBM::Particle * Pa = Particles[i];
            for (size_t n=std::max(0.0,double(Pa->X(0)-Pa->R-Lat[0].dx)/Lat[0].dx);n<=std::min(double(Lat[0].Ndim(0)-1),double(Pa->X(0)+Pa->R+Lat[0].dx)/Lat[0].dx);n++)
            for (size_t m=std::max(0.0,double(Pa->X(1)-Pa->R-Lat[0].dx)/Lat[0].dx);m<=std::min(double(Lat[0].Ndim(1)-1),double(Pa->X(1)+Pa->R+Lat[0].dx)/Lat[0].dx);m++)
            for (size_t l=std::max(0.0,double(Pa->X(2)-Pa->R-Lat[0].dx)/Lat[0].dx);l<=std::min(double(Lat[0].Ndim(2)-1),double(Pa->X(2)+Pa->R+Lat[0].dx)/Lat[0].dx);l++)
            {
                Cell  * cell = Lat[0].GetCell(iVec3_t(n,m,l));
                double x     = Lat[0].dx*cell->Index(0);
                double y     = Lat[0].dx*cell->Index(1);
                double z     = Lat[0].dx*cell->Index(2);
                Vec3_t  C(x,y,z);
                Array<Vec3_t> P(8);

                P[0] = C - 0.5*Lat[0].dx*OrthoSys::e0 - 0.5*Lat[0].dx*OrthoSys::e1 + 0.5*Lat[0].dx*OrthoSys::e2; 
                P[1] = C + 0.5*Lat[0].dx*OrthoSys::e0 - 0.5*Lat[0].dx*OrthoSys::e1 + 0.5*Lat[0].dx*OrthoSys::e2;
                P[2] = C + 0.5*Lat[0].dx*OrthoSys::e0 + 0.5*Lat[0].dx*OrthoSys::e1 + 0.5*Lat[0].dx*OrthoSys::e2;
                P[3] = C - 0.5*Lat[0].dx*OrthoSys::e0 + 0.5*Lat[0].dx*OrthoSys::e1 + 0.5*Lat[0].dx*OrthoSys::e2;
                P[4] = C - 0.5*Lat[0].dx*OrthoSys::e0 - 0.5*Lat[0].dx*OrthoSys::e1 - 0.5*Lat[0].dx*OrthoSys::e2; 
                P[5] = C + 0.5*Lat[0].dx*OrthoSys::e0 - 0.5*Lat[0].dx*OrthoSys::e1 - 0.5*Lat[0].dx*OrthoSys::e2;
                P[6] = C + 0.5*Lat[0].dx*OrthoSys::e0 + 0.5*Lat[0].dx*OrthoSys::e1 - 0.5*Lat[0].dx*OrthoSys::e2;
                P[7] = C - 0.5*Lat[0].dx*OrthoSys::e0 + 0.5*Lat[0].dx*OrthoSys::e1 - 0.5*Lat[0].dx*OrthoSys::e2;

                double len = 0.0;
                
                for (size_t j=0;j<4;j++)
                {
                    Vec3_t D;
                    double a; 
                    double b; 
                    double c; 
                    D = P[(j+1)%4] - P[j];
                    a = dot(D,D);
                    b = 2*dot(P[j]-Pa->X,D);
                    c = dot(P[j]-Pa->X,P[j]-Pa->X) - Pa->R*Pa->R;
                    if (b*b-4*a*c>0.0)
                    {
                        double ta = (-b - sqrt(b*b-4*a*c))/(2*a);
                        double tb = (-b + sqrt(b*b-4*a*c))/(2*a);
                        if (ta>1.0&&tb>1.0) continue;
                        if (ta<0.0&&tb<0.0) continue;
                        if (ta<0.0) ta = 0.0;
                        if (tb>1.0) tb = 1.0;
                        len += norm((tb-ta)*D);
                    }
                    D = P[(j+1)%4 + 4] - P[j + 4];
                    a = dot(D,D);
                    b = 2*dot(P[j + 4]-Pa->X,D);
                    c = dot(P[j + 4]-Pa->X,P[j + 4]-Pa->X) - Pa->R*Pa->R;
                    if (b*b-4*a*c>0.0)
                    {
                        double ta = (-b - sqrt(b*b-4*a*c))/(2*a);
                        double tb = (-b + sqrt(b*b-4*a*c))/(2*a);
                        if (ta>1.0&&tb>1.0) continue;
                        if (ta<0.0&&tb<0.0) continue;
                        if (ta<0.0) ta = 0.0;
                        if (tb>1.0) tb = 1.0;
                        len += norm((tb-ta)*D);
                    }
                    D = P[j+4] - P[j];
                    a = dot(D,D);
                    b = 2*dot(P[j]-Pa->X,D);
                    c = dot(P[j]-Pa->X,P[j]-Pa->X) - Pa->R*Pa->R;
                    if (b*b-4*a*c>0.0)
                    {
                        double ta = (-b - sqrt(b*b-4*a*c))/(2*a);
                        double tb = (-b + sqrt(b*b-4*a*c))/(2*a);
                        if (ta>1.0&&tb>1.0) continue;
                        if (ta<0.0&&tb<0.0) continue;
                        if (ta<0.0) ta = 0.0;
                        if (tb>1.0) tb = 1.0;
                        len += norm((tb-ta)*D);
                    }
                }

                if (len>0.0)
                {
                    for (size_t j=0;j<Lat.Size();j++)
                    {
                        cell = Lat[j].GetCell(iVec3_t(n,m,l));
                        cell->Gamma   = std::max(len/(12.0*Lat[0].dx),cell->Gamma);
                        if (fabs(cell->Gamma-1.0)<1.0e-12&&fabs(Lat[0].G)>1.0e-12) continue;
                        Vec3_t B      = C - Pa->X;
                        Vec3_t VelP   = Pa->V + cross(Pa->W,B);
                        double rho = cell->Rho;
                        double Bn  = (cell->Gamma*(cell->Tau-0.5))/((1.0-cell->Gamma)+(cell->Tau-0.5));
                        for (size_t k=0;k<cell->Nneigh;k++)
                        {
                            double Fvpp    = cell->Feq(cell->Op[k],VelP,rho);
                            double Fvp     = cell->Feq(k          ,VelP,rho);
                            cell->Omeis[k] = cell->F[cell->Op[k]] - Fvpp - (cell->F[k] - Fvp);
                            Vec3_t Flbm    = -Bn*cell->Omeis[k]*cell->C[k];
                            Pa->F          += Flbm;
                            Pa->T          += cross(B,Flbm);
                        }
                    }
                }
            }
        }
    }
}

inline void Domain::ResetDisplacements()
{
    for (size_t i=0; i<Particles.Size(); i++)
    {
        Particles[i]->X0 = Particles[i]->X;
    }
}

inline double Domain::MaxDisplacement()
{
    double md = 0.0;
    for (size_t i=0; i<Particles.Size(); i++)
    {
        double mdp = norm(Particles[i]->X - Particles[i]->X0);
        if (mdp>md) md = mdp;
    }
    return md;
}

inline void Domain::ResetContacts()
{
    if (Particles.Size()==0) return;
    for (size_t i=0; i<Particles.Size()-1; i++)
    {
        bool pi_has_vf = !Particles[i]->IsFree();
        for (size_t j=i+1; j<Particles.Size(); j++)
        {
            bool pj_has_vf = !Particles[j]->IsFree();

            bool close = (norm(Particles[i]->X-Particles[j]->X)<=Particles[i]->R+Particles[j]->R+2*Alpha);
            if ((pi_has_vf && pj_has_vf) || !close) continue;
            
            // checking if the interacton exist for that pair of particles
            set<pair<Particle *, Particle *> >::iterator it = Listofpairs.find(make_pair(Particles[i],Particles[j]));
            if (it != Listofpairs.end())
            {
                continue;
            }
            Listofpairs.insert(make_pair(Particles[i],Particles[j]));
            CInteractons.Push (new Interacton(Particles[i],Particles[j]));
        }
    }

    Interactons.Resize(0);
    for (size_t i=0; i<CInteractons.Size(); i++)
    {
        if(CInteractons[i]->UpdateContacts(Alpha)) Interactons.Push(CInteractons[i]);
    }
}

inline void Domain::BoundingBox(Vec3_t & minX, Vec3_t & maxX)
{
    minX = Vec3_t(Particles[0]->X(0) - Particles[0]->R, Particles[0]->X(1) - Particles[0]->R, Particles[0]->X(2) - Particles[0]->R);
    maxX = Vec3_t(Particles[0]->X(0) + Particles[0]->R, Particles[0]->X(1) + Particles[0]->R, Particles[0]->X(2) + Particles[0]->R);
    for (size_t i=1; i<Particles.Size(); i++)
    {
        if (minX(0)>(Particles[i]->X(0) - Particles[i]->R)&&Particles[i]->IsFree()) minX(0) = Particles[i]->X(0) - Particles[i]->R;
        if (minX(1)>(Particles[i]->X(1) - Particles[i]->R)&&Particles[i]->IsFree()) minX(1) = Particles[i]->X(1) - Particles[i]->R;
        if (minX(2)>(Particles[i]->X(2) - Particles[i]->R)&&Particles[i]->IsFree()) minX(2) = Particles[i]->X(2) - Particles[i]->R;
        if (maxX(0)<(Particles[i]->X(0) + Particles[i]->R)&&Particles[i]->IsFree()) maxX(0) = Particles[i]->X(0) + Particles[i]->R;
        if (maxX(1)<(Particles[i]->X(1) + Particles[i]->R)&&Particles[i]->IsFree()) maxX(1) = Particles[i]->X(1) + Particles[i]->R;
        if (maxX(2)<(Particles[i]->X(2) + Particles[i]->R)&&Particles[i]->IsFree()) maxX(2) = Particles[i]->X(2) + Particles[i]->R;
    }
}

inline void Domain::AddDisk(int TheTag, Vec3_t const & TheX, Vec3_t const & TheV, Vec3_t const & TheW, double Therho, double TheR, double dt)
{
    Particles.Push(new Disk(TheTag,TheX,TheV,TheW,Therho,TheR,dt));
}

inline void Domain::AddSphere(int TheTag, Vec3_t const & TheX, Vec3_t const & TheV, Vec3_t const & TheW, double Therho, double TheR, double dt)
{
    Particles.Push(new Sphere(TheTag,TheX,TheV,TheW,Therho,TheR,dt));
}

inline void Domain::Solve(double Tf, double dtOut, ptDFun_t ptSetup, ptDFun_t ptReport,
                          char const * FileKey, bool RenderVideo, size_t Nproc)
{
    // info
    Util::Stopwatch stopwatch;
    printf("\n%s--- Solving ---------------------------------------------------------------------%s\n",TERM_CLR1,TERM_RST);

    idx_out     = 0;
    //Connect particles and lattice
    //
    
    //for(size_t i=0;i<Particles.Size();i++) Particles[i]->ImprintDisk(Lat[0]);
    ImprintLattice();
    ResetContacts();
    ResetDisplacements();
    
    // Creates pair of cells to speed up body force calculation
    for (size_t i=0;i<Lat[0].Cells.Size();i++)
    {
        Cell * c = Lat[0].Cells[i];
        for (size_t j=1;j<c->Nneigh;j++)
        {
            Cell * nb = Lat[0].Cells[c->Neighs[j]];
            if (nb->ID>c->ID) 
            {
                if (!c->IsSolid||!nb->IsSolid) CellPairs.Push(iVec3_t(i,nb->ID,j));
            }
        }
    }
#ifdef USE_THREAD
    MtData MTD[Nproc];
    for (size_t i=0;i<Nproc;i++)
    {
        MTD[i].N_Proc   = Nproc;
        MTD[i].ProcRank = i;
        MTD[i].Dom      = this;
    }
    pthread_t thrs[Nproc];   
#endif
    double tout = Time;
    while (Time < Tf)
    {
        if (ptSetup!=NULL) (*ptSetup) ((*this), UserData);
        if (Time >= tout)
        {
            if (FileKey!=NULL)
            {
                String fn;
                fn.Printf    ("%s_%04d", FileKey, idx_out);
                if ( RenderVideo) 
                {
                    #ifdef USE_HDF5
                    WriteXDMF(fn.CStr());
                    #else
                    //WriteVTK (fn.CStr());
                    #endif
                }
                if (ptReport!=NULL) (*ptReport) ((*this), UserData);
            }
            tout += dtOut;
            idx_out++;
        }

        //Assigning a vlaue of zero to the particles forces and torques
        for(size_t i=0;i<Particles.Size();i++)
        {
            Particles[i]->F = Particles[i]->Ff;
            Particles[i]->T = Particles[i]->Tf;
        }

#ifdef USE_THREAD
        for (size_t i=0;i<Nproc;i++)
        {
            pthread_create(&thrs[i], NULL, GlobalSetZeroGamma, &MTD[i]);
        }
        for (size_t i=0;i<Nproc;i++)
        {
            pthread_join(thrs[i], NULL);
        }
#else 
        //Set Gamma values of the lattice cell to zero
        for(size_t j=0;j<Lat.Size();j++)
        {
            Lat[j].SetZeroGamma();
        }
#endif
        //Connect particles and lattice
        ImprintLattice();

        //Move Particles
        for(size_t i=0;i<Interactons.Size();i++) Interactons[i]->CalcForce(dt);
        for(size_t i=0;i<Particles.Size()  ;i++) Particles[i]->Translate(dt);

        //Move fluid
#ifdef USE_THREAD
        if (Lat.Size()>1||fabs(Lat[0].G)>1.0e-12)
        {
            for (size_t i=0;i<Nproc;i++)
            {
                pthread_create(&thrs[i], NULL, GlobalApplyForce, &MTD[i]);
            }
            for (size_t i=0;i<Nproc;i++)
            {
                pthread_join(thrs[i], NULL);
            }
        }
        for (size_t i=0;i<Nproc;i++)
        {
            pthread_create(&thrs[i], NULL, GlobalCollide, &MTD[i]);
        }
        for (size_t i=0;i<Nproc;i++)
        {
            pthread_join(thrs[i], NULL);
        }
        for (size_t i=0;i<Nproc;i++)
        {
            pthread_create(&thrs[i], NULL, GlobalBounceBack, &MTD[i]);
        }
        for (size_t i=0;i<Nproc;i++)
        {
            pthread_join(thrs[i], NULL);
        }
        for (size_t i=0;i<Nproc;i++)
        {
            pthread_create(&thrs[i], NULL, GlobalStream1, &MTD[i]);
        }
        for (size_t i=0;i<Nproc;i++)
        {
            pthread_join(thrs[i], NULL);
        }
        for (size_t i=0;i<Nproc;i++)
        {
            pthread_create(&thrs[i], NULL, GlobalStream2, &MTD[i]);
        }
        for (size_t i=0;i<Nproc;i++)
        {
            pthread_join(thrs[i], NULL);
        }

#else 
        if (Lat.Size()>1||fabs(Lat[0].G)>1.0e-12) ApplyForce();
        Collide();
        for(size_t j=0;j<Lat.Size();j++)
        {
            Lat[j].BounceBack();
            Lat[j].Stream();
        }
#endif

        if (MaxDisplacement()>Alpha)
        {
            ResetContacts();
            ResetDisplacements();
        }

        Time += dt;
    }
    printf("%s  Final CPU time       = %s\n",TERM_CLR2, TERM_RST);
}
}


#endif