diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 6c096ecd0..ef6004109 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -611,9 +611,9 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, S6, Fi, ipc, ip, e KinematicsLoop: do k = 1_pInt, phase_Nkinematics(material_phase(ipc,ip,el)) kinematicsType: select case (phase_kinematics(k,material_phase(ipc,ip,el))) case (KINEMATICS_cleavage_opening_ID) kinematicsType - call kinematics_cleavage_opening_LiAndItsTangent(my_Li, my_dLi_dS, S6, ipc, ip, el) + call kinematics_cleavage_opening_LiAndItsTangent(my_Li, my_dLi_dS, math_6toSym33(S6), ipc, ip, el) case (KINEMATICS_slipplane_opening_ID) kinematicsType - call kinematics_slipplane_opening_LiAndItsTangent(my_Li, my_dLi_dS, S6, ipc, ip, el) + call kinematics_slipplane_opening_LiAndItsTangent(my_Li, my_dLi_dS, math_6toSym33(S6), ipc, ip, el) case (KINEMATICS_thermal_expansion_ID) kinematicsType call kinematics_thermal_expansion_LiAndItsTangent(my_Li, my_dLi_dS, ipc, ip, el) case default kinematicsType @@ -901,7 +901,7 @@ subroutine constitutive_collectDotState(S6, FeArray, Fi, FpArray, subdt, subfrac sourceType: select case (phase_source(s,material_phase(ipc,ip,el))) case (SOURCE_damage_anisoBrittle_ID) sourceType - call source_damage_anisoBrittle_dotState (S6, ipc, ip, el) !< correct stress? + call source_damage_anisoBrittle_dotState (math_6toSym33(S6), ipc, ip, el) !< correct stress? case (SOURCE_damage_isoDuctile_ID) sourceType call source_damage_isoDuctile_dotState ( ipc, ip, el) diff --git a/src/kinematics_cleavage_opening.f90 b/src/kinematics_cleavage_opening.f90 index 89d9dcd68..89c2f6ff0 100644 --- a/src/kinematics_cleavage_opening.f90 +++ b/src/kinematics_cleavage_opening.f90 @@ -113,10 +113,10 @@ subroutine kinematics_cleavage_opening_init() tempInt = config_phase(p)%getInts('ncleavage') kinematics_cleavage_opening_Ncleavage(1:size(tempInt),instance) = tempInt - tempFloat = config_phase(p)%getFloats('anisobrittle_criticaldisplacement',requiredShape=shape(tempInt)) + tempFloat = config_phase(p)%getFloats('anisobrittle_criticaldisplacement',requiredSize=size(tempInt)) kinematics_cleavage_opening_critDisp(1:size(tempInt),instance) = tempFloat - tempFloat = config_phase(p)%getFloats('anisobrittle_criticalload',requiredShape=shape(tempInt)) + tempFloat = config_phase(p)%getFloats('anisobrittle_criticalload',requiredSize=size(tempInt)) kinematics_cleavage_opening_critLoad(1:size(tempInt),instance) = tempFloat kinematics_cleavage_opening_Ncleavage(1:lattice_maxNcleavageFamily,instance) = & @@ -138,9 +138,11 @@ end subroutine kinematics_cleavage_opening_init !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the velocity gradient !-------------------------------------------------------------------------------------------------- -subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, Tstar_v, ipc, ip, el) +subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, ip, el) use prec, only: & tol_math_check + use math, only: & + math_mul33xx33 use material, only: & material_phase, & material_homog, & @@ -148,7 +150,6 @@ subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, Tstar_v, damageMapping use lattice, only: & lattice_Scleavage, & - lattice_Scleavage_v, & lattice_maxNcleavageFamily, & lattice_NcleavageSystem @@ -157,8 +158,8 @@ subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, Tstar_v, ipc, & !< grain number ip, & !< integration point number el !< element number - real(pReal), intent(in), dimension(6) :: & - Tstar_v !< 2nd Piola-Kirchhoff stress + real(pReal), intent(in), dimension(3,3) :: & + S real(pReal), intent(out), dimension(3,3) :: & Ld !< damage velocity gradient real(pReal), intent(out), dimension(3,3,3,3) :: & @@ -181,9 +182,9 @@ subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, Tstar_v, do f = 1_pInt,lattice_maxNcleavageFamily index_myFamily = sum(lattice_NcleavageSystem(1:f-1_pInt,phase)) ! at which index starts my family do i = 1_pInt,kinematics_cleavage_opening_Ncleavage(f,instance) ! process each (active) cleavage system in family - traction_d = dot_product(Tstar_v,lattice_Scleavage_v(1:6,1,index_myFamily+i,phase)) - traction_t = dot_product(Tstar_v,lattice_Scleavage_v(1:6,2,index_myFamily+i,phase)) - traction_n = dot_product(Tstar_v,lattice_Scleavage_v(1:6,3,index_myFamily+i,phase)) + traction_d = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,1,index_myFamily+i,phase)) + traction_t = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,2,index_myFamily+i,phase)) + traction_n = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,3,index_myFamily+i,phase)) traction_crit = kinematics_cleavage_opening_critLoad(f,instance)* & damage(homog)%p(damageOffset)*damage(homog)%p(damageOffset) udotd = & diff --git a/src/kinematics_slipplane_opening.f90 b/src/kinematics_slipplane_opening.f90 index 573fe7d78..33714d573 100644 --- a/src/kinematics_slipplane_opening.f90 +++ b/src/kinematics_slipplane_opening.f90 @@ -113,10 +113,10 @@ subroutine kinematics_slipplane_opening_init() tempInt = config_phase(p)%getInts('ncleavage') kinematics_slipplane_opening_Nslip(1:size(tempInt),instance) = tempInt - tempFloat = config_phase(p)%getFloats('anisoductile_criticalplasticstrain',requiredShape=shape(tempInt)) + tempFloat = config_phase(p)%getFloats('anisoductile_criticalplasticstrain',requiredSize=size(tempInt)) kinematics_slipplane_opening_critPlasticStrain(1:size(tempInt),instance) = tempFloat - tempFloat = config_phase(p)%getFloats('anisoductile_criticalload',requiredShape=shape(tempInt)) + tempFloat = config_phase(p)%getFloats('anisoductile_criticalload',requiredSize=size(tempInt)) kinematics_slipplane_opening_critLoad(1:size(tempInt),instance) = tempFloat kinematics_slipplane_opening_Nslip(1:lattice_maxNslipFamily,instance) = & @@ -136,9 +136,11 @@ end subroutine kinematics_slipplane_opening_init !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the velocity gradient !-------------------------------------------------------------------------------------------------- -subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, Tstar_v, ipc, ip, el) +subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, ip, el) use prec, only: & tol_math_check + use math, only: & + math_mul33xx33 use lattice, only: & lattice_maxNslipFamily, & lattice_NslipSystem, & @@ -151,9 +153,6 @@ subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, Tstar_v, damage, & damageMapping use math, only: & - math_Plain3333to99, & - math_symmetric33, & - math_Mandel33to6, & math_tensorproduct33 implicit none @@ -161,16 +160,14 @@ subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, Tstar_v, ipc, & !< grain number ip, & !< integration point number el !< element number - real(pReal), intent(in), dimension(6) :: & - Tstar_v !< 2nd Piola-Kirchhoff stress + real(pReal), intent(in), dimension(3,3) :: & + S real(pReal), intent(out), dimension(3,3) :: & Ld !< damage velocity gradient real(pReal), intent(out), dimension(3,3,3,3) :: & dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor) real(pReal), dimension(3,3) :: & projection_d, projection_t, projection_n !< projection modes 3x3 tensor - real(pReal), dimension(6) :: & - projection_d_v, projection_t_v, projection_n_v !< projection modes 3x3 vector integer(pInt) :: & instance, phase, & homog, damageOffset, & @@ -196,13 +193,10 @@ subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, Tstar_v, projection_n = math_tensorproduct33(lattice_sn(1:3,index_myFamily+i,phase),& lattice_sn(1:3,index_myFamily+i,phase)) - projection_d_v(1:6) = math_Mandel33to6(math_symmetric33(projection_d(1:3,1:3))) - projection_t_v(1:6) = math_Mandel33to6(math_symmetric33(projection_t(1:3,1:3))) - projection_n_v(1:6) = math_Mandel33to6(math_symmetric33(projection_n(1:3,1:3))) - traction_d = dot_product(Tstar_v,projection_d_v(1:6)) - traction_t = dot_product(Tstar_v,projection_t_v(1:6)) - traction_n = dot_product(Tstar_v,projection_n_v(1:6)) + traction_d = math_mul33xx33(S,projection_d) + traction_t = math_mul33xx33(S,projection_t) + traction_n = math_mul33xx33(S,projection_n) traction_crit = kinematics_slipplane_opening_critLoad(f,instance)* & damage(homog)%p(damageOffset) ! degrading critical load carrying capacity by damage diff --git a/src/lattice.f90 b/src/lattice.f90 index 9be30a5d3..410c14628 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -31,8 +31,7 @@ module lattice lattice_Scleavage !< Schmid matrices for cleavage systems real(pReal), allocatable, dimension(:,:,:,:), protected, public :: & - lattice_Sslip_v, & !< Mandel notation of lattice_Sslip - lattice_Scleavage_v !< Mandel notation of lattice_Scleavege + lattice_Sslip_v !< Mandel notation of lattice_Sslip real(pReal), allocatable, dimension(:,:,:), protected, public :: & lattice_sn, & !< normal direction of slip system @@ -776,7 +775,6 @@ subroutine lattice_init allocate(lattice_interactionSlipSlip(lattice_maxNslip,lattice_maxNslip,Nphases),source=0_pInt) ! other:me allocate(lattice_Scleavage(3,3,3,lattice_maxNslip,Nphases),source=0.0_pReal) - allocate(lattice_Scleavage_v(6,3,lattice_maxNslip,Nphases),source=0.0_pReal) allocate(lattice_NcleavageSystem(lattice_maxNcleavageFamily,Nphases),source=0_pInt) allocate(CoverA(Nphases),source=0.0_pReal) @@ -1060,13 +1058,6 @@ subroutine lattice_initializeStructure(myPhase,CoverA) enddo enddo - do i = 1_pInt,myNcleavage ! store slip system vectors and Schmid matrix for my structure - do j = 1_pInt,3_pInt - lattice_Scleavage_v(1:6,j,i,myPhase) = & - math_sym33to6(math_symmetric33(lattice_Scleavage(1:3,1:3,j,i,myPhase))) - enddo - enddo - end subroutine lattice_initializeStructure diff --git a/src/source_damage_anisoBrittle.f90 b/src/source_damage_anisoBrittle.f90 index b8bd3246d..eabf43799 100644 --- a/src/source_damage_anisoBrittle.f90 +++ b/src/source_damage_anisoBrittle.f90 @@ -309,7 +309,9 @@ end subroutine source_damage_anisoBrittle_init !-------------------------------------------------------------------------------------------------- !> @brief calculates derived quantities from state !-------------------------------------------------------------------------------------------------- -subroutine source_damage_anisoBrittle_dotState(Tstar_v, ipc, ip, el) +subroutine source_damage_anisoBrittle_dotState(S, ipc, ip, el) + use math, only: & + math_mul33xx33 use material, only: & phaseAt, phasememberAt, & sourceState, & @@ -317,7 +319,7 @@ subroutine source_damage_anisoBrittle_dotState(Tstar_v, ipc, ip, el) damage, & damageMapping use lattice, only: & - lattice_Scleavage_v, & + lattice_Scleavage, & lattice_maxNcleavageFamily, & lattice_NcleavageSystem @@ -326,8 +328,8 @@ subroutine source_damage_anisoBrittle_dotState(Tstar_v, ipc, ip, el) ipc, & !< component-ID of integration point ip, & !< integration point el !< element - real(pReal), intent(in), dimension(6) :: & - Tstar_v !< 2nd Piola Kirchhoff stress tensor (Mandel) + real(pReal), intent(in), dimension(3,3) :: & + S integer(pInt) :: & phase, & constituent, & @@ -350,9 +352,9 @@ subroutine source_damage_anisoBrittle_dotState(Tstar_v, ipc, ip, el) do f = 1_pInt,lattice_maxNcleavageFamily index_myFamily = sum(lattice_NcleavageSystem(1:f-1_pInt,phase)) ! at which index starts my family do i = 1_pInt,source_damage_anisoBrittle_Ncleavage(f,instance) ! process each (active) cleavage system in family - traction_d = dot_product(Tstar_v,lattice_Scleavage_v(1:6,1,index_myFamily+i,phase)) - traction_t = dot_product(Tstar_v,lattice_Scleavage_v(1:6,2,index_myFamily+i,phase)) - traction_n = dot_product(Tstar_v,lattice_Scleavage_v(1:6,3,index_myFamily+i,phase)) + traction_d = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,1,index_myFamily+i,phase)) + traction_t = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,2,index_myFamily+i,phase)) + traction_n = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,3,index_myFamily+i,phase)) traction_crit = source_damage_anisoBrittle_critLoad(f,instance)* & damage(homog)%p(damageOffset)*damage(homog)%p(damageOffset)