Dynamic heart (#29)

* change pressure BC to be define by body part

* SpringNormalOnSurfaceParticles with body part

* option to use inner or outer surface for spring

* zero position of particles fixed

* adding function to avoid LevelSetComplexShape

* "virtual" added to comply with github CI test

* both options for normalDirectionCalculation

* "virtual" removed since no change in CI test

* additional findNormalDirection dor 2D_build

* comment to see if geometry_level_set is called

* accelaration output added to vtu

* adding displacement output & removing density

* Vec3d to Vecd

* adding normal vectors to vtu-File

* von Mises Strain implemented

* CMake name change to differentiate BUILD_TEST

* test for strain calculation

* correction of von Mises strain calculation

* displ., normal, strain output calc. added in 2D

Co-authored-by: JohnVirtonomy <[email protected]>

CMakeLists.txt

@@ -12,13 +12,13 @@ option(STATIC_BUILD "STATIC_BUILD" 0)
 option(ONLY_3D "ONLY_3D" 0)
 option(BUILD_WITH_SIMBODY "BUILD_WITH_SIMBODY" 0)
 option(EMSCRIPTEN "EMSCRIPTEN" 0)
-option(BUILD_TESTS "BUILD_TESTS" 1)
+option(BUILD_SPHINXSYS_TESTS "BUILD_SPHINXSYS_TESTS" 1)
 
 if(EMSCRIPTEN)
     set(STATIC_BUILD 1)
     set(BUILD_WITH_SIMBODY 1)
     set(ONLY_3D 1)
-    set(BUILD_TESTS 1)
+    set(BUILD_SPHINXSYS_TESTS 1)
     add_definitions(-D__SIMBODY_WITHOUT_LAPACK__)
     add_definitions(-D__EIGEN__)
     add_definitions(-D__EMSCRIPTEN__)
@@ -126,7 +126,7 @@ endif()
 
 add_subdirectory(SPHINXsys)
 
-if(BUILD_TESTS)
+if(BUILD_SPHINXSYS_TESTS)
     add_subdirectory(tests)
 endif()
 

SPHINXsys/src/for_2D_build/geometries/geometry.cpp

@@ -259,6 +259,22 @@ namespace SPH
 		return is_contain ? direction_to_surface : -1.0 * direction_to_surface;
 	}
 	//=================================================================================================//
+	Vec2d ComplexShape::findNormalDirectionComplexShape(const Vec2d &input_pnt) //function to differentiate from LevelSetComplexShape::findNormalDirection
+	{
+		bool is_contain = checkContain(input_pnt);
+		Vecd displacement_to_surface = findClosestPoint(input_pnt) - input_pnt;
+		while (displacement_to_surface.norm() < Eps)
+		{
+			Vecd jittered = input_pnt; //jittering
+			for (int l = 0; l != input_pnt.size(); ++l)
+				jittered[l] = input_pnt[l] + (((Real)rand() / (RAND_MAX)) - 0.5) * 100.0 * Eps;
+			if (checkContain(jittered) == is_contain)
+				displacement_to_surface = findClosestPoint(jittered) - jittered;
+		}
+		Vecd direction_to_surface = displacement_to_surface.normalize();
+		return is_contain ? direction_to_surface : -1.0 * direction_to_surface;
+	}
+	//=================================================================================================//
 	BoundingBox ComplexShape::findBounds()
 	{
 		return multi_ploygen_.findBounds();

SPHINXsys/src/for_2D_build/geometries/geometry.h

@@ -99,6 +99,7 @@ namespace SPH
 		virtual Real findSignedDistance(const Vec2d &input_pnt);
 		/** Normal direction point toward outside of the complex shape. */
 		virtual Vec2d findNormalDirection(const Vec2d &input_pnt);
+		Vec2d findNormalDirectionComplexShape(const Vec2d &input_pnt);
 
 	protected:
 		MultiPolygon multi_ploygen_;

SPHINXsys/src/for_2D_build/particles/solid_particles_supplementary.cpp

@@ -32,4 +32,34 @@ namespace SPH {
 			+ 3.0 * sigmaxy * sigmaxy);
 	}
 	//=================================================================================================//
