Xiangyu/composite solid

src/shared/materials/complex_solid.cpp

@@ -0,0 +1,57 @@
+/**
+ * @file 	complex_solid.cpp
+ * @brief 	These are classes for define complex solid materials.
+ * @author	Yaru Ren and Xiangyu Hu
+ */
+
+#include "complex_solid.h"
+
+#include "base_local_dynamics.h"
+#include "base_particles.hpp"
+
+namespace SPH
+{
+//=================================================================================================//
+CompositeSolid::CompositeSolid(Real rho0)
+    : ElasticSolid(rho0)
+{
+    material_type_name_ = "CompositeSolid";
+}
+//=================================================================================================//
+Matd CompositeSolid::StressPK2(Matd &deformation, size_t index_i)
+{
+    return composite_materials_[material_id_[index_i]]->StressPK2(deformation, index_i);
+}
+//=================================================================================================//
+Matd CompositeSolid::StressPK1(Matd &deformation, size_t index_i)
+{
+    return composite_materials_[material_id_[index_i]]->StressPK1(deformation, index_i);
+}
+//=================================================================================================//
+Matd CompositeSolid::StressCauchy(Matd &almansi_strain, Matd &F, size_t index_i)
+{
+    return Matd::Identity();
+}
+//=================================================================================================//
+Real CompositeSolid::CompositeDensity(size_t index_i)
+{
+    return composite_materials_[material_id_[index_i]]->ReferenceDensity();
+}
+//=================================================================================================//
+void CompositeSolid::initializeLocalParameters(BaseParticles *base_particles)
+{
+    ElasticSolid::initializeLocalParameters(base_particles);
+    base_particles->registerVariable(material_id_, "MaterialID");
+
+    for (size_t i = 0; i < composite_materials_.size(); ++i)
+    {
+        composite_materials_[i]->initializeLocalParameters(base_particles);
+    }
+}
+//=================================================================================================//
+MaterialIdInitialization::MaterialIdInitialization(SPHBody &sph_body)
+    : LocalDynamics(sph_body), GeneralDataDelegateSimple(sph_body),
+      material_id_(*particles_->getVariableByName<int>("MaterialID")),
+      pos0_(*particles_->getVariableByName<Vecd>("InitialPosition")){};
+//=================================================================================================//
+} // namespace SPH

src/shared/materials/complex_solid.h

@@ -26,9 +26,11 @@
  * @author	Chi ZHang and Xiangyu Hu
  */
 
-#pragma once
+#ifndef COMPLEX_SOLID_H
+#define COMPLEX_SOLID_H
 
 #include "elastic_solid.h"
+#include <general_dynamics.h>
 
 namespace SPH
 {
@@ -53,4 +55,55 @@ class ActiveMuscle : public MuscleType
     virtual Matd StressPK2(Matd &deformation, size_t index_i) override;
     virtual ActiveMuscle<MuscleType> *ThisObjectPtr() override { return this; };
 };
+
+/**
+ * @class CompositeSolid
+ */
+class CompositeSolid : public ElasticSolid
+{
+    StdLargeVec<int> material_id_;
+
+  protected:
+    UniquePtrsKeeper<ElasticSolid> composite_ptrs_keeper_;
+    StdVec<ElasticSolid *> composite_materials_;
+
+  public:
+    explicit CompositeSolid(Real rho0);
+    virtual ~CompositeSolid(){};
+
+    virtual void initializeLocalParameters(BaseParticles *base_particles) override;
+    virtual Matd StressPK2(Matd &deformation, size_t index_i) override;
+    virtual Matd StressPK1(Matd &deformation, size_t index_i) override;
+    virtual Matd StressCauchy(Matd &almansi_strain, Matd &F, size_t index_i) override;
+    virtual Real VolumetricKirchhoff(Real J) override { return 0.0; };
+    virtual std::string getRelevantStressMeasureName() override { return "PK2"; };
+
+    Real CompositeDensity(size_t index_i);
+
+    template <class ElasticSolidType, typename... Args>
+    void add(Args &&...args)
+    {
+        ElasticSolid *added_material =
+            composite_ptrs_keeper_.createPtr<ElasticSolidType>(std::forward<Args>(args)...);
+        composite_materials_.push_back(added_material);
+        c0_ = SMAX(c0_, added_material->ReferenceSoundSpeed());
+        setContactStiffness(c0_);
+    };
+};
+
+/**
+ * @class MaterialIdInitialization
+ */
+class MaterialIdInitialization
+    : public LocalDynamics,
+      public GeneralDataDelegateSimple
+{
+  public:
+    explicit MaterialIdInitialization(SPHBody &sph_body);
+
+  protected:
+    StdLargeVec<int> &material_id_;
+    StdLargeVec<Vecd> &pos0_;
+};
 } // namespace SPH
