diff --git a/src/homogenization_damage.f90 b/src/homogenization_damage.f90
index 8a4e596a8..b8aecccae 100644
--- a/src/homogenization_damage.f90
+++ b/src/homogenization_damage.f90
@@ -114,7 +114,7 @@ module function homogenization_mu_phi(ce) result(mu)
   real(pReal) :: mu
 
 
-  mu = lattice_mu_phi(material_phaseID(1,ce))
+  mu = phase_mu_phi(1,ce)
 
 end function homogenization_mu_phi
 
@@ -128,7 +128,7 @@ module function homogenization_K_phi(ce) result(K)
   real(pReal), dimension(3,3) :: K
 
 
-  K = crystallite_push33ToRef(1,ce,lattice_K_phi(1:3,1:3,material_phaseID(1,ce))) \
+  K = phase_K_phi(1,ce) &
     * num_damage%charLength**2
 
 end function homogenization_K_phi
diff --git a/src/homogenization_thermal.f90 b/src/homogenization_thermal.f90
index 27988dfb2..0effe09ed 100644
--- a/src/homogenization_thermal.f90
+++ b/src/homogenization_thermal.f90
@@ -102,8 +102,15 @@ module function homogenization_mu_T(ce) result(mu)
   integer, intent(in) :: ce
   real(pReal) :: mu
 
+  integer :: co
 
-  mu = c_P(ce) * rho(ce)
+
+  mu = phase_mu_T(1,ce)
+  do co = 2, homogenization_Nconstituents(material_homogenizationID(ce))
+    mu = mu + phase_mu_T(co,ce)
+  enddo
+
+  mu = mu / real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
 
 end function homogenization_mu_T
 
@@ -119,9 +126,9 @@ module function homogenization_K_T(ce) result(K)
   integer :: co
 
 
-  K = crystallite_push33ToRef(1,ce,lattice_K_T(:,:,material_phaseID(1,ce)))
+  K = phase_K_T(1,ce)
   do co = 2, homogenization_Nconstituents(material_homogenizationID(ce))
-    K = K + crystallite_push33ToRef(co,ce,lattice_K_T(:,:,material_phaseID(co,ce)))
+    K = K + phase_K_T(co,ce)
   enddo
 
   K = K / real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
@@ -187,46 +194,4 @@ module subroutine thermal_results(ho,group)
 
 end subroutine thermal_results
 
-
-!--------------------------------------------------------------------------------------------------
-!> @brief Homogenize specific heat capacity.
-!--------------------------------------------------------------------------------------------------
-function c_P(ce)
-
-  integer, intent(in) :: ce
-  real(pReal) :: c_P
-
-  integer :: co
-
-
-  c_P = lattice_c_p(material_phaseID(1,ce))
-  do co = 2, homogenization_Nconstituents(material_homogenizationID(ce))
-    c_P = c_P + lattice_c_p(material_phaseID(co,ce))
-  enddo
-
-  c_P = c_P / real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
-
-end function c_P
-
-
-!--------------------------------------------------------------------------------------------------
-!> @brief Homogenize mass density.
-!--------------------------------------------------------------------------------------------------
-function rho(ce)
-
-  integer, intent(in) :: ce
-  real(pReal) :: rho
-
-  integer :: co
-
-
-  rho = lattice_rho(material_phaseID(1,ce))
-  do co = 2, homogenization_Nconstituents(material_homogenizationID(ce))
-    rho = rho + lattice_rho(material_phaseID(co,ce))
-  enddo
-
-  rho = rho / real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
-
-end function rho
-
 end submodule thermal
diff --git a/src/phase.f90 b/src/phase.f90
index 346772ce4..54116f3c9 100644
--- a/src/phase.f90
+++ b/src/phase.f90
@@ -184,7 +184,7 @@ module phase
 
     module subroutine phase_thermal_setField(T,dot_T, co,ce)
       real(pReal), intent(in) :: T, dot_T
-      integer, intent(in) :: ce, co
+      integer, intent(in) :: co, ce
     end subroutine phase_thermal_setField
 
     module subroutine phase_set_phi(phi,co,ce)
