//#####################################################################
// 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
//#####################################################################
#include <Hybrid_Methods/Examples_And_Drivers/MPM_TRANSFERS.h>
#include <Core/Diff_Tools/Diff_Test.h>
#include <Core/Math_Tools/pow.h>
#include <Core/Random_Numbers/RANDOM_NUMBERS.h>
namespace PhysBAM{
//#####################################################################
// Function Transfer_P2G
//#####################################################################
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)
{
    for(auto it(it_);it.Valid();it.Next())
    {
        T w=it.Weight();
        TV V=Vp;
        auto idx=it.Index();
        if(Cp)
        {
            TV Y=grid.Node_X(idx)-Xp;
            V+=*Cp*Y;
            if(Hp) V+=Contract<2>(*Hp,Y)*Y*(T).5;
        }
        m(idx)+=w*mp;
        m_u(idx)+=w*mp*V;
    }
}
//#####################################################################
// Function Transfer_P2G_PH
//#####################################################################
template<class TV,class T> void
Transfer_P2G_PH(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)
{
    for(auto it(it_);it.Valid();it.Next())
    {
        T w=it.Weight();
        TV V=Vp;
        auto idx=it.Index();
        if(Cp)
        {
            TV Y=grid.Node_X(idx)-Xp;
            V+=*Cp*Y;
            if(Hp) V+=Contract<2>(*Hp,Y)*Y*(T).5;
        }
        m(idx)+=w*mp;
        m_u(idx)+=w*mp*V;
    }
}
//#####################################################################
// Function Transfer_G2P    (velocity)
//#####################################################################
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_,//v_save
    const ARRAYS_ND_BASE<TV,VECTOR<int,TV::m> >& Vi_u_,//v_friction_save
    const ARRAYS_ND_BASE<TV,VECTOR<int,TV::m> >& Vi_x_,//v
    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)
{
    TV V_x,V_n,V_u;
    MATRIX<T,TV::m> B;
    SYMMETRIC_TENSOR<T,0,TV::m> H;
    for(auto it(it_);it.Valid();it.Next())
    {
        T w=it.Weight();
        auto idx=it.Index();
        auto Vi_n=Vi_n_(idx);
        auto Vi_u=Vi_u_(idx);
        auto vi_x=Vi_x_(idx);
        V_x+=w*vi_x;
        V_n+=w*Vi_n;
        V_u+=w*Vi_u;
        TV V_grid=opt.use_midpoint?(T).5*(vi_x+Vi_n):vi_x;
        if(opt.use_affine)
        {
            TV Z=grid.Node_X(it.Index());
            if(opt.use_midpoint && !opt.use_polypic)
            {
                // TV V_grid=(T).5*(Vi_u+Vi_n);
                TV xi_new=Z+dt*V_grid;
                B+=w/2*(Outer_Product(Vi_u,Z+xi_new)+Outer_Product(Z-xi_new,Vi_u));
            }
            else B+=Outer_Product(w*Vi_u,Z-Xp_n);
            if(opt.use_polypic)
                H+=Tensor_Product<0>(Outer_Product(Z-Xp_n),w*Vi_u);
        }
    }
    if(opt.use_midpoint) V_x=(T).5*(V_n+V_x);
    TV xp_new=Xp_n+dt*V_x;
    Vp=V_u+opt.flip*(Vp_n-V_n);
    if(opt.use_polypic)
    {
        T d=1/Dp_inv;
        auto HH=H/(2*d*d);
        TV e=it_.scratch.m3,f=it_.scratch.m4;
        TV g=e*e-d*f+d*d*d;
        for(int a=0;a<TV::m;a++)
        {
            for(int b=0;b<TV::m;b++)
                HH(a,b,b)=(d*d*V_u(a)-d*H(a,b,b)+e(b)*B(a,b))/g(b);
            for(int b=0;b<TV::m;b++)
                B(a,b)=((d*d-f(b))*B(a,b)-d*e(b)*V_u(a)+e(b)*H(a,b,b))/g(b);
        }
        Vp-=d*Contract<1,2>(HH);
        *Cp=B;
        *Hp=HH*2;
    }
    else if(opt.use_affine)
    {
        if(opt.use_midpoint)
            B-=(T).5*(Outer_Product(V_u,Xp_n+xp_new)
                +Outer_Product(Xp_n-xp_new,V_u));
        *Cp=B*Dp_inv;
    }
}
//#####################################################################
// Function Transfer_G2P  (update position)
//#####################################################################
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)
{
    TV V_x,V_n;
    for(auto it(it_);it.Valid();it.Next())
    {
        T w=it.Weight();
        auto idx=it.Index();
        auto Vi_n=Vi_n_(idx);
        auto Vi_x=Vi_x_(idx);
        V_x+=w*Vi_x;
        V_n+=w*Vi_n;
    }
    if(opt.use_midpoint) V_x=(T).5*(V_n+V_x);
    Xp=Xp_n+dt*V_x;
}
//#####################################################################
// Function Transfer_G2P
//#####################################################################
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)
{
    MATRIX<T,TV::m> grad_Vp=Gather(opt,it_,Vi_x_);
    if(opt.use_midpoint)
        grad_Vp=(T).5*(grad_Vp+Gather(opt,it_,Vi_n_));
    return grad_Vp*dt;
}
//#####################################################################
// Function Transfer_G2P       (Defomration gradient)
//#####################################################################
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)
{
    if(opt.dilation_only)
    {
        T a=dt_grad_u.Trace(),b=a+1;
        if(opt.use_quad_F) b+=(T).5*a*a;
        T J=Fp.Determinant()*b;
        Fp=MATRIX<T,d>()+pow<1,d>(J);
    }
    else
    {
        MATRIX<T,d> A=dt_grad_u+1;
        if(opt.use_S)
        {
            T k=dt*opt.inv_Wi;
            *Sp=(Conjugate(A,*Sp)+k)/(1+k);
        }
        if(opt.use_quad_F) A+=(T).5*sqr(dt_grad_u);
        Fp=A*Fp;
    }
}
//#####################################################################
// Function Transfer_G2P_PS
//#####################################################################
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)
{
    
    MATRIX<T,d> A=dt_grad_u+1;
    if(opt.use_S)
    {
        T k=dt*opt.inv_Wi;
        *Sp=(Conjugate(A,*Sp)+k)/(1+k);
    }
    if(opt.use_quad_F) A+=(T).5*sqr(dt_grad_u);
    Fp=A*Fp;
}
//#####################################################################
// Function Gather
//#####################################################################
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)
{
    MATRIX<T,TV::m> M;
    for(auto it(it_);it.Valid();it.Next())
        M+=Outer_Product(u(it.Index()),it.Gradient(opt.use_cpic));
    return M;
}
//#####################################################################
// Function Scatter
//#####################################################################
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)
{
    for(auto it(it_);it.Valid();it.Next())
        f(it.Index())+=K*it.Gradient(opt.use_cpic);
}
//#####################################################################
// Function Diff_Times_G2P_F
//#####################################################################
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)
{
    if(opt.dilation_only)
    {
        T dt_div_Vp=dt_grad_u.Trace(),ddiv_Vp=grad_du.Trace();
        T a=dt_div_Vp,b=a+1,da=ddiv_Vp,db=da;
        if(opt.use_quad_F)
        {
            b+=(T).5*a*a;
            db+=a*da;
        }
        T det=Fp.Determinant(),J=det*b,dJ=det*db;
        return MATRIX<T,d>()+pow<1,d>(J)*dJ/(d*J);
    }
    MATRIX<T,d> dA=grad_du;
    if(opt.use_quad_F)
        dA+=(T).5*(grad_du*dt_grad_u+dt_grad_u*grad_du);
    return dA*Fp;
}