+#endif // COMPLEX_SOLID_H

src/shared/materials/elastic_solid.cpp

@@ -53,13 +53,18 @@ Real LinearElasticSolid::getLambda(Real youngs_modulus, Real poisson_ratio)
     return nu_ * youngs_modulus / (1.0 + poisson_ratio) / (1.0 - 2.0 * poisson_ratio);
 }
 //=================================================================================================//
-Matd LinearElasticSolid::StressPK2(Matd &F, size_t particle_index_i)
+Matd LinearElasticSolid::StressPK1(Matd &F, size_t index_i)
+{
+    return F * StressPK2(F, index_i);
+}
+//=================================================================================================//
+Matd LinearElasticSolid::StressPK2(Matd &F, size_t index_i)
 {
     Matd strain = 0.5 * (F.transpose() + F) - Matd::Identity();
     return lambda0_ * strain.trace() * Matd::Identity() + 2.0 * G0_ * strain;
 }
 //=================================================================================================//
-Matd LinearElasticSolid::StressCauchy(Matd &almansi_strain, Matd &F, size_t particle_index_i)
+Matd LinearElasticSolid::StressCauchy(Matd &almansi_strain, Matd &F, size_t index_i)
 {
     return lambda0_ * almansi_strain.trace() * Matd::Identity() + 2.0 * G0_ * almansi_strain;
 }
@@ -69,13 +74,13 @@ Real LinearElasticSolid::VolumetricKirchhoff(Real J)
     return K0_ * J * (J - 1);
 }
 //=================================================================================================//
-Matd SaintVenantKirchhoffSolid::StressPK2(Matd &F, size_t particle_index_i)
+Matd SaintVenantKirchhoffSolid::StressPK2(Matd &F, size_t index_i)
 {
     Matd strain = 0.5 * (F.transpose() * F - Matd::Identity());
     return lambda0_ * strain.trace() * Matd::Identity() + 2.0 * G0_ * strain;
 }
 //=================================================================================================//
-Matd NeoHookeanSolid::StressPK2(Matd &F, size_t particle_index_i)
+Matd NeoHookeanSolid::StressPK2(Matd &F, size_t index_i)
 {
     // This formulation allows negative determinant of F. Please refer Eq. (12) in
     // Smith et al. (2018) Stable Neo-Hookean Flesh Simulation.
@@ -85,7 +90,7 @@ Matd NeoHookeanSolid::StressPK2(Matd &F, size_t particle_index_i)
     return G0_ * Matd::Identity() + (lambda0_ * (J - 1.0) - G0_) * J * right_cauchy.inverse();
 }
 //=================================================================================================//
-Matd NeoHookeanSolid::StressCauchy(Matd &almansi_strain, Matd &F, size_t particle_index_i)
+Matd NeoHookeanSolid::StressCauchy(Matd &almansi_strain, Matd &F, size_t index_i)
 {
     Real J = F.determinant();
     Matd B = (-2.0 * almansi_strain + Matd::Identity()).inverse();
@@ -100,7 +105,7 @@ Real NeoHookeanSolid::VolumetricKirchhoff(Real J)
     return 0.5 * K0_ * (J * J - 1);
 }
 //=================================================================================================//