@@ -192,6 +192,28 @@ module phase
       integer, intent(in) :: co, ce
     end subroutine phase_set_phi
 
+
+    module function phase_mu_phi(co,ce) result(mu)
+      integer, intent(in) :: co, ce
+      real(pReal) :: mu
+    end function phase_mu_phi
+
+    module function phase_K_phi(co,ce) result(K)
+      integer, intent(in) :: co, ce
+      real(pReal), dimension(3,3) :: K
+    end function phase_K_phi
+
+
+    module function phase_mu_T(co,ce) result(mu)
+      integer, intent(in) :: co, ce
+      real(pReal) :: mu
+    end function phase_mu_T
+
+    module function phase_K_T(co,ce) result(K)
+      integer, intent(in) :: co, ce
+      real(pReal), dimension(3,3) :: K
+    end function phase_K_T
+
 ! == cleaned:end ===================================================================================
 
     module function thermal_stress(Delta_t,ph,me) result(converged_)
@@ -297,6 +319,10 @@ module phase
     phase_homogenizedC, &
     phase_f_phi, &
     phase_f_T, &
+    phase_K_phi, &
+    phase_K_T, &
+    phase_mu_phi, &
+    phase_mu_T, &
     phase_results, &
     phase_allocateState, &
     phase_forward, &
diff --git a/src/phase_damage.f90 b/src/phase_damage.f90
index 5f6d65ab8..9a20030d2 100644
--- a/src/phase_damage.f90
+++ b/src/phase_damage.f90
@@ -345,6 +345,33 @@ function phase_damage_collectDotState(ph,me) result(broken)
 end function phase_damage_collectDotState
 
 
+!--------------------------------------------------------------------------------------------------
+!> @brief Damage viscosity.
+!--------------------------------------------------------------------------------------------------
+module function phase_mu_phi(co,ce) result(mu)
+
+  integer, intent(in) :: co, ce
+  real(pReal) :: mu
+
+
+  mu = lattice_mu_phi(material_phaseID(co,ce))
+
+end function phase_mu_phi
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief Damage conductivity/diffusivity in reference configuration.
+!--------------------------------------------------------------------------------------------------
+module function phase_K_phi(co,ce) result(K)
+
+  integer, intent(in) :: co, ce
+  real(pReal), dimension(3,3) :: K
+
+
+  K = crystallite_push33ToRef(co,ce,lattice_K_phi(1:3,1:3,material_phaseID(co,ce)))
+
+end function phase_K_phi
+
 
 !--------------------------------------------------------------------------------------------------
 !> @brief for constitutive models having an instantaneous change of state
diff --git a/src/phase_thermal.f90 b/src/phase_thermal.f90
index c2b4b34b2..a0e6adb60 100644
--- a/src/phase_thermal.f90
+++ b/src/phase_thermal.f90
@@ -184,6 +184,35 @@ function phase_thermal_collectDotState(ph,me) result(broken)
 end function phase_thermal_collectDotState
 
 
+!--------------------------------------------------------------------------------------------------
+!> @brief Damage viscosity.
+!--------------------------------------------------------------------------------------------------
+module function phase_mu_T(co,ce) result(mu)
+
+  integer, intent(in) :: co, ce
+  real(pReal) :: mu
+
+
+  mu = lattice_c_p(material_phaseID(co,ce)) &
+     * lattice_rho(material_phaseID(co,ce))
+
+end function phase_mu_T
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief Damage conductivity/diffusivity in reference configuration.
+!--------------------------------------------------------------------------------------------------
+module function phase_K_T(co,ce) result(K)
+
+  integer, intent(in) :: co, ce
+  real(pReal), dimension(3,3) :: K
+
+
+  K = crystallite_push33ToRef(co,ce,lattice_K_T(1:3,1:3,material_phaseID(co,ce)))
+
+end function phase_K_T
+
+
 module function thermal_stress(Delta_t,ph,me) result(converged_)           ! ?? why is this called "stress" when it seems closer to "updateState" ??
 
   real(pReal), intent(in) :: Delta_t