Updating branch with time changes

FEBioMech/FE3FieldElasticSolidDomain.cpp

@@ -49,6 +49,8 @@ void FE3FieldElasticSolidDomain::ELEM_DATA::Serialize(DumpStream& ar)
 	ar & Lk;
 	ar & eJt;
 	ar & eJp;
+	ar & eJpp;
+	ar & eJ_star;
 }
 
 //-----------------------------------------------------------------------------
@@ -92,7 +94,7 @@ bool FE3FieldElasticSolidDomain::Init()
 	for (int i=0; i<NE; ++i)
 	{
 		ELEM_DATA& d = m_Data[i];
-		d.eJ = d.eJt = d.eJp = 1.0;
+		d.eJ = d.eJt = d.eJp = d.eJ_star = d.eJpp = 1.0;
 		d.ep = 0.0;
 		d.Lk = 0.0;
 	}
@@ -109,7 +111,7 @@ void FE3FieldElasticSolidDomain::Reset()
 	for (size_t i=0; i<NE; ++i)
 	{
 		ELEM_DATA& d = m_Data[i];
-        d.eJ = d.eJt = d.eJp = 1.0;
+        d.eJ = d.eJt = d.eJp = d.eJ_star = d.eJpp = 1.0;
         d.ep = 0.0;
         d.Lk = 0.0;
 	}
@@ -240,12 +242,12 @@ void FE3FieldElasticSolidDomain::ElementDilatationalStiffness(FEModel& fem, int
 	}
 
 	// get effective modulus
-	double k = pmi->UJJ(ed.eJ);
+	double k = pmi->UJJ(ed.eJ, ed.eJ_star);
 
 	// next, we add the Lagrangian contribution
 	// note that this term will always be zero if the material does not
 	// use the augmented lagrangian
-	k += ed.Lk*pmi->hpp(ed.eJ);
+	k += ed.Lk*pmi->hpp(ed.eJ, ed.eJ_star);
 
 	// divide by initial volume
 	k /= Ve;
@@ -448,6 +450,8 @@ void FE3FieldElasticSolidDomain::ElementGeometricalStiffness(int iel, matrix &ke
 //! This function loops over all elements and updates the stress
 void FE3FieldElasticSolidDomain::Update(const FETimeInfo& tp)
 {
+
+	feLog("FE3FieldElasticSolidDomain Update ----- \n");
 	bool berr = false;
 	int NE = (int) m_Elem.size();
 	#pragma omp parallel for shared(NE, berr)
@@ -507,23 +511,28 @@ void FE3FieldElasticSolidDomain::UpdateElementStress(int iel, const FETimeInfo&
 	}
 
 	// calculate the average dilatation and pressure
-	double v = 0, vt = 0, V = 0;
+	double v = 0, vt = 0, V = 0, v_star = 0;
 	for (int n=0; n<nint; ++n)
 	{
+		FEMaterialPoint& mp = *el.GetMaterialPoint(n);
+		FEElasticMaterialPoint& pt = *(mp.ExtractData<FEElasticMaterialPoint>());
+
 		v += detJt(el, n, m_alphaf)*gw[n];
         vt+= detJt(el, n)*gw[n];
 		V += detJ0(el, n)*gw[n];
+		v_star += pt.m_J_star;
 	}
 
 	// calculate volume ratio
+	ed.eJ_star = v_star / nint;
 	ed.eJ = v / V;
     ed.eJt = vt / V;
-    double eUt = mat.U(ed.eJt);
-    double eUp = mat.U(ed.eJp);
+    double eUt = mat.U(ed.eJt, ed.eJ_star);
+    double eUp = mat.U(ed.eJp, ed.eJ_star);
 
 	// Calculate pressure. This is a sum of a Lagrangian term and a penalty term
 	//      <--- Lag. mult. -->  <-- penalty -->
-	ed.ep = ed.Lk*mat.hp(ed.eJ) + mat.UJ(ed.eJ);
+	ed.ep = ed.Lk*mat.hp(ed.eJ, ed.eJ_star) + mat.UJ(ed.eJ, ed.eJ_star);
 //	ed.ep = mat.UJ(ed.eJ);
 
 	// loop over the integration points and calculate
@@ -581,9 +590,9 @@ void FE3FieldElasticSolidDomain::UpdateElementStress(int iel, const FETimeInfo&
 			double Wp = pt.m_Wp;
             mat3ds D = pt.RateOfDeformation();
             double D2 = D.dotdot(D);
-            if (D2 > std::numeric_limits<double>::epsilon())
+            if (D2 > 0)
                 pt.m_s += D*(((Wt-Wp)/(dt*pt.m_J) - pt.m_s.dotdot(D))/D2);
-            if (fabs(ed.eJt - ed.eJp) > std::numeric_limits<double>::epsilon())
+            if (ed.eJt != ed.eJp)
                 pt.m_s += mat3dd((eUt-eUp)/(ed.eJ*(ed.eJt-ed.eJp)));
         }
         else
@@ -615,7 +624,7 @@ bool FE3FieldElasticSolidDomain::Augment(int naug)
 		L0 = ed.Lk;
 		normL0 += L0*L0;
 
-		L1 = L0 + k*pmi->h(ed.eJ);
+		L1 = L0 + k*pmi->h(ed.eJ, ed.eJ_star);
 		normL1 += L1*L1;
 	}
 
@@ -643,9 +652,9 @@ bool FE3FieldElasticSolidDomain::Augment(int naug)
 		{
 			ELEM_DATA& ed = m_Data[n];
 
-			double hi = pmi->h(ed.eJ);
-			ed.Lk += k*pmi->h(ed.eJ);
-			ed.ep = ed.Lk*pmi->hp(ed.eJ) + k*log(ed.eJ)/ed.eJ;
+			double hi = pmi->h(ed.eJ, ed.eJ_star);
+			ed.Lk += k*pmi->h(ed.eJ, ed.eJ_star);
+			ed.ep = ed.Lk*pmi->hp(ed.eJ, ed.eJ_star) + k*log(ed.eJ/ed.eJ_star)/ed.eJ;
 		}
 	}
 

