//#####################################################################
// Copyright 2003-2007, Ron Fedkiw, Geoffrey Irving, Igor Neverov, Craig Schroeder, Andrew Selle, Tamar Shinar, 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 <Deformables/Forces/FINITE_VOLUME.h>
#include <Core/Arrays/ARRAYS_ND.h>
#include <Core/Arrays/INDIRECT_ARRAY.h>
#include <Core/Data_Structures/SPARSE_UNION_FIND.h>
#include <Core/Implicit_Objects/IMPLICIT_OBJECT.h>
#include <Core/Log/LOG.h>
#include <Core/Math_Tools/Robust_Arithmetic.h>
#include <Core/Math_Tools/Robust_Functions.h>
#include <Core/Matrices/DIAGONAL_MATRIX.h>
#include <Core/Matrices/MATRIX.h>
#include <Core/Matrices/SYMMETRIC_MATRIX.h>
#include <Core/Meshes/SEGMENT_MESH.h>
#include <Core/Meshes/TETRAHEDRALIZED_VOLUME.h>
#include <Core/Meshes/TRIANGULATED_AREA.h>
#include <Core/Meshes/TRIANGULATED_SURFACE.h>
#include <Deformables/Constitutive_Models/ANISOTROPIC_CONSTITUTIVE_MODEL.h>
#include <Deformables/Constitutive_Models/DIAGONALIZED_ISOTROPIC_STRESS_DERIVATIVE.h>
#include <Deformables/Constitutive_Models/DIAGONALIZED_STRESS_DERIVATIVE.h>
#include <Deformables/Constitutive_Models/ISOTROPIC_CONSTITUTIVE_MODEL.h>
#include <Deformables/Constitutive_Models/PLASTICITY_MODEL.h>
#include <Deformables/Forces/DIAGONALIZED_SEMI_IMPLICIT_ELEMENT.h>
#include <Deformables/Forces/STRAIN_MEASURE.h>
#include <Deformables/Particles/DEFORMABLE_PARTICLES.h>
namespace PhysBAM{
//#####################################################################
// Constructor
//#####################################################################
template<class TV,int d> FINITE_VOLUME<TV,d>::
FINITE_VOLUME(const bool use_uniform_density_input,STRAIN_MEASURE<TV,d>& strain_measure,CONSTITUTIVE_MODEL<T,d>& constitutive_model,PLASTICITY_MODEL<T,d>* plasticity_model)
    :DEFORMABLES_FORCES<TV>(dynamic_cast<DEFORMABLE_PARTICLES<TV>&>(strain_measure.particles)),strain_measure(strain_measure),constitutive_model(constitutive_model),
    plasticity_model(plasticity_model),dPi_dFe(0),dP_dFe(0),Be_scales_save(0),V(0),twice_max_strain_per_time_step(0),semi_implicit_data(0),use_uniform_density(use_uniform_density_input),
    density_list(0),implicit_surface(0),node_stiffness(0),edge_stiffness(0),force_segments(0),destroy_data(false)
{
    Update_Be_Scales();
    isotropic_model=dynamic_cast<ISOTROPIC_CONSTITUTIVE_MODEL<T,d>*>(&constitutive_model);
    anisotropic_model=dynamic_cast<ANISOTROPIC_CONSTITUTIVE_MODEL<T,d>*>(&constitutive_model);
    if(anisotropic_model) Save_V();
    strain_measure.mesh.Initialize_Incident_Elements();
    if(use_uniform_density){
        T total_mass=0;
        ARRAY<int> mesh_particles=Get_Unique(strain_measure.mesh_object.mesh.elements.Flattened());
        for(int i=0;i<mesh_particles.m;i++) total_mass+=particles.mass(mesh_particles(i));
        density=total_mass/strain_measure.mesh_object.Total_Size();
        if(density==0) density=TV::m==1?1:TV::m==2?100:1000;}
    else{
        density_list=new ARRAY<T>(strain_measure.mesh_object.mesh.elements.m);
        for(int i=0;i<strain_measure.mesh.elements.m;i++){
            const VECTOR<int,d+1>& nodes=strain_measure.mesh.elements(i);
            T volume=strain_measure.mesh_object.Signed_Size(i);
            for(int j=0;j<nodes.m;j++) (*density_list)(i)+=particles.mass(nodes(j))/(*strain_measure.mesh.incident_elements)(nodes(j)).m/volume;
            if((*density_list)(i)==0) (*density_list)(i)=TV::m==1?1:TV::m==2?100:1000;}}
}
//#####################################################################
// Destructor
//#####################################################################
template<class TV,int d> FINITE_VOLUME<TV,d>::
~FINITE_VOLUME()
{
    delete dPi_dFe;delete dP_dFe;delete Be_scales_save;delete V;delete semi_implicit_data;delete node_stiffness;delete edge_stiffness;delete density_list;
    if(destroy_data){delete &strain_measure;delete &constitutive_model;delete plasticity_model;}
}
//#####################################################################
// Function Create
//#####################################################################
template<class TV,int d> FINITE_VOLUME<TV,d>* FINITE_VOLUME<TV,d>::
Create(T_OBJECT& object,CONSTITUTIVE_MODEL<T,d>* constitutive_model,PLASTICITY_MODEL<T,d>* plasticity_model,const bool verbose,const bool use_uniform_density)
{
    STRAIN_MEASURE<TV,d>* strain_measure=new STRAIN_MEASURE<TV,d>(object);
    if(verbose) strain_measure->Print_Altitude_Statistics();
    FINITE_VOLUME* fvm=new FINITE_VOLUME(use_uniform_density,*strain_measure,*constitutive_model,plasticity_model);
    fvm->destroy_data=true;return fvm;
}
//#####################################################################
// Function Save_V
//#####################################################################
template<class TV,int d> void FINITE_VOLUME<TV,d>::
Save_V()
{
    if(!V) V=new ARRAY<MATRIX<T,d> >();
}
//#####################################################################
// Function Save_Stress_Derivative
//#####################################################################
template<class TV,int d> void FINITE_VOLUME<TV,d>::
Save_Stress_Derivative()
{
    if(isotropic_model || anisotropic_model->use_isotropic_component_of_stress_derivative_only){
        if(!dPi_dFe) dPi_dFe=new ARRAY<DIAGONALIZED_ISOTROPIC_STRESS_DERIVATIVE<VECTOR<T,d> > >();
        delete dP_dFe;
        dP_dFe=0;}
    else{
        if(!dP_dFe) dP_dFe=new ARRAY<DIAGONALIZED_STRESS_DERIVATIVE<T,d> >();
        delete dPi_dFe;
        dPi_dFe=0;}
}
//#####################################################################
// Function Use_Quasistatics
//#####################################################################
template<class TV,int d> void FINITE_VOLUME<TV,d>::
Use_Quasistatics()
{
    Save_V();Save_Stress_Derivative();
}
//#####################################################################
// Function Use_Stiffness_Matrix
//#####################################################################
template<class TV,int d> void FINITE_VOLUME<TV,d>::
Use_Stiffness_Matrix()
{
    node_stiffness=new ARRAY<SYMMETRIC_MATRIX<T,TV::m> >();edge_stiffness=new ARRAY<MATRIX<T,TV::m> >();
    if(!strain_measure.mesh.segment_mesh) strain_measure.mesh.Initialize_Segment_Mesh();
    if(!strain_measure.mesh.element_edges) strain_measure.mesh.Initialize_Element_Edges();
    if(force_segments) delete force_segments;
    force_segments=new ARRAY<int>();
}
//#####################################################################
// Function Update_Be_Scales
//#####################################################################
template<class TV,int d> void FINITE_VOLUME<TV,d>::
Update_Be_Scales()
{
    Be_scales.Resize(strain_measure.Dm_inverse.m,no_init);
    for(int t=0;t<Be_scales.m;t++) Be_scales(t)=-(T)1/factorial(d)/strain_measure.Dm_inverse(t).Determinant();
}
//#####################################################################
// Function Initialize_Be_Scales_Save
//#####################################################################
template<class TV,int d> void FINITE_VOLUME<TV,d>::
Initialize_Be_Scales_Save()
{
    Be_scales_save=new ARRAY<T>(Be_scales);
}
//#####################################################################
// Function Copy_Be_Scales_Save_Into_Be_Scales
//#####################################################################
template<class TV,int d> void FINITE_VOLUME<TV,d>::
Copy_Be_Scales_Save_Into_Be_Scales(const ARRAY<int>& map)
{
    Be_scales=Be_scales_save->Subset(map);
}
//#####################################################################
// Function Add_Dependencies
//#####################################################################
template<class TV,int d> void FINITE_VOLUME<TV,d>::
Add_Dependencies(SEGMENT_MESH& dependency_mesh) const
{
    strain_measure.mesh.Add_Dependencies(dependency_mesh);
}
//#####################################################################
// Function Update_Mpi
//#####################################################################
template<class TV,int d> void FINITE_VOLUME<TV,d>::
Update_Mpi(const ARRAY<bool>& particle_is_simulated)
{
    if(!strain_measure.mesh.segment_mesh) strain_measure.mesh.Initialize_Segment_Mesh();
    Update_Force_Elements(force_elements,strain_measure.mesh.elements,particle_is_simulated);
    if(force_segments) Update_Force_Elements(*force_segments,strain_measure.mesh.segment_mesh->elements,particle_is_simulated);
    Update_Force_Particles(force_particles,
        strain_measure.mesh.elements.Flattened(),false,
        particle_is_simulated,true);
}
//#####################################################################
// Function Update_Position_Based_State
//#####################################################################
namespace{
template<class T,int d> void Set_Inversion_Based_On_Implicit_Surface(FINITE_VOLUME<VECTOR<T,d>,d>&,const IMPLICIT_OBJECT<VECTOR<T,d> >&,const int)
{}
template<class T> void Set_Inversion_Based_On_Implicit_Surface(FINITE_VOLUME<VECTOR<T,3>,2>& fvm,const IMPLICIT_OBJECT<VECTOR<T,3> >& implicit_surface,const int triangle)
{
    VECTOR<T,3> normal=implicit_surface.Normal(fvm.strain_measure.mesh_object.Centroid(triangle));
    if(normal.Dot(Weighted_Normal(fvm.U(triangle)))<0){
        fvm.Fe_hat(triangle).x.y*=-1;fvm.U(triangle).Set_Column(1,-fvm.U(triangle).Column(1));}
}}
template<class TV,int d> void FINITE_VOLUME<TV,d>::
Update_Position_Based_State(const T time,const bool is_position_update,const bool update_hessian)
{
    Save_Stress_Derivative();
    if(anisotropic_model && !V) PHYSBAM_FATAL_ERROR();
    int elements=strain_measure.Dm_inverse.m;
    U.Resize(elements,no_init);
    De_inverse_hat.Resize(elements,no_init);
    Fe_hat.Resize(elements,no_init);
    if(dPi_dFe) dPi_dFe->Resize(elements,no_init);
    if(dP_dFe) dP_dFe->Resize(elements,no_init);
    if(V) V->Resize(elements,no_init);
    if(node_stiffness)
        node_stiffness->Resize(strain_measure.mesh_object.particles.Size(),init_all);
    if(edge_stiffness)
        edge_stiffness->Resize(strain_measure.mesh.segment_mesh->elements.m,init_all);
    if(!plasticity_model){
        MATRIX<MATRIX<T,TV::m>,d+1,d+1> dfdx;
        for(int t:force_elements){
            auto [U_,F_,V_local]=strain_measure.F(t).Singular_Value_Decomposition();
            U(t)=U_;
            Fe_hat(t)=F_;
            if(implicit_surface) Set_Inversion_Based_On_Implicit_Surface(*this,*implicit_surface,t);
            if(dPi_dFe) constitutive_model.Isotropic_Stress_Derivative(Fe_hat(t),(*dPi_dFe)(t),t);
            if(dP_dFe) anisotropic_model->Stress_Derivative(Fe_hat(t),V_local,(*dP_dFe)(t),t);
            De_inverse_hat(t)=strain_measure.Dm_inverse(t)*V_local;
            if(V) (*V)(t)=V_local;
            if(node_stiffness && edge_stiffness){
                T_MATRIX U_scale=U(t)*Be_scales(t);
                for(int k=0;k<TV::m;k++)
                    for(int l=0;l<d;l++){
                        T_MATRIX dDs,dG;dDs(k,l)=(T)1;
                        if(!dPi_dFe && !dP_dFe) // precompute damping matrix
                            dG=U_scale*constitutive_model.P_From_Strain_Rate(Fe_hat(t),U(t).Transpose_Times(dDs)*De_inverse_hat(t),t).Times_Transpose(De_inverse_hat(t));
                        else if(isotropic_model) // precompute stiffness matrix
                            dG=U_scale*(*dPi_dFe)(t).Differential(U(t).Transpose_Times(dDs)*De_inverse_hat(t)).Times_Transpose(De_inverse_hat(t));
                        else if(dP_dFe)
                            dG=U_scale*anisotropic_model->dP_From_dF(U(t).Transpose_Times(dDs)*De_inverse_hat(t),Fe_hat(t),(*V)(t),(*dP_dFe)(t),t).Times_Transpose(De_inverse_hat(t));
                        else
                            dG=U_scale*anisotropic_model->dP_From_dF(U(t).Transpose_Times(dDs)*De_inverse_hat(t),Fe_hat(t),(*V)(t),(*dPi_dFe)(t),t).Times_Transpose(De_inverse_hat(t));
                        for(int i=0;i<TV::m;i++) for(int j=0;j<d;j++) dfdx(l+1,j+1)(k,i)=dG(i,j);}
                for(int i=1;i<d+1;i++){dfdx(i,0)=MATRIX<T,TV::m>();for(int j=1;j<d+1;j++) dfdx(i,0)-=dfdx(i,j);}
                for(int j=0;j<d+1;j++){dfdx(0,j)=MATRIX<T,TV::m>();for(int i=1;i<d+1;i++) dfdx(0,j)-=dfdx(i,j);}
                VECTOR<int,d+1> nodes=strain_measure.mesh.elements(t);
                for(int v=0;v<nodes.m;v++) (*node_stiffness)(nodes[v])+=dfdx(v,v).Symmetric_Part();
                VECTOR<int,d*(d+1)/2> edges=(*strain_measure.mesh.element_edges)(t);
                for(int e=0;e<edges.m;e++){
                    int i,j;strain_measure.mesh.segment_mesh->elements(edges[e]).Get(i,j);
                    int m=nodes.Find(i),n=nodes.Find(j);assert(m!=n);
                    (*edge_stiffness)(edges[e])+=dfdx(m,n);}}}}
    else{
        if(implicit_surface) PHYSBAM_NOT_IMPLEMENTED();
        if(anisotropic_model) PHYSBAM_NOT_IMPLEMENTED();
        DIAGONAL_MATRIX<T,d> Fe_project_hat;
        for(int t:force_elements){
            T_MATRIX F=strain_measure.F(t);
            auto [U_,F_,V_local]=(F*plasticity_model->Fp_inverse(t)).Singular_Value_Decomposition();
            U(t)=U_;
            Fe_hat(t)=F_;
            if(implicit_surface) Set_Inversion_Based_On_Implicit_Surface(*this,*implicit_surface,t);
            if(plasticity_model->Project_Fe(Fe_hat(t),Fe_project_hat)){
                plasticity_model->Project_Fp(t,Fe_project_hat.Inverse()*U(t).Transpose_Times(F));
                auto [U_,F_,V_]=(F*plasticity_model->Fp_inverse(t)).Singular_Value_Decomposition();
                U(t)=U_;
                Fe_hat(t)=F_;
                V_local=V_;}
            De_inverse_hat(t)=strain_measure.Dm_inverse(t)*plasticity_model->Fp_inverse(t)*V_local;
            Be_scales(t)=-(T)1/factorial(d)/De_inverse_hat(t).Determinant();}}
}
//#####################################################################
// Function Add_Forces
//#####################################################################
template<class TV,int d> void FINITE_VOLUME<TV,d>::
Add_Forces(ARRAY_VIEW<TV> F,const T time) const
{
    if(anisotropic_model)
        for(int t:force_elements){
            T_MATRIX forces=U(t)*(Be_scales(t)*anisotropic_model->P_From_Strain(Fe_hat(t),(*V)(t),t)).Times_Transpose(De_inverse_hat(t));
            strain_measure.Distribute_Force(F,t,forces);}
    else
        for(int t:force_elements){
            T_MATRIX forces=U(t)*(Be_scales(t)*isotropic_model->P_From_Strain(Fe_hat(t),t)).Times_Transpose(De_inverse_hat(t));
            strain_measure.Distribute_Force(F,t,forces);}
}
//#####################################################################
// Function Add_Velocity_Dependent_Forces
//#####################################################################
template<class TV,int d> void FINITE_VOLUME<TV,d>::
Add_Velocity_Dependent_Forces(ARRAY_VIEW<const TV> V,ARRAY_VIEW<TV> F,const T time) const
{
    //force_elements.Print();
    if(node_stiffness && edge_stiffness && !dPi_dFe && !dP_dFe){
        if(dP_dFe || dPi_dFe) PHYSBAM_FATAL_ERROR();
        for(int p:force_particles){
            F(p)+=(*node_stiffness)(p)*V(p);}
        for(int e:*force_segments){
            int m,n;strain_measure.mesh.segment_mesh->elements(e).Get(m,n);
            F(m)+=(*edge_stiffness)(e)*V(n);F(n)+=(*edge_stiffness)(e).Transpose_Times(V(m));}}
    else{
        for(int t:force_elements){
            MATRIX<T,d> Fe_dot_hat=U(t).Transpose_Times(strain_measure.Ds(V,t))*De_inverse_hat(t);
            T_MATRIX forces=U(t)*(Be_scales(t)*constitutive_model.P_From_Strain_Rate(Fe_hat(t),Fe_dot_hat,t)).Times_Transpose(De_inverse_hat(t));
            strain_measure.Distribute_Force(F,t,forces);}}
}
//#####################################################################
// Function Enforce_Definiteness
//#####################################################################
template<class TV,int d> void FINITE_VOLUME<TV,d>::
Enforce_Definiteness(const bool enforce_definiteness_input)
{
    constitutive_model.enforce_definiteness=enforce_definiteness_input;
}
//#####################################################################
// Function Intialize_CFL
//#####################################################################
template<class TV,int d> void FINITE_VOLUME<TV,d>::
Initialize_CFL(ARRAY_VIEW<FREQUENCY_DATA> frequency)
{
    T one_over_cfl_number=1/cfl_number,one_over_cfl_number_squared=sqr(one_over_cfl_number);
    HASHTABLE<int,FREQUENCY_DATA> particle_frequency;
    for(int t:force_elements){
        T one_over_altitude_squared_and_density=(T)1/(sqr(strain_measure.Rest_Altitude(t))*(use_uniform_density?density:(*density_list)(t)));
        T elastic_squared=constitutive_model.Maximum_Elastic_Stiffness(t)*one_over_altitude_squared_and_density*one_over_cfl_number_squared;
        T damping=constitutive_model.Maximum_Damping_Stiffness(t)*one_over_altitude_squared_and_density*one_over_cfl_number;
        const VECTOR<int,d+1>& nodes=strain_measure.mesh.elements(t);
        for(int j=0;j<nodes.m;j++){FREQUENCY_DATA& frequency_data=particle_frequency.Get_Or_Insert(nodes[j]);
            frequency_data.elastic_squared=max(frequency_data.elastic_squared,elastic_squared);
            frequency_data.damping=max(frequency_data.damping,damping);}}
    for(auto& p:particle_frequency){
        frequency(p.first).elastic_squared+=p.second.elastic_squared;
        frequency(p.first).damping+=p.second.damping;}
}
//#####################################################################
// Function CFL_Strain_Rate
//#####################################################################
template<class TV,int d> typename TV::SCALAR FINITE_VOLUME<TV,d>::
CFL_Strain_Rate() const
{
    T max_strain_rate=0;
    for(int t:force_elements){
        max_strain_rate=max(max_strain_rate,strain_measure.Velocity_Gradient(t).Max_Abs());}
    return Robust_Divide(max_strain_per_time_step,max_strain_rate);
}
//#####################################################################
// Function Potential_Energy
//#####################################################################
template<class TV,int d> typename TV::SCALAR FINITE_VOLUME<TV,d>::
Potential_Energy(const T time) const
{
    T potential_energy=0;
    if(anisotropic_model) PHYSBAM_NOT_IMPLEMENTED();
    else
        for(int t:force_elements){
            potential_energy-=Be_scales(t)*isotropic_model->Energy_Density(Fe_hat(t),t);}
    return potential_energy;
}
//#####################################################################
// Function Add_Force_Derivatives
//#####################################################################
template<class TV,int d> void FINITE_VOLUME<TV,d>::
Add_Force_Derivatives(ARRAY_VIEW<const TV> VV,ARRAY_VIEW<TV> F,const T time,bool transpose) const
{
    if(node_stiffness && edge_stiffness){
        for(int p:force_particles){
            F(p)+=(*node_stiffness)(p)*VV(p);}
        for(int e:*force_segments){
            int m,n;strain_measure.mesh.segment_mesh->elements(e).Get(m,n);
            F(m)+=(*edge_stiffness)(e)*VV(n);F(n)+=(*edge_stiffness)(e).Transpose_Times(VV(m));}}
    else if(anisotropic_model){
        if(!dPi_dFe && !dP_dFe) PHYSBAM_FATAL_ERROR();
        for(int t:force_elements){
            T_MATRIX dDs=strain_measure.Ds(VV,t),dG;
            if(dP_dFe) dG=U(t)*(Be_scales(t)*anisotropic_model->dP_From_dF(U(t).Transpose_Times(dDs)*De_inverse_hat(t),Fe_hat(t),(*V)(t),(*dP_dFe)(t),t)).Times_Transpose(De_inverse_hat(t));
            else dG=U(t)*(Be_scales(t)*anisotropic_model->dP_From_dF(U(t).Transpose_Times(dDs)*De_inverse_hat(t),Fe_hat(t),(*V)(t),(*dPi_dFe)(t),t)).Times_Transpose(De_inverse_hat(t));
            strain_measure.Distribute_Force(F,t,dG);}}
    else{
        if(!dPi_dFe) PHYSBAM_FATAL_ERROR();
        for(int t:force_elements){
            T_MATRIX dDs=strain_measure.Ds(VV,t);
            T_MATRIX dG=U(t)*(Be_scales(t)*(*dPi_dFe)(t).Differential(U(t).Transpose_Times(dDs)*De_inverse_hat(t))).Times_Transpose(De_inverse_hat(t));
            strain_measure.Distribute_Force(F,t,dG);}}
}
//#####################################################################
// Function Compute_Sound_Speed
//#####################################################################
template<class TV,int d> auto FINITE_VOLUME<TV,d>::
Compute_Sound_Speed() const -> T
{
    assert(!anisotropic_model && isotropic_model);
    ARRAY<int> counts(particles.number);
    counts.Subset(strain_measure.mesh.elements.Flattened())+=1;
    T speed=0;
    for(int t:force_elements)
    {
        T mass=0;
        for(int p:strain_measure.mesh.elements(t))
            mass+=particles.mass(p)/counts(p);
        T density=mass/strain_measure.mesh_object.Element_Size(t);
        speed=max(speed,isotropic_model->Sound_Speed(Fe_hat(t),density,t));
    }
    return speed;
}

template<class T_OBJECT,class T,int d> FINITE_VOLUME<typename T_OBJECT::VECTOR_T,d>*
Create_Finite_Volume(T_OBJECT& object,CONSTITUTIVE_MODEL<T,d>* constitutive_model,const bool limit_time_step_by_strain_rate,const T max_strain_per_time_step,
    const bool use_rest_state_for_strain_rate,const bool precompute_stiffness_matrix,const bool verbose,PLASTICITY_MODEL<T,d>* plasticity_model,const bool use_uniform_density)
{
    typedef typename T_OBJECT::VECTOR_T TV;
    FINITE_VOLUME<TV,d>* fvm=FINITE_VOLUME<TV,d>::Create(object,constitutive_model,plasticity_model,verbose,use_uniform_density);
    fvm->Limit_Time_Step_By_Strain_Rate(limit_time_step_by_strain_rate,max_strain_per_time_step);
    fvm->Use_Rest_State_For_Strain_Rate(use_rest_state_for_strain_rate);
    if(precompute_stiffness_matrix) fvm->Use_Stiffness_Matrix();
    return fvm;
}

template<class T_OBJECT,class T,int d> FINITE_VOLUME<typename T_OBJECT::VECTOR_T,d>*
Create_Quasistatic_Finite_Volume(T_OBJECT& object,CONSTITUTIVE_MODEL<T,d>* constitutive_model,const bool precompute_stiffness_matrix,const bool verbose,const bool use_uniform_density)
{
    typedef typename T_OBJECT::VECTOR_T TV;
    FINITE_VOLUME<TV,d>* fvm=Create_Finite_Volume(object,constitutive_model,true,(T).1,true,precompute_stiffness_matrix,verbose,(PLASTICITY_MODEL<T,d>*)0,use_uniform_density);
    fvm->Use_Quasistatics();if(precompute_stiffness_matrix) fvm->Use_Stiffness_Matrix();return fvm;
}
//#####################################################################
template class FINITE_VOLUME<VECTOR<float,2>,2>;
template class FINITE_VOLUME<VECTOR<float,3>,2>;
template class FINITE_VOLUME<VECTOR<float,3>,3>;
template class FINITE_VOLUME<VECTOR<double,2>,2>;
template class FINITE_VOLUME<VECTOR<double,3>,2>;
template class FINITE_VOLUME<VECTOR<double,3>,3>;
template FINITE_VOLUME<TETRAHEDRALIZED_VOLUME<double>::VECTOR_T,3>* Create_Finite_Volume<TETRAHEDRALIZED_VOLUME<double>,double,3>(TETRAHEDRALIZED_VOLUME<double>&,CONSTITUTIVE_MODEL<double,3>*,bool,double,bool,bool,bool,PLASTICITY_MODEL<double,3>*,bool);
template FINITE_VOLUME<TETRAHEDRALIZED_VOLUME<double>::VECTOR_T,3>* Create_Quasistatic_Finite_Volume<TETRAHEDRALIZED_VOLUME<double>,double,3>(TETRAHEDRALIZED_VOLUME<double>&,CONSTITUTIVE_MODEL<double,3>*,bool,bool,bool);
template FINITE_VOLUME<TETRAHEDRALIZED_VOLUME<float>::VECTOR_T,3>* Create_Finite_Volume<TETRAHEDRALIZED_VOLUME<float>,float,3>(TETRAHEDRALIZED_VOLUME<float>&,CONSTITUTIVE_MODEL<float,3>*,bool,float,bool,bool,bool,PLASTICITY_MODEL<float,3>*,bool);
template FINITE_VOLUME<TETRAHEDRALIZED_VOLUME<float>::VECTOR_T,3>* Create_Quasistatic_Finite_Volume<TETRAHEDRALIZED_VOLUME<float>,float,3>(TETRAHEDRALIZED_VOLUME<float>&,CONSTITUTIVE_MODEL<float,3>*,bool,bool,bool);
template FINITE_VOLUME<TRIANGULATED_AREA<double>::VECTOR_T,2>* Create_Finite_Volume<TRIANGULATED_AREA<double>,double,2>(TRIANGULATED_AREA<double>&,CONSTITUTIVE_MODEL<double,2>*,bool,double,bool,bool,bool,PLASTICITY_MODEL<double,2>*,bool);
template FINITE_VOLUME<TRIANGULATED_AREA<double>::VECTOR_T,2>* Create_Quasistatic_Finite_Volume<TRIANGULATED_AREA<double>,double,2>(TRIANGULATED_AREA<double>&,CONSTITUTIVE_MODEL<double,2>*,bool,bool,bool);
template FINITE_VOLUME<TRIANGULATED_AREA<float>::VECTOR_T,2>* Create_Finite_Volume<TRIANGULATED_AREA<float>,float,2>(TRIANGULATED_AREA<float>&,CONSTITUTIVE_MODEL<float,2>*,bool,float,bool,bool,bool,PLASTICITY_MODEL<float,2>*,bool);
template FINITE_VOLUME<TRIANGULATED_AREA<float>::VECTOR_T,2>* Create_Quasistatic_Finite_Volume<TRIANGULATED_AREA<float>,float,2>(TRIANGULATED_AREA<float>&,CONSTITUTIVE_MODEL<float,2>*,bool,bool,bool);
template FINITE_VOLUME<TRIANGULATED_SURFACE<double>::VECTOR_T,2>* Create_Finite_Volume<TRIANGULATED_SURFACE<double>,double,2>(TRIANGULATED_SURFACE<double>&,CONSTITUTIVE_MODEL<double,2>*,bool,double,bool,bool,bool,PLASTICITY_MODEL<double,2>*,bool);
template FINITE_VOLUME<TRIANGULATED_SURFACE<float>::VECTOR_T,2>* Create_Finite_Volume<TRIANGULATED_SURFACE<float>,float,2>(TRIANGULATED_SURFACE<float>&,CONSTITUTIVE_MODEL<float,2>*,bool,float,bool,bool,bool,PLASTICITY_MODEL<float,2>*,bool);
}
