//#####################################################################
// Copyright 2002-2006, Ronald Fedkiw, Sergey Koltakov, Eran Guendelman, Geoffrey Irving, Neil Molino, Andrew Selle, Rachel Weinstein.
// This file is part of PhysBAM whose distribution is governed by the license contained in the accompanying file PHYSBAM_COPYRIGHT.txt.
//#####################################################################
#include <Core/Implicit_Objects/IMPLICIT_OBJECT.h>
#include <Core/Basic_Geometry/BOUNDED_HORIZONTAL_PLANE.h>
#include <Core/Basic_Geometry/BOWL.h>
#include <Core/Basic_Geometry/CYLINDER.h>
#include <Core/Basic_Geometry/PLANE.h>
#include <Core/Basic_Geometry/RANGE_TESSELLATION.h>
#include <Core/Basic_Geometry/RAY.h>
#include <Core/Basic_Geometry/RING.h>
#include <Core/Basic_Geometry/SMOOTH_GEAR.h>
#include <Core/Basic_Geometry/SPHERE.h>
#include <Core/Basic_Geometry/TORUS.h>
#include <Core/Diff_Tools/Diff_Test.h>
#include <Core/Grids_Uniform_Computations/MARCHING_CUBES.h>
#include <Core/Implicit_Objects/ANALYTIC_IMPLICIT_OBJECT.h>
#include <Core/Implicit_Objects/IMPLICIT_OBJECT_COMBINED.h>
#include <Core/Implicit_Objects/IMPLICIT_OBJECT_TRANSFORMED.h>
#include <Core/Implicit_Objects/LEVELSET_IMPLICIT_OBJECT.h>
#include <Core/Implicit_Objects/MULTIBODY_LEVELSET_IMPLICIT_OBJECT.h>
#include <Core/Level_Sets/LEVELSET_UTILITIES.h>
#include <Core/Meshes/TRIANGULATED_SURFACE.h>
#include <Core/Nonlinear_Equations/ITERATIVE_SOLVER.h>
#include <Core/Polynomials/QUADRATIC.h>
#include <Core/Random_Numbers/RANDOM_NUMBERS.h>
namespace PhysBAM{
//#####################################################################
// Constructor
//#####################################################################
template<class TV> IMPLICIT_OBJECT<TV>::
IMPLICIT_OBJECT()
{
    Use_Secondary_Interpolation(false);
}
//#####################################################################
// Destructor
//#####################################################################
template<class TV> IMPLICIT_OBJECT<TV>::
~IMPLICIT_OBJECT()
{}
//#####################################################################
// Function Iterative_Solver_Tolerance
//#####################################################################
namespace{
//template<class T> inline T Iterative_Solver_Tolerance(){static_assert((T)false);}
//template<> inline float Iterative_Solver_Tolerance(){return (float).01;}
//template<> inline double Iterative_Solver_Tolerance(){return .001;}
}
//#####################################################################
// Function Box
//#####################################################################
template<class TV> RANGE<TV>& IMPLICIT_OBJECT<TV>::
Box()
{
    return box;
}
//#####################################################################
// Function Closest_Point_On_Boundary
//#####################################################################
template<class TV> TV IMPLICIT_OBJECT<TV>::
Closest_Point_On_Boundary(const TV& location,const T tolerance,const int max_iterations,T* distance) const
{
    T dist=(*this)(location);
    TV new_location=location-dist*Normal(location);
    if(!tolerance){
        if(distance) *distance=dist;
        return new_location;} // only take one iteration
    for(int i=1;i<max_iterations;i++){
        T dist=(*this)(new_location);
        if(abs(dist)<=tolerance) break;
        new_location-=dist*Normal(new_location);}
    if(distance) *distance=(new_location-location).Magnitude();
    return new_location;
}
//#####################################################################
// Function Signed_Distance
//#####################################################################
template<class TV> typename TV::SCALAR IMPLICIT_OBJECT<TV>::
Signed_Distance(const TV& location) const
{
    return (*this)(location);
}
//#####################################################################
// Function Inside
//#####################################################################
template<class TV> bool IMPLICIT_OBJECT<TV>::
Inside(const TV& location,const T thickness_over_two) const
{
    return box.Inside(location,thickness_over_two) && (*this)(location)<=-thickness_over_two;
}
//#####################################################################
// Function Outside
//#####################################################################
template<class TV> bool IMPLICIT_OBJECT<TV>::
Outside(const TV& location,const T thickness_over_two) const
{
    return box.Outside(location,thickness_over_two) || (*this)(location)>=thickness_over_two;
}
//#####################################################################
// Function Boundary
//#####################################################################
template<class TV> bool IMPLICIT_OBJECT<TV>::
Boundary(const TV& location,const T thickness_over_two) const
{
    return !Inside(location,thickness_over_two) && !Outside(location,thickness_over_two);
}
//#####################################################################
// Function Test_Diff
//#####################################################################
template<class TV> void IMPLICIT_OBJECT<TV>::
Test_Diff(const RANGE<TV>& range,bool test_hess) const
{
    T eps=(T)1e-6;
    RANDOM_NUMBERS<T> rand;
    TV X,dX;
    rand.Fill_Uniform(X,range);
    rand.Fill_Uniform(dX,-eps,eps);

    T p0=Signed_Distance(X),p1=Signed_Distance(X+dX);
    TV n0=Normal(X),n1=Normal(X+dX);
    Diff_Test("IMPLICIT_OBJECT diff",p0,X,n0,p1,X+dX,n1);
    if(test_hess)
        Diff_Test("IMPLICIT_OBJECT hess",n0,X,Hessian(X),n1,X+dX,Hessian(X+dX));
}
//#####################################################################
// Function Principal_Curvatures_Helper
//#####################################################################
template<class T> VECTOR<T,0> inline
Principal_Curvatures_Helper(const VECTOR<T,1>& n,const SYMMETRIC_MATRIX<T,1>& H)
{
    return VECTOR<T,0>();
}
//#####################################################################
// Function Principal_Curvatures_Helper
//#####################################################################
template<class T> VECTOR<T,1> inline
Principal_Curvatures_Helper(const VECTOR<T,2>& n,const SYMMETRIC_MATRIX<T,2>& H)
{
    VECTOR<T,2> t=n.Perpendicular();
    return VECTOR<T,1>(t.Dot(H*t));
}
//#####################################################################
// Function Principal_Curvatures_Helper
//#####################################################################
template<class T> VECTOR<T,2> inline
Principal_Curvatures_Helper(const VECTOR<T,3>& n,const SYMMETRIC_MATRIX<T,3>& H)
{
    SYMMETRIC_MATRIX<T,3> P=(T)1-Outer_Product(n),M=Conjugate(P,H);
    T trace=M.Trace();
    QUADRATIC<T> quadratic(-1,trace,sqr(M(0,2))+sqr(M(0,1))+sqr(M(1,2))-M(1,1)*M(0,0)-M(2,2)*M(0,0)-M(2,2)*M(1,1));
    quadratic.Compute_Roots();
    if(quadratic.roots == 0) return (T).5*VECTOR<T,2>(trace,trace);
    if(quadratic.roots == 1) return VECTOR<T,2>(quadratic.root[0],quadratic.root[0]);
    return VECTOR<T,2>(quadratic.root[0],quadratic.root[1]);
}
//#####################################################################
// Function Compute_Principal_Curvatures
//#####################################################################
template<class T,int d> VECTOR<T,d-1>
Compute_Principal_Curvatures(const VECTOR<T,d>& n,const SYMMETRIC_MATRIX<T,d>& H)
{
    return Principal_Curvatures_Helper(n,H);
}
//#####################################################################
// Function Generate_Triangles_Helper
//#####################################################################
namespace {
template<class T> TRIANGULATED_SURFACE<T>* Generate_Triangles_Helper(const LEVELSET_IMPLICIT_OBJECT<VECTOR<T,3> >* implicit)
{ // TODO(jontg): There should be a better implementation for this...
    TRIANGULATED_SURFACE<T>* surface=TRIANGULATED_SURFACE<T>::Create();
    MARCHING_CUBES<VECTOR<T,3> >::Create_Surface(*surface,implicit->levelset.Native_Grid(),at_nodes,implicit->levelset.phi);
    return surface;
}

template<class T> TRIANGULATED_SURFACE<T>* Generate_Triangles_Helper(const MULTIBODY_LEVELSET_IMPLICIT_OBJECT<VECTOR<T,3> >* implicit)
{
    ARRAY<TRIANGULATED_SURFACE<T>*> triangulated_surfaces(implicit->levelsets->m);
    for(int i=0;i<implicit->levelsets->m;i++) triangulated_surfaces(i)=Generate_Triangles(*(*implicit->levelsets)(i));
    ARRAY<FRAME<VECTOR<T,3> > > identity_frames(implicit->levelsets->m);
    TRIANGULATED_SURFACE<T>* union_object=TRIANGULATED_SURFACE<T>::Union_Mesh_Objects_Relatively(triangulated_surfaces,identity_frames);
    triangulated_surfaces.Delete_Pointers_And_Clean_Memory();
    return union_object;
}

template<class T> TRIANGULATED_SURFACE<T>* Generate_Triangles_Helper(const IMPLICIT_OBJECT_COMBINED<VECTOR<T,3> >* implicit)
{
    TRIANGULATED_SURFACE<T>* triangles_1=Generate_Triangles(*implicit->implicit_object1);
    TRIANGULATED_SURFACE<T>* triangles_2=Generate_Triangles(*implicit->implicit_object2);
    assert(triangles_1->particles.Size()==triangles_2->particles.Size());
    for(int i=0;i<triangles_1->particles.Size();i++)
        triangles_1->particles.X(i)=(1-implicit->alpha)*triangles_1->particles.X(i)+implicit->alpha*triangles_2->particles.X(i);
    return triangles_1;
}

template<class T> TRIANGULATED_SURFACE<T>* Generate_Triangles_Helper(const IMPLICIT_OBJECT_TRANSFORMED<VECTOR<T,3>,FRAME<VECTOR<T,3> > >* implicit)
{ // TODO(jontg): Better way to templatize this?
    TRIANGULATED_SURFACE<T>* surface=Generate_Triangles(*implicit->object_space_implicit_object);
    for(int i=0;i<surface->particles.Size();i++) surface->particles.X(i)=implicit->World_Space_Point(surface->particles.X(i));
    return surface;
}
}
//#####################################################################
// Function Generate_Triangles
//#####################################################################
template<class T> POINT_SIMPLICES_1D<T>* Generate_Triangles(IMPLICIT_OBJECT<VECTOR<T,1> > const& implicit) {PHYSBAM_FATAL_ERROR();}
template<class T> SEGMENTED_CURVE<VECTOR<T,2> >* Generate_Triangles(IMPLICIT_OBJECT<VECTOR<T,2> > const& implicit) {PHYSBAM_FATAL_ERROR();}
template<class T> TRIANGULATED_SURFACE<T>* Generate_Triangles(IMPLICIT_OBJECT<VECTOR<T,3> > const& implicit_input)
{
    typedef VECTOR<T,3> TV;
    if(const LEVELSET_IMPLICIT_OBJECT<TV>* implicit=dynamic_cast<const LEVELSET_IMPLICIT_OBJECT<TV>*>(&implicit_input)) return Generate_Triangles_Helper(implicit);
    else if(const MULTIBODY_LEVELSET_IMPLICIT_OBJECT<TV>* implicit=dynamic_cast<const MULTIBODY_LEVELSET_IMPLICIT_OBJECT<TV>*>(&implicit_input)) return Generate_Triangles_Helper(implicit);
    else if(const IMPLICIT_OBJECT_COMBINED<TV>* implicit=dynamic_cast<const IMPLICIT_OBJECT_COMBINED<TV>*>(&implicit_input)) return Generate_Triangles_Helper(implicit); 
    else if(const IMPLICIT_OBJECT_TRANSFORMED<TV,FRAME<TV> >* implicit=dynamic_cast<const IMPLICIT_OBJECT_TRANSFORMED<TV,FRAME<TV> >*>(&implicit_input)) return Generate_Triangles_Helper(implicit); 
    else if(const ANALYTIC_IMPLICIT_OBJECT<PLANE<T> >* implicit=dynamic_cast<const ANALYTIC_IMPLICIT_OBJECT<PLANE<T> >*>(&implicit_input)) return Generate_Triangles(implicit->analytic);
    else if(const ANALYTIC_IMPLICIT_OBJECT<RING<T> >* implicit=dynamic_cast<const ANALYTIC_IMPLICIT_OBJECT<RING<T> >*>(&implicit_input)) return Generate_Triangles(implicit->analytic);
    else if(const ANALYTIC_IMPLICIT_OBJECT<BOWL<T> >* implicit=dynamic_cast<const ANALYTIC_IMPLICIT_OBJECT<BOWL<T> >*>(&implicit_input)) return Generate_Triangles(implicit->analytic);
    else if(const ANALYTIC_IMPLICIT_OBJECT<SPHERE<TV> >* implicit=dynamic_cast<const ANALYTIC_IMPLICIT_OBJECT<SPHERE<TV> >*>(&implicit_input)) return Generate_Triangles(implicit->analytic);
    else if(const ANALYTIC_IMPLICIT_OBJECT<SMOOTH_GEAR<TV> >* implicit=dynamic_cast<const ANALYTIC_IMPLICIT_OBJECT<SMOOTH_GEAR<TV> >*>(&implicit_input)) return Generate_Triangles(implicit->analytic);
    else if(const ANALYTIC_IMPLICIT_OBJECT<BOUNDED_HORIZONTAL_PLANE<TV> >* implicit=dynamic_cast<const ANALYTIC_IMPLICIT_OBJECT<BOUNDED_HORIZONTAL_PLANE<TV> >*>(&implicit_input)) return Generate_Triangles(implicit->analytic);
    else if(const ANALYTIC_IMPLICIT_OBJECT<ORIENTED_BOX<TV> >* implicit=dynamic_cast<const ANALYTIC_IMPLICIT_OBJECT<ORIENTED_BOX<TV> >*>(&implicit_input)) return Generate_Triangles(implicit->analytic);
    else if(const ANALYTIC_IMPLICIT_OBJECT<RANGE<TV> >* implicit=dynamic_cast<const ANALYTIC_IMPLICIT_OBJECT<RANGE<TV> >*>(&implicit_input)) return Generate_Triangles(implicit->analytic);
    else if(const ANALYTIC_IMPLICIT_OBJECT<TORUS<T> >* implicit=dynamic_cast<const ANALYTIC_IMPLICIT_OBJECT<TORUS<T> >*>(&implicit_input)) return Generate_Triangles(implicit->analytic);
    else if(const ANALYTIC_IMPLICIT_OBJECT<CYLINDER<T> >* implicit=dynamic_cast<const ANALYTIC_IMPLICIT_OBJECT<CYLINDER<T> >*>(&implicit_input)) return Generate_Triangles(implicit->analytic);
    else{
        LOG::cout<<"Trying to generate triangles on an object of type:\n\t"<<typeid(implicit_input).name()<<std::endl;
        PHYSBAM_FATAL_ERROR();}
}
//#####################################################################
// Undefined Functions
//#####################################################################
template class IMPLICIT_OBJECT<VECTOR<double,1> >;
template class IMPLICIT_OBJECT<VECTOR<double,2> >;
template class IMPLICIT_OBJECT<VECTOR<double,3> >;
template class IMPLICIT_OBJECT<VECTOR<float,1> >;
template class IMPLICIT_OBJECT<VECTOR<float,2> >;
template class IMPLICIT_OBJECT<VECTOR<float,3> >;
template<class TV> void IMPLICIT_OBJECT<TV>::Update_Box(){PHYSBAM_FUNCTION_IS_NOT_DEFINED();}
template<class TV> void IMPLICIT_OBJECT<TV>::Update_Minimum_Cell_Size(const int maximum_depth){PHYSBAM_FUNCTION_IS_NOT_DEFINED();}
template<class TV> typename TV::SCALAR IMPLICIT_OBJECT<TV>::Minimum_Cell_Size_Within_Box(const RANGE<TV>& box) const {PHYSBAM_FUNCTION_IS_NOT_DEFINED();}
template<class TV> typename TV::SCALAR IMPLICIT_OBJECT<TV>::operator()(const TV& location) const {PHYSBAM_FUNCTION_IS_NOT_DEFINED();}
template<class TV> typename TV::SCALAR IMPLICIT_OBJECT<TV>::Extended_Phi(const TV& location) const {PHYSBAM_FUNCTION_IS_NOT_DEFINED();}
template<class TV> typename TV::SCALAR IMPLICIT_OBJECT<TV>::Phi_Secondary(const TV& location) const {PHYSBAM_FUNCTION_IS_NOT_DEFINED();}
template<class TV> TV IMPLICIT_OBJECT<TV>::Normal(const TV& location,const int aggregate) const {PHYSBAM_FUNCTION_IS_NOT_DEFINED();}
template<class TV> TV IMPLICIT_OBJECT<TV>::Extended_Normal(const TV& location,const int aggregate) const {PHYSBAM_FUNCTION_IS_NOT_DEFINED();}
template<class TV> void IMPLICIT_OBJECT<TV>::Compute_Normals(){PHYSBAM_FUNCTION_IS_NOT_DEFINED();}
template<class TV> void IMPLICIT_OBJECT<TV>::Compute_Cell_Minimum_And_Maximum(const bool recompute_if_exists){PHYSBAM_FUNCTION_IS_NOT_DEFINED();}
template<class TV> void IMPLICIT_OBJECT<TV>::Rescale(const T scaling_factor) {PHYSBAM_FUNCTION_IS_NOT_DEFINED();}
template<class TV> void IMPLICIT_OBJECT<TV>::Translate(const TV& translation){PHYSBAM_FUNCTION_IS_NOT_DEFINED();}
template<class TV> void IMPLICIT_OBJECT<TV>::Inflate(const T inflation_distance){PHYSBAM_FUNCTION_IS_NOT_DEFINED();}
template<class TV> bool IMPLICIT_OBJECT<TV>::Lazy_Inside(const TV& location,const T contour_value) const {PHYSBAM_FUNCTION_IS_NOT_DEFINED();}
template<class TV> bool IMPLICIT_OBJECT<TV>::Lazy_Inside_And_Value(const TV& location,T& phi_value,const T contour_value) const {PHYSBAM_FUNCTION_IS_NOT_DEFINED();}
template<class TV> bool IMPLICIT_OBJECT<TV>::Lazy_Inside_Extended_Levelset(const TV& unclamped_X,const T contour_value) const {PHYSBAM_FUNCTION_IS_NOT_DEFINED();}
template<class TV> bool IMPLICIT_OBJECT<TV>::Lazy_Inside_Extended_Levelset_And_Value(const TV& unclamped_X,T& phi_value,const T contour_value) const {PHYSBAM_FUNCTION_IS_NOT_DEFINED();}
template<class TV> bool IMPLICIT_OBJECT<TV>::Lazy_Outside(const TV& location,const T contour_value) const {PHYSBAM_FUNCTION_IS_NOT_DEFINED();}
template<class TV> bool IMPLICIT_OBJECT<TV>::Lazy_Outside_Extended_Levelset(const TV& unclamped_X,const T contour_value) const {PHYSBAM_FUNCTION_IS_NOT_DEFINED();}
template<class TV> bool IMPLICIT_OBJECT<TV>::Lazy_Outside_Extended_Levelset_And_Value(const TV& unclamped_X,T& phi_value,const T contour_value) const {PHYSBAM_FUNCTION_IS_NOT_DEFINED();}
template<class TV> SYMMETRIC_MATRIX<typename TV::SCALAR,TV::m> IMPLICIT_OBJECT<TV>::Hessian(const TV& X) const {PHYSBAM_FUNCTION_IS_NOT_DEFINED();}
template<class TV> typename TV::SCALAR IMPLICIT_OBJECT<TV>::Integration_Step(const T phi) const {PHYSBAM_FUNCTION_IS_NOT_DEFINED();}
template<class TV> typename TV::SCALAR IMPLICIT_OBJECT<TV>::Minimum_Cell_Size() const {PHYSBAM_FUNCTION_IS_NOT_DEFINED();} 
//#####################################################################
template VECTOR<double,0> Compute_Principal_Curvatures<double,1>(VECTOR<double,1> const&,SYMMETRIC_MATRIX<double,1> const&);
template VECTOR<double,1> Compute_Principal_Curvatures<double,2>(VECTOR<double,2> const&,SYMMETRIC_MATRIX<double,2> const&);
template VECTOR<double,2> Compute_Principal_Curvatures<double,3>(VECTOR<double,3> const&,SYMMETRIC_MATRIX<double,3> const&);
template VECTOR<float,0> Compute_Principal_Curvatures<float,1>(VECTOR<float,1> const&,SYMMETRIC_MATRIX<float,1> const&);
template VECTOR<float,1> Compute_Principal_Curvatures<float,2>(VECTOR<float,2> const&,SYMMETRIC_MATRIX<float,2> const&);
template VECTOR<float,2> Compute_Principal_Curvatures<float,3>(VECTOR<float,3> const&,SYMMETRIC_MATRIX<float,3> const&);
template POINT_SIMPLICES_1D<float>* Generate_Triangles(const IMPLICIT_OBJECT<VECTOR<float,1> >&);
template SEGMENTED_CURVE<VECTOR<float,2> >* Generate_Triangles(const IMPLICIT_OBJECT<VECTOR<float,2> >&);
template TRIANGULATED_SURFACE<float>* Generate_Triangles(const IMPLICIT_OBJECT<VECTOR<float,3> >&); 
template POINT_SIMPLICES_1D<double>* Generate_Triangles(const IMPLICIT_OBJECT<VECTOR<double,1> >&);
template SEGMENTED_CURVE<VECTOR<double,2> >* Generate_Triangles(const IMPLICIT_OBJECT<VECTOR<double,2> >&);
template TRIANGULATED_SURFACE<double>* Generate_Triangles(const IMPLICIT_OBJECT<VECTOR<double,3> >&); 
}
