diff --git a/VERSION b/VERSION index e78e2c589..fc6f90a9e 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.3-2881-gc07efe84 +v2.0.3-2929-gfb02c69a diff --git a/src/commercialFEM_fileList.f90 b/src/commercialFEM_fileList.f90 index 942363a9f..e98631fe4 100644 --- a/src/commercialFEM_fileList.f90 +++ b/src/commercialFEM_fileList.f90 @@ -26,16 +26,8 @@ #endif #include "material.f90" #include "lattice.f90" -#include "source_thermal_dissipation.f90" -#include "source_thermal_externalheat.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 "constitutive.f90" +#include "constitutive_plastic.f90" #include "constitutive_plastic_none.f90" #include "constitutive_plastic_isotropic.f90" #include "constitutive_plastic_phenopowerlaw.f90" @@ -43,6 +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 "crystallite.f90" #include "thermal_isothermal.f90" #include "thermal_adiabatic.f90" diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 128b45d53..feff4af86 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 @@ -17,15 +17,6 @@ module constitutive use discretization use geometry_plastic_nonlocal, only: & geometry_plastic_nonlocal_disable - 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 @@ -36,128 +27,14 @@ 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 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_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) :: & - instance, & - of - end subroutine plastic_isotropic_LiAndItsTangent + module subroutine thermal_init + end subroutine thermal_init module subroutine plastic_isotropic_dotState(Mp,instance,of) @@ -205,8 +82,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 @@ -221,31 +98,136 @@ module constitutive 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 plastic_dislotwin_dependentState(T,instance,of) - integer, intent(in) :: & - instance, & + 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 - real(pReal), intent(in) :: & + 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 - end subroutine plastic_dislotwin_dependentState + 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 subroutine plastic_disloUCLA_dependentState(instance,of) - integer, intent(in) :: & - instance, & - of - end subroutine plastic_disloUCLA_dependentState + 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_dependentState(F, Fp, instance, of, ip, el) - real(pReal), dimension(3,3), intent(in) :: & - F, & - Fp + 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, & - of, & - ip, & - el - end subroutine plastic_nonlocal_dependentState + i, & + e + type(rotation), dimension(1,discretization_nIP,discretization_nElem), intent(in) :: & + orientation !< crystal orientation + end subroutine plastic_nonlocal_updateCompatibility + + 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) :: & + instance, & + 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_kinehardening_deltaState(Mp,instance,of) @@ -266,58 +248,60 @@ 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) :: & + 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 - end function plastic_dislotwin_homogenizedC + 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_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 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 plastic_results + end subroutine plastic_results + + module subroutine damage_results + end subroutine damage_results end interface + interface constitutive_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 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 + type :: tDebugOptions logical :: & @@ -333,23 +317,25 @@ 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 !-------------------------------------------------------------------------------------------------- -!> @brief allocates arrays pointing to array of the various constitutive modules +!> @brief Initialze constitutive models for individual physics !-------------------------------------------------------------------------------------------------- subroutine constitutive_init @@ -367,33 +353,10 @@ subroutine constitutive_init debugConstitutive%ip = debug_root%get_asInt('integrationpoint',defaultVal = 1) debugConstitutive%grain = debug_root%get_asInt('grain',defaultVal = 1) -!-------------------------------------------------------------------------------------------------- -! 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) @@ -407,8 +370,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 @@ -416,134 +378,29 @@ 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) :: & + real(pReal), dimension(6,6) :: & + constitutive_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 - constitutive_homogenizedC = plastic_dislotwin_homogenizedC(ipc,ip,el) + constitutive_homogenizedC = plastic_dislotwin_homogenizedC(ipc,ip,el) case default plasticityType - constitutive_homogenizedC = lattice_C66(1:6,1:6,material_phaseAt(ipc,el)) + constitutive_homogenizedC = lattice_C66(1:6,1:6,material_phaseAt(ipc,el)) end select plasticityType 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? @@ -739,7 +596,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) :: & @@ -794,7 +651,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) @@ -897,37 +754,8 @@ end function constitutive_deltaState !-------------------------------------------------------------------------------------------------- subroutine constitutive_results - integer :: p - character(len=pStringLen) :: group - 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 + call plastic_results + call damage_results end subroutine constitutive_results diff --git a/src/constitutive_damage.f90 b/src/constitutive_damage.f90 new file mode 100644 index 000000000..9e0c686b0 --- /dev/null +++ b/src/constitutive_damage.f90 @@ -0,0 +1,202 @@ +!---------------------------------------------------------------------------------------------------- +!> @brief internal microstructure state for all damage sources and kinematics constitutive models +!---------------------------------------------------------------------------------------------------- +submodule(constitutive) constitutive_damage + + 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, & !< phase ID of element + constituent !< position of element within its phase instance + real(pReal), intent(in) :: & + phi !< damage parameter + 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, & !< phase ID of element + constituent !< position of element within its phase instance + real(pReal), intent(in) :: & + phi !< damage parameter + 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, & !< phase ID of element + constituent !< position of element within its phase instance + real(pReal), intent(in) :: & + phi !< damage parameter + 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, & !< phase ID of element + constituent !< position of element within its phase instance + real(pReal), intent(in) :: & + phi !< damage parameter + real(pReal), intent(out) :: & + localphiDot, & + dLocalphiDot_dPhi + end subroutine source_damage_isoDuctile_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 + +!---------------------------------------------------------------------------------------------- +!< @brief initialize damage sources and kinematics mechanism +!---------------------------------------------------------------------------------------------- +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 + + +!---------------------------------------------------------------------------------------------- +!< @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 !< damage parameter + 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 constitutive_damage_getRateAndItsTangents + + +!---------------------------------------------------------------------------------------------- +!< @brief writes damage sources resultsvto HDF5 output file +!---------------------------------------------------------------------------------------------- +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' + call results_closeGroup(results_addGroup(group)) + + sourceType: select case (phase_source(i,p)) + + case (SOURCE_damage_anisoBrittle_ID) sourceType + call source_damage_anisoBrittle_results(p,group) + case (SOURCE_damage_anisoDuctile_ID) sourceType + call source_damage_anisoDuctile_results(p,group) + case (SOURCE_damage_isoBrittle_ID) sourceType + call source_damage_isoBrittle_results(p,group) + case (SOURCE_damage_isoDuctile_ID) sourceType + call source_damage_isoDuctile_results(p,group) + end select sourceType + + enddo SourceLoop + enddo + +end subroutine damage_results + + +end submodule constitutive_damage diff --git a/src/constitutive_plastic.f90 b/src/constitutive_plastic.f90 new file mode 100644 index 000000000..f7b1569d2 --- /dev/null +++ b/src/constitutive_plastic.f90 @@ -0,0 +1,348 @@ +!---------------------------------------------------------------------------------------------------- +!> @brief internal microstructure state for all plasticity constitutive models +!---------------------------------------------------------------------------------------------------- +submodule(constitutive) constitutive_plastic + + 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, & !< deformation gradient + Fp !< plastic deformation gradient + integer, intent(in) :: & + instance, & + of, & + ip, & !< current integration point + el !< current element number + 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 + + +contains + + +!-------------------------------------------------------------------------------------------------- +!> @brief Initialize constitutive models for plasticity +!-------------------------------------------------------------------------------------------------- +module subroutine plastic_init + + 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 plasticity constitutive models +!-------------------------------------------------------------------------------------------------- +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 + + 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_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 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 + 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_plastic_LpAndItsTangents + + +!-------------------------------------------------------------------------------------------- +!> @brief writes plasticity constitutive results to HDF5 output file +!-------------------------------------------------------------------------------------------- +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_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 dfd9dbaac..4971b48f5 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..3b3398ce2 --- /dev/null +++ b/src/constitutive_thermal.f90 @@ -0,0 +1,116 @@ +!---------------------------------------------------------------------------------------------------- +!> @brief internal microstructure state for all thermal sources and kinematics constitutive models +!---------------------------------------------------------------------------------------------------- +submodule(constitutive) constitutive_thermal + + 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 !< phase ID of element + real(pReal), intent(in), dimension(3,3) :: & + Tstar !< 2nd Piola Kirchoff stress tensor for a given element + real(pReal), intent(in), dimension(3,3) :: & + 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 + +contains + +!---------------------------------------------------------------------------------------------- +!< @brief initializes thermal sources and kinematics mechanism +!---------------------------------------------------------------------------------------------- +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 + + +!---------------------------------------------------------------------------------------------- +!< @brief calculates thermal dissipation rate +!---------------------------------------------------------------------------------------------- +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 Kirchoff stress + Lp !< plastic velocity gradient + 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, & + S(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 constitutive_thermal_getRateAndItsTangents + + +end submodule constitutive_thermal 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 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..8fc2dc704 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,22 +66,22 @@ 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) +pure module function kinematics_thermal_expansion_initialStrain(homog,phase,offset) result(initialStrain) - integer, intent(in) :: & - phase, & - homog, & - offset + integer, intent(in) :: & + phase, & + homog, & + offset - real(pReal), dimension(3,3) :: & - kinematics_thermal_expansion_initialStrain !< initial thermal strain (should be small strain, though) + real(pReal), dimension(3,3) :: & + 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 + 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 @@ -102,7 +89,7 @@ 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/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 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)