diff --git a/src/coordinate_routines.F90 b/src/coordinate_routines.F90
index aa09519b..a0f71893 100644
--- a/src/coordinate_routines.F90
+++ b/src/coordinate_routines.F90
@@ -3885,10 +3885,7 @@ SUBROUTINE Coordinates_MaterialSystemCalculate(geometricInterpPointMetrics,fibre
       ELSE
         !No fibre field
         numberOfNuDimensions=0
-        DO xiIdx=1,numberOfXiDimensions
-          dNudXiTemp(1:numberOfXDimensions,xiIdx)=geometricInterpPointMetrics%DX_DXI(1:numberOfXDimensions,xiIdx)
-          CALL Normalise(dNudXiTemp(1:numberOfXDimensions,xiIdx),dXdNu(1:numberOfXDimensions,xiIdx),err,error,*999)
-        ENDDO !xiIdx
+        CALL IdentityMatrix(dXdNU,err,error,*999)
       ENDIF
       !Calculate dNu/dX the inverse of dX/dNu (same as transpose due to orthogonality)
       CALL MatrixTranspose(dXdNu(1:numberOfXDimensions,1:numberOfXDimensions),dNudX(1:numberOfXDimensions,1: &
diff --git a/src/equations_routines.F90 b/src/equations_routines.F90
index 365c8e21..6b3a950e 100644
--- a/src/equations_routines.F90
+++ b/src/equations_routines.F90
@@ -660,7 +660,8 @@ SUBROUTINE Equations_InterpolationInitialise(equations,err,error,*)
       CALL Field_InterpolatedPointsInitialise(equations%interpolation%independentInterpParameters,  &
         & equations%interpolation%independentInterpPoint,err,error,*999)
       IF(equations%interpolation%independentField%type==FIELD_GEOMETRIC_TYPE.OR. &
-        & equations%interpolation%independentField%type==FIELD_FIBRE_TYPE) THEN
+         & equations%interpolation%independentField%type==FIELD_FIBRE_TYPE.OR. &
+         & equations%interpolation%independentField%type==FIELD_GEOMETRIC_GENERAL_TYPE) THEN
         CALL Field_InterpolatedPointsMetricsInitialise(equations%interpolation%independentInterpPoint,  &
           &  equations%interpolation%independentInterpPointMetrics,err,error,*999)
       END IF
diff --git a/src/equations_set_constants.F90 b/src/equations_set_constants.F90
index 53a9bbac..1f1c21e4 100644
--- a/src/equations_set_constants.F90
+++ b/src/equations_set_constants.F90
@@ -121,6 +121,8 @@ MODULE EquationsSetConstants
   INTEGER(INTG), PARAMETER :: EQUATIONS_SET_THREE_DIMENSIONAL_SUBTYPE=4
   INTEGER(INTG), PARAMETER :: EQUATIONS_SET_PLATE_SUBTYPE=5
   INTEGER(INTG), PARAMETER :: EQUATIONS_SET_SHELL_SUBTYPE=6
+  INTEGER(INTG), PARAMETER :: EQUATIONS_SET_ONE_DIM_STOKES_DAMPING_SUBTYPE=7
+  INTEGER(INTG), PARAMETER :: EQUATIONS_SET_TWO_DIM_PLANE_STRESS_STOKES_DAMPING_SUBTYPE=8
   !  Finite elasticity
   INTEGER(INTG), PARAMETER :: EQUATIONS_SET_MOONEY_RIVLIN_SUBTYPE=1
   INTEGER(INTG), PARAMETER :: EQUATIONS_SET_ISOTROPIC_EXPONENTIAL_SUBTYPE=2
diff --git a/src/finite_elasticity_routines.F90 b/src/finite_elasticity_routines.F90
index c5b8c450..6356612e 100644
--- a/src/finite_elasticity_routines.F90
+++ b/src/finite_elasticity_routines.F90
@@ -3757,7 +3757,7 @@ SUBROUTINE FiniteElasticity_StressStrainCalculate(equationsSet,derivedType,field
     CALL FieldVariable_FieldGet(fieldVariable,field,err,error,*999)
     fieldVarType=fieldVariable%VARIABLE_TYPE
    
-    CALL Field_NumberOfComponentsCheck(field,fieldVarType,6,err,error,*999)
+    CALL Field_NumberOfComponentsCheck(field,fieldVarType,numberOfComponents,err,error,*999)
     CALL Field_ComponentInterpolationGet(field,fieldVarType,1,fieldInterpolation,err,error,*999)
     !Check the interpolation type
     SELECT CASE(fieldInterpolation)
@@ -4260,7 +4260,8 @@ SUBROUTINE FiniteElasticity_StressStrainPoint(equationsSet,evaluateType,numberOf
       
       SELECT CASE(equationsSet%specification(3))
       CASE(EQUATIONS_SET_MOONEY_RIVLIN_SUBTYPE, &
-        & EQUATIONS_SET_MR_AND_GROWTH_LAW_IN_CELLML_SUBTYPE)
+        & EQUATIONS_SET_MR_AND_GROWTH_LAW_IN_CELLML_SUBTYPE, &
+        & EQUATIONS_SET_NEARLY_INCOMPRESSIBLE_MOONEY_RIVLIN_SUBTYPE)
         !Calculate the Cauchy stress tensor (in Voigt form) at the gauss point.
         Jznu=dependentInterpolatedPointMetrics%JACOBIAN/geometricInterpolatedPointMetrics%JACOBIAN
         ! Note that some problems, e.g. active contraction, require additonal fields to be evaluated at Gauss points. This is
@@ -6995,17 +6996,20 @@ SUBROUTINE FINITE_ELASTICITY_GAUSS_STRESS_TENSOR(EQUATIONS_SET,DEPENDENT_INTERPO
     REAL(DP), INTENT(IN) :: Jznu !Determinant of deformation gradient tensor (AZL)
     INTEGER(INTG), INTENT(IN) :: ELEMENT_NUMBER,GAUSS_POINT_NUMBER !<Element/Gauss point number
     INTEGER(INTG), INTENT(OUT) :: err !<The error code
+    INTEGER(INTG) :: i
     TYPE(VARYING_STRING), INTENT(OUT) :: error !<The error string
     !Local Variables
     INTEGER(INTG) :: PRESSURE_COMPONENT,component_idx,dof_idx
     REAL(DP) :: P
-    REAL(DP) :: I1 !Invariants, if needed
+    REAL(DP) :: I1,I2,I3 !Invariants, if needed
     REAL(DP) :: TEMPTERM1,TEMPTERM2,VALUE !Temporary variables
     REAL(DP) :: ONETHIRD_TRACE
     TYPE(VARYING_STRING) :: LOCAL_ERROR
     TYPE(FIELD_VARIABLE_TYPE), POINTER :: FIELD_VARIABLE
-    REAL(DP) :: MOD_DZDNU(3,3),MOD_DZDNUT(3,3),AZL(3,3)
+    REAL(DP) :: MOD_DZDNU(3,3),MOD_DZDNUT(3,3),AZL(3,3),AZU(3,3),PIOLA_TENSOR(3,3),IDENTITY(3,3)
+    REAL(DP) :: AZL_SQUARED(3,3)
     REAL(DP) :: B(6),E(6),DQ_DE(6)
+    REAL(DP) :: WV_PRIME
     REAL(DP), POINTER :: C(:) !Parameters for constitutive laws
 
     ENTERS("FINITE_ELASTICITY_GAUSS_STRESS_TENSOR",err,error,*999)
@@ -7021,6 +7025,38 @@ SUBROUTINE FINITE_ELASTICITY_GAUSS_STRESS_TENSOR(EQUATIONS_SET,DEPENDENT_INTERPO
     C=>MATERIALS_INTERPOLATED_POINT%VALUES(:,NO_PART_DERIV)
 
     SELECT CASE(EQUATIONS_SET%specification(3))
+    CASE(EQUATIONS_SET_NEARLY_INCOMPRESSIBLE_MOONEY_RIVLIN_SUBTYPE)
+
+       CALL INVERT(AZL,AZU,I3,ERR,ERROR,*999)
+       
+       IDENTITY=0.0_DP
+       DO i=1,3
+          IDENTITY(i,i)=1.0_DP
+       ENDDO
+
+      !Form of constitutive model is:
+      ! W_hat=c1*(I1_hat-3)+c2*(I2_hat-3)+p*J*C^(-1) + W^v(J)
+      ! take W^v(J) = 1/2 * kappa * (J-1)^2
+      WV_PRIME = C(3)*(Jznu - 1.0_DP)
+      !compute the invariants, I3 a few lines up
+      I1 = AZL(1,1) + AZL(2,2) + AZL(3,3)
+      CALL MatrixProduct(AZL,AZL,AZL_SQUARED,err,error,*999)
+      I2 = 0.5_DP * (I1**2 - AZL_SQUARED(1,1) - AZL_SQUARED(2,2) - AZL_SQUARED(3,3))
+
+      PIOLA_TENSOR=2.0_DP*Jznu**(-2.0_DP/3.0_DP)*((C(1)+C(2)*I1)*IDENTITY-C(2)*AZL &
+        & -(C(1)*I1+2.0_DP*C(2)*I2-1.5_DP*WV_PRIME*Jznu**(5.0_DP/3.0_DP))/3.0_DP*AZU)
+
+      !Calculate isochoric fictitious 2nd Piola tensor (in Voigt form)
+      STRESS_TENSOR(1)=PIOLA_TENSOR(1,1)
+      STRESS_TENSOR(2)=PIOLA_TENSOR(2,2)
+      STRESS_TENSOR(3)=PIOLA_TENSOR(3,3)
+      STRESS_TENSOR(4)=PIOLA_TENSOR(2,1)
+      STRESS_TENSOR(5)=PIOLA_TENSOR(3,1)
+      STRESS_TENSOR(6)=PIOLA_TENSOR(3,2)
+
+      !Do push-forward of 2nd Piola tensor. 
+      CALL FINITE_ELASTICITY_PUSH_STRESS_TENSOR(STRESS_TENSOR,MOD_DZDNU,Jznu,err,error,*999)
+
     CASE(EQUATIONS_SET_MOONEY_RIVLIN_ACTIVECONTRACTION_SUBTYPE, &
       & EQUATIONS_SET_MOONEY_RIVLIN_SUBTYPE, &
       & EQUATIONS_SET_MR_AND_GROWTH_LAW_IN_CELLML_SUBTYPE)
diff --git a/src/linear_elasticity_routines.F90 b/src/linear_elasticity_routines.F90
index 8bdd42a1..2dc14bac 100644
--- a/src/linear_elasticity_routines.F90
+++ b/src/linear_elasticity_routines.F90
@@ -1116,24 +1116,27 @@ SUBROUTINE LINEAR_ELASTICITY_FINITE_ELEMENT_CALCULATE(EQUATIONS_SET,ELEMENT_NUMB
     TYPE(VARYING_STRING), INTENT(OUT) :: ERROR !<The error string
 
     !Local Variables
-    INTEGER(INTG) :: FIELD_VAR_TYPE,ng,ni,xi,mi,ns,nhs,ms,mhs,TOTAL_DEPENDENT_BASIS_EP,MESH_COMPONENT
+    INTEGER(INTG) :: FIELD_VAR_TYPE,ng,ni,xi,mi,ns,nhs,ms,mhs,mh,nh,TOTAL_DEPENDENT_BASIS_EP,MESH_COMPONENT
     INTEGER(INTG) :: NUMBER_OF_DEPENDENT_COMPONENTS,NUMBER_OF_XI !,NUMBER_OF_GEOMETRIC_COMPONENTS
     INTEGER(INTG) :: OFF_DIAG_COMP(3),OFF_DIAG_DEP_VAR(2,2,3),DIAG_SUB_MAT_LOC(3),OFF_DIAG_SUB_MAT_LOC(2,3)
     INTEGER(INTG) :: DEPENDENT_BASES_EP(3) !,GEOMETRIC_BASES_EP(:)
     REAL(DP) :: JRWG,C(6,6),JRWG_DIAG_C(3,3),JRWG_OFF_DIAG_C(2,3)
-    REAL(DP):: SF(64*3)
+    REAL(DP) :: SF(64*3)
 
     !LOGICAL :: SAME_BASIS
     TYPE(EquationsType), POINTER :: equations
     TYPE(EquationsMappingVectorType), POINTER :: vectorMapping
+    TYPE(EquationsMappingDynamicType), POINTER :: dynamicMapping
     TYPE(EquationsMappingLinearType), POINTER :: linearMapping
     TYPE(EquationsMatricesVectorType), POINTER :: vectorMatrices
+    TYPE(EquationsMatricesDynamicType), POINTER :: dynamicMatrices
+    TYPE(EquationsMatricesRHSType), POINTER :: rhsVector
+    TYPE(EquationsMatrixType), POINTER :: dampingMatrix,stiffnessMatrix
     TYPE(EquationsMatrixType), POINTER :: equationsMatrix
     TYPE(EquationsVectorType), POINTER :: vectorEquations
     TYPE(FIELD_TYPE), POINTER :: dependentField,geometricField
     TYPE(BASIS_PTR_TYPE) :: DEPENDENT_BASES(3) !,GEOMETRIC_BASES(:)
     TYPE(QUADRATURE_SCHEME_PTR_TYPE) :: QUADRATURE_SCHEMES(3)
-    TYPE(EquationsMatricesRHSType), POINTER :: rhsVector
     TYPE(FIELD_INTERPOLATION_PARAMETERS_TYPE), POINTER :: GEOMETRIC_INTERPOLATION_PARAMETERS
     TYPE(FIELD_INTERPOLATION_PARAMETERS_TYPE), POINTER :: DEPENDENT_INTERPOLATION_PARAMETERS
     !TYPE(FIELD_INTERPOLATION_PARAMETERS_TYPE), POINTER :: FIBRE_INTERPOLATION_PARAMETERS
@@ -1170,11 +1173,23 @@ SUBROUTINE LINEAR_ELASTICITY_FINITE_ELEMENT_CALCULATE(EQUATIONS_SET,ELEMENT_NUMB
         geometricField=>equations%interpolation%geometricField
         !fibreField=>equations%interpolation%fibreField
         vectorMatrices=>vectorEquations%vectorMatrices
-        equationsMatrix=>vectorMatrices%linearMatrices%matrices(1)%ptr
         rhsVector=>vectorMatrices%rhsVector
         vectorMapping=>vectorEquations%vectorMapping
-        linearMapping=>vectorMapping%linearMapping
-        FIELD_VARIABLE=>linearMapping%equationsMatrixToVarMaps(1)%variable
+        SELECT CASE(EQUATIONS_SET%SPECIFICATION(3))
+        CASE(EQUATIONS_SET_THREE_DIMENSIONAL_SUBTYPE,EQUATIONS_SET_TWO_DIMENSIONAL_PLANE_STRESS_SUBTYPE, &
+             & EQUATIONS_SET_TWO_DIMENSIONAL_PLANE_STRAIN_SUBTYPE,EQUATIONS_SET_ONE_DIMENSIONAL_SUBTYPE)
+           linearMapping=>vectorMapping%linearMapping
+           equationsMatrix=>vectorMatrices%linearMatrices%matrices(1)%ptr
+           FIELD_VARIABLE=>linearMapping%equationsMatrixToVarMaps(1)%variable
+        CASE(EQUATIONS_SET_ONE_DIM_STOKES_DAMPING_SUBTYPE, &
+             & EQUATIONS_SET_TWO_DIM_PLANE_STRESS_STOKES_DAMPING_SUBTYPE)
+           dynamicMatrices=>vectorMatrices%dynamicMatrices
+           dynamicMapping=>vectorMapping%dynamicMapping
+           stiffnessMatrix=>dynamicMatrices%matrices(1)%ptr
+           dampingMatrix=>dynamicMatrices%matrices(2)%ptr
+           FIELD_VARIABLE=>dynamicMapping%equationsMatrixToVarMaps(1)%variable
+        END SELECT
+
         FIELD_VAR_TYPE=FIELD_VARIABLE%VARIABLE_TYPE
         !Note that the dimension/number of components for FIELD_U_VARIABLE_TYPE has to be the same as the FIELD_DELUDELN_VARIABLE_TYPE
         !& Number of Geometric field components = number of dependent field component
@@ -1378,8 +1393,176 @@ SUBROUTINE LINEAR_ELASTICITY_FINITE_ELEMENT_CALCULATE(EQUATIONS_SET,ELEMENT_NUMB
                 equationsMatrix%elementMatrix%matrix(mhs,nhs) = equationsMatrix%elementMatrix%matrix(nhs,mhs)
               ENDDO !nhs
             ENDDO !mhs
+         ENDIF
+         
+      CASE(EQUATIONS_SET_ONE_DIM_STOKES_DAMPING_SUBTYPE, &
+           & EQUATIONS_SET_TWO_DIM_PLANE_STRESS_STOKES_DAMPING_SUBTYPE)
+         
+          !//
+          !Parameters for number of off diagonal stress/strain terms for a given number of xi directions and order of calculation for shear terms
+          !These parameters do not change for 1D,2D,3D Linear Elasticity
+          OFF_DIAG_COMP = [0,1,3]
+          OFF_DIAG_DEP_VAR(1,1,:) = [1,1,2]
+          OFF_DIAG_DEP_VAR(1,2,:) = [2,3,3]
+          OFF_DIAG_DEP_VAR(2,1,:) = OFF_DIAG_DEP_VAR(1,2,:)
+          OFF_DIAG_DEP_VAR(2,2,:) = OFF_DIAG_DEP_VAR(1,1,:)
+          !//
+          DIAG_SUB_MAT_LOC(:) = [0,DEPENDENT_BASES_EP(1),SUM(DEPENDENT_BASES_EP(1:2))]
+          OFF_DIAG_SUB_MAT_LOC(1,:) = [0,0,DEPENDENT_BASES_EP(1)]
+          OFF_DIAG_SUB_MAT_LOC(2,:) = [DEPENDENT_BASES_EP(1),DIAG_SUB_MAT_LOC(3),DIAG_SUB_MAT_LOC(3)]
+
+          
+          IF(stiffnessMatrix%updateMatrix.OR.dampingMatrix%updateMatrix.OR.rhsVector%updateVector) THEN
+            !Get the geometric, fibre & material field interpolation parameters
+            GEOMETRIC_INTERPOLATION_PARAMETERS=>equations%interpolation%geometricInterpParameters(FIELD_U_VARIABLE_TYPE)%ptr
+            !FIBRE_INTERPOLATION_PARAMETERS=>equations%interpolation%FIBRE_INTERP_PARAMETERS
+            MATERIALS_INTERPOLATION_PARAMETERS=>equations%interpolation%materialsInterpParameters(FIELD_U_VARIABLE_TYPE)%ptr
+            CALL FIELD_INTERPOLATION_PARAMETERS_ELEMENT_GET(FIELD_VALUES_SET_TYPE,ELEMENT_NUMBER, &
+              & GEOMETRIC_INTERPOLATION_PARAMETERS,err,error,*999)
+            !CALL FIELD_INTERPOLATION_PARAMETERS_ELEMENT_GET(FIELD_VALUES_SET_TYPE,ELEMENT_NUMBER, &
+            !  & FIBRE_INTERPOLATION_PARAMETERS,err,error,*999)
+            CALL FIELD_INTERPOLATION_PARAMETERS_ELEMENT_GET(FIELD_VALUES_SET_TYPE,ELEMENT_NUMBER, &
+              & MATERIALS_INTERPOLATION_PARAMETERS,err,error,*999)
+            geometricInterpPointMetrics=>equations%interpolation%geometricInterpPointMetrics(FIELD_U_VARIABLE_TYPE)%ptr
+            materialsInterpPoint=>equations%interpolation%materialsInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr
+            !Loop over gauss points & integrate upper triangular portion of Stiffness matrix
+            DO ng=1,QUADRATURE_SCHEMES(1)%ptr%NUMBER_OF_GAUSS !Gauss point index
+              !Interpolate geometric, fibre and material fields at gauss points & calculate geometric field metrics
+              CALL FIELD_INTERPOLATE_GAUSS(FIRST_PART_DERIV,BASIS_DEFAULT_QUADRATURE_SCHEME,ng,equations%interpolation% &
+                & geometricInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr,err,error,*999)
+              !CALL FIELD_INTERPOLATE_GAUSS(FIRST_PART_DERIV,BASIS_DEFAULT_QUADRATURE_SCHEME,ng,equations%interpolation% &
+              !  & FIBRE_INTERP_POINT,err,error,*999)
+              CALL FIELD_INTERPOLATE_GAUSS(NO_PART_DERIV,BASIS_DEFAULT_QUADRATURE_SCHEME,ng,materialsInterpPoint,err,error,*999)
+              !!TODO:: Add option to only evaluate required metrics
+              CALL FIELD_INTERPOLATED_POINT_METRICS_CALCULATE(NUMBER_OF_XI,geometricInterpPointMetrics,err,error,*999)
+              !Calculate JRWG.
+              JRWG=QUADRATURE_SCHEMES(1)%ptr%GAUSS_WEIGHTS(ng)*geometricInterpPointMetrics%JACOBIAN
+              DO xi=1,NUMBER_OF_XI
+                DPHIDX_COMP(xi)%DPHIDX = 0.0_DP
+                DO ns=1,DEPENDENT_BASES_EP(xi)
+                  DO mi=1,NUMBER_OF_XI
+                    DO ni=1,NUMBER_OF_XI
+                      DPHIDX_COMP(xi)%DPHIDX(ns,mi) = DPHIDX_COMP(xi)%DPHIDX(ns,mi)+geometricInterpPointMetrics%DXI_DX(mi,ni)* &
+                        & QUADRATURE_SCHEMES(xi)%ptr%GAUSS_BASIS_FNS(ns,PARTIAL_DERIVATIVE_FIRST_DERIVATIVE_MAP(ni),ng)
+                    ENDDO !ni
+                  ENDDO !mi
+                ENDDO !ns
+              ENDDO !xi
+              !CALL COORDINATE_MATERIAL_COORDINATE_SYSTEM_CALCULATE(equations%interpolation%geometricInterpPoint, &
+              !  equations%interpolation%FIBRE_INTERP_POINT,DXDNU,err,error,*999)
+              !Create Linear Elasticity Tensor C
+              CALL LINEAR_ELASTICITY_TENSOR(EQUATIONS_SET%SPECIFICATION(3),materialsInterpPoint,C,err,error,*999)
+              !Store Elasticity Tensor diagonal & off diagonal stress coefficients
+              JRWG_DIAG_C(3,:) = [JRWG*C(4,4),JRWG*C(5,5),JRWG*C(3,3)]
+              JRWG_DIAG_C(2,:) = [JRWG*C(6,6),JRWG*C(2,2),JRWG_DIAG_C(3,2)]
+              JRWG_DIAG_C(1,:) = [JRWG*C(1,1),JRWG_DIAG_C(2,1),JRWG_DIAG_C(3,1)]
+              JRWG_OFF_DIAG_C(1,:) = [JRWG*C(1,2),JRWG*C(1,3),JRWG*C(2,3)]
+              JRWG_OFF_DIAG_C(2,:) = [JRWG*C(6,6),JRWG*C(4,4),JRWG*C(5,5)]             
+
+              !Construct Stiffness Matrix 
+              IF(stiffnessMatrix%updateMatrix) THEN
+
+                 !Construct Stiffness Matrix Diagonal Terms
+                 DO xi=1,NUMBER_OF_XI
+                    DO ns=1,DEPENDENT_BASES_EP(xi)
+                       DO ms=ns,DEPENDENT_BASES_EP(xi)
+                          stiffnessMatrix%elementMatrix%matrix(DIAG_SUB_MAT_LOC(xi)+ns,DIAG_SUB_MAT_LOC(xi)+ms) = &
+                               & stiffnessMatrix%elementMatrix%matrix(DIAG_SUB_MAT_LOC(xi)+ns,DIAG_SUB_MAT_LOC(xi)+ms) + &
+                               & DOT_PRODUCT(DPHIDX_COMP(xi)%DPHIDX(ns,1:NUMBER_OF_XI)*DPHIDX_COMP(xi)%DPHIDX(ms,1:NUMBER_OF_XI), &
+                               & JRWG_DIAG_C(xi,1:NUMBER_OF_XI)) 
+                       ENDDO !ms
+                    ENDDO !ns
+                 ENDDO !xi
+                 
+                 !Construct Stiffness Matrix Off-Diagonal Terms
+                 DO xi=1,OFF_DIAG_COMP(NUMBER_OF_XI)
+                    DO ns=1,DEPENDENT_BASES_EP(OFF_DIAG_DEP_VAR(1,1,xi))
+                       DO ms=1,DEPENDENT_BASES_EP(OFF_DIAG_DEP_VAR(1,2,xi))
+                          stiffnessMatrix%elementMatrix%matrix(OFF_DIAG_SUB_MAT_LOC(1,xi)+ns,OFF_DIAG_SUB_MAT_LOC(2,xi)+ms) = &
+                               & stiffnessMatrix%elementMatrix%matrix(OFF_DIAG_SUB_MAT_LOC(1,xi)+ns,OFF_DIAG_SUB_MAT_LOC(2,xi)+ms)& 
+                               & +DOT_PRODUCT(DPHIDX_COMP(OFF_DIAG_DEP_VAR(1,1,xi))%DPHIDX(ns,OFF_DIAG_DEP_VAR(1,:,xi))* &
+                               &  DPHIDX_COMP(OFF_DIAG_DEP_VAR(1,2,xi))%DPHIDX(ms,OFF_DIAG_DEP_VAR(2,:,xi)),JRWG_OFF_DIAG_C(:,xi))
+                       ENDDO !ms
+                    ENDDO !ns
+                 ENDDO !xi
+              ENDIF            
+            
+              !Construct Damping Matrix
+              IF(dampingMatrix%updateMatrix) THEN
+                 !Construct Damping Matrix
+                 mhs=0          
+                 DO mh=1,FIELD_VARIABLE%NUMBER_OF_COMPONENTS
+                    !Loop over element rows
+                    DO ms=1,DEPENDENT_BASES_EP(mh)
+                       mhs=mhs+1
+                       nhs=0
+                       !Loop over element columsn
+                       DO nh=1,FIELD_VARIABLE%NUMBER_OF_COMPONENTS
+                          DO ns=1,DEPENDENT_BASES_EP(nh)
+                             nhs=nhs+1
+                             dampingMatrix%elementMatrix%matrix(mhs,nhs)=dampingMatrix%elementMatrix%matrix(mhs,nhs)+ &
+                                  & QUADRATURE_SCHEMES(mh)%ptr%GAUSS_BASIS_FNS(ms,NO_PART_DERIV,ng)* &
+                                  & QUADRATURE_SCHEMES(nh)%ptr%GAUSS_BASIS_FNS(ns,NO_PART_DERIV,ng)* &
+                                  & materialsInterpPoint%values(2+FIELD_VARIABLE%NUMBER_OF_COMPONENTS,1)*JRWG
+                          ENDDO !ns
+                       ENDDO !nh
+                       !IF(rhsVector%updateVector) rhsVector%elementVector%vector(mhs)=0.0_DP
+                    ENDDO !ms
+                 ENDDO !mh
+                 !dampingMatrix%elementMatrix%matrix(:,:)=dampingMatrix%elementMatrix%matrix(:,:)* &
+                 !     & materialsInterpPoint%values(2+FIELD_VARIABLE%NUMBER_OF_COMPONENTS,1)
+              ENDIF
+                      
+            ENDDO !ng
+            !If Plane Stress/Strain problem multiply equation matrix by thickness
+
+            DO mhs=1,TOTAL_DEPENDENT_BASIS_EP
+               DO nhs=mhs,TOTAL_DEPENDENT_BASIS_EP
+                  !!TODO::Bring 2D plane stress/strain element thickness in through a field - element constant when it can be exported by field i/o. Currently brought in through material field (Temporary)
+                  stiffnessMatrix%elementMatrix%matrix(mhs,nhs)=stiffnessMatrix%elementMatrix%matrix(mhs,nhs)* &
+                       & materialsInterpPoint%values(1,1)
+               ENDDO !nhs
+            ENDDO !mhs
           ENDIF
 
+          IF(stiffnessMatrix%updateMatrix) THEN
+            !Transpose upper triangular portion of Stiffness matrix to give lower triangular portion. Has to be done after scale factors are applied
+            !!TODO:: Use symmetric linear equation solver or alternatively traspose to give full matrix when asemmbling or when creating solver matrices
+            !!TODO:: Better to use SIZE(equationsMatrix%elementMatrix%matrix,1) as apposed to TOTAL_DEPENDENT_BASIS_EP? Is the size re-calculated at end of every loop?
+            DO mhs=2,TOTAL_DEPENDENT_BASIS_EP
+              DO nhs=1,mhs-1
+                stiffnessMatrix%elementMatrix%matrix(mhs,nhs) = stiffnessMatrix%elementMatrix%matrix(nhs,mhs)
+              ENDDO !nhs
+            ENDDO !mhs
+         ENDIF
+         
+
+         !Construct Damping Matrix
+         IF(dampingMatrix%updateMatrix) THEN
+            !Construct Damping Matrix
+            mhs=0          
+            DO mh=1,FIELD_VARIABLE%NUMBER_OF_COMPONENTS
+               !Loop over element rows
+               DO ms=1,DEPENDENT_BASES_EP(mh)
+                  mhs=mhs+1
+                  nhs=0
+                  !Loop over element columsn
+                  DO nh=1,FIELD_VARIABLE%NUMBER_OF_COMPONENTS
+                     DO ns=1,DEPENDENT_BASES_EP(nh)
+                        nhs=nhs+1
+                        
+                        dampingMatrix%elementMatrix%matrix(mhs,nhs)=dampingMatrix%elementMatrix%matrix(mhs,nhs)+ &
+                             & materialsInterpPoint%values(3+FIELD_VARIABLE%NUMBER_OF_COMPONENTS,1)* &
+                             & stiffnessMatrix%elementMatrix%matrix(mhs,nhs)
+                     ENDDO !ns
+                  ENDDO !nh
+                  !IF(rhsVector%updateVector) rhsVector%elementVector%vector(mhs)=0.0_DP
+               ENDDO !ms
+            ENDDO !mh
+            !dampingMatrix%elementMatrix%matrix(:,:)=dampingMatrix%elementMatrix%matrix(:,:)* &
+            !     & materialsInterpPoint%values(2+FIELD_VARIABLE%NUMBER_OF_COMPONENTS,1)
+         ENDIF
+
         CASE(EQUATIONS_SET_PLATE_SUBTYPE)
           CALL FlagError("Not implemented.",err,error,*999)
         CASE(EQUATIONS_SET_SHELL_SUBTYPE)
@@ -1460,10 +1643,28 @@ SUBROUTINE LINEAR_ELASTICITY_TENSOR(EQUATIONS_SET_SUBTYPE,MATERIALS_INTERPOLATED
       E1 = MATERIALS_INTERPOLATED_POINT%values(2,1)
       v12 = MATERIALS_INTERPOLATED_POINT%values(3,1)
       v21 = v12
-      gama = 1.0_DP/(1.0_DP-v12*v21)
-      C11 = E1*gama
+      gama = 1.0_DP/(1.0_DP-v12-v21)
+      C11 = E1*gama*(1.0_DP-v12)/(1.0_DP+v12)
+      C22 = C11
+      C12 = E1*v21/((1.0_DP+v21)*(1.0_DP-2.0_DP*v21))
+      !C12 = v12*C11
+      C21 = C12
+      C66 = E1/(2.0_DP*(1.0_DP+v12)) != G12
+      ELASTICITY_TENSOR(1,1)=C11
+      ELASTICITY_TENSOR(1,2)=C21
+      ELASTICITY_TENSOR(2,1)=C21
+      ELASTICITY_TENSOR(2,2)=C22
+      ELASTICITY_TENSOR(6,6)=C66
+    CASE(EQUATIONS_SET_TWO_DIM_PLANE_STRESS_STOKES_DAMPING_SUBTYPE)
+      !Plane Stress Isotropic Elasticity Tensor
+      E1 = MATERIALS_INTERPOLATED_POINT%values(2,1)
+      v12 = MATERIALS_INTERPOLATED_POINT%values(3,1)
+      v21 = v12
+      gama = 1.0_DP/(1.0_DP-v12-v21)
+      C11 = E1*gama*(1.0_DP-v12)/(1.0_DP+v12)
       C22 = C11
-      C12 = C11*v21
+      C12 = E1*v21/((1.0_DP+v21)*(1.0_DP-2.0_DP*v21))
+      !C12 = v12*C11
       C21 = C12
       C66 = E1/(2.0_DP*(1.0_DP+v12)) != G12
       ELASTICITY_TENSOR(1,1)=C11
@@ -1493,6 +1694,11 @@ SUBROUTINE LINEAR_ELASTICITY_TENSOR(EQUATIONS_SET_SUBTYPE,MATERIALS_INTERPOLATED
       E1 = MATERIALS_INTERPOLATED_POINT%values(2,1)
       C11 = E1
       ELASTICITY_TENSOR(1,1)=C11
+   CASE(EQUATIONS_SET_ONE_DIM_STOKES_DAMPING_SUBTYPE)
+      !Plane Strain Isotropic Linear Elasticity Tensor
+      E1 = MATERIALS_INTERPOLATED_POINT%values(2,1)
+      C11 = E1
+      ELASTICITY_TENSOR(1,1)=C11
     CASE(EQUATIONS_SET_PLATE_SUBTYPE)
       CALL FlagError("Not implemented.",err,error,*999)
     CASE(EQUATIONS_SET_SHELL_SUBTYPE)
@@ -1571,6 +1777,9 @@ SUBROUTINE LINEAR_ELASTICITY_EQUATIONS_SET_SETUP(EQUATIONS_SET,EQUATIONS_SET_SET
           END SELECT
         CASE(EQUATIONS_SET_SETUP_GEOMETRY_TYPE)
           !Do nothing???
+        !-----------------------------------------------------------------
+        ! D e p e n d e n t   f i e l d
+        !-----------------------------------------------------------------
         CASE(EQUATIONS_SET_SETUP_DEPENDENT_TYPE)
           SELECT CASE(EQUATIONS_SET_SETUP%ACTION_TYPE)
           CASE(EQUATIONS_SET_SETUP_START_ACTION)
@@ -1690,6 +1899,9 @@ SUBROUTINE LINEAR_ELASTICITY_EQUATIONS_SET_SETUP(EQUATIONS_SET,EQUATIONS_SET_SET
               & " is invalid for a linear elasticity equation"
             CALL FlagError(localError,err,error,*999)
           END SELECT
+        !-----------------------------------------------------------------
+        ! M a t e r i a l   f i e l d
+        !-----------------------------------------------------------------          
         CASE(EQUATIONS_SET_SETUP_MATERIALS_TYPE)
           SELECT CASE(EQUATIONS_SET_SETUP%ACTION_TYPE)
           CASE(EQUATIONS_SET_SETUP_START_ACTION)
@@ -1767,6 +1979,9 @@ SUBROUTINE LINEAR_ELASTICITY_EQUATIONS_SET_SETUP(EQUATIONS_SET,EQUATIONS_SET_SET
               & " is invalid for a linear elasticity equation."
             CALL FlagError(localError,err,error,*999)
           END SELECT
+        !-----------------------------------------------------------------
+        ! S o u r c e   f i e l d
+        !----------------------------------------------------------------- 
         CASE(EQUATIONS_SET_SETUP_SOURCE_TYPE)
           SELECT CASE(EQUATIONS_SET_SETUP%ACTION_TYPE)
           CASE(EQUATIONS_SET_SETUP_START_ACTION)
@@ -1921,7 +2136,8 @@ SUBROUTINE LINEAR_ELASTICITY_EQUATIONS_SET_SETUP(EQUATIONS_SET,EQUATIONS_SET_SET
       !
       ! TWO DIMENSIONAL PLANE STRESS & PLANE STRAIN ELASTICITY
       !
-      CASE(EQUATIONS_SET_TWO_DIMENSIONAL_PLANE_STRESS_SUBTYPE,EQUATIONS_SET_TWO_DIMENSIONAL_PLANE_STRAIN_SUBTYPE)
+     CASE(EQUATIONS_SET_TWO_DIMENSIONAL_PLANE_STRESS_SUBTYPE,EQUATIONS_SET_TWO_DIMENSIONAL_PLANE_STRAIN_SUBTYPE, &
+        & EQUATIONS_SET_TWO_DIM_PLANE_STRESS_STOKES_DAMPING_SUBTYPE)
         SELECT CASE(EQUATIONS_SET_SETUP%SETUP_TYPE)
         CASE(EQUATIONS_SET_SETUP_INITIAL_TYPE)
           SELECT CASE(EQUATIONS_SET_SETUP%ACTION_TYPE)
@@ -1939,6 +2155,9 @@ SUBROUTINE LINEAR_ELASTICITY_EQUATIONS_SET_SETUP(EQUATIONS_SET,EQUATIONS_SET_SET
           END SELECT
         CASE(EQUATIONS_SET_SETUP_GEOMETRY_TYPE)
           !Do nothing???
+        !-----------------------------------------------------------------
+        ! D e p e n d e n t   f i e l d
+        !-----------------------------------------------------------------           
         CASE(EQUATIONS_SET_SETUP_DEPENDENT_TYPE)
           SELECT CASE(EQUATIONS_SET_SETUP%ACTION_TYPE)
           CASE(EQUATIONS_SET_SETUP_START_ACTION)
@@ -2058,6 +2277,9 @@ SUBROUTINE LINEAR_ELASTICITY_EQUATIONS_SET_SETUP(EQUATIONS_SET,EQUATIONS_SET_SET
               & " is invalid for a linear elasticity equation"
             CALL FlagError(localError,err,error,*999)
           END SELECT
+        !-----------------------------------------------------------------
+        ! M a t e r i a l   f i e l d
+        !-----------------------------------------------------------------          
         CASE(EQUATIONS_SET_SETUP_MATERIALS_TYPE)
           SELECT CASE(EQUATIONS_SET_SETUP%ACTION_TYPE)
           CASE(EQUATIONS_SET_SETUP_START_ACTION)
@@ -2107,12 +2329,17 @@ SUBROUTINE LINEAR_ELASTICITY_EQUATIONS_SET_SETUP(EQUATIONS_SET,EQUATIONS_SET_SET
                 CALL FIELD_DATA_TYPE_CHECK(EQUATIONS_SET_SETUP%FIELD,FIELD_U_VARIABLE_TYPE,FIELD_DP_TYPE,err,error,*999)
                 !2 Components for 2D Isotropic Linear Elasticity
                 !TODO:: Temporarily set to 3 to allow thickness to passed in. Remove once a thickness, element constant field is defined and can be exported/viewed by cmgui
-                CALL FIELD_NUMBER_OF_COMPONENTS_CHECK(EQUATIONS_SET_SETUP%FIELD,FIELD_U_VARIABLE_TYPE,3,err,error,*999) 
+                SELECT CASE(EQUATIONS_SET%SPECIFICATION(3))
+                CASE(EQUATIONS_SET_TWO_DIMENSIONAL_PLANE_STRESS_SUBTYPE)
+                   CALL FIELD_NUMBER_OF_COMPONENTS_CHECK(EQUATIONS_SET_SETUP%FIELD,FIELD_U_VARIABLE_TYPE,3,err,error,*999)
+                CASE(EQUATIONS_SET_TWO_DIM_PLANE_STRESS_STOKES_DAMPING_SUBTYPE)
+                   CALL FIELD_NUMBER_OF_COMPONENTS_CHECK(EQUATIONS_SET_SETUP%FIELD,FIELD_U_VARIABLE_TYPE,5,err,error,*999)
+                END SELECT
               ENDIF
             ELSE
               CALL FlagError("Equations materials is not associated.",err,error,*999)
             ENDIF
-          CASE(EQUATIONS_SET_SETUP_FINISH_ACTION)
+         CASE(EQUATIONS_SET_SETUP_FINISH_ACTION)
             EQUATIONS_MATERIALS=>EQUATIONS_SET%MATERIALS
             IF(ASSOCIATED(EQUATIONS_MATERIALS)) THEN
               IF(EQUATIONS_MATERIALS%MATERIALS_FIELD_AUTO_CREATED) THEN
@@ -2120,9 +2347,17 @@ SUBROUTINE LINEAR_ELASTICITY_EQUATIONS_SET_SETUP(EQUATIONS_SET,EQUATIONS_SET_SET
                 CALL FIELD_CREATE_FINISH(EQUATIONS_MATERIALS%MATERIALS_FIELD,err,error,*999)
                 !Set the default values for the materials field
                 CALL FIELD_COMPONENT_VALUES_INITIALISE(EQUATIONS_MATERIALS%MATERIALS_FIELD,FIELD_U_VARIABLE_TYPE, &
-                & FIELD_VALUES_SET_TYPE,1,30.0E6_DP,err,error,*999) !Young's Modulus
+                     & FIELD_VALUES_SET_TYPE,1,1.0_DP,err,error,*999) !Thickness
                 CALL FIELD_COMPONENT_VALUES_INITIALISE(EQUATIONS_MATERIALS%MATERIALS_FIELD,FIELD_U_VARIABLE_TYPE, &
-                & FIELD_VALUES_SET_TYPE,2,0.25_DP,err,error,*999) !Poisson's Ratio
+                     & FIELD_VALUES_SET_TYPE,2,30.0E6_DP,err,error,*999) !Young's Modulus
+                CALL FIELD_COMPONENT_VALUES_INITIALISE(EQUATIONS_MATERIALS%MATERIALS_FIELD,FIELD_U_VARIABLE_TYPE, &
+                     & FIELD_VALUES_SET_TYPE,3,0.25_DP,err,error,*999) !Poisson's Ratio
+                IF(EQUATIONS_SET%SPECIFICATION(3) == EQUATIONS_SET_TWO_DIM_PLANE_STRESS_STOKES_DAMPING_SUBTYPE) THEN
+                   CALL FIELD_COMPONENT_VALUES_INITIALISE(EQUATIONS_MATERIALS%MATERIALS_FIELD,FIELD_U_VARIABLE_TYPE, &
+                        & FIELD_VALUES_SET_TYPE,4,30.0E6_DP,err,error,*999) !Stokes Damping Coefficient
+                   CALL FIELD_COMPONENT_VALUES_INITIALISE(EQUATIONS_MATERIALS%MATERIALS_FIELD,FIELD_U_VARIABLE_TYPE, &
+                        & FIELD_VALUES_SET_TYPE,5,1.0_DP,err,error,*999) !Viscoelastic (Rayleigh) Coefficient
+                ENDIF
               ENDIF
             ELSE
               CALL FlagError("Equations set materials is not associated.",err,error,*999)
@@ -2133,6 +2368,9 @@ SUBROUTINE LINEAR_ELASTICITY_EQUATIONS_SET_SETUP(EQUATIONS_SET,EQUATIONS_SET_SET
               & " is invalid for a linear elasticity equation."
             CALL FlagError(localError,err,error,*999)
           END SELECT
+        !-----------------------------------------------------------------
+        ! S o u r c e   f i e l d
+        !----------------------------------------------------------------- 
         CASE(EQUATIONS_SET_SETUP_SOURCE_TYPE)
           SELECT CASE(EQUATIONS_SET_SETUP%ACTION_TYPE)
           CASE(EQUATIONS_SET_SETUP_START_ACTION)
@@ -2213,7 +2451,12 @@ SUBROUTINE LINEAR_ELASTICITY_EQUATIONS_SET_SETUP(EQUATIONS_SET,EQUATIONS_SET_SET
               !Create the equations
               CALL Equations_CreateStart(EQUATIONS_SET,equations,err,error,*999)
               CALL Equations_LinearityTypeSet(equations,EQUATIONS_LINEAR,err,error,*999)
-              CALL Equations_TimeDependenceTypeSet(equations,EQUATIONS_STATIC,err,error,*999)
+              SELECT CASE(EQUATIONS_SET%SPECIFICATION(3))
+              CASE(EQUATIONS_SET_TWO_DIMENSIONAL_PLANE_STRESS_SUBTYPE,EQUATIONS_SET_TWO_DIMENSIONAL_PLANE_STRAIN_SUBTYPE)
+                 CALL Equations_TimeDependenceTypeSet(equations,EQUATIONS_STATIC,err,error,*999)
+              CASE(EQUATIONS_SET_TWO_DIM_PLANE_STRESS_STOKES_DAMPING_SUBTYPE)
+                 CALL Equations_TimeDependenceTypeSet(equations,EQUATIONS_FIRST_ORDER_DYNAMIC,err,error,*999)
+              END SELECT
             ELSE
               CALL FlagError("Equations set dependent field has not bee finished.",err,error,*999)
             ENDIF
@@ -2225,27 +2468,69 @@ SUBROUTINE LINEAR_ELASTICITY_EQUATIONS_SET_SETUP(EQUATIONS_SET,EQUATIONS_SET_SET
               CALL Equations_CreateFinish(equations,err,error,*999)
               NULLIFY(vectorEquations)
               CALL Equations_VectorEquationsGet(equations,vectorEquations,err,error,*999)
-              !Create the equations mapping.
-              CALL EquationsMapping_VectorCreateStart(vectorEquations,FIELD_DELUDELN_VARIABLE_TYPE,vectorMapping,err,error,*999)
-              CALL EquationsMapping_LinearMatricesNumberSet(vectorMapping,1,err,error,*999)
-              CALL EquationsMapping_LinearMatricesVariableTypesSet(vectorMapping,[FIELD_U_VARIABLE_TYPE], &
-                & err,error,*999)
-              CALL EquationsMapping_RHSVariableTypeSet(vectorMapping,FIELD_DELUDELN_VARIABLE_TYPE,err,error,*999)
-              CALL EquationsMapping_VectorCreateFinish(vectorMapping,err,error,*999)
-              !Create the equations matrices
-              CALL EquationsMatrices_VectorCreateStart(vectorEquations,vectorMatrices,err,error,*999)
-              SELECT CASE(equations%sparsityType)
-              CASE(EQUATIONS_MATRICES_FULL_MATRICES)
-                CALL EquationsMatrices_LinearStorageTypeSet(vectorMatrices,[MATRIX_BLOCK_STORAGE_TYPE],err,error,*999)
-              CASE(EQUATIONS_MATRICES_SPARSE_MATRICES)
-                CALL EquationsMatrices_LinearStorageTypeSet(vectorMatrices,[MATRIX_COMPRESSED_ROW_STORAGE_TYPE],err,error,*999)
-                CALL EquationsMatrices_LinearStructureTypeSet(vectorMatrices,[EQUATIONS_MATRIX_FEM_STRUCTURE],err,error,*999)
-              CASE DEFAULT
-                localError="The equations matrices sparsity type of "// &
-                  & TRIM(NumberToVString(equations%sparsityType,"*",err,error))//" is invalid."
-                CALL FlagError(localError,err,error,*999)
+              SELECT CASE(EQUATIONS_SET%SPECIFICATION(3))
+              CASE(EQUATIONS_SET_TWO_DIMENSIONAL_PLANE_STRESS_SUBTYPE,EQUATIONS_SET_TWO_DIMENSIONAL_PLANE_STRAIN_SUBTYPE)
+                 !Create the equations mapping.
+                 CALL EquationsMapping_VectorCreateStart(vectorEquations,FIELD_DELUDELN_VARIABLE_TYPE,vectorMapping,err,error,*999)
+                 CALL EquationsMapping_LinearMatricesNumberSet(vectorMapping,1,err,error,*999)
+                 CALL EquationsMapping_LinearMatricesVariableTypesSet(vectorMapping,[FIELD_U_VARIABLE_TYPE], &
+                      & err,error,*999)
+                 CALL EquationsMapping_RHSVariableTypeSet(vectorMapping,FIELD_DELUDELN_VARIABLE_TYPE,err,error,*999)
+                 CALL EquationsMapping_VectorCreateFinish(vectorMapping,err,error,*999)
+                 !Create the equations matrices
+                 CALL EquationsMatrices_VectorCreateStart(vectorEquations,vectorMatrices,err,error,*999)
+                 SELECT CASE(equations%sparsityType)
+                 CASE(EQUATIONS_MATRICES_FULL_MATRICES)
+                    CALL EquationsMatrices_LinearStorageTypeSet(vectorMatrices,[MATRIX_BLOCK_STORAGE_TYPE],err,error,*999)
+                 CASE(EQUATIONS_MATRICES_SPARSE_MATRICES)
+                    CALL EquationsMatrices_LinearStorageTypeSet(vectorMatrices,[MATRIX_COMPRESSED_ROW_STORAGE_TYPE],err,error,*999)
+                    CALL EquationsMatrices_LinearStructureTypeSet(vectorMatrices,[EQUATIONS_MATRIX_FEM_STRUCTURE],err,error,*999)
+                 CASE DEFAULT
+                    localError="The equations matrices sparsity type of "// &
+                         & TRIM(NumberToVString(equations%sparsityType,"*",err,error))//" is invalid."
+                    CALL FlagError(localError,err,error,*999)
+                 END SELECT
+                 CALL EquationsMatrices_VectorCreateFinish(vectorMatrices,err,error,*999)
+              CASE(EQUATIONS_SET_TWO_DIM_PLANE_STRESS_STOKES_DAMPING_SUBTYPE)
+                 !Create the equations mapping.
+                 CALL EquationsMapping_VectorCreateStart(vectorEquations,FIELD_DELUDELN_VARIABLE_TYPE,vectorMapping,err,error,*999)
+                 CALL EquationsMapping_DynamicMatricesSet(vectorMapping,.TRUE.,.TRUE.,err,error,*999)
+                 CALL EquationsMapping_DynamicVariableTypeSet(vectorMapping,FIELD_U_VARIABLE_TYPE,err,error,*999)
+                 CALL EquationsMapping_RHSVariableTypeSet(vectorMapping,FIELD_DELUDELN_VARIABLE_TYPE,err,error,*999)
+                 
+                 !CALL EquationsMapping_SourceVariableTypeSet(vectorMapping,FIELD_U_VARIABLE_TYPE,err,error,*999)
+                 CALL EquationsMapping_VectorCreateFinish(vectorMapping,err,error,*999)
+                 !Create the equations matrices
+                 CALL EquationsMatrices_VectorCreateStart(vectorEquations,vectorMatrices,err,error,*999)
+                 !Set up matrix storage and structure
+                 IF(EQUATIONS%lumpingType==EQUATIONS_LUMPED_MATRICES) THEN
+                    !Set up lumping
+                    CALL EquationsMatrices_DynamicLumpingTypeSet(vectorMatrices, &
+                         & [EQUATIONS_MATRIX_UNLUMPED,EQUATIONS_MATRIX_LUMPED],err,error,*999)
+                    CALL EquationsMatrices_DynamicStorageTypeSet(vectorMatrices, &
+                         & [DISTRIBUTED_MATRIX_COMPRESSED_ROW_STORAGE_TYPE,DISTRIBUTED_MATRIX_DIAGONAL_STORAGE_TYPE], &
+                         & err,error,*999)
+                    CALL EquationsMatrices_DynamicStructureTypeSet(vectorMatrices, &
+                         [EQUATIONS_MATRIX_FEM_STRUCTURE,EQUATIONS_MATRIX_DIAGONAL_STRUCTURE],err,error,*999)
+                 ELSE
+                    SELECT CASE(EQUATIONS%sparsityType)
+                    CASE(EQUATIONS_MATRICES_FULL_MATRICES) 
+                       CALL EquationsMatrices_DynamicStorageTypeSet(vectorMatrices, &
+                            & [DISTRIBUTED_MATRIX_BLOCK_STORAGE_TYPE,DISTRIBUTED_MATRIX_BLOCK_STORAGE_TYPE],err,error,*999)
+                    CASE(EQUATIONS_MATRICES_SPARSE_MATRICES)
+                       CALL EquationsMatrices_DynamicStorageTypeSet(vectorMatrices, &
+                            & [DISTRIBUTED_MATRIX_COMPRESSED_ROW_STORAGE_TYPE,DISTRIBUTED_MATRIX_COMPRESSED_ROW_STORAGE_TYPE], &
+                            & err,error,*999)
+                       CALL EquationsMatrices_DynamicStructureTypeSet(vectorMatrices, &
+                            [EQUATIONS_MATRIX_FEM_STRUCTURE,EQUATIONS_MATRIX_FEM_STRUCTURE],err,error,*999)                  
+                    CASE DEFAULT
+                       localError="The equations matrices sparsity type of "// &
+                            & TRIM(NumberToVString(EQUATIONS%sparsityType,"*",err,error))//" is invalid."
+                       CALL FlagError(localError,err,error,*999)
+                    END SELECT
+                 ENDIF
+                 CALL EquationsMatrices_VectorCreateFinish(vectorMatrices,err,error,*999)
               END SELECT
-              CALL EquationsMatrices_VectorCreateFinish(vectorMatrices,err,error,*999)
             CASE(EQUATIONS_SET_BEM_SOLUTION_METHOD)
               CALL FlagError("Not implemented.",err,error,*999)
             CASE(EQUATIONS_SET_FD_SOLUTION_METHOD)
@@ -2276,7 +2561,7 @@ SUBROUTINE LINEAR_ELASTICITY_EQUATIONS_SET_SETUP(EQUATIONS_SET,EQUATIONS_SET_SET
       !
       ! ONE DIMENSIONAL ELASTICITY
       !
-      CASE(EQUATIONS_SET_ONE_DIMENSIONAL_SUBTYPE)
+      CASE(EQUATIONS_SET_ONE_DIMENSIONAL_SUBTYPE,EQUATIONS_SET_ONE_DIM_STOKES_DAMPING_SUBTYPE)
         SELECT CASE(EQUATIONS_SET_SETUP%SETUP_TYPE)
         CASE(EQUATIONS_SET_SETUP_INITIAL_TYPE)
           SELECT CASE(EQUATIONS_SET_SETUP%ACTION_TYPE)
@@ -2294,6 +2579,9 @@ SUBROUTINE LINEAR_ELASTICITY_EQUATIONS_SET_SETUP(EQUATIONS_SET,EQUATIONS_SET_SET
           END SELECT
         CASE(EQUATIONS_SET_SETUP_GEOMETRY_TYPE)
           !Do nothing???
+        !-----------------------------------------------------------------
+        ! D e p e n d e n t   f i e l d
+        !-----------------------------------------------------------------           
         CASE(EQUATIONS_SET_SETUP_DEPENDENT_TYPE)
           SELECT CASE(EQUATIONS_SET_SETUP%ACTION_TYPE)
           CASE(EQUATIONS_SET_SETUP_START_ACTION)
@@ -2413,6 +2701,9 @@ SUBROUTINE LINEAR_ELASTICITY_EQUATIONS_SET_SETUP(EQUATIONS_SET,EQUATIONS_SET_SET
               & " is invalid for a linear elasticity equation"
             CALL FlagError(localError,err,error,*999)
           END SELECT
+        !-----------------------------------------------------------------
+        ! M a t e r i a l   f i e l d
+        !-----------------------------------------------------------------          
         CASE(EQUATIONS_SET_SETUP_MATERIALS_TYPE)
           SELECT CASE(EQUATIONS_SET_SETUP%ACTION_TYPE)
           CASE(EQUATIONS_SET_SETUP_START_ACTION)
@@ -2461,7 +2752,12 @@ SUBROUTINE LINEAR_ELASTICITY_EQUATIONS_SET_SETUP(EQUATIONS_SET,EQUATIONS_SET_SET
                   & err,error,*999)
                 CALL FIELD_DATA_TYPE_CHECK(EQUATIONS_SET_SETUP%FIELD,FIELD_U_VARIABLE_TYPE,FIELD_DP_TYPE,err,error,*999)
                 !2 Components for 2D Isotropic Linear Elasticity
-                CALL FIELD_NUMBER_OF_COMPONENTS_CHECK(EQUATIONS_SET_SETUP%FIELD,FIELD_U_VARIABLE_TYPE,2,err,error,*999) 
+                SELECT CASE(EQUATIONS_SET%SPECIFICATION(3))
+                CASE(EQUATIONS_SET_ONE_DIMENSIONAL_SUBTYPE)
+                  CALL FIELD_NUMBER_OF_COMPONENTS_CHECK(EQUATIONS_SET_SETUP%FIELD,FIELD_U_VARIABLE_TYPE,2,err,error,*999)
+                CASE(EQUATIONS_SET_ONE_DIM_STOKES_DAMPING_SUBTYPE)
+                  CALL FIELD_NUMBER_OF_COMPONENTS_CHECK(EQUATIONS_SET_SETUP%FIELD,FIELD_U_VARIABLE_TYPE,4,err,error,*999)
+                END SELECT
               ENDIF
             ELSE
               CALL FlagError("Equations materials is not associated.",err,error,*999)
@@ -2477,6 +2773,12 @@ SUBROUTINE LINEAR_ELASTICITY_EQUATIONS_SET_SETUP(EQUATIONS_SET,EQUATIONS_SET_SET
                 & FIELD_VALUES_SET_TYPE,1,30.0E6_DP,err,error,*999) !Young's Modulus
                 CALL FIELD_COMPONENT_VALUES_INITIALISE(EQUATIONS_MATERIALS%MATERIALS_FIELD,FIELD_U_VARIABLE_TYPE, &
                 & FIELD_VALUES_SET_TYPE,2,0.25_DP,err,error,*999) !Poisson's Ratio
+                IF(EQUATIONS_SET%SPECIFICATION(3) == EQUATIONS_SET_ONE_DIM_STOKES_DAMPING_SUBTYPE) THEN
+                   CALL FIELD_COMPONENT_VALUES_INITIALISE(EQUATIONS_MATERIALS%MATERIALS_FIELD,FIELD_U_VARIABLE_TYPE, &
+                        & FIELD_VALUES_SET_TYPE,3,30.0E6_DP,err,error,*999) !Stokes Damping Coefficient
+                   CALL FIELD_COMPONENT_VALUES_INITIALISE(EQUATIONS_MATERIALS%MATERIALS_FIELD,FIELD_U_VARIABLE_TYPE, &
+                        & FIELD_VALUES_SET_TYPE,4,1.0_DP,err,error,*999) !Viscoelastic (Rayleigh) Coefficient
+                ENDIF                
               ENDIF
             ELSE
               CALL FlagError("Equations set materials is not associated.",err,error,*999)
@@ -2487,6 +2789,9 @@ SUBROUTINE LINEAR_ELASTICITY_EQUATIONS_SET_SETUP(EQUATIONS_SET,EQUATIONS_SET_SET
               & " is invalid for a linear elasticity equation."
             CALL FlagError(localError,err,error,*999)
           END SELECT
+        !-----------------------------------------------------------------
+        ! S o u r c e   f i e l d
+        !----------------------------------------------------------------- 
         CASE(EQUATIONS_SET_SETUP_SOURCE_TYPE)
           SELECT CASE(EQUATIONS_SET_SETUP%ACTION_TYPE)
           CASE(EQUATIONS_SET_SETUP_START_ACTION)
@@ -2565,11 +2870,16 @@ SUBROUTINE LINEAR_ELASTICITY_EQUATIONS_SET_SETUP(EQUATIONS_SET,EQUATIONS_SET_SET
               !Create the equations
               CALL Equations_CreateStart(EQUATIONS_SET,equations,err,error,*999)
               CALL Equations_LinearityTypeSet(equations,EQUATIONS_LINEAR,err,error,*999)
-              CALL Equations_TimeDependenceTypeSet(equations,EQUATIONS_STATIC,err,error,*999)
+              SELECT CASE(EQUATIONS_SET%SPECIFICATION(3))
+              CASE(EQUATIONS_SET_ONE_DIMENSIONAL_SUBTYPE)
+                 CALL Equations_TimeDependenceTypeSet(equations,EQUATIONS_STATIC,err,error,*999)                 
+              CASE(EQUATIONS_SET_ONE_DIM_STOKES_DAMPING_SUBTYPE)
+                 CALL Equations_TimeDependenceTypeSet(equations,EQUATIONS_FIRST_ORDER_DYNAMIC,err,error,*999) 
+              END SELECT
             ELSE
-              CALL FlagError("Equations set dependent field has not bee finished.",err,error,*999)
+              CALL FlagError("Equations set dependent field has not been finished.",err,error,*999)
             ENDIF
-          CASE(EQUATIONS_SET_SETUP_FINISH_ACTION)
+          CASE(EQUATIONS_SET_SETUP_FINISH_ACTION)         
             SELECT CASE(EQUATIONS_SET%SOLUTION_METHOD)
             CASE(EQUATIONS_SET_FEM_SOLUTION_METHOD)
               !Finish the equations creation
@@ -2577,27 +2887,70 @@ SUBROUTINE LINEAR_ELASTICITY_EQUATIONS_SET_SETUP(EQUATIONS_SET,EQUATIONS_SET_SET
               CALL Equations_CreateFinish(equations,err,error,*999)
               NULLIFY(vectorEquations)
               CALL Equations_VectorEquationsGet(equations,vectorEquations,err,error,*999)
-              !Create the equations mapping.
-              CALL EquationsMapping_VectorCreateStart(vectorEquations,FIELD_DELUDELN_VARIABLE_TYPE,vectorMapping,err,error,*999)
-              CALL EquationsMapping_LinearMatricesNumberSet(vectorMapping,1,err,error,*999)
-              CALL EquationsMapping_LinearMatricesVariableTypesSet(vectorMapping,[FIELD_U_VARIABLE_TYPE], &
-                & err,error,*999)
-              CALL EquationsMapping_RHSVariableTypeSet(vectorMapping,FIELD_DELUDELN_VARIABLE_TYPE,err,error,*999)
-              CALL EquationsMapping_VectorCreateFinish(vectorMapping,err,error,*999)
-              !Create the equations matrices
-              CALL EquationsMatrices_VectorCreateStart(vectorEquations,vectorMatrices,err,error,*999)
-              SELECT CASE(equations%sparsityType)
-              CASE(EQUATIONS_MATRICES_FULL_MATRICES)
-                CALL EquationsMatrices_LinearStorageTypeSet(vectorMatrices,[MATRIX_BLOCK_STORAGE_TYPE],err,error,*999)
-              CASE(EQUATIONS_MATRICES_SPARSE_MATRICES)
-                CALL EquationsMatrices_LinearStorageTypeSet(vectorMatrices,[MATRIX_COMPRESSED_ROW_STORAGE_TYPE],err,error,*999)
-                CALL EquationsMatrices_LinearStructureTypeSet(vectorMatrices,[EQUATIONS_MATRIX_FEM_STRUCTURE],err,error,*999)
-              CASE DEFAULT
-                localError="The equations matrices sparsity type of "// &
-                  & TRIM(NumberToVString(equations%sparsityType,"*",err,error))//" is invalid."
-                CALL FlagError(localError,err,error,*999)
+              SELECT CASE(EQUATIONS_SET%SPECIFICATION(3))                 
+              CASE(EQUATIONS_SET_ONE_DIMENSIONAL_SUBTYPE)
+                 !Create the equations mapping.
+                 CALL EquationsMapping_VectorCreateStart(vectorEquations,FIELD_DELUDELN_VARIABLE_TYPE,vectorMapping,err,error,*999)
+                 CALL EquationsMapping_LinearMatricesNumberSet(vectorMapping,1,err,error,*999)
+                 CALL EquationsMapping_LinearMatricesVariableTypesSet(vectorMapping,[FIELD_U_VARIABLE_TYPE], &
+                      & err,error,*999)
+                 CALL EquationsMapping_RHSVariableTypeSet(vectorMapping,FIELD_DELUDELN_VARIABLE_TYPE,err,error,*999)
+                 CALL EquationsMapping_VectorCreateFinish(vectorMapping,err,error,*999)
+                 !Create the equations matrices
+                 CALL EquationsMatrices_VectorCreateStart(vectorEquations,vectorMatrices,err,error,*999)
+                 SELECT CASE(equations%sparsityType)
+                 CASE(EQUATIONS_MATRICES_FULL_MATRICES)
+                    CALL EquationsMatrices_LinearStorageTypeSet(vectorMatrices,[MATRIX_BLOCK_STORAGE_TYPE],err,error,*999)
+                 CASE(EQUATIONS_MATRICES_SPARSE_MATRICES)
+                    CALL EquationsMatrices_LinearStorageTypeSet(vectorMatrices,[MATRIX_COMPRESSED_ROW_STORAGE_TYPE],err,error,*999)
+                    CALL EquationsMatrices_LinearStructureTypeSet(vectorMatrices,[EQUATIONS_MATRIX_FEM_STRUCTURE],err,error,*999)
+                 CASE DEFAULT
+                    localError="The equations matrices sparsity type of "// &
+                         & TRIM(NumberToVString(equations%sparsityType,"*",err,error))//" is invalid."
+                    CALL FlagError(localError,err,error,*999)
+                 END SELECT
+                 CALL EquationsMatrices_VectorCreateFinish(vectorMatrices,err,error,*999)
+              CASE(EQUATIONS_SET_ONE_DIM_STOKES_DAMPING_SUBTYPE)
+                 !Finish the creation of the equations
+                 !Create the equations mapping.
+                 CALL EquationsMapping_VectorCreateStart(vectorEquations,FIELD_DELUDELN_VARIABLE_TYPE,vectorMapping,err,error,*999)
+                 CALL EquationsMapping_DynamicMatricesSet(vectorMapping,.TRUE.,.TRUE.,err,error,*999)
+                 CALL EquationsMapping_DynamicVariableTypeSet(vectorMapping,FIELD_U_VARIABLE_TYPE,err,error,*999)
+                 CALL EquationsMapping_RHSVariableTypeSet(vectorMapping,FIELD_DELUDELN_VARIABLE_TYPE,err,error,*999)
+                 
+                 !CALL EquationsMapping_SourceVariableTypeSet(vectorMapping,FIELD_U_VARIABLE_TYPE,err,error,*999)
+                 CALL EquationsMapping_VectorCreateFinish(vectorMapping,err,error,*999)
+                 !Create the equations matrices
+                 CALL EquationsMatrices_VectorCreateStart(vectorEquations,vectorMatrices,err,error,*999)
+                 !Set up matrix storage and structure
+                 IF(EQUATIONS%lumpingType==EQUATIONS_LUMPED_MATRICES) THEN
+                    !Set up lumping
+                    CALL EquationsMatrices_DynamicLumpingTypeSet(vectorMatrices, &
+                         & [EQUATIONS_MATRIX_UNLUMPED,EQUATIONS_MATRIX_LUMPED],err,error,*999)
+                    CALL EquationsMatrices_DynamicStorageTypeSet(vectorMatrices, &
+                         & [DISTRIBUTED_MATRIX_COMPRESSED_ROW_STORAGE_TYPE,DISTRIBUTED_MATRIX_DIAGONAL_STORAGE_TYPE], &
+                         & err,error,*999)
+                    CALL EquationsMatrices_DynamicStructureTypeSet(vectorMatrices, &
+                         [EQUATIONS_MATRIX_FEM_STRUCTURE,EQUATIONS_MATRIX_DIAGONAL_STRUCTURE],err,error,*999)
+                 ELSE
+                    SELECT CASE(EQUATIONS%sparsityType)
+                    CASE(EQUATIONS_MATRICES_FULL_MATRICES) 
+                       CALL EquationsMatrices_DynamicStorageTypeSet(vectorMatrices, &
+                            & [DISTRIBUTED_MATRIX_BLOCK_STORAGE_TYPE,DISTRIBUTED_MATRIX_BLOCK_STORAGE_TYPE],err,error,*999)
+                    CASE(EQUATIONS_MATRICES_SPARSE_MATRICES)
+                       CALL EquationsMatrices_DynamicStorageTypeSet(vectorMatrices, &
+                            & [DISTRIBUTED_MATRIX_COMPRESSED_ROW_STORAGE_TYPE,DISTRIBUTED_MATRIX_COMPRESSED_ROW_STORAGE_TYPE], &
+                            & err,error,*999)
+                       CALL EquationsMatrices_DynamicStructureTypeSet(vectorMatrices, &
+                            [EQUATIONS_MATRIX_FEM_STRUCTURE,EQUATIONS_MATRIX_FEM_STRUCTURE],err,error,*999)                  
+                    CASE DEFAULT
+                       localError="The equations matrices sparsity type of "// &
+                            & TRIM(NumberToVString(EQUATIONS%sparsityType,"*",err,error))//" is invalid."
+                       CALL FlagError(localError,err,error,*999)
+                    END SELECT
+                 ENDIF
+                 CALL EquationsMatrices_VectorCreateFinish(vectorMatrices,err,error,*999)                 
               END SELECT
-              CALL EquationsMatrices_VectorCreateFinish(vectorMatrices,err,error,*999)
             CASE(EQUATIONS_SET_BEM_SOLUTION_METHOD)
               CALL FlagError("Not implemented.",err,error,*999)
             CASE(EQUATIONS_SET_FD_SOLUTION_METHOD)
@@ -2704,6 +3057,24 @@ SUBROUTINE LinearElasticity_EquationsSetSolutionMethodSet(EQUATIONS_SET,SOLUTION
           localError="The specified solution method of "//TRIM(NumberToVString(SOLUTION_METHOD,"*",err,error))//" is invalid."
           CALL FlagError(localError,err,error,*999)
         END SELECT
+      CASE(EQUATIONS_SET_TWO_DIM_PLANE_STRESS_STOKES_DAMPING_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
+          localError="The specified solution method of "//TRIM(NumberToVString(SOLUTION_METHOD,"*",err,error))//" is invalid."
+          CALL FlagError(localError,err,error,*999)
+        END SELECT
       CASE(EQUATIONS_SET_TWO_DIMENSIONAL_PLANE_STRAIN_SUBTYPE)
         SELECT CASE(SOLUTION_METHOD)
         CASE(EQUATIONS_SET_FEM_SOLUTION_METHOD)
@@ -2740,6 +3111,24 @@ SUBROUTINE LinearElasticity_EquationsSetSolutionMethodSet(EQUATIONS_SET,SOLUTION
           localError="The specified solution method of "//TRIM(NumberToVString(SOLUTION_METHOD,"*",err,error))//" is invalid."
           CALL FlagError(localError,err,error,*999)
         END SELECT
+      CASE(EQUATIONS_SET_ONE_DIM_STOKES_DAMPING_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
+          localError="The specified solution method of "//TRIM(NumberToVString(SOLUTION_METHOD,"*",err,error))//" is invalid."
+          CALL FlagError(localError,err,error,*999)
+        END SELECT
       CASE(EQUATIONS_SET_PLATE_SUBTYPE)
         CALL FlagError("Not implemented.",err,error,*999)
       CASE(EQUATIONS_SET_SHELL_SUBTYPE)
@@ -2787,6 +3176,8 @@ SUBROUTINE LinearElasticity_EquationsSetSpecificationSet(equationsSet,specificat
       CASE(EQUATIONS_SET_THREE_DIMENSIONAL_SUBTYPE, &
           & EQUATIONS_SET_TWO_DIMENSIONAL_PLANE_STRESS_SUBTYPE, &
           & EQUATIONS_SET_TWO_DIMENSIONAL_PLANE_STRAIN_SUBTYPE, &
+          & EQUATIONS_SET_TWO_DIM_PLANE_STRESS_STOKES_DAMPING_SUBTYPE, &
+          & EQUATIONS_SET_ONE_DIM_STOKES_DAMPING_SUBTYPE, &          
           & EQUATIONS_SET_ONE_DIMENSIONAL_SUBTYPE)
         !ok
       CASE(EQUATIONS_SET_PLATE_SUBTYPE)
@@ -2849,77 +3240,110 @@ SUBROUTINE LINEAR_ELASTICITY_PROBLEM_SETUP(PROBLEM,PROBLEM_SETUP,err,error,*)
       ELSE IF(SIZE(PROBLEM%SPECIFICATION,1)<3) THEN
         CALL FlagError("Problem specification must have three entries for a linear elasticity problem.",err,error,*999)
       END IF
-      SELECT CASE(PROBLEM%SPECIFICATION(3))
-      CASE(PROBLEM_NO_SUBTYPE)
-        SELECT CASE(PROBLEM_SETUP%SETUP_TYPE)
-        CASE(PROBLEM_SETUP_INITIAL_TYPE)
-          SELECT CASE(PROBLEM_SETUP%ACTION_TYPE)
-          CASE(PROBLEM_SETUP_START_ACTION)
+      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)
+         CASE(PROBLEM_SETUP_FINISH_ACTION)
             !Do nothing????
-          CASE DEFAULT
+         CASE DEFAULT
             localError="The action type of "//TRIM(NumberToVString(PROBLEM_SETUP%ACTION_TYPE,"*",err,error))// &
               & " for a setup type of "//TRIM(NumberToVString(PROBLEM_SETUP%SETUP_TYPE,"*",err,error))// &
               & " is invalid for a linear elasticity problem."
             CALL FlagError(localError,err,error,*999)
-          END SELECT
-        CASE(PROBLEM_SETUP_CONTROL_TYPE)
-          SELECT CASE(PROBLEM_SETUP%ACTION_TYPE)
-          CASE(PROBLEM_SETUP_START_ACTION)
+         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)
-          CASE(PROBLEM_SETUP_FINISH_ACTION)
+            IF(PROBLEM%SPECIFICATION(3) == PROBLEM_STOKES_DAMPING_LINEAR_ELASTICITY_SUBTYPE) THEN
+               CALL CONTROL_LOOP_TYPE_SET(CONTROL_LOOP,PROBLEM_CONTROL_TIME_LOOP_TYPE,err,error,*999)
+            ENDIF
+         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
+         CASE DEFAULT
             localError="The action type of "//TRIM(NumberToVString(PROBLEM_SETUP%ACTION_TYPE,"*",err,error))// &
               & " for a setup type of "//TRIM(NumberToVString(PROBLEM_SETUP%SETUP_TYPE,"*",err,error))// &
               & " is invalid for a linear elasticity problem."
             CALL FlagError(localError,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)
+         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 linear solver
-            CALL SOLVERS_SOLVER_GET(SOLVERS,1,SOLVER,err,error,*999)
-            CALL SOLVER_TYPE_SET(SOLVER,SOLVER_LINEAR_TYPE,err,error,*999)
-            !Set solver defaults
-            CALL SOLVER_LIBRARY_TYPE_SET(SOLVER,SOLVER_PETSC_LIBRARY,err,error,*999)
-          CASE(PROBLEM_SETUP_FINISH_ACTION)
+            SELECT CASE(PROBLEM%SPECIFICATION(3))
+            !Standard Linear Elasticity
+            CASE(PROBLEM_NO_SUBTYPE) 
+               CALL SOLVERS_NUMBER_SET(SOLVERS,1,err,error,*999)
+               !Set the solver to be a linear solver
+               CALL SOLVERS_SOLVER_GET(SOLVERS,1,SOLVER,err,error,*999)
+               CALL SOLVER_TYPE_SET(SOLVER,SOLVER_LINEAR_TYPE,err,error,*999)
+               !Set solver defaults
+               CALL SOLVER_LIBRARY_TYPE_SET(SOLVER,SOLVER_PETSC_LIBRARY,err,error,*999)
+            !Linear Elasticity with Stokes Damping
+            CASE(PROBLEM_STOKES_DAMPING_LINEAR_ELASTICITY_SUBTYPE)
+               CALL SOLVERS_NUMBER_SET(SOLVERS,1,err,error,*999)
+               !Set the solver to be a linear 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)
+               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 DEFAULT
+               localError="Problem subtype "//TRIM(NumberToVString(PROBLEM%SPECIFICATION(3),"*",err,error))// &
+                    & " is not valid for a linear elasticity type of an elasticity problem class."
+               CALL FlagError(localError,err,error,*999)
+            END SELECT
+            
+         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
+         CASE DEFAULT
             localError="The action type of "//TRIM(NumberToVString(PROBLEM_SETUP%ACTION_TYPE,"*",err,error))// &
               & " for a setup type of "//TRIM(NumberToVString(PROBLEM_SETUP%SETUP_TYPE,"*",err,error))// &
               & " is invalid for a linear elasticity problem."
             CALL FlagError(localError,err,error,*999)
-          END SELECT
-        CASE(PROBLEM_SETUP_SOLVER_EQUATIONS_TYPE)
-          SELECT CASE(PROBLEM_SETUP%ACTION_TYPE)
-          CASE(PROBLEM_SETUP_START_ACTION)
+         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_STATIC,err,error,*999)
-            CALL SOLVER_EQUATIONS_SPARSITY_TYPE_SET(SOLVER_EQUATIONS,SOLVER_SPARSE_MATRICES,err,error,*999)
-          CASE(PROBLEM_SETUP_FINISH_ACTION)
+            SELECT CASE(PROBLEM%SPECIFICATION(3))
+            CASE(PROBLEM_NO_SUBTYPE)
+               !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_STATIC,err,error,*999)
+               CALL SOLVER_EQUATIONS_SPARSITY_TYPE_SET(SOLVER_EQUATIONS,SOLVER_SPARSE_MATRICES,err,error,*999)
+            CASE(PROBLEM_STOKES_DAMPING_LINEAR_ELASTICITY_SUBTYPE)
+               !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 DEFAULT
+               localError="Problem subtype "//TRIM(NumberToVString(PROBLEM%SPECIFICATION(3),"*",err,error))// &
+                    & " is not valid for a linear elasticity type of an elasticity problem class."
+               CALL FlagError(localError,err,error,*999)
+            END SELECT
+         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)
@@ -2929,24 +3353,19 @@ SUBROUTINE LINEAR_ELASTICITY_PROBLEM_SETUP(PROBLEM,PROBLEM_SETUP,err,error,*)
             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
+         CASE DEFAULT
             localError="The action type of "//TRIM(NumberToVString(PROBLEM_SETUP%ACTION_TYPE,"*",err,error))// &
               & " for a setup type of "//TRIM(NumberToVString(PROBLEM_SETUP%SETUP_TYPE,"*",err,error))// &
               & " is invalid for a linear elasticity problem."
             CALL FlagError(localError,err,error,*999)
-          END SELECT
-        CASE DEFAULT
-          localError="The setup type of "//TRIM(NumberToVString(PROBLEM_SETUP%SETUP_TYPE,"*",err,error))// &
-            & " is invalid for a linear elasticity problem."
-          CALL FlagError(localError,err,error,*999)
-        END SELECT
+         END SELECT
       CASE DEFAULT
-        localError="Problem subtype "//TRIM(NumberToVString(PROBLEM%SPECIFICATION(3),"*",err,error))// &
-          & " is not valid for a linear elasticity type of an elasticity problem class."
-        CALL FlagError(localError,err,error,*999)
+         localError="The setup type of "//TRIM(NumberToVString(PROBLEM_SETUP%SETUP_TYPE,"*",err,error))// &
+            & " is invalid for a linear elasticity problem."
+         CALL FlagError(localError,err,error,*999)
       END SELECT
     ELSE
-      CALL FlagError("Problem is not associated.",err,error,*999)
+       CALL FlagError("Problem is not associated.",err,error,*999)
     ENDIF
 
     EXITS("LINEAR_ELASTICITY_PROBLEM_SETUP")
@@ -2977,7 +3396,8 @@ SUBROUTINE LinearElasticity_ProblemSpecificationSet(problem,problemSpecification
       IF(SIZE(problemSpecification,1)==3) THEN
         problemSubtype=problemSpecification(3)
         SELECT CASE(problemSubtype)
-        CASE(PROBLEM_NO_SUBTYPE)
+        CASE(PROBLEM_NO_SUBTYPE, &
+             & PROBLEM_STOKES_DAMPING_LINEAR_ELASTICITY_SUBTYPE)
           !ok
         CASE DEFAULT
           localError="The third problem specification of "//TRIM(NumberToVstring(problemSubtype,"*",err,error))// &
diff --git a/src/opencmiss_iron.F90 b/src/opencmiss_iron.F90
index 300cf38f..5dace590 100644
--- a/src/opencmiss_iron.F90
+++ b/src/opencmiss_iron.F90
@@ -2515,6 +2515,10 @@ MODULE OpenCMISS_Iron
   INTEGER(INTG), PARAMETER :: CMFE_EQUATIONS_SET_ONE_DIMENSIONAL_SUBTYPE = EQUATIONS_SET_ONE_DIMENSIONAL_SUBTYPE !<One dimensional linear elasticity equations set subtype \see OpenCMISS_EquationsSetSubtypes,OpenCMISS
   INTEGER(INTG), PARAMETER :: CMFE_EQUATIONS_SET_PLATE_SUBTYPE = EQUATIONS_SET_PLATE_SUBTYPE !<Plate linear elasticity equations set subtype \see OpenCMISS_EquationsSetSubtypes,OpenCMISS
   INTEGER(INTG), PARAMETER :: CMFE_EQUATIONS_SET_SHELL_SUBTYPE = EQUATIONS_SET_SHELL_SUBTYPE !<Shell linear elasticity equations set subtype \see OpenCMISS_EquationsSetSubtypes,OpenCMISS
+  INTEGER(INTG), PARAMETER :: CMFE_EQUATIONS_SET_ONE_DIM_STOKES_DAMPING_SUBTYPE = &
+    & EQUATIONS_SET_ONE_DIM_STOKES_DAMPING_SUBTYPE !<Plane stress linear elasticity equations set subtype \see OpenCMISS_EquationsSetSubtypes,OpenCMISS
+  INTEGER(INTG), PARAMETER :: CMFE_EQUATIONS_SET_TWO_DIM_PLANE_STRESS_STOKES_DAMPING_SUBTYPE = &
+    & EQUATIONS_SET_TWO_DIM_PLANE_STRESS_STOKES_DAMPING_SUBTYPE !<Plane stress linear elasticity equations set subtype \see OpenCMISS_EquationsSetSubtypes,OpenCMISS
   INTEGER(INTG), PARAMETER :: CMFE_EQUATIONS_SET_MOONEY_RIVLIN_SUBTYPE = EQUATIONS_SET_MOONEY_RIVLIN_SUBTYPE !< Mooney-Rivlin constitutive law for finite elasticity equations set subtype \see OpenCMISS_EquationsSetSubtypes,OpenCMISS
   INTEGER(INTG), PARAMETER :: CMFE_EQUATIONS_SET_INCOMPRESSIBLE_MOONEY_RIVLIN_SUBTYPE = &
   & EQUATIONS_SET_INCOMPRESSIBLE_MOONEY_RIVLIN_SUBTYPE !< Incompressible Mooney-Rivlin constitutive law for finite elasticity equations set subtype \see OpenCMISS_EquationsSetSubtypes,OpenCMISS
@@ -3081,6 +3085,8 @@ MODULE OpenCMISS_Iron
     & CMFE_EQUATIONS_SET_TWO_DIMENSIONAL_PLANE_STRAIN_SUBTYPE,CMFE_EQUATIONS_SET_ONE_DIMENSIONAL_SUBTYPE, &
     & CMFE_EQUATIONS_SET_PLATE_SUBTYPE, &
     & CMFE_EQUATIONS_SET_SHELL_SUBTYPE, &
+    & CMFE_EQUATIONS_SET_ONE_DIM_STOKES_DAMPING_SUBTYPE, &
+    & CMFE_EQUATIONS_SET_TWO_DIM_PLANE_STRESS_STOKES_DAMPING_SUBTYPE, &
     & CMFE_EQUATIONS_SET_INCOMPRESSIBLE_MOONEY_RIVLIN_SUBTYPE,CMFE_EQUATIONS_SET_NEARLY_INCOMPRESSIBLE_MOONEY_RIVLIN_SUBTYPE, &
     & CMFE_EQUATIONS_SET_MOONEY_RIVLIN_SUBTYPE, &
     & CMFE_EQUATIONS_SET_REFERENCE_STATE_MOONEY_RIVLIN_SUBTYPE, CMFE_EQUATIONS_SET_ISOTROPIC_EXPONENTIAL_SUBTYPE, &
@@ -5979,7 +5985,8 @@ MODULE OpenCMISS_Iron
     & PROBLEM_GROWTH_FINITE_ELASTICITY_NAVIER_STOKES_ALE_SUBTYPE !<Coupled Finite Elasticity with growth Navier Stokes moving mesh subtype \see OpenCMISS_ProblemSubtypes,OpenCMISS
   INTEGER(INTG), PARAMETER :: CMFE_PROBLEM_GROWTH_FINITE_ELAST_RBS_NAVIER_STOKES_ALE_SUBTYPE = &
     & PROBLEM_GROWTH_FINITE_ELASTICITY_RBS_NAVIER_STOKES_ALE_SUBTYPE !<Coupled Finite Elasticity with growth RBS Navier Stokes moving mesh subtype \see OpenCMISS_ProblemSubtypes,OpenCMISS
-
+  INTEGER(INTG), PARAMETER :: CMFE_PROBLEM_STOKES_DAMPING_LINEAR_ELASTICITY_SUBTYPE = &
+    & PROBLEM_STOKES_DAMPING_LINEAR_ELASTICITY_SUBTYPE !<Dynamica Linear Elasticity subtype \see OpenCMISS_ProblemSubtypes,OpenCMISS
   INTEGER(INTG), PARAMETER :: CMFE_PROBLEM_STATIC_FINITE_ELASTICITY_SUBTYPE = PROBLEM_STATIC_FINITE_ELASTICITY_SUBTYPE !<Static finite elasticity subtype \see OpenCMISS_ProblemSubtypes,OpenCMISS
   INTEGER(INTG), PARAMETER :: CMFE_PROBLEM_QUASISTATIC_FINITE_ELASTICITY_SUBTYPE = PROBLEM_QUASISTATIC_FINITE_ELASTICITY_SUBTYPE !<Quasistatic finite elasticity subtype \see OpenCMISS_ProblemSubtypes,OpenCMISS
   INTEGER(INTG), PARAMETER :: CMFE_PROBLEM_QUASISTATIC_FINITE_ELASTICITY_WITH_GROWTH_SUBTYPE = &
@@ -6143,6 +6150,8 @@ MODULE OpenCMISS_Iron
 
   PUBLIC CMFE_PROBLEM_FINITE_ELASTICITY_WITH_ACTIVE_SUBTYPE,CMFE_PROBLEM_FINITE_ELASTICITY_WITH_CELLML_SUBTYPE, &
     & CMFE_PROBLEM_FINITE_ELASTICITY_WITH_GROWTH_CELLML_SUBTYPE
+  
+  PUBLIC CMFE_PROBLEM_STOKES_DAMPING_LINEAR_ELASTICITY_SUBTYPE
 
 !!==================================================================================================================================
 !!
diff --git a/src/problem_constants.F90 b/src/problem_constants.F90
index 31ae1d7b..388469c3 100644
--- a/src/problem_constants.F90
+++ b/src/problem_constants.F90
@@ -112,7 +112,8 @@ MODULE PROBLEM_CONSTANTS
   !Problem subtypes
   INTEGER(INTG), PARAMETER :: PROBLEM_NO_SUBTYPE=0
   !Elasticity class
-  !  Linear elasticity - uses PROBLEM_NO_SUBTYPE=0
+  !  Linear elasticity
+  INTEGER(INTG), PARAMETER :: PROBLEM_STOKES_DAMPING_LINEAR_ELASTICITY_SUBTYPE=1
   !  Finite elasticity
   INTEGER(INTG), PARAMETER :: PROBLEM_STATIC_FINITE_ELASTICITY_SUBTYPE=1
   INTEGER(INTG), PARAMETER :: PROBLEM_QUASISTATIC_FINITE_ELASTICITY_SUBTYPE=2
diff --git a/src/problem_routines.F90 b/src/problem_routines.F90
index 5121a2ea..dae71fe8 100644
--- a/src/problem_routines.F90
+++ b/src/problem_routines.F90
@@ -461,9 +461,9 @@ RECURSIVE SUBROUTINE PROBLEM_CONTROL_LOOP_SOLVE(CONTROL_LOOP,err,error,*)
             
             !Precompute the number of iterations from total time span and time increment if it was not specified explicitely 
             IF(TIME_LOOP%NUMBER_OF_ITERATIONS==0) THEN
-              TIME_LOOP%NUMBER_OF_ITERATIONS=CEILING((TIME_LOOP%STOP_TIME-TIME_LOOP%START_TIME)/TIME_LOOP%TIME_INCREMENT)
+              TIME_LOOP%NUMBER_OF_ITERATIONS=NINT((TIME_LOOP%STOP_TIME-TIME_LOOP%START_TIME)/TIME_LOOP%TIME_INCREMENT)
             !If number of iterations was specified but does not match TIME_INCREMENT, e.g. TIME_INCREMENT is still at the default value, compute correct TIME_INCREMENT
-            ELSE IF(CEILING((TIME_LOOP%STOP_TIME-TIME_LOOP%START_TIME)/TIME_LOOP%TIME_INCREMENT) &
+            ELSE IF(NINT((TIME_LOOP%STOP_TIME-TIME_LOOP%START_TIME)/TIME_LOOP%TIME_INCREMENT) &
               & /= TIME_LOOP%NUMBER_OF_ITERATIONS) THEN
               TIME_LOOP%TIME_INCREMENT = (TIME_LOOP%STOP_TIME-TIME_LOOP%START_TIME)/TIME_LOOP%NUMBER_OF_ITERATIONS
             ENDIF
diff --git a/src/reaction_diffusion_equation_routines.F90 b/src/reaction_diffusion_equation_routines.F90
index 2f2af265..0e88f8bd 100644
--- a/src/reaction_diffusion_equation_routines.F90
+++ b/src/reaction_diffusion_equation_routines.F90
@@ -52,6 +52,7 @@ MODULE REACTION_DIFFUSION_EQUATION_ROUTINES
   USE Constants
   USE CONTROL_LOOP_ROUTINES
   USE ControlLoopAccessRoutines
+  USE COORDINATE_ROUTINES
   USE DistributedMatrixVector
   USE DOMAIN_MAPPINGS
   USE EquationsRoutines
@@ -66,6 +67,7 @@ MODULE REACTION_DIFFUSION_EQUATION_ROUTINES
   USE ISO_VARYING_STRING
   USE Kinds
   USE MatrixVector
+  USE Maths
 #ifndef NOMPIMOD
   USE MPI
 #endif
@@ -114,6 +116,8 @@ MODULE REACTION_DIFFUSION_EQUATION_ROUTINES
 
   PUBLIC REACTION_DIFFUSION_CONTROL_LOOP_POST_LOOP
 
+  PUBLIC ReactionDiffusion_GaussDeformationGradientTensor
+
 
 CONTAINS
 
@@ -132,6 +136,7 @@ SUBROUTINE ReactionDiffusion_EquationsSetSetup(EQUATIONS_SET,EQUATIONS_SET_SETUP
     !Local Variables
     INTEGER(INTG) :: component_idx,DIMENSION_MULTIPLIER,GEOMETRIC_COMPONENT_NUMBER,GEOMETRIC_SCALING_TYPE, &
       & NUMBER_OF_DIMENSIONS,NUMBER_OF_MATERIALS_COMPONENTS,GEOMETRIC_MESH_COMPONENT
+    INTEGER(INTG) :: INDEPENDENT_FIELD_NUMBER_OF_VARIABLES,INDEPENDENT_FIELD_NUMBER_OF_COMPONENTS
     TYPE(DECOMPOSITION_TYPE), POINTER :: GEOMETRIC_DECOMPOSITION
     TYPE(EquationsType), POINTER :: equations
     TYPE(EquationsMappingVectorType), POINTER :: vectorMapping
@@ -170,6 +175,9 @@ SUBROUTINE ReactionDiffusion_EquationsSetSetup(EQUATIONS_SET,EQUATIONS_SET_SETUP
         END SELECT
       CASE(EQUATIONS_SET_SETUP_GEOMETRY_TYPE)
         !\todo Check geometric dimension
+      !-----------------------------------------------------------------
+      ! D e p e n d e n t   f i e l d
+      !-----------------------------------------------------------------
       CASE(EQUATIONS_SET_SETUP_DEPENDENT_TYPE)
         SELECT CASE(EQUATIONS_SET_SETUP%ACTION_TYPE)
         CASE(EQUATIONS_SET_SETUP_START_ACTION)
@@ -288,6 +296,9 @@ SUBROUTINE ReactionDiffusion_EquationsSetSetup(EQUATIONS_SET,EQUATIONS_SET_SETUP
             & " is invalid for a reaction diffusion equation"
           CALL FlagError(localError,err,error,*999)
         END SELECT
+      !-----------------------------------------------------------------
+      ! M a t e r i a l   f i e l d
+      !----------------------------------------------------------------- 
       CASE(EQUATIONS_SET_SETUP_MATERIALS_TYPE)
         SELECT CASE(EQUATIONS_SET_SETUP%ACTION_TYPE)
         CASE(EQUATIONS_SET_SETUP_START_ACTION)
@@ -335,7 +346,7 @@ SUBROUTINE ReactionDiffusion_EquationsSetSetup(EQUATIONS_SET,EQUATIONS_SET_SETUP
                   CALL FIELD_COMPONENT_MESH_COMPONENT_SET(EQUATIONS_MATERIALS%MATERIALS_FIELD,FIELD_U_VARIABLE_TYPE, &
                     & component_idx,GEOMETRIC_COMPONENT_NUMBER,err,error,*999)
                   CALL FIELD_COMPONENT_INTERPOLATION_SET(EQUATIONS_MATERIALS%MATERIALS_FIELD,FIELD_U_VARIABLE_TYPE, &
-                    & component_idx,FIELD_CONSTANT_INTERPOLATION,err,error,*999)
+                    & component_idx,FIELD_NODE_BASED_INTERPOLATION,err,error,*999)
                 ENDDO !components_idx
                 !Default the storage co-efficient to the first geometric interpolation setup with constant interpolation
                 IF(EQUATIONS_SET%SPECIFICATION(3)==EQUATIONS_SET_CELLML_REAC_SPLIT_REAC_DIFF_SUBTYPE) THEN
@@ -345,7 +356,7 @@ SUBROUTINE ReactionDiffusion_EquationsSetSetup(EQUATIONS_SET,EQUATIONS_SET_SETUP
                   CALL FIELD_COMPONENT_MESH_COMPONENT_SET(EQUATIONS_MATERIALS%MATERIALS_FIELD,FIELD_U_VARIABLE_TYPE, &
                     & component_idx,GEOMETRIC_COMPONENT_NUMBER,err,error,*999)
                   CALL FIELD_COMPONENT_INTERPOLATION_SET(EQUATIONS_MATERIALS%MATERIALS_FIELD,FIELD_U_VARIABLE_TYPE, &
-                    & component_idx,FIELD_CONSTANT_INTERPOLATION,err,error,*999)
+                    & component_idx,FIELD_NODE_BASED_INTERPOLATION,err,error,*999)
                 ENDIF
                 !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)
@@ -408,6 +419,9 @@ SUBROUTINE ReactionDiffusion_EquationsSetSetup(EQUATIONS_SET,EQUATIONS_SET_SETUP
             & " is invalid for a reaction diffusion equation."
           CALL FlagError(localError,err,error,*999)
         END SELECT
+      !-----------------------------------------------------------------
+      ! S o u r c e   f i e l d
+      !----------------------------------------------------------------- 
       CASE(EQUATIONS_SET_SETUP_SOURCE_TYPE)
         SELECT CASE(EQUATIONS_SET_SETUP%ACTION_TYPE)
         CASE(EQUATIONS_SET_SETUP_START_ACTION)
@@ -526,6 +540,124 @@ SUBROUTINE ReactionDiffusion_EquationsSetSetup(EQUATIONS_SET,EQUATIONS_SET_SETUP
             & " is invalid for a reaction diffusion equation."
           CALL FlagError(localError,err,error,*999)
         END SELECT
+      !-----------------------------------------------------------------
+      ! I n d e p e n d e n t   f i e l d
+      !-----------------------------------------------------------------
+      CASE(EQUATIONS_SET_SETUP_INDEPENDENT_TYPE)
+        SELECT CASE(EQUATIONS_SET_SETUP%ACTION_TYPE)
+        !Set start action
+        CASE(EQUATIONS_SET_SETUP_START_ACTION)
+          INDEPENDENT_FIELD_NUMBER_OF_VARIABLES=5
+          !Create the auto created independent field
+          IF(EQUATIONS_SET%INDEPENDENT%INDEPENDENT_FIELD_AUTO_CREATED) THEN
+            !start field creation with name 'INDEPENDENT_FIELD'
+            CALL FIELD_CREATE_START(EQUATIONS_SET_SETUP%FIELD_USER_NUMBER,EQUATIONS_SET%REGION, &
+              & EQUATIONS_SET%INDEPENDENT%INDEPENDENT_FIELD,err,error,*999)
+            !start creation of a new field
+            CALL FIELD_TYPE_SET_AND_LOCK(EQUATIONS_SET%INDEPENDENT%INDEPENDENT_FIELD,FIELD_GEOMETRIC_GENERAL_TYPE,err,error,*999)
+            !label the field
+            CALL FIELD_LABEL_SET(EQUATIONS_SET%INDEPENDENT%INDEPENDENT_FIELD,"Independent Field",err,error,*999)
+            !define new created field to be independent
+            CALL FIELD_DEPENDENT_TYPE_SET_AND_LOCK(EQUATIONS_SET%INDEPENDENT%INDEPENDENT_FIELD, &
+              & FIELD_INDEPENDENT_TYPE,err,error,*999)
+            !look for decomposition rule already defined
+            CALL FIELD_MESH_DECOMPOSITION_GET(EQUATIONS_SET%GEOMETRY%GEOMETRIC_FIELD,GEOMETRIC_DECOMPOSITION,err,error,*999)
+            !apply decomposition rule found on new created field
+            CALL FIELD_MESH_DECOMPOSITION_SET_AND_LOCK(EQUATIONS_SET%INDEPENDENT%INDEPENDENT_FIELD, &
+              & GEOMETRIC_DECOMPOSITION,err,error,*999)
+            !point new field to geometric field
+            CALL FIELD_GEOMETRIC_FIELD_SET_AND_LOCK(EQUATIONS_SET%INDEPENDENT%INDEPENDENT_FIELD,EQUATIONS_SET% & 
+              & GEOMETRY%GEOMETRIC_FIELD,err,error,*999)
+            !set number of variables to 1 (1 for U)
+            CALL FIELD_NUMBER_OF_VARIABLES_SET_AND_LOCK(EQUATIONS_SET%INDEPENDENT%INDEPENDENT_FIELD, &
+              & INDEPENDENT_FIELD_NUMBER_OF_VARIABLES,err,error,*999)
+            CALL FIELD_VARIABLE_TYPES_SET_AND_LOCK(EQUATIONS_SET%INDEPENDENT%INDEPENDENT_FIELD,[FIELD_U_VARIABLE_TYPE, &
+                 & FIELD_DELUDELN_VARIABLE_TYPE, FIELD_V_VARIABLE_TYPE, &
+                 & FIELD_U1_VARIABLE_TYPE, FIELD_U2_VARIABLE_TYPE],err,error,*999)
+            !Velocity Field
+            CALL FIELD_DIMENSION_SET_AND_LOCK(EQUATIONS_SET%INDEPENDENT%INDEPENDENT_FIELD,FIELD_U_VARIABLE_TYPE, &
+              & FIELD_VECTOR_DIMENSION_TYPE,err,error,*999)
+            !Velocity Gradient Field
+            CALL FIELD_DIMENSION_SET_AND_LOCK(EQUATIONS_SET%INDEPENDENT%INDEPENDENT_FIELD,FIELD_DELUDELN_VARIABLE_TYPE, &
+              & FIELD_VECTOR_DIMENSION_TYPE,err,error,*999)
+            !TODO Temperary Field to store time step. Need to delete once time step is accessed in Finite Element Calculate
+            CALL FIELD_DIMENSION_SET_AND_LOCK(EQUATIONS_SET%INDEPENDENT%INDEPENDENT_FIELD,FIELD_V_VARIABLE_TYPE, &
+              & FIELD_VECTOR_DIMENSION_TYPE,err,error,*999)
+            !Stores Deformation Field of Previous Time Step, Required for Temporal Derivative of Jacobian
+            CALL FIELD_DIMENSION_SET_AND_LOCK(EQUATIONS_SET%INDEPENDENT%INDEPENDENT_FIELD,FIELD_U1_VARIABLE_TYPE, &
+              & FIELD_VECTOR_DIMENSION_TYPE,err,error,*999)
+            !Stores Deformation Field of Current Time Step
+            CALL FIELD_DIMENSION_SET_AND_LOCK(EQUATIONS_SET%INDEPENDENT%INDEPENDENT_FIELD,FIELD_U2_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_DATA_TYPE_SET_AND_LOCK(EQUATIONS_SET%INDEPENDENT%INDEPENDENT_FIELD,FIELD_DELUDELN_VARIABLE_TYPE, &
+              & FIELD_DP_TYPE,err,error,*999)
+            CALL FIELD_DATA_TYPE_SET_AND_LOCK(EQUATIONS_SET%INDEPENDENT%INDEPENDENT_FIELD,FIELD_V_VARIABLE_TYPE, &
+              & FIELD_DP_TYPE,err,error,*999)
+            CALL FIELD_DATA_TYPE_SET_AND_LOCK(EQUATIONS_SET%INDEPENDENT%INDEPENDENT_FIELD,FIELD_U1_VARIABLE_TYPE, &
+              & FIELD_DP_TYPE,err,error,*999)
+            CALL FIELD_DATA_TYPE_SET_AND_LOCK(EQUATIONS_SET%INDEPENDENT%INDEPENDENT_FIELD,FIELD_U2_VARIABLE_TYPE, &
+              & FIELD_DP_TYPE,err,error,*999)            
+            !calculate number of components with one component for each dimension
+            CALL FIELD_NUMBER_OF_COMPONENTS_GET(EQUATIONS_SET%GEOMETRY%GEOMETRIC_FIELD,FIELD_U_VARIABLE_TYPE, &
+              & NUMBER_OF_DIMENSIONS,err,error,*999)
+            INDEPENDENT_FIELD_NUMBER_OF_COMPONENTS=NUMBER_OF_DIMENSIONS
+            CALL FIELD_NUMBER_OF_COMPONENTS_SET_AND_LOCK(EQUATIONS_SET%INDEPENDENT%INDEPENDENT_FIELD, & 
+              & FIELD_U_VARIABLE_TYPE,INDEPENDENT_FIELD_NUMBER_OF_COMPONENTS,err,error,*999)
+            CALL FIELD_NUMBER_OF_COMPONENTS_SET_AND_LOCK(EQUATIONS_SET%INDEPENDENT%INDEPENDENT_FIELD, & 
+              & FIELD_DELUDELN_VARIABLE_TYPE,INDEPENDENT_FIELD_NUMBER_OF_COMPONENTS,err,error,*999)
+            CALL FIELD_NUMBER_OF_COMPONENTS_SET_AND_LOCK(EQUATIONS_SET%INDEPENDENT%INDEPENDENT_FIELD, & 
+              & FIELD_V_VARIABLE_TYPE,INDEPENDENT_FIELD_NUMBER_OF_COMPONENTS,err,error,*999)
+            CALL FIELD_NUMBER_OF_COMPONENTS_SET_AND_LOCK(EQUATIONS_SET%INDEPENDENT%INDEPENDENT_FIELD, & 
+              & FIELD_U1_VARIABLE_TYPE,INDEPENDENT_FIELD_NUMBER_OF_COMPONENTS,err,error,*999)
+            CALL FIELD_NUMBER_OF_COMPONENTS_SET_AND_LOCK(EQUATIONS_SET%INDEPENDENT%INDEPENDENT_FIELD, & 
+              & FIELD_U2_VARIABLE_TYPE,INDEPENDENT_FIELD_NUMBER_OF_COMPONENTS,err,error,*999)            
+            DO component_idx=1,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_SET%INDEPENDENT%INDEPENDENT_FIELD,FIELD_U_VARIABLE_TYPE, &
+                 & component_idx,GEOMETRIC_MESH_COMPONENT,err,error,*999)
+              CALL FIELD_COMPONENT_MESH_COMPONENT_SET(EQUATIONS_SET%INDEPENDENT%INDEPENDENT_FIELD,FIELD_DELUDELN_VARIABLE_TYPE, &
+                & component_idx,GEOMETRIC_MESH_COMPONENT,err,error,*999)
+              CALL FIELD_COMPONENT_MESH_COMPONENT_SET(EQUATIONS_SET%INDEPENDENT%INDEPENDENT_FIELD,FIELD_V_VARIABLE_TYPE, &
+                 & component_idx,GEOMETRIC_MESH_COMPONENT,err,error,*999)
+              CALL FIELD_COMPONENT_MESH_COMPONENT_SET(EQUATIONS_SET%INDEPENDENT%INDEPENDENT_FIELD,FIELD_U1_VARIABLE_TYPE, &
+                 & component_idx,GEOMETRIC_MESH_COMPONENT,err,error,*999)
+              CALL FIELD_COMPONENT_MESH_COMPONENT_SET(EQUATIONS_SET%INDEPENDENT%INDEPENDENT_FIELD,FIELD_U2_VARIABLE_TYPE, &
+                 & component_idx,GEOMETRIC_MESH_COMPONENT,err,error,*999)
+              CALL FIELD_COMPONENT_INTERPOLATION_SET_AND_LOCK(EQUATIONS_SET%INDEPENDENT%INDEPENDENT_FIELD, &
+                & FIELD_U_VARIABLE_TYPE,component_idx,FIELD_NODE_BASED_INTERPOLATION,err,error,*999)
+              CALL FIELD_COMPONENT_INTERPOLATION_SET_AND_LOCK(EQUATIONS_SET%INDEPENDENT%INDEPENDENT_FIELD, &
+                & FIELD_DELUDELN_VARIABLE_TYPE,component_idx,FIELD_NODE_BASED_INTERPOLATION,err,error,*999)
+              CALL FIELD_COMPONENT_INTERPOLATION_SET_AND_LOCK(EQUATIONS_SET%INDEPENDENT%INDEPENDENT_FIELD, &
+                & FIELD_V_VARIABLE_TYPE,component_idx,FIELD_NODE_BASED_INTERPOLATION,err,error,*999)
+              CALL FIELD_COMPONENT_INTERPOLATION_SET_AND_LOCK(EQUATIONS_SET%INDEPENDENT%INDEPENDENT_FIELD, &
+                & FIELD_U1_VARIABLE_TYPE,component_idx,FIELD_NODE_BASED_INTERPOLATION,err,error,*999)
+              CALL FIELD_COMPONENT_INTERPOLATION_SET_AND_LOCK(EQUATIONS_SET%INDEPENDENT%INDEPENDENT_FIELD, &
+                & FIELD_U2_VARIABLE_TYPE,component_idx,FIELD_NODE_BASED_INTERPOLATION,err,error,*999)
+            ENDDO !component_idx
+              !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_SET%INDEPENDENT%INDEPENDENT_FIELD,GEOMETRIC_SCALING_TYPE,err,error,*999)
+          ELSE
+            !Check the user specified field- Characteristic equation
+            CALL FIELD_TYPE_CHECK(EQUATIONS_SET_SETUP%FIELD,FIELD_GEOMETRIC_GENERAL_TYPE,err,error,*999)
+            CALL FIELD_DIMENSION_CHECK(EQUATIONS_SET_SETUP%FIELD,FIELD_U_VARIABLE_TYPE,FIELD_VECTOR_DIMENSION_TYPE, &
+              & err,error,*999)
+            CALL FIELD_DATA_TYPE_CHECK(EQUATIONS_SET_SETUP%FIELD,FIELD_U_VARIABLE_TYPE,FIELD_DP_TYPE,err,error,*999)
+          ENDIF    
+        !Specify finish action
+        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)
+          ENDIF
+        CASE DEFAULT
+          localError="The action type of "//TRIM(NumberToVString(EQUATIONS_SET_SETUP%ACTION_TYPE,"*",err,error))// &
+            & " for a setup type of "//TRIM(NumberToVString(EQUATIONS_SET_SETUP%SETUP_TYPE,"*",err,error))// &
+            & " is invalid for a reaction diffusion equation."
+          CALL FlagError(localError,err,error,*999)
+        END SELECT
       CASE(EQUATIONS_SET_SETUP_ANALYTIC_TYPE)
         SELECT CASE(EQUATIONS_SET_SETUP%ACTION_TYPE)
         CASE(EQUATIONS_SET_SETUP_START_ACTION)
@@ -780,8 +912,14 @@ SUBROUTINE ReactionDiffusion_FiniteElementCalculate(EQUATIONS_SET,ELEMENT_NUMBER
     !Local Variables
     INTEGER(INTG) FIELD_VAR_TYPE,mh,mhs,ms,ng,nh,nhs,ni,nj,ns,component_idx
     LOGICAL :: USE_FIBRES
-    REAL(DP) :: DIFFUSIVITY(3,3),DPHIDX(3,64),RWG,SUM,STORAGE_COEFFICIENT,C_PARAM
+    REAL(DP) :: DIFFUSIVITY(3,3),DPHIDX(3,64),PHI(3,64),DZDNU(3,3),PREVDZDNU(3,3),RWG,SUM,STORAGE_COEFFICIENT,C_PARAM
+    REAL(DP) :: Timestep, Jznu, prevJznu, dJdt, VELOCITY, VELOCITYGRAD, Jznu_ref, prevJznu_ref
+    REAL(DP) :: DZDNUT(3,3),DZDNUT_DZDNU(3,3),DZDNUT_DZDNU_INV(3,3),prevDZDNUT(3,3),prevDZDNUT_DZDNU(3,3),prevDZDNUT_DZDNU_INV(3,3)
     TYPE(BASIS_TYPE), POINTER :: DEPENDENT_BASIS,GEOMETRIC_BASIS,FIBRE_BASIS
+    TYPE(FIELD_INTERPOLATED_POINT_TYPE), POINTER :: geometricInterpPoint, prevIndependentInterpPoint,independentInterpPoint
+    TYPE(FIELD_INTERPOLATED_POINT_TYPE), POINTER :: fibreInterpolatedPoint
+    TYPE(FIELD_INTERPOLATED_POINT_METRICS_TYPE), POINTER :: geometricInterpPointMetrics
+    TYPE(FIELD_INTERPOLATED_POINT_METRICS_TYPE), POINTER :: prevIndependentInterpPointMetrics, independentInterpPointMetrics    
     TYPE(EquationsType), POINTER :: equations
     TYPE(EquationsMappingVectorType), POINTER :: vectorMapping
     TYPE(EquationsMappingDynamicType), POINTER :: dynamicMapping
@@ -791,11 +929,29 @@ SUBROUTINE ReactionDiffusion_FiniteElementCalculate(EQUATIONS_SET,ELEMENT_NUMBER
     TYPE(EquationsMatricesSourceType), POINTER :: sourceVector
     TYPE(EquationsMatrixType), POINTER :: dampingMatrix,stiffnessMatrix
     TYPE(EquationsVectorType), POINTER :: vectorEquations
-    TYPE(FIELD_TYPE), POINTER :: dependentField,geometricField,fibreField,materialsField,sourceField
+    TYPE(FIELD_TYPE), POINTER :: dependentField,geometricField,fibreField,materialsField,sourceField,independentField
     TYPE(FIELD_VARIABLE_TYPE), POINTER :: FIELD_VARIABLE,GEOMETRIC_VARIABLE
     TYPE(QUADRATURE_SCHEME_TYPE), POINTER :: QUADRATURE_SCHEME
-    
+    !TYPE(SOLVER_TYPE), POINTER :: solver !<A pointer to the solver to perform the pre-solve actions for.
+    !REAL(DP) :: currentTime,timeIncrement
+    !TYPE(CONTROL_LOOP_TYPE), POINTER :: controlLoop
+
     ENTERS("ReactionDiffusion_FiniteElementCalculate",err,error,*999)
+    
+    NULLIFY(fibreInterpolatedPoint)
+    NULLIFY(geometricInterpPoint)
+    NULLIFY(geometricInterpPointMetrics)
+    NULLIFY(independentInterpPoint)
+    NULLIFY(independentInterpPointMetrics)
+    NULLIFY(prevIndependentInterpPoint)
+    NULLIFY(prevIndependentInterpPointMetrics)
+
+    !TODO Find out how to access time step interval value
+    !IF(.NOT.ASSOCIATED(solver)) CALL FlagError("Solver is not associated.",err,error,*999)
+    !NULLIFY(controlLoop)
+    !CALL Solver_ControlLoopGet(solver,controlLoop,err,error,*999)
+    !CALL ControlLoop_CurrentTimesGet(controlLoop,currentTime,timeIncrement,err,error,*999)    
+    
 
     IF(ASSOCIATED(EQUATIONS_SET)) THEN
       IF(.NOT.ALLOCATED(EQUATIONS_SET%SPECIFICATION)) THEN
@@ -811,6 +967,7 @@ SUBROUTINE ReactionDiffusion_FiniteElementCalculate(EQUATIONS_SET,ELEMENT_NUMBER
         dependentField=>equations%interpolation%dependentField
         geometricField=>equations%interpolation%geometricField
         materialsField=>equations%interpolation%materialsField
+        independentField=>equations%interpolation%independentField
         IF(EQUATIONS_SET%SPECIFICATION(3)==EQUATIONS_SET_CELLML_REAC_NO_SPLIT_REAC_DIFF_SUBTYPE .OR. &
             & EQUATIONS_SET%SPECIFICATION(3)==EQUATIONS_SET_CONSTANT_REAC_DIFF_SUBTYPE) THEN
            sourceField=>equations%interpolation%sourceField
@@ -831,6 +988,16 @@ SUBROUTINE ReactionDiffusion_FiniteElementCalculate(EQUATIONS_SET,ELEMENT_NUMBER
           & geometricInterpParameters(FIELD_U_VARIABLE_TYPE)%ptr,err,error,*999)
         CALL FIELD_INTERPOLATION_PARAMETERS_ELEMENT_GET(FIELD_VALUES_SET_TYPE,ELEMENT_NUMBER,equations%interpolation% &
           & materialsInterpParameters(FIELD_U_VARIABLE_TYPE)%ptr,err,error,*999)
+        CALL FIELD_INTERPOLATION_PARAMETERS_ELEMENT_GET(FIELD_VALUES_SET_TYPE,ELEMENT_NUMBER,equations%interpolation% &
+          & independentInterpParameters(FIELD_U_VARIABLE_TYPE)%ptr,err,error,*999)
+        CALL FIELD_INTERPOLATION_PARAMETERS_ELEMENT_GET(FIELD_VALUES_SET_TYPE,ELEMENT_NUMBER,equations%interpolation% &
+          & independentInterpParameters(FIELD_DELUDELN_VARIABLE_TYPE)%ptr,err,error,*999)
+        CALL FIELD_INTERPOLATION_PARAMETERS_ELEMENT_GET(FIELD_VALUES_SET_TYPE,ELEMENT_NUMBER,equations%interpolation% &
+          & independentInterpParameters(FIELD_V_VARIABLE_TYPE)%ptr,err,error,*999)
+        CALL FIELD_INTERPOLATION_PARAMETERS_ELEMENT_GET(FIELD_VALUES_SET_TYPE,ELEMENT_NUMBER,equations%interpolation% &
+          & independentInterpParameters(FIELD_U1_VARIABLE_TYPE)%ptr,err,error,*999)
+        CALL FIELD_INTERPOLATION_PARAMETERS_ELEMENT_GET(FIELD_VALUES_SET_TYPE,ELEMENT_NUMBER,equations%interpolation% &
+          & independentInterpParameters(FIELD_U2_VARIABLE_TYPE)%ptr,err,error,*999)
         IF(USE_FIBRES) CALL FIELD_INTERPOLATION_PARAMETERS_ELEMENT_GET(FIELD_VALUES_SET_TYPE,ELEMENT_NUMBER,equations% &
           & interpolation%fibreInterpParameters(FIELD_U_VARIABLE_TYPE)%ptr,err,error,*999)
         IF(EQUATIONS_SET%SPECIFICATION(3)==EQUATIONS_SET_CELLML_REAC_NO_SPLIT_REAC_DIFF_SUBTYPE .OR. &
@@ -849,16 +1016,39 @@ SUBROUTINE ReactionDiffusion_FiniteElementCalculate(EQUATIONS_SET,ELEMENT_NUMBER
         dynamicMapping=>vectorMapping%dynamicMapping
         FIELD_VARIABLE=>dynamicMapping%equationsMatrixToVarMaps(1)%variable
         FIELD_VAR_TYPE=FIELD_VARIABLE%VARIABLE_TYPE
+        geometricInterpPoint=>equations%interpolation%geometricInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr
+        geometricInterpPointMetrics=>equations%interpolation%geometricInterpPointMetrics(FIELD_U_VARIABLE_TYPE)%ptr
+        prevIndependentInterpPoint=>equations%interpolation%independentInterpPoint(FIELD_U1_VARIABLE_TYPE)%ptr
+        independentInterpPoint=>equations%interpolation%independentInterpPoint(FIELD_U2_VARIABLE_TYPE)%ptr        
+        prevIndependentInterpPointMetrics=>equations%interpolation%independentInterpPointMetrics(FIELD_U1_VARIABLE_TYPE)%ptr
+        independentInterpPointMetrics=>equations%interpolation%independentInterpPointMetrics(FIELD_U2_VARIABLE_TYPE)%ptr
         IF(stiffnessMatrix%updateMatrix.OR.dampingMatrix%updateMatrix.OR.rhsVector%updateVector) THEN
           !Loop over gauss points
           DO ng=1,QUADRATURE_SCHEME%NUMBER_OF_GAUSS
             CALL FIELD_INTERPOLATE_GAUSS(FIRST_PART_DERIV,BASIS_DEFAULT_QUADRATURE_SCHEME,ng,equations%interpolation% &
               & geometricInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr,err,error,*999)
-            CALL FIELD_INTERPOLATED_POINT_METRICS_CALCULATE(GEOMETRIC_BASIS%NUMBER_OF_XI,equations%interpolation% &
-              & geometricInterpPointMetrics(FIELD_U_VARIABLE_TYPE)%ptr,err,error,*999)
             CALL FIELD_INTERPOLATE_GAUSS(NO_PART_DERIV,BASIS_DEFAULT_QUADRATURE_SCHEME,ng,equations%interpolation% &
               & materialsInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr,err,error,*999)
+            CALL FIELD_INTERPOLATE_GAUSS(NO_PART_DERIV,BASIS_DEFAULT_QUADRATURE_SCHEME,ng,equations%interpolation% &
+              & independentInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr,err,error,*999)
+            CALL FIELD_INTERPOLATE_GAUSS(NO_PART_DERIV,BASIS_DEFAULT_QUADRATURE_SCHEME,ng,equations%interpolation% &
+              & independentInterpPoint(FIELD_DELUDELN_VARIABLE_TYPE)%ptr,err,error,*999)
+            CALL FIELD_INTERPOLATE_GAUSS(NO_PART_DERIV,BASIS_DEFAULT_QUADRATURE_SCHEME,ng,equations%interpolation% &
+              & independentInterpPoint(FIELD_V_VARIABLE_TYPE)%ptr,err,error,*999)
+            CALL FIELD_INTERPOLATE_GAUSS(FIRST_PART_DERIV,BASIS_DEFAULT_QUADRATURE_SCHEME,ng,equations%interpolation% &
+              & independentInterpPoint(FIELD_U1_VARIABLE_TYPE)%ptr,err,error,*999)
+            CALL FIELD_INTERPOLATE_GAUSS(FIRST_PART_DERIV,BASIS_DEFAULT_QUADRATURE_SCHEME,ng,equations%interpolation% &
+              & independentInterpPoint(FIELD_U2_VARIABLE_TYPE)%ptr,err,error,*999)
+            CALL FIELD_INTERPOLATED_POINT_METRICS_CALCULATE(GEOMETRIC_BASIS%NUMBER_OF_XI, &
+              & geometricInterpPointMetrics,err,error,*999)
+            CALL FIELD_INTERPOLATED_POINT_METRICS_CALCULATE(GEOMETRIC_BASIS%NUMBER_OF_XI, &
+              & prevIndependentInterpPointMetrics,err,error,*999)
+            CALL FIELD_INTERPOLATED_POINT_METRICS_CALCULATE(GEOMETRIC_BASIS%NUMBER_OF_XI, &
+              & IndependentInterpPointMetrics,err,error,*999)
+      
+            
             IF(USE_FIBRES) THEN
+              fibreInterpolatedPoint=>equations%interpolation%fibreInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr
               CALL FIELD_INTERPOLATE_GAUSS(FIRST_PART_DERIV,BASIS_DEFAULT_QUADRATURE_SCHEME,ng,equations%interpolation% &
                 & fibreInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr,err,error,*999)
               CALL FIELD_INTERPOLATED_POINT_METRICS_CALCULATE(FIBRE_BASIS%NUMBER_OF_XI,equations%interpolation% &
@@ -869,6 +1059,27 @@ SUBROUTINE ReactionDiffusion_FiniteElementCalculate(EQUATIONS_SET,ELEMENT_NUMBER
               CALL FIELD_INTERPOLATE_GAUSS(NO_PART_DERIV,BASIS_DEFAULT_QUADRATURE_SCHEME,ng,equations%interpolation% &
                 & sourceInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr,err,error,*999)
             ENDIF
+            nj = GEOMETRIC_VARIABLE%NUMBER_OF_COMPONENTS
+            CALL ReactionDiffusion_GaussDeformationGradientTensor(independentInterpPointMetrics, &
+                 & geometricInterpPointMetrics, &
+                 & fibreInterpolatedPoint,dzdnu,ELEMENT_NUMBER,err,error,*999)
+            CALL MatrixTranspose(dZdNu(1:nj,1:nj), dZdNuT(1:nj,1:nj),err,error,*999)
+            CALL MatrixProduct(dZdNuT(1:nj,1:nj),dZdNu(1:nj,1:nj),dZdNuT_dZdNu(1:nj,1:nj),err,error,*999)
+            CALL Invert(dZdNuT_dZdNu(1:nj,1:nj),dZdNuT_dZdNu_Inv(1:nj,1:nj),Jznu,err,error,*999)
+            CALL Determinant(dZdNu(1:nj,1:nj),Jznu,err,error,*999)
+
+
+            !dZdNu and Jznu Previous
+            CALL ReactionDiffusion_GaussDeformationGradientTensor(prevIndependentInterpPointMetrics, &
+                   & geometricInterpPointMetrics, &
+                   & fibreInterpolatedPoint,prevdZdNU,ELEMENT_NUMBER,err,error,*999)
+            CALL Determinant(prevdZdNu(1:nj,1:nj),prevJznu,err,error,*999)
+           
+            !TODO This was a solution for 1D. Working on Generalizing this to 2D/3D using deformation Gradient Tensors
+            Jznu_ref = independentInterpPointMetrics%JACOBIAN/geometricInterpPointMetrics%JACOBIAN
+            prevJznu_ref = prevIndependentInterpPointMetrics%JACOBIAN/geometricInterpPointMetrics%JACOBIAN
+            Timestep=equations%interpolation%independentInterpPoint(FIELD_V_VARIABLE_TYPE)%ptr%VALUES(nj,1)
+            dJdt= ((Jznu - prevJznu) / (Timestep))
             !Calculate RWG.
             RWG=equations%interpolation%geometricInterpPointMetrics(FIELD_U_VARIABLE_TYPE)%ptr%JACOBIAN* &
               & QUADRATURE_SCHEME%GAUSS_WEIGHTS(ng)
@@ -876,7 +1087,10 @@ SUBROUTINE ReactionDiffusion_FiniteElementCalculate(EQUATIONS_SET,ELEMENT_NUMBER
             DIFFUSIVITY=0.0_DP
             IF(USE_FIBRES) THEN
               !Calculate the diffusivity tensor in fibre coordinates
-              CALL FlagError("Not implemented.",err,error,*999)
+              !TODO Implement
+              DO nj=1,GEOMETRIC_VARIABLE%NUMBER_OF_COMPONENTS !first three components of material field are the diffusivities
+                DIFFUSIVITY(nj,nj)=equations%interpolation%materialsInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr%VALUES(nj,1)
+              ENDDO !nj
             ELSE
               !Use the diffusivity tensor in geometric coordinates
               DO nj=1,GEOMETRIC_VARIABLE%NUMBER_OF_COMPONENTS !first three components of material field are the diffusivities
@@ -890,15 +1104,25 @@ SUBROUTINE ReactionDiffusion_FiniteElementCalculate(EQUATIONS_SET,ELEMENT_NUMBER
             DO nj=1,GEOMETRIC_VARIABLE%NUMBER_OF_COMPONENTS
               DO ms=1,DEPENDENT_BASIS%NUMBER_OF_ELEMENT_PARAMETERS
                 DPHIDX(nj,ms)=0.0_DP
+                PHI(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%geometricInterpPointMetrics(FIELD_U_VARIABLE_TYPE)%ptr%DXI_DX(ni,nj)
-                ENDDO !ni
+               ENDDO !ni
+               PHI(nj,ms)=PHI(nj,ms) + &
+                    & QUADRATURE_SCHEME%GAUSS_BASIS_FNS(ms,NO_PART_DERIV,ng)
               ENDDO !ms
             ENDDO !nj            
             !Loop over field components
-            mhs=0          
+            !WRITE(*,*) ng
+            !WRITE(*,*) "PHI mh 1, ng 1", QUADRATURE_SCHEME%GAUSS_BASIS_FNS(1,NO_PART_DERIV,1)
+            !WRITE(*,*) "PHI mh 2, ng 1", QUADRATURE_SCHEME%GAUSS_BASIS_FNS(2,NO_PART_DERIV,1)
+            !WRITE(*,*) "PHI mh 1, ng 2", QUADRATURE_SCHEME%GAUSS_BASIS_FNS(1,NO_PART_DERIV,2)
+            !WRITE(*,*) "PHI mh 2, ng 2", QUADRATURE_SCHEME%GAUSS_BASIS_FNS(2,NO_PART_DERIV,2)
+            mhs=0
+            stiffnessMatrix%updateMatrix = .TRUE.
+            dampingMatrix%updateMatrix = .TRUE.
             DO mh=1,FIELD_VARIABLE%NUMBER_OF_COMPONENTS
               !Loop over element rows
               DO ms=1,DEPENDENT_BASIS%NUMBER_OF_ELEMENT_PARAMETERS
@@ -910,11 +1134,23 @@ SUBROUTINE ReactionDiffusion_FiniteElementCalculate(EQUATIONS_SET,ELEMENT_NUMBER
                     nhs=nhs+1
                     SUM=0.0_DP
                     IF(stiffnessMatrix%updateMatrix) THEN
-                      DO ni=1,GEOMETRIC_VARIABLE%NUMBER_OF_COMPONENTS
-                        DO nj=1,GEOMETRIC_VARIABLE%NUMBER_OF_COMPONENTS
-                          SUM=SUM+DIFFUSIVITY(ni,nj)*DPHIDX(ni,mhs)*DPHIDX(nj,nhs)
+                       DO ni=1,GEOMETRIC_VARIABLE%NUMBER_OF_COMPONENTS
+                         DO nj=1,GEOMETRIC_VARIABLE%NUMBER_OF_COMPONENTS
+                           !Get the velocity (Used velocity for timestep, !TODO update names)
+                           !TODO Implement Advection in Separate Equation Set SubType
+                           VELOCITY=equations%interpolation%independentInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr%VALUES(nj,1)
+                           VELOCITYGRAD= &
+                             & equations%interpolation%independentInterpPoint(FIELD_DELUDELN_VARIABLE_TYPE)%ptr%VALUES(nj,1)
+                           SUM=SUM-VELOCITY*PHI(ni,mhs)*DPHIDX(nj,nhs)/(Jznu)
+
+                           != u * nabla(w), u = function and w = deformation velocity. Captures rate of expanding element
+                           SUM=SUM+DIFFUSIVITY(ni,nj)*DPHIDX(ni,mhs)*DPHIDX(nj,nhs)*dZdNuT_dZdNu_Inv(ni,nj)
+                     
                        ENDDO !nj
-                      ENDDO !ni
+                    ENDDO !ni
+                      SUM=SUM+dJdt*(1/Jznu)* &
+                           & QUADRATURE_SCHEME%GAUSS_BASIS_FNS(mhs,NO_PART_DERIV,ng)* &
+                           & QUADRATURE_SCHEME%GAUSS_BASIS_FNS(nhs,NO_PART_DERIV,ng)
                       stiffnessMatrix%elementMatrix%matrix(mhs,nhs)=stiffnessMatrix%elementMatrix%matrix(mhs,nhs)+(SUM*RWG)
                     ENDIF
                     IF(dampingMatrix%updateMatrix) THEN
@@ -1654,8 +1890,8 @@ SUBROUTINE REACTION_DIFFUSION_CONTROL_LOOP_POST_LOOP(CONTROL_LOOP,err,error,*)
                     dynamicMatrices=>vectorMatrices%dynamicMatrices
                     stiffnessMatrix=>dynamicMatrices%matrices(1)%ptr
                     dampingMatrix=>dynamicMatrices%matrices(2)%ptr
-                    stiffnessMatrix%updateMatrix = .FALSE.
-                    dampingMatrix%updateMatrix = .FALSE.
+                    stiffnessMatrix%updateMatrix = .TRUE.
+                    dampingMatrix%updateMatrix = .TRUE.
                   ELSE
                     CALL FlagError("Equations not associated.",err,error,*999)
                   ENDIF
@@ -1688,4 +1924,90 @@ SUBROUTINE REACTION_DIFFUSION_CONTROL_LOOP_POST_LOOP(CONTROL_LOOP,err,error,*)
 999 ERRORSEXITS("REACTION_DIFFUSION_CONTROL_LOOP_POST_LOOP",err,error)
     RETURN 1
   END SUBROUTINE REACTION_DIFFUSION_CONTROL_LOOP_POST_LOOP
+
+
+  
+  !>Evaluates the deformation gradient tensor at a given Gauss point
+  SUBROUTINE ReactionDiffusion_GaussDeformationGradientTensor(dependentInterpPointMetrics,geometricInterpPointMetrics,&
+    & fibreInterpolatedPoint,dZdNu,ELEM_NUM,err,error,*)
+
+    !Argument variables
+    TYPE(FIELD_INTERPOLATED_POINT_METRICS_TYPE), POINTER :: dependentInterpPointMetrics,geometricInterpPointMetrics
+    TYPE(FIELD_INTERPOLATED_POINT_TYPE), POINTER :: fibreInterpolatedPoint
+    REAL(DP), INTENT(OUT) :: dZdNu(3,3) !<dZdNu(coordinateIdx,coordianteIdx). On return, the deformation gradient tensor
+    INTEGER(INTG), INTENT(IN) :: ELEM_NUM
+    INTEGER(INTG), INTENT(OUT) :: err   !<The error code
+    TYPE(VARYING_STRING), INTENT(OUT) :: error !<The error string
+    !Local Variables
+    INTEGER(INTG) :: numberOfXDimensions,numberOfXiDimensions,numberOfZDimensions,i,j
+    REAL(DP) ::detdZdX,detdZdNu,dNudX(3,3),dXdNu(3,3),dNuDXi(3,3),dXidNu(3,3)
+
+    ENTERS("ReactionDiffusion_GaussDeformationGradientTensor",err,error,*999)
+    IF(ASSOCIATED(dependentInterpPointMetrics)) THEN
+      IF(ASSOCIATED(geometricInterpPointMetrics)) THEN
+        numberOfXDimensions=geometricInterpPointMetrics%NUMBER_OF_X_DIMENSIONS
+        numberOfXiDimensions=geometricInterpPointMetrics%NUMBER_OF_XI_DIMENSIONS
+        numberOfZDimensions=dependentInterpPointMetrics%NUMBER_OF_X_DIMENSIONS
+
+        CALL Coordinates_MaterialSystemCalculate(geometricInterpPointMetrics,fibreInterpolatedPoint,dNudX,dXdNu, &
+          & dNudXi(1:numberOfXDimensions,1:numberOfXiDimensions), &
+          & dXidNu(1:numberOfXiDimensions,1:numberOfXDimensions),err,error,*999)
+
+        !dZ/dNu = dZ/dXi * dXi/dNu  (deformation gradient tensor, F)
+        CALL MatrixProduct(dependentInterpPointMetrics%DX_DXI(1:numberOfZDimensions,1:numberOfXiDimensions), &
+          & dXiDNu(1:numberOfXiDimensions,1:numberOfXDimensions),dZdNu(1:numberOfZDimensions,1:numberOfXDimensions), &
+          & err,error,*999)
+
+        !TODO Debugging
+        !WRITE(*,*) "DZ_DXI", dependentInterpPointMetrics%DX_DXI(1:numberOfZDimensions,1:numberOfXiDimensions)
+        !WRITE(*,*) "DXI_DNU", dXiDNu(1:numberOfXiDimensions,1:numberOfXDimensions)
+        !WRITE(*,*) "DZ_DNU", dZdNu(1:numberOfZDimensions,1:numberOfXDimensions)
+        IF(numberOfZDimensions == 1) THEN
+          dZdNu(:,2) = [0.0_DP,1.0_DP,0.0_DP]
+          dZdNu(:,3) = [0.0_DP,0.0_DP,1.0_DP]
+          dZdNu(2:3,1) = 0.0_DP
+        ELSEIF(numberOfZDimensions == 2) THEN
+          dZdNu(:,3) = [0.0_DP,0.0_DP,1.0_DP]
+          dZdNu(3,1:2) = 0.0_DP
+        ENDIF
+
+        !IF(ELEM_NUM == 1) THEN
+           !WRITE(*,*) "dXi_dX", dXidNu(1:numberOfXiDimensions,1:numberOfXDimensions)
+           !WRITE(*,*) "dX_dXi", dNudXi(1:numberOfXDimensions,1:numberOfXiDimensions)
+           !WRITE(*,*) "dZ_dXi", dependentInterpPointMetrics%DX_DXI(1:numberOfZDimensions,1:numberOfXiDimensions)
+           !WRITE(*,*) "dZ_dX", dZdNu(1:numberOfZDimensions,1:numberOfXDimensions)
+        !ENDIF
+
+        IF(DIAGNOSTICS1) THEN
+          CALL WriteString(DIAGNOSTIC_OUTPUT_TYPE,"",err,error,*999)
+          CALL WriteString(DIAGNOSTIC_OUTPUT_TYPE,"Calculated deformation gradient tensor:",err,error,*999)
+          CALL WriteStringValue(DIAGNOSTIC_OUTPUT_TYPE,"  Number of Z dimensions  = ",numberOfZDimensions,err,error,*999)
+          CALL WriteStringValue(DIAGNOSTIC_OUTPUT_TYPE,"  Number of Xi dimensions = ",numberOfXiDimensions,err,error,*999)
+          CALL WriteString(DIAGNOSTIC_OUTPUT_TYPE,"  Derivative of X wrt to Nu coordinates:",err,error,*999)
+          CALL WriteStringMatrix(DIAGNOSTIC_OUTPUT_TYPE,1,1,numberOfXiDimensions,1,1,numberOfXDimensions, &
+            & numberOfXDimensions,numberOfXDimensions,dXidNu,WRITE_STRING_MATRIX_NAME_AND_INDICES, &
+            & '("    dX_dNu','(",I1,",:)','   :",3(X,E13.6))','(19X,3(X,E13.6))',err,error,*999)
+          CALL WriteString(DIAGNOSTIC_OUTPUT_TYPE,"  Deformation gradient tensor wrt Nu coordinates:",err,error,*999)
+          CALL WriteStringMatrix(DIAGNOSTIC_OUTPUT_TYPE,1,1,numberOfZDimensions,1,1,numberOfXDimensions, &
+            & numberOfXDimensions,numberOfXDimensions,dZdNu,WRITE_STRING_MATRIX_NAME_AND_INDICES, &
+            & '("    dZ_dNu','(",I1,",:)','   :",3(X,E13.6))','(19X,3(X,E13.6))',err,error,*999)
+          CALL Determinant(dZdNu,detdZdNu,err,error,*999)
+          CALL WriteStringValue(DIAGNOSTIC_OUTPUT_TYPE,"  Determinant dZ_dNu  = ",detdZdNu,err,error,*999)
+        ENDIF
+
+      ELSE
+        CALL FlagError("Geometric interpolated point metrics is not associated.",err,error,*999)
+      ENDIF
+    ELSE
+      CALL FlagError("Dependent interpolated point metrics is not associated.",err,error,*999)
+    ENDIF
+
+    EXITS("ReactionDiffusion_GaussDeformationGradientTensor")
+    RETURN
+999 ERRORS("ReactionDiffusion_GaussDeformationGradientTensor",err,error)
+    EXITS("ReactionDiffusion_GaussDeformationGradientTensor")
+    RETURN 1
+
+  END SUBROUTINE ReactionDiffusion_GaussDeformationGradientTensor
+
 END MODULE REACTION_DIFFUSION_EQUATION_ROUTINES