diff --git a/src/commercialFEM_fileList.f90 b/src/commercialFEM_fileList.f90
index 942363a9f..f02d7bc92 100644
--- a/src/commercialFEM_fileList.f90
+++ b/src/commercialFEM_fileList.f90
@@ -26,6 +26,17 @@
 #endif
 #include "material.f90"
 #include "lattice.f90"
+#include "constitutive.f90"
+#include "constitutive_plastic.f90"
+#include "constitutive_damage.f90"
+#include "constitutive_thermal.f90"
+#include "constitutive_plastic_none.f90"
+#include "constitutive_plastic_isotropic.f90"
+#include "constitutive_plastic_phenopowerlaw.f90"
+#include "constitutive_plastic_kinehardening.f90"
+#include "constitutive_plastic_dislotwin.f90"
+#include "constitutive_plastic_disloUCLA.f90"
+#include "constitutive_plastic_nonlocal.f90"
 #include "source_thermal_dissipation.f90"
 #include "source_thermal_externalheat.f90"
 #include "source_damage_isoBrittle.f90"
@@ -35,14 +46,6 @@
 #include "kinematics_cleavage_opening.f90"
 #include "kinematics_slipplane_opening.f90"
 #include "kinematics_thermal_expansion.f90"
-#include "constitutive.f90"
-#include "constitutive_plastic_none.f90"
-#include "constitutive_plastic_isotropic.f90"
-#include "constitutive_plastic_phenopowerlaw.f90"
-#include "constitutive_plastic_kinehardening.f90"
-#include "constitutive_plastic_dislotwin.f90"
-#include "constitutive_plastic_disloUCLA.f90"
-#include "constitutive_plastic_nonlocal.f90"
 #include "crystallite.f90"
 #include "thermal_isothermal.f90"
 #include "thermal_adiabatic.f90"
diff --git a/src/constitutive.f90 b/src/constitutive.f90
index 5f509b0ba..2af434df5 100644
--- a/src/constitutive.f90
+++ b/src/constitutive.f90
@@ -35,13 +35,16 @@ module constitutive
     module subroutine thermal_init
     end subroutine thermal_init
 
-!   module pure function kinematics_thermal_expansion_initialStrain(homog,phase,offset)
-!
-!     integer, intent(in) :: &
-!       phase, &
-!       homog, &
-!       offset
-!   end function kinematics_thermal_expansion_initialStrain
+    pure module function kinematics_thermal_expansion_initialStrain(homog,phase,offset) result(initialStrain)
+
+     integer, intent(in) :: &
+       phase, &
+       homog, &
+       offset
+     real(pReal), dimension(3,3) :: &
+       initialStrain
+
+    end function kinematics_thermal_expansion_initialStrain
 
     module function plastic_dislotwin_homogenizedC(ipc,ip,el) result(homogenizedC)
       real(pReal), dimension(6,6) :: &
@@ -591,7 +594,7 @@ pure function constitutive_initialFi(ipc, ip, el)
         homog = material_homogenizationAt(el)
         offset = thermalMapping(homog)%p(ip,el)
         constitutive_initialFi = &
-          constitutive_initialFi !+ kinematics_thermal_expansion_initialStrain(homog,phase,offset)
+          constitutive_initialFi + kinematics_thermal_expansion_initialStrain(homog,phase,offset)
     end select kinematicsType
   enddo KinematicsLoop
 
diff --git a/src/kinematics_thermal_expansion.f90 b/src/kinematics_thermal_expansion.f90
index 00871e94f..8fc2dc704 100644
--- a/src/kinematics_thermal_expansion.f90
+++ b/src/kinematics_thermal_expansion.f90
@@ -66,24 +66,24 @@ end subroutine kinematics_thermal_expansion_init
 !--------------------------------------------------------------------------------------------------
 !> @brief  report initial thermal strain based on current temperature deviation from reference
 !--------------------------------------------------------------------------------------------------
-!pure module function kinematics_thermal_expansion_initialStrain(homog,phase,offset)
-!
-! integer, intent(in) :: &
-!   phase, &
-!   homog, &
-!   offset
-!
-! real(pReal), dimension(3,3) :: &
-!   kinematics_thermal_expansion_initialStrain                                                      !< initial thermal strain (should be small strain, though)
-!
-! associate(prm => param(kinematics_thermal_expansion_instance(phase)))
-! kinematics_thermal_expansion_initialStrain = &
-!   (temperature(homog)%p(offset) - prm%T_ref)**1 / 1. * prm%expansion(1:3,1:3,1) + &               ! constant  coefficient
-!   (temperature(homog)%p(offset) - prm%T_ref)**2 / 2. * prm%expansion(1:3,1:3,2) + &               ! linear    coefficient
-!   (temperature(homog)%p(offset) - prm%T_ref)**3 / 3. * prm%expansion(1:3,1:3,3)                   ! quadratic coefficient
-! end associate
-!
-!end function kinematics_thermal_expansion_initialStrain
+pure module function kinematics_thermal_expansion_initialStrain(homog,phase,offset) result(initialStrain)
+
+ integer, intent(in) :: &
+   phase, &
+   homog, &
+   offset
+
+ real(pReal), dimension(3,3) :: &
+   initialStrain                                                                                   !< initial thermal strain (should be small strain, though)
+
+ associate(prm => param(kinematics_thermal_expansion_instance(phase)))
+ initialStrain = &
+   (temperature(homog)%p(offset) - prm%T_ref)**1 / 1. * prm%expansion(1:3,1:3,1) + &               ! constant  coefficient
+   (temperature(homog)%p(offset) - prm%T_ref)**2 / 2. * prm%expansion(1:3,1:3,2) + &               ! linear    coefficient
+   (temperature(homog)%p(offset) - prm%T_ref)**3 / 3. * prm%expansion(1:3,1:3,3)                   ! quadratic coefficient
+ end associate
+
+end function kinematics_thermal_expansion_initialStrain
 
 
 !--------------------------------------------------------------------------------------------------