-Matd NeoHookeanSolidIncompressible::StressPK2(Matd &F, size_t particle_index_i)
+Matd NeoHookeanSolidIncompressible::StressPK2(Matd &F, size_t index_i)
 {
     Matd right_cauchy = F.transpose() * F;
     Real I_1 = right_cauchy.trace();       // first strain invariant
@@ -109,7 +114,7 @@ Matd NeoHookeanSolidIncompressible::StressPK2(Matd &F, size_t particle_index_i)
 }
 //=================================================================================================//
 Matd NeoHookeanSolidIncompressible::
-    StressCauchy(Matd &almansi_strain, Matd &F, size_t particle_index_i)
+    StressCauchy(Matd &almansi_strain, Matd &F, size_t index_i)
 {
     // TODO: implement
     return {};
@@ -136,7 +141,7 @@ OrthotropicSolid::OrthotropicSolid(Real rho_0, std::array<Vecd, Dimensions> a, s
     CalculateAllLambda();
 };
 //=================================================================================================//
-Matd OrthotropicSolid::StressPK2(Matd &F, size_t particle_index_i)
+Matd OrthotropicSolid::StressPK2(Matd &F, size_t index_i)
 {
     Matd strain = 0.5 * (F.transpose() * F - Matd::Identity());
     Matd stress_PK2 = Matd::Zero();
@@ -166,7 +171,7 @@ void OrthotropicSolid::CalculateA0()
         A_[i] = a_[i] * a_[i].transpose();
 }
 //=================================================================================================//
-Matd FeneNeoHookeanSolid::StressPK2(Matd &F, size_t particle_index_i)
+Matd FeneNeoHookeanSolid::StressPK2(Matd &F, size_t index_i)
 {
     Matd right_cauchy = F.transpose() * F;
     Matd strain = 0.5 * (right_cauchy - Matd::Identity());

src/shared/materials/elastic_solid.h

@@ -70,6 +70,8 @@ class ElasticSolid : public Solid
     Real BulkModulus() { return K0_; };
     Real PoissonRatio() { return nu_; };
 
+    /** 1st Piola-Kirchhoff stress through deformation. */
+    virtual Matd StressPK1(Matd &deformation, size_t particle_index_i) = 0;
     /** 2nd Piola-Kirchhoff stress through deformation. */
     virtual Matd StressPK2(Matd &deformation, size_t particle_index_i) = 0;
     /** Cauchy stress through Eulerian Almansi strain tensor. */
@@ -115,6 +117,7 @@ class LinearElasticSolid : public ElasticSolid
     explicit LinearElasticSolid(Real rho0, Real youngs_modulus, Real poisson_ratio);
     virtual ~LinearElasticSolid(){};
 
+    virtual Matd StressPK1(Matd &deformation, size_t particle_index_i) override;
     virtual Matd StressPK2(Matd &deformation, size_t particle_index_i) override;
     virtual Matd StressCauchy(Matd &almansi_strain, Matd &F, size_t particle_index_i) override;
     /** Volumetric Kirchhoff stress from determinate */

src/shared/particle_dynamics/solid_dynamics/elastic_dynamics.cpp

@@ -70,8 +70,7 @@ Integration1stHalf::
     numerical_dissipation_factor_ = 0.25;
 }
 //=================================================================================================//
-Integration1stHalfPK2::
-    Integration1stHalfPK2(BaseInnerRelation &inner_relation)
+Integration1stHalfPK2::Integration1stHalfPK2(BaseInnerRelation &inner_relation)
     : Integration1stHalf(inner_relation){};
 //=================================================================================================//
 void Integration1stHalfPK2::initialization(size_t index_i, Real dt)
@@ -81,7 +80,7 @@ void Integration1stHalfPK2::initialization(size_t index_i, Real dt)
     rho_[index_i] = rho0_ / F_[index_i].determinant();
     // obtain the first Piola-Kirchhoff stress from the second Piola-Kirchhoff stress
     // it seems using reproducing correction here increases convergence rate near the free surface
-    stress_PK1_B_[index_i] = F_[index_i] * elastic_solid_.StressPK2(F_[index_i], index_i) * B_[index_i];
+    stress_PK1_B_[index_i] = elastic_solid_.StressPK1(F_[index_i], index_i) * B_[index_i];
 }
 //=================================================================================================//
 Integration1stHalfKirchhoff::

