cardioModel_mod.f90

Go to the documentation of this file.

module cardiomodel_mod

  use choral_constants
  use choral_variables
  use real_type, only: rp
  use basic_tools
  use funclib
  use abstract_interfaces, only: r3tor3
  use algebra_lin
  use fespace_mod
  use quadmesh_mod
  use diffusion
  use csr_mod
 
  implicit none
  private


  !       %----------------------------------------%
  !       |                                        |
  !       |          PUBLIC DATA                   |
  !       |                                        |
  !       %----------------------------------------%
  public :: cardiomodel

  public :: print
  public :: monodomain_assemble
  public :: scale_conductivities


  !       %----------------------------------------%
  !       |                                        |
  !       |          DERIVED TYPE                  |
  !       |                                        |
  !       %----------------------------------------%
  type :: cardiomodel

     !! data provided by the constructor
     !!
     procedure(r3tor3) , pointer, nopass :: fibre => vector_field_e_x

     real(RP)        :: am

     !! Data defined by the code  
     !!
     real(RP) :: gil
     real(RP) :: git
     real(RP) :: gel
     real(RP) :: get

  end type cardiomodel


  !       %----------------------------------------%
  !       |                                        |
  !       |          GENERIC INTERFACES            |
  !       |                                        |
  !       %----------------------------------------%
  interface print
     module procedure cardiomodel_print
  end interface print

  interface cardiomodel
     module procedure cardiomodel_create
  end interface cardiomodel

contains

  function cardiomodel_create(fibre, Am, conduc_type) result(cm)

    type(cardiomodel)             :: cm
    procedure(r3tor3)             :: fibre
    real(RP)         , intent(in) :: am
    integer          , intent(in) :: conduc_type

    if (choral_verb>0) write(*,*) &
         & "cardioModel_mod : create"

    cm%fibre => fibre

    cm%Am = am

    call set_conductivities(cm%GIL, &
         &cm%GIT, cm%GEL, cm%GET, conduc_type)

  end function cardiomodel_create


  subroutine set_conductivities(gil, git, gel, get, type)
    real(RP), intent(out) :: gil, git,gel,get
    integer , intent(in)  :: type

    if (choral_verb>1) write(*,"(A36,$)") &
         & " cardioModel_mod : set_conduct.   = "
    select case(type)
    case(le_guyader)
       if (choral_verb>1)  write(*,*) "LE_GUYADER"
       gil = 1.741_rp
       git = 0.1934_rp
       gel = 3.906_rp
       get = 1.970_rp

    case(le_guyader_2)
       if (choral_verb>1)  write(*,*) "LE_GUYADER_2"
       gil = 1.741_rp
       git = 0.475_rp
       gel = 3.906_rp
       get = 1.970_rp

    case(cond_iso_1)
       if (choral_verb>1)  write(*,*) "COND_ISO_1"
       gil = 1._rp
       git = 1._rp
       gel = 1._rp
       get = 1._rp

    case default
       call quit("cardioModel_mod: set_conductivities: 1")

    end select

  end subroutine set_conductivities


  subroutine scale_conductivities(cm, sc)
    type(cardiomodel), intent(inout) :: cm
    real(RP)         , intent(in)    :: sc

    cm%GIL = cm%GIL * sc
    cm%GIT = cm%GIT * sc
    cm%GEL = cm%GEL * sc
    cm%GET = cm%GET * sc

  end subroutine scale_conductivities


  subroutine cardiomodel_print(cm)
    type(cardiomodel), intent(in) :: cm

    write(*,*)"cardioModel_mod : print "
    write(*,*)"  Am                             =", cm%am
    write(*,*)"  Intra conduc. GIL              =", cm%GIL
    write(*,*)"                GIT              =", cm%GIT
    write(*,*)"  Extra conduc. GEL              =", cm%GEL
    write(*,*)"                GET              =", cm%GET

  end subroutine cardiomodel_print


  subroutine monodomain_assemble(M, S, cm, X_h, qd_M, qd_S)
    type(csr)         , intent(out):: M, S
    class(cardiomodel), intent(in) :: cm
    type(fespace)     , intent(in) :: X_h
    integer           , intent(in) :: qd_M, qd_S

    type(quadmesh) :: qdm
    real(RP) :: gl, gt

    if (choral_verb>0) write(*,*) &
         & "cardioModel_mod : monodomain"

    ! monodomain conductivities
    gl = 2._rp / (1._rp/cm%GIL + 1._rp/cm%GEL )
    gt = 2._rp / (1._rp/cm%GIT + 1._rp/cm%GET )
    if (choral_verb>1) write(*,*) &
         & "  conduc. long/trans             =", &
         & real(gl,SP), real(gt,SP) 

    ! mass matrix
    qdm = quadmesh(x_h%mesh)
    call set(qdm, qd_m)
    call assemble(qdm)
    call diffusion_massmat(m, one_r3, x_h, qdm)

    ! stiffness matrix 
    qdm = quadmesh(x_h%mesh)
    call set(qdm, qd_s)
    call assemble(qdm)
    call diffusion_stiffmat(s, b, x_h, qdm)    
    call scale(s, 1._rp/cm%Am)
    
  contains 

    function b(x, w1, w2) 
      real(RP)                           :: b
      real(RP), dimension(3), intent(in) :: x, w1, w2

      real(RP), dimension(3) :: fd

      fd = cm%fibre(x)
      b  = usigmav(w1, w2, fd, gl, gt)

    end function b

  end subroutine monodomain_assemble



  function usigmav(u, v, fd, gl, gt) result(r)
    real(RP)  :: r
    real(RP), dimension(3), intent(in)    :: u, v, fd
    real(RP)              , intent(in)    :: gl, gt 

    real(RP), dimension(3) :: e1, e2, e3
    real(RP) :: nrm
    
    nrm = sqrt(sum(fd*fd))
    e1  = fd / nrm

    e2(1) = e1(3) - e1(2)
    e2(2) = e1(1) - e1(3)
    e2(3) = e1(2) - e1(1)

    nrm = sqrt(sum(e2*e2))
    e2  = e2/nrm

    e3  = crossprod_3d(e1,e2)

    r =     sum(u*e1) * sum(v*e1) * gl
    r = r + sum(u*e2) * sum(v*e2) * gt
    r = r + sum(u*e3) * sum(v*e3) * gt

  end function usigmav


end module cardiomodel_mod