AbstractCardiacMechanicsSolver.hpp

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 #ifndef ABSTRACTCARDIACMECHANICSSOLVER_HPP_
00030 #define ABSTRACTCARDIACMECHANICSSOLVER_HPP_
00031 
00032 #include "NonlinearElasticitySolver.hpp"
00033 #include "NashHunterPoleZeroLaw.hpp"
00034 #include "QuadraticBasisFunction.hpp" // not included in NonlinearElasticitySolver.hpp, just the cpp
00035 #include "LinearBasisFunction.hpp"
00036 #include "AbstractContractionModel.hpp"
00037 #include "FibreReader.hpp"
00038 
00039 
00048 template<unsigned DIM>
00049 class AbstractCardiacMechanicsSolver : public NonlinearElasticitySolver<DIM>
00050 {
00051 protected:
00052     static const unsigned STENCIL_SIZE = NonlinearElasticitySolver<DIM>::STENCIL_SIZE;
00053     static const unsigned NUM_NODES_PER_ELEMENT = NonlinearElasticitySolver<DIM>::NUM_NODES_PER_ELEMENT;
00054     static const unsigned NUM_VERTICES_PER_ELEMENT = NonlinearElasticitySolver<DIM>::NUM_VERTICES_PER_ELEMENT;
00055 
00061     std::vector<AbstractContractionModel*> mContractionModelSystems;
00062 
00068     std::vector<double> mStretches;
00069 
00071     unsigned mTotalQuadPoints;
00072 
00074     bool mAllocatedMaterialLawMemory;
00075 
00077     double mCurrentTime;
00079     double mNextTime;
00081     double mOdeTimestep;
00082 
00084     c_matrix<double,DIM,DIM> mConstantFibreSheetDirections;
00085 
00090     std::vector<c_matrix<double,DIM,DIM> >* mpVariableFibreSheetDirections;
00091     
00096     bool mFibreSheetDirectionsDefinedByQuadraturePoint;
00097 
00099     c_matrix<double,DIM,DIM>* mpCurrentElementFibreSheetMatrix;
00101     c_vector<double,DIM> mCurrentElementFibreDirection;
00102 
00107     virtual bool IsImplicitSolver()=0;
00108 
00127     void ComputeStressAndStressDerivative(AbstractIncompressibleMaterialLaw<DIM>* pMaterialLaw,
00128                                           c_matrix<double,DIM,DIM>& rC, 
00129                                           c_matrix<double,DIM,DIM>& rInvC, 
00130                                           double pressure,
00131                                           unsigned elementIndex,
00132                                           unsigned currentQuadPointGlobalIndex,
00133                                           c_matrix<double,DIM,DIM>& rT,
00134                                           FourthOrderTensor<DIM,DIM,DIM,DIM>& rDTdE,
00135                                           bool computeDTdE);
00136     
00137 
00138 
00139 
00152     virtual void GetActiveTensionAndTensionDerivs(double currentFibreStretch,
00153                                                   unsigned currentQuadPointGlobalIndex,
00154                                                   bool assembleJacobian,
00155                                                   double& rActiveTension,
00156                                                   double& rDerivActiveTensionWrtLambda,
00157                                                   double& rDerivActiveTensionWrtDLambdaDt)=0;
00158 public:
00168     AbstractCardiacMechanicsSolver(QuadraticMesh<DIM>* pQuadMesh,
00169                                    std::string outputDirectory,
00170                                    std::vector<unsigned>& rFixedNodes,
00171                                    AbstractIncompressibleMaterialLaw<DIM>* pMaterialLaw);
00172 
00173 
00177     ~AbstractCardiacMechanicsSolver();
00178     
00179     
00181     unsigned GetTotalNumQuadPoints()
00182     {
00183         return mTotalQuadPoints;
00184     }
00185 
00187     virtual GaussianQuadratureRule<DIM>* GetQuadratureRule()
00188     {
00189         return this->mpQuadratureRule;
00190     }
00191 
00192 
00199     void SetConstantFibreSheetDirections(const c_matrix<double,DIM,DIM>& rFibreSheetMatrix);
00200 
00212     void SetVariableFibreSheetDirections(std::string orthoFile, bool definedPerQuadraturePoint);
00213 
00214 
00215 
00228     void SetCalciumAndVoltage(std::vector<double>& rCalciumConcentrations,
00229                               std::vector<double>& rVoltages);
00230 
00238     virtual void Solve(double time, double nextTime, double odeTimestep)=0;
00239     
00240     
00241 
00261     void ComputeDeformationGradientAndStretchInEachElement(std::vector<c_matrix<double,DIM,DIM> >& rDeformationGradients,
00262                                                            std::vector<double>& rStretches);
00263 };
00264 
00265 
00266 template<unsigned DIM>
00267 AbstractCardiacMechanicsSolver<DIM>::AbstractCardiacMechanicsSolver(QuadraticMesh<DIM>* pQuadMesh,
00268                                                                     std::string outputDirectory,
00269                                                                     std::vector<unsigned>& rFixedNodes,
00270                                                                     AbstractIncompressibleMaterialLaw<DIM>* pMaterialLaw)
00271    : NonlinearElasticitySolver<DIM>(pQuadMesh,
00272                                     pMaterialLaw!=NULL ? pMaterialLaw : new NashHunterPoleZeroLaw<DIM>,
00273                                     zero_vector<double>(DIM),
00274                                     DOUBLE_UNSET,
00275                                     outputDirectory,
00276                                     rFixedNodes),
00277                                     mCurrentTime(DBL_MAX),
00278                                     mNextTime(DBL_MAX),
00279                                     mOdeTimestep(DBL_MAX)
00280 {
00281     // compute total num quad points
00282     mTotalQuadPoints = pQuadMesh->GetNumElements()*this->mpQuadratureRule->GetNumQuadPoints();
00283 
00284     // note that if pMaterialLaw is NULL a new NashHunter law was sent to the
00285     // NonlinElas constuctor (see above)
00286     mAllocatedMaterialLawMemory = (pMaterialLaw==NULL);
00287 
00288     // initialise the store of fibre stretches
00289     mStretches.resize(mTotalQuadPoints, 1.0);
00290 
00291     // initialise fibre/sheet direction matrix to be the identity, fibres in X-direction, and sheet in XY-plane
00292     mConstantFibreSheetDirections = zero_matrix<double>(DIM,DIM);
00293     for(unsigned i=0; i<DIM; i++)
00294     {
00295         mConstantFibreSheetDirections(i,i) = 1.0;
00296     }
00297 
00298     mpVariableFibreSheetDirections = NULL;
00299 }
00300 
00301 
00302 template<unsigned DIM>
00303 AbstractCardiacMechanicsSolver<DIM>::~AbstractCardiacMechanicsSolver()
00304 {
00305     if(mAllocatedMaterialLawMemory)
00306     {
00307         assert(this->mMaterialLaws.size()==1); 
00308         delete this->mMaterialLaws[0];
00309     }
00310 
00311     if(mpVariableFibreSheetDirections)
00312     {
00313         delete mpVariableFibreSheetDirections;
00314     }
00315 }
00316 
00317 
00318 
00319 template<unsigned DIM>
00320 void AbstractCardiacMechanicsSolver<DIM>::SetCalciumAndVoltage(std::vector<double>& rCalciumConcentrations,
00321                                                                std::vector<double>& rVoltages)
00322 
00323 {
00324     assert(rCalciumConcentrations.size() == this->mTotalQuadPoints);
00325     assert(rVoltages.size() == this->mTotalQuadPoints);
00326 
00327     ContractionModelInputParameters input_parameters;
00328 
00329     for(unsigned i=0; i<rCalciumConcentrations.size(); i++)
00330     {
00331         input_parameters.intracellularCalciumConcentration = rCalciumConcentrations[i];
00332         input_parameters.voltage = rVoltages[i];
00333 
00334         mContractionModelSystems[i]->SetInputParameters(input_parameters);
00335     }
00336 }
00337 
00338 template<unsigned DIM>
00339 void AbstractCardiacMechanicsSolver<DIM>::ComputeStressAndStressDerivative(AbstractIncompressibleMaterialLaw<DIM>* pMaterialLaw,
00340                                                                            c_matrix<double,DIM,DIM>& rC, 
00341                                                                            c_matrix<double,DIM,DIM>& rInvC, 
00342                                                                            double pressure, 
00343                                                                            unsigned elementIndex,
00344                                                                            unsigned currentQuadPointGlobalIndex,
00345                                                                            c_matrix<double,DIM,DIM>& rT,
00346                                                                            FourthOrderTensor<DIM,DIM,DIM,DIM>& rDTdE,
00347                                                                            bool assembleJacobian)
00348 {
00349     if(!mpVariableFibreSheetDirections) // constant fibre directions
00350     {
00351         mpCurrentElementFibreSheetMatrix = &mConstantFibreSheetDirections;
00352     }
00353     else if(!mFibreSheetDirectionsDefinedByQuadraturePoint) // fibre directions defined for each mechanics mesh element
00354     {
00355         mpCurrentElementFibreSheetMatrix = &(*mpVariableFibreSheetDirections)[elementIndex];
00356     }
00357     else // fibre directions defined for each mechanics mesh quadrature point
00358     {
00359         mpCurrentElementFibreSheetMatrix = &(*mpVariableFibreSheetDirections)[currentQuadPointGlobalIndex];
00360     }
00361         
00362     for(unsigned i=0; i<DIM; i++)
00363     {
00364         mCurrentElementFibreDirection(i) = (*mpCurrentElementFibreSheetMatrix)(i,0);
00365     }
00366 
00367     // 1. Compute T and dTdE for the PASSIVE part of the strain energy.
00368     pMaterialLaw->SetChangeOfBasisMatrix(*mpCurrentElementFibreSheetMatrix);
00369     pMaterialLaw->ComputeStressAndStressDerivative(rC,rInvC,pressure,rT,rDTdE,assembleJacobian);
00370 
00371     // 2. Compute the active tension and add to the stress and stress-derivative
00372     double I4 = inner_prod(mCurrentElementFibreDirection, prod(rC, mCurrentElementFibreDirection));
00373     double lambda = sqrt(I4);
00374 
00375     double active_tension = 0;
00376     double d_act_tension_dlam = 0.0;     // Set and used if assembleJacobian==true
00377     double d_act_tension_d_dlamdt = 0.0; // Set and used if assembleJacobian==true
00378 
00379     GetActiveTensionAndTensionDerivs(lambda, currentQuadPointGlobalIndex, assembleJacobian,
00380                                      active_tension, d_act_tension_dlam, d_act_tension_d_dlamdt);
00381 
00382     // amend the stress and dTdE using the active tension
00383     double dTdE_coeff = -2*active_tension/(I4*I4); // note: I4*I4 = lam^4
00384     if(IsImplicitSolver())
00385     {
00386         double dt = mNextTime-mCurrentTime;
00387         dTdE_coeff += (d_act_tension_dlam + d_act_tension_d_dlamdt/dt)/(lambda*I4); // note: I4*lam = lam^3
00388     }
00389 
00390     rT += (active_tension/I4)*outer_prod(mCurrentElementFibreDirection,mCurrentElementFibreDirection);
00391 
00392     if(assembleJacobian)
00393     {
00394         for (unsigned M=0; M<DIM; M++)
00395         {
00396             for (unsigned N=0; N<DIM; N++)
00397             {
00398                 for (unsigned P=0; P<DIM; P++)
00399                 {
00400                     for (unsigned Q=0; Q<DIM; Q++)
00401                     {
00402                         rDTdE(M,N,P,Q) +=  dTdE_coeff * mCurrentElementFibreDirection(M)
00403                                                       * mCurrentElementFibreDirection(N)
00404                                                       * mCurrentElementFibreDirection(P)
00405                                                       * mCurrentElementFibreDirection(Q);
00406                     }
00407                 }
00408             }
00409         }
00410     }
00411 }
00412 
00413 
00414 
00415 
00416 template<unsigned DIM>
00417 void AbstractCardiacMechanicsSolver<DIM>::ComputeDeformationGradientAndStretchInEachElement(
00418     std::vector<c_matrix<double,DIM,DIM> >& rDeformationGradients,
00419     std::vector<double>& rStretches)
00420 {
00421     assert(rStretches.size()==this->mpQuadMesh->GetNumElements());
00422     
00423     // this will only work currently if the coarse mesh fibre info is defined per element, not per quad point 
00424     assert(!mpVariableFibreSheetDirections || !mFibreSheetDirectionsDefinedByQuadraturePoint);
00425    
00426     static c_matrix<double,DIM,NUM_VERTICES_PER_ELEMENT> element_current_displacements;
00427     static c_matrix<double,DIM,NUM_VERTICES_PER_ELEMENT> grad_lin_phi;
00428     static c_matrix<double,DIM,DIM> F;      // the deformation gradient, F = dx/dX, F_{iM} = dx_i/dX_M
00429 
00430     static c_matrix<double,DIM,DIM> jacobian;
00431     static c_matrix<double,DIM,DIM> inverse_jacobian;
00432     double jacobian_determinant;
00433     ChastePoint<DIM> quadrature_point; // not needed, but has to be passed in  
00434     
00435     // loop over all the elements
00436     for(unsigned elem_index=0; elem_index<this->mpQuadMesh->GetNumElements(); elem_index++)
00437     {
00438         Element<DIM,DIM>* p_elem = this->mpQuadMesh->GetElement(elem_index);
00439 
00440         // get the fibre direction for this element
00441         mpCurrentElementFibreSheetMatrix = mpVariableFibreSheetDirections ? &(*mpVariableFibreSheetDirections)[elem_index] : &mConstantFibreSheetDirections;
00442         for(unsigned i=0; i<DIM; i++)
00443         {
00444             mCurrentElementFibreDirection(i) = (*mpCurrentElementFibreSheetMatrix)(i,0);
00445         }
00446         
00447         // get the displacement at this element
00448         for (unsigned II=0; II<NUM_VERTICES_PER_ELEMENT; II++)
00449         {
00450             for (unsigned JJ=0; JJ<DIM; JJ++)
00451             {
00452                 element_current_displacements(JJ,II) = this->mCurrentSolution[DIM*p_elem->GetNodeGlobalIndex(II) + JJ];
00453             }
00454         }
00455 
00456         // set up the linear basis functions
00457         this->mpQuadMesh->GetInverseJacobianForElement(elem_index, jacobian, jacobian_determinant, inverse_jacobian);
00458         LinearBasisFunction<DIM>::ComputeTransformedBasisFunctionDerivatives(quadrature_point, inverse_jacobian, grad_lin_phi);
00459         
00460         F = identity_matrix<double>(DIM,DIM); 
00461 
00462         // loop over the vertices and interpolate F, the deformation gradient
00463         // (note: could use matrix-mult if this becomes inefficient
00464         for (unsigned node_index=0; node_index<NUM_VERTICES_PER_ELEMENT; node_index++)
00465         {
00466             for (unsigned i=0; i<DIM; i++)
00467             {
00468                 for (unsigned M=0; M<DIM; M++)
00469                 {
00470                     F(i,M) += grad_lin_phi(M,node_index)*element_current_displacements(i,node_index);
00471                 }
00472             }
00473         }
00474 
00475         rDeformationGradients[elem_index] = F;
00476         
00477         // compute and save the stretch: m^T C m = ||Fm||
00478         c_vector<double,DIM> deformed_fibre = prod(F, mCurrentElementFibreDirection);
00479         rStretches[elem_index] = norm_2(deformed_fibre);
00480     }
00481 }
00482 
00483 
00484 
00485 
00486 template<unsigned DIM>
00487 void AbstractCardiacMechanicsSolver<DIM>::SetVariableFibreSheetDirections(std::string orthoFile, bool definedPerQuadraturePoint)
00488 {
00489     mFibreSheetDirectionsDefinedByQuadraturePoint = definedPerQuadraturePoint;
00490     
00491     FileFinder finder(orthoFile, RelativeTo::ChasteSourceRoot);
00492     FibreReader<DIM> reader(finder, ORTHO);
00493     
00494     unsigned num_entries = reader.GetNumLinesOfData();
00495 
00496     if(!mFibreSheetDirectionsDefinedByQuadraturePoint && (num_entries!=this->mpQuadMesh->GetNumElements()) )
00497     {
00498         std::stringstream ss;
00499         ss << "Number of entries defined at top of file " << orthoFile << " does not match number of elements in the mesh, "
00500            << "found " <<  num_entries << ", expected " << this->mpQuadMesh->GetNumElements();
00501         EXCEPTION(ss.str());
00502     }
00503 
00504     if(mFibreSheetDirectionsDefinedByQuadraturePoint && (num_entries!=mTotalQuadPoints) )
00505     {
00506         std::stringstream ss;
00507         ss << "Number of entries defined at top of file " << orthoFile << " does not match number of quadrature points defined, "
00508            << "found " <<  num_entries << ", expected " << mTotalQuadPoints;
00509         EXCEPTION(ss.str());
00510     }
00511     
00512     mpVariableFibreSheetDirections = new std::vector<c_matrix<double,DIM,DIM> >(num_entries, zero_matrix<double>(DIM,DIM));
00513     for(unsigned index=0; index<num_entries; index++)
00514     {
00515         reader.GetNextFibreSheetAndNormalMatrix( (*mpVariableFibreSheetDirections)[index] );
00516     }
00517 }
00518 
00519 
00520 
00521 template<unsigned DIM>
00522 void AbstractCardiacMechanicsSolver<DIM>::SetConstantFibreSheetDirections(const c_matrix<double,DIM,DIM>& rFibreSheetMatrix)
00523 {
00524     mConstantFibreSheetDirections = rFibreSheetMatrix;
00525     // check orthogonality
00526     c_matrix<double,DIM,DIM>  temp = prod(trans(rFibreSheetMatrix),rFibreSheetMatrix);
00527     for(unsigned i=0; i<DIM; i++)
00528     {
00529         for(unsigned j=0; j<DIM; j++)
00530         {
00531             double val = (i==j ? 1.0 : 0.0);
00532             if(fabs(temp(i,j)-val)>1e-4)
00533             {
00534                 std::stringstream string_stream;
00535                 string_stream << "The given fibre-sheet matrix, " << rFibreSheetMatrix << ", is not orthogonal";
00536                 EXCEPTION(string_stream.str());
00537             }
00538         }
00539     }
00540 }
00541 
00542 
00543 #endif /*ABSTRACTCARDIACMECHANICSSOLVER_HPP_*/

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