//#####################################################################
// Copyright 2015-2020, Craig Schroeder.
// This file is part of PhysBAM whose distribution is governed by the license contained in the accompanying file PHYSBAM_COPYRIGHT.txt.
//#####################################################################
#include <Hybrid_Methods/Examples_And_Drivers/MPM_HELPERS.h>
#include <Core/Grids/GRID.h>
#include <Hybrid_Methods/Examples_And_Drivers/MPM_PARTICLES.h>
namespace PhysBAM{
//#####################################################################
// Constructor
//#####################################################################
template<class TV> PARTICLE_BOUNDARY_CONDITIONS<TV>::
PARTICLE_BOUNDARY_CONDITIONS(MPM_PARTICLES<TV>& particles,const GRID<TV>& grid,
    const VECTOR<MPM::BC_TYPE,TV::m*2>& side_bc_type)
    :side_bc_type(side_bc_type),range(grid.domain),grid_domain(grid.domain),
    particles(particles)
{
    // avoid being on the edge and rasterizing to a grid cell outside the domain.
    TV eps=std::numeric_limits<T>::epsilon()*4*range.Edge_Lengths();
    for(int i=0;i<TV::m;i++)
        if(side_bc_type(i)!=MPM::BC_PERIODIC)
        {
            range.min_corner(i)+=eps(i);
            range.max_corner(i)-=eps(i);
        }
    periodic_shift=range.Edge_Lengths();
}
//#####################################################################
// Function Move_Particle
//#####################################################################
template<class TV> void PARTICLE_BOUNDARY_CONDITIONS<TV>::
Enforce(int p)
{
    TV& X=particles.X(p);
    // TODO: particles.C?
    if(range.Lazy_Inside(X)) return;
    bool has_noslip=false;
    bool only_periodic=true;
    auto* C=particles.store_C?&particles.C(p):nullptr;
//    C*(X1-X0)=V1-V0;

    TV Y=X; // Closest feasible point
    TV& V=particles.V(p);

    VECTOR<MPM::BC_TYPE,2*TV::m> actual_bc_type=side_bc_type;
    for(int a=0;a<TV::m;a++)
    {
        if(X(a)<range.min_corner(a))
        {
            if(actual_bc_type(2*a)==MPM::BC_MIXED)
                actual_bc_type(2*a)=bc_mixed_type(2*a)(grid_domain.Clamp(X));
            switch(actual_bc_type(2*a))
            {
                case MPM::BC_FREE:
                    particles.valid(p)=false;
                    return;
                case MPM::BC_PERIODIC:
                    X(a)=wrap(X(a),grid_domain.min_corner(a),grid_domain.max_corner(a));
                    Y(a)=X(a);
                    break;
                case MPM::BC_NOSLIP:
                    has_noslip=true;
                default:
                    only_periodic=false;
                    Y(a)=range.min_corner(a);
                    break;
            }
        }
        else if(X(a)>range.max_corner(a))
        {
            if(actual_bc_type(2*a+1)==MPM::BC_MIXED)
                actual_bc_type(2*a+1)=bc_mixed_type(2*a+1)(grid_domain.Clamp(X));
            switch(actual_bc_type(2*a+1))
            {
                case MPM::BC_FREE:
                    particles.valid(p)=false;
                    return;
                case MPM::BC_PERIODIC:
                    X(a)=wrap(X(a),grid_domain.min_corner(a),grid_domain.max_corner(a));
                    Y(a)=X(a);
                    break;
                case MPM::BC_NOSLIP:
                    has_noslip=true;
                default:
                    only_periodic=false;
                    Y(a)=range.max_corner(a);
                    break;
            }
        }
    }
    if(only_periodic) return;

    TV U;
    if(bc_v) U=bc_v(Y);
    if(has_noslip)
    {
        X=Y;
        V=2*U-V;
        // C is unchanged
        return;
    }
    // Only BC_SEP and BC_SLIP remain.

    for(int a=0;a<TV::m;a++)
    {
        if(X(a)<range.min_corner(a))
        {
            if(C) C->Set_Column(a,-C->Column(a));
            if(actual_bc_type(2*a)==MPM::BC_SEP && V(a)>=U(a)) continue;
            V(a)=2*U(a)-V(a);
            if(C) C->Set_Row(a,-C->Row(a));
        }
        else if(X(a)>range.max_corner(a))
        {
            if(C) C->Set_Column(a,-C->Column(a));
            if(actual_bc_type(2*a+1)==MPM::BC_SEP && V(a)<=U(a)) continue;
            V(a)=2*U(a)-V(a);
            if(C) C->Set_Row(a,-C->Row(a));
        }
    }
    X=Y;
}
template class PARTICLE_BOUNDARY_CONDITIONS<VECTOR<double,2> >;
template class PARTICLE_BOUNDARY_CONDITIONS<VECTOR<double,3> >;
template class PARTICLE_BOUNDARY_CONDITIONS<VECTOR<float,2> >;
template class PARTICLE_BOUNDARY_CONDITIONS<VECTOR<float,3> >;
}
