monodomain_equations_routines.f90

Go to the documentation of this file.

MODULE monodomain_equations_routines

  USE base_routines
  USE basis_routines
  USE boundary_conditions_routines
  USE constants
  USE control_loop_routines
  USE distributed_matrix_vector
  USE domain_mappings
  USE electrophysiology_cell_routines
  USE equations_routines
  USE equations_mapping_routines
  USE equations_matrices_routines
  USE equations_set_constants
  USE field_io_routines
  USE field_routines
  USE fitting_routines
  USE input_output
  USE iso_varying_string
  USE kinds
  USE matrix_vector
  USE node_routines
  USE problem_constants
  USE strings
  USE solver_routines
  USE timer
  USE types

#include "macros.h"  


  IMPLICIT NONE

  PRIVATE

  !Module parameters

  !Module types

  !Module variables

  !Interfaces

  PUBLIC monodomain_control_loop_post_loop
 
  PUBLIC monodomain_equation_equations_set_setup
  
  PUBLIC monodomain_finiteelementcalculate

  PUBLIC monodomain_equationssetsolutionmethodset
  
  PUBLIC monodomain_equationssetspecificationset

  PUBLIC monodomain_problemspecificationset
  
  PUBLIC monodomain_equation_problem_setup
  
  PUBLIC monodomain_pre_solve,monodomain_post_solve
  
CONTAINS

  !
  !================================================================================================================================
  !

  SUBROUTINE monodomain_control_loop_post_loop(CONTROL_LOOP,ERR,ERROR,*)

    !Argument variables
    TYPE(control_loop_type), POINTER :: CONTROL_LOOP
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) :: equations_set_idx
    TYPE(control_loop_time_type), POINTER :: TIME_LOOP,TIME_LOOP_PARENT
    TYPE(control_loop_type), POINTER :: PARENT_LOOP
    TYPE(equations_set_type), POINTER :: EQUATIONS_SET
    TYPE(field_type), POINTER :: DEPENDENT_FIELD
    TYPE(problem_type), POINTER :: PROBLEM
    TYPE(region_type), POINTER :: DEPENDENT_REGION   
    TYPE(solver_type), POINTER :: SOLVER
    TYPE(solver_equations_type), POINTER :: SOLVER_EQUATIONS
    TYPE(solver_mapping_type), POINTER :: SOLVER_MAPPING
    TYPE(solvers_type), POINTER :: SOLVERS
    TYPE(varying_string) :: FILENAME,LOCAL_ERROR,METHOD
    INTEGER(INTG) :: OUTPUT_ITERATION_NUMBER,CURRENT_LOOP_ITERATION

    enters("MONODOMAIN_CONTROL_LOOP_POST_LOOP",err,error,*999)

    IF(ASSOCIATED(control_loop)) THEN
      IF(control_loop%OUTPUT_TYPE>=control_loop_progress_output) THEN
        SELECT CASE(control_loop%LOOP_TYPE)
        CASE(problem_control_simple_type)
          !do nothing
        CASE(problem_control_fixed_loop_type)
          !do nothing
        CASE(problem_control_time_loop_type)
          !Export the dependent field for this time step
          time_loop=>control_loop%TIME_LOOP
          IF(ASSOCIATED(time_loop)) THEN
            problem=>control_loop%PROBLEM
            IF(ASSOCIATED(problem)) THEN
              NULLIFY(solvers)
              NULLIFY(solver)
              !Get the solver. 
              CALL control_loop_solvers_get(control_loop,solvers,err,error,*999)            
              CALL solvers_solver_get(solvers,1,solver,err,error,*999)
              !Loop over the equations sets associated with the solver
              solver_equations=>solver%SOLVER_EQUATIONS
              IF(ASSOCIATED(solver_equations)) THEN
                solver_mapping=>solver_equations%SOLVER_MAPPING
                IF(ASSOCIATED(solver_mapping)) THEN
                  DO equations_set_idx=1,solver_mapping%NUMBER_OF_EQUATIONS_SETS
                    equations_set=>solver_mapping%EQUATIONS_SETS(equations_set_idx)%PTR
                    IF(ASSOCIATED(equations_set)) THEN
                      dependent_field=>equations_set%DEPENDENT%DEPENDENT_FIELD
                      NULLIFY(dependent_region)
                      CALL field_region_get(dependent_field,dependent_region,err,error,*999)
                      NULLIFY(parent_loop)
                      parent_loop=>control_loop%PARENT_LOOP
                      IF(ASSOCIATED(parent_loop)) THEN
                        !add the iteration number of the parent loop to the filename
                        NULLIFY(time_loop_parent)
                        time_loop_parent=>parent_loop%TIME_LOOP
                        IF(ASSOCIATED(time_loop_parent)) THEN
                          output_iteration_number=time_loop_parent%OUTPUT_NUMBER
                          current_loop_iteration=time_loop_parent%GLOBAL_ITERATION_NUMBER
                          filename="Time_"//trim(number_to_vstring(dependent_region%USER_NUMBER,"*",err,error))// &
                            & "_"//trim(number_to_vstring(time_loop_parent%GLOBAL_ITERATION_NUMBER,"*",err,error))// &
                            & "_"//trim(number_to_vstring(time_loop%ITERATION_NUMBER,"*",err,error))
                        ELSE
                          output_iteration_number=time_loop%OUTPUT_NUMBER
                          current_loop_iteration=time_loop%GLOBAL_ITERATION_NUMBER
                          filename="Time_"//trim(number_to_vstring(dependent_region%USER_NUMBER,"*",err,error))// &
                            & "_"//trim(number_to_vstring(time_loop%GLOBAL_ITERATION_NUMBER,"*",err,error))
                        ENDIF
                      ELSE
                        output_iteration_number=time_loop%OUTPUT_NUMBER
                        current_loop_iteration=time_loop%GLOBAL_ITERATION_NUMBER
                        filename="Time_"//trim(number_to_vstring(dependent_region%USER_NUMBER,"*",err,error))// &
                          & "_"//trim(number_to_vstring(time_loop%GLOBAL_ITERATION_NUMBER,"*",err,error))
                      ENDIF
                      method="FORTRAN"
                      IF(output_iteration_number/=0.AND.mod(current_loop_iteration,output_iteration_number)==0) THEN
                        CALL field_io_nodes_export(dependent_region%FIELDS,filename,method,err,error,*999)
                      ENDIF
                    ELSE
                      local_error="Equations set is not associated for equations set index "// &
                        & trim(number_to_vstring(equations_set_idx,"*",err,error))// &
                        & " in the solver mapping."
                      CALL flagerror(local_error,err,error,*999)
                    ENDIF
                  ENDDO !equations_set_idx
                ELSE
                  CALL flagerror("Solver equations solver mapping is not associated.",err,error,*999)
                ENDIF
              ELSE
                CALL flagerror("Solver solver equations are not associated.",err,error,*999)
              ENDIF
            ELSE
              CALL flagerror("Control loop problem is not associated.",err,error,*999)
            ENDIF
          ELSE
            CALL flagerror("Time loop is not associated.",err,error,*999)
          ENDIF
        CASE(problem_control_while_loop_type)
          !do nothing
        CASE(problem_control_load_increment_loop_type)
          !do nothing
        CASE DEFAULT
          local_error="The control loop type of "//trim(number_to_vstring(control_loop%LOOP_TYPE,"*",err,error))// &
            & " is invalid."
          CALL flagerror(local_error,err,error,*999)
        END SELECT
      ENDIF
    ELSE
      CALL flagerror("Control loop is not associated.",err,error,*999)
    ENDIF

    exits("MONODOMAIN_CONTROL_LOOP_POST_LOOP")
    RETURN
999 errorsexits("MONODOMAIN_CONTROL_LOOP_POST_LOOP",err,error)
    RETURN 1
    
  END SUBROUTINE monodomain_control_loop_post_loop
  !
  !================================================================================================================================
  !
  !

  SUBROUTINE monodomain_equationssetspecificationset(equationsSet,specification,err,error,*)

    !Argument variables
    TYPE(equations_set_type), POINTER :: equationsSet
    INTEGER(INTG), INTENT(IN) :: specification(:)
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables
    INTEGER(INTG) :: subtype
    TYPE(varying_string) :: localError

    enters("Monodomain_EquationsSetSpecificationSet",err,error,*999)

    IF(ASSOCIATED(equationsset)) THEN
      IF(SIZE(specification,1)<3) THEN
        CALL flagerror("Equations set specification must have at least three entries for a monodomain class equations set.", &
          & err,error,*999)
      END IF
      SELECT CASE(specification(2))
      CASE(equations_set_monodomain_strang_splitting_equation_type)
        subtype=specification(3)
        SELECT CASE(subtype)
        CASE(equations_set_monodomain_buenoorovio_subtype, &
          & equations_set_monodomain_tentusscher06_subtype)
          !ok
        CASE DEFAULT
          localerror="Equations set subtype "//trim(numbertovstring(subtype,"*",err,error))// &
            & " is not valid for a Monodomain equation type of a Strang splitting equations set class."
          CALL flagerror(localerror,err,error,*999)
        END SELECT
        !Set full specification
        IF(ALLOCATED(equationsset%specification)) THEN
          CALL flagerror("Equations set specification is already allocated.",err,error,*999)
        ELSE
          ALLOCATE(equationsset%specification(3),stat=err)
          IF(err/=0) CALL flagerror("Could not allocate equations set specification.",err,error,*999)
        END IF
        equationsset%specification(1:3)=[equations_set_bioelectrics_class, &
          & equations_set_monodomain_strang_splitting_equation_type,subtype]
      CASE DEFAULT
        localerror="Equations set equation type "//trim(numbertovstring(specification(2),"*",err,error))// &
          & " is not valid for a monodomain equations set class."
      END SELECT
    ELSE
      CALL flagerror("Equations set is not associated",err,error,*999)
    END IF

    CALL exits("Monodomain_EquationsSetSpecificationSet")
    RETURN