//#####################################################################
// Function Times_Diff_G2P_F
//#####################################################################
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)
{
    MATRIX<T,d> K;
    if(opt.dilation_only)
    {
        T det=Fp.Determinant(),a=dt_grad_u.Trace(),b=a+1;
        if(opt.use_quad_F) b+=(T).5*a*a;
        T J=det*b,k=det*pow<1,d>(J)/(J*d)*dM.Trace();
        if(opt.use_quad_F) k*=a+1;
        K+=k;
    }
    else
    {
        K=dM.Times_Transpose(Fp);
        if(opt.use_quad_F)
        {
            auto C=(T).5*dt_grad_u;
            K+=C.Transpose_Times(K)+K.Times_Transpose(C);
        }
    }
    return K;
}

//#####################################################################
// Function Times_Hess_Times_G2P_F
//#####################################################################
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)
{
    MATRIX<T,d> dK;
    if(opt.dilation_only)
    {
        T det=Fp.Determinant(),a=dt_grad_u.Trace(),da=grad_du.Trace();
        T b=a+1,db=da;
        if(opt.use_quad_F)
        {
            db+=a*da;
            b+=(T).5*a*a;
        }
        T J=det*b,dJ=det*db;
        T k=det*pow<1,d>(J)/(J*d)*dM.Trace();
        T dk=k*(1-d)/(J*d)*dJ;
        if(opt.use_quad_F)
        {
            dk=dk*(a+1)+k*da;
            k*=a+1;
        }
        dK+=dk;
    }
    else
    {
        if(opt.use_quad_F)
        {
            auto K=dM.Times_Transpose(Fp);
            auto dC=(T).5*grad_du;
            dK=dC.Transpose_Times(K)+K.Times_Transpose(dC);
        }
    }
    return dK;
}

//#####################################################################
// Function Diff_Times_G2P_S
//#####################################################################
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)
{
    if(!opt.use_S) return {};
    PHYSBAM_ASSERT(!opt.dilation_only);

    T c=1/(1+dt*opt.inv_Wi);
    auto A=grad_du.Times_Transpose((dt_grad_u+1)*(c*Sp));
    return A.Twice_Symmetric_Part();
}

//#####################################################################
// Function Times_Diff_G2P_S
//#####################################################################
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)
{
    if(!opt.use_S) return {};
    PHYSBAM_ASSERT(!opt.dilation_only);

    T c=1/(1+dt*opt.inv_Wi);
    return dM.Twice_Symmetric_Part()*(dt_grad_u+1)*(c*Sp);
}

//#####################################################################
// Function Times_Hess_Times_G2P_S
//#####################################################################
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)
{
    if(!opt.use_S) return {};
    PHYSBAM_ASSERT(!opt.dilation_only);

    T c=1/(1+dt*opt.inv_Wi);
    return dM.Twice_Symmetric_Part()*grad_du*(c*Sp);
}

// Lower s if necessary so that |t*(a+t*b)|<=bound for all 0<=t<=s
template<class T>
void Enforce_Limit_Max(T& s,T bound,T a,T b)
{
    if(a<0){a=-a;b=-b;}
    T p=s*(a+s*b),d=a*a+4*b*bound;
    if(d>=0)
    {
        if(p>bound) s=2*bound/(a+sqrt(d));
    }
    else
    {
        if(p<-bound) s=-(a+sqrt(a*a-4*b*bound))/(2*b);
    }
    PHYSBAM_ASSERT(s>=0);
}

