FELinearSystem now takes FEModel as constructor parameter instead of …

…FESolver.
febiosoftware · Oct 15, 2024 · 622ca03 · 622ca03
1 parent 1306d47
commit 622ca03
Show file tree

Hide file tree

Showing 27 changed files with 49 additions and 49 deletions.
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/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/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

diff --git a/FEBioTest/FEContactDiagnostic.cpp b/FEBioTest/FEContactDiagnostic.cpp
@@ -101,7 +101,7 @@ bool FEContactDiagnostic::Run()
 	int neq = solver.m_neq;
 	vector<double> Fd(neq, 0.0);
 	vector<double> ui(neq, 0.0);
-	FELinearSystem LS(&solver, K, Fd, ui, true);
+	FELinearSystem LS(&fem, K, Fd, ui, true);
 
 	// build the stiffness matrix
 	K0.Zero();

diff --git a/FEBioTest/FEContactDiagnosticBiphasic.cpp b/FEBioTest/FEContactDiagnosticBiphasic.cpp
@@ -155,7 +155,7 @@ bool FEContactDiagnosticBiphasic::Run()
 	int neq = solver.m_neq;
 	vector<double> Fd(neq, 0.0);
 	vector<double> ui(neq, 0.0);
-	FELinearSystem LS(&solver, K, Fd, ui, true);
+	FELinearSystem LS(&fem, K, Fd, ui, true);
 
     // build the stiffness matrix
     K.Zero();

diff --git a/FECore/FELinearSolver.cpp b/FECore/FELinearSolver.cpp
@@ -251,7 +251,7 @@ bool FELinearSolver::ReformStiffness()
 	{
 		TRACK_TIME(TimerID::Timer_Stiffness);
 
-		FELinearSystem K(this, *m_pK, m_R, m_u, (m_msymm == REAL_SYMMETRIC));
+		FELinearSystem K(GetFEModel(), *m_pK, m_R, m_u, (m_msymm == REAL_SYMMETRIC));
 		if (!StiffnessMatrix(K)) return false;
 
 		// do call back

diff --git a/FECore/FELinearSystem.cpp b/FECore/FELinearSystem.cpp
@@ -32,7 +32,7 @@ SOFTWARE.*/
 #include "FEModel.h"
 
 //-----------------------------------------------------------------------------
-FELinearSystem::FELinearSystem(FESolver* solver, FEGlobalMatrix& K, vector<double>& F, vector<double>& u, bool bsymm) : m_K(K), m_F(F), m_u(u), m_solver(solver)
+FELinearSystem::FELinearSystem(FEModel* fem, FEGlobalMatrix& K, vector<double>& F, vector<double>& u, bool bsymm) : m_K(K), m_F(F), m_u(u), m_fem(fem)
 {
 	m_bsymm = bsymm;
 }
@@ -50,14 +50,6 @@ bool FELinearSystem::IsSymmetric() const
 	return m_bsymm;
 }
 
-//-----------------------------------------------------------------------------
-// Get the solver that is using this linear system
-FESolver* FELinearSystem::GetSolver()
-{
-	return m_solver;
-}
-
-//-----------------------------------------------------------------------------
 //! assemble global stiffness matrix
 void FELinearSystem::Assemble(const FEElementMatrix& ke)
 {
@@ -98,16 +90,16 @@ void FELinearSystem::Assemble(const FEElementMatrix& ke)
 		}
 	}
 
-#pragma omp critical (LC_assemble)
-	{
-	FEModel* fem = m_solver->GetFEModel();
-	FELinearConstraintManager& LCM = fem->GetLinearConstraintManager();
-	if (LCM.LinearConstraints())
+	if (m_fem) 
 	{
-		const vector<int>& en = ke.Nodes();
-		LCM.AssembleStiffness(m_K, m_F, m_u, en, lmi, lmj, ke);
+		FELinearConstraintManager& LCM = m_fem->GetLinearConstraintManager();
+		if (LCM.LinearConstraints())
+		{
+			const vector<int>& en = ke.Nodes();
+			#pragma omp critical (LC_assemble)
+			LCM.AssembleStiffness(m_K, m_F, m_u, en, lmi, lmj, ke);
+		}
 	}
-	} // omp critical
 }
 
 //-----------------------------------------------------------------------------

