AbstractCardiacMechanicsSolver.hpp

/*

Copyright (C) University of Oxford, 2005-2011

University of Oxford means the Chancellor, Masters and Scholars of the
University of Oxford, having an administrative office at Wellington
Square, Oxford OX1 2JD, UK.

This file is part of Chaste.

Chaste is free software: you can redistribute it and/or modify it
under the terms of the GNU Lesser General Public License as published
by the Free Software Foundation, either version 2.1 of the License, or
(at your option) any later version.

Chaste is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
License for more details. The offer of Chaste under the terms of the
License is subject to the License being interpreted in accordance with
English Law and subject to any action against the University of Oxford
being under the jurisdiction of the English Courts.

You should have received a copy of the GNU Lesser General Public License
along with Chaste. If not, see <http://www.gnu.org/licenses/>.

*/

#ifndef ABSTRACTCARDIACMECHANICSSOLVER_HPP_
#define ABSTRACTCARDIACMECHANICSSOLVER_HPP_

#include "NonlinearElasticitySolver.hpp"
#include "NashHunterPoleZeroLaw.hpp"
#include "QuadraticBasisFunction.hpp" // not included in NonlinearElasticitySolver.hpp, just the cpp
#include "LinearBasisFunction.hpp"
#include "AbstractContractionModel.hpp"
#include "FibreReader.hpp"

#include "AbstractCardiacMechanicsSolverInterface.hpp"

template<class ELASTICITY_SOLVER, unsigned DIM>
class AbstractCardiacMechanicsSolver : public ELASTICITY_SOLVER, public AbstractCardiacMechanicsSolverInterface<DIM>
{
protected:
    static const unsigned NUM_VERTICES_PER_ELEMENT = ELASTICITY_SOLVER::NUM_VERTICES_PER_ELEMENT; 
    std::vector<AbstractContractionModel*> mContractionModelSystems;

    std::vector<double> mStretches;

    unsigned mTotalQuadPoints;

    double mCurrentTime;
    double mNextTime;
    double mOdeTimestep;

    c_matrix<double,DIM,DIM> mConstantFibreSheetDirections;

    std::vector<c_matrix<double,DIM,DIM> >* mpVariableFibreSheetDirections;
    
    bool mFibreSheetDirectionsDefinedByQuadraturePoint;

    c_matrix<double,DIM,DIM>* mpCurrentElementFibreSheetMatrix;
    c_vector<double,DIM> mCurrentElementFibreDirection;

    virtual bool IsImplicitSolver()=0;

    void ComputeStressAndStressDerivative(AbstractMaterialLaw<DIM>* pMaterialLaw,
                                          c_matrix<double,DIM,DIM>& rC, 
                                          c_matrix<double,DIM,DIM>& rInvC, 
                                          double pressure,
                                          unsigned elementIndex,
                                          unsigned currentQuadPointGlobalIndex,
                                          c_matrix<double,DIM,DIM>& rT,
                                          FourthOrderTensor<DIM,DIM,DIM,DIM>& rDTdE,
                                          bool computeDTdE);
    



    virtual void GetActiveTensionAndTensionDerivs(double currentFibreStretch,
                                                  unsigned currentQuadPointGlobalIndex,
                                                  bool assembleJacobian,
                                                  double& rActiveTension,
                                                  double& rDerivActiveTensionWrtLambda,
                                                  double& rDerivActiveTensionWrtDLambdaDt)=0;
public:
    AbstractCardiacMechanicsSolver(QuadraticMesh<DIM>* pQuadMesh,
                                   std::string outputDirectory,
                                   std::vector<unsigned>& rFixedNodes,
                                   AbstractMaterialLaw<DIM>* pMaterialLaw);


    ~AbstractCardiacMechanicsSolver();
    
    
    unsigned GetTotalNumQuadPoints()
    {
        return mTotalQuadPoints;
    }

    virtual GaussianQuadratureRule<DIM>* GetQuadratureRule()
    {
        return this->mpQuadratureRule;
    }


    void SetConstantFibreSheetDirections(const c_matrix<double,DIM,DIM>& rFibreSheetMatrix);

    void SetVariableFibreSheetDirections(std::string orthoFile, bool definedPerQuadraturePoint);



    void SetCalciumAndVoltage(std::vector<double>& rCalciumConcentrations,
                              std::vector<double>& rVoltages);