//#####################################################################
// Function Transfer_G2P_Limit
//#####################################################################
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)
{
    PHYSBAM_ASSERT(!opt.use_quad_F);
    PHYSBAM_ASSERT(!opt.use_midpoint);
    TV V_x,V_x_s;
    MATRIX<T,TV::m> grad_Vp,grad_Vp_s;
    for(auto it(it_);it.Valid();it.Next())
    {
        T w=it.Weight();
        auto dw=it.Gradient(opt.use_cpic);
        auto idx=it.Index();
        auto Vi_x=Vi_n_(idx);
        auto Vi_x_s=Vi_x_(idx)-Vi_x;
        V_x+=w*Vi_x;
        V_x_s+=w*Vi_x_s;
        grad_Vp+=Outer_Product(Vi_x,dw);
        grad_Vp_s+=Outer_Product(Vi_x_s,dw);
    }
    if(cfl)
        for(int i=0;i<TV::m;i++)
            Enforce_Limit_Max(s,cfl*grid.dX(i),dt*V_x(i),dt*V_x_s(i));
    if(cfl_F)
    {
        if(opt.dilation_only)
            Enforce_Limit_Max(s,cfl_F,dt*grad_Vp.Trace(),dt*grad_Vp_s.Trace());
        else
            for(int i=0;i<TV::m;i++)
                for(int j=0;j<TV::m;j++)
                    Enforce_Limit_Max(s,cfl_F,dt*grad_Vp(i,j),dt*grad_Vp_s(i,j));
    }
    return s;
}

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)
{
    typedef VECTOR<int,TV::m> IV;
    T eps=(T)1e-6;
    T mp=opt.use_midpoint?(T).5:1;

    ARRAY<TV,IV> dVi_x(Vi_x.domain),Vi_y(Vi_x.domain);
    rand.Fill_Uniform(dVi_x.array,-eps,eps);
    Vi_y.array=Vi_x.array+dVi_x.array;

    MATRIX<T,TV::m> Fp0,Fp1,dt_grad_u0,dt_grad_u1,M;
    rand.Fill_Uniform(Fp0,-1,1);
    Fp1=Fp0;
    auto Fn=Fp0;
    rand.Fill_Uniform(M,-1,1);
    SYMMETRIC_MATRIX<T,TV::m> Sp0,Sp1,S;
    rand.Fill_Uniform(Sp0,-1,1);
    rand.Fill_Uniform(S,-1,1);
    Sp1=Sp0;
    auto Sn=Sp0;
    dt_grad_u0=Compute_dt_Grad_U(grid,opt,it,Vi_n,Vi_x,dt);
    dt_grad_u1=Compute_dt_Grad_U(grid,opt,it,Vi_n,Vi_y,dt);
    Transfer_G2P(opt,dt,Fp0,&Sp0,dt_grad_u0);
    Transfer_G2P(opt,dt,Fp1,&Sp1,dt_grad_u1);
    
    auto K=Gather(opt,it,dVi_x);
    auto dFp0=Diff_Times_G2P_F(opt,Fn,dt_grad_u0,K)*dt*mp;
    auto dFp1=Diff_Times_G2P_F(opt,Fn,dt_grad_u1,K)*dt*mp;

    Diff_Times_Test("Fp",Fp0,dFp0,Vi_x.array,Fp1,dFp1,Vi_y.array);

    ARRAY<TV,IV> tu(Vi_x.domain),tu0(Vi_x.domain),tu1(Vi_x.domain);
    Scatter(opt,it,tu,Times_Diff_G2P_F(opt,Fn,M,dt_grad_u0));
    auto diff_times=M.Double_Contract(dFp0)/eps;
    auto times_diff=tu.array.Dot(dVi_x.array)/eps*dt*mp;
    LOG::printf("Fp vs Fp^T  %P %P -> %P\n",diff_times,times_diff,diff_times-times_diff);

    Scatter(opt,it,tu0,Times_Hess_Times_G2P_F(opt,Fn,dt*mp*dt*mp/eps*M,dt_grad_u0,K));
    Scatter(opt,it,tu1,Times_Hess_Times_G2P_F(opt,Fn,dt*mp*dt*mp/eps*M,dt_grad_u1,K));

    auto t0=M.Double_Contract(dFp0)/eps,t1=M.Double_Contract(dFp1)/eps;
    Diff_Test("dFp",t0,Vi_x.array,tu0.array,t1,Vi_y.array,tu1.array);

    if(opt.use_S)
    {
        auto dSp0=Diff_Times_G2P_S(opt,Sn,dt_grad_u0,K,dt)*dt*mp;
        auto dSp1=Diff_Times_G2P_S(opt,Sn,dt_grad_u1,K,dt)*dt*mp;

        Diff_Times_Test("Sp",Sp0,dSp0,Vi_x.array,Sp1,dSp1,Vi_y.array);

        ARRAY<TV,IV> tu(Vi_x.domain),tu0(Vi_x.domain),tu1(Vi_x.domain);
        Scatter(opt,it,tu,Times_Diff_G2P_S(opt,Sn,S,dt_grad_u0,dt));
        auto diff_times=S.Double_Contract(dSp0)/eps;
        auto times_diff=tu.array.Dot(dVi_x.array)/eps*dt*mp;
        LOG::printf("Sp vs Sp^T  %P %P -> %P\n",diff_times,times_diff,diff_times-times_diff);

        Scatter(opt,it,tu0,Times_Hess_Times_G2P_S(opt,Sn,dt*mp*dt*mp/eps*S,dt_grad_u0,K,dt));
        Scatter(opt,it,tu1,Times_Hess_Times_G2P_S(opt,Sn,dt*mp*dt*mp/eps*S,dt_grad_u1,K,dt));

        auto t0=S.Double_Contract(dSp0)/eps,t1=S.Double_Contract(dSp1)/eps;
        Diff_Test("dSp",t0,Vi_x.array,tu0.array,t1,Vi_y.array,tu1.array);
    }
}

