| //##################################################################### |
| // Copyright 2003-2004, Ron Fedkiw, Geoffrey Irving, Neil Molino, 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. |
| //##################################################################### |
| // Class DIAGONALIZED_FINITE_VOLUME_3D |
| //##################################################################### |
| #ifndef __DIAGONALIZED_FINITE_VOLUME_3D__ |
| #define __DIAGONALIZED_FINITE_VOLUME_3D__ |
| |
| #include "SOLIDS_FORCES.h" |
| #include "../Constitutive_Models/STRAIN_MEASURE_3D.h" |
| #include "../Constitutive_Models/DIAGONALIZED_CONSTITUTIVE_MODEL_3D.h" |
| #include "DIAGONALIZED_SEMI_IMPLICIT_ELEMENT_3D.h" |
| #include "../Grids/SEGMENT_MESH.h" |
| #include "../Thread_Utilities/THREAD_DIVISION_PARAMETERS.h" |
| |
| //#define READ_THREADING_AUXILIARY_STRUCTURES_FROM_SNAPSHOT |
| |
| namespace PhysBAM |
| { |
| |
| template<class T> |
| class DIAGONALIZED_FINITE_VOLUME_3D_THREADING_AUXILIARY_STRUCTURES |
| { |
| public: |
| LIST_ARRAYS<int>* extended_edges; |
| LIST_ARRAYS<int>* extended_tetrahedrons; |
| LIST_ARRAYS<int>* extended_tetrahedron_extended_edges; |
| LIST_ARRAY<MATRIX_3X3<T> >* extended_edge_stiffness; |
| ARRAY<VECTOR_2D<int> >* node_ranges; |
| ARRAY<VECTOR_2D<int> >* internal_edge_ranges; |
| ARRAY<VECTOR_2D<int> >* external_edge_ranges; |
| ARRAY<VECTOR_2D<int> >* extended_tetrahedron_ranges; |
| LIST_ARRAY<int>* extended_tetrahedron_parents; |
| LIST_ARRAY<MATRIX_3X3<T> >* extended_U; |
| LIST_ARRAY<MATRIX_3X3<T> >* extended_De_inverse_hat; |
| LIST_ARRAY<DIAGONAL_MATRIX_3X3<T> >* extended_Fe_hat; |
| LIST_ARRAY<DIAGONALIZED_ISOTROPIC_STRESS_DERIVATIVE_3D<T> >* extended_dP_dFe; |
| LIST_ARRAY<MATRIX_3X3<T> >* extended_V; |
| |
| DIAGONALIZED_FINITE_VOLUME_3D_THREADING_AUXILIARY_STRUCTURES() |
| : extended_edges (0), extended_tetrahedrons (0), extended_tetrahedron_extended_edges (0), extended_edge_stiffness (0), node_ranges (0), internal_edge_ranges (0), external_edge_ranges (0), |
| extended_tetrahedron_ranges (0), extended_tetrahedron_parents (0), extended_U (0), extended_De_inverse_hat (0), extended_Fe_hat (0), extended_dP_dFe (0), extended_V (0) |
| {} |
| |
| ~DIAGONALIZED_FINITE_VOLUME_3D_THREADING_AUXILIARY_STRUCTURES() |
| { |
| delete extended_edges; |
| delete extended_tetrahedrons; |
| delete extended_tetrahedron_extended_edges; |
| delete extended_edge_stiffness; |
| delete node_ranges; |
| delete internal_edge_ranges; |
| delete external_edge_ranges; |
| delete extended_tetrahedron_ranges; |
| delete extended_tetrahedron_parents; |
| delete extended_U; |
| delete extended_De_inverse_hat; |
| delete extended_Fe_hat; |
| delete extended_dP_dFe; |
| delete extended_V; |
| } |
| |
| void Allocate() |
| { |
| extended_edges = new LIST_ARRAYS<int>; |
| extended_tetrahedrons = new LIST_ARRAYS<int>; |
| extended_tetrahedron_extended_edges = new LIST_ARRAYS<int>; |
| extended_edge_stiffness = new LIST_ARRAY<MATRIX_3X3<T> >; |
| node_ranges = new ARRAY<VECTOR_2D<int> >; |
| internal_edge_ranges = new ARRAY<VECTOR_2D<int> >; |
| external_edge_ranges = new ARRAY<VECTOR_2D<int> >; |
| extended_tetrahedron_ranges = new ARRAY<VECTOR_2D<int> >; |
| extended_tetrahedron_parents = new LIST_ARRAY<int>; |
| extended_U = new LIST_ARRAY<MATRIX_3X3<T> >; |
| extended_De_inverse_hat = new LIST_ARRAY<MATRIX_3X3<T> >; |
| extended_Fe_hat = new LIST_ARRAY<DIAGONAL_MATRIX_3X3<T> >; |
| extended_dP_dFe = new LIST_ARRAY<DIAGONALIZED_ISOTROPIC_STRESS_DERIVATIVE_3D<T> >; |
| extended_V = new LIST_ARRAY<MATRIX_3X3<T> >; |
| } |
| |
| template<class RW> |
| void Read (std::istream& input_stream) |
| { |
| extended_edges->template Read<RW> (input_stream); |
| extended_tetrahedrons->template Read<RW> (input_stream); |
| extended_tetrahedron_extended_edges->template Read<RW> (input_stream); |
| node_ranges->template Read<RW> (input_stream); |
| internal_edge_ranges->template Read<RW> (input_stream); |
| external_edge_ranges->template Read<RW> (input_stream); |
| extended_tetrahedron_ranges->template Read<RW> (input_stream); |
| extended_tetrahedron_parents->template Read<RW> (input_stream); |
| } |
| |
| template<class RW> |
| void Write (std::ostream& output_stream) const |
| { |
| extended_edges->template Write<RW> (output_stream); |
| extended_tetrahedrons->template Write<RW> (output_stream); |
| extended_tetrahedron_extended_edges->template Write<RW> (output_stream); |
| node_ranges->template Write<RW> (output_stream); |
| internal_edge_ranges->template Write<RW> (output_stream); |
| external_edge_ranges->template Write<RW> (output_stream); |
| extended_tetrahedron_ranges->template Write<RW> (output_stream); |
| extended_tetrahedron_parents->template Write<RW> (output_stream); |
| } |
| }; |
| |
| template<class T> |
| class DIAGONALIZED_FINITE_VOLUME_3D: public SOLIDS_FORCES<T, VECTOR_3D<T> > |
| { |
| public: |
| using SOLIDS_FORCES<T, VECTOR_3D<T> >::particles; |
| using SOLIDS_FORCES<T, VECTOR_3D<T> >::CFL_initialized; |
| using SOLIDS_FORCES<T, VECTOR_3D<T> >::CFL_elastic_time_step; |
| using SOLIDS_FORCES<T, VECTOR_3D<T> >::CFL_damping_time_step; |
| using SOLIDS_FORCES<T, VECTOR_3D<T> >::max_strain_per_time_step; |
| |
| STRAIN_MEASURE_3D<T>& strain_measure; |
| DIAGONALIZED_CONSTITUTIVE_MODEL_3D<T>& constitutive_model; |
| LIST_ARRAY<T> Be_scales; |
| LIST_ARRAY<MATRIX_3X3<T> > U; |
| LIST_ARRAY<MATRIX_3X3<T> > De_inverse_hat; |
| LIST_ARRAY<DIAGONAL_MATRIX_3X3<T> > Fe_hat; |
| LIST_ARRAY<DIAGONALIZED_ISOTROPIC_STRESS_DERIVATIVE_3D<T> >* dP_dFe; |
| LIST_ARRAY<T>* Be_scales_save; |
| LIST_ARRAY<MATRIX_3X3<T> >* V; |
| T twice_max_strain_per_time_step; // for asynchronous |
| LIST_ARRAY<DIAGONALIZED_SEMI_IMPLICIT_ELEMENT_3D<T> >* semi_implicit_data; |
| LIST_ARRAY<SYMMETRIC_MATRIX_3X3<T> >* node_stiffness; |
| LIST_ARRAY<MATRIX_3X3<T> >* edge_stiffness; |
| DIAGONALIZED_FINITE_VOLUME_3D_THREADING_AUXILIARY_STRUCTURES<T>* threading_auxiliary_structures; |
| |
| DIAGONALIZED_FINITE_VOLUME_3D (STRAIN_MEASURE_3D<T>& strain_measure_input, DIAGONALIZED_CONSTITUTIVE_MODEL_3D<T>& constitutive_model_input) |
| : SOLIDS_FORCES<T, VECTOR_3D<T> > (strain_measure_input.tetrahedralized_volume.particles), strain_measure (strain_measure_input), constitutive_model (constitutive_model_input), |
| dP_dFe (0), Be_scales_save (0), V (0), semi_implicit_data (0), node_stiffness (0), edge_stiffness (0), threading_auxiliary_structures (0) |
| { |
| Update_Be_Scales(); |
| |
| if (constitutive_model.anisotropic) Save_V(); |
| } |
| |
| ~DIAGONALIZED_FINITE_VOLUME_3D() |
| { |
| delete dP_dFe; |
| delete Be_scales_save; |
| delete V; |
| delete semi_implicit_data; |
| delete node_stiffness; |
| delete edge_stiffness; |
| delete threading_auxiliary_structures; |
| } |
| |
| void Save_V() |
| { |
| if (!V) V = new LIST_ARRAY<MATRIX_3X3<T> >; |
| } |
| |
| void Save_Stress_Derivative() |
| { |
| if (!dP_dFe) dP_dFe = new LIST_ARRAY<DIAGONALIZED_ISOTROPIC_STRESS_DERIVATIVE_3D<T> >; |
| } |
| |
| void Use_Quasistatics() |
| { |
| Save_V(); |
| Save_Stress_Derivative(); |
| } |
| |
| void Use_Stiffness_Matrix() |
| { |
| node_stiffness = new LIST_ARRAY<SYMMETRIC_MATRIX_3X3<T> >; |
| |
| if (!THREAD_DIVISION_PARAMETERS<T>::Thread_Divisions_Enabled()) edge_stiffness = new LIST_ARRAY<MATRIX_3X3<T> >; |
| |
| if (!strain_measure.tetrahedron_mesh.segment_mesh) strain_measure.tetrahedron_mesh.Initialize_Segment_Mesh(); |
| |
| if (!strain_measure.tetrahedron_mesh.tetrahedron_edges) strain_measure.tetrahedron_mesh.Initialize_Tetrahedron_Edges(); |
| |
| if (THREAD_DIVISION_PARAMETERS<T>::Thread_Divisions_Enabled()) |
| #ifdef READ_THREADING_AUXILIARY_STRUCTURES_FROM_SNAPSHOT |
| Read_Threading_Auxiliary_Structures(); |
| |
| #else |
| Update_Threading_Auxiliary_Structures(); |
| #endif |
| } |
| |
| void Enforce_Definiteness (const bool enforce_definiteness_input) |
| { |
| constitutive_model.enforce_definiteness = enforce_definiteness_input; |
| } |
| |
| void Update_Be_Scales() |
| { |
| Be_scales.Resize_Array (strain_measure.Dm_inverse.m); |
| |
| for (int t = 1; t <= Be_scales.m; t++) Be_scales (t) = - (T) one_sixth / strain_measure.Dm_inverse (t).Determinant(); |
| } |
| |
| void Initialize_Be_Scales_Save() |
| { |
| Be_scales_save = new LIST_ARRAY<T> (Be_scales.m); |
| LIST_ARRAY<T>::copy (Be_scales, *Be_scales_save); |
| } |
| |
| void Copy_Be_Scales_Save_Into_Be_Scales (const LIST_ARRAY<int>& map) |
| { |
| Be_scales.Resize_Array (map.m); |
| |
| for (int t = 1; t <= map.m; t++) Be_scales (t) = (*Be_scales_save) (map (t)); |
| } |
| |
| //##################################################################### |
| void Update_Threading_Auxiliary_Structures(); |
| void Read_Threading_Auxiliary_Structures(); |
| void Update_Position_Based_State(); |
| static void Update_Position_Based_State_Helper (long thread_id, void* helper_raw); |
| void Update_Position_Based_State_Parallel(); |
| void Update_Position_Based_State_Serial(); |
| void Delete_Position_Based_State(); |
| void Add_Velocity_Independent_Forces (ARRAY<VECTOR_3D<T> >& F) const; |
| static void Add_Velocity_Independent_Forces_Helper (long thread_id, void* helper_raw); |
| void Add_Velocity_Independent_Forces_Parallel (ARRAY<VECTOR_3D<T> >& F) const; |
| void Add_Velocity_Independent_Forces_Serial (ARRAY<VECTOR_3D<T> >& F) const; |
| void Add_Velocity_Dependent_Forces (ARRAY<VECTOR_3D<T> >& F) const; |
| void Add_Force_Differential (const ARRAY<VECTOR_3D<T> >& dX, ARRAY<VECTOR_3D<T> >& dF) const; |
| static void Add_Force_Differential_Helper (long thread_id, void* helper_raw); |
| void Add_Force_Differential_Parallel (const ARRAY<VECTOR_3D<T> >& dX, ARRAY<VECTOR_3D<T> >& dF) const; |
| void Add_Force_Differential (const ARRAY<VECTOR_3D<T> >& dX_full, ARRAY<VECTOR_3D<T> >& dF_full, const int partition_id) const; |
| void Add_Force_Differential_Serial (const ARRAY<VECTOR_3D<T> >& dX, ARRAY<VECTOR_3D<T> >& dF) const; |
| void Initialize_CFL(); |
| T CFL_Strain_Rate() const; |
| void Semi_Implicit_Impulse_Precomputation (const T time, const T cfl, const T max_dt, ARRAY<T>* time_plus_dt, const bool verbose); |
| void Semi_Implicit_Recompute_Dt (const int element, T& time_plus_dt); |
| void Add_Semi_Implicit_Impulse (const int element, const T dt, T* time_plus_dt); |
| //##################################################################### |
| }; |
| } |
| #endif |