ImplicitCardiacMechanicsAssembler.cpp

00001 /*
00002 
00003 Copyright (C) University of Oxford, 2005-2010
00004 
00005 University of Oxford means the Chancellor, Masters and Scholars of the
00006 University of Oxford, having an administrative office at Wellington
00007 Square, Oxford OX1 2JD, UK.
00008 
00009 This file is part of Chaste.
00010 
00011 Chaste is free software: you can redistribute it and/or modify it
00012 under the terms of the GNU Lesser General Public License as published
00013 by the Free Software Foundation, either version 2.1 of the License, or
00014 (at your option) any later version.
00015 
00016 Chaste is distributed in the hope that it will be useful, but WITHOUT
00017 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
00019 License for more details. The offer of Chaste under the terms of the
00020 License is subject to the License being interpreted in accordance with
00021 English Law and subject to any action against the University of Oxford
00022 being under the jurisdiction of the English Courts.
00023 
00024 You should have received a copy of the GNU Lesser General Public License
00025 along with Chaste. If not, see <http://www.gnu.org/licenses/>.
00026 
00027 */
00028 
00029 #include "ImplicitCardiacMechanicsAssembler.hpp"
00030 #include "Kerchoffs2003ContractionModel.hpp"
00031 #include "NhsModelWithBackwardSolver.hpp"
00032 #include "NonPhysiologicalContractionModel.hpp"
00033 
00034 template<unsigned DIM>
00035 ImplicitCardiacMechanicsAssembler<DIM>::ImplicitCardiacMechanicsAssembler(
00036                                   ContractionModel contractionModel,
00037                                   QuadraticMesh<DIM>* pQuadMesh,
00038                                   std::string outputDirectory,
00039                                   std::vector<unsigned>& rFixedNodes,
00040                                   AbstractIncompressibleMaterialLaw<DIM>* pMaterialLaw)
00041     : AbstractCardiacMechanicsAssembler<DIM>(pQuadMesh,
00042                                              outputDirectory,
00043                                              rFixedNodes,
00044                                              pMaterialLaw)
00045 {
00046     switch(contractionModel)
00047     {
00048         case NONPHYSIOL1:
00049         case NONPHYSIOL2:
00050         case NONPHYSIOL3:
00051         {
00052             unsigned option = (contractionModel==NONPHYSIOL1 ? 1 : (contractionModel==NONPHYSIOL2? 2 : 3));
00053             for(unsigned i=0; i<this->mTotalQuadPoints; i++)
00054             {
00055                 this->mContractionModelSystems.push_back(new NonPhysiologicalContractionModel(option));
00056             }
00057             break;
00058         }
00059         case NHS:
00060         {
00061             for(unsigned i=0; i<this->mTotalQuadPoints; i++)
00062             {
00063                 this->mContractionModelSystems.push_back(new NhsModelWithBackwardSolver);
00064             }
00065             break;
00066         }
00067         case KERCHOFFS2003: //stretch dependent
00068         {
00069             for(unsigned i=0; i<this->mTotalQuadPoints; i++)
00070             {
00071                 this->mContractionModelSystems.push_back(new Kerchoffs2003ContractionModel());
00072             }
00073             break;
00074         }
00075         default:
00076         {
00077             #define COVERAGE_IGNORE // currently all available contraction models are acceptable for implicit
00078             EXCEPTION("Unknown or disallowed contraction model");
00079             #undef COVERAGE_IGNORE
00080         }
00081     }
00082 
00083     // initialise stores
00084     mStretchesLastTimeStep.resize(this->mTotalQuadPoints, 1.0);
00085 }
00086 
00087 template<unsigned DIM>
00088 ImplicitCardiacMechanicsAssembler<DIM>::~ImplicitCardiacMechanicsAssembler()
00089 {
00090     for(unsigned i=0; i<this->mContractionModelSystems.size(); i++)
00091     {
00092         delete this->mContractionModelSystems[i];
00093     }
00094 }
00095 
00096 
00097 template<unsigned DIM>
00098 std::vector<double>& ImplicitCardiacMechanicsAssembler<DIM>::rGetFibreStretches()
00099 {
00100     return this->mStretches;
00101 }
00102 
00103 
00104 template<unsigned DIM>
00105 void ImplicitCardiacMechanicsAssembler<DIM>::Solve(double time, double nextTime, double odeTimestep)
00106 {
00107     // set the times, which are used in AssembleOnElement
00108     assert(time < nextTime);
00109     this->mCurrentTime = time;
00110     this->mNextTime = nextTime;
00111     this->mOdeTimestep = odeTimestep;
00112 
00113     // solve
00114     NonlinearElasticityAssembler<DIM>::Solve();
00115 
00116     // assemble residual again (to solve the cell models implicitly again
00117     // using the correct value of the deformation x (in case this wasn't the
00118     // last thing that was done
00119     if(this->GetNumNewtonIterations() > 0)
00120     {
00121         this->AssembleSystem(true,false);
00122     }
00123 
00124     // now update state variables, and set lambda at last timestep. Note
00125     // lambda was set in AssembleOnElement
00126     for(unsigned i=0; i<this->mContractionModelSystems.size(); i++)
00127     {
00128          this->mContractionModelSystems[i]->UpdateStateVariables();
00129          mStretchesLastTimeStep[i] = this->mStretches[i];
00130     }
00131 }
00132 
00133 
00134 
00135 template<unsigned DIM>
00136 void ImplicitCardiacMechanicsAssembler<DIM>::GetActiveTensionAndTensionDerivs(double currentFibreStretch,
00137                                                                               unsigned currentQuadPointGlobalIndex,
00138                                                                               bool assembleJacobian,
00139                                                                               double& rActiveTension,
00140                                                                               double& rDerivActiveTensionWrtLambda,
00141                                                                               double& rDerivActiveTensionWrtDLambdaDt)
00142 {
00143     // save this fibre stretch
00144     this->mStretches[currentQuadPointGlobalIndex] = currentFibreStretch;
00145 
00146     // compute dlam/dt
00147     double dlam_dt = (currentFibreStretch-mStretchesLastTimeStep[currentQuadPointGlobalIndex])/(this->mNextTime-this->mCurrentTime);
00148 
00149     AbstractContractionModel* p_contraction_model = this->mContractionModelSystems[currentQuadPointGlobalIndex];
00150 
00151     // Set this stretch and stretch rate
00152     p_contraction_model->SetStretchAndStretchRate(currentFibreStretch, dlam_dt);
00153 
00154     // Call RunDoNotUpdate() on the contraction model to solve it using this stretch, and get the active tension
00155     try
00156     {
00157         p_contraction_model->RunDoNotUpdate(this->mCurrentTime,this->mNextTime,this->mOdeTimestep);
00158         rActiveTension = p_contraction_model->GetNextActiveTension();
00159     }
00160     catch (Exception& e)
00161     {
00162         #define COVERAGE_IGNORE
00163         // if this failed during assembling the Jacobian this is a fatal error.
00164         if(assembleJacobian)
00165         {
00166             // probably shouldn't be able to get here
00167             EXCEPTION("Failure in solving contraction models using current stretches for assembling Jacobian");
00168         }
00169         // if this failed during assembling the residual, the stretches are too large, so we just
00170         // set the active tension to infinity so that the residual will be infinite
00171         rActiveTension = DBL_MAX;
00172         std::cout << "WARNING: could not solve contraction model with this stretch and stretch rate. "
00173                   << "Setting active tension to infinity (DBL_MAX) so that the residual(-norm) is also infinite\n" << std::flush;
00174         assert(0); // just to see if we ever get here, can be removed..
00175         return;
00176         #undef COVERAGE_IGNORE
00177     }
00178 
00179     // if assembling the Jacobian, numerically evaluate dTa/dlam & dTa/d(lamdot)
00180     if(assembleJacobian)
00181     {
00182         // get active tension for (lam+h,dlamdt)
00183         double h1 = std::max(1e-6, currentFibreStretch/100);
00184         p_contraction_model->SetStretchAndStretchRate(currentFibreStretch+h1, dlam_dt);
00185         p_contraction_model->RunDoNotUpdate(this->mCurrentTime,this->mNextTime,this->mOdeTimestep);
00186         double active_tension_at_lam_plus_h = p_contraction_model->GetNextActiveTension();
00187 
00188         // get active tension for (lam,dlamdt+h)
00189         double h2 = std::max(1e-6, dlam_dt/100);
00190         p_contraction_model->SetStretchAndStretchRate(currentFibreStretch, dlam_dt+h2);
00191         p_contraction_model->RunDoNotUpdate(this->mCurrentTime,this->mNextTime,this->mOdeTimestep);
00192         double active_tension_at_dlamdt_plus_h = p_contraction_model->GetNextActiveTension();
00193 
00194         rDerivActiveTensionWrtLambda = (active_tension_at_lam_plus_h - rActiveTension)/h1;
00195         rDerivActiveTensionWrtDLambdaDt = (active_tension_at_dlamdt_plus_h - rActiveTension)/h2;
00196     }
00197 
00198     // Re-set the stretch and stretch rate and recompute the active tension so that
00199     // if this guess turns out to the solution, we can just update the state variables
00200     //  -- not needed as AssembleSystem(true,false) [ie assemble residual] is
00201     //     called in ImplicitCardiacMechanicsAssembler<DIM>::Solve() above
00202     //     after the solve and before the update.
00203     //  -- The SetActiveTensionInitialGuess() would make this very fast
00204     //     (compared to AssembleSystem(true,false) above), but the NHS class uses the last
00205     //     active tension as the initial guess anyway..
00206     //p_contraction_model->SetStretchAndStretchRate(currentFibreStretch, dlam_dt);
00207     //p_contraction_model->SetActiveTensionInitialGuess(rActiveTension);   // this line would now need a cast
00208     //p_contraction_model->RunDoNotUpdate(this->mCurrentTime,this->mNextTime,this->mOdeTimestep);
00209     //assert( fabs(p_contraction_model->GetNextActiveTension()-rActiveTension)<1e-8);
00210 }
00211 
00212 
00213 
00214 template class ImplicitCardiacMechanicsAssembler<2>;
00215 template class ImplicitCardiacMechanicsAssembler<3>;
00216 
00217 

Generated by  doxygen 1.6.2