From b5a10f23868f00a7fefc0914b4d97b9de49b5639 Mon Sep 17 00:00:00 2001 From: Sharan Roongta Date: Thu, 9 Jul 2020 01:01:08 +0200 Subject: [PATCH 01/21] sources and kinematics modules under submodules --- src/constitutive.f90 | 592 ++++++++++----------- src/constitutive_damage.f90 | 163 ++++++ src/constitutive_plastic.f90 | 274 ++++++++++ src/constitutive_plastic_disloUCLA.f90 | 2 +- src/constitutive_plastic_dislotwin.f90 | 2 +- src/constitutive_plastic_isotropic.f90 | 2 +- src/constitutive_plastic_kinehardening.f90 | 2 +- src/constitutive_plastic_none.f90 | 2 +- src/constitutive_plastic_nonlocal.f90 | 2 +- src/constitutive_plastic_phenopowerlaw.f90 | 2 +- src/constitutive_thermal.f90 | 112 ++++ src/damage_local.f90 | 40 +- src/damage_nonlocal.f90 | 41 +- src/kinematics_cleavage_opening.f90 | 20 +- src/kinematics_slipplane_opening.f90 | 20 +- src/kinematics_thermal_expansion.f90 | 57 +- src/source_damage_anisoBrittle.f90 | 29 +- src/source_damage_anisoDuctile.f90 | 33 +- src/source_damage_isoBrittle.f90 | 31 +- src/source_damage_isoDuctile.f90 | 31 +- src/source_thermal_dissipation.f90 | 19 +- src/source_thermal_externalheat.f90 | 21 +- src/thermal_adiabatic.f90 | 44 +- src/thermal_conduction.f90 | 43 +- 24 files changed, 937 insertions(+), 647 deletions(-) create mode 100644 src/constitutive_damage.f90 create mode 100644 src/constitutive_plastic.f90 create mode 100644 src/constitutive_thermal.f90 diff --git a/src/constitutive.f90 b/src/constitutive.f90 index d920c481f..5f509b0ba 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -16,15 +16,6 @@ module constitutive use lattice use discretization use geometry_plastic_nonlocal - use source_thermal_dissipation - use source_thermal_externalheat - use source_damage_isoBrittle - use source_damage_isoDuctile - use source_damage_anisoBrittle - use source_damage_anisoDuctile - use kinematics_cleavage_opening - use kinematics_slipplane_opening - use kinematics_thermal_expansion implicit none private @@ -35,114 +26,31 @@ module constitutive interface - module subroutine plastic_none_init - end subroutine plastic_none_init + module subroutine plastic_init + end subroutine plastic_init - module subroutine plastic_isotropic_init - end subroutine plastic_isotropic_init + module subroutine damage_init + end subroutine damage_init - module subroutine plastic_phenopowerlaw_init - end subroutine plastic_phenopowerlaw_init + module subroutine thermal_init + end subroutine thermal_init - module subroutine plastic_kinehardening_init - end subroutine plastic_kinehardening_init +! module pure function kinematics_thermal_expansion_initialStrain(homog,phase,offset) +! +! integer, intent(in) :: & +! phase, & +! homog, & +! offset +! end function kinematics_thermal_expansion_initialStrain - module subroutine plastic_dislotwin_init - end subroutine plastic_dislotwin_init - - module subroutine plastic_disloUCLA_init - end subroutine plastic_disloUCLA_init - - module subroutine plastic_nonlocal_init - end subroutine plastic_nonlocal_init - - - module subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) - real(pReal), dimension(3,3), intent(out) :: & - Lp !< plastic velocity gradient - real(pReal), dimension(3,3,3,3), intent(out) :: & - dLp_dMp !< derivative of Lp with respect to the Mandel stress - - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - integer, intent(in) :: & - instance, & - of - end subroutine plastic_isotropic_LpAndItsTangent - - pure module subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) - real(pReal), dimension(3,3), intent(out) :: & - Lp !< plastic velocity gradient - real(pReal), dimension(3,3,3,3), intent(out) :: & - dLp_dMp !< derivative of Lp with respect to the Mandel stress - - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - integer, intent(in) :: & - instance, & - of - end subroutine plastic_phenopowerlaw_LpAndItsTangent - - pure module subroutine plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) - real(pReal), dimension(3,3), intent(out) :: & - Lp !< plastic velocity gradient - real(pReal), dimension(3,3,3,3), intent(out) :: & - dLp_dMp !< derivative of Lp with respect to the Mandel stress - - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - integer, intent(in) :: & - instance, & - of - end subroutine plastic_kinehardening_LpAndItsTangent - - module subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of) - real(pReal), dimension(3,3), intent(out) :: & - Lp !< plastic velocity gradient - real(pReal), dimension(3,3,3,3), intent(out) :: & - dLp_dMp !< derivative of Lp with respect to the Mandel stress - - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - real(pReal), intent(in) :: & - T - integer, intent(in) :: & - instance, & - of - end subroutine plastic_dislotwin_LpAndItsTangent - - pure module subroutine plastic_disloUCLA_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of) - real(pReal), dimension(3,3), intent(out) :: & - Lp !< plastic velocity gradient - real(pReal), dimension(3,3,3,3), intent(out) :: & - dLp_dMp !< derivative of Lp with respect to the Mandel stress - - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - real(pReal), intent(in) :: & - T - integer, intent(in) :: & - instance, & - of - end subroutine plastic_disloUCLA_LpAndItsTangent - - module subroutine plastic_nonlocal_LpAndItsTangent(Lp,dLp_dMp, & - Mp,Temperature,instance,of,ip,el) - real(pReal), dimension(3,3), intent(out) :: & - Lp !< plastic velocity gradient - real(pReal), dimension(3,3,3,3), intent(out) :: & - dLp_dMp !< derivative of Lp with respect to the Mandel stress - - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - real(pReal), intent(in) :: & - Temperature - integer, intent(in) :: & - instance, & - of, & - ip, & !< current integration point - el !< current element number - end subroutine plastic_nonlocal_LpAndItsTangent + module function plastic_dislotwin_homogenizedC(ipc,ip,el) result(homogenizedC) + real(pReal), dimension(6,6) :: & + homogenizedC + integer, intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + end function plastic_dislotwin_homogenizedC module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,of) @@ -158,6 +66,45 @@ module constitutive of end subroutine plastic_isotropic_LiAndItsTangent + module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, ip, el) + + integer, intent(in) :: & + ipc, & !< grain number + ip, & !< integration point number + el !< element number + 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) + end subroutine kinematics_cleavage_opening_LiAndItsTangent + + module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, ip, el) + + integer, intent(in) :: & + ipc, & !< grain number + ip, & !< integration point number + el !< element number + 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) + end subroutine kinematics_slipplane_opening_LiAndItsTangent + + module subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, ipc, ip, el) + + integer, intent(in) :: & + ipc, & !< grain number + ip, & !< integration point number + el !< element number + real(pReal), intent(out), dimension(3,3) :: & + Li !< thermal velocity gradient + real(pReal), intent(out), dimension(3,3,3,3) :: & + dLi_dTstar !< derivative of Li with respect to Tstar (4th-order tensor defined to be zero) + end subroutine kinematics_thermal_expansion_LiAndItsTangent module subroutine plastic_isotropic_dotState(Mp,instance,of) real(pReal), dimension(3,3), intent(in) :: & @@ -220,31 +167,39 @@ module constitutive el !< current element number end subroutine plastic_nonlocal_dotState + + module subroutine source_damage_anisoBrittle_dotState(S, ipc, ip, el) - module subroutine plastic_dislotwin_dependentState(T,instance,of) - integer, intent(in) :: & - instance, & - of - real(pReal), intent(in) :: & - T - end subroutine plastic_dislotwin_dependentState - - module subroutine plastic_disloUCLA_dependentState(instance,of) - integer, intent(in) :: & - instance, & - of - end subroutine plastic_disloUCLA_dependentState - - module subroutine plastic_nonlocal_dependentState(F, Fp, instance, of, ip, el) - real(pReal), dimension(3,3), intent(in) :: & - F, & - Fp integer, intent(in) :: & - instance, & - of, & - ip, & - el - end subroutine plastic_nonlocal_dependentState + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + real(pReal), intent(in), dimension(3,3) :: & + S + end subroutine source_damage_anisoBrittle_dotState + + module subroutine source_damage_anisoDuctile_dotState(ipc, ip, el) + + integer, intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + end subroutine source_damage_anisoDuctile_dotState + + module subroutine source_damage_isoDuctile_dotState(ipc, ip, el) + + integer, intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + end subroutine source_damage_isoDuctile_dotState + + module subroutine source_thermal_externalheat_dotState(phase, of) + + integer, intent(in) :: & + phase, & + of + end subroutine source_thermal_externalheat_dotState module subroutine plastic_kinehardening_deltaState(Mp,instance,of) @@ -265,24 +220,17 @@ module constitutive el end subroutine plastic_nonlocal_deltaState - - module function plastic_dislotwin_homogenizedC(ipc,ip,el) result(homogenizedC) - real(pReal), dimension(6,6) :: & - homogenizedC - integer, intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element - end function plastic_dislotwin_homogenizedC - - module subroutine plastic_nonlocal_updateCompatibility(orientation,instance,i,e) + module subroutine source_damage_isoBrittle_deltaState(C, Fe, ipc, ip, el) integer, intent(in) :: & - instance, & - i, & - e - type(rotation), dimension(1,discretization_nIP,discretization_nElem), intent(in) :: & - orientation !< crystal orientation - end subroutine plastic_nonlocal_updateCompatibility + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + real(pReal), intent(in), dimension(3,3) :: & + Fe + real(pReal), intent(in), dimension(6,6) :: & + C + end subroutine source_damage_isoBrittle_deltaState + module subroutine plastic_isotropic_results(instance,group) @@ -315,6 +263,88 @@ module constitutive character(len=*), intent(in) :: group end subroutine plastic_nonlocal_results + module subroutine source_damage_anisoBrittle_results(phase,group) + integer, intent(in) :: phase + character(len=*), intent(in) :: group + end subroutine source_damage_anisoBrittle_results + + module subroutine source_damage_anisoDuctile_results(phase,group) + integer, intent(in) :: phase + character(len=*), intent(in) :: group + end subroutine source_damage_anisoDuctile_results + + module subroutine source_damage_isoBrittle_results(phase,group) + integer, intent(in) :: phase + character(len=*), intent(in) :: group + end subroutine source_damage_isoBrittle_results + + module subroutine source_damage_isoDuctile_results(phase,group) + integer, intent(in) :: phase + character(len=*), intent(in) :: group + end subroutine source_damage_isoDuctile_results + + module subroutine plastic_dependentState(F, Fp, ipc, ip, el) + + integer, intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + real(pReal), intent(in), dimension(3,3) :: & + F, & !< elastic deformation gradient + Fp !< plastic deformation gradient + end subroutine plastic_dependentState + + module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & + S, Fi, ipc, ip, el) + integer, intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + real(pReal), intent(in), dimension(3,3) :: & + S, & !< 2nd Piola-Kirchhoff stress + Fi !< intermediate deformation gradient + real(pReal), intent(out), dimension(3,3) :: & + Lp !< plastic velocity gradient + real(pReal), intent(out), dimension(3,3,3,3) :: & + dLp_dS, & + dLp_dFi !< derivative of Lp with respect to Fi + + end subroutine plastic_LpAndItsTangents + + module subroutine plastic_nonlocal_updateCompatibility(orientation,instance,i,e) + integer, intent(in) :: & + instance, & + i, & + e + type(rotation), dimension(1,discretization_nIP,discretization_nElem), intent(in) :: & + orientation !< crystal orientation + end subroutine plastic_nonlocal_updateCompatibility + + module subroutine damage_source_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el) + integer, intent(in) :: & + ip, & !< integration point number + el !< element number + real(pReal), intent(in) :: & + phi + real(pReal), intent(inout) :: & + phiDot, & + dPhiDot_dPhi + end subroutine damage_source_getRateAndItsTangents + + module subroutine thermal_source_getRateAndItsTangents(TDot, dTDot_dT, T, Tstar, Lp, ip, el) + integer, intent(in) :: & + ip, & !< integration point number + el !< element number + real(pReal), intent(in) :: & + T + real(pReal), intent(in), dimension(:,:,:,:,:) :: & + Tstar, & + Lp + real(pReal), intent(inout) :: & + TDot, & + dTDot_dT + end subroutine thermal_source_getRateAndItsTangents + end interface @@ -332,16 +362,18 @@ module constitutive type(tDebugOptions) :: debugConstitutive public :: & - plastic_nonlocal_updateCompatibility, & constitutive_init, & constitutive_homogenizedC, & - constitutive_dependentState, & constitutive_LpAndItsTangents, & + constitutive_dependentState, & constitutive_LiAndItsTangents, & constitutive_initialFi, & constitutive_SandItsTangents, & constitutive_collectDotState, & constitutive_deltaState, & + plastic_nonlocal_updateCompatibility, & + constitutive_damage_getRateAndItsTangents, & + constitutive_thermal_getRateAndItsTangents, & constitutive_results contains @@ -368,31 +400,9 @@ subroutine constitutive_init !-------------------------------------------------------------------------------------------------- ! initialized plasticity - if (any(phase_plasticity == PLASTICITY_NONE_ID)) call plastic_none_init - if (any(phase_plasticity == PLASTICITY_ISOTROPIC_ID)) call plastic_isotropic_init - if (any(phase_plasticity == PLASTICITY_PHENOPOWERLAW_ID)) call plastic_phenopowerlaw_init - if (any(phase_plasticity == PLASTICITY_KINEHARDENING_ID)) call plastic_kinehardening_init - if (any(phase_plasticity == PLASTICITY_DISLOTWIN_ID)) call plastic_dislotwin_init - if (any(phase_plasticity == PLASTICITY_DISLOUCLA_ID)) call plastic_disloucla_init - if (any(phase_plasticity == PLASTICITY_NONLOCAL_ID)) then - call plastic_nonlocal_init - else - call geometry_plastic_nonlocal_disable - endif -!-------------------------------------------------------------------------------------------------- -! initialize source mechanisms - if (any(phase_source == SOURCE_thermal_dissipation_ID)) call source_thermal_dissipation_init - if (any(phase_source == SOURCE_thermal_externalheat_ID)) call source_thermal_externalheat_init - if (any(phase_source == SOURCE_damage_isoBrittle_ID)) call source_damage_isoBrittle_init - if (any(phase_source == SOURCE_damage_isoDuctile_ID)) call source_damage_isoDuctile_init - if (any(phase_source == SOURCE_damage_anisoBrittle_ID)) call source_damage_anisoBrittle_init - if (any(phase_source == SOURCE_damage_anisoDuctile_ID)) call source_damage_anisoDuctile_init - -!-------------------------------------------------------------------------------------------------- -! initialize kinematic mechanisms - if (any(phase_kinematics == KINEMATICS_cleavage_opening_ID)) call kinematics_cleavage_opening_init - if (any(phase_kinematics == KINEMATICS_slipplane_opening_ID)) call kinematics_slipplane_opening_init - if (any(phase_kinematics == KINEMATICS_thermal_expansion_ID)) call kinematics_thermal_expansion_init + call plastic_init + call damage_init + call thermal_init write(6,'(/,a)') ' <<<+- constitutive init -+>>>'; flush(6) @@ -416,6 +426,43 @@ subroutine constitutive_init end subroutine constitutive_init +subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & + S, Fi, ipc, ip, el) + integer, intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + real(pReal), intent(in), dimension(3,3) :: & + S, & !< 2nd Piola-Kirchhoff stress + Fi !< intermediate deformation gradient + real(pReal), intent(out), dimension(3,3) :: & + Lp !< plastic velocity gradient + real(pReal), intent(out), dimension(3,3,3,3) :: & + dLp_dS, & + dLp_dFi !< derivative of Lp with respect to Fi + + call plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, S, Fi, ipc, ip, el) + +end subroutine constitutive_LpAndItsTangents + +!-------------------------------------------------------------------------------------------------- +!> @brief calls microstructure function of the different constitutive models +!-------------------------------------------------------------------------------------------------- +subroutine constitutive_dependentState(F, Fp, ipc, ip, el) + + integer, intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + real(pReal), intent(in), dimension(3,3) :: & + F, & !< elastic deformation gradient + Fp !< plastic deformation gradient + + call plastic_dependentState(F,Fp,ipc,ip,el) + +end subroutine constitutive_dependentState + + !-------------------------------------------------------------------------------------------------- !> @brief returns the homogenize elasticity matrix !> ToDo: homogenizedC66 would be more consistent @@ -438,111 +485,6 @@ function constitutive_homogenizedC(ipc,ip,el) end function constitutive_homogenizedC -!-------------------------------------------------------------------------------------------------- -!> @brief calls microstructure function of the different constitutive models -!-------------------------------------------------------------------------------------------------- -subroutine constitutive_dependentState(F, Fp, ipc, ip, el) - - integer, intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element - real(pReal), intent(in), dimension(3,3) :: & - F, & !< elastic deformation gradient - Fp !< plastic deformation gradient - integer :: & - ho, & !< homogenization - tme, & !< thermal member position - instance, of - - ho = material_homogenizationAt(el) - tme = thermalMapping(ho)%p(ip,el) - of = material_phasememberAt(ipc,ip,el) - instance = phase_plasticityInstance(material_phaseAt(ipc,el)) - - plasticityType: select case (phase_plasticity(material_phaseAt(ipc,el))) - case (PLASTICITY_DISLOTWIN_ID) plasticityType - call plastic_dislotwin_dependentState(temperature(ho)%p(tme),instance,of) - case (PLASTICITY_DISLOUCLA_ID) plasticityType - call plastic_disloUCLA_dependentState(instance,of) - case (PLASTICITY_NONLOCAL_ID) plasticityType - call plastic_nonlocal_dependentState (F,Fp,instance,of,ip,el) - end select plasticityType - -end subroutine constitutive_dependentState - - -!-------------------------------------------------------------------------------------------------- -!> @brief contains the constitutive equation for calculating the velocity gradient -! ToDo: Discuss whether it makes sense if crystallite handles the configuration conversion, i.e. -! Mp in, dLp_dMp out -!-------------------------------------------------------------------------------------------------- -subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & - S, Fi, ipc, ip, el) - integer, intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element - real(pReal), intent(in), dimension(3,3) :: & - S, & !< 2nd Piola-Kirchhoff stress - Fi !< intermediate deformation gradient - real(pReal), intent(out), dimension(3,3) :: & - Lp !< plastic velocity gradient - real(pReal), intent(out), dimension(3,3,3,3) :: & - dLp_dS, & - dLp_dFi !< derivative of Lp with respect to Fi - real(pReal), dimension(3,3,3,3) :: & - dLp_dMp !< derivative of Lp with respect to Mandel stress - real(pReal), dimension(3,3) :: & - Mp !< Mandel stress work conjugate with Lp - integer :: & - ho, & !< homogenization - tme !< thermal member position - integer :: & - i, j, instance, of - - ho = material_homogenizationAt(el) - tme = thermalMapping(ho)%p(ip,el) - - Mp = matmul(matmul(transpose(Fi),Fi),S) - of = material_phasememberAt(ipc,ip,el) - instance = phase_plasticityInstance(material_phaseAt(ipc,el)) - - plasticityType: select case (phase_plasticity(material_phaseAt(ipc,el))) - - case (PLASTICITY_NONE_ID) plasticityType - Lp = 0.0_pReal - dLp_dMp = 0.0_pReal - - case (PLASTICITY_ISOTROPIC_ID) plasticityType - call plastic_isotropic_LpAndItsTangent (Lp,dLp_dMp,Mp,instance,of) - - case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType - call plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) - - case (PLASTICITY_KINEHARDENING_ID) plasticityType - call plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) - - case (PLASTICITY_NONLOCAL_ID) plasticityType - call plastic_nonlocal_LpAndItsTangent (Lp,dLp_dMp,Mp, temperature(ho)%p(tme),instance,of,ip,el) - - case (PLASTICITY_DISLOTWIN_ID) plasticityType - call plastic_dislotwin_LpAndItsTangent (Lp,dLp_dMp,Mp,temperature(ho)%p(tme),instance,of) - - case (PLASTICITY_DISLOUCLA_ID) plasticityType - call plastic_disloucla_LpAndItsTangent (Lp,dLp_dMp,Mp,temperature(ho)%p(tme),instance,of) - - end select plasticityType - - do i=1,3; do j=1,3 - dLp_dFi(i,j,1:3,1:3) = matmul(matmul(Fi,S),transpose(dLp_dMp(i,j,1:3,1:3))) + & - matmul(matmul(Fi,dLp_dMp(i,j,1:3,1:3)),S) - dLp_dS(i,j,1:3,1:3) = matmul(matmul(transpose(Fi),Fi),dLp_dMp(i,j,1:3,1:3)) ! ToDo: @PS: why not: dLp_dMp:(FiT Fi) - enddo; enddo - -end subroutine constitutive_LpAndItsTangents - - !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the velocity gradient ! ToDo: MD: S is Mi? @@ -649,7 +591,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 @@ -890,43 +832,97 @@ function constitutive_deltaState(S, Fe, Fi, ipc, ip, el, phase, of) result(broke end function constitutive_deltaState +subroutine constitutive_thermal_getRateAndItsTangents(Tdot, dTDot_dT, T, Tstar, Lp, ip, el) + + integer, intent(in) :: & + ip, & !< integration point number + el !< element number + real(pReal), intent(in) :: & + T + real(pReal), intent(in), dimension(:,:,:,:,:) :: & + Tstar, & + Lp + real(pReal), intent(inout) :: & + Tdot, & + dTdot_dT + + call thermal_source_getRateAndItsTangents(Tdot, dTdot_dT, T, Tstar, Lp, ip, el) + +end subroutine constitutive_thermal_getRateAndItsTangents + +subroutine constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el) + + integer, intent(in) :: & + ip, & !< integration point number + el !< element number + real(pReal), intent(in) :: & + phi + real(pReal), intent(inout) :: & + phiDot, & + dPhiDot_dPhi + + call damage_source_getRateAndItsTangents(phiDot,dPhiDot_dPhi,phi,ip,el) + +end subroutine constitutive_damage_getRateAndItsTangents + !-------------------------------------------------------------------------------------------------- !> @brief writes constitutive results to HDF5 output file !-------------------------------------------------------------------------------------------------- subroutine constitutive_results - integer :: p - character(len=pStringLen) :: group - do p=1,size(config_name_phase) + integer :: p,i + character(len=pStringLen) :: group,group_plastic,group_sources + + plasticityLoop: do p=1,size(config_name_phase) group = trim('current/constituent')//'/'//trim(config_name_phase(p)) call results_closeGroup(results_addGroup(group)) - group = trim(group)//'/plastic' + group_plastic = trim(group)//'/plastic' - call results_closeGroup(results_addGroup(group)) + call results_closeGroup(results_addGroup(group_plastic)) select case(phase_plasticity(p)) case(PLASTICITY_ISOTROPIC_ID) - call plastic_isotropic_results(phase_plasticityInstance(p),group) + call plastic_isotropic_results(phase_plasticityInstance(p),group_plastic) case(PLASTICITY_PHENOPOWERLAW_ID) - call plastic_phenopowerlaw_results(phase_plasticityInstance(p),group) + call plastic_phenopowerlaw_results(phase_plasticityInstance(p),group_plastic) case(PLASTICITY_KINEHARDENING_ID) - call plastic_kinehardening_results(phase_plasticityInstance(p),group) + call plastic_kinehardening_results(phase_plasticityInstance(p),group_plastic) case(PLASTICITY_DISLOTWIN_ID) - call plastic_dislotwin_results(phase_plasticityInstance(p),group) + call plastic_dislotwin_results(phase_plasticityInstance(p),group_plastic) case(PLASTICITY_DISLOUCLA_ID) - call plastic_disloUCLA_results(phase_plasticityInstance(p),group) + call plastic_disloUCLA_results(phase_plasticityInstance(p),group_plastic) case(PLASTICITY_NONLOCAL_ID) - call plastic_nonlocal_results(phase_plasticityInstance(p),group) + call plastic_nonlocal_results(phase_plasticityInstance(p),group_plastic) end select - enddo + sourceLoop: do i = 1, phase_Nsources(p) + group_sources = trim(group)//'/sources' + + call results_closeGroup(results_addGroup(group_sources)) + + sourceType: select case (phase_source(i,p)) + + case (SOURCE_damage_anisoBrittle_ID) sourceType + call source_damage_anisoBrittle_results(p,group_sources) + case (SOURCE_damage_anisoDuctile_ID) sourceType + call source_damage_anisoDuctile_results(p,group_sources) + case (SOURCE_damage_isoBrittle_ID) sourceType + call source_damage_isoBrittle_results(p,group_sources) + case (SOURCE_damage_isoDuctile_ID) sourceType + call source_damage_isoDuctile_results(p,group_sources) + end select sourceType + + enddo SourceLoop + + enddo plasticityLoop + end subroutine constitutive_results diff --git a/src/constitutive_damage.f90 b/src/constitutive_damage.f90 new file mode 100644 index 000000000..cb81c252d --- /dev/null +++ b/src/constitutive_damage.f90 @@ -0,0 +1,163 @@ +submodule(constitutive) constitutive_damage + + implicit none + + interface + + module subroutine source_damage_anisoBrittle_init + end subroutine source_damage_anisoBrittle_init + + module subroutine source_damage_anisoDuctile_init + end subroutine source_damage_anisoDuctile_init + + module subroutine source_damage_isoBrittle_init + end subroutine source_damage_isoBrittle_init + + module subroutine source_damage_isoDuctile_init + end subroutine source_damage_isoDuctile_init + + module subroutine kinematics_cleavage_opening_init + end subroutine kinematics_cleavage_opening_init + + module subroutine kinematics_slipplane_opening_init + end subroutine kinematics_slipplane_opening_init + + module subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) + + integer, intent(in) :: & + phase, & + constituent + real(pReal), intent(in) :: & + phi + real(pReal), intent(out) :: & + localphiDot, & + dLocalphiDot_dPhi + + end subroutine source_damage_anisoBrittle_getRateAndItsTangent + + module subroutine source_damage_anisoDuctile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) + + integer, intent(in) :: & + phase, & + constituent + real(pReal), intent(in) :: & + phi + real(pReal), intent(out) :: & + localphiDot, & + dLocalphiDot_dPhi + + end subroutine source_damage_anisoDuctile_getRateAndItsTangent + + module subroutine source_damage_isoBrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) + + integer, intent(in) :: & + phase, & + constituent + real(pReal), intent(in) :: & + phi + real(pReal), intent(out) :: & + localphiDot, & + dLocalphiDot_dPhi + + end subroutine source_damage_isoBrittle_getRateAndItsTangent + + module subroutine source_damage_isoDuctile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) + + integer, intent(in) :: & + phase, & + constituent + real(pReal), intent(in) :: & + phi + real(pReal), intent(out) :: & + localphiDot, & + dLocalphiDot_dPhi + + end subroutine source_damage_isoDuctile_getRateAndItsTangent + + module subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDot_dT, Tstar, Lp, phase) + + integer, intent(in) :: & + phase + real(pReal), intent(in), dimension(3,3) :: & + Tstar + real(pReal), intent(in), dimension(3,3) :: & + Lp + + real(pReal), intent(out) :: & + TDot, & + dTDot_dT + + end subroutine source_thermal_dissipation_getRateAndItsTangent + end interface + +contains + +module subroutine damage_init + +! initialize source mechanisms + if (any(phase_source == SOURCE_damage_isoBrittle_ID)) call source_damage_isoBrittle_init + if (any(phase_source == SOURCE_damage_isoDuctile_ID)) call source_damage_isoDuctile_init + if (any(phase_source == SOURCE_damage_anisoBrittle_ID)) call source_damage_anisoBrittle_init + if (any(phase_source == SOURCE_damage_anisoDuctile_ID)) call source_damage_anisoDuctile_init + +!-------------------------------------------------------------------------------------------------- +! initialize kinematic mechanisms + if (any(phase_kinematics == KINEMATICS_cleavage_opening_ID)) call kinematics_cleavage_opening_init + if (any(phase_kinematics == KINEMATICS_slipplane_opening_ID)) call kinematics_slipplane_opening_init + +end subroutine damage_init + + +module subroutine damage_source_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el) + + integer, intent(in) :: & + ip, & !< integration point number + el !< element number + real(pReal), intent(in) :: & + phi + real(pReal), intent(inout) :: & + phiDot, & + dPhiDot_dPhi + + real(pReal) :: & + localphiDot, & + dLocalphiDot_dPhi + integer :: & + phase, & + grain, & + source, & + constituent + + phiDot = 0.0_pReal + dPhiDot_dPhi = 0.0_pReal + + do grain = 1, homogenization_Ngrains(material_homogenizationAt(el)) + phase = material_phaseAt(grain,el) + constituent = material_phasememberAt(grain,ip,el) + do source = 1, phase_Nsources(phase) + select case(phase_source(source,phase)) + case (SOURCE_damage_isoBrittle_ID) + call source_damage_isobrittle_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) + + case (SOURCE_damage_isoDuctile_ID) + call source_damage_isoductile_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) + + case (SOURCE_damage_anisoBrittle_ID) + call source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) + + case (SOURCE_damage_anisoDuctile_ID) + call source_damage_anisoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) + + case default + localphiDot = 0.0_pReal + dLocalphiDot_dPhi = 0.0_pReal + + end select + phiDot = phiDot + localphiDot + dPhiDot_dPhi = dPhiDot_dPhi + dLocalphiDot_dPhi + enddo + enddo + +end subroutine damage_source_getRateAndItsTangents + +end submodule diff --git a/src/constitutive_plastic.f90 b/src/constitutive_plastic.f90 new file mode 100644 index 000000000..d4d0f19a2 --- /dev/null +++ b/src/constitutive_plastic.f90 @@ -0,0 +1,274 @@ +submodule(constitutive) constitutive_plastic + + implicit none + + interface + + module subroutine plastic_none_init + end subroutine plastic_none_init + + module subroutine plastic_isotropic_init + end subroutine plastic_isotropic_init + + module subroutine plastic_phenopowerlaw_init + end subroutine plastic_phenopowerlaw_init + + module subroutine plastic_kinehardening_init + end subroutine plastic_kinehardening_init + + module subroutine plastic_dislotwin_init + end subroutine plastic_dislotwin_init + + module subroutine plastic_disloUCLA_init + end subroutine plastic_disloUCLA_init + + module subroutine plastic_nonlocal_init + end subroutine plastic_nonlocal_init + + + module subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) + real(pReal), dimension(3,3), intent(out) :: & + Lp !< plastic velocity gradient + real(pReal), dimension(3,3,3,3), intent(out) :: & + dLp_dMp !< derivative of Lp with respect to the Mandel stress + + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + integer, intent(in) :: & + instance, & + of + end subroutine plastic_isotropic_LpAndItsTangent + + pure module subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) + real(pReal), dimension(3,3), intent(out) :: & + Lp !< plastic velocity gradient + real(pReal), dimension(3,3,3,3), intent(out) :: & + dLp_dMp !< derivative of Lp with respect to the Mandel stress + + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + integer, intent(in) :: & + instance, & + of + end subroutine plastic_phenopowerlaw_LpAndItsTangent + + pure module subroutine plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) + real(pReal), dimension(3,3), intent(out) :: & + Lp !< plastic velocity gradient + real(pReal), dimension(3,3,3,3), intent(out) :: & + dLp_dMp !< derivative of Lp with respect to the Mandel stress + + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + integer, intent(in) :: & + instance, & + of + end subroutine plastic_kinehardening_LpAndItsTangent + + module subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of) + real(pReal), dimension(3,3), intent(out) :: & + Lp !< plastic velocity gradient + real(pReal), dimension(3,3,3,3), intent(out) :: & + dLp_dMp !< derivative of Lp with respect to the Mandel stress + + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + real(pReal), intent(in) :: & + T + integer, intent(in) :: & + instance, & + of + end subroutine plastic_dislotwin_LpAndItsTangent + + pure module subroutine plastic_disloUCLA_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of) + real(pReal), dimension(3,3), intent(out) :: & + Lp !< plastic velocity gradient + real(pReal), dimension(3,3,3,3), intent(out) :: & + dLp_dMp !< derivative of Lp with respect to the Mandel stress + + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + real(pReal), intent(in) :: & + T + integer, intent(in) :: & + instance, & + of + end subroutine plastic_disloUCLA_LpAndItsTangent + + module subroutine plastic_nonlocal_LpAndItsTangent(Lp,dLp_dMp, & + Mp,Temperature,instance,of,ip,el) + real(pReal), dimension(3,3), intent(out) :: & + Lp !< plastic velocity gradient + real(pReal), dimension(3,3,3,3), intent(out) :: & + dLp_dMp !< derivative of Lp with respect to the Mandel stress + + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + real(pReal), intent(in) :: & + Temperature + integer, intent(in) :: & + instance, & + of, & + ip, & !< current integration point + el !< current element number + end subroutine plastic_nonlocal_LpAndItsTangent + + + module subroutine plastic_dislotwin_dependentState(T,instance,of) + integer, intent(in) :: & + instance, & + of + real(pReal), intent(in) :: & + T + end subroutine plastic_dislotwin_dependentState + + module subroutine plastic_disloUCLA_dependentState(instance,of) + integer, intent(in) :: & + instance, & + of + end subroutine plastic_disloUCLA_dependentState + + module subroutine plastic_nonlocal_dependentState(F, Fp, instance, of, ip, el) + real(pReal), dimension(3,3), intent(in) :: & + F, & + Fp + integer, intent(in) :: & + instance, & + of, & + ip, & + el + end subroutine plastic_nonlocal_dependentState + + end interface + + +contains + + +!-------------------------------------------------------------------------------------------------- +!> @brief allocates arrays pointing to array of the various constitutive modules +!-------------------------------------------------------------------------------------------------- +module subroutine plastic_init + +!-------------------------------------------------------------------------------------------------- +! initialized plasticity + if (any(phase_plasticity == PLASTICITY_NONE_ID)) call plastic_none_init + if (any(phase_plasticity == PLASTICITY_ISOTROPIC_ID)) call plastic_isotropic_init + if (any(phase_plasticity == PLASTICITY_PHENOPOWERLAW_ID)) call plastic_phenopowerlaw_init + if (any(phase_plasticity == PLASTICITY_KINEHARDENING_ID)) call plastic_kinehardening_init + if (any(phase_plasticity == PLASTICITY_DISLOTWIN_ID)) call plastic_dislotwin_init + if (any(phase_plasticity == PLASTICITY_DISLOUCLA_ID)) call plastic_disloucla_init + if (any(phase_plasticity == PLASTICITY_NONLOCAL_ID)) then + call plastic_nonlocal_init + else + call geometry_plastic_nonlocal_disable + endif + +end subroutine plastic_init + + +!-------------------------------------------------------------------------------------------------- +!> @brief calls microstructure function of the different constitutive models +!-------------------------------------------------------------------------------------------------- +module subroutine plastic_dependentState(F, Fp, ipc, ip, el) + + integer, intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + real(pReal), intent(in), dimension(3,3) :: & + F, & !< elastic deformation gradient + Fp !< plastic deformation gradient + integer :: & + ho, & !< homogenization + tme, & !< thermal member position + instance, of + + ho = material_homogenizationAt(el) + tme = thermalMapping(ho)%p(ip,el) + of = material_phasememberAt(ipc,ip,el) + instance = phase_plasticityInstance(material_phaseAt(ipc,el)) + + plasticityType: select case (phase_plasticity(material_phaseAt(ipc,el))) + case (PLASTICITY_DISLOTWIN_ID) plasticityType + call plastic_dislotwin_dependentState(temperature(ho)%p(tme),instance,of) + case (PLASTICITY_DISLOUCLA_ID) plasticityType + call plastic_disloUCLA_dependentState(instance,of) + case (PLASTICITY_NONLOCAL_ID) plasticityType + call plastic_nonlocal_dependentState (F,Fp,instance,of,ip,el) + end select plasticityType + +end subroutine plastic_dependentState + +!-------------------------------------------------------------------------------------------------- +!> @brief contains the constitutive equation for calculating the velocity gradient +! ToDo: Discuss whether it makes sense if crystallite handles the configuration conversion, i.e. +! Mp in, dLp_dMp out +!-------------------------------------------------------------------------------------------------- +module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & + S, Fi, ipc, ip, el) + integer, intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + real(pReal), intent(in), dimension(3,3) :: & + S, & !< 2nd Piola-Kirchhoff stress + Fi !< intermediate deformation gradient + real(pReal), intent(out), dimension(3,3) :: & + Lp !< plastic velocity gradient + real(pReal), intent(out), dimension(3,3,3,3) :: & + dLp_dS, & + dLp_dFi !< derivative of Lp with respect to Fi + real(pReal), dimension(3,3,3,3) :: & + dLp_dMp !< derivative of Lp with respect to Mandel stress + real(pReal), dimension(3,3) :: & + Mp !< Mandel stress work conjugate with Lp + integer :: & + ho, & !< homogenization + tme !< thermal member position + integer :: & + i, j, instance, of + + ho = material_homogenizationAt(el) + tme = thermalMapping(ho)%p(ip,el) + + Mp = matmul(matmul(transpose(Fi),Fi),S) + of = material_phasememberAt(ipc,ip,el) + instance = phase_plasticityInstance(material_phaseAt(ipc,el)) + + plasticityType: select case (phase_plasticity(material_phaseAt(ipc,el))) + + case (PLASTICITY_NONE_ID) plasticityType + Lp = 0.0_pReal + dLp_dMp = 0.0_pReal + + case (PLASTICITY_ISOTROPIC_ID) plasticityType + call plastic_isotropic_LpAndItsTangent (Lp,dLp_dMp,Mp,instance,of) + + case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType + call plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) + + case (PLASTICITY_KINEHARDENING_ID) plasticityType + call plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) + + case (PLASTICITY_NONLOCAL_ID) plasticityType + call plastic_nonlocal_LpAndItsTangent (Lp,dLp_dMp,Mp, temperature(ho)%p(tme),instance,of,ip,el) + + case (PLASTICITY_DISLOTWIN_ID) plasticityType + call plastic_dislotwin_LpAndItsTangent (Lp,dLp_dMp,Mp,temperature(ho)%p(tme),instance,of) + + case (PLASTICITY_DISLOUCLA_ID) plasticityType + call plastic_disloucla_LpAndItsTangent (Lp,dLp_dMp,Mp,temperature(ho)%p(tme),instance,of) + + end select plasticityType + + do i=1,3; do j=1,3 + dLp_dFi(i,j,1:3,1:3) = matmul(matmul(Fi,S),transpose(dLp_dMp(i,j,1:3,1:3))) + & + matmul(matmul(Fi,dLp_dMp(i,j,1:3,1:3)),S) + dLp_dS(i,j,1:3,1:3) = matmul(matmul(transpose(Fi),Fi),dLp_dMp(i,j,1:3,1:3)) ! ToDo: @PS: why not: dLp_dMp:(FiT Fi) + enddo; enddo + +end subroutine plastic_LpAndItsTangents + +end submodule constitutive_plastic + diff --git a/src/constitutive_plastic_disloUCLA.f90 b/src/constitutive_plastic_disloUCLA.f90 index 8c94817b9..71f8f4a27 100644 --- a/src/constitutive_plastic_disloUCLA.f90 +++ b/src/constitutive_plastic_disloUCLA.f90 @@ -5,7 +5,7 @@ !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @brief crystal plasticity model for bcc metals, especially Tungsten !-------------------------------------------------------------------------------------------------- -submodule(constitutive) plastic_disloUCLA +submodule(constitutive:constitutive_plastic) plastic_disloUCLA real(pReal), parameter :: & kB = 1.38e-23_pReal !< Boltzmann constant in J/Kelvin diff --git a/src/constitutive_plastic_dislotwin.f90 b/src/constitutive_plastic_dislotwin.f90 index 09294ecd9..004238817 100644 --- a/src/constitutive_plastic_dislotwin.f90 +++ b/src/constitutive_plastic_dislotwin.f90 @@ -7,7 +7,7 @@ !> @brief material subroutine incoprorating dislocation and twinning physics !> @details to be done !-------------------------------------------------------------------------------------------------- -submodule(constitutive) plastic_dislotwin +submodule(constitutive:constitutive_plastic) plastic_dislotwin real(pReal), parameter :: & kB = 1.38e-23_pReal !< Boltzmann constant in J/Kelvin diff --git a/src/constitutive_plastic_isotropic.f90 b/src/constitutive_plastic_isotropic.f90 index 60f40fe37..3484b6f64 100644 --- a/src/constitutive_plastic_isotropic.f90 +++ b/src/constitutive_plastic_isotropic.f90 @@ -7,7 +7,7 @@ !! resolving the stress on the slip systems. Will give the response of phenopowerlaw for an !! untextured polycrystal !-------------------------------------------------------------------------------------------------- -submodule(constitutive) plastic_isotropic +submodule(constitutive:constitutive_plastic) plastic_isotropic type :: tParameters real(pReal) :: & diff --git a/src/constitutive_plastic_kinehardening.f90 b/src/constitutive_plastic_kinehardening.f90 index d3e84cfe2..163b5794a 100644 --- a/src/constitutive_plastic_kinehardening.f90 +++ b/src/constitutive_plastic_kinehardening.f90 @@ -5,7 +5,7 @@ !> @brief Phenomenological crystal plasticity using a power law formulation for the shear rates !! and a Voce-type kinematic hardening rule !-------------------------------------------------------------------------------------------------- -submodule(constitutive) plastic_kinehardening +submodule(constitutive:constitutive_plastic) plastic_kinehardening type :: tParameters real(pReal) :: & diff --git a/src/constitutive_plastic_none.f90 b/src/constitutive_plastic_none.f90 index da3ee9796..4e6033499 100644 --- a/src/constitutive_plastic_none.f90 +++ b/src/constitutive_plastic_none.f90 @@ -4,7 +4,7 @@ !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @brief Dummy plasticity for purely elastic material !-------------------------------------------------------------------------------------------------- -submodule(constitutive) plastic_none +submodule(constitutive:constitutive_plastic) plastic_none contains diff --git a/src/constitutive_plastic_nonlocal.f90 b/src/constitutive_plastic_nonlocal.f90 index ea55024b9..53b173db3 100644 --- a/src/constitutive_plastic_nonlocal.f90 +++ b/src/constitutive_plastic_nonlocal.f90 @@ -4,7 +4,7 @@ !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH !> @brief material subroutine for plasticity including dislocation flux !-------------------------------------------------------------------------------------------------- -submodule(constitutive) plastic_nonlocal +submodule(constitutive:constitutive_plastic) plastic_nonlocal use geometry_plastic_nonlocal, only: & nIPneighbors => geometry_plastic_nonlocal_nIPneighbors, & IPneighborhood => geometry_plastic_nonlocal_IPneighborhood, & diff --git a/src/constitutive_plastic_phenopowerlaw.f90 b/src/constitutive_plastic_phenopowerlaw.f90 index 53e55b319..3f6ade8d2 100644 --- a/src/constitutive_plastic_phenopowerlaw.f90 +++ b/src/constitutive_plastic_phenopowerlaw.f90 @@ -4,7 +4,7 @@ !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @brief phenomenological crystal plasticity formulation using a powerlaw fitting !-------------------------------------------------------------------------------------------------- -submodule(constitutive) plastic_phenopowerlaw +submodule(constitutive:constitutive_plastic) plastic_phenopowerlaw type :: tParameters real(pReal) :: & diff --git a/src/constitutive_thermal.f90 b/src/constitutive_thermal.f90 new file mode 100644 index 000000000..f1cd12d7f --- /dev/null +++ b/src/constitutive_thermal.f90 @@ -0,0 +1,112 @@ +submodule(constitutive) constitutive_thermal + + implicit none + + interface + + module subroutine source_thermal_dissipation_init + end subroutine source_thermal_dissipation_init + + module subroutine source_thermal_externalheat_init + end subroutine source_thermal_externalheat_init + + module subroutine kinematics_thermal_expansion_init + end subroutine kinematics_thermal_expansion_init + + module subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDot_dT, Tstar, Lp, phase) + + integer, intent(in) :: & + phase + real(pReal), intent(in), dimension(3,3) :: & + Tstar + real(pReal), intent(in), dimension(3,3) :: & + Lp + + real(pReal), intent(out) :: & + TDot, & + dTDot_dT + + end subroutine source_thermal_dissipation_getRateAndItsTangent + + module subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, dTDot_dT, phase, of) + + integer, intent(in) :: & + phase, & + of + real(pReal), intent(out) :: & + TDot, & + dTDot_dT + + end subroutine source_thermal_externalheat_getRateAndItsTangent + + end interface +contains + +module subroutine thermal_init + +! initialize source mechanisms + if (any(phase_source == SOURCE_thermal_dissipation_ID)) call source_thermal_dissipation_init + if (any(phase_source == SOURCE_thermal_externalheat_ID)) call source_thermal_externalheat_init + +!-------------------------------------------------------------------------------------------------- +!initialize kinematic mechanisms + if (any(phase_kinematics == KINEMATICS_thermal_expansion_ID)) call kinematics_thermal_expansion_init + +end subroutine thermal_init + + +module subroutine thermal_source_getRateAndItsTangents(Tdot, dTdot_dT, T, Tstar, Lp, ip, el) + + integer, intent(in) :: & + ip, & !< integration point number + el !< element number + real(pReal), intent(in) :: & + T + real(pReal), intent(in), dimension(:,:,:,:,:) :: & + Tstar, & + Lp + real(pReal), intent(inout) :: & + Tdot, & + dTdot_dT + + real(pReal) :: & + my_Tdot, & + my_dTdot_dT + integer :: & + phase, & + homog, & + instance, & + grain, & + source, & + constituent + + homog = material_homogenizationAt(el) + instance = thermal_typeInstance(homog) + + do grain = 1, homogenization_Ngrains(homog) + phase = material_phaseAt(grain,el) + constituent = material_phasememberAt(grain,ip,el) + do source = 1, phase_Nsources(phase) + select case(phase_source(source,phase)) + case (SOURCE_thermal_dissipation_ID) + call source_thermal_dissipation_getRateAndItsTangent(my_Tdot, my_dTdot_dT, & + Tstar(1:3,1:3,grain,ip,el), & + Lp(1:3,1:3,grain,ip,el), & + phase) + + case (SOURCE_thermal_externalheat_ID) + call source_thermal_externalheat_getRateAndItsTangent(my_Tdot, my_dTdot_dT, & + phase, constituent) + + case default + my_Tdot = 0.0_pReal + my_dTdot_dT = 0.0_pReal + end select + Tdot = Tdot + my_Tdot + dTdot_dT = dTdot_dT + my_dTdot_dT + enddo + enddo + +end subroutine thermal_source_getRateAndItsTangents + +end submodule diff --git a/src/damage_local.f90 b/src/damage_local.f90 index 85701acee..66b50064b 100644 --- a/src/damage_local.f90 +++ b/src/damage_local.f90 @@ -9,10 +9,7 @@ module damage_local use config use numerics use YAML_types - use source_damage_isoBrittle - use source_damage_isoDuctile - use source_damage_anisoBrittle - use source_damage_anisoDuctile + use constitutive use results implicit none @@ -128,42 +125,13 @@ subroutine damage_local_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el el !< element number real(pReal), intent(in) :: & phi - integer :: & - phase, & - grain, & - source, & - constituent real(pReal) :: & - phiDot, dPhiDot_dPhi, localphiDot, dLocalphiDot_dPhi + phiDot, dPhiDot_dPhi phiDot = 0.0_pReal dPhiDot_dPhi = 0.0_pReal - do grain = 1, homogenization_Ngrains(material_homogenizationAt(el)) - phase = material_phaseAt(grain,el) - constituent = material_phasememberAt(grain,ip,el) - do source = 1, phase_Nsources(phase) - select case(phase_source(source,phase)) - case (SOURCE_damage_isoBrittle_ID) - call source_damage_isobrittle_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) - - case (SOURCE_damage_isoDuctile_ID) - call source_damage_isoductile_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) - - case (SOURCE_damage_anisoBrittle_ID) - call source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) - - case (SOURCE_damage_anisoDuctile_ID) - call source_damage_anisoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) - - case default - localphiDot = 0.0_pReal - dLocalphiDot_dPhi = 0.0_pReal - - end select - phiDot = phiDot + localphiDot - dPhiDot_dPhi = dPhiDot_dPhi + dLocalphiDot_dPhi - enddo - enddo + + call constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el) phiDot = phiDot/real(homogenization_Ngrains(material_homogenizationAt(el)),pReal) dPhiDot_dPhi = dPhiDot_dPhi/real(homogenization_Ngrains(material_homogenizationAt(el)),pReal) diff --git a/src/damage_nonlocal.f90 b/src/damage_nonlocal.f90 index d3c3c07a2..914133de7 100644 --- a/src/damage_nonlocal.f90 +++ b/src/damage_nonlocal.f90 @@ -10,10 +10,7 @@ module damage_nonlocal use YAML_types use crystallite use lattice - use source_damage_isoBrittle - use source_damage_isoDuctile - use source_damage_anisoBrittle - use source_damage_anisoDuctile + use constitutive use results implicit none @@ -97,43 +94,13 @@ subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el !< element number real(pReal), intent(in) :: & phi - integer :: & - phase, & - grain, & - source, & - constituent real(pReal) :: & - phiDot, dPhiDot_dPhi, localphiDot, dLocalphiDot_dPhi + phiDot, dPhiDot_dPhi phiDot = 0.0_pReal dPhiDot_dPhi = 0.0_pReal - do grain = 1, homogenization_Ngrains(material_homogenizationAt(el)) - phase = material_phaseAt(grain,el) - constituent = material_phasememberAt(grain,ip,el) - do source = 1, phase_Nsources(phase) - select case(phase_source(source,phase)) - case (SOURCE_damage_isoBrittle_ID) - call source_damage_isobrittle_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) - - case (SOURCE_damage_isoDuctile_ID) - call source_damage_isoductile_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) - - case (SOURCE_damage_anisoBrittle_ID) - call source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) - - case (SOURCE_damage_anisoDuctile_ID) - call source_damage_anisoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) - - case default - localphiDot = 0.0_pReal - dLocalphiDot_dPhi = 0.0_pReal - - end select - phiDot = phiDot + localphiDot - dPhiDot_dPhi = dPhiDot_dPhi + dLocalphiDot_dPhi - enddo - enddo - + + call constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el) phiDot = phiDot/real(homogenization_Ngrains(material_homogenizationAt(el)),pReal) dPhiDot_dPhi = dPhiDot_dPhi/real(homogenization_Ngrains(material_homogenizationAt(el)),pReal) diff --git a/src/kinematics_cleavage_opening.f90 b/src/kinematics_cleavage_opening.f90 index e35f37e0e..2e314024c 100644 --- a/src/kinematics_cleavage_opening.f90 +++ b/src/kinematics_cleavage_opening.f90 @@ -4,16 +4,7 @@ !> @brief material subroutine incorporating kinematics resulting from opening of cleavage planes !> @details to be done !-------------------------------------------------------------------------------------------------- -module kinematics_cleavage_opening - use prec - use IO - use config - use math - use lattice - use material - - implicit none - private +submodule(constitutive:constitutive_damage) kinematics_cleavage_opening integer, dimension(:), allocatable :: kinematics_cleavage_opening_instance @@ -31,9 +22,6 @@ module kinematics_cleavage_opening type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance) - public :: & - kinematics_cleavage_opening_init, & - kinematics_cleavage_opening_LiAndItsTangent contains @@ -42,7 +30,7 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -subroutine kinematics_cleavage_opening_init +module subroutine kinematics_cleavage_opening_init integer :: Ninstance,p integer, dimension(:), allocatable :: N_cl !< active number of cleavage systems per family @@ -95,7 +83,7 @@ end subroutine kinematics_cleavage_opening_init !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the velocity gradient !-------------------------------------------------------------------------------------------------- -subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, ip, el) +module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, ip, el) integer, intent(in) :: & ipc, & !< grain number @@ -158,4 +146,4 @@ subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, i end subroutine kinematics_cleavage_opening_LiAndItsTangent -end module kinematics_cleavage_opening +end submodule kinematics_cleavage_opening diff --git a/src/kinematics_slipplane_opening.f90 b/src/kinematics_slipplane_opening.f90 index 847dc6c72..b15d206d3 100644 --- a/src/kinematics_slipplane_opening.f90 +++ b/src/kinematics_slipplane_opening.f90 @@ -4,16 +4,7 @@ !> @brief material subroutine incorporating kinematics resulting from opening of slip planes !> @details to be done !-------------------------------------------------------------------------------------------------- -module kinematics_slipplane_opening - use prec - use config - use IO - use math - use lattice - use material - - implicit none - private +submodule(constitutive:constitutive_damage) kinematics_slipplane_opening integer, dimension(:), allocatable :: kinematics_slipplane_opening_instance @@ -33,9 +24,6 @@ module kinematics_slipplane_opening type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance) - public :: & - kinematics_slipplane_opening_init, & - kinematics_slipplane_opening_LiAndItsTangent contains @@ -44,7 +32,7 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -subroutine kinematics_slipplane_opening_init +module subroutine kinematics_slipplane_opening_init integer :: Ninstance,p,i character(len=pStringLen) :: extmsg = '' @@ -107,7 +95,7 @@ end subroutine kinematics_slipplane_opening_init !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the velocity gradient !-------------------------------------------------------------------------------------------------- -subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, ip, el) +module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, ip, el) integer, intent(in) :: & ipc, & !< grain number @@ -183,4 +171,4 @@ subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, end subroutine kinematics_slipplane_opening_LiAndItsTangent -end module kinematics_slipplane_opening +end submodule kinematics_slipplane_opening diff --git a/src/kinematics_thermal_expansion.f90 b/src/kinematics_thermal_expansion.f90 index 63da6eb51..00871e94f 100644 --- a/src/kinematics_thermal_expansion.f90 +++ b/src/kinematics_thermal_expansion.f90 @@ -3,16 +3,7 @@ !> @brief material subroutine incorporating kinematics resulting from thermal expansion !> @details to be done !-------------------------------------------------------------------------------------------------- -module kinematics_thermal_expansion - use prec - use IO - use config - use math - use lattice - use material - - implicit none - private +submodule(constitutive:constitutive_thermal) kinematics_thermal_expansion integer, dimension(:), allocatable :: kinematics_thermal_expansion_instance @@ -25,10 +16,6 @@ module kinematics_thermal_expansion type(tParameters), dimension(:), allocatable :: param - public :: & - kinematics_thermal_expansion_init, & - kinematics_thermal_expansion_initialStrain, & - kinematics_thermal_expansion_LiAndItsTangent contains @@ -37,7 +24,7 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -subroutine kinematics_thermal_expansion_init +module subroutine kinematics_thermal_expansion_init integer :: Ninstance,p,i real(pReal), dimension(:), allocatable :: temp @@ -79,30 +66,30 @@ end subroutine kinematics_thermal_expansion_init !-------------------------------------------------------------------------------------------------- !> @brief report initial thermal strain based on current temperature deviation from reference !-------------------------------------------------------------------------------------------------- -pure 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) +! +! 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 !-------------------------------------------------------------------------------------------------- !> @brief constitutive equation for calculating the velocity gradient !-------------------------------------------------------------------------------------------------- -subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, ipc, ip, el) +module subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, ipc, ip, el) integer, intent(in) :: & ipc, & !< grain number @@ -140,4 +127,4 @@ subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, ipc, ip, end subroutine kinematics_thermal_expansion_LiAndItsTangent -end module kinematics_thermal_expansion +end submodule kinematics_thermal_expansion diff --git a/src/source_damage_anisoBrittle.f90 b/src/source_damage_anisoBrittle.f90 index c5ff1f906..20cf3a914 100644 --- a/src/source_damage_anisoBrittle.f90 +++ b/src/source_damage_anisoBrittle.f90 @@ -4,18 +4,7 @@ !> @brief material subroutine incorporating anisotropic brittle damage source mechanism !> @details to be done !-------------------------------------------------------------------------------------------------- -module source_damage_anisoBrittle - use prec - use IO - use math - use discretization - use material - use config - use lattice - use results - - implicit none - private +submodule (constitutive:constitutive_damage) source_damage_anisoBrittle integer, dimension(:), allocatable :: & source_damage_anisoBrittle_offset, & !< which source is my current source mechanism? @@ -39,12 +28,6 @@ module source_damage_anisoBrittle type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance) - public :: & - source_damage_anisoBrittle_init, & - source_damage_anisoBrittle_dotState, & - source_damage_anisobrittle_getRateAndItsTangent, & - source_damage_anisoBrittle_results - contains @@ -52,7 +35,7 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -subroutine source_damage_anisoBrittle_init +module subroutine source_damage_anisoBrittle_init integer :: Ninstance,sourceOffset,NipcMyPhase,p integer, dimension(:), allocatable :: N_cl @@ -123,7 +106,7 @@ end subroutine source_damage_anisoBrittle_init !-------------------------------------------------------------------------------------------------- !> @brief calculates derived quantities from state !-------------------------------------------------------------------------------------------------- -subroutine source_damage_anisoBrittle_dotState(S, ipc, ip, el) +module subroutine source_damage_anisoBrittle_dotState(S, ipc, ip, el) integer, intent(in) :: & ipc, & !< component-ID of integration point @@ -172,7 +155,7 @@ end subroutine source_damage_anisoBrittle_dotState !-------------------------------------------------------------------------------------------------- !> @brief returns local part of nonlocal damage driving force !-------------------------------------------------------------------------------------------------- -subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) +module subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) integer, intent(in) :: & phase, & @@ -199,7 +182,7 @@ end subroutine source_damage_anisoBrittle_getRateAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief writes results to HDF5 output file !-------------------------------------------------------------------------------------------------- -subroutine source_damage_anisoBrittle_results(phase,group) +module subroutine source_damage_anisoBrittle_results(phase,group) integer, intent(in) :: phase character(len=*), intent(in) :: group @@ -218,4 +201,4 @@ subroutine source_damage_anisoBrittle_results(phase,group) end subroutine source_damage_anisoBrittle_results -end module source_damage_anisoBrittle +end submodule source_damage_anisoBrittle diff --git a/src/source_damage_anisoDuctile.f90 b/src/source_damage_anisoDuctile.f90 index e8a76dc3a..3722f8b1c 100644 --- a/src/source_damage_anisoDuctile.f90 +++ b/src/source_damage_anisoDuctile.f90 @@ -4,23 +4,13 @@ !> @brief material subroutine incorporating anisotropic ductile damage source mechanism !> @details to be done !-------------------------------------------------------------------------------------------------- -module source_damage_anisoDuctile - use prec - use IO - use math - use discretization - use material - use config - use results - - implicit none - private +submodule(constitutive:constitutive_damage) source_damage_anisoDuctile integer, dimension(:), allocatable :: & source_damage_anisoDuctile_offset, & !< which source is my current damage mechanism? source_damage_anisoDuctile_instance !< instance of damage source mechanism - type, private :: tParameters !< container type for internal constitutive parameters + type :: tParameters !< container type for internal constitutive parameters real(pReal) :: & n real(pReal), dimension(:), allocatable :: & @@ -29,14 +19,7 @@ module source_damage_anisoDuctile output end type tParameters - type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance) - - - public :: & - source_damage_anisoDuctile_init, & - source_damage_anisoDuctile_dotState, & - source_damage_anisoDuctile_getRateAndItsTangent, & - source_damage_anisoDuctile_results + type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance) contains @@ -45,7 +28,7 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -subroutine source_damage_anisoDuctile_init +module subroutine source_damage_anisoDuctile_init integer :: Ninstance,sourceOffset,NipcMyPhase,p integer, dimension(:), allocatable :: N_sl @@ -105,7 +88,7 @@ end subroutine source_damage_anisoDuctile_init !-------------------------------------------------------------------------------------------------- !> @brief calculates derived quantities from state !-------------------------------------------------------------------------------------------------- -subroutine source_damage_anisoDuctile_dotState(ipc, ip, el) +module subroutine source_damage_anisoDuctile_dotState(ipc, ip, el) integer, intent(in) :: & ipc, & !< component-ID of integration point @@ -136,7 +119,7 @@ end subroutine source_damage_anisoDuctile_dotState !-------------------------------------------------------------------------------------------------- !> @brief returns local part of nonlocal damage driving force !-------------------------------------------------------------------------------------------------- -subroutine source_damage_anisoDuctile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) +module subroutine source_damage_anisoDuctile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) integer, intent(in) :: & phase, & @@ -163,7 +146,7 @@ end subroutine source_damage_anisoDuctile_getRateAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief writes results to HDF5 output file !-------------------------------------------------------------------------------------------------- -subroutine source_damage_anisoDuctile_results(phase,group) +module subroutine source_damage_anisoDuctile_results(phase,group) integer, intent(in) :: phase character(len=*), intent(in) :: group @@ -182,4 +165,4 @@ subroutine source_damage_anisoDuctile_results(phase,group) end subroutine source_damage_anisoDuctile_results -end module source_damage_anisoDuctile +end submodule source_damage_anisoDuctile diff --git a/src/source_damage_isoBrittle.f90 b/src/source_damage_isoBrittle.f90 index c4c4c72a4..00704fe26 100644 --- a/src/source_damage_isoBrittle.f90 +++ b/src/source_damage_isoBrittle.f90 @@ -4,23 +4,13 @@ !> @brief material subroutine incoprorating isotropic brittle damage source mechanism !> @details to be done !-------------------------------------------------------------------------------------------------- -module source_damage_isoBrittle - use prec - use IO - use math - use discretization - use material - use config - use results - - implicit none - private +submodule(constitutive:constitutive_damage) source_damage_isoBrittle integer, dimension(:), allocatable :: & source_damage_isoBrittle_offset, & source_damage_isoBrittle_instance - type, private :: tParameters !< container type for internal constitutive parameters + type :: tParameters !< container type for internal constitutive parameters real(pReal) :: & critStrainEnergy, & N @@ -30,13 +20,6 @@ module source_damage_isoBrittle type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance) - - public :: & - source_damage_isoBrittle_init, & - source_damage_isoBrittle_deltaState, & - source_damage_isoBrittle_getRateAndItsTangent, & - source_damage_isoBrittle_results - contains @@ -44,7 +27,7 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -subroutine source_damage_isoBrittle_init +module subroutine source_damage_isoBrittle_init integer :: Ninstance,sourceOffset,NipcMyPhase,p character(len=pStringLen) :: extmsg = '' @@ -99,7 +82,7 @@ end subroutine source_damage_isoBrittle_init !-------------------------------------------------------------------------------------------------- !> @brief calculates derived quantities from state !-------------------------------------------------------------------------------------------------- -subroutine source_damage_isoBrittle_deltaState(C, Fe, ipc, ip, el) +module subroutine source_damage_isoBrittle_deltaState(C, Fe, ipc, ip, el) integer, intent(in) :: & ipc, & !< component-ID of integration point @@ -145,7 +128,7 @@ end subroutine source_damage_isoBrittle_deltaState !-------------------------------------------------------------------------------------------------- !> @brief returns local part of nonlocal damage driving force !-------------------------------------------------------------------------------------------------- -subroutine source_damage_isoBrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) +module subroutine source_damage_isoBrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) integer, intent(in) :: & phase, & @@ -174,7 +157,7 @@ end subroutine source_damage_isoBrittle_getRateAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief writes results to HDF5 output file !-------------------------------------------------------------------------------------------------- -subroutine source_damage_isoBrittle_results(phase,group) +module subroutine source_damage_isoBrittle_results(phase,group) integer, intent(in) :: phase character(len=*), intent(in) :: group @@ -193,4 +176,4 @@ subroutine source_damage_isoBrittle_results(phase,group) end subroutine source_damage_isoBrittle_results -end module source_damage_isoBrittle +end submodule source_damage_isoBrittle diff --git a/src/source_damage_isoDuctile.f90 b/src/source_damage_isoDuctile.f90 index cf392fdc4..517332316 100644 --- a/src/source_damage_isoDuctile.f90 +++ b/src/source_damage_isoDuctile.f90 @@ -4,22 +4,13 @@ !> @brief material subroutine incorporating isotropic ductile damage source mechanism !> @details to be done !-------------------------------------------------------------------------------------------------- -module source_damage_isoDuctile - use prec - use IO - use discretization - use material - use config - use results - - implicit none - private +submodule (constitutive:constitutive_damage) source_damage_isoDuctile integer, dimension(:), allocatable :: & source_damage_isoDuctile_offset, & !< which source is my current damage mechanism? source_damage_isoDuctile_instance !< instance of damage source mechanism - type, private :: tParameters !< container type for internal constitutive parameters + type:: tParameters !< container type for internal constitutive parameters real(pReal) :: & critPlasticStrain, & N @@ -27,15 +18,9 @@ module source_damage_isoDuctile output end type tParameters - type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance) + type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance) - public :: & - source_damage_isoDuctile_init, & - source_damage_isoDuctile_dotState, & - source_damage_isoDuctile_getRateAndItsTangent, & - source_damage_isoDuctile_Results - contains @@ -43,7 +28,7 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -subroutine source_damage_isoDuctile_init +module subroutine source_damage_isoDuctile_init integer :: Ninstance,sourceOffset,NipcMyPhase,p character(len=pStringLen) :: extmsg = '' @@ -98,7 +83,7 @@ end subroutine source_damage_isoDuctile_init !-------------------------------------------------------------------------------------------------- !> @brief calculates derived quantities from state !-------------------------------------------------------------------------------------------------- -subroutine source_damage_isoDuctile_dotState(ipc, ip, el) +module subroutine source_damage_isoDuctile_dotState(ipc, ip, el) integer, intent(in) :: & ipc, & !< component-ID of integration point @@ -129,7 +114,7 @@ end subroutine source_damage_isoDuctile_dotState !-------------------------------------------------------------------------------------------------- !> @brief returns local part of nonlocal damage driving force !-------------------------------------------------------------------------------------------------- -subroutine source_damage_isoDuctile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) +module subroutine source_damage_isoDuctile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) integer, intent(in) :: & phase, & @@ -156,7 +141,7 @@ end subroutine source_damage_isoDuctile_getRateAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief writes results to HDF5 output file !-------------------------------------------------------------------------------------------------- -subroutine source_damage_isoDuctile_results(phase,group) +module subroutine source_damage_isoDuctile_results(phase,group) integer, intent(in) :: phase character(len=*), intent(in) :: group @@ -175,4 +160,4 @@ subroutine source_damage_isoDuctile_results(phase,group) end subroutine source_damage_isoDuctile_results -end module source_damage_isoDuctile +end submodule source_damage_isoDuctile diff --git a/src/source_thermal_dissipation.f90 b/src/source_thermal_dissipation.f90 index 8fa5ae0c7..58a0c6b3c 100644 --- a/src/source_thermal_dissipation.f90 +++ b/src/source_thermal_dissipation.f90 @@ -4,14 +4,7 @@ !> @brief material subroutine for thermal source due to plastic dissipation !> @details to be done !-------------------------------------------------------------------------------------------------- -module source_thermal_dissipation - use prec - use discretization - use material - use config - - implicit none - private +submodule(constitutive:constitutive_thermal) source_thermal_dissipation integer, dimension(:), allocatable :: & source_thermal_dissipation_offset, & !< which source is my current thermal dissipation mechanism? @@ -25,10 +18,6 @@ module source_thermal_dissipation type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance) - public :: & - source_thermal_dissipation_init, & - source_thermal_dissipation_getRateAndItsTangent - contains @@ -36,7 +25,7 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -subroutine source_thermal_dissipation_init +module subroutine source_thermal_dissipation_init integer :: Ninstance,sourceOffset,NipcMyPhase,p @@ -76,7 +65,7 @@ end subroutine source_thermal_dissipation_init !-------------------------------------------------------------------------------------------------- !> @brief Ninstances dissipation rate !-------------------------------------------------------------------------------------------------- -subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDot_dT, Tstar, Lp, phase) +module subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDot_dT, Tstar, Lp, phase) integer, intent(in) :: & phase @@ -96,4 +85,4 @@ subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDot_dT, Tstar end subroutine source_thermal_dissipation_getRateAndItsTangent -end module source_thermal_dissipation +end submodule source_thermal_dissipation diff --git a/src/source_thermal_externalheat.f90 b/src/source_thermal_externalheat.f90 index b83b29596..8ffc7a4fb 100644 --- a/src/source_thermal_externalheat.f90 +++ b/src/source_thermal_externalheat.f90 @@ -4,14 +4,8 @@ !> @author Philip Eisenlohr, Michigan State University !> @brief material subroutine for variable heat source !-------------------------------------------------------------------------------------------------- -module source_thermal_externalheat - use prec - use discretization - use material - use config +submodule(constitutive:constitutive_thermal) source_thermal_externalheat - implicit none - private integer, dimension(:), allocatable :: & source_thermal_externalheat_offset, & !< which source is my current thermal dissipation mechanism? @@ -28,11 +22,6 @@ module source_thermal_externalheat type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance) - public :: & - source_thermal_externalheat_init, & - source_thermal_externalheat_dotState, & - source_thermal_externalheat_getRateAndItsTangent - contains @@ -40,7 +29,7 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -subroutine source_thermal_externalheat_init +module subroutine source_thermal_externalheat_init integer :: Ninstance,sourceOffset,NipcMyPhase,p @@ -84,7 +73,7 @@ end subroutine source_thermal_externalheat_init !> @brief rate of change of state !> @details state only contains current time to linearly interpolate given heat powers !-------------------------------------------------------------------------------------------------- -subroutine source_thermal_externalheat_dotState(phase, of) +module subroutine source_thermal_externalheat_dotState(phase, of) integer, intent(in) :: & phase, & @@ -103,7 +92,7 @@ end subroutine source_thermal_externalheat_dotState !-------------------------------------------------------------------------------------------------- !> @brief returns local heat generation rate !-------------------------------------------------------------------------------------------------- -subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, dTDot_dT, phase, of) +module subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, dTDot_dT, phase, of) integer, intent(in) :: & phase, & @@ -135,4 +124,4 @@ subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, dTDot_dT, phas end subroutine source_thermal_externalheat_getRateAndItsTangent -end module source_thermal_externalheat +end submodule source_thermal_externalheat diff --git a/src/thermal_adiabatic.f90 b/src/thermal_adiabatic.f90 index a11e7b0a1..c52f0a3d0 100644 --- a/src/thermal_adiabatic.f90 +++ b/src/thermal_adiabatic.f90 @@ -8,8 +8,7 @@ module thermal_adiabatic use numerics use material use results - use source_thermal_dissipation - use source_thermal_externalheat + use constitutive use crystallite use lattice @@ -125,46 +124,15 @@ subroutine thermal_adiabatic_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el) T real(pReal), intent(out) :: & Tdot, dTdot_dT - - real(pReal) :: & - my_Tdot, my_dTdot_dT integer :: & - phase, & - homog, & - instance, & - grain, & - source, & - constituent - - homog = material_homogenizationAt(el) - instance = thermal_typeInstance(homog) - + homog + Tdot = 0.0_pReal dTdot_dT = 0.0_pReal - do grain = 1, homogenization_Ngrains(homog) - phase = material_phaseAt(grain,el) - constituent = material_phasememberAt(grain,ip,el) - do source = 1, phase_Nsources(phase) - select case(phase_source(source,phase)) - case (SOURCE_thermal_dissipation_ID) - call source_thermal_dissipation_getRateAndItsTangent(my_Tdot, my_dTdot_dT, & - crystallite_S(1:3,1:3,grain,ip,el), & - crystallite_Lp(1:3,1:3,grain,ip,el), & - phase) - case (SOURCE_thermal_externalheat_ID) - call source_thermal_externalheat_getRateAndItsTangent(my_Tdot, my_dTdot_dT, & - phase, constituent) - - case default - my_Tdot = 0.0_pReal - my_dTdot_dT = 0.0_pReal - end select - Tdot = Tdot + my_Tdot - dTdot_dT = dTdot_dT + my_dTdot_dT - enddo - enddo - + homog = material_homogenizationAt(el) + call constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, crystallite_S, crystallite_Lp, ip, el) + Tdot = Tdot/real(homogenization_Ngrains(homog),pReal) dTdot_dT = dTdot_dT/real(homogenization_Ngrains(homog),pReal) diff --git a/src/thermal_conduction.f90 b/src/thermal_conduction.f90 index 82eed4cc4..60766710d 100644 --- a/src/thermal_conduction.f90 +++ b/src/thermal_conduction.f90 @@ -9,8 +9,7 @@ module thermal_conduction use lattice use results use crystallite - use source_thermal_dissipation - use source_thermal_externalheat + use constitutive implicit none private @@ -84,46 +83,14 @@ subroutine thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el) T real(pReal), intent(out) :: & Tdot, dTdot_dT - real(pReal) :: & - my_Tdot, my_dTdot_dT integer :: & - phase, & - homog, & - offset, & - instance, & - grain, & - source, & - constituent - - homog = material_homogenizationAt(el) - offset = material_homogenizationMemberAt(ip,el) - instance = thermal_typeInstance(homog) - + homog + Tdot = 0.0_pReal dTdot_dT = 0.0_pReal - do grain = 1, homogenization_Ngrains(homog) - phase = material_phaseAt(grain,el) - constituent = material_phasememberAt(grain,ip,el) - do source = 1, phase_Nsources(phase) - select case(phase_source(source,phase)) - case (SOURCE_thermal_dissipation_ID) - call source_thermal_dissipation_getRateAndItsTangent(my_Tdot, my_dTdot_dT, & - crystallite_S(1:3,1:3,grain,ip,el), & - crystallite_Lp(1:3,1:3,grain,ip,el), & - phase) - case (SOURCE_thermal_externalheat_ID) - call source_thermal_externalheat_getRateAndItsTangent(my_Tdot, my_dTdot_dT, & - phase, constituent) - case default - my_Tdot = 0.0_pReal - my_dTdot_dT = 0.0_pReal - - end select - Tdot = Tdot + my_Tdot - dTdot_dT = dTdot_dT + my_dTdot_dT - enddo - enddo + homog = material_homogenizationAt(el) + call constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, crystallite_S,crystallite_Lp ,ip, el) Tdot = Tdot/real(homogenization_Ngrains(homog),pReal) dTdot_dT = dTdot_dT/real(homogenization_Ngrains(homog),pReal) From bc1d73c03b5bdbd11f7a833eaddcda91e0aa2d50 Mon Sep 17 00:00:00 2001 From: Sharan Roongta Date: Thu, 9 Jul 2020 01:49:48 +0200 Subject: [PATCH 02/21] trying new structure for all constitutive modules --- src/commercialFEM_fileList.f90 | 19 ++++++++------- src/constitutive.f90 | 19 ++++++++------- src/kinematics_thermal_expansion.f90 | 36 ++++++++++++++-------------- 3 files changed, 40 insertions(+), 34 deletions(-) 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 !-------------------------------------------------------------------------------------------------- From fd7110ce4516c4a87b345f6b095bf86a07ee00a9 Mon Sep 17 00:00:00 2001 From: Sharan Roongta Date: Fri, 10 Jul 2020 14:59:07 +0200 Subject: [PATCH 03/21] probably a more readable structure --- src/constitutive.f90 | 218 +++++++++++------------------------ src/constitutive_damage.f90 | 71 ++++++++++++ src/constitutive_plastic.f90 | 153 ++++++++++++++++++++++++ src/constitutive_thermal.f90 | 45 ++++++++ 4 files changed, 335 insertions(+), 152 deletions(-) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 2af434df5..e1e6ae47d 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -46,15 +46,15 @@ module constitutive end function kinematics_thermal_expansion_initialStrain - module function plastic_dislotwin_homogenizedC(ipc,ip,el) result(homogenizedC) + module function plastic_homogenizedC(ipc,ip,el) result(homogenizedC) real(pReal), dimension(6,6) :: & homogenizedC integer, intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point el !< element - end function plastic_dislotwin_homogenizedC - + end function plastic_homogenizedC + module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,of) real(pReal), dimension(3,3), intent(out) :: & @@ -109,101 +109,65 @@ module constitutive dLi_dTstar !< derivative of Li with respect to Tstar (4th-order tensor defined to be zero) end subroutine kinematics_thermal_expansion_LiAndItsTangent - module subroutine plastic_isotropic_dotState(Mp,instance,of) - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - integer, intent(in) :: & - instance, & - of - end subroutine plastic_isotropic_dotState - - module subroutine plastic_phenopowerlaw_dotState(Mp,instance,of) - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - integer, intent(in) :: & - instance, & - of - end subroutine plastic_phenopowerlaw_dotState - - module subroutine plastic_kinehardening_dotState(Mp,instance,of) - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - integer, intent(in) :: & - instance, & - of - end subroutine plastic_kinehardening_dotState - - module subroutine plastic_dislotwin_dotState(Mp,T,instance,of) - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - real(pReal), intent(in) :: & - T - integer, intent(in) :: & - instance, & - of - end subroutine plastic_dislotwin_dotState - - module subroutine plastic_disloUCLA_dotState(Mp,T,instance,of) - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - real(pReal), intent(in) :: & - T - integer, intent(in) :: & - instance, & - of - end subroutine plastic_disloUCLA_dotState - - module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature,timestep, & - instance,of,ip,el) - real(pReal), dimension(3,3), intent(in) ::& - Mp !< MandelStress - real(pReal), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem), intent(in) :: & - F, & !< deformation gradient - Fp !< plastic deformation gradient - real(pReal), intent(in) :: & - Temperature, & !< temperature - timestep !< substepped crystallite time increment - integer, intent(in) :: & - instance, & - of, & - ip, & !< current integration point - el !< current element number - end subroutine plastic_nonlocal_dotState - - - module subroutine source_damage_anisoBrittle_dotState(S, ipc, ip, el) + module function plastic_dotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el,phase,of) result(broken_plastic) integer, intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point - el !< element - real(pReal), intent(in), dimension(3,3) :: & - S - end subroutine source_damage_anisoBrittle_dotState - - module subroutine source_damage_anisoDuctile_dotState(ipc, ip, el) - - integer, intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element - end subroutine source_damage_anisoDuctile_dotState - - module subroutine source_damage_isoDuctile_dotState(ipc, ip, el) - - integer, intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element - end subroutine source_damage_isoDuctile_dotState - - module subroutine source_thermal_externalheat_dotState(phase, of) - - integer, intent(in) :: & + el, & !< element phase, & of - end subroutine source_thermal_externalheat_dotState + real(pReal), intent(in) :: & + subdt !< timestep + real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: & + FArray, & !< elastic deformation gradient + FpArray !< plastic deformation gradient + real(pReal), intent(in), dimension(3,3) :: & + Fi !< intermediate deformation gradient + real(pReal), intent(in), dimension(3,3) :: & + S !< 2nd Piola Kirchhoff stress (vector notation) + logical :: broken_plastic + end function plastic_dotState + module function damage_dotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el,phase,of) result(broken_damage) + + integer, intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el, & !< element + phase, & + of + real(pReal), intent(in) :: & + subdt !< timestep + real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: & + FArray, & !< elastic deformation gradient + FpArray !< plastic deformation gradient + real(pReal), intent(in), dimension(3,3) :: & + Fi !< intermediate deformation gradient + real(pReal), intent(in), dimension(3,3) :: & + S !< 2nd Piola Kirchhoff stress (vector notation) + logical :: broken_damage + end function damage_dotState + + module function thermal_dotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el,phase,of) result(broken_thermal) + + integer, intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el, & !< element + phase, & + of + real(pReal), intent(in) :: & + subdt !< timestep + real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: & + FArray, & !< elastic deformation gradient + FpArray !< plastic deformation gradient + real(pReal), intent(in), dimension(3,3) :: & + Fi !< intermediate deformation gradient + real(pReal), intent(in), dimension(3,3) :: & + S !< 2nd Piola Kirchhoff stress (vector notation) + logical :: broken_thermal + end function thermal_dotState module subroutine plastic_kinehardening_deltaState(Mp,instance,of) real(pReal), dimension(3,3), intent(in) :: & @@ -478,12 +442,7 @@ function constitutive_homogenizedC(ipc,ip,el) ip, & !< integration point el !< element - plasticityType: select case (phase_plasticity(material_phaseAt(ipc,el))) - case (PLASTICITY_DISLOTWIN_ID) plasticityType - constitutive_homogenizedC = plastic_dislotwin_homogenizedC(ipc,ip,el) - case default plasticityType - constitutive_homogenizedC = lattice_C66(1:6,1:6,material_phaseAt(ipc,el)) - end select plasticityType + constitutive_homogenizedC = plastic_homogenizedC(ipc,ip,el) end function constitutive_homogenizedC @@ -695,65 +654,20 @@ function constitutive_collectDotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el Fi !< intermediate deformation gradient real(pReal), intent(in), dimension(3,3) :: & S !< 2nd Piola Kirchhoff stress (vector notation) - real(pReal), dimension(3,3) :: & - Mp - integer :: & - ho, & !< homogenization - tme, & !< thermal member position - i, & !< counter in source loop - instance logical :: broken + logical :: & + broken_plastic, & + broken_thermal, & + broken_damage - ho = material_homogenizationAt(el) - tme = thermalMapping(ho)%p(ip,el) - instance = phase_plasticityInstance(phase) - Mp = matmul(matmul(transpose(Fi),Fi),S) + broken_plastic = plastic_dotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el,phase,of) + broken = broken_plastic + broken_damage = damage_dotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el,phase,of) + broken = broken .or. broken_damage + broken_thermal = thermal_dotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el,phase,of) + broken = broken .or. broken_thermal - plasticityType: select case (phase_plasticity(phase)) - - case (PLASTICITY_ISOTROPIC_ID) plasticityType - call plastic_isotropic_dotState (Mp,instance,of) - - case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType - call plastic_phenopowerlaw_dotState(Mp,instance,of) - - case (PLASTICITY_KINEHARDENING_ID) plasticityType - call plastic_kinehardening_dotState(Mp,instance,of) - - case (PLASTICITY_DISLOTWIN_ID) plasticityType - call plastic_dislotwin_dotState (Mp,temperature(ho)%p(tme),instance,of) - - case (PLASTICITY_DISLOUCLA_ID) plasticityType - call plastic_disloucla_dotState (Mp,temperature(ho)%p(tme),instance,of) - - case (PLASTICITY_NONLOCAL_ID) plasticityType - call plastic_nonlocal_dotState (Mp,FArray,FpArray,temperature(ho)%p(tme),subdt, & - instance,of,ip,el) - end select plasticityType - broken = any(IEEE_is_NaN(plasticState(phase)%dotState(:,of))) - - SourceLoop: do i = 1, phase_Nsources(phase) - - sourceType: select case (phase_source(i,phase)) - - case (SOURCE_damage_anisoBrittle_ID) sourceType - call source_damage_anisoBrittle_dotState (S, ipc, ip, el) !< correct stress? - - case (SOURCE_damage_isoDuctile_ID) sourceType - call source_damage_isoDuctile_dotState ( ipc, ip, el) - - case (SOURCE_damage_anisoDuctile_ID) sourceType - call source_damage_anisoDuctile_dotState ( ipc, ip, el) - - case (SOURCE_thermal_externalheat_ID) sourceType - call source_thermal_externalheat_dotState(phase,of) - - end select sourceType - - broken = broken .or. any(IEEE_is_NaN(sourceState(phase)%p(i)%dotState(:,of))) - - enddo SourceLoop end function constitutive_collectDotState diff --git a/src/constitutive_damage.f90 b/src/constitutive_damage.f90 index cb81c252d..a6264d57d 100644 --- a/src/constitutive_damage.f90 +++ b/src/constitutive_damage.f90 @@ -22,6 +22,33 @@ submodule(constitutive) constitutive_damage module subroutine kinematics_slipplane_opening_init end subroutine kinematics_slipplane_opening_init + module subroutine source_damage_anisoBrittle_dotState(S, ipc, ip, el) + + integer, intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + real(pReal), intent(in), dimension(3,3) :: & + S + end subroutine source_damage_anisoBrittle_dotState + + module subroutine source_damage_anisoDuctile_dotState(ipc, ip, el) + + integer, intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + end subroutine source_damage_anisoDuctile_dotState + + module subroutine source_damage_isoDuctile_dotState(ipc, ip, el) + + integer, intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + end subroutine source_damage_isoDuctile_dotState + + module subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) integer, intent(in) :: & @@ -107,6 +134,50 @@ module subroutine damage_init end subroutine damage_init +!-------------------------------------------------------------------------------------------------- +!> @brief contains the constitutive equation for calculating the rate of change of microstructure +!-------------------------------------------------------------------------------------------------- +module function damage_dotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el,phase,of) result(broken_damage) + + integer, intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el, & !< element + phase, & + of + real(pReal), intent(in) :: & + subdt !< timestep + real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: & + FArray, & !< elastic deformation gradient + FpArray !< plastic deformation gradient + real(pReal), intent(in), dimension(3,3) :: & + Fi !< intermediate deformation gradient + real(pReal), intent(in), dimension(3,3) :: & + S !< 2nd Piola Kirchhoff stress (vector notation) + logical :: broken_damage + integer :: i + + SourceLoop: do i = 1, phase_Nsources(phase) + + sourceType: select case (phase_source(i,phase)) + + case (SOURCE_damage_anisoBrittle_ID) sourceType + call source_damage_anisoBrittle_dotState (S, ipc, ip, el) !< correct stress? + + case (SOURCE_damage_isoDuctile_ID) sourceType + call source_damage_isoDuctile_dotState ( ipc, ip, el) + + case (SOURCE_damage_anisoDuctile_ID) sourceType + call source_damage_anisoDuctile_dotState ( ipc, ip, el) + + end select sourceType + + enddo sourceLoop + + broken_damage = any(IEEE_is_NaN(sourceState(phase)%p(i)%dotState(:,of))) + +end function damage_dotState + module subroutine damage_source_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el) diff --git a/src/constitutive_plastic.f90 b/src/constitutive_plastic.f90 index d4d0f19a2..5d190758a 100644 --- a/src/constitutive_plastic.f90 +++ b/src/constitutive_plastic.f90 @@ -25,6 +25,67 @@ submodule(constitutive) constitutive_plastic module subroutine plastic_nonlocal_init end subroutine plastic_nonlocal_init + module subroutine plastic_isotropic_dotState(Mp,instance,of) + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + integer, intent(in) :: & + instance, & + of + end subroutine plastic_isotropic_dotState + + module subroutine plastic_phenopowerlaw_dotState(Mp,instance,of) + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + integer, intent(in) :: & + instance, & + of + end subroutine plastic_phenopowerlaw_dotState + + module subroutine plastic_kinehardening_dotState(Mp,instance,of) + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + integer, intent(in) :: & + instance, & + of + end subroutine plastic_kinehardening_dotState + + module subroutine plastic_dislotwin_dotState(Mp,T,instance,of) + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + real(pReal), intent(in) :: & + T + integer, intent(in) :: & + instance, & + of + end subroutine plastic_dislotwin_dotState + + module subroutine plastic_disloUCLA_dotState(Mp,T,instance,of) + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + real(pReal), intent(in) :: & + T + integer, intent(in) :: & + instance, & + of + end subroutine plastic_disloUCLA_dotState + + module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature,timestep, & + instance,of,ip,el) + real(pReal), dimension(3,3), intent(in) ::& + Mp !< MandelStress + real(pReal), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem), intent(in) :: & + F, & !< deformation gradient + Fp !< plastic deformation gradient + real(pReal), intent(in) :: & + Temperature, & !< temperature + timestep !< substepped crystallite time increment + integer, intent(in) :: & + instance, & + of, & + ip, & !< current integration point + el !< current element number + end subroutine plastic_nonlocal_dotState + module subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) real(pReal), dimension(3,3), intent(out) :: & @@ -113,6 +174,15 @@ submodule(constitutive) constitutive_plastic el !< current element number end subroutine plastic_nonlocal_LpAndItsTangent + module function plastic_dislotwin_homogenizedC(ipc,ip,el) result(homogenizedC) + real(pReal), dimension(6,6) :: & + homogenizedC + integer, intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + end function plastic_dislotwin_homogenizedC + module subroutine plastic_dislotwin_dependentState(T,instance,of) integer, intent(in) :: & @@ -166,6 +236,89 @@ module subroutine plastic_init end subroutine plastic_init +!-------------------------------------------------------------------------------------------------- +!> @brief contains the constitutive equation for calculating the rate of change of microstructure +!-------------------------------------------------------------------------------------------------- +module function plastic_dotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el,phase,of) result(broken_plastic) + + integer, intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el, & !< element + phase, & + of + real(pReal), intent(in) :: & + subdt !< timestep + real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: & + FArray, & !< elastic deformation gradient + FpArray !< plastic deformation gradient + real(pReal), intent(in), dimension(3,3) :: & + Fi !< intermediate deformation gradient + real(pReal), intent(in), dimension(3,3) :: & + S !< 2nd Piola Kirchhoff stress (vector notation) + real(pReal), dimension(3,3) :: & + Mp + integer :: & + ho, & !< homogenization + tme, & !< thermal member position + i, & !< counter in source loop + instance + logical :: broken_plastic + + + ho = material_homogenizationAt(el) + tme = thermalMapping(ho)%p(ip,el) + instance = phase_plasticityInstance(phase) + + Mp = matmul(matmul(transpose(Fi),Fi),S) + + plasticityType: select case (phase_plasticity(phase)) + + case (PLASTICITY_ISOTROPIC_ID) plasticityType + call plastic_isotropic_dotState (Mp,instance,of) + + case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType + call plastic_phenopowerlaw_dotState(Mp,instance,of) + + case (PLASTICITY_KINEHARDENING_ID) plasticityType + call plastic_kinehardening_dotState(Mp,instance,of) + + case (PLASTICITY_DISLOTWIN_ID) plasticityType + call plastic_dislotwin_dotState (Mp,temperature(ho)%p(tme),instance,of) + + case (PLASTICITY_DISLOUCLA_ID) plasticityType + call plastic_disloucla_dotState (Mp,temperature(ho)%p(tme),instance,of) + + case (PLASTICITY_NONLOCAL_ID) plasticityType + call plastic_nonlocal_dotState (Mp,FArray,FpArray,temperature(ho)%p(tme),subdt, & + instance,of,ip,el) + end select plasticityType + broken_plastic = any(IEEE_is_NaN(plasticState(phase)%dotState(:,of))) + +end function plastic_dotState + + +!-------------------------------------------------------------------------------------------------- +!> @brief returns the homogenize elasticity matrix +!> ToDo: homogenizedC66 would be more consistent +!-------------------------------------------------------------------------------------------------- +module function plastic_homogenizedC(ipc,ip,el) result(homogenizedC) + + real(pReal), dimension(6,6) :: homogenizedC + integer, intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + + plasticityType: select case (phase_plasticity(material_phaseAt(ipc,el))) + case (PLASTICITY_DISLOTWIN_ID) plasticityType + homogenizedC = plastic_dislotwin_homogenizedC(ipc,ip,el) + case default plasticityType + homogenizedC = lattice_C66(1:6,1:6,material_phaseAt(ipc,el)) + end select plasticityType + +end function plastic_homogenizedC + !-------------------------------------------------------------------------------------------------- !> @brief calls microstructure function of the different constitutive models diff --git a/src/constitutive_thermal.f90 b/src/constitutive_thermal.f90 index f1cd12d7f..7891ebba4 100644 --- a/src/constitutive_thermal.f90 +++ b/src/constitutive_thermal.f90 @@ -13,6 +13,13 @@ submodule(constitutive) constitutive_thermal module subroutine kinematics_thermal_expansion_init end subroutine kinematics_thermal_expansion_init + module subroutine source_thermal_externalheat_dotState(phase, of) + integer, intent(in) :: & + phase, & + of + end subroutine source_thermal_externalheat_dotState + + module subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDot_dT, Tstar, Lp, phase) integer, intent(in) :: & @@ -54,6 +61,44 @@ module subroutine thermal_init end subroutine thermal_init +!-------------------------------------------------------------------------------------------------- +!> @brief contains the constitutive equation for calculating the rate of change of microstructure +!-------------------------------------------------------------------------------------------------- +module function thermal_dotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el,phase,of) result(broken_thermal) + + integer, intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el, & !< element + phase, & + of + real(pReal), intent(in) :: & + subdt !< timestep + real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: & + FArray, & !< elastic deformation gradient + FpArray !< plastic deformation gradient + real(pReal), intent(in), dimension(3,3) :: & + Fi !< intermediate deformation gradient + real(pReal), intent(in), dimension(3,3) :: & + S !< 2nd Piola Kirchhoff stress (vector notation) + logical :: broken_thermal + integer :: i + + SourceLoop: do i = 1, phase_Nsources(phase) + + sourceType: select case (phase_source(i,phase)) + + case (SOURCE_thermal_externalheat_ID) sourceType + call source_thermal_externalheat_dotState(phase,of) + + end select sourceType + + broken_thermal = any(IEEE_is_NaN(sourceState(phase)%p(i)%dotState(:,of))) + + enddo sourceLoop + +end function thermal_dotState + module subroutine thermal_source_getRateAndItsTangents(Tdot, dTdot_dT, T, Tstar, Lp, ip, el) From 957c51fb0783094f34d48e412ecf72ae31fd7ba1 Mon Sep 17 00:00:00 2001 From: Sharan Roongta Date: Fri, 10 Jul 2020 15:13:56 +0200 Subject: [PATCH 04/21] cleaner --- src/constitutive_damage.f90 | 33 +++---------------- src/constitutive_plastic.f90 | 61 ++++++------------------------------ src/constitutive_thermal.f90 | 40 ++++------------------- 3 files changed, 19 insertions(+), 115 deletions(-) diff --git a/src/constitutive_damage.f90 b/src/constitutive_damage.f90 index a6264d57d..86ac47aab 100644 --- a/src/constitutive_damage.f90 +++ b/src/constitutive_damage.f90 @@ -137,24 +137,8 @@ end subroutine damage_init !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the rate of change of microstructure !-------------------------------------------------------------------------------------------------- -module function damage_dotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el,phase,of) result(broken_damage) +module procedure damage_dotState - integer, intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el, & !< element - phase, & - of - real(pReal), intent(in) :: & - subdt !< timestep - real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: & - FArray, & !< elastic deformation gradient - FpArray !< plastic deformation gradient - real(pReal), intent(in), dimension(3,3) :: & - Fi !< intermediate deformation gradient - real(pReal), intent(in), dimension(3,3) :: & - S !< 2nd Piola Kirchhoff stress (vector notation) - logical :: broken_damage integer :: i SourceLoop: do i = 1, phase_Nsources(phase) @@ -176,19 +160,10 @@ module function damage_dotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el,phase broken_damage = any(IEEE_is_NaN(sourceState(phase)%p(i)%dotState(:,of))) -end function damage_dotState +end procedure damage_dotState -module subroutine damage_source_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el) - - integer, intent(in) :: & - ip, & !< integration point number - el !< element number - real(pReal), intent(in) :: & - phi - real(pReal), intent(inout) :: & - phiDot, & - dPhiDot_dPhi +module procedure damage_source_getRateAndItsTangents real(pReal) :: & localphiDot, & @@ -229,6 +204,6 @@ module subroutine damage_source_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, enddo enddo -end subroutine damage_source_getRateAndItsTangents +end procedure damage_source_getRateAndItsTangents end submodule diff --git a/src/constitutive_plastic.f90 b/src/constitutive_plastic.f90 index 5d190758a..f33e7b32f 100644 --- a/src/constitutive_plastic.f90 +++ b/src/constitutive_plastic.f90 @@ -239,32 +239,14 @@ end subroutine plastic_init !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the rate of change of microstructure !-------------------------------------------------------------------------------------------------- -module function plastic_dotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el,phase,of) result(broken_plastic) +module procedure plastic_dotState - integer, intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el, & !< element - phase, & - of - real(pReal), intent(in) :: & - subdt !< timestep - real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: & - FArray, & !< elastic deformation gradient - FpArray !< plastic deformation gradient - real(pReal), intent(in), dimension(3,3) :: & - Fi !< intermediate deformation gradient - real(pReal), intent(in), dimension(3,3) :: & - S !< 2nd Piola Kirchhoff stress (vector notation) real(pReal), dimension(3,3) :: & Mp integer :: & ho, & !< homogenization tme, & !< thermal member position - i, & !< counter in source loop instance - logical :: broken_plastic - ho = material_homogenizationAt(el) tme = thermalMapping(ho)%p(ip,el) @@ -295,20 +277,14 @@ module function plastic_dotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el,phas end select plasticityType broken_plastic = any(IEEE_is_NaN(plasticState(phase)%dotState(:,of))) -end function plastic_dotState +end procedure plastic_dotState !-------------------------------------------------------------------------------------------------- !> @brief returns the homogenize elasticity matrix !> ToDo: homogenizedC66 would be more consistent !-------------------------------------------------------------------------------------------------- -module function plastic_homogenizedC(ipc,ip,el) result(homogenizedC) - - real(pReal), dimension(6,6) :: homogenizedC - integer, intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element +module procedure plastic_homogenizedC plasticityType: select case (phase_plasticity(material_phaseAt(ipc,el))) case (PLASTICITY_DISLOTWIN_ID) plasticityType @@ -317,21 +293,14 @@ module function plastic_homogenizedC(ipc,ip,el) result(homogenizedC) homogenizedC = lattice_C66(1:6,1:6,material_phaseAt(ipc,el)) end select plasticityType -end function plastic_homogenizedC +end procedure plastic_homogenizedC !-------------------------------------------------------------------------------------------------- !> @brief calls microstructure function of the different constitutive models !-------------------------------------------------------------------------------------------------- -module subroutine plastic_dependentState(F, Fp, ipc, ip, el) +module procedure plastic_dependentState - integer, intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element - real(pReal), intent(in), dimension(3,3) :: & - F, & !< elastic deformation gradient - Fp !< plastic deformation gradient integer :: & ho, & !< homogenization tme, & !< thermal member position @@ -351,27 +320,15 @@ module subroutine plastic_dependentState(F, Fp, ipc, ip, el) call plastic_nonlocal_dependentState (F,Fp,instance,of,ip,el) end select plasticityType -end subroutine plastic_dependentState +end procedure plastic_dependentState !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the velocity gradient ! ToDo: Discuss whether it makes sense if crystallite handles the configuration conversion, i.e. ! Mp in, dLp_dMp out !-------------------------------------------------------------------------------------------------- -module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & - S, Fi, ipc, ip, el) - integer, intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element - real(pReal), intent(in), dimension(3,3) :: & - S, & !< 2nd Piola-Kirchhoff stress - Fi !< intermediate deformation gradient - real(pReal), intent(out), dimension(3,3) :: & - Lp !< plastic velocity gradient - real(pReal), intent(out), dimension(3,3,3,3) :: & - dLp_dS, & - dLp_dFi !< derivative of Lp with respect to Fi +module procedure plastic_LpAndItsTangents + real(pReal), dimension(3,3,3,3) :: & dLp_dMp !< derivative of Lp with respect to Mandel stress real(pReal), dimension(3,3) :: & @@ -421,7 +378,7 @@ module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & dLp_dS(i,j,1:3,1:3) = matmul(matmul(transpose(Fi),Fi),dLp_dMp(i,j,1:3,1:3)) ! ToDo: @PS: why not: dLp_dMp:(FiT Fi) enddo; enddo -end subroutine plastic_LpAndItsTangents +end procedure plastic_LpAndItsTangents end submodule constitutive_plastic diff --git a/src/constitutive_thermal.f90 b/src/constitutive_thermal.f90 index 7891ebba4..4a34ac224 100644 --- a/src/constitutive_thermal.f90 +++ b/src/constitutive_thermal.f90 @@ -64,24 +64,8 @@ end subroutine thermal_init !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the rate of change of microstructure !-------------------------------------------------------------------------------------------------- -module function thermal_dotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el,phase,of) result(broken_thermal) +module procedure thermal_dotState - integer, intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el, & !< element - phase, & - of - real(pReal), intent(in) :: & - subdt !< timestep - real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: & - FArray, & !< elastic deformation gradient - FpArray !< plastic deformation gradient - real(pReal), intent(in), dimension(3,3) :: & - Fi !< intermediate deformation gradient - real(pReal), intent(in), dimension(3,3) :: & - S !< 2nd Piola Kirchhoff stress (vector notation) - logical :: broken_thermal integer :: i SourceLoop: do i = 1, phase_Nsources(phase) @@ -93,26 +77,14 @@ module function thermal_dotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el,phas end select sourceType - broken_thermal = any(IEEE_is_NaN(sourceState(phase)%p(i)%dotState(:,of))) - enddo sourceLoop -end function thermal_dotState + broken_thermal = any(IEEE_is_NaN(sourceState(phase)%p(i)%dotState(:,of))) + +end procedure thermal_dotState -module subroutine thermal_source_getRateAndItsTangents(Tdot, dTdot_dT, T, Tstar, Lp, ip, el) - - integer, intent(in) :: & - ip, & !< integration point number - el !< element number - real(pReal), intent(in) :: & - T - real(pReal), intent(in), dimension(:,:,:,:,:) :: & - Tstar, & - Lp - real(pReal), intent(inout) :: & - Tdot, & - dTdot_dT +module procedure thermal_source_getRateAndItsTangents real(pReal) :: & my_Tdot, & @@ -152,6 +124,6 @@ module subroutine thermal_source_getRateAndItsTangents(Tdot, dTdot_dT, T, Tstar, enddo enddo -end subroutine thermal_source_getRateAndItsTangents +end procedure thermal_source_getRateAndItsTangents end submodule From 80fb571fb4fbcbb30d01e76b235c7a363c40ef4a Mon Sep 17 00:00:00 2001 From: Sharan Roongta Date: Fri, 10 Jul 2020 17:10:23 +0200 Subject: [PATCH 05/21] common functions to be clubbed together --- src/constitutive.f90 | 268 +++++++++++++++++++++-------------- src/constitutive_damage.f90 | 54 ------- src/constitutive_plastic.f90 | 112 +-------------- src/constitutive_thermal.f90 | 28 ---- src/crystallite.f90 | 4 +- 5 files changed, 164 insertions(+), 302 deletions(-) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index e1e6ae47d..a8070bd52 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -35,6 +35,103 @@ module constitutive module subroutine thermal_init end subroutine thermal_init + + module subroutine plastic_isotropic_dotState(Mp,instance,of) + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + integer, intent(in) :: & + instance, & + of + end subroutine plastic_isotropic_dotState + + module subroutine plastic_phenopowerlaw_dotState(Mp,instance,of) + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + integer, intent(in) :: & + instance, & + of + end subroutine plastic_phenopowerlaw_dotState + + module subroutine plastic_kinehardening_dotState(Mp,instance,of) + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + integer, intent(in) :: & + instance, & + of + end subroutine plastic_kinehardening_dotState + + module subroutine plastic_dislotwin_dotState(Mp,T,instance,of) + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + real(pReal), intent(in) :: & + T + integer, intent(in) :: & + instance, & + of + end subroutine plastic_dislotwin_dotState + + module subroutine plastic_disloUCLA_dotState(Mp,T,instance,of) + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + real(pReal), intent(in) :: & + T + integer, intent(in) :: & + instance, & + of + end subroutine plastic_disloUCLA_dotState + + module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature,timestep, & + instance,of,ip,el) + real(pReal), dimension(3,3), intent(in) ::& + Mp !< MandelStress + real(pReal), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem), intent(in) :: & + F, & !< deformation gradient + Fp !< plastic deformation gradient + real(pReal), intent(in) :: & + Temperature, & !< temperature + timestep !< substepped crystallite time increment + integer, intent(in) :: & + instance, & + of, & + ip, & !< current integration point + el !< current element number + end subroutine plastic_nonlocal_dotState + + + module subroutine source_damage_anisoBrittle_dotState(S, ipc, ip, el) + + integer, intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + real(pReal), intent(in), dimension(3,3) :: & + S + end subroutine source_damage_anisoBrittle_dotState + + module subroutine source_damage_anisoDuctile_dotState(ipc, ip, el) + + integer, intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + end subroutine source_damage_anisoDuctile_dotState + + module subroutine source_damage_isoDuctile_dotState(ipc, ip, el) + + integer, intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + end subroutine source_damage_isoDuctile_dotState + + module subroutine source_thermal_externalheat_dotState(phase, of) + + integer, intent(in) :: & + phase, & + of + end subroutine source_thermal_externalheat_dotState + + pure module function kinematics_thermal_expansion_initialStrain(homog,phase,offset) result(initialStrain) integer, intent(in) :: & @@ -46,14 +143,14 @@ module constitutive end function kinematics_thermal_expansion_initialStrain - module function plastic_homogenizedC(ipc,ip,el) result(homogenizedC) + module function constitutive_homogenizedC(ipc,ip,el) result(homogenizedC) real(pReal), dimension(6,6) :: & homogenizedC integer, intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point el !< element - end function plastic_homogenizedC + end function constitutive_homogenizedC module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,of) @@ -109,65 +206,6 @@ module constitutive dLi_dTstar !< derivative of Li with respect to Tstar (4th-order tensor defined to be zero) end subroutine kinematics_thermal_expansion_LiAndItsTangent - module function plastic_dotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el,phase,of) result(broken_plastic) - - integer, intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el, & !< element - phase, & - of - real(pReal), intent(in) :: & - subdt !< timestep - real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: & - FArray, & !< elastic deformation gradient - FpArray !< plastic deformation gradient - real(pReal), intent(in), dimension(3,3) :: & - Fi !< intermediate deformation gradient - real(pReal), intent(in), dimension(3,3) :: & - S !< 2nd Piola Kirchhoff stress (vector notation) - logical :: broken_plastic - end function plastic_dotState - - module function damage_dotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el,phase,of) result(broken_damage) - - integer, intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el, & !< element - phase, & - of - real(pReal), intent(in) :: & - subdt !< timestep - real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: & - FArray, & !< elastic deformation gradient - FpArray !< plastic deformation gradient - real(pReal), intent(in), dimension(3,3) :: & - Fi !< intermediate deformation gradient - real(pReal), intent(in), dimension(3,3) :: & - S !< 2nd Piola Kirchhoff stress (vector notation) - logical :: broken_damage - end function damage_dotState - - module function thermal_dotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el,phase,of) result(broken_thermal) - - integer, intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el, & !< element - phase, & - of - real(pReal), intent(in) :: & - subdt !< timestep - real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: & - FArray, & !< elastic deformation gradient - FpArray !< plastic deformation gradient - real(pReal), intent(in), dimension(3,3) :: & - Fi !< intermediate deformation gradient - real(pReal), intent(in), dimension(3,3) :: & - S !< 2nd Piola Kirchhoff stress (vector notation) - logical :: broken_thermal - end function thermal_dotState module subroutine plastic_kinehardening_deltaState(Mp,instance,of) real(pReal), dimension(3,3), intent(in) :: & @@ -261,7 +299,7 @@ module constitutive Fp !< plastic deformation gradient end subroutine plastic_dependentState - module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & + module subroutine constitutive_plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & S, Fi, ipc, ip, el) integer, intent(in) :: & ipc, & !< component-ID of integration point @@ -276,7 +314,7 @@ module constitutive dLp_dS, & dLp_dFi !< derivative of Lp with respect to Fi - end subroutine plastic_LpAndItsTangents + end subroutine constitutive_plastic_LpAndItsTangents module subroutine plastic_nonlocal_updateCompatibility(orientation,instance,i,e) integer, intent(in) :: & @@ -331,7 +369,7 @@ module constitutive public :: & constitutive_init, & constitutive_homogenizedC, & - constitutive_LpAndItsTangents, & + constitutive_plastic_LpAndItsTangents, & constitutive_dependentState, & constitutive_LiAndItsTangents, & constitutive_initialFi, & @@ -392,25 +430,6 @@ subroutine constitutive_init end subroutine constitutive_init - -subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & - S, Fi, ipc, ip, el) - integer, intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element - real(pReal), intent(in), dimension(3,3) :: & - S, & !< 2nd Piola-Kirchhoff stress - Fi !< intermediate deformation gradient - real(pReal), intent(out), dimension(3,3) :: & - Lp !< plastic velocity gradient - real(pReal), intent(out), dimension(3,3,3,3) :: & - dLp_dS, & - dLp_dFi !< derivative of Lp with respect to Fi - - call plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, S, Fi, ipc, ip, el) - -end subroutine constitutive_LpAndItsTangents !-------------------------------------------------------------------------------------------------- !> @brief calls microstructure function of the different constitutive models @@ -430,23 +449,6 @@ subroutine constitutive_dependentState(F, Fp, ipc, ip, el) end subroutine constitutive_dependentState -!-------------------------------------------------------------------------------------------------- -!> @brief returns the homogenize elasticity matrix -!> ToDo: homogenizedC66 would be more consistent -!-------------------------------------------------------------------------------------------------- -function constitutive_homogenizedC(ipc,ip,el) - - real(pReal), dimension(6,6) :: constitutive_homogenizedC - integer, intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element - - constitutive_homogenizedC = plastic_homogenizedC(ipc,ip,el) - -end function constitutive_homogenizedC - - !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the velocity gradient ! ToDo: MD: S is Mi? @@ -654,19 +656,65 @@ function constitutive_collectDotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el Fi !< intermediate deformation gradient real(pReal), intent(in), dimension(3,3) :: & S !< 2nd Piola Kirchhoff stress (vector notation) + real(pReal), dimension(3,3) :: & + Mp + integer :: & + ho, & + tme, & + i, & + instance logical :: broken - logical :: & - broken_plastic, & - broken_thermal, & - broken_damage + ho = material_homogenizationAt(el) + tme = thermalMapping(ho)%p(ip,el) + instance = phase_plasticityInstance(phase) - broken_plastic = plastic_dotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el,phase,of) - broken = broken_plastic - broken_damage = damage_dotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el,phase,of) - broken = broken .or. broken_damage - broken_thermal = thermal_dotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el,phase,of) - broken = broken .or. broken_thermal + Mp = matmul(matmul(transpose(Fi),Fi),S) + + plasticityType: select case (phase_plasticity(phase)) + + case (PLASTICITY_ISOTROPIC_ID) plasticityType + call plastic_isotropic_dotState (Mp,instance,of) + + case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType + call plastic_phenopowerlaw_dotState(Mp,instance,of) + + case (PLASTICITY_KINEHARDENING_ID) plasticityType + call plastic_kinehardening_dotState(Mp,instance,of) + + case (PLASTICITY_DISLOTWIN_ID) plasticityType + call plastic_dislotwin_dotState (Mp,temperature(ho)%p(tme),instance,of) + + case (PLASTICITY_DISLOUCLA_ID) plasticityType + call plastic_disloucla_dotState (Mp,temperature(ho)%p(tme),instance,of) + + case (PLASTICITY_NONLOCAL_ID) plasticityType + call plastic_nonlocal_dotState (Mp,FArray,FpArray,temperature(ho)%p(tme),subdt, & + instance,of,ip,el) + end select plasticityType + broken = any(IEEE_is_NaN(plasticState(phase)%dotState(:,of))) + + SourceLoop: do i = 1, phase_Nsources(phase) + + sourceType: select case (phase_source(i,phase)) + + case (SOURCE_damage_anisoBrittle_ID) sourceType + call source_damage_anisoBrittle_dotState (S, ipc, ip, el) !< correct stress? + + case (SOURCE_damage_isoDuctile_ID) sourceType + call source_damage_isoDuctile_dotState ( ipc, ip, el) + + case (SOURCE_damage_anisoDuctile_ID) sourceType + call source_damage_anisoDuctile_dotState ( ipc, ip, el) + + case (SOURCE_thermal_externalheat_ID) sourceType + call source_thermal_externalheat_dotState(phase,of) + + end select sourceType + + broken = broken .or. any(IEEE_is_NaN(sourceState(phase)%p(i)%dotState(:,of))) + + enddo SourceLoop end function constitutive_collectDotState diff --git a/src/constitutive_damage.f90 b/src/constitutive_damage.f90 index 86ac47aab..b7d1632f8 100644 --- a/src/constitutive_damage.f90 +++ b/src/constitutive_damage.f90 @@ -22,32 +22,6 @@ submodule(constitutive) constitutive_damage module subroutine kinematics_slipplane_opening_init end subroutine kinematics_slipplane_opening_init - module subroutine source_damage_anisoBrittle_dotState(S, ipc, ip, el) - - integer, intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element - real(pReal), intent(in), dimension(3,3) :: & - S - end subroutine source_damage_anisoBrittle_dotState - - module subroutine source_damage_anisoDuctile_dotState(ipc, ip, el) - - integer, intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element - end subroutine source_damage_anisoDuctile_dotState - - module subroutine source_damage_isoDuctile_dotState(ipc, ip, el) - - integer, intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element - end subroutine source_damage_isoDuctile_dotState - module subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) @@ -134,34 +108,6 @@ module subroutine damage_init end subroutine damage_init -!-------------------------------------------------------------------------------------------------- -!> @brief contains the constitutive equation for calculating the rate of change of microstructure -!-------------------------------------------------------------------------------------------------- -module procedure damage_dotState - - integer :: i - - SourceLoop: do i = 1, phase_Nsources(phase) - - sourceType: select case (phase_source(i,phase)) - - case (SOURCE_damage_anisoBrittle_ID) sourceType - call source_damage_anisoBrittle_dotState (S, ipc, ip, el) !< correct stress? - - case (SOURCE_damage_isoDuctile_ID) sourceType - call source_damage_isoDuctile_dotState ( ipc, ip, el) - - case (SOURCE_damage_anisoDuctile_ID) sourceType - call source_damage_anisoDuctile_dotState ( ipc, ip, el) - - end select sourceType - - enddo sourceLoop - - broken_damage = any(IEEE_is_NaN(sourceState(phase)%p(i)%dotState(:,of))) - -end procedure damage_dotState - module procedure damage_source_getRateAndItsTangents diff --git a/src/constitutive_plastic.f90 b/src/constitutive_plastic.f90 index f33e7b32f..c91db80ec 100644 --- a/src/constitutive_plastic.f90 +++ b/src/constitutive_plastic.f90 @@ -25,67 +25,6 @@ submodule(constitutive) constitutive_plastic module subroutine plastic_nonlocal_init end subroutine plastic_nonlocal_init - module subroutine plastic_isotropic_dotState(Mp,instance,of) - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - integer, intent(in) :: & - instance, & - of - end subroutine plastic_isotropic_dotState - - module subroutine plastic_phenopowerlaw_dotState(Mp,instance,of) - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - integer, intent(in) :: & - instance, & - of - end subroutine plastic_phenopowerlaw_dotState - - module subroutine plastic_kinehardening_dotState(Mp,instance,of) - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - integer, intent(in) :: & - instance, & - of - end subroutine plastic_kinehardening_dotState - - module subroutine plastic_dislotwin_dotState(Mp,T,instance,of) - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - real(pReal), intent(in) :: & - T - integer, intent(in) :: & - instance, & - of - end subroutine plastic_dislotwin_dotState - - module subroutine plastic_disloUCLA_dotState(Mp,T,instance,of) - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - real(pReal), intent(in) :: & - T - integer, intent(in) :: & - instance, & - of - end subroutine plastic_disloUCLA_dotState - - module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature,timestep, & - instance,of,ip,el) - real(pReal), dimension(3,3), intent(in) ::& - Mp !< MandelStress - real(pReal), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem), intent(in) :: & - F, & !< deformation gradient - Fp !< plastic deformation gradient - real(pReal), intent(in) :: & - Temperature, & !< temperature - timestep !< substepped crystallite time increment - integer, intent(in) :: & - instance, & - of, & - ip, & !< current integration point - el !< current element number - end subroutine plastic_nonlocal_dotState - module subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) real(pReal), dimension(3,3), intent(out) :: & @@ -236,55 +175,12 @@ module subroutine plastic_init end subroutine plastic_init -!-------------------------------------------------------------------------------------------------- -!> @brief contains the constitutive equation for calculating the rate of change of microstructure -!-------------------------------------------------------------------------------------------------- -module procedure plastic_dotState - - real(pReal), dimension(3,3) :: & - Mp - integer :: & - ho, & !< homogenization - tme, & !< thermal member position - instance - - ho = material_homogenizationAt(el) - tme = thermalMapping(ho)%p(ip,el) - instance = phase_plasticityInstance(phase) - - Mp = matmul(matmul(transpose(Fi),Fi),S) - - plasticityType: select case (phase_plasticity(phase)) - - case (PLASTICITY_ISOTROPIC_ID) plasticityType - call plastic_isotropic_dotState (Mp,instance,of) - - case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType - call plastic_phenopowerlaw_dotState(Mp,instance,of) - - case (PLASTICITY_KINEHARDENING_ID) plasticityType - call plastic_kinehardening_dotState(Mp,instance,of) - - case (PLASTICITY_DISLOTWIN_ID) plasticityType - call plastic_dislotwin_dotState (Mp,temperature(ho)%p(tme),instance,of) - - case (PLASTICITY_DISLOUCLA_ID) plasticityType - call plastic_disloucla_dotState (Mp,temperature(ho)%p(tme),instance,of) - - case (PLASTICITY_NONLOCAL_ID) plasticityType - call plastic_nonlocal_dotState (Mp,FArray,FpArray,temperature(ho)%p(tme),subdt, & - instance,of,ip,el) - end select plasticityType - broken_plastic = any(IEEE_is_NaN(plasticState(phase)%dotState(:,of))) - -end procedure plastic_dotState - !-------------------------------------------------------------------------------------------------- !> @brief returns the homogenize elasticity matrix !> ToDo: homogenizedC66 would be more consistent !-------------------------------------------------------------------------------------------------- -module procedure plastic_homogenizedC +module procedure constitutive_homogenizedC plasticityType: select case (phase_plasticity(material_phaseAt(ipc,el))) case (PLASTICITY_DISLOTWIN_ID) plasticityType @@ -293,7 +189,7 @@ module procedure plastic_homogenizedC homogenizedC = lattice_C66(1:6,1:6,material_phaseAt(ipc,el)) end select plasticityType -end procedure plastic_homogenizedC +end procedure constitutive_homogenizedC !-------------------------------------------------------------------------------------------------- @@ -327,7 +223,7 @@ end procedure plastic_dependentState ! ToDo: Discuss whether it makes sense if crystallite handles the configuration conversion, i.e. ! Mp in, dLp_dMp out !-------------------------------------------------------------------------------------------------- -module procedure plastic_LpAndItsTangents +module procedure constitutive_plastic_LpAndItsTangents real(pReal), dimension(3,3,3,3) :: & dLp_dMp !< derivative of Lp with respect to Mandel stress @@ -378,7 +274,7 @@ module procedure plastic_LpAndItsTangents dLp_dS(i,j,1:3,1:3) = matmul(matmul(transpose(Fi),Fi),dLp_dMp(i,j,1:3,1:3)) ! ToDo: @PS: why not: dLp_dMp:(FiT Fi) enddo; enddo -end procedure plastic_LpAndItsTangents +end procedure constitutive_plastic_LpAndItsTangents end submodule constitutive_plastic diff --git a/src/constitutive_thermal.f90 b/src/constitutive_thermal.f90 index 4a34ac224..48908b997 100644 --- a/src/constitutive_thermal.f90 +++ b/src/constitutive_thermal.f90 @@ -13,12 +13,6 @@ submodule(constitutive) constitutive_thermal module subroutine kinematics_thermal_expansion_init end subroutine kinematics_thermal_expansion_init - module subroutine source_thermal_externalheat_dotState(phase, of) - integer, intent(in) :: & - phase, & - of - end subroutine source_thermal_externalheat_dotState - module subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDot_dT, Tstar, Lp, phase) @@ -61,28 +55,6 @@ module subroutine thermal_init end subroutine thermal_init -!-------------------------------------------------------------------------------------------------- -!> @brief contains the constitutive equation for calculating the rate of change of microstructure -!-------------------------------------------------------------------------------------------------- -module procedure thermal_dotState - - integer :: i - - SourceLoop: do i = 1, phase_Nsources(phase) - - sourceType: select case (phase_source(i,phase)) - - case (SOURCE_thermal_externalheat_ID) sourceType - call source_thermal_externalheat_dotState(phase,of) - - end select sourceType - - enddo sourceLoop - - broken_thermal = any(IEEE_is_NaN(sourceState(phase)%p(i)%dotState(:,of))) - -end procedure thermal_dotState - module procedure thermal_source_getRateAndItsTangents diff --git a/src/crystallite.f90 b/src/crystallite.f90 index fbce3ab47..eaa01270d 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -554,7 +554,7 @@ subroutine crystallite_stressTangent dLidS = math_mul3333xx3333(dLidFi,dFidS) + dLidS endif - call constitutive_LpAndItsTangents(devNull,dLpdS,dLpdFi, & + call constitutive_plastic_LpAndItsTangents(devNull,dLpdS,dLpdFi, & crystallite_S (1:3,1:3,c,i,e), & crystallite_Fi(1:3,1:3,c,i,e),c,i,e) ! call constitutive law to calculate Lp tangent in lattice configuration dLpdS = math_mul3333xx3333(dLpdFi,dFidS) + dLpdS @@ -915,7 +915,7 @@ function integrateStress(ipc,ip,el,timeFraction) result(broken) call constitutive_SandItsTangents(S, dS_dFe, dS_dFi, & Fe, Fi_new, ipc, ip, el) - call constitutive_LpAndItsTangents(Lp_constitutive, dLp_dS, dLp_dFi, & + call constitutive_plastic_LpAndItsTangents(Lp_constitutive, dLp_dS, dLp_dFi, & S, Fi_new, ipc, ip, el) !* update current residuum and check for convergence of loop From 4145ac90d785fba1e7dd8ea789915d4982b38f55 Mon Sep 17 00:00:00 2001 From: Sharan Roongta Date: Fri, 10 Jul 2020 18:19:07 +0200 Subject: [PATCH 06/21] more cleaning --- src/constitutive.f90 | 46 +++++++++++------------------------- src/constitutive_plastic.f90 | 4 ++-- src/crystallite.f90 | 4 ++-- 3 files changed, 18 insertions(+), 36 deletions(-) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index a8070bd52..c44c07de1 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -132,17 +132,6 @@ module constitutive end subroutine source_thermal_externalheat_dotState - 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 constitutive_homogenizedC(ipc,ip,el) result(homogenizedC) real(pReal), dimension(6,6) :: & homogenizedC @@ -288,7 +277,7 @@ module constitutive character(len=*), intent(in) :: group end subroutine source_damage_isoDuctile_results - module subroutine plastic_dependentState(F, Fp, ipc, ip, el) + module subroutine constitutive_plastic_dependentState(F, Fp, ipc, ip, el) integer, intent(in) :: & ipc, & !< component-ID of integration point @@ -297,7 +286,7 @@ module constitutive real(pReal), intent(in), dimension(3,3) :: & F, & !< elastic deformation gradient Fp !< plastic deformation gradient - end subroutine plastic_dependentState + end subroutine constitutive_plastic_dependentState module subroutine constitutive_plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & S, Fi, ipc, ip, el) @@ -315,6 +304,17 @@ module constitutive dLp_dFi !< derivative of Lp with respect to Fi end subroutine constitutive_plastic_LpAndItsTangents + + 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 subroutine plastic_nonlocal_updateCompatibility(orientation,instance,i,e) integer, intent(in) :: & @@ -370,7 +370,7 @@ module constitutive constitutive_init, & constitutive_homogenizedC, & constitutive_plastic_LpAndItsTangents, & - constitutive_dependentState, & + constitutive_plastic_dependentState, & constitutive_LiAndItsTangents, & constitutive_initialFi, & constitutive_SandItsTangents, & @@ -430,24 +430,6 @@ subroutine constitutive_init end subroutine constitutive_init - -!-------------------------------------------------------------------------------------------------- -!> @brief calls microstructure function of the different constitutive models -!-------------------------------------------------------------------------------------------------- -subroutine constitutive_dependentState(F, Fp, ipc, ip, el) - - integer, intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element - real(pReal), intent(in), dimension(3,3) :: & - F, & !< elastic deformation gradient - Fp !< plastic deformation gradient - - call plastic_dependentState(F,Fp,ipc,ip,el) - -end subroutine constitutive_dependentState - !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the velocity gradient diff --git a/src/constitutive_plastic.f90 b/src/constitutive_plastic.f90 index c91db80ec..bfab42833 100644 --- a/src/constitutive_plastic.f90 +++ b/src/constitutive_plastic.f90 @@ -195,7 +195,7 @@ end procedure constitutive_homogenizedC !-------------------------------------------------------------------------------------------------- !> @brief calls microstructure function of the different constitutive models !-------------------------------------------------------------------------------------------------- -module procedure plastic_dependentState +module procedure constitutive_plastic_dependentState integer :: & ho, & !< homogenization @@ -216,7 +216,7 @@ module procedure plastic_dependentState call plastic_nonlocal_dependentState (F,Fp,instance,of,ip,el) end select plasticityType -end procedure plastic_dependentState +end procedure constitutive_plastic_dependentState !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the velocity gradient diff --git a/src/crystallite.f90 b/src/crystallite.f90 index eaa01270d..2ca9e8d00 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -279,7 +279,7 @@ subroutine crystallite_init do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1),FEsolving_execIP(2) do c = 1,homogenization_Ngrains(material_homogenizationAt(e)) - call constitutive_dependentState(crystallite_partionedF0(1:3,1:3,c,i,e), & + call constitutive_plastic_dependentState(crystallite_partionedF0(1:3,1:3,c,i,e), & crystallite_partionedFp0(1:3,1:3,c,i,e), & c,i,e) ! update dependent state variables to be consistent with basic states enddo @@ -874,7 +874,7 @@ function integrateStress(ipc,ip,el,timeFraction) result(broken) F = crystallite_subF(1:3,1:3,ipc,ip,el) endif - call constitutive_dependentState(crystallite_partionedF(1:3,1:3,ipc,ip,el), & + call constitutive_plastic_dependentState(crystallite_partionedF(1:3,1:3,ipc,ip,el), & crystallite_Fp(1:3,1:3,ipc,ip,el),ipc,ip,el) Lpguess = crystallite_Lp(1:3,1:3,ipc,ip,el) ! take as first guess From 3563bce6cbeb0db53a524bf4c90bf68c16538e27 Mon Sep 17 00:00:00 2001 From: Sharan Roongta Date: Fri, 10 Jul 2020 18:29:36 +0200 Subject: [PATCH 07/21] better --- src/constitutive.f90 | 181 ++++++++++++++--------------------- src/constitutive_damage.f90 | 4 +- src/constitutive_thermal.f90 | 4 +- 3 files changed, 78 insertions(+), 111 deletions(-) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index c44c07de1..de735fe55 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -141,6 +141,80 @@ module constitutive el !< element end function constitutive_homogenizedC + module subroutine constitutive_plastic_dependentState(F, Fp, ipc, ip, el) + + integer, intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + real(pReal), intent(in), dimension(3,3) :: & + F, & !< elastic deformation gradient + Fp !< plastic deformation gradient + end subroutine constitutive_plastic_dependentState + + module subroutine constitutive_plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & + S, Fi, ipc, ip, el) + integer, intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + real(pReal), intent(in), dimension(3,3) :: & + S, & !< 2nd Piola-Kirchhoff stress + Fi !< intermediate deformation gradient + real(pReal), intent(out), dimension(3,3) :: & + Lp !< plastic velocity gradient + real(pReal), intent(out), dimension(3,3,3,3) :: & + dLp_dS, & + dLp_dFi !< derivative of Lp with respect to Fi + + end subroutine constitutive_plastic_LpAndItsTangents + + 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 subroutine plastic_nonlocal_updateCompatibility(orientation,instance,i,e) + integer, intent(in) :: & + instance, & + i, & + e + type(rotation), dimension(1,discretization_nIP,discretization_nElem), intent(in) :: & + orientation !< crystal orientation + end subroutine plastic_nonlocal_updateCompatibility + + module subroutine constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el) + integer, intent(in) :: & + ip, & !< integration point number + el !< element number + real(pReal), intent(in) :: & + phi + real(pReal), intent(inout) :: & + phiDot, & + dPhiDot_dPhi + + end subroutine constitutive_damage_getRateAndItsTangents + + module subroutine constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, Tstar, Lp, ip, el) + integer, intent(in) :: & + ip, & !< integration point number + el !< element number + real(pReal), intent(in) :: & + T + real(pReal), intent(in), dimension(:,:,:,:,:) :: & + Tstar, & + Lp + real(pReal), intent(inout) :: & + TDot, & + dTDot_dT + end subroutine constitutive_thermal_getRateAndItsTangents + module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,of) real(pReal), dimension(3,3), intent(out) :: & @@ -226,7 +300,6 @@ module constitutive end subroutine source_damage_isoBrittle_deltaState - module subroutine plastic_isotropic_results(instance,group) integer, intent(in) :: instance character(len=*), intent(in) :: group @@ -277,79 +350,6 @@ module constitutive character(len=*), intent(in) :: group end subroutine source_damage_isoDuctile_results - module subroutine constitutive_plastic_dependentState(F, Fp, ipc, ip, el) - - integer, intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element - real(pReal), intent(in), dimension(3,3) :: & - F, & !< elastic deformation gradient - Fp !< plastic deformation gradient - end subroutine constitutive_plastic_dependentState - - module subroutine constitutive_plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & - S, Fi, ipc, ip, el) - integer, intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element - real(pReal), intent(in), dimension(3,3) :: & - S, & !< 2nd Piola-Kirchhoff stress - Fi !< intermediate deformation gradient - real(pReal), intent(out), dimension(3,3) :: & - Lp !< plastic velocity gradient - real(pReal), intent(out), dimension(3,3,3,3) :: & - dLp_dS, & - dLp_dFi !< derivative of Lp with respect to Fi - - end subroutine constitutive_plastic_LpAndItsTangents - - 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 subroutine plastic_nonlocal_updateCompatibility(orientation,instance,i,e) - integer, intent(in) :: & - instance, & - i, & - e - type(rotation), dimension(1,discretization_nIP,discretization_nElem), intent(in) :: & - orientation !< crystal orientation - end subroutine plastic_nonlocal_updateCompatibility - - module subroutine damage_source_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el) - integer, intent(in) :: & - ip, & !< integration point number - el !< element number - real(pReal), intent(in) :: & - phi - real(pReal), intent(inout) :: & - phiDot, & - dPhiDot_dPhi - end subroutine damage_source_getRateAndItsTangents - - module subroutine thermal_source_getRateAndItsTangents(TDot, dTDot_dT, T, Tstar, Lp, ip, el) - integer, intent(in) :: & - ip, & !< integration point number - el !< element number - real(pReal), intent(in) :: & - T - real(pReal), intent(in), dimension(:,:,:,:,:) :: & - Tstar, & - Lp - real(pReal), intent(inout) :: & - TDot, & - dTDot_dT - end subroutine thermal_source_getRateAndItsTangents - end interface @@ -779,39 +779,6 @@ function constitutive_deltaState(S, Fe, Fi, ipc, ip, el, phase, of) result(broke end function constitutive_deltaState -subroutine constitutive_thermal_getRateAndItsTangents(Tdot, dTDot_dT, T, Tstar, Lp, ip, el) - - integer, intent(in) :: & - ip, & !< integration point number - el !< element number - real(pReal), intent(in) :: & - T - real(pReal), intent(in), dimension(:,:,:,:,:) :: & - Tstar, & - Lp - real(pReal), intent(inout) :: & - Tdot, & - dTdot_dT - - call thermal_source_getRateAndItsTangents(Tdot, dTdot_dT, T, Tstar, Lp, ip, el) - -end subroutine constitutive_thermal_getRateAndItsTangents - -subroutine constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el) - - integer, intent(in) :: & - ip, & !< integration point number - el !< element number - real(pReal), intent(in) :: & - phi - real(pReal), intent(inout) :: & - phiDot, & - dPhiDot_dPhi - - call damage_source_getRateAndItsTangents(phiDot,dPhiDot_dPhi,phi,ip,el) - -end subroutine constitutive_damage_getRateAndItsTangents - !-------------------------------------------------------------------------------------------------- !> @brief writes constitutive results to HDF5 output file diff --git a/src/constitutive_damage.f90 b/src/constitutive_damage.f90 index b7d1632f8..6dae6e298 100644 --- a/src/constitutive_damage.f90 +++ b/src/constitutive_damage.f90 @@ -109,7 +109,7 @@ module subroutine damage_init end subroutine damage_init -module procedure damage_source_getRateAndItsTangents +module procedure constitutive_damage_getRateAndItsTangents real(pReal) :: & localphiDot, & @@ -150,6 +150,6 @@ module procedure damage_source_getRateAndItsTangents enddo enddo -end procedure damage_source_getRateAndItsTangents +end procedure constitutive_damage_getRateAndItsTangents end submodule diff --git a/src/constitutive_thermal.f90 b/src/constitutive_thermal.f90 index 48908b997..5526fb1d3 100644 --- a/src/constitutive_thermal.f90 +++ b/src/constitutive_thermal.f90 @@ -56,7 +56,7 @@ module subroutine thermal_init end subroutine thermal_init -module procedure thermal_source_getRateAndItsTangents +module procedure constitutive_thermal_getRateAndItsTangents real(pReal) :: & my_Tdot, & @@ -96,6 +96,6 @@ module procedure thermal_source_getRateAndItsTangents enddo enddo -end procedure thermal_source_getRateAndItsTangents +end procedure constitutive_thermal_getRateAndItsTangents end submodule From 77567bd3981094683f4329a5d9338b84ba022151 Mon Sep 17 00:00:00 2001 From: Sharan Roongta Date: Fri, 10 Jul 2020 23:41:56 +0200 Subject: [PATCH 08/21] To circumvent Marc internal compiler error --- src/constitutive.f90 | 29 +++++++++++++++++++++++++---- src/constitutive_plastic.f90 | 24 ------------------------ 2 files changed, 25 insertions(+), 28 deletions(-) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index de735fe55..8ecfb1759 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -131,16 +131,15 @@ module constitutive of end subroutine source_thermal_externalheat_dotState - - module function constitutive_homogenizedC(ipc,ip,el) result(homogenizedC) + module function plastic_dislotwin_homogenizedC(ipc,ip,el) result(homogenizedC) real(pReal), dimension(6,6) :: & homogenizedC integer, intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point el !< element - end function constitutive_homogenizedC - + end function plastic_dislotwin_homogenizedC + module subroutine constitutive_plastic_dependentState(F, Fp, ipc, ip, el) integer, intent(in) :: & @@ -430,6 +429,28 @@ subroutine constitutive_init end subroutine constitutive_init +!-------------------------------------------------------------------------------------------------- +!> @brief returns the homogenize elasticity matrix +!> ToDo: homogenizedC66 would be more consistent +!-------------------------------------------------------------------------------------------------- +function constitutive_homogenizedC(ipc,ip,el) + + real(pReal) , dimension(6,6) :: & + constitutive_homogenizedC + integer, intent(in) :: & + ipc, & + ip, & + el + + plasticityType: select case (phase_plasticity(material_phaseAt(ipc,el))) + case (PLASTICITY_DISLOTWIN_ID) plasticityType + constitutive_homogenizedC = plastic_dislotwin_homogenizedC(ipc,ip,el) + case default plasticityType + constitutive_homogenizedC = lattice_C66(1:6,1:6,material_phaseAt(ipc,el)) + end select plasticityType + +end function constitutive_homogenizedC + !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the velocity gradient diff --git a/src/constitutive_plastic.f90 b/src/constitutive_plastic.f90 index bfab42833..ad7e5a54b 100644 --- a/src/constitutive_plastic.f90 +++ b/src/constitutive_plastic.f90 @@ -113,15 +113,6 @@ submodule(constitutive) constitutive_plastic el !< current element number end subroutine plastic_nonlocal_LpAndItsTangent - module function plastic_dislotwin_homogenizedC(ipc,ip,el) result(homogenizedC) - real(pReal), dimension(6,6) :: & - homogenizedC - integer, intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element - end function plastic_dislotwin_homogenizedC - module subroutine plastic_dislotwin_dependentState(T,instance,of) integer, intent(in) :: & @@ -176,21 +167,6 @@ module subroutine plastic_init end subroutine plastic_init -!-------------------------------------------------------------------------------------------------- -!> @brief returns the homogenize elasticity matrix -!> ToDo: homogenizedC66 would be more consistent -!-------------------------------------------------------------------------------------------------- -module procedure constitutive_homogenizedC - - plasticityType: select case (phase_plasticity(material_phaseAt(ipc,el))) - case (PLASTICITY_DISLOTWIN_ID) plasticityType - homogenizedC = plastic_dislotwin_homogenizedC(ipc,ip,el) - case default plasticityType - homogenizedC = lattice_C66(1:6,1:6,material_phaseAt(ipc,el)) - end select plasticityType - -end procedure constitutive_homogenizedC - !-------------------------------------------------------------------------------------------------- !> @brief calls microstructure function of the different constitutive models From 3a5e3b36c10b0c9da51e5fbc554a13a804e708bd Mon Sep 17 00:00:00 2001 From: Sharan Roongta Date: Sun, 12 Jul 2020 13:04:26 +0200 Subject: [PATCH 09/21] better function name, crystallite should not know which physics is involved --- src/constitutive.f90 | 12 ++++++------ src/constitutive_plastic.f90 | 8 ++++---- src/crystallite.f90 | 8 ++++---- 3 files changed, 14 insertions(+), 14 deletions(-) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 8ecfb1759..62e9c661f 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -140,7 +140,7 @@ module constitutive el !< element end function plastic_dislotwin_homogenizedC - module subroutine constitutive_plastic_dependentState(F, Fp, ipc, ip, el) + module subroutine constitutive_dependentState(F, Fp, ipc, ip, el) integer, intent(in) :: & ipc, & !< component-ID of integration point @@ -149,9 +149,9 @@ module constitutive real(pReal), intent(in), dimension(3,3) :: & F, & !< elastic deformation gradient Fp !< plastic deformation gradient - end subroutine constitutive_plastic_dependentState + end subroutine constitutive_dependentState - module subroutine constitutive_plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & + module subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & S, Fi, ipc, ip, el) integer, intent(in) :: & ipc, & !< component-ID of integration point @@ -166,7 +166,7 @@ module constitutive dLp_dS, & dLp_dFi !< derivative of Lp with respect to Fi - end subroutine constitutive_plastic_LpAndItsTangents + end subroutine constitutive_LpAndItsTangents pure module function kinematics_thermal_expansion_initialStrain(homog,phase,offset) result(initialStrain) @@ -368,8 +368,8 @@ module constitutive public :: & constitutive_init, & constitutive_homogenizedC, & - constitutive_plastic_LpAndItsTangents, & - constitutive_plastic_dependentState, & + constitutive_LpAndItsTangents, & + constitutive_dependentState, & constitutive_LiAndItsTangents, & constitutive_initialFi, & constitutive_SandItsTangents, & diff --git a/src/constitutive_plastic.f90 b/src/constitutive_plastic.f90 index ad7e5a54b..7ede355e1 100644 --- a/src/constitutive_plastic.f90 +++ b/src/constitutive_plastic.f90 @@ -171,7 +171,7 @@ end subroutine plastic_init !-------------------------------------------------------------------------------------------------- !> @brief calls microstructure function of the different constitutive models !-------------------------------------------------------------------------------------------------- -module procedure constitutive_plastic_dependentState +module procedure constitutive_dependentState integer :: & ho, & !< homogenization @@ -192,14 +192,14 @@ module procedure constitutive_plastic_dependentState call plastic_nonlocal_dependentState (F,Fp,instance,of,ip,el) end select plasticityType -end procedure constitutive_plastic_dependentState +end procedure constitutive_dependentState !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the velocity gradient ! ToDo: Discuss whether it makes sense if crystallite handles the configuration conversion, i.e. ! Mp in, dLp_dMp out !-------------------------------------------------------------------------------------------------- -module procedure constitutive_plastic_LpAndItsTangents +module procedure constitutive_LpAndItsTangents real(pReal), dimension(3,3,3,3) :: & dLp_dMp !< derivative of Lp with respect to Mandel stress @@ -250,7 +250,7 @@ module procedure constitutive_plastic_LpAndItsTangents dLp_dS(i,j,1:3,1:3) = matmul(matmul(transpose(Fi),Fi),dLp_dMp(i,j,1:3,1:3)) ! ToDo: @PS: why not: dLp_dMp:(FiT Fi) enddo; enddo -end procedure constitutive_plastic_LpAndItsTangents +end procedure constitutive_LpAndItsTangents end submodule constitutive_plastic diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 2ca9e8d00..fbce3ab47 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -279,7 +279,7 @@ subroutine crystallite_init do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1),FEsolving_execIP(2) do c = 1,homogenization_Ngrains(material_homogenizationAt(e)) - call constitutive_plastic_dependentState(crystallite_partionedF0(1:3,1:3,c,i,e), & + call constitutive_dependentState(crystallite_partionedF0(1:3,1:3,c,i,e), & crystallite_partionedFp0(1:3,1:3,c,i,e), & c,i,e) ! update dependent state variables to be consistent with basic states enddo @@ -554,7 +554,7 @@ subroutine crystallite_stressTangent dLidS = math_mul3333xx3333(dLidFi,dFidS) + dLidS endif - call constitutive_plastic_LpAndItsTangents(devNull,dLpdS,dLpdFi, & + call constitutive_LpAndItsTangents(devNull,dLpdS,dLpdFi, & crystallite_S (1:3,1:3,c,i,e), & crystallite_Fi(1:3,1:3,c,i,e),c,i,e) ! call constitutive law to calculate Lp tangent in lattice configuration dLpdS = math_mul3333xx3333(dLpdFi,dFidS) + dLpdS @@ -874,7 +874,7 @@ function integrateStress(ipc,ip,el,timeFraction) result(broken) F = crystallite_subF(1:3,1:3,ipc,ip,el) endif - call constitutive_plastic_dependentState(crystallite_partionedF(1:3,1:3,ipc,ip,el), & + call constitutive_dependentState(crystallite_partionedF(1:3,1:3,ipc,ip,el), & crystallite_Fp(1:3,1:3,ipc,ip,el),ipc,ip,el) Lpguess = crystallite_Lp(1:3,1:3,ipc,ip,el) ! take as first guess @@ -915,7 +915,7 @@ function integrateStress(ipc,ip,el,timeFraction) result(broken) call constitutive_SandItsTangents(S, dS_dFe, dS_dFi, & Fe, Fi_new, ipc, ip, el) - call constitutive_plastic_LpAndItsTangents(Lp_constitutive, dLp_dS, dLp_dFi, & + call constitutive_LpAndItsTangents(Lp_constitutive, dLp_dS, dLp_dFi, & S, Fi_new, ipc, ip, el) !* update current residuum and check for convergence of loop From 5602abe690b48c76b60c3c0d8d9fa841519e924c Mon Sep 17 00:00:00 2001 From: Sharan Roongta Date: Sun, 12 Jul 2020 13:27:28 +0200 Subject: [PATCH 10/21] generic interfaces makes sense --- src/constitutive.f90 | 23 +++++++++++++++++------ src/constitutive_plastic.f90 | 8 ++++---- src/damage_local.f90 | 2 +- src/damage_nonlocal.f90 | 2 +- src/thermal_adiabatic.f90 | 2 +- src/thermal_conduction.f90 | 2 +- 6 files changed, 25 insertions(+), 14 deletions(-) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 62e9c661f..cb970957f 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -140,7 +140,7 @@ module constitutive el !< element end function plastic_dislotwin_homogenizedC - module subroutine constitutive_dependentState(F, Fp, ipc, ip, el) + module subroutine constitutive_plastic_dependentState(F, Fp, ipc, ip, el) integer, intent(in) :: & ipc, & !< component-ID of integration point @@ -149,9 +149,9 @@ module constitutive real(pReal), intent(in), dimension(3,3) :: & F, & !< elastic deformation gradient Fp !< plastic deformation gradient - end subroutine constitutive_dependentState + end subroutine constitutive_plastic_dependentState - module subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & + module subroutine constitutive_plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & S, Fi, ipc, ip, el) integer, intent(in) :: & ipc, & !< component-ID of integration point @@ -166,7 +166,7 @@ module constitutive dLp_dS, & dLp_dFi !< derivative of Lp with respect to Fi - end subroutine constitutive_LpAndItsTangents + end subroutine constitutive_plastic_LpAndItsTangents pure module function kinematics_thermal_expansion_initialStrain(homog,phase,offset) result(initialStrain) @@ -351,6 +351,18 @@ module constitutive end interface + interface constitutive_LpAndItsTangents + module procedure :: constitutive_plastic_LpAndItsTangents + end interface constitutive_LpAndItsTangents + + interface constitutive_dependentState + module procedure :: constitutive_plastic_dependentState + end interface constitutive_dependentState + + interface constitutive_getRateAndItsTangents + module procedure :: constitutive_damage_getRateAndItsTangents , & + constitutive_thermal_getRateAndItsTangents + end interface constitutive_getRateAndItsTangents type :: tDebugOptions logical :: & @@ -376,8 +388,7 @@ module constitutive constitutive_collectDotState, & constitutive_deltaState, & plastic_nonlocal_updateCompatibility, & - constitutive_damage_getRateAndItsTangents, & - constitutive_thermal_getRateAndItsTangents, & + constitutive_getRateAndItsTangents, & constitutive_results contains diff --git a/src/constitutive_plastic.f90 b/src/constitutive_plastic.f90 index 7ede355e1..ad7e5a54b 100644 --- a/src/constitutive_plastic.f90 +++ b/src/constitutive_plastic.f90 @@ -171,7 +171,7 @@ end subroutine plastic_init !-------------------------------------------------------------------------------------------------- !> @brief calls microstructure function of the different constitutive models !-------------------------------------------------------------------------------------------------- -module procedure constitutive_dependentState +module procedure constitutive_plastic_dependentState integer :: & ho, & !< homogenization @@ -192,14 +192,14 @@ module procedure constitutive_dependentState call plastic_nonlocal_dependentState (F,Fp,instance,of,ip,el) end select plasticityType -end procedure constitutive_dependentState +end procedure constitutive_plastic_dependentState !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the velocity gradient ! ToDo: Discuss whether it makes sense if crystallite handles the configuration conversion, i.e. ! Mp in, dLp_dMp out !-------------------------------------------------------------------------------------------------- -module procedure constitutive_LpAndItsTangents +module procedure constitutive_plastic_LpAndItsTangents real(pReal), dimension(3,3,3,3) :: & dLp_dMp !< derivative of Lp with respect to Mandel stress @@ -250,7 +250,7 @@ module procedure constitutive_LpAndItsTangents dLp_dS(i,j,1:3,1:3) = matmul(matmul(transpose(Fi),Fi),dLp_dMp(i,j,1:3,1:3)) ! ToDo: @PS: why not: dLp_dMp:(FiT Fi) enddo; enddo -end procedure constitutive_LpAndItsTangents +end procedure constitutive_plastic_LpAndItsTangents end submodule constitutive_plastic diff --git a/src/damage_local.f90 b/src/damage_local.f90 index 66b50064b..3d448f7a4 100644 --- a/src/damage_local.f90 +++ b/src/damage_local.f90 @@ -131,7 +131,7 @@ subroutine damage_local_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el phiDot = 0.0_pReal dPhiDot_dPhi = 0.0_pReal - call constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el) + call constitutive_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el) phiDot = phiDot/real(homogenization_Ngrains(material_homogenizationAt(el)),pReal) dPhiDot_dPhi = dPhiDot_dPhi/real(homogenization_Ngrains(material_homogenizationAt(el)),pReal) diff --git a/src/damage_nonlocal.f90 b/src/damage_nonlocal.f90 index 914133de7..530a79899 100644 --- a/src/damage_nonlocal.f90 +++ b/src/damage_nonlocal.f90 @@ -100,7 +100,7 @@ subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, phiDot = 0.0_pReal dPhiDot_dPhi = 0.0_pReal - call constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el) + call constitutive_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el) phiDot = phiDot/real(homogenization_Ngrains(material_homogenizationAt(el)),pReal) dPhiDot_dPhi = dPhiDot_dPhi/real(homogenization_Ngrains(material_homogenizationAt(el)),pReal) diff --git a/src/thermal_adiabatic.f90 b/src/thermal_adiabatic.f90 index c52f0a3d0..9e4710a1e 100644 --- a/src/thermal_adiabatic.f90 +++ b/src/thermal_adiabatic.f90 @@ -131,7 +131,7 @@ subroutine thermal_adiabatic_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el) dTdot_dT = 0.0_pReal homog = material_homogenizationAt(el) - call constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, crystallite_S, crystallite_Lp, ip, el) + call constitutive_getRateAndItsTangents(TDot, dTDot_dT, T, crystallite_S, crystallite_Lp, ip, el) Tdot = Tdot/real(homogenization_Ngrains(homog),pReal) dTdot_dT = dTdot_dT/real(homogenization_Ngrains(homog),pReal) diff --git a/src/thermal_conduction.f90 b/src/thermal_conduction.f90 index 60766710d..2afe411f1 100644 --- a/src/thermal_conduction.f90 +++ b/src/thermal_conduction.f90 @@ -90,7 +90,7 @@ subroutine thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el) dTdot_dT = 0.0_pReal homog = material_homogenizationAt(el) - call constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, crystallite_S,crystallite_Lp ,ip, el) + call constitutive_getRateAndItsTangents(TDot, dTDot_dT, T, crystallite_S,crystallite_Lp ,ip, el) Tdot = Tdot/real(homogenization_Ngrains(homog),pReal) dTdot_dT = dTdot_dT/real(homogenization_Ngrains(homog),pReal) From debe096a538acf10f5d38ff718d07e78311a8888 Mon Sep 17 00:00:00 2001 From: Sharan Roongta Date: Sun, 12 Jul 2020 15:22:40 +0200 Subject: [PATCH 11/21] results placed where it belongs; cleaning --- src/constitutive.f90 | 109 +++-------------------------------- src/constitutive_damage.f90 | 70 ++++++++++++++++++---- src/constitutive_plastic.f90 | 72 ++++++++++++++++++++++- src/constitutive_thermal.f90 | 2 - 4 files changed, 135 insertions(+), 118 deletions(-) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index cb970957f..040858e3f 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -298,56 +298,11 @@ module constitutive C end subroutine source_damage_isoBrittle_deltaState - - module subroutine plastic_isotropic_results(instance,group) - integer, intent(in) :: instance - character(len=*), intent(in) :: group - end subroutine plastic_isotropic_results - - module subroutine plastic_phenopowerlaw_results(instance,group) - integer, intent(in) :: instance - character(len=*), intent(in) :: group - end subroutine plastic_phenopowerlaw_results - - module subroutine plastic_kinehardening_results(instance,group) - integer, intent(in) :: instance - character(len=*), intent(in) :: group - end subroutine plastic_kinehardening_results - - module subroutine plastic_dislotwin_results(instance,group) - integer, intent(in) :: instance - character(len=*), intent(in) :: group - end subroutine plastic_dislotwin_results - - module subroutine plastic_disloUCLA_results(instance,group) - integer, intent(in) :: instance - character(len=*), intent(in) :: group - end subroutine plastic_disloUCLA_results - - module subroutine plastic_nonlocal_results(instance,group) - integer, intent(in) :: instance - character(len=*), intent(in) :: group - end subroutine plastic_nonlocal_results - - module subroutine source_damage_anisoBrittle_results(phase,group) - integer, intent(in) :: phase - character(len=*), intent(in) :: group - end subroutine source_damage_anisoBrittle_results - - module subroutine source_damage_anisoDuctile_results(phase,group) - integer, intent(in) :: phase - character(len=*), intent(in) :: group - end subroutine source_damage_anisoDuctile_results - - module subroutine source_damage_isoBrittle_results(phase,group) - integer, intent(in) :: phase - character(len=*), intent(in) :: group - end subroutine source_damage_isoBrittle_results - - module subroutine source_damage_isoDuctile_results(phase,group) - integer, intent(in) :: phase - character(len=*), intent(in) :: group - end subroutine source_damage_isoDuctile_results + module subroutine plastic_results + end subroutine plastic_results + + module subroutine damage_results + end subroutine damage_results end interface @@ -817,58 +772,8 @@ end function constitutive_deltaState !-------------------------------------------------------------------------------------------------- subroutine constitutive_results - integer :: p,i - character(len=pStringLen) :: group,group_plastic,group_sources - - plasticityLoop: do p=1,size(config_name_phase) - group = trim('current/constituent')//'/'//trim(config_name_phase(p)) - call results_closeGroup(results_addGroup(group)) - - group_plastic = trim(group)//'/plastic' - - call results_closeGroup(results_addGroup(group_plastic)) - select case(phase_plasticity(p)) - - case(PLASTICITY_ISOTROPIC_ID) - call plastic_isotropic_results(phase_plasticityInstance(p),group_plastic) - - case(PLASTICITY_PHENOPOWERLAW_ID) - call plastic_phenopowerlaw_results(phase_plasticityInstance(p),group_plastic) - - case(PLASTICITY_KINEHARDENING_ID) - call plastic_kinehardening_results(phase_plasticityInstance(p),group_plastic) - - case(PLASTICITY_DISLOTWIN_ID) - call plastic_dislotwin_results(phase_plasticityInstance(p),group_plastic) - - case(PLASTICITY_DISLOUCLA_ID) - call plastic_disloUCLA_results(phase_plasticityInstance(p),group_plastic) - - case(PLASTICITY_NONLOCAL_ID) - call plastic_nonlocal_results(phase_plasticityInstance(p),group_plastic) - end select - - sourceLoop: do i = 1, phase_Nsources(p) - group_sources = trim(group)//'/sources' - - call results_closeGroup(results_addGroup(group_sources)) - - sourceType: select case (phase_source(i,p)) - - case (SOURCE_damage_anisoBrittle_ID) sourceType - call source_damage_anisoBrittle_results(p,group_sources) - case (SOURCE_damage_anisoDuctile_ID) sourceType - call source_damage_anisoDuctile_results(p,group_sources) - case (SOURCE_damage_isoBrittle_ID) sourceType - call source_damage_isoBrittle_results(p,group_sources) - case (SOURCE_damage_isoDuctile_ID) sourceType - call source_damage_isoDuctile_results(p,group_sources) - end select sourceType - - enddo SourceLoop - - enddo plasticityLoop - + call plastic_results + call damage_results end subroutine constitutive_results diff --git a/src/constitutive_damage.f90 b/src/constitutive_damage.f90 index 6dae6e298..fe9eb4e41 100644 --- a/src/constitutive_damage.f90 +++ b/src/constitutive_damage.f90 @@ -1,7 +1,5 @@ submodule(constitutive) constitutive_damage - implicit none - interface module subroutine source_damage_anisoBrittle_init @@ -24,7 +22,6 @@ submodule(constitutive) constitutive_damage module subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) - integer, intent(in) :: & phase, & constituent @@ -33,11 +30,9 @@ submodule(constitutive) constitutive_damage real(pReal), intent(out) :: & localphiDot, & dLocalphiDot_dPhi - end subroutine source_damage_anisoBrittle_getRateAndItsTangent module subroutine source_damage_anisoDuctile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) - integer, intent(in) :: & phase, & constituent @@ -46,11 +41,9 @@ submodule(constitutive) constitutive_damage real(pReal), intent(out) :: & localphiDot, & dLocalphiDot_dPhi - end subroutine source_damage_anisoDuctile_getRateAndItsTangent module subroutine source_damage_isoBrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) - integer, intent(in) :: & phase, & constituent @@ -59,11 +52,9 @@ submodule(constitutive) constitutive_damage real(pReal), intent(out) :: & localphiDot, & dLocalphiDot_dPhi - end subroutine source_damage_isoBrittle_getRateAndItsTangent module subroutine source_damage_isoDuctile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) - integer, intent(in) :: & phase, & constituent @@ -72,11 +63,9 @@ submodule(constitutive) constitutive_damage real(pReal), intent(out) :: & localphiDot, & dLocalphiDot_dPhi - end subroutine source_damage_isoDuctile_getRateAndItsTangent module subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDot_dT, Tstar, Lp, phase) - integer, intent(in) :: & phase real(pReal), intent(in), dimension(3,3) :: & @@ -87,8 +76,28 @@ submodule(constitutive) constitutive_damage real(pReal), intent(out) :: & TDot, & dTDot_dT - end subroutine source_thermal_dissipation_getRateAndItsTangent + + module subroutine source_damage_anisoBrittle_results(phase,group) + integer, intent(in) :: phase + character(len=*), intent(in) :: group + end subroutine source_damage_anisoBrittle_results + + module subroutine source_damage_anisoDuctile_results(phase,group) + integer, intent(in) :: phase + character(len=*), intent(in) :: group + end subroutine source_damage_anisoDuctile_results + + module subroutine source_damage_isoBrittle_results(phase,group) + integer, intent(in) :: phase + character(len=*), intent(in) :: group + end subroutine source_damage_isoBrittle_results + + module subroutine source_damage_isoDuctile_results(phase,group) + integer, intent(in) :: phase + character(len=*), intent(in) :: group + end subroutine source_damage_isoDuctile_results + end interface contains @@ -152,4 +161,41 @@ module procedure constitutive_damage_getRateAndItsTangents end procedure constitutive_damage_getRateAndItsTangents + +module subroutine damage_results + + integer :: p,i + character(len=pStringLen) :: group + + do p = 1, size(config_name_phase) + sourceLoop: do i = 1, phase_Nsources(p) + group = trim('current/constituent')//'/'//trim(config_name_phase(p)) + group = trim(group)//'/sources' + + sourceType: select case (phase_source(i,p)) + + case (SOURCE_damage_anisoBrittle_ID) sourceType + group = trim(group)//'/damage_anisoBrittle' + call results_closeGroup(results_addGroup(group)) + call source_damage_anisoBrittle_results(p,group) + case (SOURCE_damage_anisoDuctile_ID) sourceType + group = trim(group)//'/damage_anisoDuctile' + call results_closeGroup(results_addGroup(group)) + call source_damage_anisoDuctile_results(p,group) + case (SOURCE_damage_isoBrittle_ID) sourceType + group = trim(group)//'/damage_isoBrittle' + call results_closeGroup(results_addGroup(group)) + call source_damage_isoBrittle_results(p,group) + case (SOURCE_damage_isoDuctile_ID) sourceType + group = trim(group)//'/damage_isoDuctile' + call results_closeGroup(results_addGroup(group)) + call source_damage_isoDuctile_results(p,group) + end select sourceType + + enddo SourceLoop + enddo + +end subroutine damage_results + + end submodule diff --git a/src/constitutive_plastic.f90 b/src/constitutive_plastic.f90 index ad7e5a54b..873887149 100644 --- a/src/constitutive_plastic.f90 +++ b/src/constitutive_plastic.f90 @@ -1,7 +1,5 @@ submodule(constitutive) constitutive_plastic - implicit none - interface module subroutine plastic_none_init @@ -139,6 +137,38 @@ submodule(constitutive) constitutive_plastic el end subroutine plastic_nonlocal_dependentState + module subroutine plastic_isotropic_results(instance,group) + integer, intent(in) :: instance + character(len=*), intent(in) :: group + end subroutine plastic_isotropic_results + + module subroutine plastic_phenopowerlaw_results(instance,group) + integer, intent(in) :: instance + character(len=*), intent(in) :: group + end subroutine plastic_phenopowerlaw_results + + module subroutine plastic_kinehardening_results(instance,group) + integer, intent(in) :: instance + character(len=*), intent(in) :: group + end subroutine plastic_kinehardening_results + + module subroutine plastic_dislotwin_results(instance,group) + integer, intent(in) :: instance + character(len=*), intent(in) :: group + end subroutine plastic_dislotwin_results + + module subroutine plastic_disloUCLA_results(instance,group) + integer, intent(in) :: instance + character(len=*), intent(in) :: group + end subroutine plastic_disloUCLA_results + + module subroutine plastic_nonlocal_results(instance,group) + integer, intent(in) :: instance + character(len=*), intent(in) :: group + end subroutine plastic_nonlocal_results + + + end interface @@ -252,5 +282,43 @@ module procedure constitutive_plastic_LpAndItsTangents end procedure constitutive_plastic_LpAndItsTangents +module subroutine plastic_results + + integer :: p + character(len=pStringLen) :: group + + plasticityLoop: do p=1,size(config_name_phase) + group = trim('current/constituent')//'/'//trim(config_name_phase(p)) + call results_closeGroup(results_addGroup(group)) + + group = trim(group)//'/plastic' + + call results_closeGroup(results_addGroup(group)) + select case(phase_plasticity(p)) + + case(PLASTICITY_ISOTROPIC_ID) + call plastic_isotropic_results(phase_plasticityInstance(p),group) + + case(PLASTICITY_PHENOPOWERLAW_ID) + call plastic_phenopowerlaw_results(phase_plasticityInstance(p),group) + + case(PLASTICITY_KINEHARDENING_ID) + call plastic_kinehardening_results(phase_plasticityInstance(p),group) + + case(PLASTICITY_DISLOTWIN_ID) + call plastic_dislotwin_results(phase_plasticityInstance(p),group) + + case(PLASTICITY_DISLOUCLA_ID) + call plastic_disloUCLA_results(phase_plasticityInstance(p),group) + + case(PLASTICITY_NONLOCAL_ID) + call plastic_nonlocal_results(phase_plasticityInstance(p),group) + end select + + enddo plasticityLoop + +end subroutine plastic_results + + end submodule constitutive_plastic diff --git a/src/constitutive_thermal.f90 b/src/constitutive_thermal.f90 index 5526fb1d3..ca40e21fa 100644 --- a/src/constitutive_thermal.f90 +++ b/src/constitutive_thermal.f90 @@ -1,7 +1,5 @@ submodule(constitutive) constitutive_thermal - implicit none - interface module subroutine source_thermal_dissipation_init From 70fb68d22490397be147df753ef169f038b160a9 Mon Sep 17 00:00:00 2001 From: Sharan Roongta Date: Sun, 12 Jul 2020 16:44:26 +0200 Subject: [PATCH 12/21] cleaning --- src/constitutive.f90 | 90 ++++++++++++++++++------------------ src/constitutive_damage.f90 | 28 +++++++++-- src/constitutive_plastic.f90 | 40 ++++++++++++---- src/constitutive_thermal.f90 | 22 ++++++++- 4 files changed, 120 insertions(+), 60 deletions(-) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 040858e3f..abff6e21d 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -101,9 +101,9 @@ module constitutive module subroutine source_damage_anisoBrittle_dotState(S, ipc, ip, el) integer, intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element real(pReal), intent(in), dimension(3,3) :: & S end subroutine source_damage_anisoBrittle_dotState @@ -111,17 +111,17 @@ module constitutive module subroutine source_damage_anisoDuctile_dotState(ipc, ip, el) integer, intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element end subroutine source_damage_anisoDuctile_dotState module subroutine source_damage_isoDuctile_dotState(ipc, ip, el) integer, intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element end subroutine source_damage_isoDuctile_dotState module subroutine source_thermal_externalheat_dotState(phase, of) @@ -143,28 +143,28 @@ module constitutive module subroutine constitutive_plastic_dependentState(F, Fp, ipc, ip, el) integer, intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element real(pReal), intent(in), dimension(3,3) :: & - F, & !< elastic deformation gradient - Fp !< plastic deformation gradient + F, & !< elastic deformation gradient + Fp !< plastic deformation gradient end subroutine constitutive_plastic_dependentState module subroutine constitutive_plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & S, Fi, ipc, ip, el) integer, intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element real(pReal), intent(in), dimension(3,3) :: & - S, & !< 2nd Piola-Kirchhoff stress - Fi !< intermediate deformation gradient + S, & !< 2nd Piola-Kirchhoff stress + Fi !< intermediate deformation gradient real(pReal), intent(out), dimension(3,3) :: & - Lp !< plastic velocity gradient + Lp !< plastic velocity gradient real(pReal), intent(out), dimension(3,3,3,3) :: & dLp_dS, & - dLp_dFi !< derivative of Lp with respect to Fi + dLp_dFi !< derivative of Lp with respect to Fi end subroutine constitutive_plastic_LpAndItsTangents @@ -190,8 +190,8 @@ module constitutive module subroutine constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el) integer, intent(in) :: & - ip, & !< integration point number - el !< element number + ip, & !< integration point number + el !< element number real(pReal), intent(in) :: & phi real(pReal), intent(inout) :: & @@ -202,8 +202,8 @@ module constitutive module subroutine constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, Tstar, Lp, ip, el) integer, intent(in) :: & - ip, & !< integration point number - el !< element number + ip, & !< integration point number + el !< element number real(pReal), intent(in) :: & T real(pReal), intent(in), dimension(:,:,:,:,:) :: & @@ -231,41 +231,41 @@ module constitutive module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, ip, el) integer, intent(in) :: & - ipc, & !< grain number - ip, & !< integration point number - el !< element number + ipc, & !< grain number + ip, & !< integration point number + el !< element number real(pReal), intent(in), dimension(3,3) :: & S real(pReal), intent(out), dimension(3,3) :: & - Ld !< damage velocity gradient + 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) + dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor) end subroutine kinematics_cleavage_opening_LiAndItsTangent module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, ip, el) integer, intent(in) :: & - ipc, & !< grain number - ip, & !< integration point number - el !< element number + ipc, & !< grain number + ip, & !< integration point number + el !< element number real(pReal), intent(in), dimension(3,3) :: & S real(pReal), intent(out), dimension(3,3) :: & - Ld !< damage velocity gradient + 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) + dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor) end subroutine kinematics_slipplane_opening_LiAndItsTangent module subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, ipc, ip, el) integer, intent(in) :: & - ipc, & !< grain number - ip, & !< integration point number - el !< element number + ipc, & !< grain number + ip, & !< integration point number + el !< element number real(pReal), intent(out), dimension(3,3) :: & - Li !< thermal velocity gradient + Li !< thermal velocity gradient real(pReal), intent(out), dimension(3,3,3,3) :: & - dLi_dTstar !< derivative of Li with respect to Tstar (4th-order tensor defined to be zero) + dLi_dTstar !< derivative of Li with respect to Tstar (4th-order tensor defined to be zero) end subroutine kinematics_thermal_expansion_LiAndItsTangent @@ -289,9 +289,9 @@ module constitutive module subroutine source_damage_isoBrittle_deltaState(C, Fe, ipc, ip, el) integer, intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element real(pReal), intent(in), dimension(3,3) :: & Fe real(pReal), intent(in), dimension(6,6) :: & @@ -613,7 +613,7 @@ function constitutive_collectDotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el integer, intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point - el, & !< element + el, & !< element phase, & of real(pReal), intent(in) :: & @@ -668,7 +668,7 @@ function constitutive_collectDotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el sourceType: select case (phase_source(i,phase)) case (SOURCE_damage_anisoBrittle_ID) sourceType - call source_damage_anisoBrittle_dotState (S, ipc, ip, el) !< correct stress? + call source_damage_anisoBrittle_dotState (S, ipc, ip, el) ! correct stress? case (SOURCE_damage_isoDuctile_ID) sourceType call source_damage_isoDuctile_dotState ( ipc, ip, el) diff --git a/src/constitutive_damage.f90 b/src/constitutive_damage.f90 index fe9eb4e41..5e1ed848d 100644 --- a/src/constitutive_damage.f90 +++ b/src/constitutive_damage.f90 @@ -23,10 +23,10 @@ submodule(constitutive) constitutive_damage module subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) integer, intent(in) :: & - phase, & - constituent + phase, & !< phase ID of element + constituent !< position of element within its phase instance real(pReal), intent(in) :: & - phi + phi !< damage value real(pReal), intent(out) :: & localphiDot, & dLocalphiDot_dPhi @@ -102,6 +102,9 @@ submodule(constitutive) constitutive_damage contains +!---------------------------------------------------------------------------------------------- +!< @brief initialize damage sources and kinematics mechanism +!---------------------------------------------------------------------------------------------- module subroutine damage_init ! initialize source mechanisms @@ -118,7 +121,19 @@ module subroutine damage_init end subroutine damage_init -module procedure constitutive_damage_getRateAndItsTangents +!---------------------------------------------------------------------------------------------- +!< @brief returns local part of nonlocal damage driving force +!---------------------------------------------------------------------------------------------- +module subroutine constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el) + + integer, intent(in) :: & + ip, & !< integration point number + el !< element number + real(pReal), intent(in) :: & + phi + real(pReal), intent(inout) :: & + phiDot, & + dPhiDot_dPhi real(pReal) :: & localphiDot, & @@ -159,9 +174,12 @@ module procedure constitutive_damage_getRateAndItsTangents enddo enddo -end procedure constitutive_damage_getRateAndItsTangents +end subroutine constitutive_damage_getRateAndItsTangents +!---------------------------------------------------------------------------------------------- +!< @brief writes damage sources resultsvto HDF5 output file +!---------------------------------------------------------------------------------------------- module subroutine damage_results integer :: p,i diff --git a/src/constitutive_plastic.f90 b/src/constitutive_plastic.f90 index 873887149..f62bca60e 100644 --- a/src/constitutive_plastic.f90 +++ b/src/constitutive_plastic.f90 @@ -168,7 +168,6 @@ submodule(constitutive) constitutive_plastic end subroutine plastic_nonlocal_results - end interface @@ -197,11 +196,18 @@ module subroutine plastic_init end subroutine plastic_init +!-------------------------------------------------------------------------------------------------- +!> @brief calls microstructure function of the different plasticity constitutive models +!-------------------------------------------------------------------------------------------------- +module subroutine constitutive_plastic_dependentState(F, Fp, ipc, ip, el) -!-------------------------------------------------------------------------------------------------- -!> @brief calls microstructure function of the different constitutive models -!-------------------------------------------------------------------------------------------------- -module procedure constitutive_plastic_dependentState + integer, intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + real(pReal), intent(in), dimension(3,3) :: & + F, & !< elastic deformation gradient + Fp !< plastic deformation gradient integer :: & ho, & !< homogenization @@ -222,14 +228,28 @@ module procedure constitutive_plastic_dependentState call plastic_nonlocal_dependentState (F,Fp,instance,of,ip,el) end select plasticityType -end procedure constitutive_plastic_dependentState +end subroutine constitutive_plastic_dependentState + !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the velocity gradient ! ToDo: Discuss whether it makes sense if crystallite handles the configuration conversion, i.e. ! Mp in, dLp_dMp out !-------------------------------------------------------------------------------------------------- -module procedure constitutive_plastic_LpAndItsTangents +module subroutine constitutive_plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & + S, Fi, ipc, ip, el) + integer, intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + real(pReal), intent(in), dimension(3,3) :: & + S, & !< 2nd Piola-Kirchhoff stress + Fi !< intermediate deformation gradient + real(pReal), intent(out), dimension(3,3) :: & + Lp !< plastic velocity gradient + real(pReal), intent(out), dimension(3,3,3,3) :: & + dLp_dS, & + dLp_dFi !< derivative of Lp with respect to Fi real(pReal), dimension(3,3,3,3) :: & dLp_dMp !< derivative of Lp with respect to Mandel stress @@ -280,8 +300,12 @@ module procedure constitutive_plastic_LpAndItsTangents dLp_dS(i,j,1:3,1:3) = matmul(matmul(transpose(Fi),Fi),dLp_dMp(i,j,1:3,1:3)) ! ToDo: @PS: why not: dLp_dMp:(FiT Fi) enddo; enddo -end procedure constitutive_plastic_LpAndItsTangents +end subroutine constitutive_plastic_LpAndItsTangents + +!-------------------------------------------------------------------------------------------- +!> @brief writes plasticity constitutive results to HDF5 output file +!-------------------------------------------------------------------------------------------- module subroutine plastic_results integer :: p diff --git a/src/constitutive_thermal.f90 b/src/constitutive_thermal.f90 index ca40e21fa..5972f5dc3 100644 --- a/src/constitutive_thermal.f90 +++ b/src/constitutive_thermal.f90 @@ -39,8 +39,12 @@ submodule(constitutive) constitutive_thermal end subroutine source_thermal_externalheat_getRateAndItsTangent end interface + contains +!---------------------------------------------------------------------------------------------- +!< @brief initializes thermal sources and kinematics mechanism +!---------------------------------------------------------------------------------------------- module subroutine thermal_init ! initialize source mechanisms @@ -54,7 +58,21 @@ module subroutine thermal_init end subroutine thermal_init -module procedure constitutive_thermal_getRateAndItsTangents +!---------------------------------------------------------------------------------------------- +!< @brief calculates thermal dissipation rate +!---------------------------------------------------------------------------------------------- +module subroutine constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, Tstar, Lp, ip, el) + integer, intent(in) :: & + ip, & !< integration point number + el !< element number + real(pReal), intent(in) :: & + T + real(pReal), intent(in), dimension(:,:,:,:,:) :: & + Tstar, & + Lp + real(pReal), intent(inout) :: & + TDot, & + dTDot_dT real(pReal) :: & my_Tdot, & @@ -94,6 +112,6 @@ module procedure constitutive_thermal_getRateAndItsTangents enddo enddo -end procedure constitutive_thermal_getRateAndItsTangents +end subroutine constitutive_thermal_getRateAndItsTangents end submodule From 0067b44a0df78b22f8976aeb1a0c8a12ab546bb2 Mon Sep 17 00:00:00 2001 From: Sharan Roongta Date: Sun, 12 Jul 2020 18:38:52 +0200 Subject: [PATCH 13/21] wrong logic --- src/constitutive_damage.f90 | 9 +-------- 1 file changed, 1 insertion(+), 8 deletions(-) diff --git a/src/constitutive_damage.f90 b/src/constitutive_damage.f90 index 5e1ed848d..8c36e4df2 100644 --- a/src/constitutive_damage.f90 +++ b/src/constitutive_damage.f90 @@ -189,24 +189,17 @@ module subroutine damage_results sourceLoop: do i = 1, phase_Nsources(p) group = trim('current/constituent')//'/'//trim(config_name_phase(p)) group = trim(group)//'/sources' + call results_closeGroup(results_addGroup(group)) sourceType: select case (phase_source(i,p)) case (SOURCE_damage_anisoBrittle_ID) sourceType - group = trim(group)//'/damage_anisoBrittle' - call results_closeGroup(results_addGroup(group)) call source_damage_anisoBrittle_results(p,group) case (SOURCE_damage_anisoDuctile_ID) sourceType - group = trim(group)//'/damage_anisoDuctile' - call results_closeGroup(results_addGroup(group)) call source_damage_anisoDuctile_results(p,group) case (SOURCE_damage_isoBrittle_ID) sourceType - group = trim(group)//'/damage_isoBrittle' - call results_closeGroup(results_addGroup(group)) call source_damage_isoBrittle_results(p,group) case (SOURCE_damage_isoDuctile_ID) sourceType - group = trim(group)//'/damage_isoDuctile' - call results_closeGroup(results_addGroup(group)) call source_damage_isoDuctile_results(p,group) end select sourceType From fcaa319f56adc6ff9bf1df48fdd304ddd9ae037e Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 13 Jul 2020 14:48:23 +0200 Subject: [PATCH 14/21] polishing --- src/commercialFEM_fileList.f90 | 6 +++--- src/constitutive.f90 | 21 +++++++++------------ src/constitutive_plastic.f90 | 4 +--- 3 files changed, 13 insertions(+), 18 deletions(-) diff --git a/src/commercialFEM_fileList.f90 b/src/commercialFEM_fileList.f90 index f02d7bc92..e98631fe4 100644 --- a/src/commercialFEM_fileList.f90 +++ b/src/commercialFEM_fileList.f90 @@ -28,8 +28,6 @@ #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" @@ -37,15 +35,17 @@ #include "constitutive_plastic_dislotwin.f90" #include "constitutive_plastic_disloUCLA.f90" #include "constitutive_plastic_nonlocal.f90" +#include "constitutive_thermal.f90" #include "source_thermal_dissipation.f90" #include "source_thermal_externalheat.f90" +#include "kinematics_thermal_expansion.f90" +#include "constitutive_damage.f90" #include "source_damage_isoBrittle.f90" #include "source_damage_isoDuctile.f90" #include "source_damage_anisoBrittle.f90" #include "source_damage_anisoDuctile.f90" #include "kinematics_cleavage_opening.f90" #include "kinematics_slipplane_opening.f90" -#include "kinematics_thermal_expansion.f90" #include "crystallite.f90" #include "thermal_isothermal.f90" #include "thermal_adiabatic.f90" diff --git a/src/constitutive.f90 b/src/constitutive.f90 index abff6e21d..2d60f7e66 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -368,8 +368,7 @@ subroutine constitutive_init debugConstitutive%ip = debug_root%get_asInt('integrationpoint',defaultVal = 1) debugConstitutive%grain = debug_root%get_asInt('grain',defaultVal = 1) -!-------------------------------------------------------------------------------------------------- -! initialized plasticity + call plastic_init call damage_init call thermal_init @@ -386,8 +385,7 @@ subroutine constitutive_init sourceState(ph)%p(s)%partionedState0 = sourceState(ph)%p(s)%state0 sourceState(ph)%p(s)%state = sourceState(ph)%p(s)%partionedState0 end forall -!-------------------------------------------------------------------------------------------------- -! determine max size of source state + constitutive_source_maxSizeDotState = max(constitutive_source_maxSizeDotState, & maxval(sourceState(ph)%p%sizeDotState)) enddo PhaseLoop2 @@ -401,12 +399,12 @@ end subroutine constitutive_init !-------------------------------------------------------------------------------------------------- function constitutive_homogenizedC(ipc,ip,el) - real(pReal) , dimension(6,6) :: & + real(pReal), dimension(6,6) :: & constitutive_homogenizedC integer, intent(in) :: & - ipc, & - ip, & - el + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element plasticityType: select case (phase_plasticity(material_phaseAt(ipc,el))) case (PLASTICITY_DISLOTWIN_ID) plasticityType @@ -628,9 +626,9 @@ function constitutive_collectDotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el real(pReal), dimension(3,3) :: & Mp integer :: & - ho, & - tme, & - i, & + ho, & !< homogenization + tme, & !< thermal member position + i, & !< counter in source loop instance logical :: broken @@ -685,7 +683,6 @@ function constitutive_collectDotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el enddo SourceLoop - end function constitutive_collectDotState diff --git a/src/constitutive_plastic.f90 b/src/constitutive_plastic.f90 index f62bca60e..acd9d5a67 100644 --- a/src/constitutive_plastic.f90 +++ b/src/constitutive_plastic.f90 @@ -175,12 +175,10 @@ contains !-------------------------------------------------------------------------------------------------- -!> @brief allocates arrays pointing to array of the various constitutive modules +!> @brief Initialize constitutive models for plasticity !-------------------------------------------------------------------------------------------------- module subroutine plastic_init -!-------------------------------------------------------------------------------------------------- -! initialized plasticity if (any(phase_plasticity == PLASTICITY_NONE_ID)) call plastic_none_init if (any(phase_plasticity == PLASTICITY_ISOTROPIC_ID)) call plastic_isotropic_init if (any(phase_plasticity == PLASTICITY_PHENOPOWERLAW_ID)) call plastic_phenopowerlaw_init From 50a7caa61a6ef39656e04a9cc5c7c527c44c96ca Mon Sep 17 00:00:00 2001 From: Sharan Roongta Date: Wed, 15 Jul 2020 14:35:21 +0200 Subject: [PATCH 15/21] cleaning --- src/constitutive.f90 | 136 ++++++++++++++++------------------- src/constitutive_damage.f90 | 52 ++++++-------- src/constitutive_plastic.f90 | 12 ++-- src/constitutive_thermal.f90 | 29 ++++---- src/damage_local.f90 | 2 +- src/damage_nonlocal.f90 | 2 +- src/thermal_adiabatic.f90 | 2 +- src/thermal_conduction.f90 | 2 +- 8 files changed, 106 insertions(+), 131 deletions(-) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 2d60f7e66..fadd9fdd6 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -1,7 +1,7 @@ !-------------------------------------------------------------------------------------------------- !> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH -!> @brief elasticity, plasticity, internal microstructure state +!> @brief elasticity, plasticity, damage & thermal internal microstructure state !-------------------------------------------------------------------------------------------------- module constitutive use prec @@ -81,8 +81,8 @@ module constitutive end subroutine plastic_disloUCLA_dotState module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature,timestep, & - instance,of,ip,el) - real(pReal), dimension(3,3), intent(in) ::& + instance,of,ip,el) + real(pReal), dimension(3,3), intent(in) :: & Mp !< MandelStress real(pReal), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem), intent(in) :: & F, & !< deformation gradient @@ -99,7 +99,6 @@ module constitutive module subroutine source_damage_anisoBrittle_dotState(S, ipc, ip, el) - integer, intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point @@ -109,7 +108,6 @@ module constitutive end subroutine source_damage_anisoBrittle_dotState module subroutine source_damage_anisoDuctile_dotState(ipc, ip, el) - integer, intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point @@ -117,7 +115,6 @@ module constitutive end subroutine source_damage_anisoDuctile_dotState module subroutine source_damage_isoDuctile_dotState(ipc, ip, el) - integer, intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point @@ -125,12 +122,36 @@ module constitutive end subroutine source_damage_isoDuctile_dotState module subroutine source_thermal_externalheat_dotState(phase, of) - integer, intent(in) :: & phase, & of end subroutine source_thermal_externalheat_dotState + module subroutine constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el) + integer, intent(in) :: & + ip, & !< integration point number + el !< element number + real(pReal), intent(in) :: & + phi !< damage parameter + real(pReal), intent(inout) :: & + phiDot, & + dPhiDot_dPhi + end subroutine constitutive_damage_getRateAndItsTangents + + module subroutine constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, S, Lp, ip, el) + integer, intent(in) :: & + ip, & !< integration point number + el !< element number + real(pReal), intent(in) :: & + T + real(pReal), intent(in), dimension(:,:,:,:,:) :: & + S, & !< current 2nd Piola Kitchoff stress vector + Lp !< plastic velocity gradient + real(pReal), intent(inout) :: & + TDot, & + dTDot_dT + end subroutine constitutive_thermal_getRateAndItsTangents + module function plastic_dislotwin_homogenizedC(ipc,ip,el) result(homogenizedC) real(pReal), dimension(6,6) :: & homogenizedC @@ -140,43 +161,13 @@ module constitutive el !< element end function plastic_dislotwin_homogenizedC - module subroutine constitutive_plastic_dependentState(F, Fp, ipc, ip, el) - - integer, intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element - real(pReal), intent(in), dimension(3,3) :: & - F, & !< elastic deformation gradient - Fp !< plastic deformation gradient - end subroutine constitutive_plastic_dependentState - - module subroutine constitutive_plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & - S, Fi, ipc, ip, el) - integer, intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element - real(pReal), intent(in), dimension(3,3) :: & - S, & !< 2nd Piola-Kirchhoff stress - Fi !< intermediate deformation gradient - real(pReal), intent(out), dimension(3,3) :: & - Lp !< plastic velocity gradient - real(pReal), intent(out), dimension(3,3,3,3) :: & - dLp_dS, & - dLp_dFi !< derivative of Lp with respect to Fi - - end subroutine constitutive_plastic_LpAndItsTangents - 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 subroutine plastic_nonlocal_updateCompatibility(orientation,instance,i,e) @@ -188,39 +179,11 @@ module constitutive orientation !< crystal orientation end subroutine plastic_nonlocal_updateCompatibility - module subroutine constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el) - integer, intent(in) :: & - ip, & !< integration point number - el !< element number - real(pReal), intent(in) :: & - phi - real(pReal), intent(inout) :: & - phiDot, & - dPhiDot_dPhi - - end subroutine constitutive_damage_getRateAndItsTangents - - module subroutine constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, Tstar, Lp, ip, el) - integer, intent(in) :: & - ip, & !< integration point number - el !< element number - real(pReal), intent(in) :: & - T - real(pReal), intent(in), dimension(:,:,:,:,:) :: & - Tstar, & - Lp - real(pReal), intent(inout) :: & - TDot, & - dTDot_dT - end subroutine constitutive_thermal_getRateAndItsTangents - - module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,of) real(pReal), dimension(3,3), intent(out) :: & Li !< inleastic velocity gradient real(pReal), dimension(3,3,3,3), intent(out) :: & dLi_dMi !< derivative of Li with respect to Mandel stress - real(pReal), dimension(3,3), intent(in) :: & Mi !< Mandel stress integer, intent(in) :: & @@ -229,7 +192,6 @@ module constitutive end subroutine plastic_isotropic_LiAndItsTangent module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, ip, el) - integer, intent(in) :: & ipc, & !< grain number ip, & !< integration point number @@ -243,7 +205,6 @@ module constitutive end subroutine kinematics_cleavage_opening_LiAndItsTangent module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, ip, el) - integer, intent(in) :: & ipc, & !< grain number ip, & !< integration point number @@ -257,7 +218,6 @@ module constitutive end subroutine kinematics_slipplane_opening_LiAndItsTangent module subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, ipc, ip, el) - integer, intent(in) :: & ipc, & !< grain number ip, & !< integration point number @@ -307,17 +267,40 @@ module constitutive end interface interface constitutive_LpAndItsTangents - module procedure :: constitutive_plastic_LpAndItsTangents + + module subroutine constitutive_plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & + S, Fi, ipc, ip, el) + integer, intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + real(pReal), intent(in), dimension(3,3) :: & + S, & !< 2nd Piola-Kirchhoff stress + Fi !< intermediate deformation gradient + real(pReal), intent(out), dimension(3,3) :: & + Lp !< plastic velocity gradient + real(pReal), intent(out), dimension(3,3,3,3) :: & + dLp_dS, & + dLp_dFi !< derivative of Lp with respect to Fi + end subroutine constitutive_plastic_LpAndItsTangents + end interface constitutive_LpAndItsTangents + interface constitutive_dependentState - module procedure :: constitutive_plastic_dependentState + + module subroutine constitutive_plastic_dependentState(F, Fp, ipc, ip, el) + integer, intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + real(pReal), intent(in), dimension(3,3) :: & + F, & !< elastic deformation gradient + Fp !< plastic deformation gradient + end subroutine constitutive_plastic_dependentState + end interface constitutive_dependentState - interface constitutive_getRateAndItsTangents - module procedure :: constitutive_damage_getRateAndItsTangents , & - constitutive_thermal_getRateAndItsTangents - end interface constitutive_getRateAndItsTangents type :: tDebugOptions logical :: & @@ -343,14 +326,15 @@ module constitutive constitutive_collectDotState, & constitutive_deltaState, & plastic_nonlocal_updateCompatibility, & - constitutive_getRateAndItsTangents, & + constitutive_damage_getRateAndItsTangents, & + constitutive_thermal_getRateAndItsTangents, & constitutive_results contains !-------------------------------------------------------------------------------------------------- -!> @brief allocates arrays pointing to array of the various constitutive modules +!> @brief Initialze constitutive models for individual physics !-------------------------------------------------------------------------------------------------- subroutine constitutive_init diff --git a/src/constitutive_damage.f90 b/src/constitutive_damage.f90 index 8c36e4df2..9e0c686b0 100644 --- a/src/constitutive_damage.f90 +++ b/src/constitutive_damage.f90 @@ -1,3 +1,6 @@ +!---------------------------------------------------------------------------------------------------- +!> @brief internal microstructure state for all damage sources and kinematics constitutive models +!---------------------------------------------------------------------------------------------------- submodule(constitutive) constitutive_damage interface @@ -26,7 +29,7 @@ submodule(constitutive) constitutive_damage phase, & !< phase ID of element constituent !< position of element within its phase instance real(pReal), intent(in) :: & - phi !< damage value + phi !< damage parameter real(pReal), intent(out) :: & localphiDot, & dLocalphiDot_dPhi @@ -34,10 +37,10 @@ submodule(constitutive) constitutive_damage module subroutine source_damage_anisoDuctile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) integer, intent(in) :: & - phase, & - constituent + phase, & !< phase ID of element + constituent !< position of element within its phase instance real(pReal), intent(in) :: & - phi + phi !< damage parameter real(pReal), intent(out) :: & localphiDot, & dLocalphiDot_dPhi @@ -45,10 +48,10 @@ submodule(constitutive) constitutive_damage module subroutine source_damage_isoBrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) integer, intent(in) :: & - phase, & - constituent + phase, & !< phase ID of element + constituent !< position of element within its phase instance real(pReal), intent(in) :: & - phi + phi !< damage parameter real(pReal), intent(out) :: & localphiDot, & dLocalphiDot_dPhi @@ -56,28 +59,15 @@ submodule(constitutive) constitutive_damage module subroutine source_damage_isoDuctile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) integer, intent(in) :: & - phase, & - constituent + phase, & !< phase ID of element + constituent !< position of element within its phase instance real(pReal), intent(in) :: & - phi + phi !< damage parameter real(pReal), intent(out) :: & localphiDot, & dLocalphiDot_dPhi end subroutine source_damage_isoDuctile_getRateAndItsTangent - module subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDot_dT, Tstar, Lp, phase) - integer, intent(in) :: & - phase - real(pReal), intent(in), dimension(3,3) :: & - Tstar - real(pReal), intent(in), dimension(3,3) :: & - Lp - - real(pReal), intent(out) :: & - TDot, & - dTDot_dT - end subroutine source_thermal_dissipation_getRateAndItsTangent - module subroutine source_damage_anisoBrittle_results(phase,group) integer, intent(in) :: phase character(len=*), intent(in) :: group @@ -108,10 +98,10 @@ contains module subroutine damage_init ! initialize source mechanisms - if (any(phase_source == SOURCE_damage_isoBrittle_ID)) call source_damage_isoBrittle_init - if (any(phase_source == SOURCE_damage_isoDuctile_ID)) call source_damage_isoDuctile_init - if (any(phase_source == SOURCE_damage_anisoBrittle_ID)) call source_damage_anisoBrittle_init - if (any(phase_source == SOURCE_damage_anisoDuctile_ID)) call source_damage_anisoDuctile_init + if (any(phase_source == SOURCE_damage_isoBrittle_ID)) call source_damage_isoBrittle_init + if (any(phase_source == SOURCE_damage_isoDuctile_ID)) call source_damage_isoDuctile_init + if (any(phase_source == SOURCE_damage_anisoBrittle_ID)) call source_damage_anisoBrittle_init + if (any(phase_source == SOURCE_damage_anisoDuctile_ID)) call source_damage_anisoDuctile_init !-------------------------------------------------------------------------------------------------- ! initialize kinematic mechanisms @@ -127,10 +117,10 @@ end subroutine damage_init module subroutine constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el) integer, intent(in) :: & - ip, & !< integration point number - el !< element number + ip, & !< integration point number + el !< element number real(pReal), intent(in) :: & - phi + phi !< damage parameter real(pReal), intent(inout) :: & phiDot, & dPhiDot_dPhi @@ -209,4 +199,4 @@ module subroutine damage_results end subroutine damage_results -end submodule +end submodule constitutive_damage diff --git a/src/constitutive_plastic.f90 b/src/constitutive_plastic.f90 index acd9d5a67..f7b1569d2 100644 --- a/src/constitutive_plastic.f90 +++ b/src/constitutive_plastic.f90 @@ -1,3 +1,6 @@ +!---------------------------------------------------------------------------------------------------- +!> @brief internal microstructure state for all plasticity constitutive models +!---------------------------------------------------------------------------------------------------- submodule(constitutive) constitutive_plastic interface @@ -42,7 +45,6 @@ submodule(constitutive) constitutive_plastic Lp !< plastic velocity gradient real(pReal), dimension(3,3,3,3), intent(out) :: & dLp_dMp !< derivative of Lp with respect to the Mandel stress - real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & @@ -128,13 +130,13 @@ submodule(constitutive) constitutive_plastic module subroutine plastic_nonlocal_dependentState(F, Fp, instance, of, ip, el) real(pReal), dimension(3,3), intent(in) :: & - F, & - Fp + F, & !< deformation gradient + Fp !< plastic deformation gradient integer, intent(in) :: & instance, & of, & - ip, & - el + ip, & !< current integration point + el !< current element number end subroutine plastic_nonlocal_dependentState module subroutine plastic_isotropic_results(instance,group) diff --git a/src/constitutive_thermal.f90 b/src/constitutive_thermal.f90 index 5972f5dc3..3b3398ce2 100644 --- a/src/constitutive_thermal.f90 +++ b/src/constitutive_thermal.f90 @@ -1,3 +1,6 @@ +!---------------------------------------------------------------------------------------------------- +!> @brief internal microstructure state for all thermal sources and kinematics constitutive models +!---------------------------------------------------------------------------------------------------- submodule(constitutive) constitutive_thermal interface @@ -13,29 +16,24 @@ submodule(constitutive) constitutive_thermal module subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDot_dT, Tstar, Lp, phase) - integer, intent(in) :: & - phase + phase !< phase ID of element real(pReal), intent(in), dimension(3,3) :: & - Tstar + Tstar !< 2nd Piola Kirchoff stress tensor for a given element real(pReal), intent(in), dimension(3,3) :: & - Lp - + Lp !< plastic velocuty gradient for a given element real(pReal), intent(out) :: & TDot, & dTDot_dT - end subroutine source_thermal_dissipation_getRateAndItsTangent module subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, dTDot_dT, phase, of) - integer, intent(in) :: & phase, & of real(pReal), intent(out) :: & TDot, & dTDot_dT - end subroutine source_thermal_externalheat_getRateAndItsTangent end interface @@ -61,15 +59,15 @@ end subroutine thermal_init !---------------------------------------------------------------------------------------------- !< @brief calculates thermal dissipation rate !---------------------------------------------------------------------------------------------- -module subroutine constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, Tstar, Lp, ip, el) +module subroutine constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, S, Lp, ip, el) integer, intent(in) :: & - ip, & !< integration point number - el !< element number + ip, & !< integration point number + el !< element number real(pReal), intent(in) :: & T real(pReal), intent(in), dimension(:,:,:,:,:) :: & - Tstar, & - Lp + S, & !< current 2nd Piola Kirchoff stress + Lp !< plastic velocity gradient real(pReal), intent(inout) :: & TDot, & dTDot_dT @@ -95,7 +93,7 @@ module subroutine constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, select case(phase_source(source,phase)) case (SOURCE_thermal_dissipation_ID) call source_thermal_dissipation_getRateAndItsTangent(my_Tdot, my_dTdot_dT, & - Tstar(1:3,1:3,grain,ip,el), & + S(1:3,1:3,grain,ip,el), & Lp(1:3,1:3,grain,ip,el), & phase) @@ -114,4 +112,5 @@ module subroutine constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, end subroutine constitutive_thermal_getRateAndItsTangents -end submodule + +end submodule constitutive_thermal diff --git a/src/damage_local.f90 b/src/damage_local.f90 index 3d448f7a4..66b50064b 100644 --- a/src/damage_local.f90 +++ b/src/damage_local.f90 @@ -131,7 +131,7 @@ subroutine damage_local_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el phiDot = 0.0_pReal dPhiDot_dPhi = 0.0_pReal - call constitutive_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el) + call constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el) phiDot = phiDot/real(homogenization_Ngrains(material_homogenizationAt(el)),pReal) dPhiDot_dPhi = dPhiDot_dPhi/real(homogenization_Ngrains(material_homogenizationAt(el)),pReal) diff --git a/src/damage_nonlocal.f90 b/src/damage_nonlocal.f90 index 530a79899..914133de7 100644 --- a/src/damage_nonlocal.f90 +++ b/src/damage_nonlocal.f90 @@ -100,7 +100,7 @@ subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, phiDot = 0.0_pReal dPhiDot_dPhi = 0.0_pReal - call constitutive_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el) + call constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el) phiDot = phiDot/real(homogenization_Ngrains(material_homogenizationAt(el)),pReal) dPhiDot_dPhi = dPhiDot_dPhi/real(homogenization_Ngrains(material_homogenizationAt(el)),pReal) diff --git a/src/thermal_adiabatic.f90 b/src/thermal_adiabatic.f90 index 9e4710a1e..c52f0a3d0 100644 --- a/src/thermal_adiabatic.f90 +++ b/src/thermal_adiabatic.f90 @@ -131,7 +131,7 @@ subroutine thermal_adiabatic_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el) dTdot_dT = 0.0_pReal homog = material_homogenizationAt(el) - call constitutive_getRateAndItsTangents(TDot, dTDot_dT, T, crystallite_S, crystallite_Lp, ip, el) + call constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, crystallite_S, crystallite_Lp, ip, el) Tdot = Tdot/real(homogenization_Ngrains(homog),pReal) dTdot_dT = dTdot_dT/real(homogenization_Ngrains(homog),pReal) diff --git a/src/thermal_conduction.f90 b/src/thermal_conduction.f90 index 2afe411f1..60766710d 100644 --- a/src/thermal_conduction.f90 +++ b/src/thermal_conduction.f90 @@ -90,7 +90,7 @@ subroutine thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el) dTdot_dT = 0.0_pReal homog = material_homogenizationAt(el) - call constitutive_getRateAndItsTangents(TDot, dTDot_dT, T, crystallite_S,crystallite_Lp ,ip, el) + call constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, crystallite_S,crystallite_Lp ,ip, el) Tdot = Tdot/real(homogenization_Ngrains(homog),pReal) dTdot_dT = dTdot_dT/real(homogenization_Ngrains(homog),pReal) From 33d42265f2ab5bb7b9a6be2255be6c75a5c9d400 Mon Sep 17 00:00:00 2001 From: Test User Date: Wed, 15 Jul 2020 22:31:40 +0200 Subject: [PATCH 16/21] [skip ci] updated version information after successful test of v2.0.3-2889-gf0fdcd0d --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index b76d77510..403eb9a88 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.3-2872-g0a2d3046 +v2.0.3-2889-gf0fdcd0d From ec56316683f17f63bed7c753b155237ec2667809 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 16 Jul 2020 15:26:00 +0200 Subject: [PATCH 17/21] IPneighborhood for MSC.Marc tested for 8 ip hexahedaron --- src/marc/discretization_marc.f90 | 88 ++++++++++++++++++++++++++------ src/math.f90 | 18 ++++--- 2 files changed, 83 insertions(+), 23 deletions(-) diff --git a/src/marc/discretization_marc.f90 b/src/marc/discretization_marc.f90 index 365bbbff7..d53dccf75 100644 --- a/src/marc/discretization_marc.f90 +++ b/src/marc/discretization_marc.f90 @@ -119,7 +119,7 @@ subroutine discretization_marc_init unscaledNormals = IPareaNormal(elem,nElems,connectivity_cell,node0_cell) call geometry_plastic_nonlocal_setIParea(norm2(unscaledNormals,1)) call geometry_plastic_nonlocal_setIPareaNormal(unscaledNormals/spread(norm2(unscaledNormals,1),1,3)) - !call geometry_plastic_nonlocal_setIPneighborhood ToDo: Support nonlocal + call geometry_plastic_nonlocal_setIPneighborhood(IPneighborhood(elem,connectivity_cell)) call geometry_plastic_nonlocal_results end subroutine discretization_marc_init @@ -733,10 +733,10 @@ end subroutine inputRead_microstructureAndHomogenization !-------------------------------------------------------------------------------------------------- -!> @brief Determine cell connectivity and definition of cell nodes +!> @brief Calculates cell node coordinates from element node coordinates !-------------------------------------------------------------------------------------------------- -subroutine buildCells(connectivity_cell,cellNodeDefinition, & - elem,connectivity_elem) +pure subroutine buildCells(connectivity_cell,cellNodeDefinition, & + elem,connectivity_elem) type(tCellNodeDefinition), dimension(:), intent(out) :: cellNodeDefinition ! definition of cell nodes for increasing number of parents integer, dimension(:,:,:),intent(out) :: connectivity_cell @@ -747,8 +747,7 @@ subroutine buildCells(connectivity_cell,cellNodeDefinition, & integer,dimension(:), allocatable :: candidates_local integer,dimension(:,:), allocatable :: parentsAndWeights,candidates_global - integer :: e,n,c,p,s,i,m,j,& - nParentNodes,nCellNode,Nelem,candidateID + integer :: e, n, c, p, s,i,m,j,nParentNodes,nCellNode,Nelem,candidateID Nelem = size(connectivity_elem,2) @@ -762,7 +761,9 @@ subroutine buildCells(connectivity_cell,cellNodeDefinition, & do e = 1, Nelem do c = 1, elem%NcellNodes realNode: if (count(elem%cellNodeParentNodeWeights(:,c) /= 0) == 1) then - where(connectivity_cell(:,:,e) == -c) connectivity_cell(:,:,e) = connectivity_elem(c,e) + where(connectivity_cell(:,:,e) == -c) + connectivity_cell(:,:,e) = connectivity_elem(c,e) + end where endif realNode enddo enddo @@ -889,10 +890,10 @@ end subroutine buildCells !-------------------------------------------------------------------------------------------------- -!> @brief Calculate cell node coordinates from element node coordinates +!> @brief Calculates cell node coordinates from element node coordinates !-------------------------------------------------------------------------------------------------- -subroutine buildCellNodes(node_cell, & - definition,node_elem) +pure subroutine buildCellNodes(node_cell, & + definition,node_elem) real(pReal), dimension(:,:), intent(out) :: node_cell !< cell node coordinates type(tCellNodeDefinition), dimension(:), intent(in) :: definition !< cell node definition (weights and parents) @@ -919,10 +920,10 @@ end subroutine buildCellNodes !-------------------------------------------------------------------------------------------------- -!> @brief Calculate IP coordinates as center of cell +!> @brief Calculates IP coordinates as center of cell !-------------------------------------------------------------------------------------------------- -subroutine buildIPcoordinates(IPcoordinates, & - connectivity_cell,node_cell) +pure subroutine buildIPcoordinates(IPcoordinates, & + connectivity_cell,node_cell) real(pReal), dimension(:,:), intent(out):: IPcoordinates !< cell-center/IP coordinates integer, dimension(:,:), intent(in) :: connectivity_cell !< connectivity for each cell @@ -947,7 +948,7 @@ end subroutine buildIPcoordinates !> @details The IP volume is calculated differently depending on the cell type. !> 2D cells assume an element depth of 1.0 !--------------------------------------------------------------------------------------------------- -function IPvolume(elem,node,connectivity) +pure function IPvolume(elem,node,connectivity) type(tElement), intent(in) :: elem real(pReal), dimension(:,:), intent(in) :: node @@ -1005,7 +1006,7 @@ end function IPvolume !-------------------------------------------------------------------------------------------------- !> @brief calculation of IP interface areas !-------------------------------------------------------------------------------------------------- -function IPareaNormal(elem,nElem,connectivity,node) +pure function IPareaNormal(elem,nElem,connectivity,node) type(tElement), intent(in) :: elem integer, intent(in) :: nElem @@ -1051,6 +1052,63 @@ function IPareaNormal(elem,nElem,connectivity,node) end function IPareaNormal +!-------------------------------------------------------------------------------------------------- +!> @brief IP neighborhood +!-------------------------------------------------------------------------------------------------- +function IPneighborhood(elem,connectivity) + + type(tElement), intent(in) :: elem ! definition of the element in use + integer, dimension(:,:,:), intent(in) :: connectivity ! cell connectivity + integer, dimension(3,size(elem%cellFace,2), & + size(connectivity,2),size(connectivity,3)) :: IPneighborhood ! neighboring IPs as [element ID, IP ID, face ID] + + integer, dimension(size(elem%cellFace,1)+3,& + size(elem%cellFace,2)*size(connectivity,2)*size(connectivity,3)) :: face + integer, dimension(size(connectivity,1)) :: myConnectivity + integer, dimension(size(elem%cellFace,1)) :: face_unordered + integer :: e,i,f,n,c,s + + c = 0 + do e = 1, size(connectivity,3) + do i = 1, size(connectivity,2) + myConnectivity = connectivity(:,i,e) + do f = 1, size(elem%cellFace,2) + c = c + 1 + face_unordered = myConnectivity(elem%cellFace(:,f)) + do n = 1, size(face_unordered) + face(n,c) = minval(face_unordered) + face_unordered(minloc(face_unordered)) = huge(face_unordered) + enddo + face(n:n+3,c) = [e,i,f] + enddo + enddo; enddo + +!-------------------------------------------------------------------------------------------------- +! sort face definitions + call math_sort(face,sortDim=1) + do c=2, size(face,1)-4 + s = 1 + e = 1 + do while (e < size(face,2)) + e = e + 1 + if(any(face(:c,s) /= face(:c,e))) then + if(e-1/=s) call math_sort(face(:,s:e-1),sortDim=c) + s = e + endif + enddo + enddo + + IPneighborhood = 0 + do c=1, size(face,2) - 1 + if(all(face(:n-1,c) == face(:n-1,c+1))) then + IPneighborhood(:,face(n+2,c+1),face(n+1,c+1),face(n+0,c+1)) = face(n:n+3,c+0) + IPneighborhood(:,face(n+2,c+0),face(n+1,c+0),face(n+0,c+0)) = face(n:n+3,c+1) + endif + enddo + +end function IPneighborhood + + !-------------------------------------------------------------------------------------------------- !> @brief return integer list corresponding to items in consecutive lines. !! First integer in array is counter diff --git a/src/math.f90 b/src/math.f90 index 7abb39424..0ef13cc3f 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -126,11 +126,12 @@ end subroutine math_init !-------------------------------------------------------------------------------------------------- -!> @brief Quicksort algorithm for two-dimensional integer arrays -! Sorting is done with respect to array(sort,:) and keeps array(/=sort,:) linked to it. -! default: sort=1 +!> @brief Sorting of two-dimensional integer arrays +!> @details Based on quicksort. +! Sorting is done with respect to array(sortDim,:) and keeps array(/=sortDim,:) linked to it. +! Default: sortDim=1 !-------------------------------------------------------------------------------------------------- -recursive subroutine math_sort(a, istart, iend, sortDim) +pure recursive subroutine math_sort(a, istart, iend, sortDim) integer, dimension(:,:), intent(inout) :: a integer, intent(in),optional :: istart,iend, sortDim @@ -155,7 +156,7 @@ recursive subroutine math_sort(a, istart, iend, sortDim) endif if (s < e) then - ipivot = qsort_partition(a,s, e, d) + call qsort_partition(a,ipivot, s,e, d) call math_sort(a, s, ipivot-1, d) call math_sort(a, ipivot+1, e, d) endif @@ -166,9 +167,10 @@ recursive subroutine math_sort(a, istart, iend, sortDim) !------------------------------------------------------------------------------------------------- !> @brief Partitioning required for quicksort !------------------------------------------------------------------------------------------------- - integer function qsort_partition(a, istart, iend, sort) + pure subroutine qsort_partition(a,p, istart, iend, sort) integer, dimension(:,:), intent(inout) :: a + integer, intent(out) :: p ! Pivot element integer, intent(in) :: istart,iend,sort integer, dimension(size(a,1)) :: tmp integer :: i,j @@ -186,7 +188,7 @@ recursive subroutine math_sort(a, istart, iend, sortDim) tmp = a(:,istart) a(:,istart) = a(:,j) a(:,j) = tmp - qsort_partition = j + p = j return else cross ! exchange values tmp = a(:,i) @@ -195,7 +197,7 @@ recursive subroutine math_sort(a, istart, iend, sortDim) endif cross enddo - end function qsort_partition + end subroutine qsort_partition end subroutine math_sort From d8e112c04262083697bcbf12a33edd973e3c2a9d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 16 Jul 2020 23:38:37 +0200 Subject: [PATCH 18/21] phaseAt needs constituent/element tuple constituent is always 1 for nonlocal --- src/crystallite.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index fbce3ab47..f33f392a2 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -641,7 +641,7 @@ subroutine crystallite_orientations if (plasticState(material_phaseAt(1,e))%nonlocal) then do i = FEsolving_execIP(1),FEsolving_execIP(2) call plastic_nonlocal_updateCompatibility(crystallite_orientation, & - phase_plasticityInstance(material_phaseAt(i,e)),i,e) + phase_plasticityInstance(material_phaseAt(1,e)),i,e) enddo endif enddo From 63cb5f1c9cdfc7361f5114efb3bd4cf1c7d6b960 Mon Sep 17 00:00:00 2001 From: Test User Date: Fri, 17 Jul 2020 13:06:46 +0200 Subject: [PATCH 19/21] [skip ci] updated version information after successful test of v2.0.3-2912-g361c35b8 --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index e78e2c589..e00afd09a 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.3-2881-gc07efe84 +v2.0.3-2912-g361c35b8 From 4b94a8c27f75827d2b83ef33237278bb5f820d2c Mon Sep 17 00:00:00 2001 From: Test User Date: Fri, 17 Jul 2020 18:43:22 +0200 Subject: [PATCH 20/21] [skip ci] updated version information after successful test of v2.0.3-2925-gd47273f5 --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index e00afd09a..60fe3e0e3 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.3-2912-g361c35b8 +v2.0.3-2925-gd47273f5 From 83504b417d339559779613ac84377223005ceb3d Mon Sep 17 00:00:00 2001 From: Test User Date: Mon, 20 Jul 2020 13:41:54 +0200 Subject: [PATCH 21/21] [skip ci] updated version information after successful test of v2.0.3-2929-gfb02c69a --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 60fe3e0e3..fc6f90a9e 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.3-2925-gd47273f5 +v2.0.3-2929-gfb02c69a