From d868a240b242948a88f414b2b6a56d6c694ef49f Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 2 Feb 2022 16:45:13 +0000 Subject: [PATCH] bugix: change of behavior --- PRIVATE | 2 +- ...hase_mechanical_eigen_thermalexpansion.f90 | 45 +++++++------------ src/phase_mechanical_elastic.f90 | 4 +- 3 files changed, 20 insertions(+), 31 deletions(-) diff --git a/PRIVATE b/PRIVATE index 5598ec489..419ba7b6f 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 5598ec4892b0921fccf63e8570f9fcd8e0365716 +Subproject commit 419ba7b6f6989a14e0c91c1f1f6621daf8215b71 diff --git a/src/phase_mechanical_eigen_thermalexpansion.f90 b/src/phase_mechanical_eigen_thermalexpansion.f90 index 70bf84cd8..a5d9868a8 100644 --- a/src/phase_mechanical_eigen_thermalexpansion.f90 +++ b/src/phase_mechanical_eigen_thermalexpansion.f90 @@ -8,10 +8,9 @@ submodule(phase:eigen) thermalexpansion integer, dimension(:), allocatable :: kinematics_thermal_expansion_instance type :: tParameters - real(pReal) :: & - T_ref - real(pReal), dimension(3,3,3) :: & - A = 0.0_pReal + type(tPolynomial) :: & + A_11, & + A_33 end type tParameters type(tParameters), dimension(:), allocatable :: param @@ -34,7 +33,7 @@ module function thermalexpansion_init(kinematics_length) result(myKinematics) phase, & mech, & kinematics, & - kinematic_type + myConfig print'(/,1x,a)', '<<<+- phase:mechanical:eigen:thermalexpansion init -+>>>' @@ -56,21 +55,13 @@ module function thermalexpansion_init(kinematics_length) result(myKinematics) do k = 1, kinematics%length if (myKinematics(k,p)) then associate(prm => param(kinematics_thermal_expansion_instance(p))) - kinematic_type => kinematics%get(k) - prm%T_ref = kinematic_type%get_asFloat('T_ref', defaultVal=T_ROOM) + myConfig => kinematics%get(k) + + prm%A_11 = polynomial(myConfig%asDict(),'A_11','T') + if (any(phase_lattice(p) == ['hP','tI'])) & + prm%A_33 = polynomial(myConfig%asDict(),'A_33','T') - prm%A(1,1,1) = kinematic_type%get_asFloat('A_11') - prm%A(1,1,2) = kinematic_type%get_asFloat('A_11,T', defaultVal=0.0_pReal) - prm%A(1,1,3) = kinematic_type%get_asFloat('A_11,T^2',defaultVal=0.0_pReal) - if (any(phase_lattice(p) == ['hP','tI'])) then - prm%A(3,3,1) = kinematic_type%get_asFloat('A_33') - prm%A(3,3,2) = kinematic_type%get_asFloat('A_33,T', defaultVal=0.0_pReal) - prm%A(3,3,3) = kinematic_type%get_asFloat('A_33,T^2',defaultVal=0.0_pReal) - end if - do i=1, size(prm%A,3) - prm%A(1:3,1:3,i) = lattice_symmetrize_33(prm%A(1:3,1:3,i),phase_lattice(p)) - end do end associate end if end do @@ -91,22 +82,20 @@ module subroutine thermalexpansion_LiAndItsTangent(Li, dLi_dTstar, ph,me) dLi_dTstar !< derivative of Li with respect to Tstar (4th-order tensor defined to be zero) real(pReal) :: T, dot_T + real(pReal), dimension(3,3) :: A T = thermal_T(ph,me) dot_T = thermal_dot_T(ph,me) associate(prm => param(kinematics_thermal_expansion_instance(ph))) - Li = dot_T * ( & - prm%A(1:3,1:3,1) & ! constant coefficient - + prm%A(1:3,1:3,2)*(T - prm%T_ref) & ! linear coefficient - + prm%A(1:3,1:3,3)*(T - prm%T_ref)**2 & ! quadratic coefficient - ) / & - (1.0_pReal & - + prm%A(1:3,1:3,1)*(T - prm%T_ref) / 1.0_pReal & - + prm%A(1:3,1:3,2)*(T - prm%T_ref)**2 / 2.0_pReal & - + prm%A(1:3,1:3,3)*(T - prm%T_ref)**3 / 3.0_pReal & - ) + + A = 0.0_pReal + A(1,1) = prm%A_11%at(T) + if (any(phase_lattice(ph) == ['hP','tI'])) A(3,3) = prm%A_33%at(T) + A = lattice_symmetrize_33(A,phase_lattice(ph)) + Li = dot_T * A + end associate dLi_dTstar = 0.0_pReal diff --git a/src/phase_mechanical_elastic.f90 b/src/phase_mechanical_elastic.f90 index c57849b1f..084cc75bc 100644 --- a/src/phase_mechanical_elastic.f90 +++ b/src/phase_mechanical_elastic.f90 @@ -48,7 +48,7 @@ module subroutine elastic_init(phases) prm%C_11 = polynomial(elastic%asDict(),'C_11','T') prm%C_12 = polynomial(elastic%asDict(),'C_12','T') prm%C_44 = polynomial(elastic%asDict(),'C_44','T') - + if (any(phase_lattice(ph) == ['hP','tI'])) then prm%C_13 = polynomial(elastic%asDict(),'C_13','T') prm%C_33 = polynomial(elastic%asDict(),'C_33','T') @@ -77,7 +77,7 @@ pure module function elastic_C66(ph,en) result(C66) associate(prm => param(ph)) - + C66 = 0.0_pReal T = thermal_T(ph,en)