999 CALL errors("Monodomain_EquationsSetSpecificationSet",err,error)
    CALL exits("Monodomain_EquationsSetSpecificationSet")
    RETURN 1
    
  END SUBROUTINE monodomain_equationssetspecificationset

  !
  !================================================================================================================================
  !

  SUBROUTINE monodomain_problemspecificationset(problem,problemSpecification,err,error,*)

    !Argument variables
    TYPE(problem_type), POINTER :: problem
    INTEGER(INTG), INTENT(IN) :: problemSpecification(:)
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables
    TYPE(varying_string) :: localError
    INTEGER(INTG) :: problemType,problemSubtype

    CALL enters("Monodomain_ProblemSpecificationSet",err,error,*999)

    IF(ASSOCIATED(problem)) THEN
      IF(SIZE(problemspecification,1)>=3) THEN
        problemtype=problemspecification(2)
        SELECT CASE(problemtype)
        CASE(problem_monodomain_strang_splitting_equation_type)
          problemsubtype=problemspecification(3)
          SELECT CASE(problemsubtype)
          CASE(problem_monodomain_buenoorovio_subtype, &
            & problem_monodomain_tentusscher06_subtype)
            !ok
          CASE DEFAULT
            localerror="Problem subtype "//trim(numbertovstring(problemsubtype,"*",err,error))// &
              & " is not valid for a Monodomain equation type of a Strang splitting  problem class."
            CALL flagerror(localerror,err,error,*999)
          END SELECT
          IF(ALLOCATED(problem%specification)) THEN
            CALL flagerror("Problem specification is already allocated.",err,error,*999)
          ELSE
            ALLOCATE(problem%specification(3),stat=err)
            IF(err/=0) CALL flagerror("Could not allocate problem specification.",err,error,*999)
          END IF
          problem%specification(1:3)=[problem_bioelectrics_class,problem_monodomain_strang_splitting_equation_type,problemsubtype]
        CASE DEFAULT
          localerror="Problem equation type "//trim(numbertovstring(problemtype,"*",err,error))// &
            & " is not valid for a monodomain problem class."
          CALL flagerror(localerror,err,error,*999)
        END SELECT
      ELSE
        CALL flagerror("Monodomain problem specification must have a type.",err,error,*999)
      END IF
    ELSE
      CALL flagerror("Problem is not associated",err,error,*999)
    END IF

    CALL exits("Monodomain_ProblemSpecificationSet")
    RETURN
999 CALL errors("Monodomain_ProblemSpecificationSet",err,error)
    CALL exits("Monodomain_ProblemSpecificationSet")
    RETURN 1
    
  END SUBROUTINE monodomain_problemspecificationset

  !
  !================================================================================================================================
  !

  SUBROUTINE monodomain_finiteelementcalculate(EQUATIONS_SET,ELEMENT_NUMBER,ERR,ERROR,*)

    !Argument variables
    TYPE(equations_set_type), POINTER :: EQUATIONS_SET
    INTEGER(INTG), INTENT(IN) :: ELEMENT_NUMBER
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) FIELD_VAR_TYPE,ng,mh,mhs,ms,nh,nhs,ni,ns,nj
    REAL(DP) :: RWG,SUM,Df, Dt, D(3,3), f(3), fnorm
    REAL(DP) :: DPHIDX(3,8) ! assumes <= 8 basis functions / DEPENDENT_BASIS%NUMBER_OF_ELEMENT_PARAMETERS <= 8
    TYPE(basis_type), POINTER :: DEPENDENT_BASIS,GEOMETRIC_BASIS
    TYPE(equations_type), POINTER :: EQUATIONS
    TYPE(equations_mapping_type), POINTER :: EQUATIONS_MAPPING
    TYPE(equations_mapping_dynamic_type), POINTER :: DYNAMIC_MAPPING
    TYPE(equations_matrices_type), POINTER :: EQUATIONS_MATRICES
    TYPE(equations_matrices_dynamic_type), POINTER :: DYNAMIC_MATRICES
    TYPE(equations_matrices_rhs_type), POINTER :: RHS_VECTOR
    TYPE(equations_matrix_type), POINTER :: DAMPING_MATRIX,STIFFNESS_MATRIX
    TYPE(field_type), POINTER :: DEPENDENT_FIELD,GEOMETRIC_FIELD,MATERIALS_FIELD
    TYPE(field_variable_type), POINTER :: FIELD_VARIABLE,GEOMETRIC_VARIABLE
    TYPE(quadrature_scheme_type), POINTER :: QUADRATURE_SCHEME
    TYPE(varying_string) :: LOCAL_ERROR
     
    enters("Monodomain_FiniteElementCalculate",err,error,*999)
    
    IF(ASSOCIATED(equations_set)) THEN
      equations=>equations_set%EQUATIONS
      IF(ASSOCIATED(equations)) THEN
        IF(.NOT.ALLOCATED(equations_set%SPECIFICATION)) THEN
          CALL flagerror("Equations set specification is not allocated.",err,error,*999)
        ELSE IF(SIZE(equations_set%SPECIFICATION,1)/=3) THEN
          CALL flagerror("Equations set specification must have three entries for a monodomain type equations set.", &
            & err,error,*999)
        END IF
        SELECT CASE(equations_set%SPECIFICATION(3))
        CASE(equations_set_monodomain_buenoorovio_subtype,equations_set_monodomain_tentusscher06_subtype)

          dependent_field=>equations%INTERPOLATION%DEPENDENT_FIELD
          geometric_field=>equations%INTERPOLATION%GEOMETRIC_FIELD
          materials_field=>equations%INTERPOLATION%MATERIALS_FIELD
          equations_matrices=>equations%EQUATIONS_MATRICES
          dynamic_matrices=>equations_matrices%DYNAMIC_MATRICES
          stiffness_matrix=>dynamic_matrices%MATRICES(1)%PTR
          damping_matrix=>dynamic_matrices%MATRICES(2)%PTR
          rhs_vector=>equations_matrices%RHS_VECTOR

          IF(.NOT.(stiffness_matrix%UPDATE_MATRIX.OR.damping_matrix%UPDATE_MATRIX.OR.rhs_vector%UPDATE_VECTOR)) RETURN ! no updates -> return immediately
 
          ! set up opencmiss interpolation and basis
          equations_mapping=>equations%EQUATIONS_MAPPING
          dynamic_mapping=>equations_mapping%DYNAMIC_MAPPING
          field_variable=>dynamic_mapping%EQUATIONS_MATRIX_TO_VAR_MAPS(1)%VARIABLE
          field_var_type=field_variable%VARIABLE_TYPE
          geometric_variable=>geometric_field%VARIABLE_TYPE_MAP(field_u_variable_type)%PTR
          dependent_basis=>dependent_field%DECOMPOSITION%DOMAIN(dependent_field%DECOMPOSITION%MESH_COMPONENT_NUMBER)%PTR% &
            & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
          geometric_basis=>geometric_field%DECOMPOSITION%DOMAIN(geometric_field%DECOMPOSITION%MESH_COMPONENT_NUMBER)%PTR% &
            & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
          quadrature_scheme=>dependent_basis%QUADRATURE%QUADRATURE_SCHEME_MAP(basis_default_quadrature_scheme)%PTR

          CALL field_interpolation_parameters_element_get(field_values_set_type,element_number,equations%INTERPOLATION% &
            & geometric_interp_parameters(field_u_variable_type)%PTR,err,error,*999)
          CALL field_interpolation_parameters_element_get(field_values_set_type,element_number,equations%INTERPOLATION% &
            & materials_interp_parameters(field_u_variable_type)%PTR,err,error,*999)


          DO ng=1,quadrature_scheme%NUMBER_OF_GAUSS
            ! get interpolated geometric and material interpolated point
            CALL field_interpolate_gauss(first_part_deriv,basis_default_quadrature_scheme,ng,equations%INTERPOLATION% &

              & geometric_interp_point(field_u_variable_type)%PTR,err,error,*999)
            CALL field_interpolated_point_metrics_calculate(geometric_basis%NUMBER_OF_XI,equations%INTERPOLATION% &
              & geometric_interp_point_metrics(field_u_variable_type)%PTR,err,error,*999)
            CALL field_interpolate_gauss(no_part_deriv,basis_default_quadrature_scheme,ng,equations%INTERPOLATION% &
              & materials_interp_point(field_u_variable_type)%PTR,err,error,*999)

            ! calculate weight = det J * gauss pt weight
            rwg=equations%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS(field_u_variable_type)%PTR%JACOBIAN* &
              & quadrature_scheme%GAUSS_WEIGHTS(ng)

            ! basis function chain rule taken out of the inner loop.
            DO ms=1,dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
            DO nj=1,geometric_variable%NUMBER_OF_COMPONENTS
              dphidx(nj,ms)=0.0_dp
              DO ni=1,dependent_basis%NUMBER_OF_XI                          
                dphidx(nj,ms)=dphidx(nj,ms) + &
                             & quadrature_scheme%GAUSS_BASIS_FNS(ms,partial_derivative_first_derivative_map(ni),ng)* &
                             & equations%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS(field_u_variable_type)%PTR%DXI_DX(ni,nj)
              ENDDO !ni
            ENDDO !nj
            ENDDO !ms

            ! Diffusion tensor  D  =  Dt I + (Df - Dt) f f^T  where Dt and Df are diffusivity/conductivity in fiber/transverse directions
            df = equations%INTERPOLATION%MATERIALS_INTERP_POINT(field_u_variable_type)%PTR%VALUES(2,no_part_deriv) ! 2 = Df
            dt = equations%INTERPOLATION%MATERIALS_INTERP_POINT(field_u_variable_type)%PTR%VALUES(3,no_part_deriv) ! 3 = Dt
            fnorm = 0.0
            DO nj=1,geometric_variable%NUMBER_OF_COMPONENTS
              f(nj) = equations%INTERPOLATION%MATERIALS_INTERP_POINT(field_u_variable_type)%PTR%VALUES(3+nj,no_part_deriv) ! 4,5[,6] = f
              fnorm = fnorm + f(nj)*f(nj)
            ENDDO
            ! normalize f, and fill in default for 0,0,0 -> 1,0,0
            fnorm = sqrt(fnorm)
            IF(fnorm < 1e-6) THEN
              f = (/ 1.0, 0.0, 0.0 /) ! default
            ELSE
              f = f / fnorm
            ENDIF
            DO ni=1,geometric_variable%NUMBER_OF_COMPONENTS
              d(ni,:)  = 0.0
              d(ni,ni) = dt
              DO nj=1,geometric_variable%NUMBER_OF_COMPONENTS
                d(ni,nj) = d(ni,nj) + (df - dt) * f(ni) * f(nj)
              ENDDO
            ENDDO

            !Loop over field components
            mhs=0          
            DO mh=1,field_variable%NUMBER_OF_COMPONENTS
              !Loop over element rows
              DO ms=1,dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
                mhs=mhs+1
                nhs=0
                IF(stiffness_matrix%UPDATE_MATRIX.OR.damping_matrix%UPDATE_MATRIX) THEN
                  !Loop over element columns. TODO: use symmetry?
                  DO nh=1,field_variable%NUMBER_OF_COMPONENTS
                    DO ns=1,dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
                      nhs=nhs+1
                      IF(stiffness_matrix%UPDATE_MATRIX) THEN
                        sum=0.0_dp
                        DO ni=1,geometric_variable%NUMBER_OF_COMPONENTS 
                        DO nj=1,geometric_variable%NUMBER_OF_COMPONENTS
                          sum=sum + d(ni,nj) * dphidx(ni,ms) * dphidx(nj,ns)
                        ENDDO !nj
                        ENDDO !ni
                        stiffness_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs)=stiffness_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs)+sum*rwg ! Aij = int D_ij * dphi_m/dx_i * dphi_n/dx_j 
                      ENDIF
 
                      IF(damping_matrix%UPDATE_MATRIX) THEN ! non mass lumped version
                         damping_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs)=damping_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs)+ &
                            & quadrature_scheme%GAUSS_BASIS_FNS(ms,no_part_deriv,ng)* &
                            & quadrature_scheme%GAUSS_BASIS_FNS(ns,no_part_deriv,ng)*rwg ! int phi_m phi_n
                      ENDIF 

                    ENDDO !ns
                  ENDDO !nh
                ENDIF
                IF(rhs_vector%UPDATE_VECTOR) rhs_vector%ELEMENT_VECTOR%VECTOR(mhs)=0.0_dp
