//#####################################################################
// Copyright 2010.
// This file is part of PhysBAM whose distribution is governed by the license contained in the accompanying file PHYSBAM_COPYRIGHT.txt.
//#####################################################################
#include <Core/Grids_Uniform_Computations/LEVELSET_MAKER_UNIFORM.h>
#include <Core/Basic_Geometry/BASIC_SIMPLEX_POLICY.h>
#include <Core/Basic_Geometry/SEGMENT_2D.h>
#include <Core/Grids/CELL_ITERATOR.h>
#include <Core/Grids/GRID.h>
#include <Core/Level_Sets/FAST_MARCHING_METHOD_UNIFORM.h>
#include <Core/Level_Sets/LEVELSET.h>
#include <Core/Log/DEBUG_SUBSTEPS.h>
#include <Core/Math_Tools/RANGE.h>
#include <Core/Meshes/SEGMENTED_CURVE.h>
#include <Core/Meshes/TRIANGULATED_SURFACE.h>
using namespace PhysBAM;
//#####################################################################
// Function Compute_Level_Set
//#####################################################################
template<class TV> void LEVELSET_MAKER_UNIFORM<TV>::
Compute_Level_Set(T_SURFACE& surface,GRID<TV>& grid,AT_CELLS,int ghost_cells,ARRAY<T,IV>& phi)
{
    PHYSBAM_ASSERT(phi.domain.min_corner.Max()<=-ghost_cells);
    if(!surface.mesh.adjacent_elements) surface.mesh.Initialize_Adjacent_Elements();
    phi.Fill(FLT_MAX);
    T dx=grid.dX.Max();

    ARRAY<IV> seed_indices;
    for(int i=0;i<surface.mesh.elements.m;i++){
        typename BASIC_SIMPLEX_POLICY<TV,TV::m-1>::SIMPLEX simplex(surface.particles.X.Subset(surface.mesh.elements(i)));
        RANGE<IV>  box(grid.Clamp_To_Cell(simplex.Bounding_Box(),ghost_cells).Thickened(1).Intersect(phi.domain));
        for(RANGE_ITERATOR<TV::m> it(box);it.Valid();it.Next()){
            TV X=grid.Cell_X(it.index);
            T dist=simplex.Distance_To_Element(X),&p=phi(it.index);
            if(dist<abs(p)+dx*(T)1e-4 && dist<dx){
                if(p==FLT_MAX) seed_indices.Append(it.index);
                bool new_sign=simplex.Signed_Distance(X)<0;
                if(abs(dist-abs(p))<dx*(T)1e-4 && new_sign!=(p<0))
                    new_sign=surface.Inside(X);
                if(abs(dist)<abs(p)) p=dist;
                p=abs(p)*(new_sign?-1:1);}}}

    ARRAY<IV> todo;
    for(int i=0;i<seed_indices.m;i++)
        if(phi(seed_indices(i))<0)
            todo.Append(seed_indices(i));
    while(todo.m){
        ARRAY<IV> next;
        for(int i=0;i<todo.m;i++)
            for(int a=0;a<TV::m;a++)
                for(int s=-1;s<2;s+=2){
                    IV index=todo(i);
                    index(a)+=s;
                    if(phi.domain.Lazy_Inside_Half_Open(index) && phi(index)==FLT_MAX){
                        phi(index)=-FLT_MAX;
                        next.Append(index);}}
        next.Exchange(todo);}

    LEVELSET<TV> levelset(grid,at_cells,phi);
    FAST_MARCHING_METHOD_UNIFORM<TV> fmm;
    fmm.seed_indices=&seed_indices;
    fmm.correct_interface_phi=false;
    fmm.Fast_Marching_Method(grid,at_cells,ghost_cells,phi,0);
}
namespace PhysBAM{
template class LEVELSET_MAKER_UNIFORM<VECTOR<float,2> >;
template class LEVELSET_MAKER_UNIFORM<VECTOR<float,3> >;
template class LEVELSET_MAKER_UNIFORM<VECTOR<double,2> >;
template class LEVELSET_MAKER_UNIFORM<VECTOR<double,3> >;
}
