Updated remodeling interface materials to calculate solid mass densit…

…y for elastic solids, solid mixtures, and multiphasic materials
febiosoftware · Apr 20, 2024 · cccf8cc · cccf8cc
1 parent 1747162
commit cccf8cc
Show file tree

Hide file tree

Showing 9 changed files with 71 additions and 3 deletions.
diff --git a/FEBioMech/FEBioMechPlot.cpp b/FEBioMech/FEBioMechPlot.cpp
@@ -1407,8 +1407,8 @@ bool FEPlotDensity::Save(FEDomain &dom, FEDataStream& a)
 	if (dom.Class() == FE_DOMAIN_SOLID)
 	{
 		FESolidDomain& bd = static_cast<FESolidDomain&>(dom);
-		FEElasticMaterial* em = dynamic_cast<FEElasticMaterial*>(bd.GetMaterial());
-		if (em == 0) return false;
+		FEElasticMaterial* em = bd.GetMaterial()->ExtractProperty<FEElasticMaterial>();
+        if (em == 0) return false;
 
 		FERemodelingElasticMaterial* rm = dynamic_cast<FERemodelingElasticMaterial*>(em);
 		if (rm)
@@ -1426,7 +1426,7 @@ bool FEPlotDensity::Save(FEDomain &dom, FEDataStream& a)
 	}
 	else if (dom.Class() == FE_DOMAIN_SHELL)
 	{
-		FEElasticMaterial* em = dynamic_cast<FEElasticMaterial*>(dom.GetMaterial());
+        FEElasticMaterial* em = dom.GetMaterial()->ExtractProperty<FEElasticMaterial>();
 		if (em == 0) return false;
 		FEDensity dens(em);
 		writeAverageElementValue<double>(dom, a, dens);

diff --git a/FEBioMech/FERemodelingElasticMaterial.cpp b/FEBioMech/FERemodelingElasticMaterial.cpp
@@ -95,6 +95,14 @@ double FERemodelingElasticMaterial::StrainEnergyDensity(FEMaterialPoint& mp)
 	return (dynamic_cast<FERemodelingInterface*>((FEElasticMaterial*)m_pBase))->StrainEnergy(mp);
 }
 
+//-----------------------------------------------------------------------------
+//! evaluate referential mass density
+double FERemodelingElasticMaterial::Density(FEMaterialPoint& mp)
+{
+    FERemodelingMaterialPoint& rpt = *(mp.ExtractData<FERemodelingMaterialPoint>());
+    return rpt.m_rhor;
+}
+
 //-----------------------------------------------------------------------------
 //! Stress function
 mat3ds FERemodelingElasticMaterial::Stress(FEMaterialPoint& mp)

diff --git a/FEBioMech/FERemodelingElasticMaterial.h b/FEBioMech/FERemodelingElasticMaterial.h
@@ -111,6 +111,9 @@ class FEBIOMECH_API FERemodelingElasticMaterial : public FEElasticMaterial
 
 	//! tangent function of stress with strain
 	tens4ds Tangent(FEMaterialPoint& pt) override;
+
+    //! evaluate referential mass density
+    double Density(FEMaterialPoint& pt) override;
 
 	//! tangent function of strain energy density with solid mass density
 	double Tangent_SE_Density(FEMaterialPoint& pt);

diff --git a/FEBioMix/FECarterHayes.cpp b/FEBioMix/FECarterHayes.cpp
@@ -190,6 +190,22 @@ tens4ds FECarterHayes::Tangent(FEMaterialPoint& mp)
 	return tens4ds(D);
 }
 
+//-----------------------------------------------------------------------------
+//! evaluate referential mass density
+double FECarterHayes::Density(FEMaterialPoint& pt)
+{
+    FERemodelingMaterialPoint* rpt = pt.ExtractData<FERemodelingMaterialPoint>();
+    if (rpt) return rpt->m_rhor;
+    else {
+        FEElasticMixtureMaterialPoint* emp = pt.ExtractData<FEElasticMixtureMaterialPoint>();
+        if (emp) {
+            rpt = emp->GetPointData(m_comp)->ExtractData<FERemodelingMaterialPoint>();
+            if (rpt) return rpt->m_rhor;
+        }
+    }
+    return 0.0;
+}
+
 //-----------------------------------------------------------------------------
 //! calculate tangent of strain energy density with mass density
 double FECarterHayes::Tangent_SE_Density(FEMaterialPoint& mp)