+	Vecd ElasticSolidParticles::displacement(size_t particle_i) //not tested in 2D
+	{
+		Vecd disp = pos_n_[particle_i] - pos_0_[particle_i];
+		return disp;
+	}
+	//=================================================================================================//
+	Vecd ElasticSolidParticles::normal(size_t particle_i) //not tested in 2D
+	{
+		Vecd normal_vec = n_[particle_i];
+		return normal_vec;
+	}
+	//=================================================================================================//
+	Real ElasticSolidParticles::von_Mises_strain(size_t particle_i) //not tested in 2D
+	{
+
+		Mat2d F = F_[particle_i];
+		Mat2d epsilon = 0.5 * (~F * F - Mat2d(1.0)); //calculation of the Green-Lagrange strain tensor
+
+
+		Real epsilonxx = epsilon(0, 0);
+		Real epsilonyy = epsilon(1, 1);
+		Real epsilonzz = 0; 			//z-components zero for 2D measures
+		Real epsilonxy = epsilon(0, 1);
+		Real epsilonxz = 0; 			//z-components zero for 2D measures
+		Real epsilonyz = 0; 			//z-components zero for 2D measures
+
+		return sqrt( (1.0 / 3.0) * (std::pow(epsilonxx - epsilonyy, 2.0) + std::pow(epsilonyy - epsilonzz, 2.0) + std::pow(epsilonzz - epsilonxx, 2.0))
+		 + 2.0 * (std::pow(epsilonxy, 2.0) + std::pow(epsilonyz, 2.0) + std::pow(epsilonxz, 2.0)));
+	}
+	//=================================================================================================//
 }

SPHINXsys/src/for_3D_build/geometries/geometry.cpp

@@ -154,6 +154,30 @@ namespace SPH
 		return complex_shape_mesh_->findSignedDistance(input_pnt);
 	}
 	//=================================================================================================//
+	Vec3d ComplexShape::findNormalDirectionComplexShape(const Vec3d &input_pnt) //function to differentiate from LevelSetComplexShape::findNormalDirection
+	{
+		bool is_contain = checkContain(input_pnt);
+		Vecd displacement_to_surface = findClosestPoint(input_pnt) - input_pnt;
+		if(input_pnt[2] < 0.01)
+		{
+			if (input_pnt[0] > 0.01 && input_pnt[0] < 0.09 &&
+			input_pnt[1] > 0.01 && input_pnt[1] < 0.09)
+			{
+				std::cout << "displacement_to_surface: " << displacement_to_surface << std::endl;
+			}
+		}
+		while (displacement_to_surface.norm() < Eps)
+		{
+			Vecd jittered = input_pnt; //jittering
+			for (int l = 0; l != input_pnt.size(); ++l)
+				jittered[l] = input_pnt[l] + (((Real)rand() / (RAND_MAX)) - 0.5) * 100.0 * Eps;
+			if (checkContain(jittered) == is_contain)
+				displacement_to_surface = findClosestPoint(jittered) - jittered;
+		}
+		Vecd direction_to_surface = displacement_to_surface.normalize();
+		return is_contain ? direction_to_surface : -1.0 * direction_to_surface;
+	}
+		//=================================================================================================//
 	Vec3d ComplexShape::findNormalDirection(const Vec3d &input_pnt)
 	{
 		return complex_shape_mesh_->findNormalDirection(input_pnt);

SPHINXsys/src/for_3D_build/geometries/geometry.h

@@ -127,6 +127,7 @@ namespace SPH
 		virtual Real findSignedDistance(const Vec3d &input_pnt);
 		/** Normal direction point toward outside of the complex shape. */
 		virtual Vec3d findNormalDirection(const Vec3d &input_pnt);
+		Vec3d findNormalDirectionComplexShape(const Vec3d &input_pnt); //function to differentiate from LevelSetComplexShape::findNormalDirection
 
 	protected:
 		/** shape container<pointer to geomtry, operation> */

SPHINXsys/src/for_3D_build/particles/solid_particles_supplementary.cpp

@@ -31,4 +31,34 @@ namespace SPH {
 			+ 3.0 * (sigmaxy * sigmaxy + sigmaxz * sigmaxz + sigmayz * sigmayz));
 	}
 	//=================================================================================================//
+	Vecd ElasticSolidParticles::displacement(size_t particle_i)
+	{
+		Vecd disp = pos_n_[particle_i] - pos_0_[particle_i];
+		return disp;
+	}
+	//=================================================================================================//
+	Vecd ElasticSolidParticles::normal(size_t particle_i)
+	{
+		Vecd normal_vec = n_[particle_i];
+		return normal_vec;
+	}
+	//=================================================================================================//
+	Real ElasticSolidParticles::von_Mises_strain(size_t particle_i)
+	{
+
+		Mat3d F = F_[particle_i];
+		Mat3d epsilon = 0.5 * (~F * F - Matd(1.0)); //calculation of the Green-Lagrange strain tensor
+
+
+		Real epsilonxx = epsilon(0, 0);
+		Real epsilonyy = epsilon(1, 1);
+		Real epsilonzz = epsilon(2, 2);
+		Real epsilonxy = epsilon(0, 1);
+		Real epsilonxz = epsilon(0, 2);
+		Real epsilonyz = epsilon(1, 2);
+
+		return sqrt( (1.0 / 3.0) * (std::pow(epsilonxx - epsilonyy, 2.0) + std::pow(epsilonyy - epsilonzz, 2.0) + std::pow(epsilonzz - epsilonxx, 2.0))
+		 + 2.0 * (std::pow(epsilonxy, 2.0) + std::pow(epsilonyz, 2.0) + std::pow(epsilonxz, 2.0)));
+	}
+	//=================================================================================================//
 }