FEBioMech/FE3FieldElasticSolidDomain.h

@@ -42,6 +42,8 @@ class FEBIOMECH_API FE3FieldElasticSolidDomain : public FEElasticSolidDomain
 		double	Lk;		// Lagrangian multiplier
         double  eJt;    // average element jacobian at current time
         double  eJp;    // average element jacobian at previous time
+        double  eJpp;    // average element jacobian at previous time
+		double	eJ_star;// average desired element jacobian at intermediate time
 
 		void Serialize(DumpStream& ar);
 	};

FEBioMech/FEElasticMaterialPoint.cpp

@@ -34,11 +34,13 @@ FEElasticMaterialPoint::FEElasticMaterialPoint(FEMaterialPointData* mp) : FEMate
 {
 	m_F.unit();
 	m_J = 1;
+	m_J_star = 1;
 	m_s.zero();
     m_v = m_a = m_gradJ = vec3d(0, 0, 0);
     m_buncoupled = false;
     m_Wt = m_Wp = 0;
     m_p = 0;
+    m_Lk = 0;
 }
 
 //-----------------------------------------------------------------------------
@@ -55,6 +57,7 @@ void FEElasticMaterialPoint::Init()
 	m_F.unit();
 
 	m_J = 1;
+	m_J_star = 1;
 
 	m_s.zero();
 
@@ -64,6 +67,7 @@ void FEElasticMaterialPoint::Init()
     m_Wt = m_Wp = 0;
 
     m_p = 0;
+    m_Lk = 0;
 
 	// don't forget to initialize the base class
     FEMaterialPointData::Init();
@@ -73,7 +77,7 @@ void FEElasticMaterialPoint::Init()
 void FEElasticMaterialPoint::Serialize(DumpStream& ar)
 {
 	FEMaterialPointData::Serialize(ar);
-    ar & m_F & m_J & m_s & m_v & m_a & m_gradJ & m_L & m_Wt & m_Wp & m_p;
+    ar & m_F & m_J & m_J_star & m_s & m_v & m_a & m_gradJ & m_L & m_Wt & m_Wp & m_p & m_Lk;
 	ar & m_buncoupled;
 }
 

FEBioMech/FEElasticMaterialPoint.h

@@ -84,6 +84,7 @@ class FEBIOMECH_API FEElasticMaterialPoint : public FEMaterialPointData
 	// deformation data at intermediate time
 	mat3d	m_F;	//!< deformation gradient
 	double	m_J;	//!< determinant of F
+	double	m_J_star;	//!< desired J
     vec3d   m_gradJ;  //!< gradient of J
     vec3d   m_v;    //!< velocity
     vec3d   m_a;    //!< acceleration
@@ -100,4 +101,7 @@ class FEBIOMECH_API FEElasticMaterialPoint : public FEMaterialPointData
 
     // previous time data
     double      m_Wp;       //!< strain energy density