tests/user_examples/test_2d_flow_stream_around_fish/2d_flow_stream_around_fish.cpp

@@ -4,7 +4,6 @@
  * @author  Yaru Ren and Xiangyu Hu
  */
 #include "2d_flow_stream_around_fish.h"
-#include "active_model.h"
 #include "sphinxsys.h"
 using namespace SPH;
 
@@ -36,7 +35,7 @@ int main(int ac, char *av[])
     SolidBody fish_body(system, makeShared<FishBody>("FishBody"));
     fish_body.defineAdaptationRatios(1.15, 2.0);
     fish_body.defineBodyLevelSetShape()->writeLevelSet(io_environment);
-    fish_body.defineParticlesAndMaterial<ElasticSolidParticles, SolidBodyMaterial>();
+    fish_body.defineParticlesAndMaterial<ElasticSolidParticles, FishBodyComposite>();
     //  Using relaxed particle distribution if needed
     (!system.RunParticleRelaxation() && system.ReloadParticles())
         ? fish_body.generateParticles<ParticleGeneratorReload>(io_environment, fish_body.getName())
@@ -151,7 +150,7 @@ int main(int ac, char *av[])
     SimpleDynamics<ImposingActiveStrain> imposing_active_strain(fish_body);
     ReduceDynamics<solid_dynamics::AcousticTimeStepSize> fish_body_computing_time_step_size(fish_body);
     /** Stress relaxation for the inserted body. */
-    Dynamics1Level<solid_dynamics::ActiveIntegration1stHalf> fish_body_stress_relaxation_first_half(fish_inner);
+    Dynamics1Level<solid_dynamics::Integration1stHalfPK2> fish_body_stress_relaxation_first_half(fish_inner);
     Dynamics1Level<solid_dynamics::Integration2ndHalf> fish_body_stress_relaxation_second_half(fish_inner);
     /** Update norm .*/
     SimpleDynamics<solid_dynamics::UpdateElasticNormalDirection> fish_body_update_normal(fish_body);

tests/user_examples/test_2d_flow_stream_around_fish/2d_flow_stream_around_fish.h

@@ -1,8 +1,8 @@
 
 #include "2d_fish_and_bones.h"
-#include "composite_material.h"
+#include "active_model.h"
 #include "sphinxsys.h"
-#define PI 3.1415926
+
 using namespace SPH;
 //----------------------------------------------------------------------
 //	Basic geometry parameters and numerical setup.
@@ -147,10 +147,10 @@ class TimeDependentAcceleration : public Gravity
 };
 
 // Material ID