!                IF(DAMPING_MATRIX%UPDATE_MATRIX) THEN  ! mass lumnped version
!                  DAMPING_MATRIX%ELEMENT_MATRIX%MATRIX(mhs,mhs)=DAMPING_MATRIX%ELEMENT_MATRIX%MATRIX(mhs,mhs)+ &
!                  & QUADRATURE_SCHEME%GAUSS_BASIS_FNS(ms,NO_PART_DERIV,ng)*RWG   !  // int phi_m
!                ENDIF
              ENDDO !ms
            ENDDO !mh
          IF(rhs_vector%UPDATE_VECTOR) rhs_vector%ELEMENT_VECTOR%VECTOR(mhs)=0.0_dp 
          ENDDO !ng
          IF(dependent_field%SCALINGS%SCALING_TYPE/=field_no_scaling) THEN
            CALL field_interpolationparametersscalefactorselementget(element_number,equations%INTERPOLATION% &
              & dependent_interp_parameters(field_var_type)%PTR,err,error,*999)
            mhs=0          
            DO mh=1,field_variable%NUMBER_OF_COMPONENTS
              !Loop over element rows
              DO ms=1,dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
                mhs=mhs+1                    
                nhs=0
                IF(stiffness_matrix%UPDATE_MATRIX.OR.damping_matrix%UPDATE_MATRIX) THEN
                  !Loop over element columns
                  DO nh=1,field_variable%NUMBER_OF_COMPONENTS
                    DO ns=1,dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
                      nhs=nhs+1
                      IF(stiffness_matrix%UPDATE_MATRIX) THEN
                        stiffness_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs)=stiffness_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs)* &
                          & equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_var_type)%PTR%SCALE_FACTORS(ms,mh)* &
                          & equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_var_type)%PTR%SCALE_FACTORS(ns,nh)
                      ENDIF
                      IF(damping_matrix%UPDATE_MATRIX) THEN ! non mass lumped version
                        damping_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs)=damping_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs)* &
                           & equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_var_type)%PTR%SCALE_FACTORS(ms,mh)* &
                           & equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_var_type)%PTR%SCALE_FACTORS(ns,nh)
                      ENDIF 
                    ENDDO !ns
                  ENDDO !nh
                ENDIF

                IF(rhs_vector%UPDATE_VECTOR) rhs_vector%ELEMENT_VECTOR%VECTOR(mhs)=rhs_vector%ELEMENT_VECTOR%VECTOR(mhs)* &
                  & equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_var_type)%PTR%SCALE_FACTORS(ms,mh) 
 !               IF(DAMPING_MATRIX%UPDATE_MATRIX) THEN ! mass lumped version
 !                  DAMPING_MATRIX%ELEMENT_MATRIX%MATRIX(mhs,nhs)=DAMPING_MATRIX%ELEMENT_MATRIX%MATRIX(mhs,nhs)* &
 !                    & EQUATIONS%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(FIELD_VAR_TYPE)%PTR%SCALE_FACTORS(ms,mh)
 !               ENDIF
              ENDDO !ms
            ENDDO !mh
          ENDIF
        CASE DEFAULT
          local_error="Equations set subtype "//trim(number_to_vstring(equations_set%SPECIFICATION(3),"*",err,error))// &
            & " is not valid for a Monodomain equation type of a Strang splitting equations set class."
          CALL flagerror(local_error,err,error,*999)
        END SELECT
      ELSE
        CALL flagerror("Equations set equations is not associated.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Equations set is not associated.",err,error,*999)
    ENDIF
 
    exits("Monodomain_FiniteElementCalculate")
    RETURN
