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;
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();
void GetActiveTensionAndTensionDerivs(double currentFibreStretch, unsigned currentQuadPointGlobalIndex, bool assembleJacobian, double &rActiveTension, double &rDerivActiveTensionWrtLambda, double &rDerivActiveTensionWrtDLambdaDt)