From 11d7f034e47d04a857aa1ef26f62aa6a9fafe461 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 22 Dec 2020 12:20:00 +0100 Subject: [PATCH] code follows modular structure --- src/constitutive.f90 | 69 +++++++++++++++++++++++++++---------- src/constitutive_damage.f90 | 69 +++++++++++++++++++------------------ 2 files changed, 86 insertions(+), 52 deletions(-) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index fe978e2ef..e68ade06f 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -152,6 +152,12 @@ module constitutive integer, intent(in) :: ph end subroutine mech_results + module subroutine damage_results(group,ph) + character(len=*), intent(in) :: group + integer, intent(in) :: ph + end subroutine damage_results + + module subroutine mech_restart_read(fileHandle) integer(HID_T), intent(in) :: fileHandle end subroutine mech_restart_read @@ -314,10 +320,6 @@ end function constitutive_deltaState C end subroutine source_damage_isoBrittle_deltaState - module subroutine damage_results - end subroutine damage_results - - module subroutine constitutive_plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & S, Fi, ipc, ip, el) @@ -468,7 +470,7 @@ end subroutine constitutive_init !-------------------------------------------------------------------------------------------------- !> @brief checks if a source mechanism is active or not !-------------------------------------------------------------------------------------------------- -module function source_active(source_label,src_length) result(active_source) +function source_active(source_label,src_length) result(active_source) character(len=*), intent(in) :: source_label !< name of source mechanism integer, intent(in) :: src_length !< max. number of sources in system @@ -499,8 +501,7 @@ end function source_active !-------------------------------------------------------------------------------------------------- !> @brief checks if a kinematic mechanism is active or not !-------------------------------------------------------------------------------------------------- - -module function kinematics_active(kinematics_label,kinematics_length) result(active_kinematics) +function kinematics_active(kinematics_label,kinematics_length) result(active_kinematics) character(len=*), intent(in) :: kinematics_label !< name of kinematic mechanism integer, intent(in) :: kinematics_length !< max. number of kinematics in system @@ -631,13 +632,10 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & end subroutine constitutive_LiAndItsTangents - - - !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the rate of change of microstructure !-------------------------------------------------------------------------------------------------- -function constitutive_source_collectDotState(S, ipc, ip, el,phase,of) result(broken) +function constitutive_damage_collectDotState(S, ipc, ip, el,phase,of) result(broken) integer, intent(in) :: & ipc, & !< component-ID of integration point @@ -667,6 +665,39 @@ function constitutive_source_collectDotState(S, ipc, ip, el,phase,of) result(bro case (SOURCE_damage_anisoDuctile_ID) sourceType call source_damage_anisoDuctile_dotState(ipc, ip, el) + end select sourceType + + broken = broken .or. any(IEEE_is_NaN(sourceState(phase)%p(i)%dotState(:,of))) + + enddo SourceLoop + +end function constitutive_damage_collectDotState + + +!-------------------------------------------------------------------------------------------------- +!> @brief contains the constitutive equation for calculating the rate of change of microstructure +!-------------------------------------------------------------------------------------------------- +function constitutive_thermal_collectDotState(S, ipc, ip, el,phase,of) result(broken) + + integer, intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el, & !< element + phase, & + of + real(pReal), intent(in), dimension(3,3) :: & + S !< 2nd Piola Kirchhoff stress (vector notation) + integer :: & + i !< counter in source loop + logical :: broken + + + broken = .false. + + SourceLoop: do i = 1, phase_Nsources(phase) + + sourceType: select case (phase_source(i,phase)) + case (SOURCE_thermal_externalheat_ID) sourceType call source_thermal_externalheat_dotState(phase,of) @@ -676,7 +707,7 @@ function constitutive_source_collectDotState(S, ipc, ip, el,phase,of) result(bro enddo SourceLoop -end function constitutive_source_collectDotState +end function constitutive_thermal_collectDotState !-------------------------------------------------------------------------------------------------- @@ -805,20 +836,18 @@ subroutine constitutive_results character(len=:), allocatable :: group - group = '/current/phase/' - call results_closeGroup(results_addGroup(group)) + call results_closeGroup(results_addGroup('/current/phase/')) do ph = 1, size(material_name_phase) - group = group//trim(material_name_phase(ph))//'/' + group = '/current/phase/'//trim(material_name_phase(ph))//'/' call results_closeGroup(results_addGroup(group)) call mech_results(group,ph) + call damage_results(group,ph) enddo - call damage_results - end subroutine constitutive_results @@ -1453,7 +1482,8 @@ subroutine integrateSourceState(g,i,e) p = material_phaseAt(g,e) c = material_phaseMemberAt(g,i,e) - broken = constitutive_source_collectDotState(crystallite_S(1:3,1:3,g,i,e), g,i,e,p,c) + broken = constitutive_thermal_collectDotState(crystallite_S(1:3,1:3,g,i,e), g,i,e,p,c) + broken = broken .or. constitutive_damage_collectDotState(crystallite_S(1:3,1:3,g,i,e), g,i,e,p,c) if(broken) return do s = 1, phase_Nsources(p) @@ -1471,7 +1501,8 @@ subroutine integrateSourceState(g,i,e) source_dotState(1:size_so(s),1,s) = sourceState(p)%p(s)%dotState(:,c) enddo - broken = constitutive_source_collectDotState(crystallite_S(1:3,1:3,g,i,e), g,i,e,p,c) + broken = constitutive_thermal_collectDotState(crystallite_S(1:3,1:3,g,i,e), g,i,e,p,c) + broken = broken .or. constitutive_damage_collectDotState(crystallite_S(1:3,1:3,g,i,e), g,i,e,p,c) if(broken) exit iteration do s = 1, phase_Nsources(p) diff --git a/src/constitutive_damage.f90 b/src/constitutive_damage.f90 index a864ca1b8..4e23b78ea 100644 --- a/src/constitutive_damage.f90 +++ b/src/constitutive_damage.f90 @@ -1,5 +1,5 @@ !---------------------------------------------------------------------------------------------------- -!> @brief internal microstructure state for all damage sources and kinematics constitutive models +!> @brief internal microstructure state for all damage sources and kinematics constitutive models !---------------------------------------------------------------------------------------------------- submodule(constitutive) constitutive_damage @@ -8,7 +8,7 @@ submodule(constitutive) constitutive_damage module function source_damage_anisoBrittle_init(source_length) result(mySources) integer, intent(in) :: source_length logical, dimension(:,:), allocatable :: mySources - end function source_damage_anisoBrittle_init + end function source_damage_anisoBrittle_init module function source_damage_anisoDuctile_init(source_length) result(mySources) integer, intent(in) :: source_length @@ -23,7 +23,7 @@ submodule(constitutive) constitutive_damage module function source_damage_isoDuctile_init(source_length) result(mySources) integer, intent(in) :: source_length logical, dimension(:,:), allocatable :: mySources - end function source_damage_isoDuctile_init + end function source_damage_isoDuctile_init module function kinematics_cleavage_opening_init(kinematics_length) result(myKinematics) integer, intent(in) :: kinematics_length @@ -39,14 +39,14 @@ submodule(constitutive) constitutive_damage 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 + constituent !< position of element within its phase instance real(pReal), intent(in) :: & - phi !< damage parameter + 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 @@ -129,7 +129,7 @@ module subroutine damage_init allocate(sourceState(ph)%p(phase_Nsources(ph))) enddo - allocate(phase_source(maxval(phase_Nsources),phases%length), source = SOURCE_undefined_ID) + allocate(phase_source(maxval(phase_Nsources),phases%length), source = SOURCE_undefined_ID) ! initialize source mechanisms if(maxval(phase_Nsources) /= 0) then @@ -141,19 +141,19 @@ module subroutine damage_init !-------------------------------------------------------------------------------------------------- ! initialize kinematic mechanisms - allocate(phase_Nkinematics(phases%length),source = 0) + allocate(phase_Nkinematics(phases%length),source = 0) do ph = 1,phases%length phase => phases%get(ph) kinematics => phase%get('kinematics',defaultVal=emptyList) phase_Nkinematics(ph) = kinematics%length enddo - - allocate(phase_kinematics(maxval(phase_Nkinematics),phases%length), source = KINEMATICS_undefined_ID) + + allocate(phase_kinematics(maxval(phase_Nkinematics),phases%length), source = KINEMATICS_undefined_ID) if(maxval(phase_Nkinematics) /= 0) then where(kinematics_cleavage_opening_init(maxval(phase_Nkinematics))) phase_kinematics = KINEMATICS_cleavage_opening_ID where(kinematics_slipplane_opening_init(maxval(phase_Nkinematics))) phase_kinematics = KINEMATICS_slipplane_opening_ID - endif + endif end subroutine damage_init @@ -167,7 +167,7 @@ module subroutine constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi ip, & !< integration point number el !< element number real(pReal), intent(in) :: & - phi !< damage parameter + phi !< damage parameter real(pReal), intent(inout) :: & phiDot, & dPhiDot_dPhi @@ -183,7 +183,7 @@ module subroutine constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi phiDot = 0.0_pReal dPhiDot_dPhi = 0.0_pReal - + do grain = 1, homogenization_Nconstituents(material_homogenizationAt(el)) phase = material_phaseAt(grain,el) constituent = material_phasememberAt(grain,ip,el) @@ -217,32 +217,35 @@ end subroutine constitutive_damage_getRateAndItsTangents !---------------------------------------------------------------------------------------------- !< @brief writes damage sources results to HDF5 output file !---------------------------------------------------------------------------------------------- -module subroutine damage_results +module subroutine damage_results(group,ph) - integer :: p,i - character(len=pStringLen) :: group + character(len=*), intent(in) :: group + integer, intent(in) :: ph - do p = 1, size(material_name_phase) + integer :: so - sourceLoop: do i = 1, phase_Nsources(p) - group = trim('current/phase')//'/'//trim(material_name_phase(p)) - group = trim(group)//'/sources' - call results_closeGroup(results_addGroup(group)) + sourceLoop: do so = 1, phase_Nsources(ph) - sourceType: select case (phase_source(i,p)) + if (phase_source(so,ph) /= SOURCE_UNDEFINED_ID) & + call results_closeGroup(results_addGroup(group//'damage/')) - 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 + sourceType: select case (phase_source(so,ph)) - enddo SourceLoop - enddo + case (SOURCE_damage_anisoBrittle_ID) sourceType + call source_damage_anisoBrittle_results(ph,group//'damage/') + + case (SOURCE_damage_anisoDuctile_ID) sourceType + call source_damage_anisoDuctile_results(ph,group//'damage/') + + case (SOURCE_damage_isoBrittle_ID) sourceType + call source_damage_isoBrittle_results(ph,group//'damage/') + + case (SOURCE_damage_isoDuctile_ID) sourceType + call source_damage_isoDuctile_results(ph,group//'damage/') + + end select sourceType + + enddo SourceLoop end subroutine damage_results