+
+    // augmented Lagrangian variable
+    double      m_Lk;       //!< strain energy density
 };

FEBioMech/FEElasticSolidDomain.cpp

@@ -38,6 +38,7 @@ SOFTWARE.*/
 #include "FEBioMech.h"
 #include <FECore/FELinearSystem.h>
 #include "FEResidualVector.h"
+#include "FEUncoupledMaterial.h"
 
 //-----------------------------------------------------------------------------
 //! constructor
@@ -53,6 +54,11 @@ FEElasticSolidDomain::FEElasticSolidDomain(FEModel* pfem) : FESolidDomain(pfem),
 	m_secant_stress = false;
 	m_secant_tangent = false;
 
+	m_blaugon = false;
+	m_augtol = 0.01;
+	m_naugmin = 0;
+	m_naugmax = 0;
+
 	// TODO: Can this be done in Init, since  there is no error checking
 	if (pfem)
 	{
@@ -119,6 +125,14 @@ void FEElasticSolidDomain::Activate()
 	}
 }
 
+//-----------------------------------------------------------------------------
+bool FEElasticSolidDomain::DoAugmentations() const
+{
+	// make sure the domain material uses an uncoupled formulation
+	if (dynamic_cast<FEUncoupledMaterial*>(m_pMat) == 0) return false;
+	return m_blaugon;
+}
+
 //-----------------------------------------------------------------------------
 //! serialization
 void FEElasticSolidDomain::Serialize(DumpStream& ar)
@@ -648,8 +662,7 @@ void FEElasticSolidDomain::UpdateElementStress(int iel, const FETimeInfo& tp)
 		pt.m_s = (m_secant_stress ? m_pMat->SecantStress(mp) : m_pMat->Stress(mp));
 
         // adjust stress for strain energy conservation
-		// (Apply only for mid-point rule)
-		if (m_alphaf == 0.5)
+        if (m_alphaf == 0.5) 
 		{
 			FEElasticMaterial* pme = dynamic_cast<FEElasticMaterial*>(m_pMat);
 
@@ -664,10 +677,8 @@ void FEElasticSolidDomain::UpdateElementStress(int iel, const FETimeInfo& tp)
 
             mat3ds D = pt.RateOfDeformation();
             double D2 = D.dotdot(D);
-			if (D2 > std::numeric_limits<double>::epsilon())
-			{
-				pt.m_s += D * (((pt.m_Wt - pt.m_Wp) / (dt * pt.m_J) - pt.m_s.dotdot(D)) / D2);
-			}
+            if (D2 > 0)
+                pt.m_s += D*(((pt.m_Wt-pt.m_Wp)/(dt*pt.m_J) - pt.m_s.dotdot(D))/D2);
         }
     }
 }
@@ -769,11 +780,101 @@ void FEElasticSolidDomain::ElementInertialForce(FESolidElement& el, vector<doubl
     }
 }
 
