should be part of eigendeformation submodule

src/commercialFEM_fileList.f90

@@ -32,10 +32,10 @@
 #include "phase_mechanics_eigendeformation.f90"
 #include "phase_mechanics_eigendeformation_cleavageopening.f90"
 #include "phase_mechanics_eigendeformation_slipplaneopening.f90"
+#include "phase_mechanics_eigendeformation_thermalexpansion.f90"
 #include "phase_thermal.f90"
 #include "phase_thermal_dissipation.f90"
 #include "phase_thermal_externalheat.f90"
-#include "phase_mechanics_eigendeformation_thermalexpansion.f90"
 #include "phase_damage.f90"
 #include "phase_damage_isobrittle.f90"
 #include "phase_damage_isoductile.f90"

src/phase.f90

@@ -186,6 +186,11 @@ module phase
       real(pReal) :: T
     end function thermal_T
 
+    module function thermal_dot_T(ph,me) result(dot_T)
+      integer, intent(in) :: ph,me
+      real(pReal) :: dot_T
+    end function thermal_dot_T
+
 
     module subroutine constitutive_mech_setF(F,co,ip,el)
       real(pReal), dimension(3,3), intent(in) :: F

src/phase_mechanics.f90

@@ -137,6 +137,11 @@ submodule(phase) mechanics
       logical, dimension(:,:), allocatable :: myKinematics
     end function kinematics_slipplane_opening_init
 
+    module function kinematics_thermal_expansion_init(kinematics_length) result(myKinematics)
+      integer, intent(in) :: kinematics_length
+      logical, dimension(:,:), allocatable :: myKinematics
+    end function kinematics_thermal_expansion_init
+
     module subroutine plastic_isotropic_results(instance,group)
       integer,          intent(in) :: instance
       character(len=*), intent(in) :: group
@@ -365,6 +370,7 @@ module subroutine mech_init(phases)
   if(maxval(phase_Nkinematics) /= 0) then
     where(kinematics_cleavage_opening_init(maxval(phase_Nkinematics)))  phase_kinematics = KINEMATICS_cleavage_opening_ID
     where(kinematics_slipplane_opening_init(maxval(phase_Nkinematics))) phase_kinematics = KINEMATICS_slipplane_opening_ID
+    where(kinematics_thermal_expansion_init(maxval(phase_Nkinematics))) phase_kinematics = KINEMATICS_thermal_expansion_ID
   endif
 
 end subroutine mech_init

src/phase_mechanics_eigendeformation_thermalexpansion.f90

@@ -3,7 +3,7 @@
 !> @brief material subroutine incorporating kinematics resulting from thermal expansion
 !> @details to be done
 !--------------------------------------------------------------------------------------------------
-submodule(phase:thermal) thermalexpansion
+submodule(phase:eigendeformation) thermalexpansion
 
   integer, dimension(:), allocatable :: kinematics_thermal_expansion_instance
 
@@ -93,21 +93,21 @@ module subroutine thermalexpansion_LiAndItsTangent(Li, dLi_dTstar, ph,me)
 
   real(pReal) :: T, dot_T
 
-  T     = current(ph)%T(me)
+  T     = thermal_T(ph,me)
-  dot_T = current(ph)%dot_T(me)
+  dot_T = thermal_dot_T(ph,me)
 
   associate(prm => param(kinematics_thermal_expansion_instance(ph)))
-  Li = dot_T * ( &
+    Li = dot_T * ( &
-                prm%A(1:3,1:3,1)*(T - prm%T_ref)**0 &                                               ! constant  coefficient
+                  prm%A(1:3,1:3,1)*(T - prm%T_ref)**0 &                                               ! constant  coefficient
-              + prm%A(1:3,1:3,2)*(T - prm%T_ref)**1 &                                               ! linear    coefficient
+                + prm%A(1:3,1:3,2)*(T - prm%T_ref)**1 &                                               ! linear    coefficient
-              + prm%A(1:3,1:3,3)*(T - prm%T_ref)**2 &                                               ! quadratic coefficient
+                + prm%A(1:3,1:3,3)*(T - prm%T_ref)**2 &                                               ! quadratic coefficient
-              ) / &
+                ) / &
-       (1.0_pReal &
+         (1.0_pReal &
-             + prm%A(1:3,1:3,1)*(T - prm%T_ref)**1 / 1. &
+               + prm%A(1:3,1:3,1)*(T - prm%T_ref)**1 / 1. &
-             + prm%A(1:3,1:3,2)*(T - prm%T_ref)**2 / 2. &
+               + prm%A(1:3,1:3,2)*(T - prm%T_ref)**2 / 2. &
-             + prm%A(1:3,1:3,3)*(T - prm%T_ref)**3 / 3. &
+               + prm%A(1:3,1:3,3)*(T - prm%T_ref)**3 / 3. &
-       )
+         )
-  end associate
+    end associate
   dLi_dTstar = 0.0_pReal
 
 end subroutine thermalexpansion_LiAndItsTangent

src/phase_thermal.f90

@@ -32,10 +32,7 @@ submodule(phase) thermal
       logical, dimension(:,:), allocatable :: mySources
     end function externalheat_init
 
-    module function kinematics_thermal_expansion_init(kinematics_length) result(myKinematics)
+
-      integer, intent(in) :: kinematics_length
-      logical, dimension(:,:), allocatable :: myKinematics
-    end function kinematics_thermal_expansion_init
 
 
     module subroutine externalheat_dotState(ph, me)
@@ -121,11 +118,6 @@ module subroutine thermal_init(phases)
                                           maxval(thermalState(ph)%p%sizeDotState))
   enddo PhaseLoop2
 
-!--------------------------------------------------------------------------------------------------
-!initialize kinematic mechanisms
-  if(maxval(phase_Nkinematics) /= 0) where(kinematics_thermal_expansion_init(maxval(phase_Nkinematics))) &
-                                           phase_kinematics = KINEMATICS_thermal_expansion_ID
-
 end subroutine thermal_init
 
 
@@ -257,6 +249,20 @@ module function thermal_T(ph,me) result(T)
 end function thermal_T
 
 
+!----------------------------------------------------------------------------------------------
+!< @brief Get rate of temperature (for use by non-thermal physics)
+!----------------------------------------------------------------------------------------------
+module function thermal_dot_T(ph,me) result(dot_T)
+
+  integer, intent(in) :: ph, me
+  real(pReal) :: dot_T
+
+
+  dot_T = current(ph)%dot_T(me)
+
+end function thermal_dot_T
+
+
 !----------------------------------------------------------------------------------------------
 !< @brief Set temperature
 !----------------------------------------------------------------------------------------------