SPHINXsys/src/shared/geometries/geometry_level_set.cpp

@@ -49,6 +49,7 @@ namespace SPH
 	//=================================================================================================//
 	Vecd LevelSetComplexShape::findNormalDirection(const Vecd &input_pnt)
 	{
+		//std::cout << "LevelSetComplexShape::findNormalDirection called" << std::endl; //to check if LevelSetComplexShape::findNormalDirection is called
 		return level_set_->probeNormalDirection(input_pnt);
 	}
 	//=================================================================================================//

SPHINXsys/src/shared/particles/base_particles.cpp

@@ -37,14 +37,14 @@ namespace SPH
 		//		add basic output particle data
 		//----------------------------------------------------------------------
 		addAVariableToWrite<indexVector, Vecd>("Velocity");
-		addAVariableToWrite<indexScalar, Real>("Density");
+		addAVariableToWrite<indexVector, Vecd>("Acceleration");
 		//----------------------------------------------------------------------
 		//		add restart output particle data
 		//----------------------------------------------------------------------
 		addAVariableNameToList<indexVector, Vecd>(variables_to_restart_, "Position");
 		addAVariableNameToList<indexVector, Vecd>(variables_to_restart_, "Velocity");
+		addAVariableNameToList<indexVector, Vecd>(variables_to_restart_, "Acceleration");
 		addAVariableNameToList<indexScalar, Real>(variables_to_restart_, "Volume");
-		addAVariableNameToList<indexScalar, Real>(variables_to_restart_, "Density");
 
 		ParticleGenerator* particle_generator = body_->particle_generator_;
 		particle_generator->initialize(body_);