+//! Do augmentation
+bool FEElasticSolidDomain::Augment(int naug)
+{
+	FEUncoupledMaterial* pmi = dynamic_cast<FEUncoupledMaterial*>(m_pMat);
+	assert(pmi);
+
+	// make sure Augmented Lagrangian flag is on
+	if (m_blaugon == false) return true;
+
+	// do the augmentation
+	int n;
+	double normL0 = 0, normL1 = 0, L0, L1;
+	double k = pmi->m_K;
+	FEMesh& mesh = *m_pMesh;
+	int NE = Elements();
+
+	for (int iel=0; iel<NE; ++iel)
+	{
+		// get the solid element
+		FESolidElement& el = m_Elem[iel];
+		int nint = el.GaussPoints();
+
+		// loop over the integration points and calculate
+		for (int n=0; n<nint; ++n)
+		{
+			FEMaterialPoint& mp = *el.GetMaterialPoint(n);
+			FEElasticMaterialPoint& pt = *(mp.ExtractData<FEElasticMaterialPoint>());
+
+			L0 = pt.m_Lk;
+			normL0 += L0*L0;
+
+			//L1 = L0 + k*pmi->h(pt.m_J, pt.m_J_star);
+			L1 = L0 + pmi->UJ(pt.m_J, pt.m_J_star);
+			normL1 += L1*L1;
+
+			//printf("%f %f %f %f %f %f\n", mp.m_r0(0), mp.m_r0(1), mp.m_r0(2), pt.m_J, pt.m_J_star, k*pmi->h(pt.m_J, pt.m_J_star));
+		}
+	}
+
+	normL0 = sqrt(normL0);
+	normL1 = sqrt(normL1);
+
+	// check convergence
+	// Note: Use absolute tolerance instead of relative?? or L0_0 at first augmentation relative to that
+	double pctn = 0;
+	if (fabs(normL1) > 1e-10) pctn = fabs((normL1 - normL0)/normL1);
+
+	feLog(" material %d\n", pmi->GetID());
+	feLog("                        CURRENT         CHANGE        REQUIRED\n");
+	feLog("   pressure norm : %15le%15le%15le\n", normL1, pctn, m_augtol);
+
+	// check convergence
+	bool bconv = true;
+	if (pctn >= m_augtol) bconv = false;
+	if (m_naugmin > naug) bconv = false;
+	if ((m_naugmax > 0) && (m_naugmax <= naug)) bconv = true;
+
+	// do the augmentation only if we have not yet converged
+	if (bconv == false)
+	{
+		for (int iel=0; iel<NE; ++iel)
+		{
+			// get the solid element
+			FESolidElement& el = m_Elem[iel];
+			int nint = el.GaussPoints();
+
+			// loop over the integration points and calculate
+			for (int n=0; n<nint; ++n)
+			{
+				FEMaterialPoint& mp = *el.GetMaterialPoint(n);
+				FEElasticMaterialPoint& pt = *(mp.ExtractData<FEElasticMaterialPoint>());
+
+				//Question: What is hi used for?
+				//double hi = pmi->h(pt.m_J, pt.m_J_star);
+				//pt.m_Lk += k*pmi->h(pt.m_J, pt.m_J_star);
+				//Question: Why reassign m_p here, doesn't it get overwritten?
+				//pt.m_p = pt.m_Lk*pmi->hp(pt.m_J, pt.m_J_star) + k*log(pt.m_J/pt.m_J_star)/pt.m_J;
+
+				pt.m_Lk += pmi->UJ(pt.m_J, pt.m_J_star);
+			}
+		}
+	}
+
+	return bconv;
+}
+
 
 //=================================================================================================
 
 BEGIN_FECORE_CLASS(FEStandardElasticSolidDomain, FEElasticSolidDomain)
 	ADD_PARAMETER(m_elemType, "elem_type", FE_PARAM_ATTRIBUTE, "$(solid_element)\0");
+	ADD_PARAMETER(m_blaugon, "laugon");
+	ADD_PARAMETER(m_augtol , "atol");
+	ADD_PARAMETER(m_naugmin, "minaug");
+	ADD_PARAMETER(m_naugmax, "maxaug");
 END_FECORE_CLASS();
 
 FEStandardElasticSolidDomain::FEStandardElasticSolidDomain(FEModel* fem) : FEElasticSolidDomain(fem)

FEBioMech/FEElasticSolidDomain.h

@@ -47,6 +47,11 @@ class FEBIOMECH_API FEElasticSolidDomain : public FESolidDomain, public FEElasti
 	//! activate
 	void Activate() override;
 
+	//! augmentation
+	bool Augment(int naug) override;
+
+	bool DoAugmentations() const;
+
 	//! initialize elements
 	void PreSolveUpdate(const FETimeInfo& timeInfo) override;
 
@@ -125,6 +130,11 @@ class FEBIOMECH_API FEElasticSolidDomain : public FESolidDomain, public FEElasti
 	bool	m_secant_stress;	//!< use secant approximation to stress
 	bool	m_secant_tangent;   //!< flag for using secant tangent
 
+	bool	m_blaugon;		//!< augmented lagrangian flag
+	double	m_augtol;		//!< augmented lagrangian tolerance
+	int		m_naugmin;		//!< minimum number of augmentations
+	int		m_naugmax;		//!< max number of augmentations
+
 protected:
 	FEDofList	m_dofU;		// displacement dofs
 	FEDofList	m_dofR;		// rigid rotation rofs

FEBioMech/FESolidSolver2.cpp

@@ -960,6 +960,9 @@ void FESolidSolver2::PrepStep()
 
 		FE3FieldElasticShellDomain* dom3fs = dynamic_cast<FE3FieldElasticShellDomain*>(dom);
 		if (dom3fs && dom3fs->DoAugmentations()) m_baugment = true;
+
+		FEElasticSolidDomain* doms = dynamic_cast<FEElasticSolidDomain*>(dom);
+		if (doms && doms->DoAugmentations()) m_baugment = true;
 	}
 
 	// see if we have to do nonlinear constraint augmentations