999 errorsexits("Monodomain_FiniteElementCalculate",err,error)
    RETURN 1
    
  END SUBROUTINE monodomain_finiteelementcalculate

  !
  !================================================================================================================================
  !

  SUBROUTINE monodomain_equation_equations_set_setup(EQUATIONS_SET,EQUATIONS_SET_SETUP,ERR,ERROR,*)

    !Argument variables
    TYPE(equations_set_type), POINTER :: EQUATIONS_SET
    TYPE(equations_set_setup_type), INTENT(INOUT) :: EQUATIONS_SET_SETUP
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    TYPE(varying_string) :: LOCAL_ERROR
    
    enters("MONODOMAIN_EQUATION_EQUATIONS_SET_SETUP",err,error,*999)

    IF(ASSOCIATED(equations_set)) THEN
      IF(.NOT.ALLOCATED(equations_set%SPECIFICATION)) THEN
        CALL flagerror("Equations set specification is not allocated.",err,error,*999)
      ELSE IF(SIZE(equations_set%SPECIFICATION,1)/=3) THEN
        CALL flagerror("Equations set specification must have three entries for a monodomain type equations set.", &
          & err,error,*999)
      END IF
      SELECT CASE(equations_set%SPECIFICATION(3))
      CASE(equations_set_monodomain_buenoorovio_subtype,equations_set_monodomain_tentusscher06_subtype)
        CALL monodomain_equationssetsubtypesetup(equations_set,equations_set_setup,err,error,*999)        
      CASE DEFAULT
        local_error="Equations set subtype "//trim(number_to_vstring(equations_set%SPECIFICATION(3),"*",err,error))// &
          & " is not valid for a Monodomain equation type of a Strang splitting equation set class."
        CALL flagerror(local_error,err,error,*999)
      END SELECT
    ELSE
      CALL flagerror("Equations set is not associated.",err,error,*999)
    ENDIF
       
    exits("MONODOMAIN_EQUATION_EQUATIONS_SET_SETUP")
    RETURN
999 errorsexits("MONODOMAIN_EQUATION_EQUATIONS_SET_SETUP",err,error)
    RETURN 1
  END SUBROUTINE monodomain_equation_equations_set_setup

  !
  !================================================================================================================================
  !

  SUBROUTINE monodomain_equationssetsolutionmethodset(EQUATIONS_SET,SOLUTION_METHOD,ERR,ERROR,*)

    !Argument variables
    TYPE(equations_set_type), POINTER :: EQUATIONS_SET
    INTEGER(INTG), INTENT(IN) :: SOLUTION_METHOD
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    TYPE(varying_string) :: LOCAL_ERROR
    
    enters("MONODOMAIN_EQUATIONS_SET_SOLUTION_METHOD_SET",err,error,*999)
    
    IF(ASSOCIATED(equations_set)) THEN
      IF(.NOT.ALLOCATED(equations_set%SPECIFICATION)) THEN
        CALL flagerror("Equations set specification is not allocated.",err,error,*999)
      ELSE IF(SIZE(equations_set%SPECIFICATION,1)/=3) THEN
        CALL flagerror("Equations set specification must have three entries for a monodomain type equations set.", &
          & err,error,*999)
      END IF
      SELECT CASE(equations_set%SPECIFICATION(3))
      CASE(equations_set_monodomain_buenoorovio_subtype,equations_set_monodomain_tentusscher06_subtype)        
        SELECT CASE(solution_method)
        CASE(equations_set_fem_solution_method)
          equations_set%SOLUTION_METHOD=equations_set_fem_solution_method
        CASE(equations_set_bem_solution_method)
          CALL flagerror("Not implemented.",err,error,*999)
        CASE(equations_set_fd_solution_method)
          CALL flagerror("Not implemented.",err,error,*999)
        CASE(equations_set_fv_solution_method)
          CALL flagerror("Not implemented.",err,error,*999)
        CASE(equations_set_gfem_solution_method)
          CALL flagerror("Not implemented.",err,error,*999)
        CASE(equations_set_gfv_solution_method)
          CALL flagerror("Not implemented.",err,error,*999)
        CASE DEFAULT
          local_error="The specified solution method of "//trim(number_to_vstring(solution_method,"*",err,error))//" is invalid."
          CALL flagerror(local_error,err,error,*999)
        END SELECT
      CASE DEFAULT
        local_error="Equations set subtype of "//trim(number_to_vstring(equations_set%SPECIFICATION(3),"*",err,error))// &
          & " is not valid for a Monodomain equation type of monodomain Strang splitting equations set class."
        CALL flagerror(local_error,err,error,*999)
      END SELECT
    ELSE
      CALL flagerror("Equations set is not associated.",err,error,*999)
    ENDIF
       
    exits("Monodomain_EquationsSetSolutionMethodSet")
    RETURN
