//#####################################################################
// Copyright 2003-2008, Christopher Allocco, Ronald Fedkiw, Geoffrey Irving, Sergey Koltakov, Neil Molino, Duc Nguyen, Andrew Selle, Eftychios Sifakis, Joseph Teran.
// This file is part of PhysBAM whose distribution is governed by the license contained in the accompanying file PHYSBAM_COPYRIGHT.txt.
//#####################################################################
#include <Core/Meshes/TRIANGULATED_AREA.h>
#include <Core/Arrays/INDIRECT_ARRAY.h>
#include <Core/Basic_Geometry/TRIANGLE_2D.h>
#include <Core/Data_Structures/HASHTABLE.h>
#include <Core/Grids/GRID.h>
#include <Core/Implicit_Objects/IMPLICIT_OBJECT.h>
#include <Core/Log/LOG.h>
#include <Core/Meshes/SEGMENTED_CURVE.h>
#include <Core/Registry/STRUCTURE_REGISTRY.h>
#include <Core/Spatial_Acceleration/MESH_HIERARCHY.h>
using namespace PhysBAM;
//#####################################################################
// Constructor
//#####################################################################
template<class T> MESH_OBJECT<VECTOR<T,2>,2>::
MESH_OBJECT()
    :MESH_OBJECT(*new TRIANGLE_MESH,*new GEOMETRY_PARTICLES<TV>)
{
    this->need_destroy_mesh=this->need_destroy_particles=true;
}
//#####################################################################
// Constructor
//#####################################################################
template<class T> MESH_OBJECT<VECTOR<T,2>,2>::
MESH_OBJECT(TRIANGLE_MESH& triangle_mesh_input,GEOMETRY_PARTICLES<TV>& particles_input)
    :MESH_OBJECT_BASE<TV,TRIANGLE_MESH>(triangle_mesh_input,particles_input),
    segmented_curve(0),nodal_areas(0)
{
}
//#####################################################################
// Destructor
//#####################################################################
template<class T> MESH_OBJECT<VECTOR<T,2>,2>::
~MESH_OBJECT()
{
    delete segmented_curve;
    delete nodal_areas;
}
//#####################################################################
// Function Clean_Memory
//#####################################################################
template<class T> void MESH_OBJECT<VECTOR<T,2>,2>::
Clean_Memory()
{
    MESH_OBJECT_BASE<TV,TRIANGLE_MESH>::Clean_Memory();
    delete segmented_curve;
    segmented_curve=0;
    delete nodal_areas;
    nodal_areas=0;
}
//#####################################################################
// Function Refresh_Auxiliary_Structures_Helper
//#####################################################################
template<class T> void MESH_OBJECT<VECTOR<T,2>,2>::
Refresh_Auxiliary_Structures_Helper()
{
    if(hierarchy) this->Initialize_Hierarchy();
}
//#####################################################################
// Funcion Initialize_Square_Mesh_And_Particles
//#####################################################################
template<class T> void MESH_OBJECT<VECTOR<T,2>,2>::
Initialize_Square_Mesh_And_Particles(const GRID<TV>& grid,const bool reverse_triangles)
{
    int m=grid.Numbers_Of_Nodes().x,n=grid.Numbers_Of_Nodes().y,particle=0;
    particles.Delete_All_Elements();mesh.Initialize_Square_Mesh(m,n,reverse_triangles);particles.Add_Elements(m*n);
    for(int j=0;j<n;j++) for(int i=0;i<m;i++) particles.X(particle++)=grid.Node_X(IV(i,j));
}
//#####################################################################
// Funcion Initialize_Circle_Mesh_And_Particles
//#####################################################################
template<class T> void MESH_OBJECT<VECTOR<T,2>,2>::
Initialize_Circle_Mesh_And_Particles(const T outer_radius,const T inner_radius,const int num_radial,const int num_tangential)
{
    int particle=0;particles.Delete_All_Elements();
    mesh.Initialize_Circle_Mesh(num_radial,num_tangential);particles.Add_Elements(num_radial*num_tangential);
    for(int j=0;j<num_radial;j++) for(int i=0;i<num_tangential;i++){
        T r=T(j)/T(num_tangential-1)*(outer_radius-inner_radius)+inner_radius,theta=T(i)/T(num_tangential)*(T)2*pi_v<T>; 
        particles.X(particle++)=TV(r*cos(theta),r*sin(theta));}
}
//#####################################################################
// Funcion Initialize_Herringbone_Mesh_And_Particles
//#####################################################################
template<class T> void MESH_OBJECT<VECTOR<T,2>,2>::
Initialize_Herring_Bone_Mesh_And_Particles(const GRID<TV>& grid)
{
    int m=grid.Numbers_Of_Nodes().x,n=grid.Numbers_Of_Nodes().y,particle=0;
    particles.Delete_All_Elements();mesh.Initialize_Herring_Bone_Mesh(m,n);particles.Add_Elements(m*n);
    for(int j=0;j<n;j++) for(int i=0;i<m;i++) particles.X(particle++)=grid.Node_X(IV(i,j));
}
//#####################################################################
// Funcion Initialize_Equilateral_Mesh_And_Particles
//#####################################################################
template<class T> void MESH_OBJECT<VECTOR<T,2>,2>::
Initialize_Equilateral_Mesh_And_Particles(const GRID<TV>& grid)
{
    int m=grid.Numbers_Of_Nodes().x,n=grid.Numbers_Of_Nodes().y,particle=0;
    particles.Delete_All_Elements();mesh.Initialize_Equilateral_Mesh(m,n);particles.Add_Elements(m*n);
    for(int j=0;j<n;j++)
        for(int i=0;i<m;i++){
            TV X=grid.Node_X(IV(i,j));
            if(j%2) X.x+=(T).5*grid.dX.x;
            particles.X(particle++)=X;}
}
//#####################################################################
// Function Inside
// note: return the first triangle that it is inside of (including boundary), otherwise returns -1
//#####################################################################
template<class T> int MESH_OBJECT<VECTOR<T,2>,2>::
Inside(const TV& location,const T thickness_over_two) const
{
    PHYSBAM_ASSERT(hierarchy);
    int found=-1;
    hierarchy->Traverse(
        [=](const RANGE<TV>& r){return !r.Outside(location,thickness_over_two);},
        [=,this,&found](int t)
        {
            int i,j,k;
            mesh.elements(t).Get(i,j,k);
            if(!TRIANGLE_2D<T>::Outside(location,particles.X(i),particles.X(j),
                    particles.X(k),thickness_over_two))
            {
                found=t;
                return false;
            }
            return true;
        });
    return found;
}
//#####################################################################
// Function Inside_Any_Simplex
//#####################################################################
template<class T> bool MESH_OBJECT<VECTOR<T,2>,2>::
Inside_Any_Simplex(const TV& location,int& triangle_id,const T thickness_over_two) const
{
    triangle_id=Inside(location,thickness_over_two);
    return triangle_id>=0;
}
//#####################################################################
// Function Check_Signed_Area_And_Make_Consistent
//#####################################################################
template<class T> void MESH_OBJECT<VECTOR<T,2>,2>::
Check_Signed_Area_And_Make_Consistent(const int triangle,const bool verbose)
{
    VECTOR<int,3>& nodes=mesh.elements(triangle);
    if(TRIANGLE_2D<T>::Signed_Area(particles.X(nodes[0]),particles.X(nodes[1]),particles.X(nodes[2])) < 0){
        if(verbose) LOG::cout<<"triangle number "<<triangle<<" is oriented improperly."<<std::endl;
        exchange(nodes[1],nodes[2]);}
}
//#####################################################################
// Function Check_Signed_Areas_And_Make_Consistent
//#####################################################################
template<class T> void MESH_OBJECT<VECTOR<T,2>,2>::
Check_Signed_Areas_And_Make_Consistent(const bool verbose)
{
    for(int t=0;t<mesh.elements.m;t++) Check_Signed_Area_And_Make_Consistent(t,verbose);
}
//#####################################################################
// Function Initialize_Segmented_Curve
//#####################################################################
template<class T> void MESH_OBJECT<VECTOR<T,2>,2>::
Initialize_Segmented_Curve()
{
    delete segmented_curve;if(!mesh.boundary_mesh) mesh.Initialize_Boundary_Mesh();
    segmented_curve=new SEGMENTED_CURVE<TV>(*mesh.boundary_mesh,particles);
}
//#####################################################################
// Function Centroid
//#####################################################################
template<class T> VECTOR<T,2> MESH_OBJECT<VECTOR<T,2>,2>::
Centroid(const int triangle) const
{
    int i,j,k;mesh.elements(triangle).Get(i,j,k);
    return ((T)1/3)*(particles.X(i)+particles.X(j)+particles.X(k));
}
//#####################################################################
// Function Rescale
//#####################################################################
template<class T> void MESH_OBJECT<VECTOR<T,2>,2>::
Rescale(const T scaling_factor)
{
    Rescale(scaling_factor,scaling_factor);
}
//#####################################################################
// Function Rescale
//#####################################################################
template<class T> void MESH_OBJECT<VECTOR<T,2>,2>::
Rescale(const T scaling_x,const T scaling_y)
{
    for(int k=0;k<particles.Size();k++) particles.X(k)*=TV(scaling_x,scaling_y);
    Refresh_Auxiliary_Structures();
}
//#####################################################################
// Function Area
//#####################################################################
template<class T> T MESH_OBJECT<VECTOR<T,2>,2>::
Area(const int triangle) const
{
    int i,j,k;mesh.elements(triangle).Get(i,j,k);
    return TRIANGLE_2D<T>::Area(particles.X(i),particles.X(j),particles.X(k));
}
//#####################################################################
// Function Signed_Area
//#####################################################################
template<class T> T MESH_OBJECT<VECTOR<T,2>,2>::
Signed_Area(const int triangle) const
{
    int i,j,k;mesh.elements(triangle).Get(i,j,k);
    return TRIANGLE_2D<T>::Signed_Area(particles.X(i),particles.X(j),particles.X(k));
}
//#####################################################################
// Funcion Minimum_Area
//#####################################################################
template<class T> T MESH_OBJECT<VECTOR<T,2>,2>::
Minimum_Area(int* index) const
{
    int k=0;T minimum=FLT_MAX;
    for(int t=0;t<mesh.elements.m;t++){
        int node1=mesh.elements(t)(0),node2=mesh.elements(t)(1),node3=mesh.elements(t)(2);
        T temp=TRIANGLE_2D<T>::Area(particles.X(node1),particles.X(node2),particles.X(node3));
        if(temp < minimum){minimum=temp;k=t;}}
    if(index) *index=k;
    return minimum;
}
//#####################################################################
// Funcion Minimum_Signed_Area
//#####################################################################
template<class T> T MESH_OBJECT<VECTOR<T,2>,2>::
Minimum_Signed_Area(int* index) const
{
    int k=0;T minimum=FLT_MAX;
    for(int t=0;t<mesh.elements.m;t++){
        int node1=mesh.elements(t)(0),node2=mesh.elements(t)(1),node3=mesh.elements(t)(2);
        T temp=TRIANGLE_2D<T>::Signed_Area(particles.X(node1),particles.X(node2),particles.X(node3));
        if(temp < minimum){minimum=temp;k=t;}}
    if(index) *index=k;
    return minimum;
}
//#####################################################################
// Funcion Total_Area
//#####################################################################
template<class T> T MESH_OBJECT<VECTOR<T,2>,2>::
Total_Area() const
{
    T area=0;
    for(int t=0;t<mesh.elements.m;t++){
        int node1=mesh.elements(t)(0),node2=mesh.elements(t)(1),node3=mesh.elements(t)(2);
        area+=TRIANGLE_2D<T>::Area(particles.X(node1),particles.X(node2),particles.X(node3));}
    return area;
}
//#####################################################################
// Funcion Minimum_Altitude
//#####################################################################
template<class T> T MESH_OBJECT<VECTOR<T,2>,2>::
Minimum_Altitude(int* index) const
{
    int k=0;T minimum=FLT_MAX;
    for(int t=0;t<mesh.elements.m;t++){
        int node1=mesh.elements(t)(0),node2=mesh.elements(t)(1),node3=mesh.elements(t)(2);
        T temp=TRIANGLE_2D<T>::Minimum_Altitude(particles.X(node1),particles.X(node2),particles.X(node3));
        if(temp < minimum){minimum=temp;k=t;}}
    if(index) *index=k;
    return minimum;
}
//#####################################################################
// Funcion Minimum_Edge_Length
//#####################################################################
template<class T> T MESH_OBJECT<VECTOR<T,2>,2>::
Minimum_Edge_Length(int* index) const
{
    int t_save=0;T minimum_squared=FLT_MAX;
    for(int t=0;t<mesh.elements.m;t++){int i,j,k;mesh.elements(t).Get(i,j,k);
        T temp=TRIANGLE_2D<T>::Minimum_Edge_Length_Squared(particles.X(i),particles.X(j),particles.X(k));if(temp < minimum_squared){minimum_squared=temp;t_save=t;}}
    if(index) *index=t_save;
    return sqrt(minimum_squared);
}
//#####################################################################
// Funcion Maximum_Edge_Length
//#####################################################################
template<class T> T MESH_OBJECT<VECTOR<T,2>,2>::
Maximum_Edge_Length(int* index) const
{
    int t_save=0;T maximum_squared=FLT_MAX;
    for(int t=0;t<mesh.elements.m;t++){int i,j,k;mesh.elements(t).Get(i,j,k);
        T temp=TRIANGLE_2D<T>::Maximum_Edge_Length_Squared(particles.X(i),particles.X(j),particles.X(k));if(temp < maximum_squared){maximum_squared=temp;t_save=t;}}
    if(index) *index=t_save;
    return sqrt(maximum_squared);
}
//#####################################################################
// Function Volume_Incident_On_A_Particle
//#####################################################################
template<class T> T MESH_OBJECT<VECTOR<T,2>,2>::
Area_Incident_On_A_Particle(const int particle_index)
{
    int incident_elements_defined=mesh.incident_elements!=0;if(!incident_elements_defined) mesh.Initialize_Incident_Elements();
    T total_incident_area=0;
    for(int t=0;t<(*mesh.incident_elements)(particle_index).m;t++){
        int i,j,k;mesh.elements((*mesh.incident_elements)(particle_index)(t)).Get(i,j,k);
        total_incident_area+=TRIANGLE_2D<T>::Area(particles.X(i),particles.X(j),particles.X(k));}
    if(!incident_elements_defined){delete mesh.incident_elements;mesh.incident_elements=0;}
    return total_incident_area;
}
//#####################################################################
namespace PhysBAM{
template class MESH_OBJECT<VECTOR<float,2>,2>;
template class MESH_OBJECT<VECTOR<double,2>,2>;
}
