Merge branch 'develop' of github.com:febiosoftware/FEBioHeat into develop

Author: michaelrossherron

FEBioHeat/FEHeatSolidDomain.cpp

@@ -1,11 +1,11 @@
 #include "FEHeatSolidDomain.h"
-#include "FECore/FEModel.h"
+#include <FECore/FEModel.h>
 #include <FECore/Integrate.h>
 using namespace std::placeholders;
 
 //-----------------------------------------------------------------------------
 //! constructor
-FEHeatSolidDomain::FEHeatSolidDomain(FEModel* pfem) : FESolidDomain(pfem), FEHeatDomain(pfem), m_dof(pfem)
+FEHeatSolidDomain::FEHeatSolidDomain(FEModel* pfem) : FESolidDomain(pfem), m_dof(pfem)
 {
 	m_pMat = 0;
 
@@ -30,49 +30,10 @@ void FEHeatSolidDomain::ConductionMatrix(FELinearSystem& ls)
 
 //-----------------------------------------------------------------------------
 // Calculate the capacitance matrix
-void FEHeatSolidDomain::CapacitanceMatrix(FELinearSystem& ls, double dt)
+void FEHeatSolidDomain::CapacitanceMatrix(FELinearSystem& ls)
 {
-	vector<int> lm;
-	vector<double> fe;
-
-	FEDofList& dofs = m_dof;
-	assert(dofs.Size()==1);
-	int dofT = dofs[0];
-
-	// get the mesh
-	FEMesh& mesh = *GetMesh();
-
-	// loop over all elements in domain
-	for (int i=0; i<(int) m_Elem.size(); ++i)
-	{
-		FESolidElement& el = m_Elem[i];
-		int ne = el.Nodes();
-
-		// element capacitance matrix
-		FEElementMatrix kc(ne, ne);
-		ElementCapacitance(el, kc, dt);
-
-		// set up the LM matrix
-		UnpackLM(el, lm);
-		kc.SetIndices(lm);
-		kc.SetNodes(el.m_node);
-
-		// assemble into global matrix
-		ls.Assemble(kc);
-
-		// we also need to assemble this in the right-hand side
-		fe.resize(ne, 0.0);
-		for (int i = 0; i<ne; ++i)
-		{
-			fe[i] = mesh.Node(el.m_node[i]).get(dofT);
-		}
-
-		// multiply with me
-		fe = kc*fe;
-
-		// assemble this vector to the right-hand side
-		ls.AssembleRHS(lm, fe);
-	}
+	auto fp = bind(&FEHeatSolidDomain::ElementCapacitance, this, _1, _2);
+	AssembleSolidDomain(*this, ls, fp);
 }
 
 //-----------------------------------------------------------------------------
@@ -92,6 +53,40 @@ void FEHeatSolidDomain::HeatSource(FEGlobalVector& R, FEHeatSource& hs)
 	}
 }
 
+void FEHeatSolidDomain::CapacitanceLoad(FEGlobalVector& R)
+{
+	auto fp = bind(&FEHeatSolidDomain::ElementCapacitanceLoad, this, _1, _2);
+	AssembleSolidDomain(*this, R, fp);
+}
+
+void FEHeatSolidDomain::ElementCapacitanceLoad(FESolidElement& el, vector<double>& fe)
+{
+	FETimeInfo& ti = GetFEModel()->GetTime();
+	double dt = ti.timeIncrement;
+
+	double c = m_pMat->Capacitance();
+	double d = m_pMat->Density();
+	double alpha = c * d / dt;
+
+	zero(fe);
+	double* w = el.GaussWeights();
+	int ne = el.Nodes();
+	int ni = el.GaussPoints();
+	for (int n = 0; n < ni; ++n)
+	{
+		FEMaterialPoint& mp = *el.GetMaterialPoint(n);
+		FEHeatMaterialPoint* pt = (mp.ExtractData<FEHeatMaterialPoint>());
+
+		double Ti = pt->m_T;
+		double* H = el.H(n);
+		double J = detJt(el, n);
+		for (int i = 0; i < ne; ++i)
+		{
+			fe[i] += alpha * H[i] * Ti * J * w[n];
+		}
+	}
+}
+
 //-----------------------------------------------------------------------------
 //! calculate element contribution to heat source term
 void FEHeatSolidDomain::ElementHeatSource(FEHeatSource& hs, FESolidElement& el, vector<double>& fe)