999 errors("Monodomain_EquationsSetSolutionMethodSet",err,error)
    exits("Monodomain_EquationsSetSolutionMethodSet")
    RETURN 1

  END SUBROUTINE monodomain_equationssetsolutionmethodset

  !
  !================================================================================================================================
  !

  SUBROUTINE monodomain_equationssetsubtypesetup(EQUATIONS_SET,EQUATIONS_SET_SETUP,ERR,ERROR,*)

    !Argument variables
    TYPE(equations_set_type), POINTER :: EQUATIONS_SET
    TYPE(equations_set_setup_type), INTENT(INOUT) :: EQUATIONS_SET_SETUP
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) :: component_idx,GEOMETRIC_MESH_COMPONENT,GEOMETRIC_SCALING_TYPE,NUMBER_OF_DIMENSIONS, &
      & NUMBER_OF_MATERIALS_COMPONENTS, NUM_COMP
    TYPE(decomposition_type), POINTER :: GEOMETRIC_DECOMPOSITION
    TYPE(equations_type), POINTER :: EQUATIONS
    TYPE(equations_mapping_type), POINTER :: EQUATIONS_MAPPING
    TYPE(equations_matrices_type), POINTER :: EQUATIONS_MATRICES
    TYPE(equations_set_materials_type), POINTER :: EQUATIONS_MATERIALS
    TYPE(varying_string) :: LOCAL_ERROR
    
    enters("MONODOMAIN_EQUATION_EQUATION_SET_SUBTYPE_SETUP",err,error,*999)
 
    NULLIFY(equations)
    NULLIFY(equations_mapping)
    NULLIFY(equations_matrices)
    NULLIFY(geometric_decomposition)


    IF(ASSOCIATED(equations_set)) THEN
      IF(.NOT.ALLOCATED(equations_set%SPECIFICATION)) THEN
        CALL flagerror("Equations set specification is not allocated.",err,error,*999)
      ELSE IF(SIZE(equations_set%SPECIFICATION,1)/=3) THEN
        CALL flagerror("Equations set specification must have three entries for a monodomain type equations set.", &
          & err,error,*999)
      END IF
      IF(equations_set%SPECIFICATION(3)==equations_set_monodomain_buenoorovio_subtype .OR. &
          & equations_set%SPECIFICATION(3)==equations_set_monodomain_tentusscher06_subtype) THEN
        SELECT CASE(equations_set_setup%SETUP_TYPE)
        CASE(equations_set_setup_initial_type)
          SELECT CASE(equations_set_setup%ACTION_TYPE)
          CASE(equations_set_setup_start_action)
            CALL monodomain_equationssetsolutionmethodset(equations_set,equations_set_fem_solution_method, &
              & err,error,*999)
          CASE(equations_set_setup_finish_action)
            !Do nothing
          CASE DEFAULT
            local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
              & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
              & " is invalid for a Monodomain equation."
            CALL flagerror(local_error,err,error,*999)
          END SELECT
        CASE(equations_set_setup_geometry_type)
          !Do nothing
        CASE(equations_set_setup_dependent_type)
          SELECT CASE(equations_set_setup%ACTION_TYPE)
          CASE(equations_set_setup_start_action)
            IF(equations_set%DEPENDENT%DEPENDENT_FIELD_AUTO_CREATED) THEN
              !Create the auto created dependent field
              CALL field_create_start(equations_set_setup%FIELD_USER_NUMBER,equations_set%REGION,equations_set%DEPENDENT% &
                & dependent_field,err,error,*999)
              CALL field_type_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_general_type,err,error,*999)
              CALL field_dependent_type_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_dependent_type,err,error,*999)
              CALL field_mesh_decomposition_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_decomposition,err,error,*999)
              CALL field_mesh_decomposition_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,geometric_decomposition, &
                & err,error,*999)
              CALL field_geometric_field_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,equations_set%GEOMETRY% &
                & geometric_field,err,error,*999)
              CALL field_number_of_variables_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,2,err,error,*999)
              CALL field_variable_types_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,(/field_u_variable_type, &
                & field_deludeln_variable_type/),err,error,*999)
              CALL field_dimension_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type, &
                & field_scalar_dimension_type,err,error,*999)
              CALL field_dimension_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_deludeln_variable_type, &
                & field_scalar_dimension_type,err,error,*999)
              CALL field_data_type_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type, &
                & field_dp_type,err,error,*999)
              CALL field_data_type_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_deludeln_variable_type, &
                & field_dp_type,err,error,*999) 

          

              CALL field_number_of_components_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type,1,&
                & err,error,*999)
              CALL field_number_of_components_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_deludeln_variable_type,&
                & 1 ,err,error,*999)

              !Default to the geometric interpolation setup
              CALL field_component_mesh_component_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type,1, &
                & geometric_mesh_component,err,error,*999)
              CALL field_component_mesh_component_set(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type,1, &
                & geometric_mesh_component,err,error,*999)
              CALL field_component_mesh_component_set(equations_set%DEPENDENT%DEPENDENT_FIELD,field_deludeln_variable_type,1, &
                & geometric_mesh_component,err,error,*999)
              SELECT CASE(equations_set%SOLUTION_METHOD)
              CASE(equations_set_fem_solution_method)
                CALL field_component_interpolation_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD, &
                  & field_u_variable_type,1,field_node_based_interpolation,err,error,*999)
                CALL field_component_interpolation_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD, &
                  & field_deludeln_variable_type,1,field_node_based_interpolation,err,error,*999)
                !Default the scaling to the geometric field scaling
                CALL field_scaling_type_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_scaling_type,err,error,*999)
                CALL field_scaling_type_set(equations_set%DEPENDENT%DEPENDENT_FIELD,geometric_scaling_type,err,error,*999)
              CASE DEFAULT
                local_error="The solution method of "//trim(number_to_vstring(equations_set%SOLUTION_METHOD,"*",err,error))// &
                  & " is invalid or not implemented"
                CALL flagerror(local_error,err,error,*999)
              END SELECT
            ELSE
              ! user specified field
              CALL flagerror("No user specified field supported!",err,error,*999)
            ENDIF
          CASE(equations_set_setup_finish_action)
            IF(equations_set%DEPENDENT%DEPENDENT_FIELD_AUTO_CREATED) THEN
              CALL field_create_finish(equations_set%DEPENDENT%DEPENDENT_FIELD,err,error,*999)
            ENDIF
          CASE DEFAULT
            local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
              & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
              & " is invalid for a Monodomain equation"
            CALL flagerror(local_error,err,error,*999)
          END SELECT

        CASE(equations_set_setup_independent_type)
          SELECT CASE(equations_set_setup%ACTION_TYPE)
          CASE(equations_set_setup_start_action)
            IF(equations_set%INDEPENDENT%INDEPENDENT_FIELD_AUTO_CREATED) THEN
              SELECT CASE(equations_set%SPECIFICATION(3))
              CASE(equations_set_monodomain_buenoorovio_subtype)
                 num_comp = 4
              CASE(equations_set_monodomain_tentusscher06_subtype)
                 num_comp = 19
              CASE DEFAULT
                CALL flagerror("Invalid cell model equations set subtype",err,error,*999)
              END SELECT    

              CALL field_create_start(equations_set_setup%FIELD_USER_NUMBER,equations_set%REGION,equations_set%INDEPENDENT% &
                & independent_field,err,error,*999)
              CALL field_type_set_and_lock(equations_set%INDEPENDENT%INDEPENDENT_FIELD,field_general_type,err,error,*999)
              CALL field_dependent_type_set_and_lock(equations_set%INDEPENDENT%INDEPENDENT_FIELD,&
                   & field_independent_type,err,error,*999)
              CALL field_mesh_decomposition_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_decomposition,err,error,*999)
              CALL field_mesh_decomposition_set_and_lock(equations_set%INDEPENDENT%INDEPENDENT_FIELD,geometric_decomposition, &
                & err,error,*999)
              CALL field_geometric_field_set_and_lock(equations_set%INDEPENDENT%INDEPENDENT_FIELD,equations_set%GEOMETRY% &
                & geometric_field,err,error,*999)
              CALL field_number_of_variables_set_and_lock(equations_set%INDEPENDENT%INDEPENDENT_FIELD,1,err,error,*999)
              CALL field_variable_types_set_and_lock(equations_set%INDEPENDENT%INDEPENDENT_FIELD,&
                   &(/field_u_variable_type/),err,error,*999)
              CALL field_dimension_set_and_lock(equations_set%INDEPENDENT%INDEPENDENT_FIELD,field_u_variable_type, &
                & field_vector_dimension_type,err,error,*999)
              CALL field_data_type_set_and_lock(equations_set%INDEPENDENT%INDEPENDENT_FIELD,field_u_variable_type, &
                & field_dp_type,err,error,*999)
         
              CALL field_number_of_components_set_and_lock(equations_set%INDEPENDENT%INDEPENDENT_FIELD,field_u_variable_type,&
                   &num_comp,err,error,*999)
              CALL field_dof_order_type_set(equations_set%INDEPENDENT%INDEPENDENT_FIELD,field_u_variable_type,&
              &  field_contiguous_component_dof_order,err,error,*999) ! dofs continuous, so first + (x-1) is x'th component index

              !Default to the geometric interpolation setup
              CALL field_component_mesh_component_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type,1, &
                & geometric_mesh_component,err,error,*999)
              CALL field_component_mesh_component_set(equations_set%INDEPENDENT%INDEPENDENT_FIELD,field_u_variable_type,1, &
                & geometric_mesh_component,err,error,*999)

              SELECT CASE(equations_set%SOLUTION_METHOD)
              CASE(equations_set_fem_solution_method)
                CALL field_component_interpolation_set_and_lock(equations_set%INDEPENDENT%INDEPENDENT_FIELD, &
                  & field_u_variable_type,1,field_node_based_interpolation,err,error,*999)
                !Default the scaling to the geometric field scaling
                CALL field_scaling_type_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_scaling_type,err,error,*999)
                CALL field_scaling_type_set(equations_set%INDEPENDENT%INDEPENDENT_FIELD,geometric_scaling_type,err,error,*999)
              END SELECT
            ELSE
              ! user specified field
              CALL flagerror("No user specified field supported!",err,error,*999)
            ENDIF

          CASE(equations_set_setup_finish_action)
            IF(equations_set%INDEPENDENT%INDEPENDENT_FIELD_AUTO_CREATED) THEN
              CALL field_create_finish(equations_set%INDEPENDENT%INDEPENDENT_FIELD,err,error,*999)
              SELECT CASE(equations_set%SPECIFICATION(3))
              CASE(equations_set_monodomain_buenoorovio_subtype)
                CALL bueno_orovio_initialize(equations_set%INDEPENDENT%INDEPENDENT_FIELD,err,error,*999) ! initialize to y0
              CASE(equations_set_monodomain_tentusscher06_subtype)
                CALL tentusscher06_initialize(equations_set%INDEPENDENT%INDEPENDENT_FIELD,err,error,*999) ! initialize to y0 
              CASE DEFAULT
                CALL flagerror("Invalid cell model equations set subtype",err,error,*999)
              END SELECT              
            ENDIF
          CASE DEFAULT
            local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
              & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
              & " is invalid for a Monodomain equation"
            CALL flagerror(local_error,err,error,*999)
          END SELECT

        CASE(equations_set_setup_materials_type)
          SELECT CASE(equations_set_setup%ACTION_TYPE)
           CASE(equations_set_setup_start_action)
            equations_materials=>equations_set%MATERIALS
             IF(ASSOCIATED(equations_materials)) THEN
               IF(equations_materials%MATERIALS_FIELD_AUTO_CREATED) THEN
                !Create the auto created materials field
                CALL field_create_start(equations_set_setup%FIELD_USER_NUMBER,equations_set%REGION,equations_materials% &
                  & materials_field,err,error,*999)
                CALL field_type_set_and_lock(equations_materials%MATERIALS_FIELD,field_material_type,err,error,*999)
                CALL field_dependent_type_set_and_lock(equations_materials%MATERIALS_FIELD,field_independent_type,err,error,*999)
                CALL field_mesh_decomposition_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_decomposition,err,error,*999)
                CALL field_mesh_decomposition_set_and_lock(equations_materials%MATERIALS_FIELD,geometric_decomposition, &
                  & err,error,*999)
                CALL field_geometric_field_set_and_lock(equations_materials%MATERIALS_FIELD,equations_set%GEOMETRY% &
                  & geometric_field,err,error,*999)
                CALL field_number_of_variables_set_and_lock(equations_materials%MATERIALS_FIELD,1,err,error,*999)
                CALL field_variable_types_set_and_lock(equations_materials%MATERIALS_FIELD,(/field_u_variable_type/), &
                  & err,error,*999)
                CALL field_dimension_set_and_lock(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
                  & field_vector_dimension_type,err,error,*999)
                CALL field_data_type_set_and_lock(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
                  & field_dp_type,err,error,*999)
                CALL field_number_of_components_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
                  & number_of_dimensions,err,error,*999)
                  !Materials field components are
                  ! 1. activation  factor (usually 0.0 or 1.0)
                  ! 2,3 for fiber/transverse conductivity   . defaults to constant interpolation 
                  ! 4,5[,6] : fiber unit vector in dimension
                  ! 7: out - activation times
                number_of_materials_components= 7 !NUMBER_OF_DIMENSIONS + 3
                 !Set the number of materials components
                CALL field_number_of_components_set_and_lock(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
                  & number_of_materials_components,err,error,*999)

                ! 1st = activation = node based, 2 3 diffusion constants
                CALL field_component_interpolation_set(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
                  & 1,field_node_based_interpolation,err,error,*999)
                CALL field_component_interpolation_set(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
                  & 2,field_constant_interpolation,err,error,*999)
                CALL field_component_interpolation_set(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
                  & 3,field_constant_interpolation,err,error,*999)
                ! 4 5 (6) fiber unit vector
                DO component_idx=1,3 !NUMBER_OF_DIMENSIONS
                  CALL field_component_mesh_component_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
                    & component_idx,geometric_mesh_component,err,error,*999)
                  CALL field_component_mesh_component_set(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
                    & component_idx+3,geometric_mesh_component,err,error,*999)
                  CALL field_component_interpolation_set(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
                    & component_idx+3,field_node_based_interpolation,err,error,*999)
                ENDDO !component_idx

                  CALL field_component_interpolation_set(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
                    & number_of_materials_components,field_node_based_interpolation,err,error,*999)

                  !Default the field scaling to that of the geometric field
                CALL field_scaling_type_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_scaling_type,err,error,*999)
                CALL field_scaling_type_set(equations_materials%MATERIALS_FIELD,geometric_scaling_type,err,error,*999)
              ELSE
                ! user specified field
                CALL flagerror("No user specified field supported!",err,error,*999)
              ENDIF
             ELSE

               CALL flagerror("Equations set materials is not associated.",err,error,*999)
             ENDIF
           CASE(equations_set_setup_finish_action)
            equations_materials=>equations_set%MATERIALS
            IF(ASSOCIATED(equations_materials)) THEN
              IF(equations_materials%MATERIALS_FIELD_AUTO_CREATED) THEN
                !Finish creating the materials field
                CALL field_create_finish(equations_materials%MATERIALS_FIELD,err,error,*999)
                !Set the default values for the materials field
                CALL field_number_of_components_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
                  & number_of_dimensions,err,error,*999)
                !Materials field components are 1 for each dimension 
                !i.e., k in div(k.grad(u(x)))
                number_of_materials_components=number_of_dimensions                             
                !First set the k values to 1.0
                DO component_idx=1,number_of_dimensions
                  CALL field_component_values_initialise(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
                    & field_values_set_type,component_idx,1.0_dp,err,error,*999)
                ENDDO !component_idx
              ENDIF
            ELSE
              CALL flagerror("Equations set materials is not associated.",err,error,*999)
            ENDIF
