36 #include "ImplicitCardiacMechanicsSolver.hpp"
38 template<
class ELASTICITY_SOLVER,
unsigned DIM>
42 std::string outputDirectory)
50 template<
class ELASTICITY_SOLVER,
unsigned DIM>
56 template<
class ELASTICITY_SOLVER,
unsigned DIM>
60 assert(time < nextTime);
61 this->mCurrentTime = time;
62 this->mNextTime = nextTime;
63 this->mOdeTimestep = odeTimestep;
66 ELASTICITY_SOLVER::Solve();
71 if (this->GetNumNewtonIterations() > 0)
73 this->AssembleSystem(
true,
false);
78 for (std::map<unsigned,DataAtQuadraturePoint>::iterator iter = this->mQuadPointToDataAtQuadPointMap.begin();
79 iter != this->mQuadPointToDataAtQuadPointMap.end();
83 iter->second.StretchLastTimeStep = iter->second.Stretch;
88 template<
class ELASTICITY_SOLVER,
unsigned DIM>
90 unsigned currentQuadPointGlobalIndex,
91 bool assembleJacobian,
92 double& rActiveTension,
93 double& rDerivActiveTensionWrtLambda,
94 double& rDerivActiveTensionWrtDLambdaDt)
98 assert(this->mMapIterator->first==currentQuadPointGlobalIndex);
103 r_data_at_quad_point.
Stretch = currentFibreStretch;
106 double dlam_dt = (currentFibreStretch-r_data_at_quad_point.
StretchLastTimeStep)/(this->mNextTime-this->mCurrentTime);
115 p_contraction_model->
RunDoNotUpdate(this->mCurrentTime,this->mNextTime,this->mOdeTimestep);
122 if (assembleJacobian)
125 EXCEPTION(
"Failure in solving contraction models using current stretches for assembling Jacobian");
129 rActiveTension = DBL_MAX;
130 std::cout <<
"WARNING: could not solve contraction model with this stretch and stretch rate. "
131 <<
"Setting active tension to infinity (DBL_MAX) so that the residual(-norm) is also infinite\n" << std::flush;
138 if (assembleJacobian)
141 double h1 = std::max(1e-6, currentFibreStretch/100);
143 p_contraction_model->
RunDoNotUpdate(this->mCurrentTime,this->mNextTime,this->mOdeTimestep);
147 double h2 = std::max(1e-6, dlam_dt/100);
149 p_contraction_model->
RunDoNotUpdate(this->mCurrentTime,this->mNextTime,this->mOdeTimestep);
152 rDerivActiveTensionWrtLambda = (active_tension_at_lam_plus_h - rActiveTension)/h1;
153 rDerivActiveTensionWrtDLambdaDt = (active_tension_at_dlamdt_plus_h - rActiveTension)/h2;
157 this->mMapIterator++;
158 if (this->mMapIterator==this->mQuadPointToDataAtQuadPointMap.end())
160 this->mMapIterator = this->mQuadPointToDataAtQuadPointMap.begin();
ImplicitCardiacMechanicsSolver(QuadraticMesh< DIM > &rQuadMesh, ElectroMechanicsProblemDefinition< DIM > &rProblemDefinition, std::string outputDirectory)
virtual ~ImplicitCardiacMechanicsSolver()
virtual void UpdateStateVariables()=0
void GetActiveTensionAndTensionDerivs(double currentFibreStretch, unsigned currentQuadPointGlobalIndex, bool assembleJacobian, double &rActiveTension, double &rDerivActiveTensionWrtLambda, double &rDerivActiveTensionWrtDLambdaDt)
#define EXCEPTION(message)
double StretchLastTimeStep
virtual void SetStretchAndStretchRate(double stretch, double stretchRate)=0
AbstractContractionModel * ContractionModel
virtual void RunDoNotUpdate(double startTime, double endTime, double timeStep)=0
virtual double GetNextActiveTension()=0
void Solve(double time, double nextTime, double odeTimestep)