@@ -135,7 +135,7 @@ namespace SPH
 	{
 		size_t total_real_particles = total_real_particles_;
 
-		//write particle positions first
+		//write current/final particle positions first
 		output_file << "   <Points>\n";
 		output_file << "    <DataArray Name=\"Position\" type=\"Float32\"  NumberOfComponents=\"3\" Format=\"ascii\">\n";
 		output_file << "    ";

SPHINXsys/src/shared/particles/solid_particles.cpp

@@ -61,10 +61,12 @@ namespace SPH {
 	//=================================================================================================//
 	void SolidParticles::initializeNormalDirectionFromGeometry()
 	{
+		std::cout << "initialize function called" << std::endl; //to check if function is called
 		ComplexShape* body_shape = body_->body_shape_;
 		for (size_t i = 0; i != total_real_particles_; ++i)
 		{
-			Vecd normal_direction = body_shape->findNormalDirection(pos_n_[i]);
+			//Vecd normal_direction = body_shape->findNormalDirection(pos_n_[i]); //calls LevelSetComplexShape::findNormalDirection
+			Vecd normal_direction = body_shape->findNormalDirectionComplexShape(pos_n_[i]); //function to differentiate from LevelSetComplexShape::findNormalDirection
 			n_[i] = normal_direction;
 			n_0_[i] = normal_direction;
 		}
@@ -120,33 +122,114 @@ namespace SPH {
 		return von_Mises_stress_max;
 	}
 	//=================================================================================================//
+	StdLargeVec<Vecd> ElasticSolidParticles::getDisplacement()
+	{
+		StdLargeVec<Vecd> displacement_vector = {};
+		for (size_t index_i = 0; index_i < pos_0_.size(); index_i++)
+		{
+			displacement_vector.push_back(displacement(index_i));
+		}
+		return displacement_vector;
+	}
+	//=================================================================================================//
+	StdLargeVec<Vecd> ElasticSolidParticles::getNormal()
+	{
+		StdLargeVec<Vecd> normal_vector = {};
+		for (size_t index_i = 0; index_i < pos_0_.size(); index_i++)
+		{
+			normal_vector.push_back(normal(index_i));
+		}
+		return normal_vector;
+	}
+	//=================================================================================================//
+	StdLargeVec<Real> ElasticSolidParticles::getVonMisesStrain()
+	{
+		StdLargeVec<Real> von_Mises_strain_vector = {};
+		for (size_t index_i = 0; index_i < pos_0_.size(); index_i++)
+		{
+			von_Mises_strain_vector.push_back(von_Mises_strain(index_i));
+		}
+		return von_Mises_strain_vector;
+	}
+	//=================================================================================================//
+	Real ElasticSolidParticles::getMaxVonMisesStrain()
+	{
+		Real von_Mises_strain_max = 0;
+		for (size_t index_i = 0; index_i < pos_0_.size(); index_i++)
+		{
+			Real von_Mises_strain_i = von_Mises_strain(index_i);
+			if (von_Mises_strain_max < von_Mises_strain_i)
+			{
+				von_Mises_strain_max = von_Mises_strain_i;
+			}
+		}
+		return von_Mises_strain_max;
+	}
+	//=================================================================================================//
 	void ElasticSolidParticles::writeParticlesToVtuFile(std::ofstream& output_file)
 	{
 		SolidParticles::writeParticlesToVtuFile(output_file);
 
 		size_t total_real_particles = total_real_particles_;
 
+		//write von Mises stress
 		output_file << "    <DataArray Name=\"von Mises stress\" type=\"Float32\" Format=\"ascii\">\n";
 		output_file << "    ";
 		for (size_t i = 0; i != total_real_particles; ++i) {
 			output_file << std::fixed << std::setprecision(9) << von_Mises_stress(i) << " ";
 		}
 		output_file << std::endl;
 		output_file << "    </DataArray>\n";
+
+		//write Displacement
+		output_file << "    <DataArray Name=\"Displacement\" type=\"Float32\" NumberOfComponents=\"3\" Format=\"ascii\">\n";
+		output_file << "    ";
+		for (size_t i = 0; i != total_real_particles; ++i) {
+			Vecd displacement_vector = displacement(i);
+			output_file << displacement_vector[0] << " " << displacement_vector[1] << " " << displacement_vector[2] << " ";
+		}
+		output_file << std::endl;
+		output_file << "    </DataArray>\n";
+
+		//write Normal Vectors
+		output_file << "    <DataArray Name=\"Normal Vector\" type=\"Float32\" NumberOfComponents=\"3\" Format=\"ascii\">\n";
+		output_file << "    ";
+		for (size_t i = 0; i != total_real_particles; ++i) {
+			Vecd normal_vector = normal(i);
+			output_file << normal_vector[0] << " " << normal_vector[1] << " " << normal_vector[2] << " ";
+		}
+		output_file << std::endl;
+		output_file << "    </DataArray>\n";
+
+		//write von Mises strain
+		output_file << "    <DataArray Name=\"von Mises strain\" type=\"Float32\" Format=\"ascii\">\n";
+		output_file << "    ";
+		for (size_t i = 0; i != total_real_particles; ++i) {
+			output_file << std::fixed << std::setprecision(9) << von_Mises_strain(i) << " ";
+		}
+		output_file << std::endl;
+		output_file << "    </DataArray>\n";
 	}
 	//=================================================================================================//
 	void ElasticSolidParticles::writePltFileHeader(std::ofstream& output_file)
 	{
 		SolidParticles::writePltFileHeader(output_file);
 
 		output_file << ",\" von Mises stress \"";
+		output_file << ",\" Displacement \"";
+		output_file << ",\" Normal Vectors \"";
+		output_file << ",\" von Mises strain \"";
 	}
 	//=================================================================================================//
 	void ElasticSolidParticles::writePltFileParticleData(std::ofstream& output_file, size_t index_i)
 	{
 		SolidParticles::writePltFileParticleData(output_file, index_i);
 
 		output_file << von_Mises_stress(index_i) << " ";
+		Vecd displacement_vector = displacement(index_i);
+		output_file << displacement_vector[0] << " " << displacement_vector[1] << " " << displacement_vector[2] << " "
+			<< index_i << " ";
+		output_file << von_Mises_strain(index_i) << " ";
 	}
 	//=============================================================================================//
 	void ActiveMuscleParticles::initializeActiveMuscleParticleData()

SPHINXsys/src/shared/particles/solid_particles.h

@@ -104,6 +104,19 @@ namespace SPH {
 		StdLargeVec<Real> getVonMisesStress();
 		Real getMaxVonMisesStress();
 
+		/**< Computing displacemnt. */
+		Vecd displacement(size_t particle_i);
+		StdLargeVec<Vecd> getDisplacement();
+
+		/**< Computing normal vector. */
+		Vecd normal (size_t particle_i);
+		StdLargeVec<Vecd> getNormal();
+
+		/**< Computing von Mises equivalent stress. */
+		Real von_Mises_strain (size_t particle_i);
+		StdLargeVec<Real> getVonMisesStrain();
+		Real getMaxVonMisesStrain();
+
 		virtual void writeParticlesToVtuFile(std::ofstream &output_file) override;
 
 		virtual ElasticSolidParticles* ThisObjectPtr() override {return this;};