//#####################################################################
// Copyright 2015, Craig Schroeder.
// This file is part of PhysBAM whose distribution is governed by the license contained in the accompanying file PHYSBAM_COPYRIGHT.txt.
//#####################################################################
// Class MPM_TRANSFERS
//#####################################################################
#ifndef __MPM_TRANSFERS__
#define __MPM_TRANSFERS__

#include <Core/Grids/GRID.h>
#include <Core/Matrices/MATRIX.h>
#include <Core/Matrices/SYMMETRIC_MATRIX.h>
#include <Core/Tensors/SYMMETRIC_TENSOR.h>
#include <Hybrid_Methods/Examples_And_Drivers/MPM_TRANSFER_OPTIONS.h>
#include <Hybrid_Methods/Iterators/PARTICLE_GRID_ITERATOR.h>
namespace PhysBAM{
template<class T> std::pair<bool,bool>
Transfer_P2G(const MPM_TRANSFER_OPTIONS<T>& opt){return {};}

template<class TV,class T> void
Transfer_P2G(const GRID<TV>& grid,const MPM_TRANSFER_OPTIONS<T>& opt,
    const PARTICLE_GRID_ITERATOR<TV>& it_,
    ARRAYS_ND_BASE<T,VECTOR<int,TV::m> >& m,
    ARRAYS_ND_BASE<TV,VECTOR<int,TV::m> >& m_u,
    T mp,const TV& Vp,const TV& Xp,const MATRIX<ND<T>,TV::m>* Cp,
    const SYMMETRIC_TENSOR<ND<T>,0,TV::m>* Hp);

template<class T> std::pair<bool,bool>
Transfer_G2P(const MPM_TRANSFER_OPTIONS<T>& opt){return {!opt.use_cpic,opt.use_cpic};}

template<class TV,class T> void
Transfer_G2P(const GRID<TV>& grid,const MPM_TRANSFER_OPTIONS<T>& opt,
    const PARTICLE_GRID_ITERATOR<TV>& it_,
    const ARRAYS_ND_BASE<TV,VECTOR<int,TV::m> >& Vi_n_,
    const ARRAYS_ND_BASE<TV,VECTOR<int,TV::m> >& Vi_u_,
    const ARRAYS_ND_BASE<TV,VECTOR<int,TV::m> >& Vi_x_,
    T Dp_inv,const TV& Xp_n,const TV& Vp_n,T dt,
    TV& Vp,MATRIX<ND<T>,TV::m>* Cp,SYMMETRIC_TENSOR<ND<T>,0,TV::m>* Hp);

template<class TV,class T> void
Transfer_G2P(const GRID<TV>& grid,const MPM_TRANSFER_OPTIONS<T>& opt,
    const PARTICLE_GRID_ITERATOR<TV>& it_,
    const ARRAYS_ND_BASE<TV,VECTOR<int,TV::m> >& Vi_n_,
    const ARRAYS_ND_BASE<TV,VECTOR<int,TV::m> >& Vi_x_,
    const TV& Xp_n,T dt,TV& Xp);

template<class TV,class T> MATRIX<T,TV::m>
Compute_dt_Grad_U(const GRID<TV>& grid,const MPM_TRANSFER_OPTIONS<T>& opt,
    const PARTICLE_GRID_ITERATOR<TV>& it_,
    const ARRAYS_ND_BASE<TV,VECTOR<int,TV::m> >& Vi_n_,
    const ARRAYS_ND_BASE<TV,VECTOR<int,TV::m> >& Vi_x_,T dt);

template<class T,int d> void
Transfer_G2P(const MPM_TRANSFER_OPTIONS<T>& opt,T dt,MATRIX<T,d>& Fp,
    SYMMETRIC_MATRIX<ND<T>,d+0>* Sp,const MATRIX<ND<T>,d>& dt_grad_u);

template<class T,int d> void
Transfer_G2P_PS(const MPM_TRANSFER_OPTIONS<T>& opt,T dt,MATRIX<T,d>& Fp,
    SYMMETRIC_MATRIX<ND<T>,d+0>* Sp,const MATRIX<ND<T>,d>& dt_grad_u);

template<class TV,class T> MATRIX<T,TV::m>
Gather(const MPM_TRANSFER_OPTIONS<T>& opt,
    const PARTICLE_GRID_ITERATOR<TV>& it_,
    const ARRAYS_ND_BASE<TV,VECTOR<int,TV::m> >& u);

template<class TV,class T> void
Scatter(const MPM_TRANSFER_OPTIONS<T>& opt,
    const PARTICLE_GRID_ITERATOR<TV>& it_,
    ARRAYS_ND_BASE<TV,VECTOR<int,TV::m> >& f,
    const MATRIX<T,TV::m>& K);

// G=Gather(opt,it,z);
// D=Diff_Times_G2P_F(opt,Fn,dt_grad_u,G);
// A=D*dt/(1+use_midpoint);
// Computes: A=sum_i diff(Fp,Vi_x(i)) * z(i)
// The sum is over grid indicies.
// z is an arbitrary grid vector field.
// dt_grad_u is obtained from Transfer_G2P.
// Fn is F_p^n (input to Transfer_G2P, not output).
template<class T,int d> MATRIX<T,d>
Diff_Times_G2P_F(const MPM_TRANSFER_OPTIONS<T>& opt,
    const MATRIX<T,d>& Fp,const MATRIX<T,d>& dt_grad_u,
    const MATRIX<T,d>& grad_du);

// D=Times_Diff_G2P_F(opt,Fn,dM,dt_grad_u);
// Scatter(opt,it,y,D);
// z=y*dt/(1+use_midpoint);
// Computes: z(i)=dM : diff(Fp,Vi_x(i))
// z is an output grid vector field.
// dt_grad_u is obtained from Transfer_G2P.
// Fn is F_p^n (input to Transfer_G2P, not output).
template<class T,int d> MATRIX<T,d>
Times_Diff_G2P_F(const MPM_TRANSFER_OPTIONS<T>& opt,
    const MATRIX<T,d>& Fp,
    const MATRIX<T,d>& dM,const MATRIX<T,d>& dt_grad_u);

// G=Gather(opt,it,z);
// D=Times_Hess_Times_G2P_F(opt,Fn,dM,dt_grad_u,G);
// Scatter(opt,it,y,D);
// w=y*sqr(dt/(1+use_midpoint));
// Computes: w(i)=sum_i dM : diff(diff(Fp,Vi_x(i)),Vi_x(i)) * z(i)
// z is an input grid vector field.
// w is an output grid vector field.
// dt_grad_u is obtained from Transfer_G2P.
// Fn is F_p^n (input to Transfer_G2P, not output).
template<class T,int d> MATRIX<T,d>
Times_Hess_Times_G2P_F(const MPM_TRANSFER_OPTIONS<T>& opt,
    const MATRIX<T,d>& Fp,
    const MATRIX<T,d>& dM,const MATRIX<T,d>& dt_grad_u,
    const MATRIX<T,d>& grad_du);

// G=Gather(opt,it,z);
// D=Diff_Times_G2P_S(opt,Sn,dt_grad_u,G);
// A=D*dt/(1+use_midpoint);
// Computes: A=sum_i diff(Sp,Vi_x(i)) * z(i)
// The sum is over grid indicies.
// z is an arbitrary grid vector field.
// dt_grad_u is obtained from Transfer_G2P.
// Sn is S_p^n (input to Transfer_G2P, not output).
template<class T,int d> SYMMETRIC_MATRIX<T,d>
Diff_Times_G2P_S(const MPM_TRANSFER_OPTIONS<T>& opt,
    const SYMMETRIC_MATRIX<T,d>& Sp,const MATRIX<T,d>& dt_grad_u,const MATRIX<T,d>& grad_du,T dt);

// D=Times_Diff_G2P_S(opt,Sn,dM,dt_grad_u);
// Scatter(opt,it,y,D);
// z=y*dt/(1+use_midpoint);
// Computes: z(i)=dM : diff(Sp,Vi_x(i))
// z is an output grid vector field.
// dt_grad_u is obtained from Transfer_G2P.
// Sn is S_p^n (input to Transfer_G2P, not output).
template<class T,int d> MATRIX<T,d>
Times_Diff_G2P_S(const MPM_TRANSFER_OPTIONS<T>& opt,
    const SYMMETRIC_MATRIX<T,d>& Sp,
    const SYMMETRIC_MATRIX<T,d>& dM,const MATRIX<T,d>& dt_grad_u,T dt);

// G=Gather(opt,it,z);
// D=Times_Hess_Times_G2P_S(opt,Sn,dM,dt_grad_u,G);
// Scatter(opt,it,y,D);
// w=y*sqr(dt/(1+use_midpoint));
// Computes: w(i)=sum_i dM : diff(diff(Sp,Vi_x(i)),Vi_x(i)) * z(i)
// z is an input grid vector field.
// w is an output grid vector field.
// dt_grad_u is obtained from Transfer_G2P.
// Sn is S_p^n (input to Transfer_G2P, not output).
template<class T,int d> MATRIX<T,d>
Times_Hess_Times_G2P_S(const MPM_TRANSFER_OPTIONS<T>& opt,
    const SYMMETRIC_MATRIX<T,d>& Sp,
    const SYMMETRIC_MATRIX<T,d>& dM,const MATRIX<T,d>& dt_grad_u,
    const MATRIX<T,d>& grad_du,T dt);

template<class TV,class T> T
Transfer_G2P_Limit(const GRID<TV>& grid,const MPM_TRANSFER_OPTIONS<T>& opt,
    const PARTICLE_GRID_ITERATOR<TV>& it_,
    const ARRAYS_ND_BASE<TV,VECTOR<int,TV::m> >& Vi_n_,
    const ARRAYS_ND_BASE<TV,VECTOR<int,TV::m> >& Vi_x_,T dt,T s,T cfl,T cfl_F);

template<class TV,class T> void
Test_Diff_Transfers(RANDOM_NUMBERS<T>& rand,const GRID<TV>& grid,
    const MPM_TRANSFER_OPTIONS<T>& opt,
    const PARTICLE_GRID_ITERATOR<TV>& it_,
    const ARRAYS_ND_BASE<TV,VECTOR<int,TV::m> >& Vi_n,
    const ARRAYS_ND_BASE<TV,VECTOR<int,TV::m> >& Vi_x,T dt);

}
#endif