diff --git a/FEBioMix/FECarterHayes.h b/FEBioMix/FECarterHayes.h
@@ -63,6 +63,9 @@ class FEBIOMIX_API FECarterHayes : public FEElasticMaterial, public FERemodeling
 	//! calculate tangent stiffness at material point
 	tens4ds Tangent(FEMaterialPoint& pt) override;
 
+    //! evaluate referential mass density
+    double Density(FEMaterialPoint& pt) override;
+
 	//! Create material point data
 	FEMaterialPointData* CreateMaterialPointData() override;
 

diff --git a/FEBioMix/FEFiberExpPowSBM.cpp b/FEBioMix/FEFiberExpPowSBM.cpp
@@ -255,6 +255,22 @@ double FEFiberExpPowSBM::StrainEnergyDensity(FEMaterialPoint& mp)
     return sed;
 }
 
+//-----------------------------------------------------------------------------
+//! evaluate referential mass density
+double FEFiberExpPowSBM::Density(FEMaterialPoint& pt)
+{
+    FERemodelingMaterialPoint* rpt = pt.ExtractData<FERemodelingMaterialPoint>();
+    if (rpt) return rpt->m_rhor;
+    else {
+        FEElasticMixtureMaterialPoint* emp = pt.ExtractData<FEElasticMixtureMaterialPoint>();
+        if (emp) {
+            rpt = emp->GetPointData(m_comp)->ExtractData<FERemodelingMaterialPoint>();
+            if (rpt) return rpt->m_rhor;
+        }
+    }
+    return 0.0;
+}
+
 //-----------------------------------------------------------------------------
 //! calculate strain energy density at material point
 double FEFiberExpPowSBM::StrainEnergy(FEMaterialPoint& mp)

diff --git a/FEBioMix/FEFiberExpPowSBM.h b/FEBioMix/FEFiberExpPowSBM.h
@@ -53,6 +53,9 @@ class FEBIOMIX_API FEFiberExpPowSBM : public FEElasticMaterial, public FERemodel
     //! Strain energy density
     double StrainEnergyDensity(FEMaterialPoint& mp) override;
 
+    //! evaluate referential mass density
+    double Density(FEMaterialPoint& pt) override;
+
     //! return fiber modulus
     double FiberModulus(double rhor) { return m_ksi0*pow(rhor/m_rho0, m_g);}
 

diff --git a/FEBioMix/FEFiberPowLinearSBM.cpp b/FEBioMix/FEFiberPowLinearSBM.cpp
@@ -265,6 +265,22 @@ double FEFiberPowLinearSBM::StrainEnergyDensity(FEMaterialPoint& mp)
     return sed;
 }
 
+//-----------------------------------------------------------------------------
+//! evaluate referential mass density
+double FEFiberPowLinearSBM::Density(FEMaterialPoint& pt)
+{
+    FERemodelingMaterialPoint* rpt = pt.ExtractData<FERemodelingMaterialPoint>();
+    if (rpt) return rpt->m_rhor;
+    else {
+        FEElasticMixtureMaterialPoint* emp = pt.ExtractData<FEElasticMixtureMaterialPoint>();
+        if (emp) {
+            rpt = emp->GetPointData(m_comp)->ExtractData<FERemodelingMaterialPoint>();
+            if (rpt) return rpt->m_rhor;
+        }
+    }
+    return 0.0;
+}
+
 //-----------------------------------------------------------------------------
 //! calculate strain energy density at material point
 double FEFiberPowLinearSBM::StrainEnergy(FEMaterialPoint& mp)

diff --git a/FEBioMix/FEFiberPowLinearSBM.h b/FEBioMix/FEFiberPowLinearSBM.h
@@ -53,6 +53,9 @@ class FEBIOMIX_API FEFiberPowLinearSBM : public FEElasticMaterial, public FERemo
     //! Strain energy density
     double StrainEnergyDensity(FEMaterialPoint& mp) override;
 
+    //! evaluate referential mass density
+    double Density(FEMaterialPoint& pt) override;
+
     //! return fiber modulus
     double FiberModulus(double rhor) { return m_E0*pow(rhor/m_rho0, m_g);}