@@ -129,8 +124,11 @@ void FEHeatSolidDomain::ElementConduction(FESolidElement& el, matrix& ke)
 }
 
 //-----------------------------------------------------------------------------
-void FEHeatSolidDomain::ElementCapacitance(FESolidElement &el, matrix &ke, double dt)
+void FEHeatSolidDomain::ElementCapacitance(FESolidElement &el, matrix &ke)
 {
+	FETimeInfo& ti = GetFEModel()->GetTime();
+	double dt = ti.timeIncrement;
+
 	// zero stiffness matrix
 	ke.zero();
 
@@ -168,6 +166,8 @@ void FEHeatSolidDomain::Update(const FETimeInfo& tp)
 			FEHeatMaterialPoint* pt = (mp.ExtractData<FEHeatMaterialPoint>());
 			assert(pt);
 
+			pt->m_T = el.Evaluate(T, n);
+
 			vec3d gradT = gradient(el, T, n);
 			mat3ds D = m_pMat->Conductivity(mp);
 			pt->m_q = -(D*gradT);

FEBioHeat/FEHeatSolidDomain.h

@@ -11,16 +11,12 @@ class FEModel;
 class FEHeatDomain
 {
 public:
-	FEHeatDomain(FEModel* pfem) : m_pfem(pfem) {}
+	FEHeatDomain() {}
 	virtual ~FEHeatDomain(){}
 	virtual void ConductionMatrix (FELinearSystem& ls) = 0;
-	virtual void CapacitanceMatrix(FELinearSystem& ls, double dt) = 0;
+	virtual void CapacitanceMatrix(FELinearSystem& ls) = 0;
 	virtual void HeatSource(FEGlobalVector& R, FEHeatSource& hs) = 0;
-
-	FEModel* GetFEModel() { return m_pfem; }
-
-protected:
-	FEModel* m_pfem;
+	virtual void CapacitanceLoad(FEGlobalVector& R) = 0;
 };
 
 //-----------------------------------------------------------------------------
@@ -41,29 +37,34 @@ class FEHeatSolidDomain : public FESolidDomain, public FEHeatDomain
 
 	//! Update state data
 	void Update(const FETimeInfo& tp) override;
+
+	const FEDofList& GetDOFList() const { return m_dof; }
 
 public: // overloaded from FEHeatDomain
 
 	//! Calculate the conduction stiffness 
 	void ConductionMatrix(FELinearSystem& ls) override;
 
 	//! Calculate capacitance stiffness matrix
-	void CapacitanceMatrix(FELinearSystem& ls, double dt) override;
+	void CapacitanceMatrix(FELinearSystem& ls) override;
 
 	// evaluate heat source
 	void HeatSource(FEGlobalVector& R, FEHeatSource& hs) override;
+
+	// evaluate capacitance load
+	void CapacitanceLoad(FEGlobalVector& R) override;
 
 protected:
 	//! calculate the conductive element stiffness matrix
 	void ElementConduction(FESolidElement& el, matrix& ke);
 
 	//! calculate the capacitance element stiffness matrix
-	void ElementCapacitance(FESolidElement& el, matrix& ke, double dt);
+	void ElementCapacitance(FESolidElement& el, matrix& ke);
 
 	//! calculate element contribution to heat source term
 	void ElementHeatSource(FEHeatSource& hs, FESolidElement& el, vector<double>& fe);
 
-	const FEDofList& GetDOFList() const { return m_dof; }
+	void ElementCapacitanceLoad(FESolidElement& el, vector<double>& fe);
 
 protected:
 	FEHeatTransferMaterial*	m_pMat;

FEBioHeat/FEHeatSolver.cpp

@@ -72,6 +72,25 @@ bool FEHeatSolver::Init()
 	return true;
 }
 
