From fd4ae7127952d2c798df249f20c89874309b9ef7 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 30 Dec 2018 21:08:48 +0100 Subject: [PATCH] takeover from 40_XX and 41_XX branch easier to focus on damage instead of doing all kinematics and sources together --- src/constitutive.f90 | 14 ++- src/damage_local.f90 | 13 +- src/damage_nonlocal.f90 | 13 +- src/kinematics_cleavage_opening.f90 | 167 ++++++++------------------ src/kinematics_slipplane_opening.f90 | 170 ++++++++------------------- src/source_damage_anisoBrittle.f90 | 85 ++++++++------ src/source_damage_anisoDuctile.f90 | 85 ++++++++------ src/source_damage_isoBrittle.f90 | 70 ++++++----- src/source_damage_isoDuctile.f90 | 70 ++++++----- 9 files changed, 302 insertions(+), 385 deletions(-) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 8294047e7..ceb396823 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -186,8 +186,8 @@ subroutine constitutive_init() !-------------------------------------------------------------------------------------------------- ! parse kinematic mechanisms from config file call IO_checkAndRewind(FILEUNIT) - if (any(phase_kinematics == KINEMATICS_cleavage_opening_ID)) call kinematics_cleavage_opening_init(FILEUNIT) - if (any(phase_kinematics == KINEMATICS_slipplane_opening_ID)) call kinematics_slipplane_opening_init(FILEUNIT) + 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(FILEUNIT) if (any(phase_kinematics == KINEMATICS_vacancy_strain_ID)) call kinematics_vacancy_strain_init(FILEUNIT) if (any(phase_kinematics == KINEMATICS_hydrogen_strain_ID)) call kinematics_hydrogen_strain_init(FILEUNIT) @@ -1173,16 +1173,18 @@ function constitutive_postResults(S6, Fi, FeArray, ipc, ip, el) SourceLoop: do s = 1_pInt, phase_Nsources(material_phase(ipc,ip,el)) startPos = endPos + 1_pInt endPos = endPos + sourceState(material_phase(ipc,ip,el))%p(s)%sizePostResults + of = phasememberAt(ipc,ip,el) sourceType: select case (phase_source(s,material_phase(ipc,ip,el))) case (SOURCE_damage_isoBrittle_ID) sourceType - constitutive_postResults(startPos:endPos) = source_damage_isoBrittle_postResults(ipc, ip, el) + constitutive_postResults(startPos:endPos) = source_damage_isoBrittle_postResults(material_phase(ipc,ip,el),of) case (SOURCE_damage_isoDuctile_ID) sourceType - constitutive_postResults(startPos:endPos) = source_damage_isoDuctile_postResults(ipc, ip, el) + constitutive_postResults(startPos:endPos) = source_damage_isoDuctile_postResults(material_phase(ipc,ip,el),of) case (SOURCE_damage_anisoBrittle_ID) sourceType - constitutive_postResults(startPos:endPos) = source_damage_anisoBrittle_postResults(ipc, ip, el) + constitutive_postResults(startPos:endPos) = source_damage_anisoBrittle_postResults(material_phase(ipc,ip,el),of) case (SOURCE_damage_anisoDuctile_ID) sourceType - constitutive_postResults(startPos:endPos) = source_damage_anisoDuctile_postResults(ipc, ip, el) + constitutive_postResults(startPos:endPos) = source_damage_anisoDuctile_postResults(material_phase(ipc,ip,el),of) end select sourceType + enddo SourceLoop end function constitutive_postResults diff --git a/src/damage_local.f90 b/src/damage_local.f90 index 74bcb00db..6569347c2 100644 --- a/src/damage_local.f90 +++ b/src/damage_local.f90 @@ -225,6 +225,7 @@ subroutine damage_local_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el homogenization_Ngrains, & mappingHomogenization, & phaseAt, & + phasememberAt, & phase_source, & phase_Nsources, & SOURCE_damage_isoBrittle_ID, & @@ -249,7 +250,8 @@ subroutine damage_local_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el integer(pInt) :: & phase, & grain, & - source + source, & + constituent real(pReal) :: & phiDot, dPhiDot_dPhi, localphiDot, dLocalphiDot_dPhi @@ -257,19 +259,20 @@ subroutine damage_local_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el dPhiDot_dPhi = 0.0_pReal do grain = 1, homogenization_Ngrains(mappingHomogenization(2,ip,el)) phase = phaseAt(grain,ip,el) + constituent = 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, grain, ip, el) + 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, grain, ip, el) + 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, grain, ip, el) + 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, grain, ip, el) + call source_damage_anisoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) case default localphiDot = 0.0_pReal diff --git a/src/damage_nonlocal.f90 b/src/damage_nonlocal.f90 index 6b9093ef1..eab808266 100644 --- a/src/damage_nonlocal.f90 +++ b/src/damage_nonlocal.f90 @@ -186,6 +186,7 @@ subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, homogenization_Ngrains, & mappingHomogenization, & phaseAt, & + phasememberAt, & phase_source, & phase_Nsources, & SOURCE_damage_isoBrittle_ID, & @@ -210,7 +211,8 @@ subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, integer(pInt) :: & phase, & grain, & - source + source, & + constituent real(pReal) :: & phiDot, dPhiDot_dPhi, localphiDot, dLocalphiDot_dPhi @@ -218,19 +220,20 @@ subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, dPhiDot_dPhi = 0.0_pReal do grain = 1, homogenization_Ngrains(mappingHomogenization(2,ip,el)) phase = phaseAt(grain,ip,el) + constituent = phasememberAt(grain,ip,el) do source = 1_pInt, phase_Nsources(phase) select case(phase_source(source,phase)) case (SOURCE_damage_isoBrittle_ID) - call source_damage_isobrittle_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, grain, ip, el) + 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, grain, ip, el) + 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, grain, ip, el) + 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, grain, ip, el) + call source_damage_anisoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) case default localphiDot = 0.0_pReal diff --git a/src/kinematics_cleavage_opening.f90 b/src/kinematics_cleavage_opening.f90 index 998b19562..89d9dcd68 100644 --- a/src/kinematics_cleavage_opening.f90 +++ b/src/kinematics_cleavage_opening.f90 @@ -11,20 +11,22 @@ module kinematics_cleavage_opening implicit none private - integer(pInt), dimension(:), allocatable, public, protected :: & - kinematics_cleavage_opening_sizePostResults, & !< cumulative size of post results - kinematics_cleavage_opening_offset, & !< which kinematics is my current damage mechanism? - kinematics_cleavage_opening_instance !< instance of damage kinematics mechanism + integer(pInt), dimension(:), allocatable, private :: kinematics_cleavage_opening_instance - integer(pInt), dimension(:,:), allocatable, target, public :: & - kinematics_cleavage_opening_sizePostResult !< size of each post result output - - character(len=64), dimension(:,:), allocatable, target, public :: & - kinematics_cleavage_opening_output !< name of each post result output - - integer(pInt), dimension(:), allocatable, target, public :: & - kinematics_cleavage_opening_Noutput !< number of outputs per instance of this damage + type, private :: tParameters !< container type for internal constitutive parameters + integer(pInt) :: & + totalNcleavage + integer(pInt), dimension(:), allocatable :: & + Ncleavage !< active number of cleavage systems per family + real(pReal) :: & + sdot0, & + n + real(pReal), dimension(:), allocatable :: & + critDip, & + critLoad + end type +! Begin Deprecated integer(pInt), dimension(:), allocatable, private :: & kinematics_cleavage_opening_totalNcleavage !< total number of cleavage systems @@ -38,6 +40,7 @@ module kinematics_cleavage_opening real(pReal), dimension(:,:), allocatable, private :: & kinematics_cleavage_opening_critDisp, & kinematics_cleavage_opening_critLoad +! End Deprecated public :: & kinematics_cleavage_opening_init, & @@ -50,7 +53,7 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -subroutine kinematics_cleavage_opening_init(fileUnit) +subroutine kinematics_cleavage_opening_init() #if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & @@ -60,41 +63,25 @@ subroutine kinematics_cleavage_opening_init(fileUnit) debug_level,& debug_constitutive,& debug_levelBasic + use config, only: & + config_phase use IO, only: & - IO_read, & - IO_lc, & - IO_getTag, & - IO_isBlank, & - IO_stringPos, & - IO_stringValue, & - IO_floatValue, & - IO_intValue, & IO_warning, & IO_error, & - IO_timeStamp, & - IO_EOF + IO_timeStamp use material, only: & phase_kinematics, & - phase_Nkinematics, & - phase_Noutput, & KINEMATICS_cleavage_opening_label, & KINEMATICS_cleavage_opening_ID - use config, only: & - material_Nphase, & - MATERIAL_partPhase use lattice, only: & lattice_maxNcleavageFamily, & lattice_NcleavageSystem implicit none - integer(pInt), intent(in) :: fileUnit + integer(pInt), allocatable, dimension(:) :: tempInt + real(pReal), allocatable, dimension(:) :: tempFloat - integer(pInt), allocatable, dimension(:) :: chunkPos - integer(pInt) :: maxNinstance,phase,instance,kinematics - integer(pInt) :: Nchunks_CleavageFamilies = 0_pInt, j - character(len=65536) :: & - tag = '', & - line = '' + integer(pInt) :: maxNinstance,p,instance,kinematics write(6,'(/,a)') ' <<<+- kinematics_'//KINEMATICS_cleavage_opening_LABEL//' init -+>>>' write(6,'(a15,a)') ' Current time: ',IO_timeStamp() @@ -106,21 +93,11 @@ subroutine kinematics_cleavage_opening_init(fileUnit) if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) & write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance - allocate(kinematics_cleavage_opening_offset(material_Nphase), source=0_pInt) - allocate(kinematics_cleavage_opening_instance(material_Nphase), source=0_pInt) - do phase = 1, material_Nphase - kinematics_cleavage_opening_instance(phase) = count(phase_kinematics(:,1:phase) == kinematics_cleavage_opening_ID) - do kinematics = 1, phase_Nkinematics(phase) - if (phase_kinematics(kinematics,phase) == kinematics_cleavage_opening_ID) & - kinematics_cleavage_opening_offset(phase) = kinematics - enddo + allocate(kinematics_cleavage_opening_instance(size(config_phase)), source=0_pInt) + do p = 1_pInt, size(config_phase) + kinematics_cleavage_opening_instance(p) = count(phase_kinematics(:,1:p) == kinematics_cleavage_opening_ID) ! ToDo: count correct? enddo - allocate(kinematics_cleavage_opening_sizePostResults(maxNinstance), source=0_pInt) - allocate(kinematics_cleavage_opening_sizePostResult(maxval(phase_Noutput),maxNinstance), source=0_pInt) - allocate(kinematics_cleavage_opening_output(maxval(phase_Noutput),maxNinstance)) - kinematics_cleavage_opening_output = '' - allocate(kinematics_cleavage_opening_Noutput(maxNinstance), source=0_pInt) allocate(kinematics_cleavage_opening_critDisp(lattice_maxNcleavageFamily,maxNinstance), source=0.0_pReal) allocate(kinematics_cleavage_opening_critLoad(lattice_maxNcleavageFamily,maxNinstance), source=0.0_pReal) allocate(kinematics_cleavage_opening_Ncleavage(lattice_maxNcleavageFamily,maxNinstance), source=0_pInt) @@ -128,84 +105,44 @@ subroutine kinematics_cleavage_opening_init(fileUnit) allocate(kinematics_cleavage_opening_sdot_0(maxNinstance), source=0.0_pReal) allocate(kinematics_cleavage_opening_N(maxNinstance), source=0.0_pReal) - rewind(fileUnit) - phase = 0_pInt - do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase) ! wind forward to - line = IO_read(fileUnit) - enddo - - parsingFile: do while (trim(line) /= IO_EOF) ! read through sections of phase part - line = IO_read(fileUnit) - if (IO_isBlank(line)) cycle ! skip empty lines - if (IO_getTag(line,'<','>') /= '') then ! stop at next part - line = IO_read(fileUnit, .true.) ! reset IO_read - exit - endif - if (IO_getTag(line,'[',']') /= '') then ! next phase section - phase = phase + 1_pInt ! advance phase section counter - cycle ! skip to next line - endif - if (phase > 0_pInt ) then; if (any(phase_kinematics(:,phase) == KINEMATICS_cleavage_opening_ID)) then ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran - instance = kinematics_cleavage_opening_instance(phase) ! which instance of my damage is present phase - chunkPos = IO_stringPos(line) - tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key - select case(tag) - case ('anisobrittle_sdot0') - kinematics_cleavage_opening_sdot_0(instance) = IO_floatValue(line,chunkPos,2_pInt) - - case ('anisobrittle_ratesensitivity') - kinematics_cleavage_opening_N(instance) = IO_floatValue(line,chunkPos,2_pInt) - - case ('ncleavage') ! - Nchunks_CleavageFamilies = chunkPos(1) - 1_pInt - do j = 1_pInt, Nchunks_CleavageFamilies - kinematics_cleavage_opening_Ncleavage(j,instance) = IO_intValue(line,chunkPos,1_pInt+j) - enddo + do p = 1_pInt, size(config_phase) + if (all(phase_kinematics(:,p) /= KINEMATICS_cleavage_opening_ID)) cycle + instance = kinematics_cleavage_opening_instance(p) + kinematics_cleavage_opening_sdot_0(instance) = config_phase(p)%getFloat('anisobrittle_sdot0') + kinematics_cleavage_opening_N(instance) = config_phase(p)%getFloat('anisobrittle_ratesensitivity') + tempInt = config_phase(p)%getInts('ncleavage') + kinematics_cleavage_opening_Ncleavage(1:size(tempInt),instance) = tempInt - case ('anisobrittle_criticaldisplacement') - do j = 1_pInt, Nchunks_CleavageFamilies - kinematics_cleavage_opening_critDisp(j,instance) = IO_floatValue(line,chunkPos,1_pInt+j) - enddo + tempFloat = config_phase(p)%getFloats('anisobrittle_criticaldisplacement',requiredShape=shape(tempInt)) + kinematics_cleavage_opening_critDisp(1:size(tempInt),instance) = tempFloat - case ('anisobrittle_criticalload') - do j = 1_pInt, Nchunks_CleavageFamilies - kinematics_cleavage_opening_critLoad(j,instance) = IO_floatValue(line,chunkPos,1_pInt+j) - enddo + tempFloat = config_phase(p)%getFloats('anisobrittle_criticalload',requiredShape=shape(tempInt)) + kinematics_cleavage_opening_critLoad(1:size(tempInt),instance) = tempFloat - end select - endif; endif - enddo parsingFile - -!-------------------------------------------------------------------------------------------------- -! sanity checks - sanityChecks: do phase = 1_pInt, material_Nphase - myPhase: if (any(phase_kinematics(:,phase) == KINEMATICS_cleavage_opening_ID)) then - instance = kinematics_cleavage_opening_instance(phase) - kinematics_cleavage_opening_Ncleavage(1:lattice_maxNcleavageFamily,instance) = & - min(lattice_NcleavageSystem(1:lattice_maxNcleavageFamily,phase),& ! limit active cleavage systems per family to min of available and requested + kinematics_cleavage_opening_Ncleavage(1:lattice_maxNcleavageFamily,instance) = & + min(lattice_NcleavageSystem(1:lattice_maxNcleavageFamily,p),& ! limit active cleavage systems per family to min of available and requested kinematics_cleavage_opening_Ncleavage(1:lattice_maxNcleavageFamily,instance)) kinematics_cleavage_opening_totalNcleavage(instance) = sum(kinematics_cleavage_opening_Ncleavage(:,instance)) ! how many cleavage systems altogether if (kinematics_cleavage_opening_sdot_0(instance) <= 0.0_pReal) & call IO_error(211_pInt,el=instance,ext_msg='sdot_0 ('//KINEMATICS_cleavage_opening_LABEL//')') - if (any(kinematics_cleavage_opening_critDisp(1:Nchunks_CleavageFamilies,instance) < 0.0_pReal)) & + if (any(kinematics_cleavage_opening_critDisp(1:size(tempInt),instance) < 0.0_pReal)) & call IO_error(211_pInt,el=instance,ext_msg='critical_displacement ('//KINEMATICS_cleavage_opening_LABEL//')') - if (any(kinematics_cleavage_opening_critLoad(1:Nchunks_CleavageFamilies,instance) < 0.0_pReal)) & + if (any(kinematics_cleavage_opening_critLoad(1:size(tempInt),instance) < 0.0_pReal)) & call IO_error(211_pInt,el=instance,ext_msg='critical_load ('//KINEMATICS_cleavage_opening_LABEL//')') if (kinematics_cleavage_opening_N(instance) <= 0.0_pReal) & call IO_error(211_pInt,el=instance,ext_msg='rate_sensitivity ('//KINEMATICS_cleavage_opening_LABEL//')') - endif myPhase - enddo sanityChecks + enddo end subroutine kinematics_cleavage_opening_init !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the velocity gradient !-------------------------------------------------------------------------------------------------- -subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar3333, Tstar_v, ipc, ip, el) +subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, Tstar_v, ipc, ip, el) use prec, only: & tol_math_check use material, only: & - phaseAt, phasememberAt, & + material_phase, & material_homog, & damage, & damageMapping @@ -225,25 +162,22 @@ subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar3333, Tstar real(pReal), intent(out), dimension(3,3) :: & Ld !< damage velocity gradient real(pReal), intent(out), dimension(3,3,3,3) :: & - dLd_dTstar3333 !< derivative of Ld with respect to Tstar (4th-order tensor) + dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor) integer(pInt) :: & - phase, & - constituent, & - instance, & + instance, phase, & homog, damageOffset, & f, i, index_myFamily, k, l, m, n real(pReal) :: & traction_d, traction_t, traction_n, traction_crit, & udotd, dudotd_dt, udott, dudott_dt, udotn, dudotn_dt - phase = phaseAt(ipc,ip,el) - constituent = phasememberAt(ipc,ip,el) + phase = material_phase(ipc,ip,el) instance = kinematics_cleavage_opening_instance(phase) homog = material_homog(ip,el) damageOffset = damageMapping(homog)%p(ip,el) Ld = 0.0_pReal - dLd_dTstar3333 = 0.0_pReal + dLd_dTstar = 0.0_pReal do f = 1_pInt,lattice_maxNcleavageFamily index_myFamily = sum(lattice_NcleavageSystem(1:f-1_pInt,phase)) ! at which index starts my family do i = 1_pInt,kinematics_cleavage_opening_Ncleavage(f,instance) ! process each (active) cleavage system in family @@ -261,7 +195,7 @@ subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar3333, Tstar dudotd_dt = sign(1.0_pReal,traction_d)*udotd*kinematics_cleavage_opening_N(instance)/ & max(0.0_pReal, abs(traction_d) - traction_crit) forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & - dLd_dTstar3333(k,l,m,n) = dLd_dTstar3333(k,l,m,n) + & + dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + & dudotd_dt*lattice_Scleavage(k,l,1,index_myFamily+i,phase)* & lattice_Scleavage(m,n,1,index_myFamily+i,phase) endif @@ -275,7 +209,7 @@ subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar3333, Tstar dudott_dt = sign(1.0_pReal,traction_t)*udott*kinematics_cleavage_opening_N(instance)/ & max(0.0_pReal, abs(traction_t) - traction_crit) forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & - dLd_dTstar3333(k,l,m,n) = dLd_dTstar3333(k,l,m,n) + & + dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + & dudott_dt*lattice_Scleavage(k,l,2,index_myFamily+i,phase)* & lattice_Scleavage(m,n,2,index_myFamily+i,phase) endif @@ -289,11 +223,10 @@ subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar3333, Tstar dudotn_dt = sign(1.0_pReal,traction_n)*udotn*kinematics_cleavage_opening_N(instance)/ & max(0.0_pReal, abs(traction_n) - traction_crit) forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & - dLd_dTstar3333(k,l,m,n) = dLd_dTstar3333(k,l,m,n) + & + dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + & dudotn_dt*lattice_Scleavage(k,l,3,index_myFamily+i,phase)* & lattice_Scleavage(m,n,3,index_myFamily+i,phase) endif - enddo enddo diff --git a/src/kinematics_slipplane_opening.f90 b/src/kinematics_slipplane_opening.f90 index 61ff84b9f..573fe7d78 100644 --- a/src/kinematics_slipplane_opening.f90 +++ b/src/kinematics_slipplane_opening.f90 @@ -11,20 +11,22 @@ module kinematics_slipplane_opening implicit none private - integer(pInt), dimension(:), allocatable, public, protected :: & - kinematics_slipplane_opening_sizePostResults, & !< cumulative size of post results - kinematics_slipplane_opening_offset, & !< which kinematics is my current damage mechanism? - kinematics_slipplane_opening_instance !< instance of damage kinematics mechanism - - integer(pInt), dimension(:,:), allocatable, target, public :: & - kinematics_slipplane_opening_sizePostResult !< size of each post result output - - character(len=64), dimension(:,:), allocatable, target, public :: & - kinematics_slipplane_opening_output !< name of each post result output - - integer(pInt), dimension(:), allocatable, target, public :: & - kinematics_slipplane_opening_Noutput !< number of outputs per instance of this damage + integer(pInt), dimension(:), allocatable, private :: kinematics_slipplane_opening_instance + type, private :: tParameters !< container type for internal constitutive parameters + integer(pInt) :: & + totalNslip + integer(pInt), dimension(:), allocatable :: & + Nslip !< active number of slip systems per family + real(pReal) :: & + sdot0, & + n + real(pReal), dimension(:), allocatable :: & + critDip, & + critPlasticStrain + end type + +! Begin Deprecated integer(pInt), dimension(:), allocatable, private :: & kinematics_slipplane_opening_totalNslip !< total number of slip systems @@ -38,6 +40,7 @@ module kinematics_slipplane_opening real(pReal), dimension(:,:), allocatable, private :: & kinematics_slipplane_opening_critPlasticStrain, & kinematics_slipplane_opening_critLoad +! End Deprecated public :: & kinematics_slipplane_opening_init, & @@ -50,7 +53,7 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -subroutine kinematics_slipplane_opening_init(fileUnit) +subroutine kinematics_slipplane_opening_init() #if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & @@ -60,41 +63,25 @@ subroutine kinematics_slipplane_opening_init(fileUnit) debug_level,& debug_constitutive,& debug_levelBasic + use config, only: & + config_phase use IO, only: & - IO_read, & - IO_lc, & - IO_getTag, & - IO_isBlank, & - IO_stringPos, & - IO_stringValue, & - IO_floatValue, & - IO_intValue, & IO_warning, & IO_error, & - IO_timeStamp, & - IO_EOF + IO_timeStamp use material, only: & phase_kinematics, & - phase_Nkinematics, & - phase_Noutput, & KINEMATICS_slipplane_opening_label, & KINEMATICS_slipplane_opening_ID - use config, only: & - material_Nphase, & - MATERIAL_partPhase use lattice, only: & lattice_maxNslipFamily, & lattice_NslipSystem implicit none - integer(pInt), intent(in) :: fileUnit + integer(pInt), allocatable, dimension(:) :: tempInt + real(pReal), allocatable, dimension(:) :: tempFloat - integer(pInt), allocatable, dimension(:) :: chunkPos - integer(pInt) :: maxNinstance,phase,instance,kinematics - integer(pInt) :: Nchunks_SlipFamilies = 0_pInt, j - character(len=65536) :: & - tag = '', & - line = '' + integer(pInt) :: maxNinstance,p,instance,kinematics write(6,'(/,a)') ' <<<+- kinematics_'//KINEMATICS_slipplane_opening_LABEL//' init -+>>>' write(6,'(a15,a)') ' Current time: ',IO_timeStamp() @@ -106,21 +93,11 @@ subroutine kinematics_slipplane_opening_init(fileUnit) if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) & write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance - allocate(kinematics_slipplane_opening_offset(material_Nphase), source=0_pInt) - allocate(kinematics_slipplane_opening_instance(material_Nphase), source=0_pInt) - do phase = 1, material_Nphase - kinematics_slipplane_opening_instance(phase) = count(phase_kinematics(:,1:phase) == kinematics_slipplane_opening_ID) - do kinematics = 1, phase_Nkinematics(phase) - if (phase_kinematics(kinematics,phase) == kinematics_slipplane_opening_ID) & - kinematics_slipplane_opening_offset(phase) = kinematics - enddo + allocate(kinematics_slipplane_opening_instance(size(config_phase)), source=0_pInt) + do p = 1_pInt, size(config_phase) + kinematics_slipplane_opening_instance(p) = count(phase_kinematics(:,1:p) == kinematics_slipplane_opening_ID) ! ToDo: count correct? enddo - allocate(kinematics_slipplane_opening_sizePostResults(maxNinstance), source=0_pInt) - allocate(kinematics_slipplane_opening_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt) - allocate(kinematics_slipplane_opening_output(maxval(phase_Noutput),maxNinstance)) - kinematics_slipplane_opening_output = '' - allocate(kinematics_slipplane_opening_Noutput(maxNinstance), source=0_pInt) allocate(kinematics_slipplane_opening_critLoad(lattice_maxNslipFamily,maxNinstance), source=0.0_pReal) allocate(kinematics_slipplane_opening_critPlasticStrain(lattice_maxNslipFamily,maxNinstance),source=0.0_pReal) allocate(kinematics_slipplane_opening_Nslip(lattice_maxNslipFamily,maxNinstance), source=0_pInt) @@ -128,61 +105,22 @@ subroutine kinematics_slipplane_opening_init(fileUnit) allocate(kinematics_slipplane_opening_N(maxNinstance), source=0.0_pReal) allocate(kinematics_slipplane_opening_sdot_0(maxNinstance), source=0.0_pReal) - rewind(fileUnit) - phase = 0_pInt - do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase) ! wind forward to - line = IO_read(fileUnit) - enddo - - parsingFile: do while (trim(line) /= IO_EOF) ! read through sections of phase part - line = IO_read(fileUnit) - if (IO_isBlank(line)) cycle ! skip empty lines - if (IO_getTag(line,'<','>') /= '') then ! stop at next part - line = IO_read(fileUnit, .true.) ! reset IO_read - exit - endif - if (IO_getTag(line,'[',']') /= '') then ! next phase section - phase = phase + 1_pInt ! advance phase section counter - cycle ! skip to next line - endif - if (phase > 0_pInt ) then; if (any(phase_kinematics(:,phase) == KINEMATICS_slipplane_opening_ID)) then ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran - instance = kinematics_slipplane_opening_instance(phase) ! which instance of my damage is present phase - chunkPos = IO_stringPos(line) - tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key - select case(tag) - case ('nslip') ! - Nchunks_SlipFamilies = chunkPos(1) - 1_pInt - do j = 1_pInt, Nchunks_SlipFamilies - kinematics_slipplane_opening_Nslip(j,instance) = IO_intValue(line,chunkPos,1_pInt+j) - enddo + do p = 1_pInt, size(config_phase) + if (all(phase_kinematics(:,p) /= KINEMATICS_slipplane_opening_ID)) cycle + instance = kinematics_slipplane_opening_instance(p) + kinematics_slipplane_opening_sdot_0(instance) = config_phase(p)%getFloat('anisoductile_sdot0') + kinematics_slipplane_opening_N(instance) = config_phase(p)%getFloat('anisoductile_ratesensitivity') + tempInt = config_phase(p)%getInts('ncleavage') + kinematics_slipplane_opening_Nslip(1:size(tempInt),instance) = tempInt - case ('anisoductile_sdot0') - kinematics_slipplane_opening_sdot_0(instance) = IO_floatValue(line,chunkPos,2_pInt) - - case ('anisoductile_criticalplasticstrain') - do j = 1_pInt, Nchunks_SlipFamilies - kinematics_slipplane_opening_critPlasticStrain(j,instance) = IO_floatValue(line,chunkPos,1_pInt+j) - enddo - - case ('anisoductile_ratesensitivity') - kinematics_slipplane_opening_N(instance) = IO_floatValue(line,chunkPos,2_pInt) + tempFloat = config_phase(p)%getFloats('anisoductile_criticalplasticstrain',requiredShape=shape(tempInt)) + kinematics_slipplane_opening_critPlasticStrain(1:size(tempInt),instance) = tempFloat - case ('anisoductile_criticalload') - do j = 1_pInt, Nchunks_SlipFamilies - kinematics_slipplane_opening_critLoad(j,instance) = IO_floatValue(line,chunkPos,1_pInt+j) - enddo - - end select - endif; endif - enddo parsingFile + tempFloat = config_phase(p)%getFloats('anisoductile_criticalload',requiredShape=shape(tempInt)) + kinematics_slipplane_opening_critLoad(1:size(tempInt),instance) = tempFloat -!-------------------------------------------------------------------------------------------------- -! sanity checks - sanityChecks: do phase = 1_pInt, material_Nphase - myPhase: if (any(phase_kinematics(:,phase) == KINEMATICS_slipplane_opening_ID)) then - instance = kinematics_slipplane_opening_instance(phase) kinematics_slipplane_opening_Nslip(1:lattice_maxNslipFamily,instance) = & - min(lattice_NslipSystem(1:lattice_maxNslipFamily,phase),& ! limit active cleavage systems per family to min of available and requested + min(lattice_NslipSystem(1:lattice_maxNslipFamily,p),& ! limit active cleavage systems per family to min of available and requested kinematics_slipplane_opening_Nslip(1:lattice_maxNslipFamily,instance)) kinematics_slipplane_opening_totalNslip(instance) = sum(kinematics_slipplane_opening_Nslip(:,instance)) if (kinematics_slipplane_opening_sdot_0(instance) <= 0.0_pReal) & @@ -191,16 +129,14 @@ subroutine kinematics_slipplane_opening_init(fileUnit) call IO_error(211_pInt,el=instance,ext_msg='criticaPlasticStrain ('//KINEMATICS_slipplane_opening_LABEL//')') if (kinematics_slipplane_opening_N(instance) <= 0.0_pReal) & call IO_error(211_pInt,el=instance,ext_msg='rate_sensitivity ('//KINEMATICS_slipplane_opening_LABEL//')') - endif myPhase - enddo sanityChecks + enddo - end subroutine kinematics_slipplane_opening_init !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the velocity gradient !-------------------------------------------------------------------------------------------------- -subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar3333, Tstar_v, ipc, ip, el) +subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, Tstar_v, ipc, ip, el) use prec, only: & tol_math_check use lattice, only: & @@ -210,19 +146,15 @@ subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar3333, Tsta lattice_st, & lattice_sn use material, only: & - phaseAt, phasememberAt, & + material_phase, & material_homog, & damage, & damageMapping use math, only: & math_Plain3333to99, & - math_I3, & - math_identity4th, & math_symmetric33, & math_Mandel33to6, & - math_tensorproduct33, & - math_det33, & - math_mul33x33 + math_tensorproduct33 implicit none integer(pInt), intent(in) :: & @@ -234,29 +166,26 @@ subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar3333, Tsta real(pReal), intent(out), dimension(3,3) :: & Ld !< damage velocity gradient real(pReal), intent(out), dimension(3,3,3,3) :: & - dLd_dTstar3333 !< derivative of Ld with respect to Tstar (4th-order tensor) + dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor) real(pReal), dimension(3,3) :: & projection_d, projection_t, projection_n !< projection modes 3x3 tensor real(pReal), dimension(6) :: & projection_d_v, projection_t_v, projection_n_v !< projection modes 3x3 vector integer(pInt) :: & - phase, & - constituent, & - instance, & + instance, phase, & homog, damageOffset, & f, i, index_myFamily, k, l, m, n real(pReal) :: & traction_d, traction_t, traction_n, traction_crit, & udotd, dudotd_dt, udott, dudott_dt, udotn, dudotn_dt - phase = phaseAt(ipc,ip,el) - constituent = phasememberAt(ipc,ip,el) + phase = material_phase(ipc,ip,el) instance = kinematics_slipplane_opening_instance(phase) homog = material_homog(ip,el) damageOffset = damageMapping(homog)%p(ip,el) Ld = 0.0_pReal - dLd_dTstar3333 = 0.0_pReal + dLd_dTstar = 0.0_pReal do f = 1_pInt,lattice_maxNslipFamily index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,phase)) ! at which index starts my family do i = 1_pInt,kinematics_slipplane_opening_Nslip(f,instance) ! process each (active) slip system in family @@ -287,7 +216,7 @@ subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar3333, Tsta Ld = Ld + udotd*projection_d dudotd_dt = udotd*kinematics_slipplane_opening_N(instance)/traction_d forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & - dLd_dTstar3333(k,l,m,n) = dLd_dTstar3333(k,l,m,n) + & + dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + & dudotd_dt*projection_d(k,l)*projection_d(m,n) endif @@ -300,9 +229,10 @@ subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar3333, Tsta Ld = Ld + udott*projection_t dudott_dt = udott*kinematics_slipplane_opening_N(instance)/traction_t forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & - dLd_dTstar3333(k,l,m,n) = dLd_dTstar3333(k,l,m,n) + & + dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + & dudott_dt*projection_t(k,l)*projection_t(m,n) endif + udotn = & kinematics_slipplane_opening_sdot_0(instance)* & (max(0.0_pReal,traction_n)/traction_crit - & @@ -311,7 +241,7 @@ subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar3333, Tsta Ld = Ld + udotn*projection_n dudotn_dt = udotn*kinematics_slipplane_opening_N(instance)/traction_n forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & - dLd_dTstar3333(k,l,m,n) = dLd_dTstar3333(k,l,m,n) + & + dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + & dudotn_dt*projection_n(k,l)*projection_n(m,n) endif enddo diff --git a/src/source_damage_anisoBrittle.f90 b/src/source_damage_anisoBrittle.f90 index 6b222c37c..b8bd3246d 100644 --- a/src/source_damage_anisoBrittle.f90 +++ b/src/source_damage_anisoBrittle.f90 @@ -49,6 +49,23 @@ module source_damage_anisoBrittle source_damage_anisoBrittle_outputID !< ID of each post result output + type, private :: tParameters !< container type for internal constitutive parameters + real(pReal) :: & + aTol, & + sdot_0, & + N + real(pReal), dimension(:), allocatable :: & + critDisp, & + critLoad + integer(pInt) :: & + totalNcleavage + integer(pInt), dimension(:), allocatable :: & + Ncleavage + end type tParameters + + type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance) + + public :: & source_damage_anisoBrittle_init, & source_damage_anisoBrittle_dotState, & @@ -94,6 +111,7 @@ subroutine source_damage_anisoBrittle_init(fileUnit) material_phase, & sourceState use config, only: & + config_phase, & material_Nphase, & MATERIAL_partPhase use numerics,only: & @@ -106,9 +124,9 @@ subroutine source_damage_anisoBrittle_init(fileUnit) integer(pInt), intent(in) :: fileUnit integer(pInt), allocatable, dimension(:) :: chunkPos - integer(pInt) :: maxNinstance,mySize=0_pInt,phase,instance,source,sourceOffset,o + integer(pInt) :: Ninstance,mySize=0_pInt,phase,instance,source,sourceOffset,o integer(pInt) :: sizeState, sizeDotState, sizeDeltaState - integer(pInt) :: NofMyPhase + integer(pInt) :: NofMyPhase,p integer(pInt) :: Nchunks_CleavageFamilies = 0_pInt, j character(len=65536) :: & tag = '', & @@ -118,11 +136,11 @@ subroutine source_damage_anisoBrittle_init(fileUnit) write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" - maxNinstance = int(count(phase_source == SOURCE_damage_anisoBrittle_ID),pInt) - if (maxNinstance == 0_pInt) return + Ninstance = int(count(phase_source == SOURCE_damage_anisoBrittle_ID),pInt) + if (Ninstance == 0_pInt) return if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) & - write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance + write(6,'(a16,1x,i5,/)') '# instances:',Ninstance allocate(source_damage_anisoBrittle_offset(material_Nphase), source=0_pInt) allocate(source_damage_anisoBrittle_instance(material_Nphase), source=0_pInt) @@ -134,19 +152,24 @@ subroutine source_damage_anisoBrittle_init(fileUnit) enddo enddo - allocate(source_damage_anisoBrittle_sizePostResults(maxNinstance), source=0_pInt) - allocate(source_damage_anisoBrittle_sizePostResult(maxval(phase_Noutput),maxNinstance), source=0_pInt) - allocate(source_damage_anisoBrittle_output(maxval(phase_Noutput),maxNinstance)) + allocate(source_damage_anisoBrittle_sizePostResults(Ninstance), source=0_pInt) + allocate(source_damage_anisoBrittle_sizePostResult(maxval(phase_Noutput),Ninstance), source=0_pInt) + allocate(source_damage_anisoBrittle_output(maxval(phase_Noutput),Ninstance)) source_damage_anisoBrittle_output = '' - allocate(source_damage_anisoBrittle_outputID(maxval(phase_Noutput),maxNinstance), source=undefined_ID) - allocate(source_damage_anisoBrittle_Noutput(maxNinstance), source=0_pInt) - allocate(source_damage_anisoBrittle_critDisp(lattice_maxNcleavageFamily,maxNinstance), source=0.0_pReal) - allocate(source_damage_anisoBrittle_critLoad(lattice_maxNcleavageFamily,maxNinstance), source=0.0_pReal) - allocate(source_damage_anisoBrittle_Ncleavage(lattice_maxNcleavageFamily,maxNinstance), source=0_pInt) - allocate(source_damage_anisoBrittle_totalNcleavage(maxNinstance), source=0_pInt) - allocate(source_damage_anisoBrittle_aTol(maxNinstance), source=0.0_pReal) - allocate(source_damage_anisoBrittle_sdot_0(maxNinstance), source=0.0_pReal) - allocate(source_damage_anisoBrittle_N(maxNinstance), source=0.0_pReal) + allocate(source_damage_anisoBrittle_outputID(maxval(phase_Noutput),Ninstance), source=undefined_ID) + allocate(source_damage_anisoBrittle_Noutput(Ninstance), source=0_pInt) + + allocate(source_damage_anisoBrittle_critDisp(lattice_maxNcleavageFamily,Ninstance), source=0.0_pReal) + allocate(source_damage_anisoBrittle_critLoad(lattice_maxNcleavageFamily,Ninstance), source=0.0_pReal) + allocate(source_damage_anisoBrittle_Ncleavage(lattice_maxNcleavageFamily,Ninstance), source=0_pInt) + allocate(source_damage_anisoBrittle_totalNcleavage(Ninstance), source=0_pInt) + allocate(source_damage_anisoBrittle_aTol(Ninstance), source=0.0_pReal) + allocate(source_damage_anisoBrittle_sdot_0(Ninstance), source=0.0_pReal) + allocate(source_damage_anisoBrittle_N(Ninstance), source=0.0_pReal) + + do p=1, size(config_phase) + if (all(phase_source(:,p) /= SOURCE_damage_anisoBrittle_ID)) cycle + enddo rewind(fileUnit) phase = 0_pInt @@ -349,26 +372,22 @@ end subroutine source_damage_anisoBrittle_dotState !-------------------------------------------------------------------------------------------------- !> @brief returns local part of nonlocal damage driving force !-------------------------------------------------------------------------------------------------- -subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ipc, ip, el) +subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) use material, only: & - phaseAt, phasememberAt, & sourceState implicit none integer(pInt), intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element + phase, & + constituent real(pReal), intent(in) :: & phi real(pReal), intent(out) :: & localphiDot, & dLocalphiDot_dPhi integer(pInt) :: & - phase, constituent, sourceOffset + sourceOffset - phase = phaseAt(ipc,ip,el) - constituent = phasememberAt(ipc,ip,el) sourceOffset = source_damage_anisoBrittle_offset(phase) localphiDot = 1.0_pReal - & @@ -381,25 +400,21 @@ end subroutine source_damage_anisobrittle_getRateAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief return array of local damage results !-------------------------------------------------------------------------------------------------- -function source_damage_anisoBrittle_postResults(ipc,ip,el) +function source_damage_anisoBrittle_postResults(phase, constituent) use material, only: & - phaseAt, phasememberAt, & sourceState implicit none - integer(pInt), intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element + integer(pInt), intent(in) :: & + phase, & + constituent real(pReal), dimension(source_damage_anisoBrittle_sizePostResults( & - source_damage_anisoBrittle_instance(phaseAt(ipc,ip,el)))) :: & + source_damage_anisoBrittle_instance(phase))) :: & source_damage_anisoBrittle_postResults integer(pInt) :: & - instance, phase, constituent, sourceOffset, o, c + instance, sourceOffset, o, c - phase = phaseAt(ipc,ip,el) - constituent = phasememberAt(ipc,ip,el) instance = source_damage_anisoBrittle_instance(phase) sourceOffset = source_damage_anisoBrittle_offset(phase) diff --git a/src/source_damage_anisoDuctile.f90 b/src/source_damage_anisoDuctile.f90 index 5978960fb..c52dd4ff4 100644 --- a/src/source_damage_anisoDuctile.f90 +++ b/src/source_damage_anisoDuctile.f90 @@ -53,6 +53,23 @@ module source_damage_anisoDuctile source_damage_anisoDuctile_outputID !< ID of each post result output + type, private :: tParameters !< container type for internal constitutive parameters + real(pReal) :: & + aTol, & + sdot_0, & + N + real(pReal), dimension(:), allocatable :: & + critPlasticStrain, & + critLoad + integer(pInt) :: & + totalNslip + integer(pInt), dimension(:), allocatable :: & + Nslip + end type tParameters + + type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance) + + public :: & source_damage_anisoDuctile_init, & source_damage_anisoDuctile_dotState, & @@ -98,6 +115,7 @@ subroutine source_damage_anisoDuctile_init(fileUnit) material_phase, & sourceState use config, only: & + config_phase, & material_Nphase, & MATERIAL_partPhase use numerics,only: & @@ -110,9 +128,9 @@ subroutine source_damage_anisoDuctile_init(fileUnit) integer(pInt), intent(in) :: fileUnit integer(pInt), allocatable, dimension(:) :: chunkPos - integer(pInt) :: maxNinstance,mySize=0_pInt,phase,instance,source,sourceOffset,o + integer(pInt) :: Ninstance,mySize=0_pInt,phase,instance,source,sourceOffset,o integer(pInt) :: sizeState, sizeDotState, sizeDeltaState - integer(pInt) :: NofMyPhase + integer(pInt) :: NofMyPhase,p integer(pInt) :: Nchunks_SlipFamilies = 0_pInt, j character(len=65536) :: & tag = '', & @@ -122,11 +140,11 @@ subroutine source_damage_anisoDuctile_init(fileUnit) write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" - maxNinstance = int(count(phase_source == SOURCE_damage_anisoDuctile_ID),pInt) - if (maxNinstance == 0_pInt) return + Ninstance = int(count(phase_source == SOURCE_damage_anisoDuctile_ID),pInt) + if (Ninstance == 0_pInt) return if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) & - write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance + write(6,'(a16,1x,i5,/)') '# instances:',Ninstance allocate(source_damage_anisoDuctile_offset(material_Nphase), source=0_pInt) allocate(source_damage_anisoDuctile_instance(material_Nphase), source=0_pInt) @@ -138,19 +156,24 @@ subroutine source_damage_anisoDuctile_init(fileUnit) enddo enddo - allocate(source_damage_anisoDuctile_sizePostResults(maxNinstance), source=0_pInt) - allocate(source_damage_anisoDuctile_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt) - allocate(source_damage_anisoDuctile_output(maxval(phase_Noutput),maxNinstance)) + allocate(source_damage_anisoDuctile_sizePostResults(Ninstance), source=0_pInt) + allocate(source_damage_anisoDuctile_sizePostResult(maxval(phase_Noutput),Ninstance),source=0_pInt) + allocate(source_damage_anisoDuctile_output(maxval(phase_Noutput),Ninstance)) source_damage_anisoDuctile_output = '' - allocate(source_damage_anisoDuctile_outputID(maxval(phase_Noutput),maxNinstance), source=undefined_ID) - allocate(source_damage_anisoDuctile_Noutput(maxNinstance), source=0_pInt) - allocate(source_damage_anisoDuctile_critLoad(lattice_maxNslipFamily,maxNinstance), source=0.0_pReal) - allocate(source_damage_anisoDuctile_critPlasticStrain(lattice_maxNslipFamily,maxNinstance),source=0.0_pReal) - allocate(source_damage_anisoDuctile_Nslip(lattice_maxNslipFamily,maxNinstance), source=0_pInt) - allocate(source_damage_anisoDuctile_totalNslip(maxNinstance), source=0_pInt) - allocate(source_damage_anisoDuctile_N(maxNinstance), source=0.0_pReal) - allocate(source_damage_anisoDuctile_sdot_0(maxNinstance), source=0.0_pReal) - allocate(source_damage_anisoDuctile_aTol(maxNinstance), source=0.0_pReal) + allocate(source_damage_anisoDuctile_outputID(maxval(phase_Noutput),Ninstance), source=undefined_ID) + allocate(source_damage_anisoDuctile_Noutput(Ninstance), source=0_pInt) + + allocate(source_damage_anisoDuctile_critLoad(lattice_maxNslipFamily,Ninstance), source=0.0_pReal) + allocate(source_damage_anisoDuctile_critPlasticStrain(lattice_maxNslipFamily,Ninstance),source=0.0_pReal) + allocate(source_damage_anisoDuctile_Nslip(lattice_maxNslipFamily,Ninstance), source=0_pInt) + allocate(source_damage_anisoDuctile_totalNslip(Ninstance), source=0_pInt) + allocate(source_damage_anisoDuctile_N(Ninstance), source=0.0_pReal) + allocate(source_damage_anisoDuctile_sdot_0(Ninstance), source=0.0_pReal) + allocate(source_damage_anisoDuctile_aTol(Ninstance), source=0.0_pReal) + + do p=1, size(config_phase) + if (all(phase_source(:,p) /= SOURCE_damage_anisoDuctile_ID)) cycle + enddo rewind(fileUnit) phase = 0_pInt @@ -338,26 +361,22 @@ end subroutine source_damage_anisoDuctile_dotState !-------------------------------------------------------------------------------------------------- !> @brief returns local part of nonlocal damage driving force !-------------------------------------------------------------------------------------------------- -subroutine source_damage_anisoDuctile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ipc, ip, el) +subroutine source_damage_anisoDuctile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) use material, only: & - phaseAt, phasememberAt, & sourceState implicit none integer(pInt), intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element + phase, & + constituent real(pReal), intent(in) :: & phi real(pReal), intent(out) :: & localphiDot, & dLocalphiDot_dPhi integer(pInt) :: & - phase, constituent, sourceOffset + sourceOffset - phase = phaseAt(ipc,ip,el) - constituent = phasememberAt(ipc,ip,el) sourceOffset = source_damage_anisoDuctile_offset(phase) localphiDot = 1.0_pReal - & @@ -371,25 +390,21 @@ end subroutine source_damage_anisoDuctile_getRateAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief return array of local damage results !-------------------------------------------------------------------------------------------------- -function source_damage_anisoDuctile_postResults(ipc,ip,el) +function source_damage_anisoDuctile_postResults(phase, constituent) use material, only: & - phaseAt, phasememberAt, & sourceState implicit none - integer(pInt), intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element + integer(pInt), intent(in) :: & + phase, & + constituent real(pReal), dimension(source_damage_anisoDuctile_sizePostResults( & - source_damage_anisoDuctile_instance(phaseAt(ipc,ip,el)))) :: & + source_damage_anisoDuctile_instance(phase))) :: & source_damage_anisoDuctile_postResults integer(pInt) :: & - instance, phase, constituent, sourceOffset, o, c + instance, sourceOffset, o, c - phase = phaseAt(ipc,ip,el) - constituent = phasememberAt(ipc,ip,el) instance = source_damage_anisoDuctile_instance(phase) sourceOffset = source_damage_anisoDuctile_offset(phase) diff --git a/src/source_damage_isoBrittle.f90 b/src/source_damage_isoBrittle.f90 index 041761afe..6f572c72b 100644 --- a/src/source_damage_isoBrittle.f90 +++ b/src/source_damage_isoBrittle.f90 @@ -39,6 +39,16 @@ module source_damage_isoBrittle source_damage_isoBrittle_outputID !< ID of each post result output + type, private :: tParameters !< container type for internal constitutive parameters + real(pReal) :: & + critStrainEnergy, & + N, & + aTol + end type tParameters + + type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance) + + public :: & source_damage_isoBrittle_init, & source_damage_isoBrittle_deltaState, & @@ -84,6 +94,7 @@ subroutine source_damage_isoBrittle_init(fileUnit) material_phase, & sourceState use config, only: & + config_phase, & material_Nphase, & MATERIAL_partPhase use numerics,only: & @@ -93,9 +104,9 @@ subroutine source_damage_isoBrittle_init(fileUnit) integer(pInt), intent(in) :: fileUnit integer(pInt), allocatable, dimension(:) :: chunkPos - integer(pInt) :: maxNinstance,mySize=0_pInt,phase,instance,source,sourceOffset,o + integer(pInt) :: Ninstance,mySize=0_pInt,phase,instance,source,sourceOffset,o integer(pInt) :: sizeState, sizeDotState, sizeDeltaState - integer(pInt) :: NofMyPhase + integer(pInt) :: NofMyPhase,p character(len=65536) :: & tag = '', & line = '' @@ -104,11 +115,11 @@ subroutine source_damage_isoBrittle_init(fileUnit) write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" - maxNinstance = int(count(phase_source == SOURCE_damage_isoBrittle_ID),pInt) - if (maxNinstance == 0_pInt) return + Ninstance = int(count(phase_source == SOURCE_damage_isoBrittle_ID),pInt) + if (Ninstance == 0_pInt) return if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) & - write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance + write(6,'(a16,1x,i5,/)') '# instances:',Ninstance allocate(source_damage_isoBrittle_offset(material_Nphase), source=0_pInt) allocate(source_damage_isoBrittle_instance(material_Nphase), source=0_pInt) @@ -120,15 +131,20 @@ subroutine source_damage_isoBrittle_init(fileUnit) enddo enddo - allocate(source_damage_isoBrittle_sizePostResults(maxNinstance), source=0_pInt) - allocate(source_damage_isoBrittle_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt) - allocate(source_damage_isoBrittle_output(maxval(phase_Noutput),maxNinstance)) + allocate(source_damage_isoBrittle_sizePostResults(Ninstance), source=0_pInt) + allocate(source_damage_isoBrittle_sizePostResult(maxval(phase_Noutput),Ninstance),source=0_pInt) + allocate(source_damage_isoBrittle_output(maxval(phase_Noutput),Ninstance)) source_damage_isoBrittle_output = '' - allocate(source_damage_isoBrittle_outputID(maxval(phase_Noutput),maxNinstance), source=undefined_ID) - allocate(source_damage_isoBrittle_Noutput(maxNinstance), source=0_pInt) - allocate(source_damage_isoBrittle_critStrainEnergy(maxNinstance), source=0.0_pReal) - allocate(source_damage_isoBrittle_N(maxNinstance), source=1.0_pReal) - allocate(source_damage_isoBrittle_aTol(maxNinstance), source=0.0_pReal) + allocate(source_damage_isoBrittle_outputID(maxval(phase_Noutput),Ninstance), source=undefined_ID) + allocate(source_damage_isoBrittle_Noutput(Ninstance), source=0_pInt) + + allocate(source_damage_isoBrittle_critStrainEnergy(Ninstance), source=0.0_pReal) + allocate(source_damage_isoBrittle_N(Ninstance), source=1.0_pReal) + allocate(source_damage_isoBrittle_aTol(Ninstance), source=0.0_pReal) + + do p=1, size(config_phase) + if (all(phase_source(:,p) /= SOURCE_damage_isoBrittle_ID)) cycle + enddo rewind(fileUnit) phase = 0_pInt @@ -306,26 +322,22 @@ end subroutine source_damage_isoBrittle_deltaState !-------------------------------------------------------------------------------------------------- !> @brief returns local part of nonlocal damage driving force !-------------------------------------------------------------------------------------------------- -subroutine source_damage_isoBrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ipc, ip, el) +subroutine source_damage_isoBrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) use material, only: & - phaseAt, phasememberAt, & sourceState implicit none integer(pInt), intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element + phase, & + constituent real(pReal), intent(in) :: & phi real(pReal), intent(out) :: & localphiDot, & dLocalphiDot_dPhi integer(pInt) :: & - phase, constituent, instance, sourceOffset + instance, sourceOffset - phase = phaseAt(ipc,ip,el) - constituent = phasememberAt(ipc,ip,el) instance = source_damage_isoBrittle_instance(phase) sourceOffset = source_damage_isoBrittle_offset(phase) @@ -340,25 +352,21 @@ end subroutine source_damage_isoBrittle_getRateAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief return array of local damage results !-------------------------------------------------------------------------------------------------- -function source_damage_isoBrittle_postResults(ipc,ip,el) +function source_damage_isoBrittle_postResults(phase, constituent) use material, only: & - phaseAt, phasememberAt, & sourceState implicit none - integer(pInt), intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element + integer(pInt), intent(in) :: & + phase, & + constituent real(pReal), dimension(source_damage_isoBrittle_sizePostResults( & - source_damage_isoBrittle_instance(phaseAt(ipc,ip,el)))) :: & + source_damage_isoBrittle_instance(phase))) :: & source_damage_isoBrittle_postResults integer(pInt) :: & - instance, phase, constituent, sourceOffset, o, c + instance, sourceOffset, o, c - phase = phaseAt(ipc,ip,el) - constituent = phasememberAt(ipc,ip,el) instance = source_damage_isoBrittle_instance(phase) sourceOffset = source_damage_isoBrittle_offset(phase) diff --git a/src/source_damage_isoDuctile.f90 b/src/source_damage_isoDuctile.f90 index e843be728..b4ecb53e4 100644 --- a/src/source_damage_isoDuctile.f90 +++ b/src/source_damage_isoDuctile.f90 @@ -39,6 +39,16 @@ module source_damage_isoDuctile source_damage_isoDuctile_outputID !< ID of each post result output + type, private :: tParameters !< container type for internal constitutive parameters + real(pReal) :: & + critPlasticStrain, & + N, & + aTol + end type tParameters + + type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance) + + public :: & source_damage_isoDuctile_init, & source_damage_isoDuctile_dotState, & @@ -84,6 +94,7 @@ subroutine source_damage_isoDuctile_init(fileUnit) material_phase, & sourceState use config, only: & + config_phase, & material_Nphase, & MATERIAL_partPhase @@ -94,9 +105,9 @@ subroutine source_damage_isoDuctile_init(fileUnit) integer(pInt), intent(in) :: fileUnit integer(pInt), allocatable, dimension(:) :: chunkPos - integer(pInt) :: maxNinstance,mySize=0_pInt,phase,instance,source,sourceOffset,o + integer(pInt) :: Ninstance,mySize=0_pInt,phase,instance,source,sourceOffset,o integer(pInt) :: sizeState, sizeDotState, sizeDeltaState - integer(pInt) :: NofMyPhase + integer(pInt) :: NofMyPhase,p character(len=65536) :: & tag = '', & line = '' @@ -105,11 +116,11 @@ subroutine source_damage_isoDuctile_init(fileUnit) write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" - maxNinstance = int(count(phase_source == SOURCE_damage_isoDuctile_ID),pInt) - if (maxNinstance == 0_pInt) return + Ninstance = int(count(phase_source == SOURCE_damage_isoDuctile_ID),pInt) + if (Ninstance == 0_pInt) return if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) & - write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance + write(6,'(a16,1x,i5,/)') '# instances:',Ninstance allocate(source_damage_isoDuctile_offset(material_Nphase), source=0_pInt) allocate(source_damage_isoDuctile_instance(material_Nphase), source=0_pInt) @@ -121,15 +132,20 @@ subroutine source_damage_isoDuctile_init(fileUnit) enddo enddo - allocate(source_damage_isoDuctile_sizePostResults(maxNinstance), source=0_pInt) - allocate(source_damage_isoDuctile_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt) - allocate(source_damage_isoDuctile_output(maxval(phase_Noutput),maxNinstance)) + allocate(source_damage_isoDuctile_sizePostResults(Ninstance), source=0_pInt) + allocate(source_damage_isoDuctile_sizePostResult(maxval(phase_Noutput),Ninstance),source=0_pInt) + allocate(source_damage_isoDuctile_output(maxval(phase_Noutput),Ninstance)) source_damage_isoDuctile_output = '' - allocate(source_damage_isoDuctile_outputID(maxval(phase_Noutput),maxNinstance), source=undefined_ID) - allocate(source_damage_isoDuctile_Noutput(maxNinstance), source=0_pInt) - allocate(source_damage_isoDuctile_critPlasticStrain(maxNinstance), source=0.0_pReal) - allocate(source_damage_isoDuctile_N(maxNinstance), source=0.0_pReal) - allocate(source_damage_isoDuctile_aTol(maxNinstance), source=0.0_pReal) + allocate(source_damage_isoDuctile_outputID(maxval(phase_Noutput),Ninstance), source=undefined_ID) + allocate(source_damage_isoDuctile_Noutput(Ninstance), source=0_pInt) + + allocate(source_damage_isoDuctile_critPlasticStrain(Ninstance), source=0.0_pReal) + allocate(source_damage_isoDuctile_N(Ninstance), source=0.0_pReal) + allocate(source_damage_isoDuctile_aTol(Ninstance), source=0.0_pReal) + + do p=1, size(config_phase) + if (all(phase_source(:,p) /= SOURCE_damage_isoDuctile_ID)) cycle + enddo rewind(fileUnit) phase = 0_pInt @@ -275,26 +291,22 @@ end subroutine source_damage_isoDuctile_dotState !-------------------------------------------------------------------------------------------------- !> @brief returns local part of nonlocal damage driving force !-------------------------------------------------------------------------------------------------- -subroutine source_damage_isoDuctile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ipc, ip, el) +subroutine source_damage_isoDuctile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) use material, only: & - phaseAt, phasememberAt, & sourceState implicit none integer(pInt), intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element + phase, & + constituent real(pReal), intent(in) :: & phi real(pReal), intent(out) :: & localphiDot, & dLocalphiDot_dPhi integer(pInt) :: & - phase, constituent, sourceOffset + sourceOffset - phase = phaseAt(ipc,ip,el) - constituent = phasememberAt(ipc,ip,el) sourceOffset = source_damage_isoDuctile_offset(phase) localphiDot = 1.0_pReal - & @@ -308,25 +320,21 @@ end subroutine source_damage_isoDuctile_getRateAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief return array of local damage results !-------------------------------------------------------------------------------------------------- -function source_damage_isoDuctile_postResults(ipc,ip,el) +function source_damage_isoDuctile_postResults(phase, constituent) use material, only: & - phaseAt, phasememberAt, & sourceState implicit none - integer(pInt), intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element + integer(pInt), intent(in) :: & + phase, & + constituent real(pReal), dimension(source_damage_isoDuctile_sizePostResults( & - source_damage_isoDuctile_instance(phaseAt(ipc,ip,el)))) :: & + source_damage_isoDuctile_instance(phase))) :: & source_damage_isoDuctile_postResults integer(pInt) :: & - instance, phase, constituent, sourceOffset, o, c + instance, sourceOffset, o, c - phase = phaseAt(ipc,ip,el) - constituent = phasememberAt(ipc,ip,el) instance = source_damage_isoDuctile_instance(phase) sourceOffset = source_damage_isoDuctile_offset(phase)