-class SolidBodyMaterial : public CompositeMaterial
+class FishBodyComposite : public CompositeSolid
 {
   public:
-    SolidBodyMaterial() : CompositeMaterial(rho0_s)
+    FishBodyComposite() : CompositeSolid(rho0_s)
     {
         add<ActiveModelSolid>(rho0_s, Youngs_modulus1, poisson);
         add<SaintVenantKirchhoffSolid>(rho0_s, Youngs_modulus2, poisson);
@@ -200,10 +200,8 @@ class ImposingActiveStrain
     explicit ImposingActiveStrain(SolidBody &solid_body)
         : solid_dynamics::ElasticDynamicsInitialCondition(solid_body),
           material_id_(*particles_->getVariableByName<int>("MaterialID")),
-          pos0_(*particles_->getVariableByName<Vecd>("InitialPosition"))
-    {
-        particles_->registerVariable(active_strain_, "ActiveStrain");
-    };
+          pos0_(*particles_->getVariableByName<Vecd>("InitialPosition")),
+          active_strain_(*particles_->getVariableByName<Matd>("ActiveStrain")){};
     virtual void update(size_t index_i, Real dt = 0.0)
     {
         if (material_id_[index_i] == 0)
@@ -230,5 +228,5 @@ class ImposingActiveStrain
   protected:
     StdLargeVec<int> &material_id_;
     StdLargeVec<Vecd> &pos0_;
-    StdLargeVec<Matd> active_strain_;
+    StdLargeVec<Matd> &active_strain_;
 };

tests/user_examples/test_2d_flow_stream_around_fish/active_model.cpp

@@ -1,40 +1,36 @@
 #include "active_model.h"
-#include <numeric>
+
+#include "base_particles.hpp"
 
 namespace SPH
 {
-	//=========================================================================================================//
-	namespace solid_dynamics
-	{
-		//=================================================================================================//
-		ActiveIntegration1stHalf::
-			ActiveIntegration1stHalf(BaseInnerRelation &inner_relation)
-			: Integration1stHalfPK2(inner_relation), active_strain_(*particles_->getVariableByName<Matd>("ActiveStrain"))
-			, material_id_(*particles_->getVariableByName<int>("MaterialID"))
-		{
-			particles_->registerVariable(F_0, "ActiveTensor");
-			particles_->registerVariable(E_e, "ElasticStrain");
-		}
-		//=================================================================================================//
-		void ActiveIntegration1stHalf::initialization(size_t index_i, Real dt)
-		{
-			Integration1stHalfPK2::initialization(index_i,dt);
-
-			// Calculate the active deformation gradient
-			F_0[index_i] = Matd::Identity();
-			F_0[index_i] = (material_id_[index_i] == 0) ? ((2.0 * active_strain_[index_i] + Matd::Identity()).llt().matrixL()) : F_0[index_i];
-
-			//  Calculate the elastic deformation 
-			Matd F_e = F_[index_i] * F_0[index_i].inverse();
-			Matd F_0_star = F_0[index_i].determinant() * (F_0[index_i].inverse()).transpose();
+//=================================================================================================//
+ActiveModelSolid::ActiveModelSolid(Real rho0, Real youngs_modulus, Real poisson_ratio)
+    : SaintVenantKirchhoffSolid(rho0, youngs_modulus, poisson_ratio)
+{
+    material_type_name_ = "ActiveModelSolid";
+}
+//=============================================================================================//
+void ActiveModelSolid::initializeLocalParameters(BaseParticles *base_particles)
+{
+    SaintVenantKirchhoffSolid::initializeLocalParameters(base_particles);
+    base_particles->registerVariable(active_strain_, "ActiveStrain");
+}
+//=================================================================================================//
+Matd ActiveModelSolid::StressPK1(Matd &F, size_t index_i)
+{
+    // Calculate the active deformation gradient
+    Matd F0 = (2.0 * active_strain_[index_i] + Matd::Identity()).llt().matrixL();
+    Matd F0_inv = F0.inverse();
 
-			// Calculate the elastic strain
-			E_e[index_i] = (material_id_[index_i] == 0) ? (0.5 * (F_[index_i].transpose() * F_[index_i] - Matd::Identity()) - active_strain_[index_i]) : F_[index_i];
+    // Calculate the elastic deformation
+    Matd F_e = F * F0_inv;
+    Matd F0_star = F0.determinant() * F0_inv.transpose();
 
-			stress_PK1_B_[index_i] = F_e * elastic_solid_.StressPK2(E_e[index_i], index_i) * F_0_star * B_[index_i];
-		}
-		//=================================================================================================//
-	}
-	//=================================================================================================//
+    // Calculate the elastic strain
+    Matd E_e = 0.5 * (F.transpose() * F - Matd::Identity()) - active_strain_[index_i];
+    return F_e * (lambda0_ * E_e.trace() * Matd::Identity() + 2.0 * G0_ * E_e) * F0_star;
 }
 //=================================================================================================//
+} // namespace SPH
+//=================================================================================================//