Merge branch 'develop' into erosion_changes

febiosoftware · Oct 24, 2024 · 29c30d8 · 29c30d8
2 parents d4174dc + 8fb2d34
commit 29c30d8
Show file tree

Hide file tree

Showing 59 changed files with 1,534 additions and 109 deletions.
diff --git a/Documentation/Doxygen/register.dox b/Documentation/Doxygen/register.dox
@@ -17,13 +17,13 @@ The following code excerpt presents an example of PluginInitialize.
 \code
 FECORE_EXPORT void PluginInitialize(FECoreKernel& fecore)
 {
-	\\ Set the kernel to match FEBio's kernel
+	// Set the kernel to match FEBio's kernel
 	FECoreKernel::SetInstance(&fecore);
 
-	\\ set the active module to which to add this new feature
+	// set the active module to which to add this new feature
 	fecore.SetActiveModule("solid");
 
-	\\ Register the classes in this plugin
+	// Register the classes in this plugin
 	REGISTER_FECORE_CLASS(MyMaterial, "My material");
 }
 \endcode

diff --git a/FEBio/FEBioApp.cpp b/FEBio/FEBioApp.cpp
@@ -175,16 +175,14 @@ int FEBioApp::RunModel()
 	{
 		// read the input file
 		if (fem.Input(m_ops.szfile) == false) nret = 1;
-		else
-		{
-			// apply configuration overrides
-			ApplyConfig(fem);
-		}
 	}
 
 	// solve the model with the task and control file
 	if (nret == 0)
 	{
+		// apply configuration overrides
+		ApplyConfig(fem);
+
 		bool bret = febio::SolveModel(fem, m_ops.sztask, m_ops.szctrl);
 
 		nret = (bret ? 0 : 1);

diff --git a/FEBioFluid/FEFluidFSISolver.cpp b/FEBioFluid/FEFluidFSISolver.cpp
@@ -1113,7 +1113,7 @@ bool FEFluidFSISolver::StiffnessMatrix()
     // get the mesh
     FEMesh& mesh = fem.GetMesh();
 
-	FESolidLinearSystem LS(this, &m_rigidSolver, *m_pK, m_Fd, m_ui, (m_msymm == REAL_SYMMETRIC), m_alphaf, m_nreq);
+	FESolidLinearSystem LS(&fem, &m_rigidSolver, *m_pK, m_Fd, m_ui, (m_msymm == REAL_SYMMETRIC), m_alphaf, m_nreq);
 
     // calculate the stiffness matrix for each domain
     for (int i=0; i<mesh.Domains(); ++i)

diff --git a/FEBioFluid/FEMultiphasicFSISolver.cpp b/FEBioFluid/FEMultiphasicFSISolver.cpp
@@ -1249,7 +1249,7 @@ bool FEMultiphasicFSISolver::StiffnessMatrix()
     // get the mesh
     FEMesh& mesh = fem.GetMesh();
 
-    FESolidLinearSystem LS(this, &m_rigidSolver, *m_pK, m_Fd, m_ui, (m_msymm == REAL_SYMMETRIC), m_alphaf, m_nreq);
+    FESolidLinearSystem LS(&fem, &m_rigidSolver, *m_pK, m_Fd, m_ui, (m_msymm == REAL_SYMMETRIC), m_alphaf, m_nreq);
 
     // calculate the stiffness matrix for each domain
     for (int i=0; i<mesh.Domains(); ++i)

diff --git a/FEBioLib/FEBioModel.cpp b/FEBioLib/FEBioModel.cpp
@@ -190,7 +190,6 @@ ModelStats FEBioModel::GetModelStats() const
 
 ModelStats FEBioModel::GetStepStats(size_t n) const
 {
-	assert((n >= 0) && (n < Steps()));
 	return m_stepStats[n];
 }
 
@@ -1336,6 +1335,28 @@ void FEBioModel::Serialize(DumpStream& ar)
 
 		// --- Save IO Data
 		SerializeIOData(ar);
+
+		if (ar.IsSaving())
+		{
+			int n = (int)m_stepStats.size();
+			ar << n;
+			for (ModelStats& s : m_stepStats)
+			{
+				ar << s.ntimeSteps << s.ntotalIters << s.ntotalReforms << s.ntotalRHS;
+			}
+		}
+		else
+		{
+			m_stepStats.clear();
+			int n = 0;
+			ar >> n;
+			for (int i = 0; i < n; ++i)
+			{
+				ModelStats s;
+				ar >> s.ntimeSteps >> s.ntotalIters >> s.ntotalReforms >> s.ntotalRHS;
+				m_stepStats.push_back(s);
+			}
+		}
 	}
 }
 