! ! ! Upto here
           CASE DEFAULT
             local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
               & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
               & " is invalid for a Monodomain equation."
             CALL flagerror(local_error,err,error,*999)
           END SELECT
        CASE(equations_set_setup_source_type)
          SELECT CASE(equations_set_setup%ACTION_TYPE)
          CASE(equations_set_setup_start_action)
            !Do nothing
          CASE(equations_set_setup_finish_action)
            !Do nothing
          CASE DEFAULT
            local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
              & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
              & " is invalid for a monodomain equation."
            CALL flagerror(local_error,err,error,*999)
          END SELECT
          CASE(equations_set_setup_equations_type)
          SELECT CASE(equations_set_setup%ACTION_TYPE)
          CASE(equations_set_setup_start_action)
            IF(equations_set%DEPENDENT%DEPENDENT_FINISHED) THEN
              CALL equations_create_start(equations_set,equations,err,error,*999)
              CALL equations_linearity_type_set(equations,equations_linear,err,error,*999)
              CALL equations_time_dependence_type_set(equations,equations_first_order_dynamic,err,error,*999)
            ELSE
              CALL flagerror("Equations set dependent field has not been finished.",err,error,*999)
            ENDIF
          CASE(equations_set_setup_finish_action)
            SELECT CASE(equations_set%SOLUTION_METHOD)
            CASE(equations_set_fem_solution_method)
              !Finish the equations creation
              CALL equations_set_equations_get(equations_set,equations,err,error,*999)
              CALL equations_create_finish(equations,err,error,*999)
              !Create the equations mapping.
              CALL equations_mapping_create_start(equations,equations_mapping,err,error,*999)
              !!! Check this 
              !CALL EQUATIONS_MAPPING_LINEAR_MATRICES_NUMBER_SET(EQUATIONS_MAPPING,1,ERR,ERROR,*999)
              CALL equations_mapping_dynamic_matrices_set(equations_mapping,.true.,.true.,err,error,*999)
              !CALL EQUATIONS_MAPPING_LINEAR_MATRICES_VARIABLE_TYPES_SET(EQUATIONS_MAPPING,(/FIELD_U_VARIABLE_TYPE/), &
               ! & ERR,ERROR,*999)
              CALL equations_mapping_dynamic_variable_type_set(equations_mapping,field_u_variable_type,err,error,*999)
              !!!! Check the above two lines
              CALL equations_mapping_rhs_variable_type_set(equations_mapping,field_deludeln_variable_type,err,error,*999)
              CALL equations_mapping_create_finish(equations_mapping,err,error,*999)
              !Create the equations matrices
              CALL equations_matrices_create_start(equations,equations_matrices,err,error,*999)
              SELECT CASE(equations%SPARSITY_TYPE)
              CASE(equations_matrices_full_matrices) 
                CALL equations_matrices_linear_storage_type_set(equations_matrices,(/matrix_block_storage_type/), &
                  & err,error,*999)
              CASE(equations_matrices_sparse_matrices) 
                !CALL EQUATIONS_MATRICES_DYNAMIC_STORAGE_TYPE_SET(EQUATIONS_MATRICES,(/MATRIX_COMPRESSED_ROW_STORAGE_TYPE/), &
                 ! & ERR,ERROR,*999)
                !CALL EquationsMatrices_DynamicStructureTypeSet(EQUATIONS_MATRICES,(/EQUATIONS_MATRIX_FEM_STRUCTURE/), &
                 ! & ERR,ERROR,*999)
                 CALL equations_matrices_dynamic_storage_type_set(equations_matrices, &
                    & (/distributed_matrix_compressed_row_storage_type,distributed_matrix_compressed_row_storage_type/), &
                    & err,error,*999)
                  CALL equationsmatrices_dynamicstructuretypeset(equations_matrices, &
                    (/equations_matrix_fem_structure,equations_matrix_fem_structure/),err,error,*999)     
              CASE DEFAULT
                local_error="The equations matrices sparsity type of "// &
                  & trim(number_to_vstring(equations%SPARSITY_TYPE,"*",err,error))//" is invalid."
                CALL flagerror(local_error,err,error,*999)
              END SELECT
              CALL equations_matrices_create_finish(equations_matrices,err,error,*999)
            CASE DEFAULT
                local_error="The solution method of "//trim(number_to_vstring(equations_set%SOLUTION_METHOD,"*",err,error))// &
                & " is invalid or not implemented."
              CALL flagerror(local_error,err,error,*999)
            END SELECT
          CASE DEFAULT
            local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
              & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
              & " is invalid for a Monodomain equation."
            CALL flagerror(local_error,err,error,*999)
          END SELECT
        CASE DEFAULT
          local_error="The setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
            & " is invalid for a Monodomain equation."
          CALL flagerror(local_error,err,error,*999)
        END SELECT
      ELSE
        local_error="The equations set subtype of "//trim(number_to_vstring(equations_set%SPECIFICATION(3),"*",err,error))// &
          & " does not equal a Monodomain equation subtype."
        CALL flagerror(local_error,err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Equations set is not associated.",err,error,*999)
    ENDIF
       
    exits("Monodomain_EquationsSetSubtypeSetup")
    RETURN
999 errorsexits("Monodomain_EquationsSetSubtypeSetup",err,error)
    RETURN 1
    
  END SUBROUTINE monodomain_equationssetsubtypesetup

  !
  !================================================================================================================================
  !
 
  SUBROUTINE monodomain_equation_problem_setup(PROBLEM,PROBLEM_SETUP,ERR,ERROR,*)
    
    !Argument variables
    TYPE(problem_type), POINTER :: PROBLEM
    TYPE(problem_setup_type), INTENT(INOUT) :: PROBLEM_SETUP
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    TYPE(varying_string) :: LOCAL_ERROR
    
    enters("MONODOMAIN_EQUATION_PROBLEM_SETUP",err,error,*999)

    IF(ASSOCIATED(problem)) THEN
      IF(.NOT.ALLOCATED(problem%SPECIFICATION)) THEN
        CALL flagerror("Problem specification is not allocated.",err,error,*999)
      ELSE IF(SIZE(problem%SPECIFICATION,1)<2) THEN
        CALL flagerror("Problem specification must have at least two entries for a monodomain equation problem.",err,error,*999)
      END IF
      SELECT CASE(problem%SPECIFICATION(2))
      CASE(problem_monodomain_strang_splitting_equation_type)
        CALL monodomain_problemstrangsplittingsetup(problem,problem_setup,err,error,*999)
      CASE DEFAULT
        local_error="Problem type "//trim(number_to_vstring(problem%SPECIFICATION(2),"*",err,error))// &
          & " is not valid for a Monodomain equation Strang splitting problem class."
        CALL flagerror(local_error,err,error,*999)
      END SELECT
    ELSE
      CALL flagerror("Problem is not associated.",err,error,*999)
    ENDIF
       
    exits("MONODOMAIN_EQUATION_PROBLEM_SETUP")
    RETURN
999 errorsexits("MONODOMAIN_EQUATION_PROBLEM_SETUP",err,error)
    RETURN 1
  END SUBROUTINE monodomain_equation_problem_setup
  
!
  !================================================================================================================================
  !
 
  SUBROUTINE monodomain_problemstrangsplittingsetup(PROBLEM,PROBLEM_SETUP,ERR,ERROR,*)
    
    !Argument variables
    TYPE(problem_type), POINTER :: PROBLEM
    TYPE(problem_setup_type), INTENT(INOUT) :: PROBLEM_SETUP
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    TYPE(varying_string) :: LOCAL_ERROR
    
    enters("Monodomain_ProblemStrangSplittingSetup",err,error,*999)

    IF(ASSOCIATED(problem)) THEN
      IF(.NOT.ALLOCATED(problem%SPECIFICATION)) THEN
        CALL flagerror("Problem specification is not allocated.",err,error,*999)
      ELSE IF(SIZE(problem%SPECIFICATION,1)<3) THEN
        CALL flagerror("Problem specification must have three entries for a monodomain Strang splitting problem.",err,error,*999)
      END IF
      SELECT CASE(problem%SPECIFICATION(3))
      CASE(problem_monodomain_buenoorovio_subtype,problem_monodomain_tentusscher06_subtype)
        CALL monodomain_equation_problem_subtype_setup(problem,problem_setup,err,error,*999)
      CASE DEFAULT
        local_error="Problem subtype "//trim(number_to_vstring(problem%SPECIFICATION(3),"*",err,error))// &
          & " is not valid for a Monodomain equation type of a Strang splitting problem class."
        CALL flagerror(local_error,err,error,*999)
      END SELECT
    ELSE
      CALL flagerror("Problem is not associated.",err,error,*999)
    ENDIF
       
    exits("Monodomain_ProblemStrangSplittingSetup")
    RETURN
999 errors("Monodomain_ProblemStrangSplittingSetup",err,error)
    exits("Monodomain_ProblemStrangSplittingSetup")
    RETURN 1
    
  END SUBROUTINE monodomain_problemstrangsplittingsetup
  
  !
  !================================================================================================================================
  !

  SUBROUTINE monodomain_equation_problem_subtype_setup(PROBLEM,PROBLEM_SETUP,ERR,ERROR,*)

    !Argument variables
    TYPE(problem_type), POINTER :: PROBLEM
    TYPE(problem_setup_type), INTENT(INOUT) :: PROBLEM_SETUP
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    TYPE(control_loop_type), POINTER :: CONTROL_LOOP,CONTROL_LOOP_ROOT
    TYPE(solver_type), POINTER :: SOLVER
    TYPE(solver_equations_type), POINTER :: SOLVER_EQUATIONS
    TYPE(solvers_type), POINTER :: SOLVERS
    TYPE(varying_string) :: LOCAL_ERROR
    
    enters("MONODOMAIN_EQUATION_PROBLEM_SUBTYPE_SETUP",err,error,*999)

    NULLIFY(control_loop)
    NULLIFY(solver)
    NULLIFY(solver_equations)
    NULLIFY(solvers)
    
    IF(ASSOCIATED(problem)) THEN
      IF(.NOT.ALLOCATED(control_loop%PROBLEM%SPECIFICATION)) THEN
        CALL flagerror("Problem specification is not allocated.",err,error,*999)
      ELSE IF(SIZE(control_loop%PROBLEM%SPECIFICATION,1)<3) THEN
        CALL flagerror("Problem specification must have three entries for a monodomain equation problem.",err,error,*999)
      END IF
      IF(problem%SPECIFICATION(3)==problem_monodomain_buenoorovio_subtype.OR. &
          & problem%SPECIFICATION(3)==problem_monodomain_tentusscher06_subtype)THEN
        SELECT CASE(problem_setup%SETUP_TYPE)
        CASE(problem_setup_initial_type)
          SELECT CASE(problem_setup%ACTION_TYPE)
          CASE(problem_setup_start_action)
            !Do nothing????
          CASE(problem_setup_finish_action)
            !Do nothing???
          CASE DEFAULT
            local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
              & " for a setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
              & " is invalid for a monodomain equation."
            CALL flagerror(local_error,err,error,*999)
          END SELECT
        CASE(problem_setup_control_type)
          SELECT CASE(problem_setup%ACTION_TYPE)
          CASE(problem_setup_start_action)
            !Set up a simple control loop
            CALL control_loop_create_start(problem,control_loop,err,error,*999)
            CALL control_loop_type_set(control_loop,problem_control_time_loop_type,err,error,*999)
          CASE(problem_setup_finish_action)
            !Finish the control loops
            control_loop_root=>problem%CONTROL_LOOP
            CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
            CALL control_loop_create_finish(control_loop,err,error,*999)            
          CASE DEFAULT
            local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
              & " for a setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
              & " is invalid for a monodomain equation."
            CALL flagerror(local_error,err,error,*999)
          END SELECT
        CASE(problem_setup_solvers_type)
           !Get the control loop
          control_loop_root=>problem%CONTROL_LOOP
          CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
          SELECT CASE(problem_setup%ACTION_TYPE)
          CASE(problem_setup_start_action)
            !Start the solvers creation
            CALL solvers_create_start(control_loop,solvers,err,error,*999)
            CALL solvers_number_set(solvers,1,err,error,*999)
            !Set the solver to be a first order dynamic solver
            CALL solvers_solver_get(solvers,1,solver,err,error,*999)
            CALL solver_type_set(solver,solver_dynamic_type,err,error,*999)
            CALL solver_dynamic_order_set(solver,solver_dynamic_first_order,err,error,*999)
            !Set solver defaults
            CALL solver_dynamic_degree_set(solver,solver_dynamic_first_degree,err,error,*999)
            CALL solver_dynamic_scheme_set(solver,solver_dynamic_crank_nicolson_scheme,err,error,*999)
            CALL solver_library_type_set(solver,solver_cmiss_library,err,error,*999)
          CASE(problem_setup_finish_action)
            !Get the solvers
            CALL control_loop_solvers_get(control_loop,solvers,err,error,*999)
            !Finish the solvers creation
            CALL solvers_create_finish(solvers,err,error,*999)
          CASE DEFAULT
            local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
              & " for a setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
                & " is invalid for a monodomain equation."
            CALL flagerror(local_error,err,error,*999)
          END SELECT
        CASE(problem_setup_solver_equations_type)
           SELECT CASE(problem_setup%ACTION_TYPE)
          CASE(problem_setup_start_action)
            !Get the control loop
            control_loop_root=>problem%CONTROL_LOOP
            CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
            !Get the solver
            CALL control_loop_solvers_get(control_loop,solvers,err,error,*999)
            CALL solvers_solver_get(solvers,1,solver,err,error,*999)
            !Create the solver equations
            CALL solver_equations_create_start(solver,solver_equations,err,error,*999)
            CALL solver_equations_linearity_type_set(solver_equations,solver_equations_linear,err,error,*999)
            CALL solver_equations_time_dependence_type_set(solver_equations,solver_equations_first_order_dynamic,err,error,*999)
            CALL solver_equations_sparsity_type_set(solver_equations,solver_sparse_matrices,err,error,*999)
          CASE(problem_setup_finish_action)
            !Get the control loop
            control_loop_root=>problem%CONTROL_LOOP
            CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
            !Get the solver equations
            CALL control_loop_solvers_get(control_loop,solvers,err,error,*999)
            CALL solvers_solver_get(solvers,1,solver,err,error,*999)
            CALL solver_solver_equations_get(solver,solver_equations,err,error,*999)
            !Finish the solver equations creation
            CALL solver_equations_create_finish(solver_equations,err,error,*999)             
          CASE DEFAULT
            local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
              & " for a setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
              & " is invalid for a monodomain equation."
            CALL flagerror(local_error,err,error,*999)
          END SELECT
        CASE DEFAULT
          local_error="The setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
            & " is invalid for a Monodomain equation."
          CALL flagerror(local_error,err,error,*999)
        END SELECT
      ELSE
        local_error="The problem subtype of "//trim(number_to_vstring(problem%SPECIFICATION(3),"*",err,error))// &
          & " does not equal a Monodomain equation subtype."
        CALL flagerror(local_error,err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Problem is not associated.",err,error,*999)
    ENDIF
    
    exits("MONODOMAIN_EQUATION_PROBLEM_SUBTYPE_SETUP")
    RETURN
999 errorsexits("MONODOMAIN_EQUATION_PROBLEM_SUBTYPE_SETUP",err,error)
    RETURN 1
  END SUBROUTINE monodomain_equation_problem_subtype_setup

  !
  !================================================================================================================================
  !

  SUBROUTINE monodomain_pre_solve(CONTROL_LOOP,SOLVER,ERR,ERROR,*)

    !Argument variables
    TYPE(control_loop_type), POINTER :: CONTROL_LOOP
    TYPE(solver_type), POINTER :: SOLVER
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    TYPE(varying_string) :: LOCAL_ERROR
    TYPE(field_type), POINTER :: DEPENDENT_FIELD, INDEPENDENT_FIELD

    enters("MONODOMAIN_PRE_SOLVE",err,error,*999)

    IF(ASSOCIATED(control_loop%PROBLEM)) THEN
      IF(.NOT.ALLOCATED(control_loop%PROBLEM%SPECIFICATION)) THEN
        CALL flagerror("Problem specification is not allocated.",err,error,*999)
      ELSE IF(SIZE(control_loop%PROBLEM%SPECIFICATION,1)<2) THEN
        CALL flagerror("Problem specification must at least two entries for a monodomain equation problem.",err,error,*999)
      END IF
      SELECT CASE(control_loop%PROBLEM%SPECIFICATION(2))
      CASE(problem_monodomain_strang_splitting_equation_type)
        dependent_field=>solver%SOLVERS%SOLVERS(1)%PTR%SOLVER_EQUATIONS%SOLVER_MAPPING% &
          & equations_sets(1)%PTR%DEPENDENT%DEPENDENT_FIELD
        independent_field=>solver%SOLVERS%SOLVERS(1)%PTR%SOLVER_EQUATIONS%SOLVER_MAPPING% &
          & equations_sets(1)%PTR%INDEPENDENT%INDEPENDENT_FIELD

        CALL field_parameterstofieldparameterscopy(independent_field,field_u_variable_type,field_values_set_type, &
          & 1,dependent_field,field_u_variable_type,field_values_set_type,1,err,error,*999)
        CALL field_parameterstofieldparameterscopy(independent_field,field_u_variable_type,field_values_set_type, &
          & 1,dependent_field,field_u_variable_type,field_previous_values_set_type,1,err,error,*999) ! also to prev.

      CASE DEFAULT
        local_error="Problem type "//trim(number_to_vstring(control_loop%PROBLEM%SPECIFICATION(2),"*",err,error))// &
          & " is not valid for a monodomain problem class."
        CALL flagerror(local_error,err,error,*999)
      END SELECT
    ELSE
      CALL flagerror("Problem is not associated.",err,error,*999)
    ENDIF
       
    exits("MONODOMAIN_PRE_SOLVE")
    RETURN
999 errorsexits("MONODOMAIN_PRE_SOLVE",err,error)
    RETURN 1
  END SUBROUTINE monodomain_pre_solve

  !
  !================================================================================================================================
  !

  SUBROUTINE monodomain_post_solve(CONTROL_LOOP,SOLVER,ERR,ERROR,*)

    !Argument variables
    TYPE(control_loop_type), POINTER :: CONTROL_LOOP
    TYPE(solver_type), POINTER :: SOLVER
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    TYPE(varying_string) :: LOCAL_ERROR
    REAL(DP) :: TMPV, TMPA
    INTEGER(INTG) :: I
    TYPE(field_type), POINTER :: DEPENDENT_FIELD, MATERIAL_FIELD,  INDEPENDENT_FIELD, GEOMETRIC_FIELD

    TYPE(equations_type), POINTER :: EQUATIONS
    TYPE(equations_matrices_type), POINTER :: EQUATIONS_MATRICES
    TYPE(equations_matrices_rhs_type), POINTER :: RHS_VECTOR
    TYPE(equations_matrix_type), POINTER :: DAMPING_MATRIX,STIFFNESS_MATRIX
    TYPE(equations_matrices_dynamic_type), POINTER :: DYNAMIC_MATRICES

    enters("MONODOMAIN_POST_SOLVE",err,error,*999)
  
    IF(ASSOCIATED(control_loop%PROBLEM)) THEN
      IF(.NOT.ALLOCATED(control_loop%PROBLEM%SPECIFICATION)) THEN
        CALL flagerror("Problem specification is not allocated.",err,error,*999)
      ELSE IF(SIZE(control_loop%PROBLEM%SPECIFICATION,1)<3) THEN
        CALL flagerror("Problem specification must have three entries for a monodomain equation problem.",err,error,*999)
      END IF
      SELECT CASE(control_loop%PROBLEM%SPECIFICATION(2))
      CASE(problem_monodomain_strang_splitting_equation_type) ! CAN NOT GET EQN TYPE?! 


       ! Don't update in next step
        equations=>solver%SOLVERS%SOLVERS(1)%PTR%SOLVER_EQUATIONS%SOLVER_MAPPING%EQUATIONS_SETS(1)%PTR%EQUATIONS
        equations_matrices=>equations%EQUATIONS_MATRICES
        dynamic_matrices=>equations_matrices%DYNAMIC_MATRICES
        stiffness_matrix=>dynamic_matrices%MATRICES(1)%PTR
        damping_matrix=>dynamic_matrices%MATRICES(2)%PTR
        rhs_vector=>equations_matrices%RHS_VECTOR
        stiffness_matrix%UPDATE_MATRIX = .false.
        damping_matrix%UPDATE_MATRIX   = .false.
        rhs_vector%UPDATE_VECTOR       = .false.

       ! integrate   cell models
        geometric_field => solver%SOLVERS%SOLVERS(1)%PTR%SOLVER_EQUATIONS%SOLVER_MAPPING% &
                           & equations_sets(1)%PTR%GEOMETRY%GEOMETRIC_FIELD
        dependent_field => solver%SOLVERS%SOLVERS(1)%PTR%SOLVER_EQUATIONS%SOLVER_MAPPING% &
                           & equations_sets(1)%PTR%DEPENDENT%DEPENDENT_FIELD
        material_field  => solver%SOLVERS%SOLVERS(1)%PTR%SOLVER_EQUATIONS%SOLVER_MAPPING% &
                           & equations_sets(1)%PTR%MATERIALS%MATERIALS_FIELD
        independent_field => solver%SOLVERS%SOLVERS(1)%PTR%SOLVER_EQUATIONS%SOLVER_MAPPING% &
                           & equations_sets(1)%PTR%INDEPENDENT%INDEPENDENT_FIELD

        CALL field_parameterstofieldparameterscopy(dependent_field,field_u_variable_type,field_values_set_type, &
         & 1, independent_field,field_u_variable_type,field_values_set_type, 1,err,error,*999) ! dependent -> independent

        SELECT CASE(control_loop%PROBLEM%SPECIFICATION(3))
        CASE(problem_monodomain_buenoorovio_subtype)
          CALL bueno_orovio_integrate(independent_field,material_field,&
          &  control_loop%TIME_LOOP%CURRENT_TIME-control_loop%TIME_LOOP%TIME_INCREMENT, control_loop%TIME_LOOP%CURRENT_TIME&  ! from t-dt to t
          & ,err,error,*999)
        CASE(problem_monodomain_tentusscher06_subtype)
           CALL tentusscher06_integrate(independent_field,material_field,&
          &  control_loop%TIME_LOOP%CURRENT_TIME-control_loop%TIME_LOOP%TIME_INCREMENT, control_loop%TIME_LOOP%CURRENT_TIME&  ! from t-dt to t
          & ,err,error,*999)
        CASE DEFAULT
          CALL flagerror("Invalid cell model subtype",err,error,*999)
        END SELECT

        DO i=1,independent_field%DECOMPOSITION%DOMAIN(1)%PTR%TOPOLOGY%NODES%NUMBER_OF_NODES
          !Default to version 1 of each derivative
          CALL field_parameter_set_get_node(independent_field,field_u_variable_type,field_values_set_type,1,1,i,1,tmpv,&  ! get local node?
               & err,error,*999)
          IF (tmpv > 0) THEN
           !Default to version 1 of each derivative
           CALL field_parameter_set_get_node(material_field,field_u_variable_type,field_values_set_type,1,1,i,7,tmpa,&
                & err,error,*999) 
            IF(abs(tmpa)<zero_tolerance) THEN
              !Default to version 1 of each derivative
              CALL field_parameter_set_update_node(material_field,field_u_variable_type,field_values_set_type,1,1,i,7,&
                   &control_loop%TIME_LOOP%CURRENT_TIME, err,error,*999)  
            ENDIF      
          ENDIF
        ENDDO

        IF(mod(control_loop%TIME_LOOP%CURRENT_TIME+1e-6,5.0)<1e-3) THEN
          WRITE(*,*)  'T=',control_loop%TIME_LOOP%CURRENT_TIME
        ENDIF

      CASE DEFAULT
        local_error="Problem type "//trim(number_to_vstring(control_loop%PROBLEM%SPECIFICATION(2),"*",err,error))// &
          & " is not valid for a monodomain problem class."
        CALL flagerror(local_error,err,error,*999)
      END SELECT
    ELSE
      CALL flagerror("Problem is not associated.",err,error,*999)
    ENDIF
       
    exits("MONODOMAIN_POST_SOLVE")
    RETURN
999 errorsexits("MONODOMAIN_POST_SOLVE",err,error)
    RETURN 1
  END SUBROUTINE monodomain_post_solve

  !
  !================================================================================================================================
  !

END MODULE monodomain_equations_routines