//#####################################################################
// Copyright 2006-2008, Geoffrey Irving, Andrew Selle, Tamar Shinar, Eftychios Sifakis.
// This file is part of PhysBAM whose distribution is governed by the license contained in the accompanying file PHYSBAM_COPYRIGHT.txt.
//#####################################################################
// Class SIMPLEX_MESH
//#####################################################################
#include <Core/Meshes/SIMPLEX_MESH_BASE.h>
#include <Core/Arrays/INDIRECT_ARRAY.h>
#include <Core/Data_Structures/HASHTABLE.h>
#include <Core/Data_Structures/UNION_FIND.h>
#include <Core/Log/LOG.h>
#include <Core/Meshes/SEGMENT_MESH.h>
#include <Core/Read_Write/READ_WRITE.h>
#include <climits>
namespace PhysBAM{
//#####################################################################
// Constructor
//#####################################################################
template<int d> SIMPLEX_MESH_BASE<d>::
SIMPLEX_MESH_BASE()
    :number_nodes(0)
{}
//#####################################################################
// Constructor
//#####################################################################
template<int d> SIMPLEX_MESH_BASE<d>::
SIMPLEX_MESH_BASE(const int number_nodes_input,const ARRAY<VECTOR<int,d+1> >& simplex_list)
{
    Initialize_Mesh(number_nodes_input,simplex_list);
}
//#####################################################################
// Constructor
//#####################################################################
template<int d> SIMPLEX_MESH_BASE<d>::
SIMPLEX_MESH_BASE(const SIMPLEX_MESH_BASE& mesh)
{
    Initialize_Mesh(mesh);
}
//#####################################################################
// Destructor
//#####################################################################
template<int d> SIMPLEX_MESH_BASE<d>::
~SIMPLEX_MESH_BASE()
{
    delete neighbor_nodes;
    delete incident_elements;
    delete adjacent_elements;
    delete neighbor_elements;
}
//#####################################################################
// Function Clean_Memory
//#####################################################################
template<int d> void SIMPLEX_MESH_BASE<d>::
Clean_Memory()
{
    elements.Clean_Memory();
    Delete_Auxiliary_Structures();
}
//#####################################################################
// Function Delete_Auxiliary_Structures
//#####################################################################
template<int d> void SIMPLEX_MESH_BASE<d>::
Delete_Auxiliary_Structures()
{
    delete neighbor_nodes;
    neighbor_nodes=0;
    delete incident_elements;
    incident_elements=0;
    delete adjacent_elements;
    adjacent_elements=0;
    delete neighbor_elements;
    neighbor_elements=0;
}
//#####################################################################
// Function Refresh_Auxiliary_Structures
//#####################################################################
template<int d> void SIMPLEX_MESH_BASE<d>::
Refresh_Auxiliary_Structures()
{
    if(neighbor_nodes) Initialize_Neighbor_Nodes();
    if(incident_elements) Initialize_Incident_Elements();
    if(adjacent_elements) Initialize_Adjacent_Elements();
    if(neighbor_elements) Initialize_Neighbor_Elements();
}
//#####################################################################
// Function Simplex
//#####################################################################
// requires incident_elements
// returns the simplex containing these nodes - returns 0 if no simplex contains these nodes
template<int d> int SIMPLEX_MESH_BASE<d>::
Simplex(const VECTOR<int,d+1>& nodes) const
{
    assert(incident_elements);
    // find shortest list of elements
    int short_list=nodes[0];VECTOR<int,d> check=nodes.Remove_Index(0);
    for(int i=0;i<d;i++)if((*incident_elements)(check[i]).m < (*incident_elements)(short_list).m) exchange(short_list,check[i]);
    // search short list for other nodes
    for(int k=0;k<(*incident_elements)(short_list).m;k++){int t=(*incident_elements)(short_list)(k);
        if(Nodes_In_Simplex(check,t)) return t;}
    return -1;
}
//#####################################################################
// Function Initialize_Neighbor_Nodes
//#####################################################################
template<int d> void SIMPLEX_MESH_BASE<d>::
Initialize_Neighbor_Nodes()
{
    delete neighbor_nodes;neighbor_nodes=new ARRAY<ARRAY<int> >(number_nodes);
    for(int t=0;t<elements.m;t++){VECTOR<int,d+1>& element=elements(t);
        for(int i=0;i<element.m;i++)for(int j=0;j<element.m;j++)if(i!=j)
            (*neighbor_nodes)(element[i]).Append_Unique(element[j]);}
    for(int p=0;p<number_nodes;p++) (*neighbor_nodes)(p).Compact(); // remove extra space
}
//#####################################################################
// Function Initialize_Incident_Elements
//#####################################################################
template<int d> void SIMPLEX_MESH_BASE<d>::
Initialize_Incident_Elements()
{
    delete incident_elements;incident_elements=new ARRAY<ARRAY<int> >(number_nodes);
    for(int t=0;t<elements.m;t++){VECTOR<int,d+1>& element=elements(t); // for each element, put it on each of its nodes' lists of simplices
        for(int i=0;i<element.m;i++) (*incident_elements)(element[i]).Append(t);}
    for(int p=0;p<number_nodes;p++) (*incident_elements)(p).Compact(); // remove extra space
}
//#####################################################################
// Function Initialize_Adjacent_Elements
//#####################################################################
template<int d> void SIMPLEX_MESH_BASE<d>::
Initialize_Adjacent_Elements()
{
    delete adjacent_elements;adjacent_elements=new ARRAY<ARRAY<int> >(elements.m);
    bool incident_elements_defined=incident_elements!=0;if(!incident_elements) Initialize_Incident_Elements();
    for(int t=0;t<elements.m;t++){VECTOR<int,d+1>& element=elements(t); // for each simplex, make the list of adjacent simplices
        for(int i=0;i<element.m;i++) Find_And_Append_Adjacent_Elements(t,element.Remove_Index(i));}
    if(!incident_elements_defined){delete incident_elements;incident_elements=0;}
    for(int t=0;t<elements.m;t++) (*adjacent_elements)(t).Compact(); // remove extra space
}
//#####################################################################
// Function Find_And_Append_Adjacent_Elements
//#####################################################################
// find adjacent simplices that contains face and append them to the adjacency list
template<int d> void SIMPLEX_MESH_BASE<d>::
Find_And_Append_Adjacent_Elements(const int simplex,const VECTOR<int,d>& face)
{
    int first_node=face[0];VECTOR<int,d-1> other_nodes=face.Remove_Index(0);
    for(int t=0;t<(*incident_elements)(first_node).m;t++){
        int simplex2=(*incident_elements)(first_node)(t);
        if(simplex2!=simplex && Nodes_In_Simplex(other_nodes,simplex2))
            (*adjacent_elements)(simplex).Append_Unique(simplex2);}
}
template<> void SIMPLEX_MESH_BASE<0>:: // specialize segment case since other_nodes would be 0 size
Find_And_Append_Adjacent_Elements(const int simplex,const VECTOR<int,0>& face)
{return;}
//#####################################################################
// Function Initialize_Neighbor_Elements
//#####################################################################
template<int d> void SIMPLEX_MESH_BASE<d>::
Initialize_Neighbor_Elements()
{
    static_assert(d<4);static const int expected_neighbors=d==1?4:d==2?20:d==3?50:0;
    delete neighbor_elements;neighbor_elements=new ARRAY<ARRAY<int> >(elements.m);
    bool incident_elements_defined=incident_elements!=0;if(!incident_elements) Initialize_Incident_Elements();
    for(int t=0;t<elements.m;t++){VECTOR<int,d+1>& element=elements(t); // for each simplex, make the list of neighbor simplices
        ARRAY<int>& neighbors=(*neighbor_elements)(t);neighbors.Preallocate(expected_neighbors);
        neighbors.Append((*incident_elements)(element[0]));
        for(int i=1;i<element.m;i++) neighbors.Append_Unique((*incident_elements)(element[i]));
        neighbors.Remove_Index_Lazy(neighbors.Find(t));neighbors.Compact();}
    if(!incident_elements_defined){delete incident_elements;incident_elements=0;}
}
//#####################################################################
// Function Add_Element_If_Not_Already_There
//#####################################################################
// will ONLY incrementally update incident_elements, assumes number_nodes is already set to accommodate the new element
template<int d> int SIMPLEX_MESH_BASE<d>::
Add_Element_If_Not_Already_There(const VECTOR<int,d+1>& nodes)
{
    int element=Simplex(nodes);
    if(element>=0) return element;
    int id=elements.Append(nodes);
    for(int i=0;i<d+1;i++) (*incident_elements)(nodes[i]).Append(id);
    return id;
}
//#####################################################################
// Function Delete_Elements_With_Missing_Nodes
//#####################################################################
template<int d> int SIMPLEX_MESH_BASE<d>::
Delete_Elements_With_Missing_Nodes()
{
    int m_save=elements.m;
    for(int t=elements.m-1;t>=0;t--)if(elements(t).Contains(-1)) elements.Remove_Index_Lazy(t);
    elements.Compact();Refresh_Auxiliary_Structures();
    return m_save-elements.m;
}
//#####################################################################
// Function Delete_Sorted_Elements
//#####################################################################
template<int d> void SIMPLEX_MESH_BASE<d>::
Delete_Sorted_Elements(const ARRAY<int>& deletion_list)
{
    elements.Remove_Sorted_Indices_Lazy(deletion_list);
    Refresh_Auxiliary_Structures();
}
//#####################################################################
// Function Delete_Sorted_Elements
// index_map. Maps from old element index to new element index. (0 if deleted, no entry if unchanged)
//#####################################################################
template<int d> void SIMPLEX_MESH_BASE<d>::
Delete_Sorted_Elements(const ARRAY<int>& deletion_list,HASHTABLE<int,int>& index_map) 
{
    int curr=0;
    HASHTABLE<int,int> reverse_index_map;
    ARRAY<int> unique_deletion_list;
    for(int k=deletion_list.m-1;k>=0;k--) if(deletion_list(k)!=curr){
        curr=deletion_list(k);int previous_index=elements.m;
        unique_deletion_list.Append(curr);
        if(curr!=elements.m){
            if(reverse_index_map.Get(elements.m,previous_index)) reverse_index_map.Delete(elements.m);
            reverse_index_map.Insert(curr,previous_index);}
        elements.Remove_Index_Lazy(curr);}
    for(auto p:reverse_index_map) index_map.Insert(p.second,p.first);
    for(int i=0;i<unique_deletion_list.m;i++) index_map.Insert(unique_deletion_list(i),0);
    elements.Compact();Refresh_Auxiliary_Structures();
}
//#####################################################################
// Function Delete_Elements
//#####################################################################
template<int d> void SIMPLEX_MESH_BASE<d>::
Delete_Elements(ARRAY<int> deletion_list)
{
    deletion_list.Sort();Delete_Sorted_Elements(deletion_list);
}
//#####################################################################
// Function Update_Adjacent_Elements_From_Incident_Elements
//#####################################################################
template<int d> void SIMPLEX_MESH_BASE<d>::
Update_Adjacent_Elements_From_Incident_Elements(const int node)
{
    if(!adjacent_elements){Initialize_Adjacent_Elements();return;}
    bool incident_elements_defined=incident_elements!=0;if(!incident_elements) Initialize_Incident_Elements();
    for(int a=0;a<(*incident_elements)(node).m;a++){
        int t=(*incident_elements)(node)(a);VECTOR<int,d+1>& element=elements(t);
        (*adjacent_elements)(t).Remove_All();
        for(int i=0;i<element.m;i++) Find_And_Append_Adjacent_Elements(t,element.Remove_Index(i));
        (*adjacent_elements)(t).Compact();}
    if(!incident_elements_defined){delete incident_elements;incident_elements=0;}
}
//#####################################################################
// Function Update_Neighbor_Nodes_From_Incident_Elements
//#####################################################################
template<int d> void SIMPLEX_MESH_BASE<d>::
Update_Neighbor_Nodes_From_Incident_Elements(const int node)
{
    if(!neighbor_nodes){Initialize_Neighbor_Nodes();return;}
    bool incident_elements_defined=incident_elements!=0;if(!incident_elements) Initialize_Incident_Elements();
    (*neighbor_nodes)(node).Remove_All();
    for(int t=0;t<(*incident_elements)(node).m;t++){
        VECTOR<int,d+1>& element=elements((*incident_elements)(node)(t));assert(element.Contains(node));
        for(int i=0;i<element.m;i++)if(element[i]!=node) (*neighbor_nodes)(node).Append_Unique(element[i]);}
    (*neighbor_nodes)(node).Compact();
    if(!incident_elements_defined){delete incident_elements;incident_elements=0;}
}
//#####################################################################
// Function Assert_Consistent
//#####################################################################
template<int d> bool SIMPLEX_MESH_BASE<d>::
Assert_Consistent() const
{
    if(neighbor_nodes) assert(neighbor_nodes->m==number_nodes);
    if(incident_elements) assert(incident_elements->m==number_nodes);
    if(adjacent_elements) assert(adjacent_elements->m==elements.m);
    if(neighbor_elements) assert(neighbor_elements->m==elements.m);
    return true;
}
//#####################################################################
// Function Set_Number_Nodes
//#####################################################################
template<int d> void SIMPLEX_MESH_BASE<d>::
Set_Number_Nodes(const int number_nodes_input)
{
    assert(Assert_Consistent());
    number_nodes=number_nodes_input;
    if(neighbor_nodes) neighbor_nodes->Resize(number_nodes);
    if(incident_elements) incident_elements->Resize(number_nodes);
}
//#####################################################################
// Function Add_Dependencies
//#####################################################################
template<int d> void SIMPLEX_MESH_BASE<d>::
Add_Dependencies(SEGMENT_MESH& dependency_mesh) const
{
    for(int t=0;t<elements.m;t++) for(int i=0;i<d;i++) for(int j=i;j<d+1;j++) dependency_mesh.Add_Element_If_Not_Already_There(VECTOR<int,2>(elements(t)[i],elements(t)[j]));
}
//#####################################################################
// Function Mark_Nodes_Referenced
//#####################################################################
template<int d> template<class T> void SIMPLEX_MESH_BASE<d>::
Mark_Nodes_Referenced(ARRAY<T>& marks,const T& mark) const
{
    for(int e=0;e<elements.m;e++) marks.Subset(elements(e)).Fill(mark);
}
//#####################################################################
// Function Simplices_On_Subsimplex
//#####################################################################
template<int d> template<int d2> void SIMPLEX_MESH_BASE<d>::
Simplices_On_Subsimplex(const VECTOR<int,d2>& subsimplex_nodes,ARRAY<int>& simplices_on_subsimplex) const
{
    assert(incident_elements);
    const ARRAY<int>& incident=(*incident_elements)(subsimplex_nodes[0]);
    VECTOR<int,d2-1> other_nodes=subsimplex_nodes.Remove_Index(0);
    for(int i=0;i<incident.m;i++){
        int simplex=incident(i);
        if(Nodes_In_Simplex(other_nodes,simplex)) simplices_on_subsimplex.Append(simplex);}
}
//#####################################################################
// Function Read
//#####################################################################
template<int d> template<class RW> void SIMPLEX_MESH_BASE<d>::
Read(UNTYPED_ISTREAM<RW>& input)
{
    Clean_Memory();
    int backward_compatible;
    Get(input,'(');
    Read_Binary(input,number_nodes,backward_compatible,elements);
    if(elements.m){
        int min_index=elements.Flattened().Min(),max_index=elements.Flattened().Min();
        if(number_nodes<0) throw READ_ERROR(LOG::sprintf("Invalid negative number_nodes = %d",number_nodes));
        if(min_index<0) throw READ_ERROR(LOG::sprintf("Invalid vertex index %d",min_index));
        if(max_index>=number_nodes) throw READ_ERROR(LOG::sprintf("Invalid vertex index %d (number_nodes = %d)",min_index,number_nodes));
    }
    Get(input,')');
}
//#####################################################################
// Function Write
//#####################################################################
template<int d> template<class RW> void SIMPLEX_MESH_BASE<d>::
Write(UNTYPED_OSTREAM<RW>& output) const
{
    Write_Binary_Group(output,number_nodes,d+1,elements);
}
//#####################################################################
template class SIMPLEX_MESH_BASE<0>;
template class SIMPLEX_MESH_BASE<1>;
template class SIMPLEX_MESH_BASE<2>;
template class SIMPLEX_MESH_BASE<3>;
template void SIMPLEX_MESH_BASE<2>::Simplices_On_Subsimplex(const VECTOR<int,2>&,ARRAY<int>&) const;
template void SIMPLEX_MESH_BASE<3>::Simplices_On_Subsimplex(const VECTOR<int,2>&,ARRAY<int>&) const;
template void SIMPLEX_MESH_BASE<3>::Simplices_On_Subsimplex(const VECTOR<int,3>&,ARRAY<int>&) const;
template void SIMPLEX_MESH_BASE<0>::Mark_Nodes_Referenced(ARRAY<int>&,const int&) const;
template void SIMPLEX_MESH_BASE<1>::Mark_Nodes_Referenced(ARRAY<int>&,const int&) const;
template void SIMPLEX_MESH_BASE<2>::Mark_Nodes_Referenced(ARRAY<int>&,const int&) const;
template void SIMPLEX_MESH_BASE<3>::Mark_Nodes_Referenced(ARRAY<int>&,const int&) const;
template void SIMPLEX_MESH_BASE<3>::Mark_Nodes_Referenced(ARRAY<bool>&,const bool&) const;
template void SIMPLEX_MESH_BASE<0>::Read(UNTYPED_ISTREAM<float>&);
template void SIMPLEX_MESH_BASE<0>::Read(UNTYPED_ISTREAM<char>&);
template void SIMPLEX_MESH_BASE<0>::Write(UNTYPED_OSTREAM<float>&) const;
template void SIMPLEX_MESH_BASE<0>::Write(UNTYPED_OSTREAM<char>&) const;
template void SIMPLEX_MESH_BASE<1>::Read(UNTYPED_ISTREAM<float>&);
template void SIMPLEX_MESH_BASE<1>::Read(UNTYPED_ISTREAM<char>&);
template void SIMPLEX_MESH_BASE<1>::Write(UNTYPED_OSTREAM<float>&) const;
template void SIMPLEX_MESH_BASE<1>::Write(UNTYPED_OSTREAM<char>&) const;
template void SIMPLEX_MESH_BASE<2>::Read(UNTYPED_ISTREAM<float>&);
template void SIMPLEX_MESH_BASE<2>::Read(UNTYPED_ISTREAM<char>&);
template void SIMPLEX_MESH_BASE<2>::Write(UNTYPED_OSTREAM<float>&) const;
template void SIMPLEX_MESH_BASE<2>::Write(UNTYPED_OSTREAM<char>&) const;
template void SIMPLEX_MESH_BASE<3>::Read(UNTYPED_ISTREAM<float>&);
template void SIMPLEX_MESH_BASE<3>::Read(UNTYPED_ISTREAM<char>&);
template void SIMPLEX_MESH_BASE<3>::Write(UNTYPED_OSTREAM<float>&) const;
template void SIMPLEX_MESH_BASE<3>::Write(UNTYPED_OSTREAM<char>&) const;
template void SIMPLEX_MESH_BASE<0>::Read(UNTYPED_ISTREAM<double>&);
template void SIMPLEX_MESH_BASE<0>::Write(UNTYPED_OSTREAM<double>&) const;
template void SIMPLEX_MESH_BASE<1>::Read(UNTYPED_ISTREAM<double>&);
template void SIMPLEX_MESH_BASE<1>::Write(UNTYPED_OSTREAM<double>&) const;
template void SIMPLEX_MESH_BASE<2>::Read(UNTYPED_ISTREAM<double>&);
template void SIMPLEX_MESH_BASE<2>::Write(UNTYPED_OSTREAM<double>&) const;
template void SIMPLEX_MESH_BASE<3>::Read(UNTYPED_ISTREAM<double>&);
template void SIMPLEX_MESH_BASE<3>::Write(UNTYPED_OSTREAM<double>&) const;
}