+void FEHeatSolver::ForceVector(FEGlobalVector& R)
+{
+	FELinearSolver::ForceVector(R);
+
+	FEModel& fem = *GetFEModel();
+	FEAnalysis* pstep = fem.GetCurrentStep();
+	bool btransient = (pstep->m_nanalysis == FEHeatAnalysis::TRANSIENT);
+
+	if (btransient)
+	{
+		// for a dynamic analysis add the capacitance load
+		for (int i = 0; i < pstep->Domains(); ++i)
+		{
+			FEHeatDomain& bd = dynamic_cast<FEHeatDomain&>(*pstep->Domain(i));
+			bd.CapacitanceLoad(R);
+		}
+	}
+}
+
 //-----------------------------------------------------------------------------
 //! Calculate the global stiffness matrix. This function simply calls 
 //! HeatStiffnessMatrix() for each domain which will calculate the
@@ -99,7 +118,7 @@ bool FEHeatSolver::StiffnessMatrix(FELinearSystem& LS)
 		// for a dynamic analysis add the capacitance matrix
 		if (btransient)
 		{
-			bd.CapacitanceMatrix(LS, tp.timeIncrement);
+			bd.CapacitanceMatrix(LS);
 		}
 	}
 

FEBioHeat/FEHeatSolver.h

@@ -21,6 +21,8 @@ class FEHeatSolver : public FELinearSolver
 
 protected: // from FELinearSolver
 
+	void ForceVector(FEGlobalVector& R) override;
+
 	//! calculate the stiffness matrix
 	bool StiffnessMatrix(FELinearSystem& LS) override; 
 

FEBioHeat/FEIsotropicFourier.cpp

@@ -1,14 +1,18 @@
 #include "FEIsotropicFourier.h"
 
-//-----------------------------------------------------------------------------
-// define the parameter list
 BEGIN_FECORE_CLASS(FEIsotropicFourier, FEMaterial)
 	ADD_PARAMETER(m_rho, FE_RANGE_GREATER(0.0), "density")->setUnits(UNIT_DENSITY);
 	ADD_PARAMETER(m_k  , FE_RANGE_GREATER(0.0), "k")->setUnits("W/L.T");
 	ADD_PARAMETER(m_c  , FE_RANGE_GREATER(0.0), "c")->setUnits("E/M.T");
 END_FECORE_CLASS();
 
-//-----------------------------------------------------------------------------
+FEIsotropicFourier::FEIsotropicFourier(FEModel* pfem) : FEHeatTransferMaterial(pfem) 
+{
+	m_rho = 1;
+	m_k = 0;
+	m_c = 0;
+}
+
 mat3ds FEIsotropicFourier::Conductivity(FEMaterialPoint& mp)
 {
 	return mat3dd(m_k);

FEBioHeat/FEIsotropicFourier.h

@@ -1,12 +1,11 @@
 #pragma once
 #include "FEHeatTransferMaterial.h"
 
-//-----------------------------------------------------------------------------
-// Isotropic Fourer heat-transfer material
+// Isotropic Fourier heat-transfer material
 class FEIsotropicFourier : public FEHeatTransferMaterial
 {
 public:
-	FEIsotropicFourier(FEModel* pfem) : FEHeatTransferMaterial(pfem) {}
+	FEIsotropicFourier(FEModel* pfem);
 
 public:
 	double	m_k;	//!< heat conductivity

README.md

@@ -0,0 +1,16 @@
+# FEBioHeat
+The FEBioHeat plugin implements a solver for the heat transfer equation that can be used with [febio](https://febio.org/).
+It currently supports steady-state and transient heat transfer analyses and assumes Fourier's law for calculating the heat flux.
+
+Supported boundary conditions are: 
+* Zero and prescribed temperature
+* Prescribed heat flux
+* Convective heat flux
+* Gap heat flux (enforces continuous heat flux across an interface)
+
+Supported body loads are:
+* Constant heat source
+* Bio-heat
+
+For transient problems, the initial temperature can be prescribed. 
+