ExplicitCardiacMechanicsSolver.cpp

00001 /*
00002 
00003 Copyright (C) University of Oxford, 2005-2011
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 "ExplicitCardiacMechanicsSolver.hpp"
00030 #include "Nash2004ContractionModel.hpp"
00031 #include "Kerchoffs2003ContractionModel.hpp"
00032 #include "NonPhysiologicalContractionModel.hpp"
00033 
00035 //#include "NhsModelWithBackSolver.hpp"
00036 
00037 template<unsigned DIM>
00038 ExplicitCardiacMechanicsSolver<DIM>::ExplicitCardiacMechanicsSolver(ContractionModel contractionModel,
00039                                                                     QuadraticMesh<DIM>* pQuadMesh,
00040                                                                     std::string outputDirectory,
00041                                                                     std::vector<unsigned>& rFixedNodes,
00042                                                                     AbstractIncompressibleMaterialLaw<DIM>* pMaterialLaw)
00043     : AbstractCardiacMechanicsSolver<DIM>(pQuadMesh,
00044                                           outputDirectory,
00045                                           rFixedNodes,
00046                                           pMaterialLaw)
00047 {
00048 //    // for showing stretch-rate-dependent models won't work with explicit (should be commented if committed)
00049 //    mStretchesLastTimestep.resize(this->mTotalQuadPoints, 1.0);
00050 
00051     switch(contractionModel)
00052     {
00053         case NONPHYSIOL1:
00054         case NONPHYSIOL2:
00055         case NONPHYSIOL3:
00056         {
00057             unsigned option = (contractionModel==NONPHYSIOL1 ? 1 : (contractionModel==NONPHYSIOL2? 2 : 3));
00058             for(unsigned i=0; i<this->mTotalQuadPoints; i++)
00059             {
00060                 this->mContractionModelSystems.push_back(new NonPhysiologicalContractionModel(option));
00061             }
00062             break;
00063         }
00064         case NASH2004: //stretch dependent, will this work with explicit??
00065         {
00066             for(unsigned i=0; i<this->mTotalQuadPoints; i++)
00067             {
00068                 this->mContractionModelSystems.push_back(new Nash2004ContractionModel);
00069             }
00070             break;
00071         }
00072         case KERCHOFFS2003: //stretch dependent, will this work with explicit? Answer: can be unstable
00073         {
00074             for(unsigned i=0; i<this->mTotalQuadPoints; i++)
00075             {
00076                 this->mContractionModelSystems.push_back(new Kerchoffs2003ContractionModel);
00077             }
00078             break;
00079         }
00080 //        // for showing stretch-rate-dependent models won't work with explicit (should be commented if committed)
00081 //        case NHS:
00082 //        {
00083 //            for(unsigned i=0; i<this->mTotalQuadPoints; i++)
00084 //            {
00085 //                this->mContractionModelSystems.push_back(new NhsModelWithBackwardSolver);
00086 //            }
00087 //            break;
00088 //        }
00089         default:
00090         {
00091             EXCEPTION("Unknown or stretch-rate-dependent contraction model");
00092         }
00093     }
00094 
00095     assert(!(this->mContractionModelSystems[0]->IsStretchRateDependent()));
00096 }
00097 
00098 template<unsigned DIM>
00099 ExplicitCardiacMechanicsSolver<DIM>::~ExplicitCardiacMechanicsSolver()
00100 {
00101     for(unsigned i=0; i<this->mContractionModelSystems.size(); i++)
00102     {
00103         delete this->mContractionModelSystems[i];
00104     }
00105 }
00106 
00107 template<unsigned DIM>
00108 void ExplicitCardiacMechanicsSolver<DIM>::GetActiveTensionAndTensionDerivs(double currentFibreStretch,
00109                                                                            unsigned currentQuadPointGlobalIndex,
00110                                                                            bool assembleJacobian,
00111                                                                            double& rActiveTension,
00112                                                                            double& rDerivActiveTensionWrtLambda,
00113                                                                            double& rDerivActiveTensionWrtDLambdaDt)
00114 {
00115     // the active tensions have already been computed for each contraction model, so can
00116     // return it straightaway..
00117     rActiveTension = this->mContractionModelSystems[currentQuadPointGlobalIndex]->GetActiveTension();
00118 
00119     // these are unset
00120     rDerivActiveTensionWrtLambda = 0.0;
00121     rDerivActiveTensionWrtDLambdaDt = 0.0;
00122 
00123     // store the value of given for this quad point, so that it can be used when computing
00124     // the active tension at the next timestep
00125     this->mStretches[currentQuadPointGlobalIndex] = currentFibreStretch;
00126 }
00127 
00128 template<unsigned DIM>
00129 void ExplicitCardiacMechanicsSolver<DIM>::Solve(double time, double nextTime, double odeTimestep)
00130 {
00131     assert(time < nextTime);
00132     this->mCurrentTime = time;
00133     this->mNextTime = nextTime;
00134     this->mOdeTimestep = odeTimestep;
00135 
00136     // assemble the residual again so that mStretches is set (in GetActiveTensionAndTensionDerivs)
00137     // using the current deformation.
00138     this->AssembleSystem(true,false);
00139 
00140     // integrate contraction models
00141     for(unsigned i=0; i<this->mContractionModelSystems.size(); i++)
00142     {
00143 //        // for showing stretch-rate-dependent models won't work with explicit (should be commented if committed)
00144 //        double dlam_dt = (this->mStretches[i] - mStretchesLastTimestep[i])/(nextTime-time);
00145 
00146         this->mContractionModelSystems[i]->SetStretchAndStretchRate(this->mStretches[i], 0.0 /*dlam_dt*/);
00147         this->mContractionModelSystems[i]->RunAndUpdate(time, nextTime, odeTimestep);
00148     }
00149 
00150     // solve
00151     NonlinearElasticitySolver<DIM>::Solve();
00152 
00153 //    // for showing stretch-rate-dependent models won't work with explicit (should be commented if committed)
00154 //    mStretchesLastTimestep = this->mStretches;
00155 }
00156 
00157 
00158 
00159 template class ExplicitCardiacMechanicsSolver<2>;
00160 template class ExplicitCardiacMechanicsSolver<3>;

Generated on Mon Apr 18 11:35:28 2011 for Chaste by  doxygen 1.5.5