    virtual void Solve(double time, double nextTime, double odeTimestep)=0;
    
    

    void ComputeDeformationGradientAndStretchInEachElement(std::vector<c_matrix<double,DIM,DIM> >& rDeformationGradients,
                                                           std::vector<double>& rStretches);
};


template<class ELASTICITY_SOLVER,unsigned DIM>
AbstractCardiacMechanicsSolver<ELASTICITY_SOLVER,DIM>::AbstractCardiacMechanicsSolver(QuadraticMesh<DIM>* pQuadMesh,
                                                                                      std::string outputDirectory,
                                                                                      std::vector<unsigned>& rFixedNodes,
                                                                                      AbstractMaterialLaw<DIM>* pMaterialLaw)
   : ELASTICITY_SOLVER(pQuadMesh,
                       pMaterialLaw,
                       zero_vector<double>(DIM),
                       DOUBLE_UNSET,
                       outputDirectory,
                       rFixedNodes),
     mCurrentTime(DBL_MAX),
     mNextTime(DBL_MAX),
     mOdeTimestep(DBL_MAX)
{
    // compute total num quad points
    mTotalQuadPoints = pQuadMesh->GetNumElements()*this->mpQuadratureRule->GetNumQuadPoints();

    // initialise the store of fibre stretches
    mStretches.resize(mTotalQuadPoints, 1.0);

    // initialise fibre/sheet direction matrix to be the identity, fibres in X-direction, and sheet in XY-plane
    mConstantFibreSheetDirections = zero_matrix<double>(DIM,DIM);
    for(unsigned i=0; i<DIM; i++)
    {
        mConstantFibreSheetDirections(i,i) = 1.0;
    }

    mpVariableFibreSheetDirections = NULL;
}


template<class ELASTICITY_SOLVER,unsigned DIM>
AbstractCardiacMechanicsSolver<ELASTICITY_SOLVER,DIM>::~AbstractCardiacMechanicsSolver()
{
    if(mpVariableFibreSheetDirections)
    {
        delete mpVariableFibreSheetDirections;
    }
}



template<class ELASTICITY_SOLVER,unsigned DIM>
void AbstractCardiacMechanicsSolver<ELASTICITY_SOLVER,DIM>::SetCalciumAndVoltage(std::vector<double>& rCalciumConcentrations,
                                                                                 std::vector<double>& rVoltages)

{
    assert(rCalciumConcentrations.size() == this->mTotalQuadPoints);
    assert(rVoltages.size() == this->mTotalQuadPoints);

    ContractionModelInputParameters input_parameters;

    for(unsigned i=0; i<rCalciumConcentrations.size(); i++)
    {
        input_parameters.intracellularCalciumConcentration = rCalciumConcentrations[i];
        input_parameters.voltage = rVoltages[i];

        mContractionModelSystems[i]->SetInputParameters(input_parameters);
    }
}