template MATRIX<double,2,2> Gather<VECTOR<double,2>,double>(MPM_TRANSFER_OPTIONS<double> const&,PARTICLE_GRID_ITERATOR<VECTOR<double,2> > const&,ARRAYS_ND_BASE<VECTOR<double,2>,VECTOR<int,2> > const&);
template MATRIX<double,1,1> Diff_Times_G2P_F<double,1>(MPM_TRANSFER_OPTIONS<double> const&,MATRIX<double,1,1> const&,MATRIX<double,1,1> const&,MATRIX<double,1,1> const&);
template MATRIX<double,1,1> Times_Diff_G2P_F<double,1>(MPM_TRANSFER_OPTIONS<double> const&,MATRIX<double,1,1> const&,MATRIX<double,1,1> const&,MATRIX<double,1,1> const&);
template MATRIX<double,1,1> Times_Hess_Times_G2P_F<double,1>(MPM_TRANSFER_OPTIONS<double> const&,MATRIX<double,1,1> const&,MATRIX<double,1,1> const&,MATRIX<double,1,1> const&,MATRIX<double,1,1> const&);
template MATRIX<double,2,2> Diff_Times_G2P_F<double,2>(MPM_TRANSFER_OPTIONS<double> const&,MATRIX<double,2,2> const&,MATRIX<double,2,2> const&,MATRIX<double,2,2> const&);
template MATRIX<double,2,2> Times_Diff_G2P_F<double,2>(MPM_TRANSFER_OPTIONS<double> const&,MATRIX<double,2,2> const&,MATRIX<double,2,2> const&,MATRIX<double,2,2> const&);
template MATRIX<double,2,2> Times_Diff_G2P_S<double,2>(MPM_TRANSFER_OPTIONS<double> const&,SYMMETRIC_MATRIX<double,2> const&,SYMMETRIC_MATRIX<double,2> const&,MATRIX<double,2,2> const&,double);
template MATRIX<double,2,2> Times_Hess_Times_G2P_F<double,2>(MPM_TRANSFER_OPTIONS<double> const&,MATRIX<double,2,2> const&,MATRIX<double,2,2> const&,MATRIX<double,2,2> const&,MATRIX<double,2,2> const&);
template MATRIX<double,2,2> Times_Hess_Times_G2P_S<double,2>(MPM_TRANSFER_OPTIONS<double> const&,SYMMETRIC_MATRIX<double,2> const&,SYMMETRIC_MATRIX<double,2> const&,MATRIX<double,2,2> const&,MATRIX<double,2,2> const&,double);
template MATRIX<double,3,3> Diff_Times_G2P_F<double,3>(MPM_TRANSFER_OPTIONS<double> const&,MATRIX<double,3,3> const&,MATRIX<double,3,3> const&,MATRIX<double,3,3> const&);
template MATRIX<double,3,3> Times_Diff_G2P_F<double,3>(MPM_TRANSFER_OPTIONS<double> const&,MATRIX<double,3,3> const&,MATRIX<double,3,3> const&,MATRIX<double,3,3> const&);
template MATRIX<double,3,3> Times_Diff_G2P_S<double,3>(MPM_TRANSFER_OPTIONS<double> const&,SYMMETRIC_MATRIX<double,3> const&,SYMMETRIC_MATRIX<double,3> const&,MATRIX<double,3,3> const&,double);
template MATRIX<double,3,3> Times_Hess_Times_G2P_F<double,3>(MPM_TRANSFER_OPTIONS<double> const&,MATRIX<double,3,3> const&,MATRIX<double,3,3> const&,MATRIX<double,3,3> const&,MATRIX<double,3,3> const&);
template MATRIX<double,3,3> Times_Hess_Times_G2P_S<double,3>(MPM_TRANSFER_OPTIONS<double> const&,SYMMETRIC_MATRIX<double,3> const&,SYMMETRIC_MATRIX<double,3> const&,MATRIX<double,3,3> const&,MATRIX<double,3,3> const&,double);
template MATRIX<double,1,1> Gather<VECTOR<double,1>,double>(MPM_TRANSFER_OPTIONS<double> const&,PARTICLE_GRID_ITERATOR<VECTOR<double,1> > const&,ARRAYS_ND_BASE<VECTOR<double,1>,VECTOR<int,1> > const&);
template MATRIX<double,3,3> Gather<VECTOR<double,3>,double>(MPM_TRANSFER_OPTIONS<double> const&,PARTICLE_GRID_ITERATOR<VECTOR<double,3> > const&,ARRAYS_ND_BASE<VECTOR<double,3>,VECTOR<int,3> > const&);
template MATRIX<float,1,1> Diff_Times_G2P_F<float,1>(MPM_TRANSFER_OPTIONS<float> const&,MATRIX<float,1,1> const&,MATRIX<float,1,1> const&,MATRIX<float,1,1> const&);
template MATRIX<float,1,1> Times_Diff_G2P_F<float,1>(MPM_TRANSFER_OPTIONS<float> const&,MATRIX<float,1,1> const&,MATRIX<float,1,1> const&,MATRIX<float,1,1> const&);
template MATRIX<float,1,1> Times_Hess_Times_G2P_F<float,1>(MPM_TRANSFER_OPTIONS<float> const&,MATRIX<float,1,1> const&,MATRIX<float,1,1> const&,MATRIX<float,1,1> const&,MATRIX<float,1,1> const&);
template MATRIX<float,2,2> Diff_Times_G2P_F<float,2>(MPM_TRANSFER_OPTIONS<float> const&,MATRIX<float,2,2> const&,MATRIX<float,2,2> const&,MATRIX<float,2,2> const&);
template MATRIX<float,2,2> Times_Diff_G2P_F<float,2>(MPM_TRANSFER_OPTIONS<float> const&,MATRIX<float,2,2> const&,MATRIX<float,2,2> const&,MATRIX<float,2,2> const&);
template MATRIX<float,2,2> Times_Diff_G2P_S<float,2>(MPM_TRANSFER_OPTIONS<float> const&,SYMMETRIC_MATRIX<float,2> const&,SYMMETRIC_MATRIX<float,2> const&,MATRIX<float,2,2> const&,float);
template MATRIX<float,2,2> Times_Hess_Times_G2P_F<float,2>(MPM_TRANSFER_OPTIONS<float> const&,MATRIX<float,2,2> const&,MATRIX<float,2,2> const&,MATRIX<float,2,2> const&,MATRIX<float,2,2> const&);
template MATRIX<float,2,2> Times_Hess_Times_G2P_S<float,2>(MPM_TRANSFER_OPTIONS<float> const&,SYMMETRIC_MATRIX<float,2> const&,SYMMETRIC_MATRIX<float,2> const&,MATRIX<float,2,2> const&,MATRIX<float,2,2> const&,float);
template MATRIX<float,3,3> Diff_Times_G2P_F<float,3>(MPM_TRANSFER_OPTIONS<float> const&,MATRIX<float,3,3> const&,MATRIX<float,3,3> const&,MATRIX<float,3,3> const&);
template MATRIX<float,3,3> Times_Diff_G2P_F<float,3>(MPM_TRANSFER_OPTIONS<float> const&,MATRIX<float,3,3> const&,MATRIX<float,3,3> const&,MATRIX<float,3,3> const&);
template MATRIX<float,3,3> Times_Diff_G2P_S<float,3>(MPM_TRANSFER_OPTIONS<float> const&,SYMMETRIC_MATRIX<float,3> const&,SYMMETRIC_MATRIX<float,3> const&,MATRIX<float,3,3> const&,float);
template MATRIX<float,3,3> Times_Hess_Times_G2P_F<float,3>(MPM_TRANSFER_OPTIONS<float> const&,MATRIX<float,3,3> const&,MATRIX<float,3,3> const&,MATRIX<float,3,3> const&,MATRIX<float,3,3> const&);
template MATRIX<float,3,3> Times_Hess_Times_G2P_S<float,3>(MPM_TRANSFER_OPTIONS<float> const&,SYMMETRIC_MATRIX<float,3> const&,SYMMETRIC_MATRIX<float,3> const&,MATRIX<float,3,3> const&,MATRIX<float,3,3> const&,float);
template MATRIX<float,1,1> Gather<VECTOR<float,1>,float>(MPM_TRANSFER_OPTIONS<float> const&,PARTICLE_GRID_ITERATOR<VECTOR<float,1> > const&,ARRAYS_ND_BASE<VECTOR<float,1>,VECTOR<int,1> > const&);
template MATRIX<float,2,2> Gather<VECTOR<float,2>,float>(MPM_TRANSFER_OPTIONS<float> const&,PARTICLE_GRID_ITERATOR<VECTOR<float,2> > const&,ARRAYS_ND_BASE<VECTOR<float,2>,VECTOR<int,2> > const&);
template MATRIX<float,3,3> Gather<VECTOR<float,3>,float>(MPM_TRANSFER_OPTIONS<float> const&,PARTICLE_GRID_ITERATOR<VECTOR<float,3> > const&,ARRAYS_ND_BASE<VECTOR<float,3>,VECTOR<int,3> > const&);
template SYMMETRIC_MATRIX<double,2> Diff_Times_G2P_S<double,2>(MPM_TRANSFER_OPTIONS<double> const&,SYMMETRIC_MATRIX<double,2> const&,MATRIX<double,2,2> const&,MATRIX<double,2,2> const&,double);
template SYMMETRIC_MATRIX<double,3> Diff_Times_G2P_S<double,3>(MPM_TRANSFER_OPTIONS<double> const&,SYMMETRIC_MATRIX<double,3> const&,MATRIX<double,3,3> const&,MATRIX<double,3,3> const&,double);
template SYMMETRIC_MATRIX<float,2> Diff_Times_G2P_S<float,2>(MPM_TRANSFER_OPTIONS<float> const&,SYMMETRIC_MATRIX<float,2> const&,MATRIX<float,2,2> const&,MATRIX<float,2,2> const&,float);
template SYMMETRIC_MATRIX<float,3> Diff_Times_G2P_S<float,3>(MPM_TRANSFER_OPTIONS<float> const&,SYMMETRIC_MATRIX<float,3> const&,MATRIX<float,3,3> const&,MATRIX<float,3,3> const&,float);
template void Scatter<VECTOR<double,1>,double>(MPM_TRANSFER_OPTIONS<double> const&,PARTICLE_GRID_ITERATOR<VECTOR<double,1> > const&,ARRAYS_ND_BASE<VECTOR<double,1>,VECTOR<int,1> >&,MATRIX<double,1,1> const&);
template void Scatter<VECTOR<double,2>,double>(MPM_TRANSFER_OPTIONS<double> const&,PARTICLE_GRID_ITERATOR<VECTOR<double,2> > const&,ARRAYS_ND_BASE<VECTOR<double,2>,VECTOR<int,2> >&,MATRIX<double,2,2> const&);
template void Scatter<VECTOR<double,3>,double>(MPM_TRANSFER_OPTIONS<double> const&,PARTICLE_GRID_ITERATOR<VECTOR<double,3> > const&,ARRAYS_ND_BASE<VECTOR<double,3>,VECTOR<int,3> >&,MATRIX<double,3,3> const&);
template void Scatter<VECTOR<float,1>,float>(MPM_TRANSFER_OPTIONS<float> const&,PARTICLE_GRID_ITERATOR<VECTOR<float,1> > const&,ARRAYS_ND_BASE<VECTOR<float,1>,VECTOR<int,1> >&,MATRIX<float,1,1> const&);
template void Scatter<VECTOR<float,2>,float>(MPM_TRANSFER_OPTIONS<float> const&,PARTICLE_GRID_ITERATOR<VECTOR<float,2> > const&,ARRAYS_ND_BASE<VECTOR<float,2>,VECTOR<int,2> >&,MATRIX<float,2,2> const&);
template void Scatter<VECTOR<float,3>,float>(MPM_TRANSFER_OPTIONS<float> const&,PARTICLE_GRID_ITERATOR<VECTOR<float,3> > const&,ARRAYS_ND_BASE<VECTOR<float,3>,VECTOR<int,3> >&,MATRIX<float,3,3> const&);
template double Transfer_G2P_Limit<VECTOR<double,2>,double>(GRID<VECTOR<double,2> > const&,MPM_TRANSFER_OPTIONS<double> const&,PARTICLE_GRID_ITERATOR<VECTOR<double,2> > const&,ARRAYS_ND_BASE<VECTOR<double,2>,VECTOR<int,2> > const&,ARRAYS_ND_BASE<VECTOR<double,2>,VECTOR<int,2> > const&,double,double,double,double);
template double Transfer_G2P_Limit<VECTOR<double,3>,double>(GRID<VECTOR<double,3> > const&,MPM_TRANSFER_OPTIONS<double> const&,PARTICLE_GRID_ITERATOR<VECTOR<double,3> > const&,ARRAYS_ND_BASE<VECTOR<double,3>,VECTOR<int,3> > const&,ARRAYS_ND_BASE<VECTOR<double,3>,VECTOR<int,3> > const&,double,double,double,double);
template float Transfer_G2P_Limit<VECTOR<float,2>,float>(GRID<VECTOR<float,2> > const&,MPM_TRANSFER_OPTIONS<float> const&,PARTICLE_GRID_ITERATOR<VECTOR<float,2> > const&,ARRAYS_ND_BASE<VECTOR<float,2>,VECTOR<int,2> > const&,ARRAYS_ND_BASE<VECTOR<float,2>,VECTOR<int,2> > const&,float,float,float,float);
template float Transfer_G2P_Limit<VECTOR<float,3>,float>(GRID<VECTOR<float,3> > const&,MPM_TRANSFER_OPTIONS<float> const&,PARTICLE_GRID_ITERATOR<VECTOR<float,3> > const&,ARRAYS_ND_BASE<VECTOR<float,3>,VECTOR<int,3> > const&,ARRAYS_ND_BASE<VECTOR<float,3>,VECTOR<int,3> > const&,float,float,float,float);
template void Test_Diff_Transfers<VECTOR<double,2>,double>(RANDOM_NUMBERS<double>&,GRID<VECTOR<double,2> > const&,MPM_TRANSFER_OPTIONS<double> const&,PARTICLE_GRID_ITERATOR<VECTOR<double,2> > const&,ARRAYS_ND_BASE<VECTOR<double,2>,VECTOR<int,2> > const&,ARRAYS_ND_BASE<VECTOR<double,2>,VECTOR<int,2> > const&,double);
template void Test_Diff_Transfers<VECTOR<double,3>,double>(RANDOM_NUMBERS<double>&,GRID<VECTOR<double,3> > const&,MPM_TRANSFER_OPTIONS<double> const&,PARTICLE_GRID_ITERATOR<VECTOR<double,3> > const&,ARRAYS_ND_BASE<VECTOR<double,3>,VECTOR<int,3> > const&,ARRAYS_ND_BASE<VECTOR<double,3>,VECTOR<int,3> > const&,double);
template void Test_Diff_Transfers<VECTOR<float,2>,float>(RANDOM_NUMBERS<float>&,GRID<VECTOR<float,2> > const&,MPM_TRANSFER_OPTIONS<float> const&,PARTICLE_GRID_ITERATOR<VECTOR<float,2> > const&,ARRAYS_ND_BASE<VECTOR<float,2>,VECTOR<int,2> > const&,ARRAYS_ND_BASE<VECTOR<float,2>,VECTOR<int,2> > const&,float);
template void Test_Diff_Transfers<VECTOR<float,3>,float>(RANDOM_NUMBERS<float>&,GRID<VECTOR<float,3> > const&,MPM_TRANSFER_OPTIONS<float> const&,PARTICLE_GRID_ITERATOR<VECTOR<float,3> > const&,ARRAYS_ND_BASE<VECTOR<float,3>,VECTOR<int,3> > const&,ARRAYS_ND_BASE<VECTOR<float,3>,VECTOR<int,3> > const&,float);
template void Transfer_P2G<VECTOR<double,2>,double>(GRID<VECTOR<double,2> > const&,MPM_TRANSFER_OPTIONS<double> const&,PARTICLE_GRID_ITERATOR<VECTOR<double,2> > const&,ARRAYS_ND_BASE<double,VECTOR<int,2> >&,ARRAYS_ND_BASE<VECTOR<double,2>,VECTOR<int,2> >&,double,VECTOR<double,2> const&,VECTOR<double,2> const&,MATRIX<double,2,2> const*,SYMMETRIC_TENSOR<double,0,2,2> const*);
template void Transfer_P2G<VECTOR<double,3>,double>(GRID<VECTOR<double,3> > const&,MPM_TRANSFER_OPTIONS<double> const&,PARTICLE_GRID_ITERATOR<VECTOR<double,3> > const&,ARRAYS_ND_BASE<double,VECTOR<int,3> >&,ARRAYS_ND_BASE<VECTOR<double,3>,VECTOR<int,3> >&,double,VECTOR<double,3> const&,VECTOR<double,3> const&,MATRIX<double,3,3> const*,SYMMETRIC_TENSOR<double,0,3,3> const*);
template void Transfer_P2G<VECTOR<float,2>,float>(GRID<VECTOR<float,2> > const&,MPM_TRANSFER_OPTIONS<float> const&,PARTICLE_GRID_ITERATOR<VECTOR<float,2> > const&,ARRAYS_ND_BASE<float,VECTOR<int,2> >&,ARRAYS_ND_BASE<VECTOR<float,2>,VECTOR<int,2> >&,float,VECTOR<float,2> const&,VECTOR<float,2> const&,MATRIX<float,2,2> const*,SYMMETRIC_TENSOR<float,0,2,2> const*);
template void Transfer_P2G<VECTOR<float,3>,float>(GRID<VECTOR<float,3> > const&,MPM_TRANSFER_OPTIONS<float> const&,PARTICLE_GRID_ITERATOR<VECTOR<float,3> > const&,ARRAYS_ND_BASE<float,VECTOR<int,3> >&,ARRAYS_ND_BASE<VECTOR<float,3>,VECTOR<int,3> >&,float,VECTOR<float,3> const&,VECTOR<float,3> const&,MATRIX<float,3,3> const*,SYMMETRIC_TENSOR<float,0,3,3> const*);
template void Transfer_G2P<VECTOR<double,2>,double>(GRID<VECTOR<double,2> > const&,MPM_TRANSFER_OPTIONS<double> const&,PARTICLE_GRID_ITERATOR<VECTOR<double,2> > const&,ARRAYS_ND_BASE<VECTOR<double,2>,VECTOR<int,2> > const&,ARRAYS_ND_BASE<VECTOR<double,2>,VECTOR<int,2> > const&,VECTOR<double,2> const&,double,VECTOR<double,2>&);
template void Transfer_G2P<VECTOR<double,2>,double>(GRID<VECTOR<double,2> > const&,MPM_TRANSFER_OPTIONS<double> const&,PARTICLE_GRID_ITERATOR<VECTOR<double,2> > const&,ARRAYS_ND_BASE<VECTOR<double,2>,VECTOR<int,2> > const&,ARRAYS_ND_BASE<VECTOR<double,2>,VECTOR<int,2> > const&,ARRAYS_ND_BASE<VECTOR<double,2>,VECTOR<int,2> > const&,double,VECTOR<double,2> const&,VECTOR<double,2> const&,double,VECTOR<double,2>&,MATRIX<double,2,2>*,SYMMETRIC_TENSOR<double,0,2,2>*);
template void Transfer_G2P<VECTOR<double,3>,double>(GRID<VECTOR<double,3> > const&,MPM_TRANSFER_OPTIONS<double> const&,PARTICLE_GRID_ITERATOR<VECTOR<double,3> > const&,ARRAYS_ND_BASE<VECTOR<double,3>,VECTOR<int,3> > const&,ARRAYS_ND_BASE<VECTOR<double,3>,VECTOR<int,3> > const&,VECTOR<double,3> const&,double,VECTOR<double,3>&);
template void Transfer_G2P<VECTOR<double,3>,double>(GRID<VECTOR<double,3> > const&,MPM_TRANSFER_OPTIONS<double> const&,PARTICLE_GRID_ITERATOR<VECTOR<double,3> > const&,ARRAYS_ND_BASE<VECTOR<double,3>,VECTOR<int,3> > const&,ARRAYS_ND_BASE<VECTOR<double,3>,VECTOR<int,3> > const&,ARRAYS_ND_BASE<VECTOR<double,3>,VECTOR<int,3> > const&,double,VECTOR<double,3> const&,VECTOR<double,3> const&,double,VECTOR<double,3>&,MATRIX<double,3,3>*,SYMMETRIC_TENSOR<double,0,3,3>*);
template void Transfer_G2P<VECTOR<float,2>,float>(GRID<VECTOR<float,2> > const&,MPM_TRANSFER_OPTIONS<float> const&,PARTICLE_GRID_ITERATOR<VECTOR<float,2> > const&,ARRAYS_ND_BASE<VECTOR<float,2>,VECTOR<int,2> > const&,ARRAYS_ND_BASE<VECTOR<float,2>,VECTOR<int,2> > const&,VECTOR<float,2> const&,float,VECTOR<float,2>&);
template void Transfer_G2P<VECTOR<float,2>,float>(GRID<VECTOR<float,2> > const&,MPM_TRANSFER_OPTIONS<float> const&,PARTICLE_GRID_ITERATOR<VECTOR<float,2> > const&,ARRAYS_ND_BASE<VECTOR<float,2>,VECTOR<int,2> > const&,ARRAYS_ND_BASE<VECTOR<float,2>,VECTOR<int,2> > const&,ARRAYS_ND_BASE<VECTOR<float,2>,VECTOR<int,2> > const&,float,VECTOR<float,2> const&,VECTOR<float,2> const&,float,VECTOR<float,2>&,MATRIX<float,2,2>*,SYMMETRIC_TENSOR<float,0,2,2>*);
template void Transfer_G2P<VECTOR<float,3>,float>(GRID<VECTOR<float,3> > const&,MPM_TRANSFER_OPTIONS<float> const&,PARTICLE_GRID_ITERATOR<VECTOR<float,3> > const&,ARRAYS_ND_BASE<VECTOR<float,3>,VECTOR<int,3> > const&,ARRAYS_ND_BASE<VECTOR<float,3>,VECTOR<int,3> > const&,VECTOR<float,3> const&,float,VECTOR<float,3>&);
template void Transfer_G2P<VECTOR<float,3>,float>(GRID<VECTOR<float,3> > const&,MPM_TRANSFER_OPTIONS<float> const&,PARTICLE_GRID_ITERATOR<VECTOR<float,3> > const&,ARRAYS_ND_BASE<VECTOR<float,3>,VECTOR<int,3> > const&,ARRAYS_ND_BASE<VECTOR<float,3>,VECTOR<int,3> > const&,ARRAYS_ND_BASE<VECTOR<float,3>,VECTOR<int,3> > const&,float,VECTOR<float,3> const&,VECTOR<float,3> const&,float,VECTOR<float,3>&,MATRIX<float,3,3>*,SYMMETRIC_TENSOR<float,0,3,3>*);
template void Transfer_G2P(MPM_TRANSFER_OPTIONS<double> const&,double,MATRIX<double,2,2>&,SYMMETRIC_MATRIX<double,2>*,MATRIX<double,2,2> const&);
template void Transfer_G2P(MPM_TRANSFER_OPTIONS<double> const&,double,MATRIX<double,3,3>&,SYMMETRIC_MATRIX<double,3>*,MATRIX<double,3,3> const&);
template void Transfer_G2P(MPM_TRANSFER_OPTIONS<float> const&,float,MATRIX<float,2,2>&,SYMMETRIC_MATRIX<float,2>*,MATRIX<float,2,2> const&);
template void Transfer_G2P(MPM_TRANSFER_OPTIONS<float> const&,float,MATRIX<float,3,3>&,SYMMETRIC_MATRIX<float,3>*,MATRIX<float,3,3> const&);
template void Transfer_G2P_PS(MPM_TRANSFER_OPTIONS<double> const&,double,MATRIX<double,2,2>&,SYMMETRIC_MATRIX<double,2>*,MATRIX<double,2,2> const&);
template void Transfer_G2P_PS(MPM_TRANSFER_OPTIONS<double> const&,double,MATRIX<double,3,3>&,SYMMETRIC_MATRIX<double,3>*,MATRIX<double,3,3> const&);
template void Transfer_G2P_PS(MPM_TRANSFER_OPTIONS<float> const&,float,MATRIX<float,2,2>&,SYMMETRIC_MATRIX<float,2>*,MATRIX<float,2,2> const&);
template void Transfer_G2P_PS(MPM_TRANSFER_OPTIONS<float> const&,float,MATRIX<float,3,3>&,SYMMETRIC_MATRIX<float,3>*,MATRIX<float,3,3> const&);
template MATRIX<double,2,2> Compute_dt_Grad_U<VECTOR<double,2>,double>(GRID<VECTOR<double,2> > const&,MPM_TRANSFER_OPTIONS<double> const&,PARTICLE_GRID_ITERATOR<VECTOR<double,2> > const&,ARRAYS_ND_BASE<VECTOR<double,2>,VECTOR<int,2> > const&,ARRAYS_ND_BASE<VECTOR<double,2>,VECTOR<int,2> > const&,double);
template MATRIX<double,3,3> Compute_dt_Grad_U<VECTOR<double,3>,double>(GRID<VECTOR<double,3> > const&,MPM_TRANSFER_OPTIONS<double> const&,PARTICLE_GRID_ITERATOR<VECTOR<double,3> > const&,ARRAYS_ND_BASE<VECTOR<double,3>,VECTOR<int,3> > const&,ARRAYS_ND_BASE<VECTOR<double,3>,VECTOR<int,3> > const&,double);
template MATRIX<float,2,2> Compute_dt_Grad_U<VECTOR<float,2>,float>(GRID<VECTOR<float,2> > const&,MPM_TRANSFER_OPTIONS<float> const&,PARTICLE_GRID_ITERATOR<VECTOR<float,2> > const&,ARRAYS_ND_BASE<VECTOR<float,2>,VECTOR<int,2> > const&,ARRAYS_ND_BASE<VECTOR<float,2>,VECTOR<int,2> > const&,float);
template MATRIX<float,3,3> Compute_dt_Grad_U<VECTOR<float,3>,float>(GRID<VECTOR<float,3> > const&,MPM_TRANSFER_OPTIONS<float> const&,PARTICLE_GRID_ITERATOR<VECTOR<float,3> > const&,ARRAYS_ND_BASE<VECTOR<float,3>,VECTOR<int,3> > const&,ARRAYS_ND_BASE<VECTOR<float,3>,VECTOR<int,3> > const&,float);
}