diff --git a/FECore/FELinearSystem.h b/FECore/FELinearSystem.h
@@ -31,7 +31,7 @@ SOFTWARE.*/
 #include "matrix.h"
 #include <vector>
 
-class FESolver;
+class FEModel;
 
 //-----------------------------------------------------------------------------
 // Experimental class to see if all the assembly operations can be moved to a class
@@ -44,17 +44,14 @@ class FECORE_API FELinearSystem
 	// and a vector F which contains the assembled contribution of the prescribed 
 	// degrees of freedom. The F vector must be added to the "force" vector. The u 
 	// vector contains the nodal values of the prescribed degrees of freedom.
-	FELinearSystem(FESolver* solver, FEGlobalMatrix& K, std::vector<double>& F, std::vector<double>& u, bool bsymm);
+	FELinearSystem(FEModel* fem, FEGlobalMatrix& K, std::vector<double>& F, std::vector<double>& u, bool bsymm);
 
 	// virtual destructor
 	virtual ~FELinearSystem();
 
 	// get symmetry flag
 	bool IsSymmetric() const;
 
-	// Get the solver that is using this linear system
-	FESolver* GetSolver();
-
 public:
 	// Assembly routine
 	// This assembles the element stiffness matrix ke into the global matrix.
@@ -70,7 +67,7 @@ class FECORE_API FELinearSystem
 
 protected:
 	bool					m_bsymm;	//!< symmetry flag
-	FESolver*				m_solver;
+	FEModel*				m_fem;
 	FEGlobalMatrix&			m_K;	//!< The global stiffness matrix
 	std::vector<double>&	m_F;	//!< Contributions from prescribed degrees of freedom
 	std::vector<double>&	m_u;	//!< the array with prescribed values

diff --git a/FECore/FEModelComponent.cpp b/FECore/FEModelComponent.cpp
@@ -29,6 +29,7 @@ SOFTWARE.*/
 #include "stdafx.h"
 #include "FEModelComponent.h"
 #include "FEModel.h"
+#include "FEAnalysis.h"
 #include "DumpStream.h"
 #include <string.h>
 
@@ -92,6 +93,13 @@ const FETimeInfo& FEModelComponent::GetTimeInfo() const
 	return GetFEModel()->GetTime();
 }
 
+FESolver* FEModelComponent::GetSolver()
+{
+	FEAnalysis* step = GetFEModel()->GetCurrentStep();
+	if (step == nullptr) return nullptr;
+	return step->GetFESolver();
+}
+
 //-----------------------------------------------------------------------------
 void FEModelComponent::AttachLoadController(const char* szparam, int lc)
 {

diff --git a/FECore/FEModelComponent.h b/FECore/FEModelComponent.h
@@ -36,6 +36,7 @@ SOFTWARE.*/
 //! to the constructor.
 class FENodeSet;
 class FEMesh;
+class FESolver;
 
 //-----------------------------------------------------------------------------
 //! This class serves as a base class for many of the FECore classes. It defines
@@ -64,6 +65,7 @@ class FECORE_API FEModelComponent : public FECoreBase
 	int GetDOFIndex(const char* szvar, int n) const;
 	int GetDOFIndex(const char* szdof) const;
 	const FETimeInfo& GetTimeInfo() const;
+	FESolver* GetSolver();
 	void AttachLoadController(const char* szparam, int lc);
 	void AttachLoadController(void* pd, int lc);
 

diff --git a/FECore/FENewtonSolver.cpp b/FECore/FENewtonSolver.cpp
@@ -1123,7 +1123,7 @@ void FENewtonSolver::UpdateModel()
 bool FENewtonSolver::StiffnessMatrix()
 {
 	// setup the linear system
-	FELinearSystem LS(this, *m_pK, m_Fd, m_ui, (m_msymm == REAL_SYMMETRIC));
+	FELinearSystem LS(GetFEModel(), *m_pK, m_Fd, m_ui, (m_msymm == REAL_SYMMETRIC));
 
 	// build the stiffness matrix
 	return StiffnessMatrix(LS);