template<class ELASTICITY_SOLVER,unsigned DIM>
void AbstractCardiacMechanicsSolver<ELASTICITY_SOLVER,DIM>::ComputeStressAndStressDerivative(AbstractMaterialLaw<DIM>* pMaterialLaw,
                                                                                             c_matrix<double,DIM,DIM>& rC,
                                                                                             c_matrix<double,DIM,DIM>& rInvC,
                                                                                             double pressure,
                                                                                             unsigned elementIndex,
                                                                                             unsigned currentQuadPointGlobalIndex,
                                                                                             c_matrix<double,DIM,DIM>& rT,
                                                                                             FourthOrderTensor<DIM,DIM,DIM,DIM>& rDTdE,
                                                                                             bool assembleJacobian)
{
    if(!mpVariableFibreSheetDirections) // constant fibre directions
    {
        mpCurrentElementFibreSheetMatrix = &mConstantFibreSheetDirections;
    }
    else if(!mFibreSheetDirectionsDefinedByQuadraturePoint) // fibre directions defined for each mechanics mesh element
    {
        mpCurrentElementFibreSheetMatrix = &(*mpVariableFibreSheetDirections)[elementIndex];
    }
    else // fibre directions defined for each mechanics mesh quadrature point
    {
        mpCurrentElementFibreSheetMatrix = &(*mpVariableFibreSheetDirections)[currentQuadPointGlobalIndex];
    }
        
    for(unsigned i=0; i<DIM; i++)
    {
        mCurrentElementFibreDirection(i) = (*mpCurrentElementFibreSheetMatrix)(i,0);
    }

    // 1. Compute T and dTdE for the PASSIVE part of the strain energy.
    pMaterialLaw->SetChangeOfBasisMatrix(*mpCurrentElementFibreSheetMatrix);
    pMaterialLaw->ComputeStressAndStressDerivative(rC,rInvC,pressure,rT,rDTdE,assembleJacobian);

    // 2. Compute the active tension and add to the stress and stress-derivative
    double I4 = inner_prod(mCurrentElementFibreDirection, prod(rC, mCurrentElementFibreDirection));
    double lambda = sqrt(I4);

    double active_tension = 0;
    double d_act_tension_dlam = 0.0;     // Set and used if assembleJacobian==true
    double d_act_tension_d_dlamdt = 0.0; // Set and used if assembleJacobian==true

    GetActiveTensionAndTensionDerivs(lambda, currentQuadPointGlobalIndex, assembleJacobian,
                                     active_tension, d_act_tension_dlam, d_act_tension_d_dlamdt);

    // amend the stress and dTdE using the active tension
    double dTdE_coeff = -2*active_tension/(I4*I4); // note: I4*I4 = lam^4
    if(IsImplicitSolver())
    {
        double dt = mNextTime-mCurrentTime;
        dTdE_coeff += (d_act_tension_dlam + d_act_tension_d_dlamdt/dt)/(lambda*I4); // note: I4*lam = lam^3
    }

    rT += (active_tension/I4)*outer_prod(mCurrentElementFibreDirection,mCurrentElementFibreDirection);

    if(assembleJacobian)
    {
        for (unsigned M=0; M<DIM; M++)
        {
            for (unsigned N=0; N<DIM; N++)
            {
                for (unsigned P=0; P<DIM; P++)
                {
                    for (unsigned Q=0; Q<DIM; Q++)
                    {
                        rDTdE(M,N,P,Q) +=  dTdE_coeff * mCurrentElementFibreDirection(M)
                                                      * mCurrentElementFibreDirection(N)
                                                      * mCurrentElementFibreDirection(P)
                                                      * mCurrentElementFibreDirection(Q);
                    }
                }
            }
        }
    }
}




template<class ELASTICITY_SOLVER,unsigned DIM>
void AbstractCardiacMechanicsSolver<ELASTICITY_SOLVER,DIM>::ComputeDeformationGradientAndStretchInEachElement(
    std::vector<c_matrix<double,DIM,DIM> >& rDeformationGradients,
    std::vector<double>& rStretches)
{
    assert(rStretches.size()==this->mpQuadMesh->GetNumElements());
    
    // this will only work currently if the coarse mesh fibre info is defined per element, not per quad point 
    assert(!mpVariableFibreSheetDirections || !mFibreSheetDirectionsDefinedByQuadraturePoint);
   
    static c_matrix<double,DIM,NUM_VERTICES_PER_ELEMENT> element_current_displacements;
    static c_matrix<double,DIM,NUM_VERTICES_PER_ELEMENT> grad_lin_phi;
    static c_matrix<double,DIM,DIM> F;      // the deformation gradient, F = dx/dX, F_{iM} = dx_i/dX_M

    static c_matrix<double,DIM,DIM> jacobian;
    static c_matrix<double,DIM,DIM> inverse_jacobian;
    double jacobian_determinant;
    ChastePoint<DIM> quadrature_point; // not needed, but has to be passed in  
    
    // loop over all the elements
    for(unsigned elem_index=0; elem_index<this->mpQuadMesh->GetNumElements(); elem_index++)
    {
        Element<DIM,DIM>* p_elem = this->mpQuadMesh->GetElement(elem_index);

        // get the fibre direction for this element
        mpCurrentElementFibreSheetMatrix = mpVariableFibreSheetDirections ? &(*mpVariableFibreSheetDirections)[elem_index] : &mConstantFibreSheetDirections;
        for(unsigned i=0; i<DIM; i++)
        {
            mCurrentElementFibreDirection(i) = (*mpCurrentElementFibreSheetMatrix)(i,0);
        }
        
        // get the displacement at this element
        for (unsigned II=0; II<NUM_VERTICES_PER_ELEMENT; II++)
        {
            for (unsigned JJ=0; JJ<DIM; JJ++)
            {
                element_current_displacements(JJ,II) = this->mCurrentSolution[DIM*p_elem->GetNodeGlobalIndex(II) + JJ];
            }
        }

        // set up the linear basis functions
        this->mpQuadMesh->GetInverseJacobianForElement(elem_index, jacobian, jacobian_determinant, inverse_jacobian);
        LinearBasisFunction<DIM>::ComputeTransformedBasisFunctionDerivatives(quadrature_point, inverse_jacobian, grad_lin_phi);
        
        F = identity_matrix<double>(DIM,DIM); 

        // loop over the vertices and interpolate F, the deformation gradient
        // (note: could use matrix-mult if this becomes inefficient
        for (unsigned node_index=0; node_index<NUM_VERTICES_PER_ELEMENT; node_index++)
        {
            for (unsigned i=0; i<DIM; i++)
            {
                for (unsigned M=0; M<DIM; M++)
                {
                    F(i,M) += grad_lin_phi(M,node_index)*element_current_displacements(i,node_index);
                }
            }
        }

        rDeformationGradients[elem_index] = F;
        
        // compute and save the stretch: m^T C m = ||Fm||
        c_vector<double,DIM> deformed_fibre = prod(F, mCurrentElementFibreDirection);
        rStretches[elem_index] = norm_2(deformed_fibre);
    }
}




