src/parsec/disk-image/parsec/parsec-benchmark/pkgs/apps/facesim/src/Public_Library/Deformable_Objects/DEFORMABLE_OBJECT_3D.h - public/gem5-resources - Git at Google

 //#####################################################################
 // Copyright 2004, Zhaosheng Bao, Ronald Fedkiw, Geoffrey Irving, Sergey Koltakov, Igor Neverov, 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 DEFORMABLE_OBJECT_3D
 //#####################################################################
 #ifndef __DEFORMABLE_OBJECT_3D__
 #define __DEFORMABLE_OBJECT_3D__

 #include "DEFORMABLE_OBJECT.h"
 #include "DEFORMABLE_OBJECT_COLLISIONS_3D.h"
 #include "../Fracture/EMBEDDED_TETRAHEDRALIZED_VOLUME_BOUNDARY_SURFACE.h"
 namespace PhysBAM
 {

 template<class T> class BODY_FORCES_3D;
 template<class T> class DIAGONALIZED_FINITE_VOLUME_3D;
 template<class T> class STRAIN_MEASURE_3D;
 template<class T> class CONSTITUTIVE_MODEL_3D;
 template<class T> class DIAGONALIZED_CONSTITUTIVE_MODEL_3D;
 template<class T> class COLLISION_PENALTY_FORCES;

 template<class T>
 class DEFORMABLE_OBJECT_3D: public DEFORMABLE_OBJECT<T, VECTOR_3D<T> >
 {
 public:
 	using DEFORMABLE_OBJECT<T, VECTOR_3D<T> >::particles;
 	using DEFORMABLE_OBJECT<T, VECTOR_3D<T> >::solids_forces;
 	using DEFORMABLE_OBJECT<T, VECTOR_3D<T> >::external_forces_and_velocities;
 	using DEFORMABLE_OBJECT<T, VECTOR_3D<T> >::id_number;
 	using DEFORMABLE_OBJECT<T, VECTOR_3D<T> >::print_diagnostics;
 	using DEFORMABLE_OBJECT<T, VECTOR_3D<T> >::print_residuals;
 	using DEFORMABLE_OBJECT<T, VECTOR_3D<T> >::Enforce_Definiteness;
 	using DEFORMABLE_OBJECT<T, VECTOR_3D<T> >::Save_Velocity;
 	using DEFORMABLE_OBJECT<T, VECTOR_3D<T> >::Restore_Velocity;
 	using DEFORMABLE_OBJECT<T, VECTOR_3D<T> >::Update_Position_Based_State;
 	using DEFORMABLE_OBJECT<T, VECTOR_3D<T> >::cfl_number;
 	using DEFORMABLE_OBJECT<T, VECTOR_3D<T> >::Add_Velocity_Independent_Forces;
 	using DEFORMABLE_OBJECT<T, VECTOR_3D<T> >::dX_full;
 	using DEFORMABLE_OBJECT<T, VECTOR_3D<T> >::R_full;


 	SEGMENTED_CURVE_3D<T>* segmented_curve;
 	TRIANGULATED_SURFACE<T>* triangulated_surface;
 	TETRAHEDRALIZED_VOLUME<T>* tetrahedralized_volume;
 	EMBEDDED_TETRAHEDRALIZED_VOLUME<T>* embedded_tetrahedralized_volume;
 	EMBEDDED_TETRAHEDRALIZED_VOLUME_BOUNDARY_SURFACE<T>* embedded_tetrahedralized_volume_boundary_surface;
 	DEFORMABLE_OBJECT_COLLISIONS_3D<T> collisions;
 	ARRAY<BODY_FORCES_3D<T>*> body_forces;
 	ARRAY<DIAGONALIZED_FINITE_VOLUME_3D<T>*> diagonalized_finite_volume_3d;
 	ARRAY<STRAIN_MEASURE_3D<T>*> strain_measure_3d;
 	ARRAY<COLLISION_PENALTY_FORCES<T>*> collision_penalty_forces;
 private:
 	ARRAY<CONSTITUTIVE_MODEL_3D<T>*> constitutive_model_3d;
 	ARRAY<DIAGONALIZED_CONSTITUTIVE_MODEL_3D<T>*> diagonalized_constitutive_model_3d;
 public:
 	bool use_nonembedded_self_collision;

 	DEFORMABLE_OBJECT_3D()
 		: segmented_curve (0), triangulated_surface (0), tetrahedralized_volume (0), embedded_tetrahedralized_volume (0),
 		  embedded_tetrahedralized_volume_boundary_surface (0), use_nonembedded_self_collision (false)
 	{}

 	virtual ~DEFORMABLE_OBJECT_3D()
 	{
 		if (segmented_curve)
 		{
 			delete &segmented_curve->segment_mesh;
 			delete segmented_curve;
 		}

 		if (triangulated_surface)
 		{
 			delete &triangulated_surface->triangle_mesh;
 			delete triangulated_surface;
 		}

 		if (tetrahedralized_volume)
 		{
 			delete &tetrahedralized_volume->tetrahedron_mesh;
 			delete tetrahedralized_volume;
 		}

 		if (embedded_tetrahedralized_volume)
 		{
 			delete &embedded_tetrahedralized_volume->tetrahedralized_volume.tetrahedron_mesh;
 			delete &embedded_tetrahedralized_volume->tetrahedralized_volume;
 			delete embedded_tetrahedralized_volume;
 		}

 		delete embedded_tetrahedralized_volume_boundary_surface;
 		Delete_Forces();
 	}

 	void Allocate_Segmented_Curve()
 	{
 		assert (!segmented_curve);
 		segmented_curve = new SEGMENTED_CURVE_3D<T> (* (new SEGMENT_MESH), particles);
 		particles.Update_Velocity();
 		particles.Store_Mass();
 	}

 	void Allocate_Triangulated_Surface()
 	{
 		assert (!triangulated_surface);
 		triangulated_surface = new TRIANGULATED_SURFACE<T> (* (new TRIANGLE_MESH), particles);
 		particles.Update_Velocity();
 		particles.Store_Mass();
 	}

 	void Allocate_Tetrahedralized_Volume()
 	{
 		assert (!tetrahedralized_volume);
 		tetrahedralized_volume = new TETRAHEDRALIZED_VOLUME<T> (* (new TETRAHEDRON_MESH), particles);
 		particles.Update_Velocity();
 		particles.Store_Mass();
 	}

 	void Allocate_Embedded_Triangulated_Surface (const int hashtable_multiplier = 15, const T interpolation_fraction_threshold = 1e-4)
 	{}

 	void Allocate_Embedded_Tetrahedralized_Volume (const int hashtable_multiplier = 15, const T interpolation_fraction_threshold = 1e-4)
 	{
 		assert (!embedded_tetrahedralized_volume);
 		assert (!embedded_tetrahedralized_volume_boundary_surface);
 		embedded_tetrahedralized_volume = new EMBEDDED_TETRAHEDRALIZED_VOLUME<T> (* (new TETRAHEDRALIZED_VOLUME<T> (* (new TETRAHEDRON_MESH), particles)));
 		particles.Update_Velocity();
 		particles.Store_Mass();
 		embedded_tetrahedralized_volume->Initialize_Parents_To_Embedded_Particles_Hash_Table (hashtable_multiplier);
 		embedded_tetrahedralized_volume->Set_Interpolation_Fraction_Threshold (interpolation_fraction_threshold);
 		embedded_tetrahedralized_volume_boundary_surface = new EMBEDDED_TETRAHEDRALIZED_VOLUME_BOUNDARY_SURFACE<T> (*embedded_tetrahedralized_volume);
 	}

 	virtual void Initialize_Collision_Geometry()
 	{
 		if (segmented_curve) collisions.Initialize_Segmented_Curve (*segmented_curve);
 		else if (triangulated_surface) collisions.Initialize_Triangulated_Surface (*triangulated_surface);
 		else if (tetrahedralized_volume)
 		{
 			if (!tetrahedralized_volume->triangulated_surface) tetrahedralized_volume->Initialize_Triangulated_Surface();

 			collisions.Initialize_Triangulated_Surface (*tetrahedralized_volume->triangulated_surface);
 		}
 		else if (embedded_tetrahedralized_volume) collisions.Initialize_Triangulated_Surface (embedded_tetrahedralized_volume_boundary_surface->boundary_surface);
 	}

 	void Initialize_Object_Collisions (const bool collide_with_interior = false, const bool enforce_tangential_collision_velocity = false)
 	{
 		if (segmented_curve) collisions.Initialize_Object_Collisions (segmented_curve->particles, enforce_tangential_collision_velocity);
 		else if (triangulated_surface) collisions.Initialize_Object_Collisions (triangulated_surface->particles, enforce_tangential_collision_velocity);
 		else if (tetrahedralized_volume)
 		{
 			collisions.Initialize_Object_Collisions (tetrahedralized_volume->particles, enforce_tangential_collision_velocity);

 			if (!collide_with_interior)
 			{
 				TETRAHEDRON_MESH& tetrahedron_mesh = tetrahedralized_volume->tetrahedron_mesh;
 				LIST_ARRAY<bool>* old_node_on_boundary = tetrahedron_mesh.node_on_boundary;
 				tetrahedron_mesh.node_on_boundary = 0;
 				tetrahedron_mesh.Initialize_Node_On_Boundary();
 				ARRAY<bool>::exchange_arrays (collisions.check_collision, tetrahedron_mesh.node_on_boundary->array);
 				delete tetrahedron_mesh.node_on_boundary;
 				tetrahedron_mesh.node_on_boundary = old_node_on_boundary;
 			}
 		}
 		else if (embedded_tetrahedralized_volume)
 		{
 			collisions.Initialize_Object_Collisions (embedded_tetrahedralized_volume->tetrahedralized_volume.particles, enforce_tangential_collision_velocity);

 			if (collide_with_interior) assert (false);
 		}
 	}

 	void Allocate_Asynchronous_Time_Stepping()
 	{}

 	// Should be accelerated by using the triangle_hierarchy to eliminate tests on parts of the surface that aren't close to collision bodies.
 	int Adjust_Nodes_For_Collision_Body_Collisions (const T dt)
 	{
 		if (embedded_tetrahedralized_volume) return collisions.Adjust_Nodes_For_Collision_Body_Collisions (*embedded_tetrahedralized_volume, *embedded_tetrahedralized_volume_boundary_surface, dt);
 		else return collisions.Adjust_Nodes_For_Collision_Body_Collisions (dt);
 	}

 	void Restore_Self_Collision_Free_State()
 	{
 		collisions.Restore_Self_Collision_Free_State();

 		if (embedded_tetrahedralized_volume) // triangle_collisions only restores the boundary mesh, need to propagate this back to the tetrahedralized volume
 		{
 			for (int p = 1; p <= embedded_tetrahedralized_volume->tetrahedralized_volume.particles.number; p++)
 			{
 				embedded_tetrahedralized_volume->tetrahedralized_volume.particles.X (p) =
 					embedded_tetrahedralized_volume_boundary_surface->boundary_particles.X (p + embedded_tetrahedralized_volume->embedded_particles.number);
 				embedded_tetrahedralized_volume->tetrahedralized_volume.particles.V (p) =
 					embedded_tetrahedralized_volume_boundary_surface->boundary_particles.V (p + embedded_tetrahedralized_volume->embedded_particles.number);
 			}

 			for (int p = 1; p <= embedded_tetrahedralized_volume->tetrahedralized_volume.particles.number; p++) // THIS IS NEW!
 			{
 				embedded_tetrahedralized_volume->embedded_particles.X (p) = embedded_tetrahedralized_volume_boundary_surface->boundary_particles.X (p);
 				embedded_tetrahedralized_volume->embedded_particles.V (p) = embedded_tetrahedralized_volume_boundary_surface->boundary_particles.V (p);
 			}
 		}
 	}

 	void Disable_Finite_Volume_Damping()
 	{
 		for (int k = 1; k <= diagonalized_finite_volume_3d.m; k++) diagonalized_finite_volume_3d (k)->Use_Velocity_Dependent_Forces (false);
 	}

 	void Disable_Spring_Elasticity()
 	{}

 	template<class RW>
 	void Read_Static_Variables (const std::string& prefix, const int frame = -1)
 	{
 		char version;
 		bool has_segmented_curve, has_triangulated_surface, has_tetrahedralized_volume, has_embedded_tetrahedralized_volume, has_embedded_triangulated_surface = false;
 		std::string f = frame == -1 ? "" : STRING_UTILITIES::string_sprintf (".%d", frame);
 		std::istream* mesh_type_input = FILE_UTILITIES::Safe_Open_Input (prefix + "mesh_type" + f);
 		Read_Binary<RW> (*mesh_type_input, version, has_segmented_curve, has_triangulated_surface, has_tetrahedralized_volume, has_embedded_tetrahedralized_volume);

 		if (version == 2) Read_Binary<RW> (*mesh_type_input, has_embedded_triangulated_surface);
 		else if (version != 1)
 		{
 			std::cerr << "Unrecognized deformable object version number " << (int) version << std::endl;
 			exit (1);
 		}

 		delete mesh_type_input;

 		if (has_segmented_curve)
 		{
 			Allocate_Segmented_Curve();
 			FILE_UTILITIES::Read_From_File<RW> (prefix + "segment_mesh" + f, segmented_curve->segment_mesh);
 		}

 		if (has_triangulated_surface)
 		{
 			Allocate_Triangulated_Surface();
 			FILE_UTILITIES::Read_From_File<RW> (prefix + "triangle_mesh" + f, triangulated_surface->triangle_mesh);
 		}

 		if (has_tetrahedralized_volume)
 		{
 			Allocate_Tetrahedralized_Volume();
 			FILE_UTILITIES::Read_From_File<RW> (prefix + "tetrahedron_mesh" + f, tetrahedralized_volume->tetrahedron_mesh);
 		}

 		if (has_embedded_tetrahedralized_volume)
 		{
 			Allocate_Embedded_Tetrahedralized_Volume();
 			FILE_UTILITIES::Read_From_File<RW> (prefix + "embedded_tetrahedralized_volume" + f, *embedded_tetrahedralized_volume);
 		}
 	}

 	template<class RW>
 	void Read_Dynamic_Variables (const std::string& prefix, const int frame)
 	{
 		particles.template Read<RW> (prefix, frame);

 		if (embedded_tetrahedralized_volume) embedded_tetrahedralized_volume->embedded_particles.template Read<RW> (prefix + "embedded_surface_", frame);
 	}

 	template<class RW>
 	void Write_Static_Variables (const std::string& prefix, const int frame = -1) const
 	{
 		char version = 2;
 		std::string f = frame == -1 ? "" : STRING_UTILITIES::string_sprintf (".%d", frame);
 		FILE_UTILITIES::Write_To_File<RW> (prefix + "mesh_type" + f, version, segmented_curve != 0, triangulated_surface != 0, tetrahedralized_volume != 0, embedded_tetrahedralized_volume != 0, false);

 		if (segmented_curve) FILE_UTILITIES::Write_To_File<RW> (prefix + "segment_mesh" + f, segmented_curve->segment_mesh);

 		if (triangulated_surface) FILE_UTILITIES::Write_To_File<RW> (prefix + "triangle_mesh" + f, triangulated_surface->triangle_mesh);

 		if (tetrahedralized_volume) FILE_UTILITIES::Write_To_File<RW> (prefix + "tetrahedron_mesh" + f, tetrahedralized_volume->tetrahedron_mesh);

 		if (embedded_tetrahedralized_volume) // writes positions
 		{
 			FILE_UTILITIES::Write_To_File<RW> (prefix + "embedded_tetrahedralized_volume" + f, *embedded_tetrahedralized_volume);
 			FILE_UTILITIES::Write_To_File<RW> (prefix + "embedded_tetrahedralized_volume_boundary_surface" + f, embedded_tetrahedralized_volume_boundary_surface->boundary_surface);
 		}
 	}

 	template<class RW>
 	void Write_Dynamic_Variables (const std::string& prefix, const int frame) const
 	{
 		particles.template Write<RW> (prefix, frame);

 		if (embedded_tetrahedralized_volume) embedded_tetrahedralized_volume->embedded_particles.template Write<RW> (prefix + "embedded_surface_", frame);
 	}

 //#####################################################################
 	void Delete_Forces();
 	void Update_Collision_Penalty_Forces_And_Derivatives();
 	void Advance_One_Time_Step (const T time, const T dt, const T cg_tolerance, const int cg_iterations, const bool perform_collision_body_collisions, const bool verbose);
 	void Advance_One_Time_Step_Semi_Implicit (const T time, const T dt, const bool perform_collision_body_collisions, const bool verbose);
 	void Advance_One_Time_Step_Quasistatic (const T time, const T dt, const T cg_tolerance, const int cg_iterations, const T newton_tolerance, const int newton_iterations, const bool use_partially_converged_result,
 						const bool verbose);
 public:
 	int Add_Quasistatic_Neo_Hookean_Elasticity (TETRAHEDRALIZED_VOLUME<T>& tetrahedralized_volume, const T youngs_modulus = 3e6, const T poissons_ratio = .475, const bool verbose = true);
 	int Add_Quasistatic_Diagonalized_Linear_Finite_Volume (TETRAHEDRALIZED_VOLUME<T>& tetrahedralized_volume, const T youngs_modulus = 3e6, const T poissons_ratio = .475, const bool verbose = true,
 			const bool precompute_stiffness_matrix = false);
 	int Add_Quasistatic_Diagonalized_Neo_Hookean_Elasticity (TETRAHEDRALIZED_VOLUME<T>& tetrahedralized_volume, const T youngs_modulus = 3e6, const T poissons_ratio = .475, const T failure_threshold = .25,
 			const bool verbose = true, const bool precompute_stiffness_matrix = false);
 	int Add_Quasistatic_Diagonalized_Linear_Elasticity (TETRAHEDRALIZED_VOLUME<T>& tetrahedralized_volume, const T youngs_modulus = 3e6, const T poissons_ratio = .475, const bool verbose = true);
 	int Add_Diagonalized_Mooney_Rivlin (TETRAHEDRALIZED_VOLUME<T>& tetrahedralized_volume, const T mu_10 = (T) 6e4, const T mu_01 = (T) 2e4, const T kappa = (T) 6e4, const T Rayleigh_coefficient = .05,
 					    const T failure_threshold = (T).25, const bool limit_time_step_by_strain_rate = true, const T max_strain_per_time_step = .1, const bool use_rest_state_for_strain_rate = true, const bool verbose = true);
 	int Add_Quasistatic_Diagonalized_Mooney_Rivlin (TETRAHEDRALIZED_VOLUME<T>& tetrahedralized_volume, const T mu_10 = (T) 6e4, const T mu_01 = (T) 2e4, const T kappa = (T) 6e4, const T failure_threshold = (T).25,
 			const bool verbose = true, const bool precompute_stiffness_matrix = false);
 	int Add_Diagonalized_Face (TETRAHEDRALIZED_VOLUME<T>& tetrahedralized_volume, const LIST_ARRAY<LIST_ARRAY<int> >& muscle_tets, const LIST_ARRAY<LIST_ARRAY<VECTOR_3D<T> > >& muscle_fibers,
 				   const LIST_ARRAY<LIST_ARRAY<T> >& muscle_densities, LIST_ARRAY<T>& muscle_activations, const LIST_ARRAY<T>* peak_isometric_stress = 0, const T mu_10 = (T) 6e4, const T mu_01 = (T) 2e4, const T kappa = (T) 6e4,
 				   const T Rayleigh_coefficient = (T).05, const T failure_threshold = (T).25, const bool limit_time_step_by_strain_rate = true, const T max_strain_per_time_step = (T).1,
 				   const bool use_rest_state_for_strain_rate = true, const bool verbose = true);
 	int Add_Quasistatic_Diagonalized_Face (TETRAHEDRALIZED_VOLUME<T>& tetrahedralized_volume, const LIST_ARRAY<LIST_ARRAY<int> >& muscle_tets, const LIST_ARRAY<LIST_ARRAY<VECTOR_3D<T> > >& muscle_fibers,
 					       const LIST_ARRAY<LIST_ARRAY<T> >& muscle_densities, LIST_ARRAY<T>& muscle_activations, int* single_activation_used_for_force_derivative = 0, const LIST_ARRAY<T>* peak_isometric_stress = 0,
 					       DIAGONALIZED_CONSTITUTIVE_MODEL_3D<T> **face_constitutive_model = 0, const T mu_10 = (T) 6e4, const T mu_01 = (T) 2e4, const T kappa = (T) 6e4, const T failure_threshold = (T).25, const bool verbose = false);
 //#####################################################################
 };
 }
 #endif
	//#####################################################################
	// Copyright 2004, Zhaosheng Bao, Ronald Fedkiw, Geoffrey Irving, Sergey Koltakov, Igor Neverov, 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 DEFORMABLE_OBJECT_3D
	//#####################################################################
	#ifndef __DEFORMABLE_OBJECT_3D__
	#define __DEFORMABLE_OBJECT_3D__

	#include "DEFORMABLE_OBJECT.h"
	#include "DEFORMABLE_OBJECT_COLLISIONS_3D.h"
	#include "../Fracture/EMBEDDED_TETRAHEDRALIZED_VOLUME_BOUNDARY_SURFACE.h"
	namespace PhysBAM
	{

	template<class T> class BODY_FORCES_3D;
	template<class T> class DIAGONALIZED_FINITE_VOLUME_3D;
	template<class T> class STRAIN_MEASURE_3D;
	template<class T> class CONSTITUTIVE_MODEL_3D;
	template<class T> class DIAGONALIZED_CONSTITUTIVE_MODEL_3D;
	template<class T> class COLLISION_PENALTY_FORCES;

	template<class T>
	class DEFORMABLE_OBJECT_3D: public DEFORMABLE_OBJECT<T, VECTOR_3D<T> >
	{
	public:
	using DEFORMABLE_OBJECT<T, VECTOR_3D<T> >::particles;
	using DEFORMABLE_OBJECT<T, VECTOR_3D<T> >::solids_forces;
	using DEFORMABLE_OBJECT<T, VECTOR_3D<T> >::external_forces_and_velocities;
	using DEFORMABLE_OBJECT<T, VECTOR_3D<T> >::id_number;
	using DEFORMABLE_OBJECT<T, VECTOR_3D<T> >::print_diagnostics;
	using DEFORMABLE_OBJECT<T, VECTOR_3D<T> >::print_residuals;
	using DEFORMABLE_OBJECT<T, VECTOR_3D<T> >::Enforce_Definiteness;
	using DEFORMABLE_OBJECT<T, VECTOR_3D<T> >::Save_Velocity;
	using DEFORMABLE_OBJECT<T, VECTOR_3D<T> >::Restore_Velocity;
	using DEFORMABLE_OBJECT<T, VECTOR_3D<T> >::Update_Position_Based_State;
	using DEFORMABLE_OBJECT<T, VECTOR_3D<T> >::cfl_number;
	using DEFORMABLE_OBJECT<T, VECTOR_3D<T> >::Add_Velocity_Independent_Forces;
	using DEFORMABLE_OBJECT<T, VECTOR_3D<T> >::dX_full;
	using DEFORMABLE_OBJECT<T, VECTOR_3D<T> >::R_full;


	SEGMENTED_CURVE_3D<T>* segmented_curve;
	TRIANGULATED_SURFACE<T>* triangulated_surface;
	TETRAHEDRALIZED_VOLUME<T>* tetrahedralized_volume;
	EMBEDDED_TETRAHEDRALIZED_VOLUME<T>* embedded_tetrahedralized_volume;
	EMBEDDED_TETRAHEDRALIZED_VOLUME_BOUNDARY_SURFACE<T>* embedded_tetrahedralized_volume_boundary_surface;
	DEFORMABLE_OBJECT_COLLISIONS_3D<T> collisions;
	ARRAY<BODY_FORCES_3D<T>*> body_forces;
	ARRAY<DIAGONALIZED_FINITE_VOLUME_3D<T>*> diagonalized_finite_volume_3d;
	ARRAY<STRAIN_MEASURE_3D<T>*> strain_measure_3d;
	ARRAY<COLLISION_PENALTY_FORCES<T>*> collision_penalty_forces;
	private:
	ARRAY<CONSTITUTIVE_MODEL_3D<T>*> constitutive_model_3d;
	ARRAY<DIAGONALIZED_CONSTITUTIVE_MODEL_3D<T>*> diagonalized_constitutive_model_3d;
	public:
	bool use_nonembedded_self_collision;

	DEFORMABLE_OBJECT_3D()
	: segmented_curve (0), triangulated_surface (0), tetrahedralized_volume (0), embedded_tetrahedralized_volume (0),
	embedded_tetrahedralized_volume_boundary_surface (0), use_nonembedded_self_collision (false)
	{}

	virtual ~DEFORMABLE_OBJECT_3D()
	{
	if (segmented_curve)
	{
	delete &segmented_curve->segment_mesh;
	delete segmented_curve;
	}

	if (triangulated_surface)
	{
	delete &triangulated_surface->triangle_mesh;
	delete triangulated_surface;
	}

	if (tetrahedralized_volume)
	{
	delete &tetrahedralized_volume->tetrahedron_mesh;
	delete tetrahedralized_volume;
	}

	if (embedded_tetrahedralized_volume)
	{
	delete &embedded_tetrahedralized_volume->tetrahedralized_volume.tetrahedron_mesh;
	delete &embedded_tetrahedralized_volume->tetrahedralized_volume;
	delete embedded_tetrahedralized_volume;
	}

	delete embedded_tetrahedralized_volume_boundary_surface;
	Delete_Forces();
	}

	void Allocate_Segmented_Curve()
	{
	assert (!segmented_curve);
	segmented_curve = new SEGMENTED_CURVE_3D<T> (* (new SEGMENT_MESH), particles);
	particles.Update_Velocity();
	particles.Store_Mass();
	}

	void Allocate_Triangulated_Surface()
	{
	assert (!triangulated_surface);
	triangulated_surface = new TRIANGULATED_SURFACE<T> (* (new TRIANGLE_MESH), particles);
	particles.Update_Velocity();
	particles.Store_Mass();
	}

	void Allocate_Tetrahedralized_Volume()
	{
	assert (!tetrahedralized_volume);
	tetrahedralized_volume = new TETRAHEDRALIZED_VOLUME<T> (* (new TETRAHEDRON_MESH), particles);
	particles.Update_Velocity();
	particles.Store_Mass();
	}

	void Allocate_Embedded_Triangulated_Surface (const int hashtable_multiplier = 15, const T interpolation_fraction_threshold = 1e-4)
	{}

	void Allocate_Embedded_Tetrahedralized_Volume (const int hashtable_multiplier = 15, const T interpolation_fraction_threshold = 1e-4)
	{
	assert (!embedded_tetrahedralized_volume);
	assert (!embedded_tetrahedralized_volume_boundary_surface);
	embedded_tetrahedralized_volume = new EMBEDDED_TETRAHEDRALIZED_VOLUME<T> (* (new TETRAHEDRALIZED_VOLUME<T> (* (new TETRAHEDRON_MESH), particles)));
	particles.Update_Velocity();
	particles.Store_Mass();
	embedded_tetrahedralized_volume->Initialize_Parents_To_Embedded_Particles_Hash_Table (hashtable_multiplier);
	embedded_tetrahedralized_volume->Set_Interpolation_Fraction_Threshold (interpolation_fraction_threshold);
	embedded_tetrahedralized_volume_boundary_surface = new EMBEDDED_TETRAHEDRALIZED_VOLUME_BOUNDARY_SURFACE<T> (*embedded_tetrahedralized_volume);
	}

	virtual void Initialize_Collision_Geometry()
	{
	if (segmented_curve) collisions.Initialize_Segmented_Curve (*segmented_curve);
	else if (triangulated_surface) collisions.Initialize_Triangulated_Surface (*triangulated_surface);
	else if (tetrahedralized_volume)
	{
	if (!tetrahedralized_volume->triangulated_surface) tetrahedralized_volume->Initialize_Triangulated_Surface();

	collisions.Initialize_Triangulated_Surface (*tetrahedralized_volume->triangulated_surface);
	}
	else if (embedded_tetrahedralized_volume) collisions.Initialize_Triangulated_Surface (embedded_tetrahedralized_volume_boundary_surface->boundary_surface);
	}

	void Initialize_Object_Collisions (const bool collide_with_interior = false, const bool enforce_tangential_collision_velocity = false)
	{
	if (segmented_curve) collisions.Initialize_Object_Collisions (segmented_curve->particles, enforce_tangential_collision_velocity);
	else if (triangulated_surface) collisions.Initialize_Object_Collisions (triangulated_surface->particles, enforce_tangential_collision_velocity);
	else if (tetrahedralized_volume)
	{
	collisions.Initialize_Object_Collisions (tetrahedralized_volume->particles, enforce_tangential_collision_velocity);

	if (!collide_with_interior)
	{
	TETRAHEDRON_MESH& tetrahedron_mesh = tetrahedralized_volume->tetrahedron_mesh;
	LIST_ARRAY<bool>* old_node_on_boundary = tetrahedron_mesh.node_on_boundary;
	tetrahedron_mesh.node_on_boundary = 0;
	tetrahedron_mesh.Initialize_Node_On_Boundary();
	ARRAY<bool>::exchange_arrays (collisions.check_collision, tetrahedron_mesh.node_on_boundary->array);
	delete tetrahedron_mesh.node_on_boundary;
	tetrahedron_mesh.node_on_boundary = old_node_on_boundary;
	}
	}
	else if (embedded_tetrahedralized_volume)
	{
	collisions.Initialize_Object_Collisions (embedded_tetrahedralized_volume->tetrahedralized_volume.particles, enforce_tangential_collision_velocity);

	if (collide_with_interior) assert (false);
	}
	}

	void Allocate_Asynchronous_Time_Stepping()
	{}

	// Should be accelerated by using the triangle_hierarchy to eliminate tests on parts of the surface that aren't close to collision bodies.
	int Adjust_Nodes_For_Collision_Body_Collisions (const T dt)
	{
	if (embedded_tetrahedralized_volume) return collisions.Adjust_Nodes_For_Collision_Body_Collisions (embedded_tetrahedralized_volume, embedded_tetrahedralized_volume_boundary_surface, dt);
	else return collisions.Adjust_Nodes_For_Collision_Body_Collisions (dt);
	}

	void Restore_Self_Collision_Free_State()
	{
	collisions.Restore_Self_Collision_Free_State();

	if (embedded_tetrahedralized_volume) // triangle_collisions only restores the boundary mesh, need to propagate this back to the tetrahedralized volume
	{
	for (int p = 1; p <= embedded_tetrahedralized_volume->tetrahedralized_volume.particles.number; p++)
	{
	embedded_tetrahedralized_volume->tetrahedralized_volume.particles.X (p) =
	embedded_tetrahedralized_volume_boundary_surface->boundary_particles.X (p + embedded_tetrahedralized_volume->embedded_particles.number);
	embedded_tetrahedralized_volume->tetrahedralized_volume.particles.V (p) =
	embedded_tetrahedralized_volume_boundary_surface->boundary_particles.V (p + embedded_tetrahedralized_volume->embedded_particles.number);
	}

	for (int p = 1; p <= embedded_tetrahedralized_volume->tetrahedralized_volume.particles.number; p++) // THIS IS NEW!
	{
	embedded_tetrahedralized_volume->embedded_particles.X (p) = embedded_tetrahedralized_volume_boundary_surface->boundary_particles.X (p);
	embedded_tetrahedralized_volume->embedded_particles.V (p) = embedded_tetrahedralized_volume_boundary_surface->boundary_particles.V (p);
	}
	}
	}

	void Disable_Finite_Volume_Damping()
	{
	for (int k = 1; k <= diagonalized_finite_volume_3d.m; k++) diagonalized_finite_volume_3d (k)->Use_Velocity_Dependent_Forces (false);
	}

	void Disable_Spring_Elasticity()
	{}

	template<class RW>
	void Read_Static_Variables (const std::string& prefix, const int frame = -1)
	{
	char version;
	bool has_segmented_curve, has_triangulated_surface, has_tetrahedralized_volume, has_embedded_tetrahedralized_volume, has_embedded_triangulated_surface = false;
	std::string f = frame == -1 ? "" : STRING_UTILITIES::string_sprintf (".%d", frame);
	std::istream* mesh_type_input = FILE_UTILITIES::Safe_Open_Input (prefix + "mesh_type" + f);
	Read_Binary<RW> (*mesh_type_input, version, has_segmented_curve, has_triangulated_surface, has_tetrahedralized_volume, has_embedded_tetrahedralized_volume);

	if (version == 2) Read_Binary<RW> (*mesh_type_input, has_embedded_triangulated_surface);
	else if (version != 1)
	{
	std::cerr << "Unrecognized deformable object version number " << (int) version << std::endl;
	exit (1);
	}

	delete mesh_type_input;

	if (has_segmented_curve)
	{
	Allocate_Segmented_Curve();
	FILE_UTILITIES::Read_From_File<RW> (prefix + "segment_mesh" + f, segmented_curve->segment_mesh);
	}

	if (has_triangulated_surface)
	{
	Allocate_Triangulated_Surface();
	FILE_UTILITIES::Read_From_File<RW> (prefix + "triangle_mesh" + f, triangulated_surface->triangle_mesh);
	}

	if (has_tetrahedralized_volume)
	{
	Allocate_Tetrahedralized_Volume();
	FILE_UTILITIES::Read_From_File<RW> (prefix + "tetrahedron_mesh" + f, tetrahedralized_volume->tetrahedron_mesh);
	}

	if (has_embedded_tetrahedralized_volume)
	{
	Allocate_Embedded_Tetrahedralized_Volume();
	FILE_UTILITIES::Read_From_File<RW> (prefix + "embedded_tetrahedralized_volume" + f, *embedded_tetrahedralized_volume);
	}
	}

	template<class RW>
	void Read_Dynamic_Variables (const std::string& prefix, const int frame)
	{
	particles.template Read<RW> (prefix, frame);

	if (embedded_tetrahedralized_volume) embedded_tetrahedralized_volume->embedded_particles.template Read<RW> (prefix + "embedded_surface_", frame);
	}

	template<class RW>
	void Write_Static_Variables (const std::string& prefix, const int frame = -1) const
	{
	char version = 2;
	std::string f = frame == -1 ? "" : STRING_UTILITIES::string_sprintf (".%d", frame);
	FILE_UTILITIES::Write_To_File<RW> (prefix + "mesh_type" + f, version, segmented_curve != 0, triangulated_surface != 0, tetrahedralized_volume != 0, embedded_tetrahedralized_volume != 0, false);

	if (segmented_curve) FILE_UTILITIES::Write_To_File<RW> (prefix + "segment_mesh" + f, segmented_curve->segment_mesh);

	if (triangulated_surface) FILE_UTILITIES::Write_To_File<RW> (prefix + "triangle_mesh" + f, triangulated_surface->triangle_mesh);

	if (tetrahedralized_volume) FILE_UTILITIES::Write_To_File<RW> (prefix + "tetrahedron_mesh" + f, tetrahedralized_volume->tetrahedron_mesh);

	if (embedded_tetrahedralized_volume) // writes positions
	{
	FILE_UTILITIES::Write_To_File<RW> (prefix + "embedded_tetrahedralized_volume" + f, *embedded_tetrahedralized_volume);
	FILE_UTILITIES::Write_To_File<RW> (prefix + "embedded_tetrahedralized_volume_boundary_surface" + f, embedded_tetrahedralized_volume_boundary_surface->boundary_surface);
	}
	}

	template<class RW>
	void Write_Dynamic_Variables (const std::string& prefix, const int frame) const
	{
	particles.template Write<RW> (prefix, frame);

	if (embedded_tetrahedralized_volume) embedded_tetrahedralized_volume->embedded_particles.template Write<RW> (prefix + "embedded_surface_", frame);
	}

	//#####################################################################
	void Delete_Forces();
	void Update_Collision_Penalty_Forces_And_Derivatives();
	void Advance_One_Time_Step (const T time, const T dt, const T cg_tolerance, const int cg_iterations, const bool perform_collision_body_collisions, const bool verbose);
	void Advance_One_Time_Step_Semi_Implicit (const T time, const T dt, const bool perform_collision_body_collisions, const bool verbose);
	void Advance_One_Time_Step_Quasistatic (const T time, const T dt, const T cg_tolerance, const int cg_iterations, const T newton_tolerance, const int newton_iterations, const bool use_partially_converged_result,
	const bool verbose);
	public:
	int Add_Quasistatic_Neo_Hookean_Elasticity (TETRAHEDRALIZED_VOLUME<T>& tetrahedralized_volume, const T youngs_modulus = 3e6, const T poissons_ratio = .475, const bool verbose = true);
	int Add_Quasistatic_Diagonalized_Linear_Finite_Volume (TETRAHEDRALIZED_VOLUME<T>& tetrahedralized_volume, const T youngs_modulus = 3e6, const T poissons_ratio = .475, const bool verbose = true,
	const bool precompute_stiffness_matrix = false);
	int Add_Quasistatic_Diagonalized_Neo_Hookean_Elasticity (TETRAHEDRALIZED_VOLUME<T>& tetrahedralized_volume, const T youngs_modulus = 3e6, const T poissons_ratio = .475, const T failure_threshold = .25,
	const bool verbose = true, const bool precompute_stiffness_matrix = false);
	int Add_Quasistatic_Diagonalized_Linear_Elasticity (TETRAHEDRALIZED_VOLUME<T>& tetrahedralized_volume, const T youngs_modulus = 3e6, const T poissons_ratio = .475, const bool verbose = true);
	int Add_Diagonalized_Mooney_Rivlin (TETRAHEDRALIZED_VOLUME<T>& tetrahedralized_volume, const T mu_10 = (T) 6e4, const T mu_01 = (T) 2e4, const T kappa = (T) 6e4, const T Rayleigh_coefficient = .05,
	const T failure_threshold = (T).25, const bool limit_time_step_by_strain_rate = true, const T max_strain_per_time_step = .1, const bool use_rest_state_for_strain_rate = true, const bool verbose = true);
	int Add_Quasistatic_Diagonalized_Mooney_Rivlin (TETRAHEDRALIZED_VOLUME<T>& tetrahedralized_volume, const T mu_10 = (T) 6e4, const T mu_01 = (T) 2e4, const T kappa = (T) 6e4, const T failure_threshold = (T).25,
	const bool verbose = true, const bool precompute_stiffness_matrix = false);
	int Add_Diagonalized_Face (TETRAHEDRALIZED_VOLUME<T>& tetrahedralized_volume, const LIST_ARRAY<LIST_ARRAY<int> >& muscle_tets, const LIST_ARRAY<LIST_ARRAY<VECTOR_3D<T> > >& muscle_fibers,
	const LIST_ARRAY<LIST_ARRAY<T> >& muscle_densities, LIST_ARRAY<T>& muscle_activations, const LIST_ARRAY<T>* peak_isometric_stress = 0, const T mu_10 = (T) 6e4, const T mu_01 = (T) 2e4, const T kappa = (T) 6e4,
	const T Rayleigh_coefficient = (T).05, const T failure_threshold = (T).25, const bool limit_time_step_by_strain_rate = true, const T max_strain_per_time_step = (T).1,
	const bool use_rest_state_for_strain_rate = true, const bool verbose = true);
	int Add_Quasistatic_Diagonalized_Face (TETRAHEDRALIZED_VOLUME<T>& tetrahedralized_volume, const LIST_ARRAY<LIST_ARRAY<int> >& muscle_tets, const LIST_ARRAY<LIST_ARRAY<VECTOR_3D<T> > >& muscle_fibers,
	const LIST_ARRAY<LIST_ARRAY<T> >& muscle_densities, LIST_ARRAY<T>& muscle_activations, int* single_activation_used_for_force_derivative = 0, const LIST_ARRAY<T>* peak_isometric_stress = 0,
	DIAGONALIZED_CONSTITUTIVE_MODEL_3D<T> **face_constitutive_model = 0, const T mu_10 = (T) 6e4, const T mu_01 = (T) 2e4, const T kappa = (T) 6e4, const T failure_threshold = (T).25, const bool verbose = false);
	//#####################################################################
	};
	}
	#endif