@@ -1546,7 +1567,7 @@ bool FEBioModel::Init()
 	}
 
 	m_report.clear();
-	m_stepStats.resize(Steps());
+	m_stepStats.clear();
 
 	FEBioPlotFile* pplt = nullptr;
 	m_lastUpdate = -1;
@@ -1640,6 +1661,9 @@ bool FEBioModel::InitLogFile()
 
 bool FEBioModel::Solve()
 {
+	// The total time is usually started on calling Input,
+	// however, in a restart Input is not called, so we start it here.
+	if (!m_TotalTime.isRunning()) m_TotalTime.start();
 	bool b = FEModel::Solve();
 	m_TotalTime.stop();
 	return b;
@@ -1693,13 +1717,7 @@ bool FEBioModel::Reset()
 	m_modelStats.ntotalIters = 0;
 	m_modelStats.ntotalRHS = 0;
 	m_modelStats.ntotalReforms = 0;
-	for (auto& s : m_stepStats)
-	{
-		s.ntimeSteps = 0;
-		s.ntotalIters = 0;
-		s.ntotalRHS = 0;
-		s.ntotalReforms = 0;
-	}
+	m_stepStats.clear();
 
 	// do the callback
 	DoCallback(CB_INIT);
@@ -1860,7 +1878,7 @@ void FEBioModel::on_cb_stepSolved()
 	stats.ntotalIters   = step->m_ntotiter;
 	stats.ntotalRHS     = step->m_ntotrhs;
 	stats.ntotalReforms = step->m_ntotref;
-	m_stepStats[GetCurrentStepIndex()] = stats;
+	m_stepStats.push_back(stats);
 	m_modelStats.ntimeSteps    += stats.ntimeSteps;
 	m_modelStats.ntotalIters   += stats.ntotalIters;
 	m_modelStats.ntotalRHS     += stats.ntotalRHS;

diff --git a/FEBioLib/FEBioRestart.cpp b/FEBioLib/FEBioRestart.cpp
@@ -30,6 +30,7 @@ SOFTWARE.*/
 #include <FEBioXML/FERestartImport.h>
 #include <FECore/DumpFile.h>
 #include <FECore/FEAnalysis.h>
+#include "FEBioModelBuilder.h"
 
 //-----------------------------------------------------------------------------
 FEBioRestart::FEBioRestart(FEModel* pfem) : FECoreTask(pfem) {}
@@ -75,6 +76,7 @@ bool FEBioRestart::Init(const char* szfile)
 
 		// the file is assumed to be a xml-text input file
 		FERestartImport file;
+		file.SetModelBuilder(new FEBioModelBuilder(fem));
 		if (file.Load(fem, szfile) == false)
 		{
 			char szerr[256];

diff --git a/FEBioMech/FEDonnanEquilibrium.cpp b/FEBioMech/FEDonnanEquilibrium.cpp
@@ -93,6 +93,12 @@ bool FEDonnanEquilibrium::Init()
 	return FEElasticMaterial::Init();
 }
 
+void FEDonnanEquilibrium::Serialize(DumpStream& ar)
+{
+	FEElasticMaterial::Serialize(ar);
+	ar & m_Rgas & m_Tabs & m_binit & m_bnew & m_phiwr;
+}
+
 //-----------------------------------------------------------------------------
 mat3ds FEDonnanEquilibrium::Stress(FEMaterialPoint& mp)
 {

diff --git a/FEBioMech/FEDonnanEquilibrium.h b/FEBioMech/FEDonnanEquilibrium.h
@@ -64,6 +64,8 @@ class FEDonnanEquilibrium : public FEElasticMaterial
 	//! Initialization routine
 	bool Init() override;
 
+	void Serialize(DumpStream& ar) override;
+
 	//! Returns the Cauchy stress
 	mat3ds Stress(FEMaterialPoint& mp) override;
 

diff --git a/FEBioMech/FEFacet2FacetSliding.cpp b/FEBioMech/FEFacet2FacetSliding.cpp
@@ -780,7 +780,7 @@ void FEFacet2FacetSliding::StiffnessMatrix(FELinearSystem& LS, const FETimeInfo&
 	FEMesh* pm = m_ss.GetMesh();
 
 	// see how many reformations we've had to do so far
-	int nref = LS.GetSolver()->m_nref;
+	int nref = GetSolver()->m_nref;
 
 	// get the "size" of the model
 	// We need this to scale the insertion distance

diff --git a/FEBioMech/FEReactiveViscoelastic.cpp b/FEBioMech/FEReactiveViscoelastic.cpp
@@ -762,3 +762,9 @@ double FEReactiveViscoelasticMaterial::Damage(FEMaterialPoint& mp)
     }
     return D;
 }
+
+void FEReactiveViscoelasticMaterial::Serialize(DumpStream& ar)
+{
+	FEElasticMaterial::Serialize(ar);
+	ar & m_nmax;
+}
diff --git a/FEBioMech/FEReactiveViscoelastic.h b/FEBioMech/FEReactiveViscoelastic.h
@@ -61,6 +61,8 @@ class FEReactiveViscoelasticMaterial :	public FEElasticMaterial
     //! data initialization
     bool Init() override;
 
+	void Serialize(DumpStream& ar) override;
+
 	//! stress function
 	mat3ds Stress(FEMaterialPoint& pt) override;
     mat3ds StressStrongBonds(FEMaterialPoint& pt);
@@ -105,7 +107,7 @@ class FEReactiveViscoelasticMaterial :	public FEElasticMaterial
 
     //! evaluate scalar strain measure (same type as trigger strain for bond breaking)
     double ScalarStrain(FEMaterialPoint& mp);
-    
+
 private:
 	FEElasticMaterial*	m_pBase;	//!< pointer to elastic solid material for strong bonds
 	FEElasticMaterial*	m_pBond;	//!< pointer to elastic solid material for reactive bonds

diff --git a/FEBioMech/FESlidingElasticInterface.cpp b/FEBioMech/FESlidingElasticInterface.cpp
@@ -1094,7 +1094,7 @@ void FESlidingElasticInterface::LoadVector(FEGlobalVector& R, const FETimeInfo&
 void FESlidingElasticInterface::StiffnessMatrix(FELinearSystem& LS, const FETimeInfo& tp)
 {
     // see how many reformations we've had to do so far
-    int nref = LS.GetSolver()->m_nref;
+    int nref = GetSolver()->m_nref;
 
     const int MN = FEElement::MAX_NODES;
 

diff --git a/FEBioMech/FESolidLinearSystem.cpp b/FEBioMech/FESolidLinearSystem.cpp
@@ -31,8 +31,9 @@ SOFTWARE.*/
 #include <FECore/FELinearConstraintManager.h>
 #include <FECore/FEModel.h>
 
-FESolidLinearSystem::FESolidLinearSystem(FESolver* solver, FERigidSolver* rigidSolver, FEGlobalMatrix& K, std::vector<double>& F, std::vector<double>& u, bool bsymm, double alpha, int nreq) : FELinearSystem(solver, K, F, u, bsymm)
+FESolidLinearSystem::FESolidLinearSystem(FEModel* fem, FERigidSolver* rigidSolver, FEGlobalMatrix& K, std::vector<double>& F, std::vector<double>& u, bool bsymm, double alpha, int nreq) : FELinearSystem(fem, K, F, u, bsymm)
 {
+	m_fem = fem;
 	m_rigidSolver = rigidSolver;
 	m_alpha = alpha;
 	m_nreq = nreq;
@@ -73,8 +74,7 @@ void FESolidLinearSystem::Assemble(const FEElementMatrix& ke)
 		vector<double>& ui = m_u;
 
 		// adjust for linear constraints
-		FEModel* fem = m_solver->GetFEModel();
-		FELinearConstraintManager& LCM = fem->GetLinearConstraintManager();
+		FELinearConstraintManager& LCM = m_fem->GetLinearConstraintManager();
 		if (LCM.LinearConstraints() > 0)
 		{
 			#pragma omp critical (LCM_assemble)

diff --git a/FEBioMech/FESolidLinearSystem.h b/FEBioMech/FESolidLinearSystem.h
@@ -35,7 +35,7 @@ class FERigidSolver;
 class FEBIOMECH_API FESolidLinearSystem : public FELinearSystem
 {
 public:
-	FESolidLinearSystem(FESolver* solver, FERigidSolver* rigidSolver, FEGlobalMatrix& K, std::vector<double>& F, std::vector<double>& u, bool bsymm, double alpha, int nreq);
+	FESolidLinearSystem(FEModel* fem, FERigidSolver* rigidSolver, FEGlobalMatrix& K, std::vector<double>& F, std::vector<double>& u, bool bsymm, double alpha, int nreq);
 
 	// Assembly routine
 	// This assembles the element stiffness matrix ke into the global matrix.
@@ -46,6 +46,7 @@ class FEBIOMECH_API FESolidLinearSystem : public FELinearSystem
 	void StiffnessAssemblyScaleFactor(double a);
 
 private:
+	FEModel* fem;
 	FERigidSolver*	m_rigidSolver;
 	double			m_alpha;
 	int				m_nreq;

diff --git a/FEBioMech/FESolidSolver.cpp b/FEBioMech/FESolidSolver.cpp
@@ -527,7 +527,7 @@ bool FESolidSolver::StiffnessMatrix()
 	const FETimeInfo& tp = fem.GetTime();
 
 	// setup the linear syster
-	FESolidLinearSystem LS(this, &m_rigidSolver, *m_pK, m_Fd, m_ui, (m_msymm == REAL_SYMMETRIC), 1.0, m_nreq);
+	FESolidLinearSystem LS(&fem, &m_rigidSolver, *m_pK, m_Fd, m_ui, (m_msymm == REAL_SYMMETRIC), 1.0, m_nreq);
 
 	// get the mesh
 	FEMesh& mesh = fem.GetMesh();

diff --git a/FEBioMech/FESolidSolver2.cpp b/FEBioMech/FESolidSolver2.cpp
@@ -287,7 +287,7 @@ bool FESolidSolver2::InitAccelerations()
 
 		// setup the linear system
 		m_pK->Zero();
-		FESolidLinearSystem LS(this, &m_rigidSolver, *m_pK, m_Fd, m_ui, (m_msymm == REAL_SYMMETRIC), 1.0, m_nreq);
+		FESolidLinearSystem LS(&fem, &m_rigidSolver, *m_pK, m_Fd, m_ui, (m_msymm == REAL_SYMMETRIC), 1.0, m_nreq);
 
 		// build the global mass matrix
 		FEMesh& mesh = fem.GetMesh();
@@ -1324,7 +1324,7 @@ bool FESolidSolver2::StiffnessMatrix()
 	FEMesh& mesh = fem.GetMesh();
 
 	// setup the linear system
-	FESolidLinearSystem LS(this, &m_rigidSolver, *m_pK, m_Fd, m_ui, (m_msymm == REAL_SYMMETRIC), m_alpha, m_nreq);
+	FESolidLinearSystem LS(&fem, &m_rigidSolver, *m_pK, m_Fd, m_ui, (m_msymm == REAL_SYMMETRIC), m_alpha, m_nreq);
 
 	// calculate the stiffness matrix for each domain
 	for (int i=0; i<mesh.Domains(); ++i) 

diff --git a/FEBioMech/FETiedElasticInterface.cpp b/FEBioMech/FETiedElasticInterface.cpp
@@ -693,7 +693,7 @@ void FETiedElasticInterface::StiffnessMatrix(FELinearSystem& LS, const FETimeInf
     FEElementMatrix ke;
 
     // see how many reformations we've had to do so far
-    int nref = LS.GetSolver()->m_nref;
+    int nref = GetSolver()->m_nref;
 
     // set higher order stiffness mutliplier
     // NOTE: this algrotihm doesn't really need this

diff --git a/FEBioMix/FEBiphasicSoluteSolver.cpp b/FEBioMix/FEBiphasicSoluteSolver.cpp
@@ -487,7 +487,7 @@ bool FEBiphasicSoluteSolver::StiffnessMatrix()
 	FEMesh& mesh = fem.GetMesh();
 
 	// setup the linear system
-	FESolidLinearSystem LS(this, &m_rigidSolver, *m_pK, m_Fd, m_ui, (m_msymm == REAL_SYMMETRIC), m_alpha, m_nreq);
+	FESolidLinearSystem LS(&fem, &m_rigidSolver, *m_pK, m_Fd, m_ui, (m_msymm == REAL_SYMMETRIC), m_alpha, m_nreq);
 
 	// calculate the stiffness matrix for each domain
 	FEAnalysis* pstep = fem.GetCurrentStep();

diff --git a/FEBioMix/FEBiphasicSolver.cpp b/FEBioMix/FEBiphasicSolver.cpp
@@ -608,7 +608,7 @@ bool FEBiphasicSolver::StiffnessMatrix()
 	FEMesh& mesh = fem.GetMesh();
 
 	// setup the linear system of equations
-	FESolidLinearSystem LS(this, &m_rigidSolver, *m_pK, m_Fd, m_ui, (m_msymm == REAL_SYMMETRIC), m_alpha, m_nreq);
+	FESolidLinearSystem LS(&fem, &m_rigidSolver, *m_pK, m_Fd, m_ui, (m_msymm == REAL_SYMMETRIC), m_alpha, m_nreq);
 
 	// calculate the stiffness matrix for each domain
 	FEAnalysis* pstep = fem.GetCurrentStep();

diff --git a/FEBioMix/FEMultiphasicSolver.cpp b/FEBioMix/FEMultiphasicSolver.cpp
@@ -752,7 +752,7 @@ bool FEMultiphasicSolver::StiffnessMatrix()
 	// get the mesh
 	FEMesh& mesh = fem.GetMesh();
 
-	FESolidLinearSystem LS(this, &m_rigidSolver, *m_pK, m_Fd, m_ui, (m_msymm == REAL_SYMMETRIC), m_alpha, m_nreq);
+	FESolidLinearSystem LS(&fem, &m_rigidSolver, *m_pK, m_Fd, m_ui, (m_msymm == REAL_SYMMETRIC), m_alpha, m_nreq);
 
 	// calculate the stiffness matrix for each domain
 	FEAnalysis* pstep = fem.GetCurrentStep();

diff --git a/FEBioMix/FESlidingInterface2.cpp b/FEBioMix/FESlidingInterface2.cpp
@@ -1209,7 +1209,7 @@ void FESlidingInterface2::StiffnessMatrix(FELinearSystem& LS, const FETimeInfo&
     double psf = GetPenaltyScaleFactor();
 
 	// see how many reformations we've had to do so far
-	int nref = LS.GetSolver()->m_nref;
+	int nref = GetSolver()->m_nref;
 
 	// set higher order stiffness mutliplier
 	// NOTE: this algrotihm doesn't really need this

diff --git a/FEBioMix/FESlidingInterface3.cpp b/FEBioMix/FESlidingInterface3.cpp
@@ -1379,7 +1379,7 @@ void FESlidingInterface3::StiffnessMatrix(FELinearSystem& LS, const FETimeInfo&
     double psf = GetPenaltyScaleFactor();
 
 	// see how many reformations we've had to do so far
-	int nref = LS.GetSolver()->m_nref;
+	int nref = GetSolver()->m_nref;
 
 	// set higher order stiffness mutliplier
 	// NOTE: this algrotihm doesn't really need this

diff --git a/FEBioMix/FESlidingInterfaceBiphasic.cpp b/FEBioMix/FESlidingInterfaceBiphasic.cpp
@@ -1590,7 +1590,7 @@ void FESlidingInterfaceBiphasic::LoadVector(FEGlobalVector& R, const FETimeInfo&
 void FESlidingInterfaceBiphasic::StiffnessMatrix(FELinearSystem& LS, const FETimeInfo& tp)
 {
     // see how many reformations we've had to do so far
-    int nref = LS.GetSolver()->m_nref;
+    int nref = GetSolver()->m_nref;
 
     const int MN = FEElement::MAX_NODES;
 

diff --git a/FEBioMix/FESlidingInterfaceBiphasicMixed.cpp b/FEBioMix/FESlidingInterfaceBiphasicMixed.cpp
@@ -1688,7 +1688,7 @@ void FESlidingInterfaceBiphasicMixed::StiffnessMatrix(FESlidingSurfaceBiphasicMi
 	int degree_p = dofs.GetVariableInterpolationOrder(ss.m_varP);
 
     // see how many reformations we've had to do so far
-    int nref = LS.GetSolver()->m_nref;
+    int nref = GetSolver()->m_nref;
 
     const int MN = FEElement::MAX_NODES;
 

diff --git a/FEBioMix/FESlidingInterfaceMP.cpp b/FEBioMix/FESlidingInterfaceMP.cpp
@@ -2011,7 +2011,7 @@ void FESlidingInterfaceMP::StiffnessMatrix(FELinearSystem& LS, const FETimeInfo&
     double psf = GetPenaltyScaleFactor();
 
     // see how many reformations we've had to do so far
-    int nref = LS.GetSolver()->m_nref;
+    int nref = GetSolver()->m_nref;
 
     // set higher order stiffness mutliplier
     // NOTE: this algorithm doesn't really need this

diff --git a/FEBioMix/FETiedBiphasicInterface.cpp b/FEBioMix/FETiedBiphasicInterface.cpp
@@ -784,7 +784,7 @@ void FETiedBiphasicInterface::StiffnessMatrix(FELinearSystem& LS, const FETimeIn
 	FEMesh* pm = m_ss.GetMesh();
 
 	// see how many reformations we've had to do so far
-	int nref = LS.GetSolver()->m_nref;
+	int nref = GetSolver()->m_nref;
 
 	// set higher order stiffness mutliplier
 	// NOTE: this algrotihm doesn't really need this

diff --git a/FEBioMix/FETiedMultiphasicInterface.cpp b/FEBioMix/FETiedMultiphasicInterface.cpp
@@ -1073,7 +1073,7 @@ void FETiedMultiphasicInterface::StiffnessMatrix(FELinearSystem& LS, const FETim
     FEModel& fem = *GetFEModel();
 
     // see how many reformations we've had to do so far
-    int nref = LS.GetSolver()->m_nref;
+    int nref = GetSolver()->m_nref;
 
     // set higher order stiffness mutliplier
     // NOTE: this algorithm doesn't really need this