template<class ELASTICITY_SOLVER,unsigned DIM>
void AbstractCardiacMechanicsSolver<ELASTICITY_SOLVER,DIM>::SetVariableFibreSheetDirections(std::string orthoFile, bool definedPerQuadraturePoint)
{
    mFibreSheetDirectionsDefinedByQuadraturePoint = definedPerQuadraturePoint;
    
    FileFinder finder(orthoFile, RelativeTo::ChasteSourceRoot);
    FibreReader<DIM> reader(finder, ORTHO);
    
    unsigned num_entries = reader.GetNumLinesOfData();

    if(!mFibreSheetDirectionsDefinedByQuadraturePoint && (num_entries!=this->mpQuadMesh->GetNumElements()) )
    {
        std::stringstream ss;
        ss << "Number of entries defined at top of file " << orthoFile << " does not match number of elements in the mesh, "
           << "found " <<  num_entries << ", expected " << this->mpQuadMesh->GetNumElements();
        EXCEPTION(ss.str());
    }

    if(mFibreSheetDirectionsDefinedByQuadraturePoint && (num_entries!=mTotalQuadPoints) )
    {
        std::stringstream ss;
        ss << "Number of entries defined at top of file " << orthoFile << " does not match number of quadrature points defined, "
           << "found " <<  num_entries << ", expected " << mTotalQuadPoints;
        EXCEPTION(ss.str());
    }
    
    mpVariableFibreSheetDirections = new std::vector<c_matrix<double,DIM,DIM> >(num_entries, zero_matrix<double>(DIM,DIM));
    for(unsigned index=0; index<num_entries; index++)
    {
        reader.GetNextFibreSheetAndNormalMatrix( (*mpVariableFibreSheetDirections)[index] );
    }
}



template<class ELASTICITY_SOLVER,unsigned DIM>
void AbstractCardiacMechanicsSolver<ELASTICITY_SOLVER,DIM>::SetConstantFibreSheetDirections(const c_matrix<double,DIM,DIM>& rFibreSheetMatrix)
{
    mConstantFibreSheetDirections = rFibreSheetMatrix;
    // check orthogonality
    c_matrix<double,DIM,DIM>  temp = prod(trans(rFibreSheetMatrix),rFibreSheetMatrix);
    for(unsigned i=0; i<DIM; i++)
    {
        for(unsigned j=0; j<DIM; j++)
        {
            double val = (i==j ? 1.0 : 0.0);
            if(fabs(temp(i,j)-val)>1e-4)
            {
                std::stringstream string_stream;
                string_stream << "The given fibre-sheet matrix, " << rFibreSheetMatrix << ", is not orthogonal";
                EXCEPTION(string_stream.str());
            }
        }
    }
}


#endif /*ABSTRACTCARDIACMECHANICSSOLVER_HPP_*/

---