diff --git a/PRIVATE b/PRIVATE index 9573ce7bd..d0d5b5a22 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 9573ce7bd2c1a7188c1aac5b83aa76d480c2bdb0 +Subproject commit d0d5b5a22be9778187b100214c782747793bb956 diff --git a/VERSION b/VERSION index 864f27f49..0ae0c6c28 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.3-1862-g0b340a6d +v2.0.3-1920-ge8c16e7f diff --git a/examples/ConfigFiles/Phase_Isotropic_AluminumIsotropic.config b/examples/ConfigFiles/Phase_Isotropic_AluminumIsotropic.config index 2a5c53ba7..b2199da85 100644 --- a/examples/ConfigFiles/Phase_Isotropic_AluminumIsotropic.config +++ b/examples/ConfigFiles/Phase_Isotropic_AluminumIsotropic.config @@ -7,7 +7,7 @@ plasticity isotropic (output) flowstress (output) strainrate -lattice_structure isotropic +lattice_structure iso c11 110.9e9 c12 58.34e9 diff --git a/examples/ConfigFiles/Phase_Isotropic_FreeSurface.config b/examples/ConfigFiles/Phase_Isotropic_FreeSurface.config index 5118c3974..ee009da3a 100644 --- a/examples/ConfigFiles/Phase_Isotropic_FreeSurface.config +++ b/examples/ConfigFiles/Phase_Isotropic_FreeSurface.config @@ -12,7 +12,7 @@ plasticity isotropic (output) flowstress (output) strainrate -lattice_structure isotropic +lattice_structure iso c11 10e9 c12 0.0 gdot0 0.001 @@ -22,4 +22,4 @@ h0 1e6 n 5 m 3 a 2 -atol_resistance 1 \ No newline at end of file +atol_resistance 1 diff --git a/examples/ConfigFiles/Phase_None_IsotropicVolumePreservation.config b/examples/ConfigFiles/Phase_None_IsotropicVolumePreservation.config index 6da2ad826..46c99985f 100644 --- a/examples/ConfigFiles/Phase_None_IsotropicVolumePreservation.config +++ b/examples/ConfigFiles/Phase_None_IsotropicVolumePreservation.config @@ -1,9 +1,8 @@ - [IsotropicVolumePreservation] elasticity hooke plasticity none ### Material parameters ### -lattice_structure isotropic +lattice_structure iso C11 100.0e9 C12 66.6666667e9 diff --git a/examples/SpectralMethod/EshelbyInclusion/material.config b/examples/SpectralMethod/EshelbyInclusion/material.config index 008c44f4b..cd926ba67 100644 --- a/examples/SpectralMethod/EshelbyInclusion/material.config +++ b/examples/SpectralMethod/EshelbyInclusion/material.config @@ -16,7 +16,7 @@ t0 330.0 #................. [isotropic matrix] -lattice_structure isotropic +lattice_structure iso plasticity none {config/elastic_isotropic.config} {config/thermal.config} @@ -45,7 +45,7 @@ plasticity none #................. [isotropic inclusion] -lattice_structure isotropic +lattice_structure iso plasticity none {config/elastic_isotropic.config} {config/thermal.config} diff --git a/src/CPFEM.f90 b/src/CPFEM.f90 index 172b79122..11b0a557f 100644 --- a/src/CPFEM.f90 +++ b/src/CPFEM.f90 @@ -27,9 +27,6 @@ module CPFEM implicit none private - real(pReal), parameter, private :: & - CPFEM_odd_stress = 1e15_pReal, & !< return value for stress in case of ping pong dummy cycle - CPFEM_odd_jacobian = 1e50_pReal !< return value for jacobian in case of ping pong dummy cycle real(pReal), dimension (:,:,:), allocatable, private :: & CPFEM_cs !< Cauchy stress real(pReal), dimension (:,:,:,:), allocatable, private :: & @@ -150,7 +147,10 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt integer(pInt) elCP, & ! crystal plasticity element number i, j, k, l, m, n, ph, homog, mySource - logical updateJaco ! flag indicating if JAcobian has to be updated + logical updateJaco ! flag indicating if Jacobian has to be updated + + real(pReal), parameter :: ODD_STRESS = 1e15_pReal, & !< return value for stress in case of ping pong dummy cycle + ODD_JACOBIAN = 1e50_pReal !< return value for jacobian in case of ping pong dummy cycle elCP = mesh_FEM2DAMASK_elem(elFE) @@ -193,8 +193,8 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt elseif (iand(mode, CPFEM_COLLECT) /= 0_pInt) then call random_number(rnd) if (rnd < 0.5_pReal) rnd = rnd - 1.0_pReal - CPFEM_cs(1:6,ip,elCP) = rnd * CPFEM_odd_stress - CPFEM_dcsde(1:6,1:6,ip,elCP) = CPFEM_odd_jacobian * math_identity2nd(6) + CPFEM_cs(1:6,ip,elCP) = rnd * ODD_STRESS + CPFEM_dcsde(1:6,1:6,ip,elCP) = ODD_JACOBIAN * math_identity2nd(6) chosenThermal2: select case (thermal_type(material_homogenizationAt(elCP))) case (THERMAL_conduction_ID) chosenThermal2 temperature(material_homogenizationAt(elCP))%p(thermalMapping(material_homogenizationAt(elCP))%p(ip,elCP)) = & @@ -226,8 +226,8 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt endif call random_number(rnd) if (rnd < 0.5_pReal) rnd = rnd - 1.0_pReal - CPFEM_cs(1:6,ip,elCP) = rnd*CPFEM_odd_stress - CPFEM_dcsde(1:6,1:6,ip,elCP) = CPFEM_odd_jacobian*math_identity2nd(6) + CPFEM_cs(1:6,ip,elCP) = ODD_STRESS * rnd + CPFEM_dcsde(1:6,1:6,ip,elCP) = ODD_JACOBIAN * math_identity2nd(6) !*** deformation gradient is not outdated @@ -257,8 +257,8 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt call random_number(rnd) if (rnd < 0.5_pReal) rnd = rnd - 1.0_pReal - CPFEM_cs(1:6,ip,elCP) = rnd * CPFEM_odd_stress - CPFEM_dcsde(1:6,1:6,ip,elCP) = CPFEM_odd_jacobian * math_identity2nd(6) + CPFEM_cs(1:6,ip,elCP) = ODD_STRESS * rnd + CPFEM_dcsde(1:6,1:6,ip,elCP) = ODD_JACOBIAN * math_identity2nd(6) else terminalIllness @@ -331,6 +331,7 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt end subroutine CPFEM_general + !-------------------------------------------------------------------------------------------------- !> @brief Forward data for new time increment. !-------------------------------------------------------------------------------------------------- diff --git a/src/DAMASK_marc.f90 b/src/DAMASK_marc.f90 index 782657da6..a3718c59d 100644 --- a/src/DAMASK_marc.f90 +++ b/src/DAMASK_marc.f90 @@ -40,7 +40,7 @@ module DAMASK_interface implicit none private - logical, public :: symmetricSolver + logical, protected, public :: symmetricSolver character(len=*), parameter, public :: INPUTFILEEXTENSION = '.dat' public :: & diff --git a/src/HDF5_utilities.f90 b/src/HDF5_utilities.f90 index b5eb2c7df..aa75d9c87 100644 --- a/src/HDF5_utilities.f90 +++ b/src/HDF5_utilities.f90 @@ -77,20 +77,6 @@ module HDF5_utilities module procedure HDF5_addAttribute_real_array end interface HDF5_addAttribute - -!-------------------------------------------------------------------------------------------------- - public :: & - HDF5_utilities_init, & - HDF5_openFile, & - HDF5_closeFile, & - HDF5_addAttribute, & - HDF5_closeGroup ,& - HDF5_openGroup, & - HDF5_addGroup, & - HDF5_read, & - HDF5_write, & - HDF5_setLink, & - HDF5_objectExists contains diff --git a/src/constitutive_plastic_dislotwin.f90 b/src/constitutive_plastic_dislotwin.f90 index 152c630ae..116abb1cb 100644 --- a/src/constitutive_plastic_dislotwin.f90 +++ b/src/constitutive_plastic_dislotwin.f90 @@ -198,10 +198,10 @@ module subroutine plastic_dislotwin_init prm%n0_sl = lattice_slip_normal(prm%N_sl,config%getString('lattice_structure'),& config%getFloat('c/a',defaultVal=0.0_pReal)) - prm%fccTwinTransNucleation = merge(.true., .false., lattice_structure(p) == LATTICE_FCC_ID) & + prm%fccTwinTransNucleation = merge(.true., .false., lattice_structure(p) == lattice_FCC_ID) & .and. (prm%N_sl(1) == 12) if(prm%fccTwinTransNucleation) & - prm%fcc_twinNucleationSlipPair = lattice_fcc_twinNucleationSlipPair + prm%fcc_twinNucleationSlipPair = lattice_FCC_TWINNUCLEATIONSLIPPAIR prm%rho_mob_0 = config%getFloats('rhoedge0', requiredSize=size(prm%N_sl)) prm%rho_dip_0 = config%getFloats('rhoedgedip0',requiredSize=size(prm%N_sl)) @@ -230,7 +230,7 @@ module subroutine plastic_dislotwin_init prm%omega = config%getFloat('omega', defaultVal = 1000.0_pReal) & * merge(12.0_pReal, & 8.0_pReal, & - lattice_structure(p) == LATTICE_FCC_ID .or. lattice_structure(p) == LATTICE_HEX_ID) + lattice_structure(p) == lattice_FCC_ID .or. lattice_structure(p) == lattice_HEX_ID) ! expand: family => system @@ -335,7 +335,7 @@ module subroutine plastic_dislotwin_init config%getFloat('a_bcc', defaultVal=0.0_pReal), & config%getFloat('a_fcc', defaultVal=0.0_pReal)) - if (lattice_structure(p) /= LATTICE_fcc_ID) then + if (lattice_structure(p) /= lattice_FCC_ID) then prm%dot_N_0_tr = config%getFloats('ndot0_trans') prm%dot_N_0_tr = math_expand(prm%dot_N_0_tr,prm%N_tr) endif diff --git a/src/damage_local.f90 b/src/damage_local.f90 index 6296631e5..1c2f51056 100644 --- a/src/damage_local.f90 +++ b/src/damage_local.f90 @@ -16,24 +16,18 @@ module damage_local implicit none private - enum, bind(c) - enumerator :: & - undefined_ID, & - damage_ID - end enum - type :: tParameters - integer(kind(undefined_ID)), dimension(:), allocatable :: & - outputID + character(len=pStringLen), allocatable, dimension(:) :: & + output end type tParameters - + type(tparameters), dimension(:), allocatable :: & param - + public :: & damage_local_init, & damage_local_updateState, & - damage_local_Results + damage_local_results contains @@ -43,41 +37,30 @@ contains !-------------------------------------------------------------------------------------------------- subroutine damage_local_init - integer :: maxNinstance,o,NofMyHomog,h - character(len=pStringLen), dimension(:), allocatable :: outputs - - write(6,'(/,a)') ' <<<+- damage_'//DAMAGE_local_label//' init -+>>>'; flush(6) + integer :: Ninstance,NofMyHomog,h - maxNinstance = count(damage_type == DAMAGE_local_ID) - if (maxNinstance == 0) return - - allocate(param(maxNinstance)) - - do h = 1, size(damage_type) + write(6,'(/,a)') ' <<<+- damage_'//DAMAGE_local_label//' init -+>>>'; flush(6) + + Ninstance = count(damage_type == DAMAGE_local_ID) + allocate(param(Ninstance)) + + do h = 1, size(config_homogenization) if (damage_type(h) /= DAMAGE_LOCAL_ID) cycle associate(prm => param(damage_typeInstance(h)),config => config_homogenization(h)) - - outputs = config%getStrings('(output)',defaultVal=emptyStringArray) - allocate(prm%outputID(0)) - - do o=1, size(outputs) - select case(outputs(o)) - case ('damage') - prm%outputID = [prm%outputID , damage_ID] - end select - enddo + + prm%output = config%getStrings('(output)',defaultVal=emptyStringArray) NofMyHomog = count(material_homogenizationAt == h) damageState(h)%sizeState = 1 allocate(damageState(h)%state0 (1,NofMyHomog), source=damage_initialPhi(h)) allocate(damageState(h)%subState0(1,NofMyHomog), source=damage_initialPhi(h)) allocate(damageState(h)%state (1,NofMyHomog), source=damage_initialPhi(h)) - + nullify(damageMapping(h)%p) damageMapping(h)%p => material_homogenizationMemberAt deallocate(damage(h)%p) damage(h)%p => damageState(h)%state(1,:) - + end associate enddo @@ -85,10 +68,10 @@ end subroutine damage_local_init !-------------------------------------------------------------------------------------------------- -!> @brief calculates local change in damage field +!> @brief calculates local change in damage field !-------------------------------------------------------------------------------------------------- function damage_local_updateState(subdt, ip, el) - + integer, intent(in) :: & ip, & !< integration point number el !< element number @@ -100,30 +83,30 @@ function damage_local_updateState(subdt, ip, el) homog, & offset real(pReal) :: & - phi, phiDot, dPhiDot_dPhi - + phi, phiDot, dPhiDot_dPhi + homog = material_homogenizationAt(el) offset = material_homogenizationMemberAt(ip,el) phi = damageState(homog)%subState0(1,offset) call damage_local_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el) phi = max(residualStiffness,min(1.0_pReal,phi + subdt*phiDot)) - + damage_local_updateState = [ abs(phi - damageState(homog)%state(1,offset)) & <= err_damage_tolAbs & .or. abs(phi - damageState(homog)%state(1,offset)) & <= err_damage_tolRel*abs(damageState(homog)%state(1,offset)), & .true.] - damageState(homog)%state(1,offset) = phi + damageState(homog)%state(1,offset) = phi end function damage_local_updateState !-------------------------------------------------------------------------------------------------- -!> @brief calculates homogenized local damage driving forces +!> @brief calculates homogenized local damage driving forces !-------------------------------------------------------------------------------------------------- subroutine damage_local_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el) - + integer, intent(in) :: & ip, & !< integration point number el !< element number @@ -135,7 +118,7 @@ subroutine damage_local_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el source, & constituent real(pReal) :: & - phiDot, dPhiDot_dPhi, localphiDot, dLocalphiDot_dPhi + phiDot, dPhiDot_dPhi, localphiDot, dLocalphiDot_dPhi phiDot = 0.0_pReal dPhiDot_dPhi = 0.0_pReal @@ -143,7 +126,7 @@ subroutine damage_local_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el phase = material_phaseAt(grain,el) constituent = material_phasememberAt(grain,ip,el) do source = 1, phase_Nsources(phase) - select case(phase_source(source,phase)) + select case(phase_source(source,phase)) case (SOURCE_damage_isoBrittle_ID) call source_damage_isobrittle_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) @@ -163,12 +146,12 @@ subroutine damage_local_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el end select phiDot = phiDot + localphiDot dPhiDot_dPhi = dPhiDot_dPhi + dLocalphiDot_dPhi - enddo + enddo enddo - + phiDot = phiDot/real(homogenization_Ngrains(material_homogenizationAt(el)),pReal) dPhiDot_dPhi = dPhiDot_dPhi/real(homogenization_Ngrains(material_homogenizationAt(el)),pReal) - + end subroutine damage_local_getSourceAndItsTangent @@ -179,14 +162,13 @@ subroutine damage_local_results(homog,group) integer, intent(in) :: homog character(len=*), intent(in) :: group - integer :: o - - associate(prm => param(damage_typeInstance(homog))) - outputsLoop: do o = 1,size(prm%outputID) - select case(prm%outputID(o)) - - case (damage_ID) + integer :: o + + associate(prm => param(damage_typeInstance(homog))) + outputsLoop: do o = 1,size(prm%output) + select case(prm%output(o)) + case ('damage') call results_writeDataset(group,damage(homog)%p,'phi',& 'damage indicator','-') end select diff --git a/src/damage_none.f90 b/src/damage_none.f90 index 239f98c4e..b8a9bd7b5 100644 --- a/src/damage_none.f90 +++ b/src/damage_none.f90 @@ -8,7 +8,7 @@ module damage_none implicit none public - + contains !-------------------------------------------------------------------------------------------------- @@ -28,10 +28,10 @@ subroutine damage_none_init allocate(damageState(h)%state0 (0,NofMyHomog)) allocate(damageState(h)%subState0(0,NofMyHomog)) allocate(damageState(h)%state (0,NofMyHomog)) - + deallocate(damage(h)%p) allocate (damage(h)%p(1), source=damage_initialPhi(h)) - + enddo end subroutine damage_none_init diff --git a/src/damage_nonlocal.f90 b/src/damage_nonlocal.f90 index 4587857aa..9644211c4 100644 --- a/src/damage_nonlocal.f90 +++ b/src/damage_nonlocal.f90 @@ -18,27 +18,21 @@ module damage_nonlocal implicit none private - enum, bind(c) - enumerator :: & - undefined_ID, & - damage_ID - end enum - type :: tParameters - integer(kind(undefined_ID)), dimension(:), allocatable :: & - outputID + character(len=pStringLen), allocatable, dimension(:) :: & + output end type tParameters - + type(tparameters), dimension(:), allocatable :: & param public :: & damage_nonlocal_init, & damage_nonlocal_getSourceAndItsTangent, & - damage_nonlocal_getDiffusion33, & + damage_nonlocal_getDiffusion, & damage_nonlocal_getMobility, & damage_nonlocal_putNonLocalDamage, & - damage_nonlocal_Results + damage_nonlocal_results contains @@ -48,29 +42,18 @@ contains !-------------------------------------------------------------------------------------------------- subroutine damage_nonlocal_init - integer :: maxNinstance,o,NofMyHomog,h - character(len=pStringLen), dimension(:), allocatable :: outputs + integer :: Ninstance,NofMyHomog,h - write(6,'(/,a)') ' <<<+- damage_'//DAMAGE_nonlocal_label//' init -+>>>'; flush(6) - - maxNinstance = count(damage_type == DAMAGE_nonlocal_ID) - if (maxNinstance == 0) return + write(6,'(/,a)') ' <<<+- damage_'//DAMAGE_nonlocal_label//' init -+>>>'; flush(6) - allocate(param(maxNinstance)) - - do h = 1, size(damage_type) + Ninstance = count(damage_type == DAMAGE_nonlocal_ID) + allocate(param(Ninstance)) + + do h = 1, size(config_homogenization) if (damage_type(h) /= DAMAGE_NONLOCAL_ID) cycle associate(prm => param(damage_typeInstance(h)),config => config_homogenization(h)) - - outputs = config%getStrings('(output)',defaultVal=emptyStringArray) - allocate(prm%outputID(0)) - - do o=1, size(outputs) - select case(outputs(o)) - case ('damage') - prm%outputID = [prm%outputID, damage_ID] - end select - enddo + + prm%output = config%getStrings('(output)',defaultVal=emptyStringArray) NofMyHomog = count(material_homogenizationAt == h) damageState(h)%sizeState = 1 @@ -82,7 +65,7 @@ subroutine damage_nonlocal_init damageMapping(h)%p => material_homogenizationMemberAt deallocate(damage(h)%p) damage(h)%p => damageState(h)%state(1,:) - + end associate enddo @@ -90,10 +73,10 @@ end subroutine damage_nonlocal_init !-------------------------------------------------------------------------------------------------- -!> @brief calculates homogenized damage driving forces +!> @brief calculates homogenized damage driving forces !-------------------------------------------------------------------------------------------------- subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el) - + integer, intent(in) :: & ip, & !< integration point number el !< element number @@ -105,7 +88,7 @@ subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, source, & constituent real(pReal) :: & - phiDot, dPhiDot_dPhi, localphiDot, dLocalphiDot_dPhi + phiDot, dPhiDot_dPhi, localphiDot, dLocalphiDot_dPhi phiDot = 0.0_pReal dPhiDot_dPhi = 0.0_pReal @@ -113,7 +96,7 @@ subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, phase = material_phaseAt(grain,el) constituent = material_phasememberAt(grain,ip,el) do source = 1, phase_Nsources(phase) - select case(phase_source(source,phase)) + select case(phase_source(source,phase)) case (SOURCE_damage_isoBrittle_ID) call source_damage_isobrittle_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) @@ -133,44 +116,44 @@ subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, end select phiDot = phiDot + localphiDot dPhiDot_dPhi = dPhiDot_dPhi + dLocalphiDot_dPhi - enddo + enddo enddo - + phiDot = phiDot/real(homogenization_Ngrains(material_homogenizationAt(el)),pReal) dPhiDot_dPhi = dPhiDot_dPhi/real(homogenization_Ngrains(material_homogenizationAt(el)),pReal) - + end subroutine damage_nonlocal_getSourceAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief returns homogenized non local damage diffusion tensor in reference configuration !-------------------------------------------------------------------------------------------------- -function damage_nonlocal_getDiffusion33(ip,el) +function damage_nonlocal_getDiffusion(ip,el) integer, intent(in) :: & ip, & !< integration point number el !< element number real(pReal), dimension(3,3) :: & - damage_nonlocal_getDiffusion33 + damage_nonlocal_getDiffusion integer :: & homog, & grain - + homog = material_homogenizationAt(el) - damage_nonlocal_getDiffusion33 = 0.0_pReal + damage_nonlocal_getDiffusion = 0.0_pReal do grain = 1, homogenization_Ngrains(homog) - damage_nonlocal_getDiffusion33 = damage_nonlocal_getDiffusion33 + & - crystallite_push33ToRef(grain,ip,el,lattice_DamageDiffusion33(1:3,1:3,material_phaseAt(grain,el))) + damage_nonlocal_getDiffusion = damage_nonlocal_getDiffusion + & + crystallite_push33ToRef(grain,ip,el,lattice_DamageDiffusion(1:3,1:3,material_phaseAt(grain,el))) enddo - damage_nonlocal_getDiffusion33 = & - charLength**2*damage_nonlocal_getDiffusion33/real(homogenization_Ngrains(homog),pReal) - -end function damage_nonlocal_getDiffusion33 - - + damage_nonlocal_getDiffusion = & + charLength**2*damage_nonlocal_getDiffusion/real(homogenization_Ngrains(homog),pReal) + +end function damage_nonlocal_getDiffusion + + !-------------------------------------------------------------------------------------------------- -!> @brief Returns homogenized nonlocal damage mobility +!> @brief Returns homogenized nonlocal damage mobility !-------------------------------------------------------------------------------------------------- real(pReal) function damage_nonlocal_getMobility(ip,el) @@ -179,9 +162,9 @@ real(pReal) function damage_nonlocal_getMobility(ip,el) el !< element number integer :: & ipc - + damage_nonlocal_getMobility = 0.0_pReal - + do ipc = 1, homogenization_Ngrains(material_homogenizationAt(el)) damage_nonlocal_getMobility = damage_nonlocal_getMobility + lattice_DamageMobility(material_phaseAt(ipc,el)) enddo @@ -205,7 +188,7 @@ subroutine damage_nonlocal_putNonLocalDamage(phi,ip,el) integer :: & homog, & offset - + homog = material_homogenizationAt(el) offset = damageMapping(homog)%p(ip,el) damage(homog)%p(offset) = phi @@ -220,14 +203,13 @@ subroutine damage_nonlocal_results(homog,group) integer, intent(in) :: homog character(len=*), intent(in) :: group - integer :: o - - associate(prm => param(damage_typeInstance(homog))) - outputsLoop: do o = 1,size(prm%outputID) - select case(prm%outputID(o)) - - case (damage_ID) + integer :: o + + associate(prm => param(damage_typeInstance(homog))) + outputsLoop: do o = 1,size(prm%output) + select case(prm%output(o)) + case ('damage') call results_writeDataset(group,damage(homog)%p,'phi',& 'damage indicator','-') end select diff --git a/src/grid/grid_damage_spectral.f90 b/src/grid/grid_damage_spectral.f90 index 5b55f3132..8460ea8be 100644 --- a/src/grid/grid_damage_spectral.f90 +++ b/src/grid/grid_damage_spectral.f90 @@ -249,8 +249,8 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr) cell = 0 do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) cell = cell + 1 - vectorField_real(1:3,i,j,k) = matmul(damage_nonlocal_getDiffusion33(1,cell) - K_ref, & - vectorField_real(1:3,i,j,k)) + vectorField_real(1:3,i,j,k) = matmul(damage_nonlocal_getDiffusion(1,cell) - K_ref, & + vectorField_real(1:3,i,j,k)) enddo; enddo; enddo call utilities_FFTvectorForward call utilities_fourierVectorDivergence !< calculate damage divergence in fourier field @@ -294,7 +294,7 @@ subroutine updateReference mu_ref = 0.0_pReal do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) cell = cell + 1 - K_ref = K_ref + damage_nonlocal_getDiffusion33(1,cell) + K_ref = K_ref + damage_nonlocal_getDiffusion(1,cell) mu_ref = mu_ref + damage_nonlocal_getMobility(1,cell) enddo; enddo; enddo K_ref = K_ref*wgt diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90 index 5dfc79e3f..1fe6aeab0 100644 --- a/src/grid/grid_thermal_spectral.f90 +++ b/src/grid/grid_thermal_spectral.f90 @@ -252,8 +252,8 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr) cell = 0 do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) cell = cell + 1 - vectorField_real(1:3,i,j,k) = matmul(thermal_conduction_getConductivity33(1,cell) - K_ref, & - vectorField_real(1:3,i,j,k)) + vectorField_real(1:3,i,j,k) = matmul(thermal_conduction_getConductivity(1,cell) - K_ref, & + vectorField_real(1:3,i,j,k)) enddo; enddo; enddo call utilities_FFTvectorForward call utilities_fourierVectorDivergence !< calculate temperature divergence in fourier field @@ -294,9 +294,8 @@ subroutine updateReference mu_ref = 0.0_pReal do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) cell = cell + 1 - K_ref = K_ref + thermal_conduction_getConductivity33(1,cell) - mu_ref = mu_ref + thermal_conduction_getMassDensity(1,cell)* & - thermal_conduction_getSpecificHeat(1,cell) + K_ref = K_ref + thermal_conduction_getConductivity(1,cell) + mu_ref = mu_ref + thermal_conduction_getMassDensity(1,cell)* thermal_conduction_getSpecificHeat(1,cell) enddo; enddo; enddo K_ref = K_ref*wgt call MPI_Allreduce(MPI_IN_PLACE,K_ref,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) diff --git a/src/kinematics_cleavage_opening.f90 b/src/kinematics_cleavage_opening.f90 index 6b060a8d9..a5f5dfae2 100644 --- a/src/kinematics_cleavage_opening.f90 +++ b/src/kinematics_cleavage_opening.f90 @@ -27,25 +27,12 @@ module kinematics_cleavage_opening sdot0, & n real(pReal), dimension(:), allocatable :: & - critDisp, & critLoad - end type + real(pReal), dimension(:,:,:,:), allocatable :: & + cleavage_systems + end type tParameters -! Begin Deprecated - integer, dimension(:), allocatable :: & - kinematics_cleavage_opening_totalNcleavage !< total number of cleavage systems - - integer, dimension(:,:), allocatable :: & - kinematics_cleavage_opening_Ncleavage !< number of cleavage systems per family - - real(pReal), dimension(:), allocatable :: & - kinematics_cleavage_opening_sdot_0, & - kinematics_cleavage_opening_N - - real(pReal), dimension(:,:), allocatable :: & - kinematics_cleavage_opening_critDisp, & - kinematics_cleavage_opening_critLoad -! End Deprecated + type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance) public :: & kinematics_cleavage_opening_init, & @@ -60,66 +47,59 @@ contains !-------------------------------------------------------------------------------------------------- subroutine kinematics_cleavage_opening_init - integer, allocatable, dimension(:) :: tempInt - real(pReal), allocatable, dimension(:) :: tempFloat - - integer :: maxNinstance,p,instance + integer :: Ninstance,p + character(len=pStringLen) :: extmsg = '' write(6,'(/,a)') ' <<<+- kinematics_'//KINEMATICS_cleavage_opening_LABEL//' init -+>>>'; flush(6) - maxNinstance = count(phase_kinematics == KINEMATICS_cleavage_opening_ID) - if (maxNinstance == 0) return - + Ninstance = count(phase_kinematics == KINEMATICS_cleavage_opening_ID) if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) & - write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance - + write(6,'(a16,1x,i5,/)') '# instances:',Ninstance + allocate(kinematics_cleavage_opening_instance(size(config_phase)), source=0) - do p = 1, 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_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) - allocate(kinematics_cleavage_opening_totalNcleavage(maxNinstance), source=0) - allocate(kinematics_cleavage_opening_sdot_0(maxNinstance), source=0.0_pReal) - allocate(kinematics_cleavage_opening_N(maxNinstance), source=0.0_pReal) + allocate(param(Ninstance)) do p = 1, size(config_phase) + kinematics_cleavage_opening_instance(p) = count(phase_kinematics(:,1:p) == kinematics_cleavage_opening_ID) 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 - tempFloat = config_phase(p)%getFloats('anisobrittle_criticaldisplacement',requiredSize=size(tempInt)) - kinematics_cleavage_opening_critDisp(1:size(tempInt),instance) = tempFloat + associate(prm => param(kinematics_cleavage_opening_instance(p)), & + config => config_phase(p)) - tempFloat = config_phase(p)%getFloats('anisobrittle_criticalload',requiredSize=size(tempInt)) - kinematics_cleavage_opening_critLoad(1:size(tempInt),instance) = tempFloat + prm%Ncleavage = config%getInts('ncleavage') + prm%totalNcleavage = sum(prm%Ncleavage) + + prm%n = config%getFloat('anisobrittle_ratesensitivity') + prm%sdot0 = config%getFloat('anisobrittle_sdot0') + + prm%critLoad = config%getFloats('anisobrittle_criticalload',requiredSize=size(prm%Ncleavage)) + + prm%cleavage_systems = lattice_SchmidMatrix_cleavage(prm%Ncleavage,config%getString('lattice_structure'),& + config%getFloat('c/a',defaultVal=0.0_pReal)) + + ! expand: family => system + prm%critLoad = math_expand(prm%critLoad, prm%Ncleavage) + + ! sanity checks + if (prm%n <= 0.0_pReal) extmsg = trim(extmsg)//' anisobrittle_n' + if (prm%sdot0 <= 0.0_pReal) extmsg = trim(extmsg)//' anisobrittle_sdot0' + if (any(prm%critLoad < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_critLoad' - 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,el=instance,ext_msg='sdot_0 ('//KINEMATICS_cleavage_opening_LABEL//')') - if (any(kinematics_cleavage_opening_critDisp(1:size(tempInt),instance) < 0.0_pReal)) & - call IO_error(211,el=instance,ext_msg='critical_displacement ('//KINEMATICS_cleavage_opening_LABEL//')') - if (any(kinematics_cleavage_opening_critLoad(1:size(tempInt),instance) < 0.0_pReal)) & - call IO_error(211,el=instance,ext_msg='critical_load ('//KINEMATICS_cleavage_opening_LABEL//')') - if (kinematics_cleavage_opening_N(instance) <= 0.0_pReal) & - call IO_error(211,el=instance,ext_msg='rate_sensitivity ('//KINEMATICS_cleavage_opening_LABEL//')') - enddo - -end subroutine kinematics_cleavage_opening_init - !-------------------------------------------------------------------------------------------------- -!> @brief contains the constitutive equation for calculating the velocity gradient +! exit if any parameter is out of range + if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'('//SOURCE_DAMAGE_ANISOBRITTLE_LABEL//')') + + end associate + enddo + +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) - + integer, intent(in) :: & ipc, & !< grain number ip, & !< integration point number @@ -130,73 +110,55 @@ subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, i 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) + integer :: & - instance, phase, & homog, damageOffset, & - f, i, index_myFamily, k, l, m, n + i, k, l, m, n real(pReal) :: & traction_d, traction_t, traction_n, traction_crit, & udotd, dudotd_dt, udott, dudott_dt, udotn, dudotn_dt - phase = material_phaseAt(ipc,el) - instance = kinematics_cleavage_opening_instance(phase) homog = material_homogenizationAt(el) damageOffset = damageMapping(homog)%p(ip,el) - + Ld = 0.0_pReal dLd_dTstar = 0.0_pReal - do f = 1,lattice_maxNcleavageFamily - index_myFamily = sum(lattice_NcleavageSystem(1:f-1,phase)) ! at which index starts my family - do i = 1,kinematics_cleavage_opening_Ncleavage(f,instance) ! process each (active) cleavage system in family - traction_d = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,1,index_myFamily+i,phase)) - traction_t = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,2,index_myFamily+i,phase)) - traction_n = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,3,index_myFamily+i,phase)) - traction_crit = kinematics_cleavage_opening_critLoad(f,instance)* & - damage(homog)%p(damageOffset)*damage(homog)%p(damageOffset) - udotd = & - sign(1.0_pReal,traction_d)* & - kinematics_cleavage_opening_sdot_0(instance)* & - (max(0.0_pReal, abs(traction_d) - traction_crit)/traction_crit)**kinematics_cleavage_opening_N(instance) - if (abs(udotd) > tol_math_check) then - Ld = Ld + udotd*lattice_Scleavage(1:3,1:3,1,index_myFamily+i,phase) - 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:3,l=1:3,m=1:3,n=1:3) & - 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 + associate(prm => param(kinematics_cleavage_opening_instance(material_phaseAt(ipc,el)))) + do i = 1,prm%totalNcleavage + traction_crit = prm%critLoad(i)* damage(homog)%p(damageOffset)**2.0_pReal - udott = & - sign(1.0_pReal,traction_t)* & - kinematics_cleavage_opening_sdot_0(instance)* & - (max(0.0_pReal, abs(traction_t) - traction_crit)/traction_crit)**kinematics_cleavage_opening_N(instance) - if (abs(udott) > tol_math_check) then - Ld = Ld + udott*lattice_Scleavage(1:3,1:3,2,index_myFamily+i,phase) - 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:3,l=1:3,m=1:3,n=1:3) & - 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 + traction_d = math_mul33xx33(S,prm%cleavage_systems(1:3,1:3,1,i)) + if (abs(traction_d) > traction_crit + tol_math_check) then + udotd = sign(1.0_pReal,traction_d)* prm%sdot0 * ((abs(traction_d) - traction_crit)/traction_crit)**prm%n + Ld = Ld + udotd*prm%cleavage_systems(1:3,1:3,1,i) + dudotd_dt = sign(1.0_pReal,traction_d)*udotd*prm%n / (abs(traction_d) - traction_crit) + forall (k=1:3,l=1:3,m=1:3,n=1:3) & + dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) & + + dudotd_dt*prm%cleavage_systems(k,l,1,i) * prm%cleavage_systems(m,n,1,i) + endif - udotn = & - sign(1.0_pReal,traction_n)* & - kinematics_cleavage_opening_sdot_0(instance)* & - (max(0.0_pReal, abs(traction_n) - traction_crit)/traction_crit)**kinematics_cleavage_opening_N(instance) - if (abs(udotn) > tol_math_check) then - Ld = Ld + udotn*lattice_Scleavage(1:3,1:3,3,index_myFamily+i,phase) - 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:3,l=1:3,m=1:3,n=1:3) & - 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 + traction_t = math_mul33xx33(S,prm%cleavage_systems(1:3,1:3,2,i)) + if (abs(traction_t) > traction_crit + tol_math_check) then + udott = sign(1.0_pReal,traction_t)* prm%sdot0 * ((abs(traction_t) - traction_crit)/traction_crit)**prm%n + Ld = Ld + udott*prm%cleavage_systems(1:3,1:3,2,i) + dudott_dt = sign(1.0_pReal,traction_t)*udott*prm%n / (abs(traction_t) - traction_crit) + forall (k=1:3,l=1:3,m=1:3,n=1:3) & + dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) & + + dudott_dt*prm%cleavage_systems(k,l,2,i) * prm%cleavage_systems(m,n,2,i) + endif + + traction_n = math_mul33xx33(S,prm%cleavage_systems(1:3,1:3,3,i)) + if (abs(traction_n) > traction_crit + tol_math_check) then + udotn = sign(1.0_pReal,traction_n)* prm%sdot0 * ((abs(traction_n) - traction_crit)/traction_crit)**prm%n + Ld = Ld + udotn*prm%cleavage_systems(1:3,1:3,3,i) + dudotn_dt = sign(1.0_pReal,traction_n)*udotn*prm%n / (abs(traction_n) - traction_crit) + forall (k=1:3,l=1:3,m=1:3,n=1:3) & + dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) & + + dudotn_dt*prm%cleavage_systems(k,l,3,i) * prm%cleavage_systems(m,n,3,i) + endif enddo - + end associate + end subroutine kinematics_cleavage_opening_LiAndItsTangent end module kinematics_cleavage_opening diff --git a/src/kinematics_slipplane_opening.f90 b/src/kinematics_slipplane_opening.f90 index 2c94448bd..ab9035918 100644 --- a/src/kinematics_slipplane_opening.f90 +++ b/src/kinematics_slipplane_opening.f90 @@ -12,12 +12,12 @@ module kinematics_slipplane_opening use math use lattice use material - + implicit none private - + integer, dimension(:), allocatable :: kinematics_slipplane_opening_instance - + type :: tParameters !< container type for internal constitutive parameters integer :: & totalNslip @@ -28,14 +28,14 @@ module kinematics_slipplane_opening n real(pReal), dimension(:), allocatable :: & critLoad - real(pReal), dimension(:,:), allocatable :: & - slip_direction, & - slip_normal, & - slip_transverse + real(pReal), dimension(:,:,:), allocatable :: & + P_d, & + P_t, & + P_n end type tParameters - + type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance) - + public :: & kinematics_slipplane_opening_init, & kinematics_slipplane_opening_LiAndItsTangent @@ -49,58 +49,66 @@ contains !-------------------------------------------------------------------------------------------------- subroutine kinematics_slipplane_opening_init - integer :: maxNinstance,p,instance + integer :: Ninstance,p,i + character(len=pStringLen) :: extmsg = '' + real(pReal), dimension(:,:), allocatable :: d,n,t write(6,'(/,a)') ' <<<+- kinematics_'//KINEMATICS_slipplane_opening_LABEL//' init -+>>>'; flush(6) - maxNinstance = count(phase_kinematics == KINEMATICS_slipplane_opening_ID) - if (maxNinstance == 0) return - + Ninstance = count(phase_kinematics == KINEMATICS_slipplane_opening_ID) if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) & - write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance - + write(6,'(a16,1x,i5,/)') '# instances:',Ninstance + allocate(kinematics_slipplane_opening_instance(size(config_phase)), source=0) + allocate(param(Ninstance)) + do p = 1, size(config_phase) - kinematics_slipplane_opening_instance(p) = count(phase_kinematics(:,1:p) == kinematics_slipplane_opening_ID) ! ToDo: count correct? - enddo - - allocate(param(maxNinstance)) - - do p = 1, size(config_phase) + kinematics_slipplane_opening_instance(p) = count(phase_kinematics(:,1:p) == kinematics_slipplane_opening_ID) if (all(phase_kinematics(:,p) /= KINEMATICS_slipplane_opening_ID)) cycle associate(prm => param(kinematics_slipplane_opening_instance(p)), & config => config_phase(p)) - instance = kinematics_slipplane_opening_instance(p) - prm%sdot0 = config_phase(p)%getFloat('anisoductile_sdot0') - prm%n = config_phase(p)%getFloat('anisoductile_ratesensitivity') - - prm%Nslip = config%getInts('nslip') - prm%critLoad = config_phase(p)%getFloats('anisoductile_criticalload',requiredSize=size(prm%Nslip )) + prm%sdot0 = config%getFloat('anisoductile_sdot0') + prm%n = config%getFloat('anisoductile_ratesensitivity') + prm%Nslip = config%getInts('nslip') + prm%totalNslip = sum(prm%Nslip) + d = lattice_slip_direction (prm%Nslip,config%getString('lattice_structure'),& + config%getFloat('c/a',defaultVal=0.0_pReal)) + t = lattice_slip_transverse(prm%Nslip,config%getString('lattice_structure'),& + config%getFloat('c/a',defaultVal=0.0_pReal)) + n = lattice_slip_normal (prm%Nslip,config%getString('lattice_structure'),& + config%getFloat('c/a',defaultVal=0.0_pReal)) + allocate(prm%P_d(3,3,size(d,2)),prm%P_t(3,3,size(t,2)),prm%P_n(3,3,size(n,2))) + + do i=1, size(n,2) + prm%P_d(1:3,1:3,i) = math_outer(d(1:3,i), n(1:3,i)) + prm%P_t(1:3,1:3,i) = math_outer(t(1:3,i), n(1:3,i)) + prm%P_n(1:3,1:3,i) = math_outer(n(1:3,i), n(1:3,i)) + enddo + + prm%critLoad = config%getFloats('anisoductile_criticalload',requiredSize=size(prm%Nslip)) + + ! expand: family => system prm%critLoad = math_expand(prm%critLoad, prm%Nslip) - - prm%slip_direction = lattice_slip_direction (prm%Nslip,config%getString('lattice_structure'),& - config%getFloat('c/a',defaultVal=0.0_pReal)) - prm%slip_normal = lattice_slip_normal (prm%Nslip,config%getString('lattice_structure'),& - config%getFloat('c/a',defaultVal=0.0_pReal)) - prm%slip_transverse = lattice_slip_transverse(prm%Nslip,config%getString('lattice_structure'),& - config%getFloat('c/a',defaultVal=0.0_pReal)) - ! if (kinematics_slipplane_opening_sdot_0(instance) <= 0.0_pReal) & - ! call IO_error(211,el=instance,ext_msg='sdot_0 ('//KINEMATICS_slipplane_opening_LABEL//')') - ! if (any(kinematics_slipplane_opening_critPlasticStrain(:,instance) < 0.0_pReal)) & - ! call IO_error(211,el=instance,ext_msg='criticaPlasticStrain ('//KINEMATICS_slipplane_opening_LABEL//')') - ! if (kinematics_slipplane_opening_N(instance) <= 0.0_pReal) & - ! call IO_error(211,el=instance,ext_msg='rate_sensitivity ('//KINEMATICS_slipplane_opening_LABEL//')') - - end associate - enddo - -end subroutine kinematics_slipplane_opening_init + ! sanity checks + if (prm%n <= 0.0_pReal) extmsg = trim(extmsg)//' anisoDuctile_n' + if (prm%sdot0 <= 0.0_pReal) extmsg = trim(extmsg)//' anisoDuctile_sdot0' + if (any(prm%critLoad < 0.0_pReal)) extmsg = trim(extmsg)//' anisoDuctile_critLoad' !-------------------------------------------------------------------------------------------------- -!> @brief contains the constitutive equation for calculating the velocity gradient +! exit if any parameter is out of range + if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'('//SOURCE_DAMAGE_ANISODUCTILE_LABEL//')') + + end associate + enddo + +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) @@ -114,8 +122,7 @@ subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, Ld !< damage velocity gradient real(pReal), intent(out), dimension(3,3,3,3) :: & dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor) - real(pReal), dimension(3,3) :: & - projection_d, projection_t, projection_n !< projection modes 3x3 tensor + integer :: & instance, phase, & homog, damageOffset, & @@ -123,26 +130,22 @@ subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, real(pReal) :: & traction_d, traction_t, traction_n, traction_crit, & udotd, dudotd_dt, udott, dudott_dt, udotn, dudotn_dt - + phase = material_phaseAt(ipc,el) instance = kinematics_slipplane_opening_instance(phase) homog = material_homogenizationAt(el) damageOffset = damageMapping(homog)%p(ip,el) - + associate(prm => param(instance)) Ld = 0.0_pReal dLd_dTstar = 0.0_pReal do i = 1, prm%totalNslip - projection_d = math_outer(prm%slip_direction(1:3,i), prm%slip_normal(1:3,i)) - projection_t = math_outer(prm%slip_transverse(1:3,i),prm%slip_normal(1:3,i)) - projection_n = math_outer(prm%slip_normal(1:3,i), prm%slip_normal(1:3,i)) - - traction_d = math_mul33xx33(S,projection_d) - traction_t = math_mul33xx33(S,projection_t) - traction_n = math_mul33xx33(S,projection_n) - - traction_crit = prm%critLoad(i)* damage(homog)%p(damageOffset) ! degrading critical load carrying capacity by damage + traction_d = math_tensordot(S,prm%P_d(1:3,1:3,i)) + traction_t = math_tensordot(S,prm%P_t(1:3,1:3,i)) + traction_n = math_tensordot(S,prm%P_n(1:3,1:3,i)) + + traction_crit = prm%critLoad(i)* damage(homog)%p(damageOffset) ! degrading critical load carrying capacity by damage udotd = sign(1.0_pReal,traction_d)* prm%sdot0* ( abs(traction_d)/traction_crit & - abs(traction_d)/prm%critLoad(i))**prm%n @@ -151,20 +154,34 @@ subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, udotn = prm%sdot0* ( max(0.0_pReal,traction_n)/traction_crit & - max(0.0_pReal,traction_n)/prm%critLoad(i))**prm%n - dudotd_dt = udotd*prm%n/traction_d - dudott_dt = udott*prm%n/traction_t - dudotn_dt = udotn*prm%n/traction_n + if (dNeq0(traction_d)) then + dudotd_dt = udotd*prm%n/traction_d + else + dudotd_dt = 0.0_pReal + endif + if (dNeq0(traction_t)) then + dudott_dt = udott*prm%n/traction_t + else + dudott_dt = 0.0_pReal + endif + if (dNeq0(traction_n)) then + dudotn_dt = udotn*prm%n/traction_n + else + dudotn_dt = 0.0_pReal + endif forall (k=1:3,l=1:3,m=1:3,n=1:3) & - dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + dudotd_dt*projection_d(k,l)*projection_d(m,n) & - + dudott_dt*projection_t(k,l)*projection_t(m,n) & - + dudotn_dt*projection_n(k,l)*projection_n(m,n) + dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) & + + dudotd_dt*prm%P_d(k,l,i)*prm%P_d(m,n,i) & + + dudott_dt*prm%P_t(k,l,i)*prm%P_t(m,n,i) & + + dudotn_dt*prm%P_n(k,l,i)*prm%P_n(m,n,i) - Ld = Ld + udotd*projection_d & - + udott*projection_t & - + udotn*projection_n + Ld = Ld & + + udotd*prm%P_d(1:3,1:3,i) & + + udott*prm%P_t(1:3,1:3,i) & + + udotn*prm%P_n(1:3,1:3,i) enddo - + end associate end subroutine kinematics_slipplane_opening_LiAndItsTangent diff --git a/src/kinematics_thermal_expansion.f90 b/src/kinematics_thermal_expansion.f90 index 7f7994959..acf3a5067 100644 --- a/src/kinematics_thermal_expansion.f90 +++ b/src/kinematics_thermal_expansion.f90 @@ -11,17 +11,21 @@ module kinematics_thermal_expansion use math use lattice use material - + implicit none private - + + integer, dimension(:), allocatable :: kinematics_thermal_expansion_instance + type :: tParameters - real(pReal), allocatable, dimension(:,:,:) :: & - expansion + real(pReal) :: & + T_ref + real(pReal), dimension(3,3,3) :: & + expansion = 0.0_pReal end type tParameters - + type(tParameters), dimension(:), allocatable :: param - + public :: & kinematics_thermal_expansion_init, & kinematics_thermal_expansion_initialStrain, & @@ -35,36 +39,40 @@ contains !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- subroutine kinematics_thermal_expansion_init - - integer :: & - Ninstance, & - p, i - real(pReal), dimension(:), allocatable :: & - temp - + + integer :: Ninstance,p,i + real(pReal), dimension(:), allocatable :: temp + write(6,'(/,a)') ' <<<+- kinematics_'//KINEMATICS_thermal_expansion_LABEL//' init -+>>>'; flush(6) - + Ninstance = count(phase_kinematics == KINEMATICS_thermal_expansion_ID) - if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) & write(6,'(a16,1x,i5,/)') '# instances:',Ninstance - + + allocate(kinematics_thermal_expansion_instance(size(config_phase)), source=0) allocate(param(Ninstance)) - - do p = 1, size(phase_kinematics) + + do p = 1, size(config_phase) + kinematics_thermal_expansion_instance(p) = count(phase_kinematics(:,1:p) == KINEMATICS_thermal_expansion_ID) if (all(phase_kinematics(:,p) /= KINEMATICS_thermal_expansion_ID)) cycle - - ! ToDo: Here we need to decide how to extend the concept of instances to - ! kinetics and sources. I would suggest that the same mechanism exists at maximum once per phase - + + associate(prm => param(kinematics_thermal_expansion_instance(p)), & + config => config_phase(p)) + + prm%T_ref = config%getFloat('reference_temperature', defaultVal=0.0_pReal) + ! read up to three parameters (constant, linear, quadratic with T) - temp = config_phase(p)%getFloats('thermal_expansion11') - !lattice_thermalExpansion33(1,1,1:size(temp),p) = temp - temp = config_phase(p)%getFloats('thermal_expansion22', & - defaultVal=[(0.0_pReal, i=1,size(temp))],requiredSize=size(temp)) - !lattice_thermalExpansion33(2,2,1:size(temp),p) = temp - temp = config_phase(p)%getFloats('thermal_expansion33', & - defaultVal=[(0.0_pReal, i=1,size(temp))],requiredSize=size(temp)) + temp = config%getFloats('thermal_expansion11') + prm%expansion(1,1,1:size(temp)) = temp + temp = config%getFloats('thermal_expansion22',defaultVal=[(0.0_pReal, i=1,size(temp))],requiredSize=size(temp)) + prm%expansion(2,2,1:size(temp)) = temp + temp = config%getFloats('thermal_expansion33',defaultVal=[(0.0_pReal, i=1,size(temp))],requiredSize=size(temp)) + prm%expansion(3,3,1:size(temp)) = temp + do i=1, size(prm%expansion,3) + prm%expansion(1:3,1:3,i) = lattice_applyLatticeSymmetry33(prm%expansion(1:3,1:3,i),config%getString('lattice_structure')) + enddo + + end associate enddo end subroutine kinematics_thermal_expansion_init @@ -74,30 +82,30 @@ end subroutine kinematics_thermal_expansion_init !> @brief report initial thermal strain based on current temperature deviation from reference !-------------------------------------------------------------------------------------------------- pure function kinematics_thermal_expansion_initialStrain(homog,phase,offset) - + integer, intent(in) :: & phase, & - homog, offset + homog, & + offset + real(pReal), dimension(3,3) :: & - kinematics_thermal_expansion_initialStrain !< initial thermal strain (should be small strain, though) - - + kinematics_thermal_expansion_initialStrain !< initial thermal strain (should be small strain, though) + + associate(prm => param(kinematics_thermal_expansion_instance(phase))) kinematics_thermal_expansion_initialStrain = & - (temperature(homog)%p(offset) - lattice_referenceTemperature(phase))**1 / 1. * & - lattice_thermalExpansion33(1:3,1:3,1,phase) + & ! constant coefficient - (temperature(homog)%p(offset) - lattice_referenceTemperature(phase))**2 / 2. * & - lattice_thermalExpansion33(1:3,1:3,2,phase) + & ! linear coefficient - (temperature(homog)%p(offset) - lattice_referenceTemperature(phase))**3 / 3. * & - lattice_thermalExpansion33(1:3,1:3,3,phase) ! quadratic coefficient - + (temperature(homog)%p(offset) - prm%T_ref)**1 / 1. * prm%expansion(1:3,1:3,1) + & ! constant coefficient + (temperature(homog)%p(offset) - prm%T_ref)**2 / 2. * prm%expansion(1:3,1:3,2) + & ! linear coefficient + (temperature(homog)%p(offset) - prm%T_ref)**3 / 3. * prm%expansion(1:3,1:3,3) ! quadratic coefficient + end associate + end function kinematics_thermal_expansion_initialStrain !-------------------------------------------------------------------------------------------------- -!> @brief contains the constitutive equation for calculating the velocity gradient +!> @brief contains the constitutive equation for calculating the velocity gradient !-------------------------------------------------------------------------------------------------- subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, ipc, ip, el) - + integer, intent(in) :: & ipc, & !< grain number ip, & !< integration point number @@ -106,31 +114,32 @@ subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, ipc, ip, 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) + integer :: & phase, & - homog, offset + homog real(pReal) :: & - T, TRef, TDot - + T, TDot + phase = material_phaseAt(ipc,el) homog = material_homogenizationAt(el) - offset = thermalMapping(homog)%p(ip,el) - T = temperature(homog)%p(offset) - TDot = temperatureRate(homog)%p(offset) - TRef = lattice_referenceTemperature(phase) - + T = temperature(homog)%p(thermalMapping(homog)%p(ip,el)) + TDot = temperatureRate(homog)%p(thermalMapping(homog)%p(ip,el)) + + associate(prm => param(kinematics_thermal_expansion_instance(phase))) Li = TDot * ( & - lattice_thermalExpansion33(1:3,1:3,1,phase)*(T - TRef)**0 & ! constant coefficient - + lattice_thermalExpansion33(1:3,1:3,2,phase)*(T - TRef)**1 & ! linear coefficient - + lattice_thermalExpansion33(1:3,1:3,3,phase)*(T - TRef)**2 & ! quadratic coefficient + prm%expansion(1:3,1:3,1)*(T - prm%T_ref)**0 & ! constant coefficient + + prm%expansion(1:3,1:3,2)*(T - prm%T_ref)**1 & ! linear coefficient + + prm%expansion(1:3,1:3,3)*(T - prm%T_ref)**2 & ! quadratic coefficient ) / & (1.0_pReal & - + lattice_thermalExpansion33(1:3,1:3,1,phase)*(T - TRef)**1 / 1. & - + lattice_thermalExpansion33(1:3,1:3,2,phase)*(T - TRef)**2 / 2. & - + lattice_thermalExpansion33(1:3,1:3,3,phase)*(T - TRef)**3 / 3. & + + prm%expansion(1:3,1:3,1)*(T - prm%T_ref)**1 / 1. & + + prm%expansion(1:3,1:3,2)*(T - prm%T_ref)**2 / 2. & + + prm%expansion(1:3,1:3,3)*(T - prm%T_ref)**3 / 3. & ) - dLi_dTstar = 0.0_pReal - + end associate + dLi_dTstar = 0.0_pReal + end subroutine kinematics_thermal_expansion_LiAndItsTangent end module kinematics_thermal_expansion diff --git a/src/lattice.f90 b/src/lattice.f90 index ab2f43284..3f5ceddb7 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -16,47 +16,35 @@ module lattice implicit none private -! BEGIN DEPRECATED - integer, parameter, public :: & - LATTICE_maxNcleavageFamily = 3 !< max # of transformation system families over lattice structures - - integer, allocatable, dimension(:,:), protected, public :: & - lattice_NcleavageSystem !< total # of transformation systems in each family - - real(pReal), allocatable, dimension(:,:,:,:,:), protected, public :: & - lattice_Scleavage !< Schmid matrices for cleavage systems -! END DEPRECATED - - !-------------------------------------------------------------------------------------------------- ! face centered cubic integer, dimension(2), parameter :: & - LATTICE_FCC_NSLIPSYSTEM = [12, 6] !< # of slip systems per family for fcc + FCC_NSLIPSYSTEM = [12, 6] !< # of slip systems per family for fcc integer, dimension(1), parameter :: & - LATTICE_FCC_NTWINSYSTEM = [12] !< # of twin systems per family for fcc + FCC_NTWINSYSTEM = [12] !< # of twin systems per family for fcc integer, dimension(1), parameter :: & - LATTICE_FCC_NTRANSSYSTEM = [12] !< # of transformation systems per family for fcc + FCC_NTRANSSYSTEM = [12] !< # of transformation systems per family for fcc - integer, dimension(2), parameter :: & - LATTICE_FCC_NCLEAVAGESYSTEM = [3, 4] !< # of cleavage systems per family for fcc + integer, dimension(1), parameter :: & + FCC_NCLEAVAGESYSTEM = [3] !< # of cleavage systems per family for fcc integer, parameter :: & #ifndef __PGI - LATTICE_FCC_NSLIP = sum(LATTICE_FCC_NSLIPSYSTEM), & !< total # of slip systems for fcc - LATTICE_FCC_NTWIN = sum(LATTICE_FCC_NTWINSYSTEM), & !< total # of twin systems for fcc - LATTICE_FCC_NTRANS = sum(LATTICE_FCC_NTRANSSYSTEM), & !< total # of transformation systems for fcc - LATTICE_FCC_NCLEAVAGE = sum(LATTICE_FCC_NCLEAVAGESYSTEM) !< total # of cleavage systems for fcc + FCC_NSLIP = sum(FCC_NSLIPSYSTEM), & !< total # of slip systems for fcc + FCC_NTWIN = sum(FCC_NTWINSYSTEM), & !< total # of twin systems for fcc + FCC_NTRANS = sum(FCC_NTRANSSYSTEM), & !< total # of transformation systems for fcc + FCC_NCLEAVAGE = sum(FCC_NCLEAVAGESYSTEM) !< total # of cleavage systems for fcc #else - LATTICE_FCC_NSLIP = 18, & - LATTICE_FCC_NTWIN = 12, & - LATTICE_FCC_NTRANS = 12, & - LATTICE_FCC_NCLEAVAGE = 7 + FCC_NSLIP = 18, & + FCC_NTWIN = 12, & + FCC_NTRANS = 12, & + FCC_NCLEAVAGE = 3 #endif - real(pReal), dimension(3+3,LATTICE_FCC_NSLIP), parameter :: & - LATTICE_FCC_SYSTEMSLIP = reshape(real([& + real(pReal), dimension(3+3,FCC_NSLIP), parameter :: & + FCC_SYSTEMSLIP = reshape(real([& ! Slip direction Plane normal ! SCHMID-BOAS notation 0, 1,-1, 1, 1, 1, & ! B2 -1, 0, 1, 1, 1, 1, & ! B4 @@ -77,10 +65,10 @@ module lattice 1, 0,-1, 1, 0, 1, & 0, 1, 1, 0, 1,-1, & 0, 1,-1, 0, 1, 1 & - ],pReal),shape(LATTICE_FCC_SYSTEMSLIP)) !< Slip system <110>{111} directions. Sorted according to Eisenlohr & Hantcherli + ],pReal),shape(FCC_SYSTEMSLIP)) !< Slip system <110>{111} directions. Sorted according to Eisenlohr & Hantcherli - real(pReal), dimension(3+3,LATTICE_FCC_NTWIN), parameter :: & - LATTICE_FCC_SYSTEMTWIN = reshape(real( [& + real(pReal), dimension(3+3,FCC_NTWIN), parameter :: & + FCC_SYSTEMTWIN = reshape(real( [& -2, 1, 1, 1, 1, 1, & 1,-2, 1, 1, 1, 1, & 1, 1,-2, 1, 1, 1, & @@ -93,10 +81,10 @@ module lattice 2, 1,-1, -1, 1,-1, & -1,-2,-1, -1, 1,-1, & -1, 1, 2, -1, 1,-1 & - ],pReal),shape(LATTICE_FCC_SYSTEMTWIN)) !< Twin system <112>{111} directions. Sorted according to Eisenlohr & Hantcherli + ],pReal),shape(FCC_SYSTEMTWIN)) !< Twin system <112>{111} directions. Sorted according to Eisenlohr & Hantcherli - integer, dimension(2,LATTICE_FCC_NTWIN), parameter, public :: & - LATTICE_FCC_TWINNUCLEATIONSLIPPAIR = reshape( [& + integer, dimension(2,FCC_NTWIN), parameter, public :: & + lattice_FCC_TWINNUCLEATIONSLIPPAIR = reshape( [& 2,3, & 1,3, & 1,2, & @@ -109,44 +97,40 @@ module lattice 11,12, & 10,12, & 10,11 & - ],shape(LATTICE_FCC_TWINNUCLEATIONSLIPPAIR)) + ],shape(lattice_FCC_TWINNUCLEATIONSLIPPAIR)) - real(pReal), dimension(3+3,LATTICE_FCC_NCLEAVAGE), parameter :: & - LATTICE_FCC_SYSTEMCLEAVAGE = reshape(real([& + real(pReal), dimension(3+3,FCC_NCLEAVAGE), parameter :: & + FCC_SYSTEMCLEAVAGE = reshape(real([& ! Cleavage direction Plane normal 0, 1, 0, 1, 0, 0, & 0, 0, 1, 0, 1, 0, & - 1, 0, 0, 0, 0, 1, & - 0, 1,-1, 1, 1, 1, & - 0,-1,-1, -1,-1, 1, & - -1, 0,-1, 1,-1,-1, & - 0, 1, 1, -1, 1,-1 & - ],pReal),shape(LATTICE_FCC_SYSTEMCLEAVAGE)) + 1, 0, 0, 0, 0, 1 & + ],pReal),shape(FCC_SYSTEMCLEAVAGE)) !-------------------------------------------------------------------------------------------------- ! body centered cubic integer, dimension(2), parameter :: & - LATTICE_BCC_NSLIPSYSTEM = [12, 12] !< # of slip systems per family for bcc + BCC_NSLIPSYSTEM = [12, 12] !< # of slip systems per family for bcc integer, dimension(1), parameter :: & - LATTICE_BCC_NTWINSYSTEM = [12] !< # of twin systems per family for bcc + BCC_NTWINSYSTEM = [12] !< # of twin systems per family for bcc - integer, dimension(2), parameter :: & - LATTICE_BCC_NCLEAVAGESYSTEM = [3, 6] !< # of cleavage systems per family for bcc + integer, dimension(1), parameter :: & + BCC_NCLEAVAGESYSTEM = [3] !< # of cleavage systems per family for bcc integer, parameter :: & #ifndef __PGI - LATTICE_BCC_NSLIP = sum(LATTICE_BCC_NSLIPSYSTEM), & !< total # of slip systems for bcc - LATTICE_BCC_NTWIN = sum(LATTICE_BCC_NTWINSYSTEM), & !< total # of twin systems for bcc - LATTICE_BCC_NCLEAVAGE = sum(LATTICE_BCC_NCLEAVAGESYSTEM) !< total # of cleavage systems for bcc + BCC_NSLIP = sum(BCC_NSLIPSYSTEM), & !< total # of slip systems for bcc + BCC_NTWIN = sum(BCC_NTWINSYSTEM), & !< total # of twin systems for bcc + BCC_NCLEAVAGE = sum(BCC_NCLEAVAGESYSTEM) !< total # of cleavage systems for bcc #else - LATTICE_BCC_NSLIP = 24, & - LATTICE_BCC_NTWIN = 12, & - LATTICE_BCC_NCLEAVAGE = 9 + BCC_NSLIP = 24, & + BCC_NTWIN = 12, & + BCC_NCLEAVAGE = 3 #endif - real(pReal), dimension(3+3,LATTICE_BCC_NSLIP), parameter :: & - LATTICE_BCC_SYSTEMSLIP = reshape(real([& + real(pReal), dimension(3+3,BCC_NSLIP), parameter :: & + BCC_SYSTEMSLIP = reshape(real([& ! Slip direction Plane normal ! Slip system <111>{110} 1,-1, 1, 0, 1, 1, & @@ -174,10 +158,10 @@ module lattice 1,-1, 1, -1, 1, 2, & -1, 1, 1, 1,-1, 2, & 1, 1, 1, 1, 1,-2 & - ],pReal),shape(LATTICE_BCC_SYSTEMSLIP)) + ],pReal),shape(BCC_SYSTEMSLIP)) - real(pReal), dimension(3+3,LATTICE_BCC_NTWIN), parameter :: & - LATTICE_BCC_SYSTEMTWIN = reshape(real([& + real(pReal), dimension(3+3,BCC_NTWIN), parameter :: & + BCC_SYSTEMTWIN = reshape(real([& ! Twin system <111>{112} -1, 1, 1, 2, 1, 1, & 1, 1, 1, -2, 1, 1, & @@ -191,46 +175,35 @@ module lattice 1,-1, 1, -1, 1, 2, & -1, 1, 1, 1,-1, 2, & 1, 1, 1, 1, 1,-2 & - ],pReal),shape(LATTICE_BCC_SYSTEMTWIN)) + ],pReal),shape(BCC_SYSTEMTWIN)) - real(pReal), dimension(3+3,LATTICE_BCC_NCLEAVAGE), parameter :: & - LATTICE_BCC_SYSTEMCLEAVAGE = reshape(real([& + real(pReal), dimension(3+3,BCC_NCLEAVAGE), parameter :: & + BCC_SYSTEMCLEAVAGE = reshape(real([& ! Cleavage direction Plane normal 0, 1, 0, 1, 0, 0, & 0, 0, 1, 0, 1, 0, & - 1, 0, 0, 0, 0, 1, & - 1,-1, 1, 0, 1, 1, & - 1, 1, 1, 0,-1, 1, & - -1, 1, 1, 1, 0, 1, & - 1, 1, 1, -1, 0, 1, & - -1, 1, 1, 1, 1, 0, & - 1, 1, 1, -1, 1, 0 & - ],pReal),shape(LATTICE_BCC_SYSTEMCLEAVAGE)) + 1, 0, 0, 0, 0, 1 & + ],pReal),shape(BCC_SYSTEMCLEAVAGE)) !-------------------------------------------------------------------------------------------------- ! hexagonal integer, dimension(6), parameter :: & - LATTICE_HEX_NSLIPSYSTEM = [3, 3, 3, 6, 12, 6] !< # of slip systems per family for hex + HEX_NSLIPSYSTEM = [3, 3, 3, 6, 12, 6] !< # of slip systems per family for hex integer, dimension(4), parameter :: & - LATTICE_HEX_NTWINSYSTEM = [6, 6, 6, 6] !< # of slip systems per family for hex - - integer, dimension(1), parameter :: & - LATTICE_HEX_NCLEAVAGESYSTEM = [3] !< # of cleavage systems per family for hex + HEX_NTWINSYSTEM = [6, 6, 6, 6] !< # of slip systems per family for hex integer, parameter :: & #ifndef __PGI - LATTICE_HEX_NSLIP = sum(LATTICE_HEX_NSLIPSYSTEM), & !< total # of slip systems for hex - LATTICE_HEX_NTWIN = sum(LATTICE_HEX_NTWINSYSTEM), & !< total # of twin systems for hex - LATTICE_HEX_NCLEAVAGE = sum(LATTICE_HEX_NCLEAVAGESYSTEM) !< total # of cleavage systems for hex + HEX_NSLIP = sum(HEX_NSLIPSYSTEM), & !< total # of slip systems for hex + HEX_NTWIN = sum(HEX_NTWINSYSTEM) !< total # of twin systems for hex #else - LATTICE_HEX_NSLIP = 33, & - LATTICE_HEX_NTWIN = 24, & - LATTICE_HEX_NCLEAVAGE = 3 + HEX_NSLIP = 33, & + HEX_NTWIN = 24 #endif - real(pReal), dimension(4+4,LATTICE_HEX_NSLIP), parameter :: & - LATTICE_HEX_SYSTEMSLIP = reshape(real([& + real(pReal), dimension(4+4,HEX_NSLIP), parameter :: & + HEX_SYSTEMSLIP = reshape(real([& ! Slip direction Plane normal ! Basal systems <-1-1.0>{00.1} (independent of c/a-ratio, Bravais notation (4 coordinate base)) 2, -1, -1, 0, 0, 0, 0, 1, & @@ -271,10 +244,10 @@ module lattice 1, 1, -2, 3, -1, -1, 2, 2, & -1, 2, -1, 3, 1, -2, 1, 2, & -2, 1, 1, 3, 2, -1, -1, 2 & - ],pReal),shape(LATTICE_HEX_SYSTEMSLIP)) !< slip systems for hex, sorted by P. Eisenlohr CCW around starting next to a_1 axis + ],pReal),shape(HEX_SYSTEMSLIP)) !< slip systems for hex, sorted by P. Eisenlohr CCW around starting next to a_1 axis - real(pReal), dimension(4+4,LATTICE_HEX_NTWIN), parameter :: & - LATTICE_HEX_SYSTEMTWIN = reshape(real([& + real(pReal), dimension(4+4,HEX_NTWIN), parameter :: & + HEX_SYSTEMTWIN = reshape(real([& ! Compression or Tension = f(twinning shear=f(c/a)) for each metal ! (according to Yoo 1981) -1, 0, 1, 1, 1, 0, -1, 2, & ! <-10.1>{10.2} shear = (3-(c/a)^2)/(sqrt(3) c/a) 0, -1, 1, 1, 0, 1, -1, 2, & @@ -303,31 +276,22 @@ module lattice -1, -1, 2, -3, -1, -1, 2, 2, & 1, -2, 1, -3, 1, -2, 1, 2, & 2, -1, -1, -3, 2, -1, -1, 2 & - ],pReal),shape(LATTICE_HEX_SYSTEMTWIN)) !< twin systems for hex, sorted by P. Eisenlohr CCW around starting next to a_1 axis - - real(pReal), dimension(4+4,LATTICE_HEX_NCLEAVAGE), parameter :: & - LATTICE_HEX_SYSTEMCLEAVAGE = reshape(real([& - ! Cleavage direction Plane normal - 2,-1,-1, 0, 0, 0, 0, 1, & - 0, 0, 0, 1, 2,-1,-1, 0, & - 0, 0, 0, 1, 0, 1,-1, 0 & - ],pReal),shape(LATTICE_HEX_SYSTEMCLEAVAGE)) - + ],pReal),shape(HEX_SYSTEMTWIN)) !< twin systems for hex, sorted by P. Eisenlohr CCW around starting next to a_1 axis !-------------------------------------------------------------------------------------------------- ! body centered tetragonal integer, dimension(13), parameter :: & - LATTICE_BCT_NSLIPSYSTEM = [2, 2, 2, 4, 2, 4, 2, 2, 4, 8, 4, 8, 8 ] !< # of slip systems per family for bct (Sn) Bieler J. Electr Mater 2009 + BCT_NSLIPSYSTEM = [2, 2, 2, 4, 2, 4, 2, 2, 4, 8, 4, 8, 8 ] !< # of slip systems per family for bct (Sn) Bieler J. Electr Mater 2009 integer, parameter :: & #ifndef __PGI - LATTICE_BCT_NSLIP = sum(LATTICE_BCT_NSLIPSYSTEM) !< total # of slip systems for bct + BCT_NSLIP = sum(BCT_NSLIPSYSTEM) !< total # of slip systems for bct #else - LATTICE_BCT_NSLIP = 52 + BCT_NSLIP = 52 #endif - real(pReal), dimension(3+3,LATTICE_BCT_NSLIP), parameter :: & - LATTICE_BCT_SYSTEMSLIP = reshape(real([& + real(pReal), dimension(3+3,BCT_NSLIP), parameter :: & + BCT_SYSTEMSLIP = reshape(real([& ! Slip direction Plane normal ! Slip family 1 {100)<001] (Bravais notation {hkl)>>' + write(6,'(/,a)') ' <<<+- lattice init -+>>>'; flush(6) Nphases = size(config_phase) - allocate(lattice_structure(Nphases),source = LATTICE_undefined_ID) + allocate(lattice_structure(Nphases),source = lattice_UNDEFINED_ID) allocate(lattice_C66(6,6,Nphases), source=0.0_pReal) - allocate(lattice_C3333(3,3,3,3,Nphases), source=0.0_pReal) - allocate(lattice_thermalExpansion33 (3,3,3,Nphases), source=0.0_pReal) ! constant, linear, quadratic coefficients - allocate(lattice_thermalConductivity33 (3,3,Nphases), source=0.0_pReal) - allocate(lattice_damageDiffusion33 (3,3,Nphases), source=0.0_pReal) - allocate(lattice_damageMobility ( Nphases), source=0.0_pReal) - allocate(lattice_massDensity ( Nphases), source=0.0_pReal) - allocate(lattice_specificHeat ( Nphases), source=0.0_pReal) - allocate(lattice_referenceTemperature ( Nphases), source=300.0_pReal) - - allocate(lattice_mu(Nphases), source=0.0_pReal) - allocate(lattice_nu(Nphases), source=0.0_pReal) - - allocate(lattice_Scleavage(3,3,3,lattice_maxNcleavage,Nphases),source=0.0_pReal) - allocate(lattice_NcleavageSystem(lattice_maxNcleavageFamily,Nphases),source=0) + allocate(lattice_thermalConductivity (3,3,Nphases), source=0.0_pReal) + allocate(lattice_damageDiffusion (3,3,Nphases), source=0.0_pReal) + allocate(lattice_damageMobility,& + lattice_massDensity,lattice_specificHeat, & + lattice_mu, lattice_nu,& + source=[(0.0_pReal,i=1,Nphases)]) do p = 1, size(config_phase) - lattice_C66(1,1,p) = config_phase(p)%getFloat('c11',defaultVal=0.0_pReal) - lattice_C66(1,2,p) = config_phase(p)%getFloat('c12',defaultVal=0.0_pReal) + lattice_C66(1,1,p) = config_phase(p)%getFloat('c11') + lattice_C66(1,2,p) = config_phase(p)%getFloat('c12') + lattice_C66(1,3,p) = config_phase(p)%getFloat('c13',defaultVal=0.0_pReal) lattice_C66(2,2,p) = config_phase(p)%getFloat('c22',defaultVal=0.0_pReal) lattice_C66(2,3,p) = config_phase(p)%getFloat('c23',defaultVal=0.0_pReal) @@ -566,36 +482,34 @@ subroutine lattice_init lattice_C66(5,5,p) = config_phase(p)%getFloat('c55',defaultVal=0.0_pReal) lattice_C66(6,6,p) = config_phase(p)%getFloat('c66',defaultVal=0.0_pReal) - CoverA = config_phase(p)%getFloat('c/a',defaultVal=0.0_pReal) - - tag = config_phase(p)%getString('lattice_structure') - select case(tag(1:3)) + structure = config_phase(p)%getString('lattice_structure') + select case(trim(structure)) case('iso') - lattice_structure(p) = LATTICE_iso_ID + lattice_structure(p) = lattice_ISO_ID case('fcc') - lattice_structure(p) = LATTICE_fcc_ID + lattice_structure(p) = lattice_FCC_ID case('bcc') - lattice_structure(p) = LATTICE_bcc_ID + lattice_structure(p) = lattice_BCC_ID case('hex') - if(CoverA < 1.0_pReal .or. CoverA > 2.0_pReal) call IO_error(131,el=p) - lattice_structure(p) = LATTICE_hex_ID + lattice_structure(p) = lattice_HEX_ID case('bct') - if(CoverA > 2.0_pReal) call IO_error(131,el=p) - lattice_structure(p) = LATTICE_bct_ID + lattice_structure(p) = lattice_BCT_ID case('ort') - lattice_structure(p) = LATTICE_ort_ID + lattice_structure(p) = lattice_ORT_ID case default - call IO_error(130,ext_msg='lattice_init') + call IO_error(130,ext_msg='lattice_init: '//trim(structure)) end select - lattice_C66(1:6,1:6,p) = lattice_symmetrizeC66(lattice_structure(p),lattice_C66(1:6,1:6,p)) + lattice_C66(1:6,1:6,p) = applyLatticeSymmetryC66(lattice_C66(1:6,1:6,p),structure) - lattice_mu(p) = 0.2_pReal *(lattice_C66(1,1,p) -lattice_C66(1,2,p) +3.0_pReal*lattice_C66(4,4,p)) ! (C11iso-C12iso)/2 with C11iso=(3*C11+2*C12+4*C44)/5 and C12iso=(C11+4*C12-2*C44)/5 + ! (C11iso-C12iso)/2 with C11iso=(3*C11+2*C12+4*C44)/5 and C12iso=(C11+4*C12-2*C44)/5 + lattice_mu(p) = 0.2_pReal *(lattice_C66(1,1,p) -lattice_C66(1,2,p) +3.0_pReal*lattice_C66(4,4,p)) + + ! C12iso/(C11iso+C12iso) with C11iso=(3*C11+2*C12+4*C44)/5 and C12iso=(C11+4*C12-2*C44)/5 lattice_nu(p) = ( lattice_C66(1,1,p) +4.0_pReal*lattice_C66(1,2,p) -2.0_pReal*lattice_C66(4,4,p)) & - / (4.0_pReal*lattice_C66(1,1,p) +6.0_pReal*lattice_C66(1,2,p) +2.0_pReal*lattice_C66(4,4,p))! C12iso/(C11iso+C12iso) with C11iso=(3*C11+2*C12+4*C44)/5 and C12iso=(C11+4*C12-2*C44)/5 + / (4.0_pReal*lattice_C66(1,1,p) +6.0_pReal*lattice_C66(1,2,p) +2.0_pReal*lattice_C66(4,4,p)) - lattice_C3333(1:3,1:3,1:3,1:3,p) = math_Voigt66to3333(lattice_C66(1:6,1:6,p)) ! Literature data is Voigt - lattice_C66(1:6,1:6,p) = math_sym3333to66(lattice_C3333(1:3,1:3,1:3,1:3,p)) ! DAMASK uses Mandel-weighting + lattice_C66(1:6,1:6,p) = math_sym3333to66(math_Voigt66to3333(lattice_C66(1:6,1:6,p))) ! Literature data is in Voigt notation do i = 1, 6 if (abs(lattice_C66(i,i,p)) @brief !!!!!!!DEPRECTATED!!!!!! -!-------------------------------------------------------------------------------------------------- -subroutine lattice_initializeStructure(myPhase,CoverA) - - integer, intent(in) :: myPhase - real(pReal), intent(in) :: & - CoverA - - integer :: & - i - - do i = 1,3 - lattice_thermalExpansion33 (1:3,1:3,i,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),& - lattice_thermalExpansion33 (1:3,1:3,i,myPhase)) - enddo - - lattice_thermalConductivity33 (1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),& - lattice_thermalConductivity33 (1:3,1:3,myPhase)) - lattice_DamageDiffusion33 (1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),& - lattice_DamageDiffusion33 (1:3,1:3,myPhase)) - - select case(lattice_structure(myPhase)) - - case (LATTICE_fcc_ID) - lattice_NcleavageSystem(1:2,myPhase) = lattice_fcc_NcleavageSystem - lattice_Scleavage(1:3,1:3,1:3,1:lattice_fcc_Ncleavage,myPhase) = & - lattice_SchmidMatrix_cleavage(lattice_fcc_ncleavageSystem,'fcc',covera) - - case (LATTICE_bcc_ID) - lattice_NcleavageSystem(1:2,myPhase) = lattice_bcc_NcleavageSystem - lattice_Scleavage(1:3,1:3,1:3,1:lattice_bcc_Ncleavage,myPhase) = & - lattice_SchmidMatrix_cleavage(lattice_bcc_ncleavagesystem,'bcc',covera) - - case (LATTICE_hex_ID) - lattice_NcleavageSystem(1:1,myPhase) = lattice_hex_NcleavageSystem - lattice_Scleavage(1:3,1:3,1:3,1:lattice_hex_Ncleavage,myPhase) = & - lattice_SchmidMatrix_cleavage(lattice_hex_ncleavagesystem,'hex',covera) - - case (LATTICE_ort_ID) - lattice_NcleavageSystem(1:3,myPhase) = lattice_ort_NcleavageSystem - lattice_Scleavage(1:3,1:3,1:3,1:lattice_ort_Ncleavage,myPhase) = & - lattice_SchmidMatrix_cleavage(lattice_ort_NcleavageSystem,'ort',covera) - - case (LATTICE_iso_ID) - lattice_NcleavageSystem(1:1,myPhase) = lattice_iso_NcleavageSystem - lattice_Scleavage(1:3,1:3,1:3,1:lattice_iso_Ncleavage,myPhase) = & - lattice_SchmidMatrix_cleavage(lattice_iso_NcleavageSystem,'iso',covera) - - end select - -end subroutine lattice_initializeStructure - - -!-------------------------------------------------------------------------------------------------- -!> @brief Symmetrizes stiffness matrix according to lattice type -!> @details J. A. Rayne and B. S. Chandrasekhar Phys. Rev. 120, 1658 Erratum Phys. Rev. 122, 1962 -!-------------------------------------------------------------------------------------------------- -pure function lattice_symmetrizeC66(struct,C66) - - integer(kind(LATTICE_undefined_ID)), intent(in) :: struct - real(pReal), dimension(6,6), intent(in) :: C66 - real(pReal), dimension(6,6) :: lattice_symmetrizeC66 - integer :: j,k - - lattice_symmetrizeC66 = 0.0_pReal - - select case(struct) - case (LATTICE_iso_ID) - do k=1,3 - do j=1,3 - lattice_symmetrizeC66(k,j) = C66(1,2) - enddo - lattice_symmetrizeC66(k,k) = C66(1,1) - lattice_symmetrizeC66(k+3,k+3) = 0.5_pReal*(C66(1,1)-C66(1,2)) - enddo - case (LATTICE_fcc_ID,LATTICE_bcc_ID) - do k=1,3 - do j=1,3 - lattice_symmetrizeC66(k,j) = C66(1,2) - enddo - lattice_symmetrizeC66(k,k) = C66(1,1) - lattice_symmetrizeC66(k+3,k+3) = C66(4,4) - enddo - case (LATTICE_hex_ID) - lattice_symmetrizeC66(1,1) = C66(1,1) - lattice_symmetrizeC66(2,2) = C66(1,1) - lattice_symmetrizeC66(3,3) = C66(3,3) - lattice_symmetrizeC66(1,2) = C66(1,2) - lattice_symmetrizeC66(2,1) = C66(1,2) - lattice_symmetrizeC66(1,3) = C66(1,3) - lattice_symmetrizeC66(3,1) = C66(1,3) - lattice_symmetrizeC66(2,3) = C66(1,3) - lattice_symmetrizeC66(3,2) = C66(1,3) - lattice_symmetrizeC66(4,4) = C66(4,4) - lattice_symmetrizeC66(5,5) = C66(4,4) - lattice_symmetrizeC66(6,6) = 0.5_pReal*(C66(1,1)-C66(1,2)) - case (LATTICE_ort_ID) - lattice_symmetrizeC66(1,1) = C66(1,1) - lattice_symmetrizeC66(2,2) = C66(2,2) - lattice_symmetrizeC66(3,3) = C66(3,3) - lattice_symmetrizeC66(1,2) = C66(1,2) - lattice_symmetrizeC66(2,1) = C66(1,2) - lattice_symmetrizeC66(1,3) = C66(1,3) - lattice_symmetrizeC66(3,1) = C66(1,3) - lattice_symmetrizeC66(2,3) = C66(2,3) - lattice_symmetrizeC66(3,2) = C66(2,3) - lattice_symmetrizeC66(4,4) = C66(4,4) - lattice_symmetrizeC66(5,5) = C66(5,5) - lattice_symmetrizeC66(6,6) = C66(6,6) - case (LATTICE_bct_ID) - lattice_symmetrizeC66(1,1) = C66(1,1) - lattice_symmetrizeC66(2,2) = C66(1,1) - lattice_symmetrizeC66(3,3) = C66(3,3) - lattice_symmetrizeC66(1,2) = C66(1,2) - lattice_symmetrizeC66(2,1) = C66(1,2) - lattice_symmetrizeC66(1,3) = C66(1,3) - lattice_symmetrizeC66(3,1) = C66(1,3) - lattice_symmetrizeC66(2,3) = C66(1,3) - lattice_symmetrizeC66(3,2) = C66(1,3) - lattice_symmetrizeC66(4,4) = C66(4,4) - lattice_symmetrizeC66(5,5) = C66(4,4) - lattice_symmetrizeC66(6,6) = C66(6,6) - case default - lattice_symmetrizeC66 = C66 - end select - -end function lattice_symmetrizeC66 - - -!-------------------------------------------------------------------------------------------------- -!> @brief Symmetrizes 2nd order tensor according to lattice type -!-------------------------------------------------------------------------------------------------- -pure function lattice_symmetrize33(struct,T33) - - integer(kind(LATTICE_undefined_ID)), intent(in) :: struct - real(pReal), dimension(3,3), intent(in) :: T33 - real(pReal), dimension(3,3) :: lattice_symmetrize33 - integer :: k - - lattice_symmetrize33 = 0.0_pReal - - select case(struct) - case (LATTICE_iso_ID,LATTICE_fcc_ID,LATTICE_bcc_ID) - do k=1,3 - lattice_symmetrize33(k,k) = T33(1,1) - enddo - case (LATTICE_hex_ID) - lattice_symmetrize33(1,1) = T33(1,1) - lattice_symmetrize33(2,2) = T33(1,1) - lattice_symmetrize33(3,3) = T33(3,3) - case (LATTICE_ort_ID,lattice_bct_ID) - lattice_symmetrize33(1,1) = T33(1,1) - lattice_symmetrize33(2,2) = T33(2,2) - lattice_symmetrize33(3,3) = T33(3,3) - case default - lattice_symmetrize33 = T33 - end select - -end function lattice_symmetrize33 - - !-------------------------------------------------------------------------------------------------- !> @brief Characteristic shear for twinning !-------------------------------------------------------------------------------------------------- @@ -807,7 +555,7 @@ function lattice_characteristicShear_Twin(Ntwin,structure,CoverA) result(charact f, & !< index of my family s !< index of my system in current family - integer, dimension(LATTICE_HEX_NTWIN), parameter :: & + integer, dimension(HEX_NTWIN), parameter :: & HEX_SHEARTWIN = reshape( [& 1, & ! <-10.1>{10.2} 1, & @@ -833,7 +581,7 @@ function lattice_characteristicShear_Twin(Ntwin,structure,CoverA) result(charact 4, & 4, & 4 & - ],[LATTICE_HEX_NTWIN]) ! indicator to formulas below + ],[HEX_NTWIN]) ! indicator to formulas below if (len_trim(structure) /= 3) & call IO_error(137,ext_msg='lattice_characteristicShear_Twin: '//trim(structure)) @@ -842,13 +590,13 @@ function lattice_characteristicShear_Twin(Ntwin,structure,CoverA) result(charact myFamilies: do f = 1,size(Ntwin,1) mySystems: do s = 1,Ntwin(f) a = a + 1 - select case(structure(1:3)) + select case(structure) case('fcc','bcc') characteristicShear(a) = 0.5_pReal*sqrt(2.0_pReal) case('hex') if (cOverA < 1.0_pReal .or. cOverA > 2.0_pReal) & call IO_error(131,ext_msg='lattice_characteristicShear_Twin') - p = sum(LATTICE_HEX_NTWINSYSTEM(1:f-1))+s + p = sum(HEX_NTWINSYSTEM(1:f-1))+s select case(HEX_SHEARTWIN(p)) ! from Christian & Mahajan 1995 p.29 case (1) ! <-10.1>{10.2} characteristicShear(a) = (3.0_pReal-cOverA**2.0_pReal)/sqrt(3.0_pReal)/CoverA @@ -886,15 +634,15 @@ function lattice_C66_twin(Ntwin,C66,structure,CoverA) if (len_trim(structure) /= 3) & call IO_error(137,ext_msg='lattice_C66_twin: '//trim(structure)) - select case(structure(1:3)) + select case(structure) case('fcc') - coordinateSystem = buildCoordinateSystem(Ntwin,LATTICE_FCC_NSLIPSYSTEM,LATTICE_FCC_SYSTEMTWIN,& + coordinateSystem = buildCoordinateSystem(Ntwin,FCC_NSLIPSYSTEM,FCC_SYSTEMTWIN,& trim(structure),0.0_pReal) case('bcc') - coordinateSystem = buildCoordinateSystem(Ntwin,LATTICE_BCC_NSLIPSYSTEM,LATTICE_BCC_SYSTEMTWIN,& + coordinateSystem = buildCoordinateSystem(Ntwin,BCC_NSLIPSYSTEM,BCC_SYSTEMTWIN,& trim(structure),0.0_pReal) case('hex') - coordinateSystem = buildCoordinateSystem(Ntwin,LATTICE_HEX_NSLIPSYSTEM,LATTICE_HEX_SYSTEMTWIN,& + coordinateSystem = buildCoordinateSystem(Ntwin,HEX_NSLIPSYSTEM,HEX_SYSTEMTWIN,& 'hex',cOverA) case default call IO_error(137,ext_msg='lattice_C66_twin: '//trim(structure)) @@ -946,7 +694,7 @@ function lattice_C66_trans(Ntrans,C_parent66,structure_target, & C_target_unrotated66(1,3) = C_bar66(1,3) C_target_unrotated66(3,3) = C_bar66(3,3) C_target_unrotated66(4,4) = C_bar66(4,4) - C_bar66(1,4)**2.0_pReal/(0.5_pReal*(C_bar66(1,1) - C_bar66(1,2))) - C_target_unrotated66 = lattice_symmetrizeC66(LATTICE_HEX_ID,C_target_unrotated66) + C_target_unrotated66 = applyLatticeSymmetryC66(C_target_unrotated66,'hex') elseif (structure_target(1:3) == 'bcc') then if (a_bcc <= 0.0_pReal .or. a_fcc <= 0.0_pReal) & call IO_error(134,ext_msg='lattice_C66_trans: '//trim(structure_target)) @@ -989,7 +737,7 @@ function lattice_nonSchmidMatrix(Nslip,nonSchmidCoefficients,sense) result(nonSc if (abs(sense) /= 1) call IO_error(0,ext_msg='lattice_nonSchmidMatrix') - coordinateSystem = buildCoordinateSystem(Nslip,LATTICE_BCC_NSLIPSYSTEM,LATTICE_BCC_SYSTEMSLIP,& + coordinateSystem = buildCoordinateSystem(Nslip,BCC_NSLIPSYSTEM,BCC_SYSTEMSLIP,& 'bcc',0.0_pReal) coordinateSystem(1:3,1,1:sum(Nslip)) = coordinateSystem(1:3,1,1:sum(Nslip)) *real(sense,pReal) ! convert unidirectional coordinate system nonSchmidMatrix = lattice_SchmidMatrix_slip(Nslip,'bcc',0.0_pReal) ! Schmid contribution @@ -1032,7 +780,7 @@ function lattice_interaction_SlipBySlip(Nslip,interactionValues,structure) resul integer, dimension(:), allocatable :: NslipMax integer, dimension(:,:), allocatable :: interactionTypes - integer, dimension(LATTICE_FCC_NSLIP,LATTICE_FCC_NSLIP), parameter :: & + integer, dimension(FCC_NSLIP,FCC_NSLIP), parameter :: & FCC_INTERACTIONSLIPSLIP = reshape( [& 1, 2, 2, 4, 6, 5, 3, 5, 5, 4, 5, 6, 9,10, 9,10,11,12, & ! -----> acting 2, 1, 2, 6, 4, 5, 5, 4, 6, 5, 3, 5, 9,10,11,12, 9,10, & ! | @@ -1067,7 +815,7 @@ function lattice_interaction_SlipBySlip(Nslip,interactionValues,structure) resul !<11: crossing btw one {110} and one {111} plane !<12: collinear btw one {110} and one {111} plane - integer, dimension(LATTICE_BCC_NSLIP,LATTICE_BCC_NSLIP), parameter :: & + integer, dimension(BCC_NSLIP,BCC_NSLIP), parameter :: & BCC_INTERACTIONSLIPSLIP = reshape( [& 1,2,6,6,5,4,4,3,4,3,5,4, 6,6,4,3,3,4,6,6,4,3,6,6, & ! -----> acting 2,1,6,6,4,3,5,4,5,4,4,3, 6,6,3,4,4,3,6,6,3,4,6,6, & ! | @@ -1102,7 +850,7 @@ function lattice_interaction_SlipBySlip(Nslip,interactionValues,structure) resul !< 5: mixed-symmetrical junction !< 6: edge junction - integer, dimension(LATTICE_HEX_NSLIP,LATTICE_HEX_NSLIP), parameter :: & + integer, dimension(HEX_NSLIP,HEX_NSLIP), parameter :: & HEX_INTERACTIONSLIPSLIP = reshape( [& 1, 2, 2, 3, 3, 3, 7, 7, 7, 13,13,13,13,13,13, 21,21,21,21,21,21,21,21,21,21,21,21, 31,31,31,31,31,31, & ! -----> acting 2, 1, 2, 3, 3, 3, 7, 7, 7, 13,13,13,13,13,13, 21,21,21,21,21,21,21,21,21,21,21,21, 31,31,31,31,31,31, & ! | @@ -1144,7 +892,7 @@ function lattice_interaction_SlipBySlip(Nslip,interactionValues,structure) resul 42,42,42, 41,41,41, 40,40,40, 39,39,39,39,39,39, 38,38,38,38,38,38,38,38,38,38,38,38, 37,37,37,37,37,36 & ],shape(HEX_INTERACTIONSLIPSLIP)) !< Slip--slip interaction types for hex (onion peel naming scheme) - integer, dimension(LATTICE_BCT_NSLIP,LATTICE_BCT_NSLIP), parameter :: & + integer, dimension(BCT_NSLIP,BCT_NSLIP), parameter :: & BCT_INTERACTIONSLIPSLIP = reshape( [& 1, 2, 3, 3, 7, 7, 13, 13, 13, 13, 21, 21, 31, 31, 31, 31, 43, 43, 57, 57, 73, 73, 73, 73, 91, 91, 91, 91, 91, 91, 91, 91, 111, 111, 111, 111, 133,133,133,133,133,133,133,133, 157,157,157,157,157,157,157,157, & ! -----> acting 2, 1, 3, 3, 7, 7, 13, 13, 13, 13, 21, 21, 31, 31, 31, 31, 43, 43, 57, 57, 73, 73, 73, 73, 91, 91, 91, 91, 91, 91, 91, 91, 111, 111, 111, 111, 133,133,133,133,133,133,133,133, 157,157,157,157,157,157,157,157, & ! | @@ -1216,19 +964,19 @@ function lattice_interaction_SlipBySlip(Nslip,interactionValues,structure) resul if (len_trim(structure) /= 3) & call IO_error(137,ext_msg='lattice_interaction_SlipBySlip: '//trim(structure)) - select case(structure(1:3)) + select case(structure) case('fcc') interactionTypes = FCC_INTERACTIONSLIPSLIP - NslipMax = LATTICE_FCC_NSLIPSYSTEM + NslipMax = FCC_NSLIPSYSTEM case('bcc') interactionTypes = BCC_INTERACTIONSLIPSLIP - NslipMax = LATTICE_BCC_NSLIPSYSTEM + NslipMax = BCC_NSLIPSYSTEM case('hex') interactionTypes = HEX_INTERACTIONSLIPSLIP - NslipMax = LATTICE_HEX_NSLIPSYSTEM + NslipMax = HEX_NSLIPSYSTEM case('bct') interactionTypes = BCT_INTERACTIONSLIPSLIP - NslipMax = LATTICE_BCT_NSLIPSYSTEM + NslipMax = BCT_NSLIPSYSTEM case default call IO_error(137,ext_msg='lattice_interaction_SlipBySlip: '//trim(structure)) end select @@ -1252,7 +1000,7 @@ function lattice_interaction_TwinByTwin(Ntwin,interactionValues,structure) resul integer, dimension(:), allocatable :: NtwinMax integer, dimension(:,:), allocatable :: interactionTypes - integer, dimension(LATTICE_FCC_NTWIN,LATTICE_FCC_NTWIN), parameter :: & + integer, dimension(FCC_NTWIN,FCC_NTWIN), parameter :: & FCC_INTERACTIONTWINTWIN = reshape( [& 1,1,1,2,2,2,2,2,2,2,2,2, & ! -----> acting 1,1,1,2,2,2,2,2,2,2,2,2, & ! | @@ -1268,7 +1016,7 @@ function lattice_interaction_TwinByTwin(Ntwin,interactionValues,structure) resul 2,2,2,2,2,2,2,2,2,1,1,1 & ],shape(FCC_INTERACTIONTWINTWIN)) !< Twin-twin interaction types for fcc - integer, dimension(LATTICE_BCC_NTWIN,LATTICE_BCC_NTWIN), parameter :: & + integer, dimension(BCC_NTWIN,BCC_NTWIN), parameter :: & BCC_INTERACTIONTWINTWIN = reshape( [& 1,3,3,3,3,3,3,2,3,3,2,3, & ! -----> acting 3,1,3,3,3,3,2,3,3,3,3,2, & ! | @@ -1286,7 +1034,7 @@ function lattice_interaction_TwinByTwin(Ntwin,interactionValues,structure) resul !< 1: self interaction !< 2: collinear interaction !< 3: other interaction - integer, dimension(LATTICE_HEX_NTWIN,LATTICE_HEX_NTWIN), parameter :: & + integer, dimension(HEX_NTWIN,HEX_NTWIN), parameter :: & HEX_INTERACTIONTWINTWIN = reshape( [& 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 7, 7, 7, 7, 7, 7, 13,13,13,13,13,13, & ! -----> acting 2, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 7, 7, 7, 7, 7, 7, 13,13,13,13,13,13, & ! | @@ -1320,16 +1068,16 @@ function lattice_interaction_TwinByTwin(Ntwin,interactionValues,structure) resul if (len_trim(structure) /= 3) & call IO_error(137,ext_msg='lattice_interaction_TwinByTwin: '//trim(structure)) - select case(structure(1:3)) + select case(structure) case('fcc') interactionTypes = FCC_INTERACTIONTWINTWIN - NtwinMax = LATTICE_FCC_NTWINSYSTEM + NtwinMax = FCC_NTWINSYSTEM case('bcc') interactionTypes = BCC_INTERACTIONTWINTWIN - NtwinMax = LATTICE_BCC_NTWINSYSTEM + NtwinMax = BCC_NTWINSYSTEM case('hex') interactionTypes = HEX_INTERACTIONTWINTWIN - NtwinMax = LATTICE_HEX_NTWINSYSTEM + NtwinMax = HEX_NTWINSYSTEM case default call IO_error(137,ext_msg='lattice_interaction_TwinByTwin: '//trim(structure)) end select @@ -1353,7 +1101,7 @@ function lattice_interaction_TransByTrans(Ntrans,interactionValues,structure) re integer, dimension(:), allocatable :: NtransMax integer, dimension(:,:), allocatable :: interactionTypes - integer, dimension(LATTICE_FCC_NTRANS,LATTICE_FCC_NTRANS), parameter :: & + integer, dimension(FCC_NTRANS,FCC_NTRANS), parameter :: & FCC_INTERACTIONTRANSTRANS = reshape( [& 1,1,1,2,2,2,2,2,2,2,2,2, & ! -----> acting 1,1,1,2,2,2,2,2,2,2,2,2, & ! | @@ -1372,9 +1120,9 @@ function lattice_interaction_TransByTrans(Ntrans,interactionValues,structure) re if (len_trim(structure) /= 3) & call IO_error(137,ext_msg='lattice_interaction_TransByTrans: '//trim(structure)) - if(structure(1:3) == 'fcc') then + if(structure == 'fcc') then interactionTypes = FCC_INTERACTIONTRANSTRANS - NtransMax = LATTICE_FCC_NTRANSSYSTEM + NtransMax = FCC_NTRANSSYSTEM else call IO_error(137,ext_msg='lattice_interaction_TransByTrans: '//trim(structure)) end if @@ -1400,7 +1148,7 @@ function lattice_interaction_SlipByTwin(Nslip,Ntwin,interactionValues,structure) NtwinMax integer, dimension(:,:), allocatable :: interactionTypes - integer, dimension(LATTICE_FCC_NTWIN,LATTICE_FCC_NSLIP), parameter :: & + integer, dimension(FCC_NTWIN,FCC_NSLIP), parameter :: & FCC_INTERACTIONSLIPTWIN = reshape( [& 1,1,1,3,3,3,2,2,2,3,3,3, & ! -----> twin (acting) 1,1,1,3,3,3,3,3,3,2,2,2, & ! | @@ -1425,7 +1173,7 @@ function lattice_interaction_SlipByTwin(Nslip,Ntwin,interactionValues,structure) !< 1: coplanar interaction !< 2: screw trace between slip system and twin habit plane (easy cross slip) !< 3: other interaction - integer, dimension(LATTICE_BCC_NTWIN,LATTICE_BCC_NSLIP), parameter :: & + integer, dimension(BCC_NTWIN,BCC_NSLIP), parameter :: & BCC_INTERACTIONSLIPTWIN = reshape( [& 3,3,3,2,2,3,3,3,3,2,3,3, & ! -----> twin (acting) 3,3,2,3,3,2,3,3,2,3,3,3, & ! | @@ -1456,7 +1204,7 @@ function lattice_interaction_SlipByTwin(Nslip,Ntwin,interactionValues,structure) !< 1: coplanar interaction !< 2: screw trace between slip system and twin habit plane (easy cross slip) !< 3: other interaction - integer, dimension(LATTICE_HEX_NTWIN,LATTICE_HEX_NSLIP), parameter :: & + integer, dimension(HEX_NTWIN,HEX_NSLIP), parameter :: & HEX_INTERACTIONSLIPTWIN = reshape( [& 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, & ! ----> twin (acting) 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, & ! | @@ -1502,19 +1250,19 @@ function lattice_interaction_SlipByTwin(Nslip,Ntwin,interactionValues,structure) if (len_trim(structure) /= 3) & call IO_error(137,ext_msg='lattice_interaction_SlipByTwin: '//trim(structure)) - select case(structure(1:3)) + select case(structure) case('fcc') interactionTypes = FCC_INTERACTIONSLIPTWIN - NslipMax = LATTICE_FCC_NSLIPSYSTEM - NtwinMax = LATTICE_FCC_NTWINSYSTEM + NslipMax = FCC_NSLIPSYSTEM + NtwinMax = FCC_NTWINSYSTEM case('bcc') interactionTypes = BCC_INTERACTIONSLIPTWIN - NslipMax = LATTICE_BCC_NSLIPSYSTEM - NtwinMax = LATTICE_BCC_NTWINSYSTEM + NslipMax = BCC_NSLIPSYSTEM + NtwinMax = BCC_NTWINSYSTEM case('hex') interactionTypes = HEX_INTERACTIONSLIPTWIN - NslipMax = LATTICE_HEX_NSLIPSYSTEM - NtwinMax = LATTICE_HEX_NTWINSYSTEM + NslipMax = HEX_NSLIPSYSTEM + NtwinMax = HEX_NTWINSYSTEM case default call IO_error(137,ext_msg='lattice_interaction_SlipByTwin: '//trim(structure)) end select @@ -1540,7 +1288,7 @@ function lattice_interaction_SlipByTrans(Nslip,Ntrans,interactionValues,structur NtransMax integer, dimension(:,:), allocatable :: interactionTypes - integer, dimension(LATTICE_FCC_NTRANS,LATTICE_FCC_NSLIP), parameter :: & + integer, dimension(FCC_NTRANS,FCC_NSLIP), parameter :: & FCC_INTERACTIONSLIPTRANS = reshape( [& 1,1,1,3,3,3,2,2,2,3,3,3, & ! -----> trans (acting) 1,1,1,3,3,3,3,3,3,2,2,2, & ! | @@ -1566,11 +1314,11 @@ function lattice_interaction_SlipByTrans(Nslip,Ntrans,interactionValues,structur if (len_trim(structure) /= 3) & call IO_error(137,ext_msg='lattice_interaction_SlipByTrans: '//trim(structure)) - select case(structure(1:3)) + select case(structure) case('fcc') interactionTypes = FCC_INTERACTIONSLIPTRANS - NslipMax = LATTICE_FCC_NSLIPSYSTEM - NtransMax = LATTICE_FCC_NTRANSSYSTEM + NslipMax = FCC_NSLIPSYSTEM + NtransMax = FCC_NTRANSSYSTEM case default call IO_error(137,ext_msg='lattice_interaction_SlipByTrans: '//trim(structure)) end select @@ -1596,13 +1344,13 @@ function lattice_interaction_TwinBySlip(Ntwin,Nslip,interactionValues,structure) NslipMax integer, dimension(:,:), allocatable :: interactionTypes - integer, dimension(LATTICE_FCC_NSLIP,LATTICE_FCC_NTWIN), parameter :: & + integer, dimension(FCC_NSLIP,FCC_NTWIN), parameter :: & FCC_INTERACTIONTWINSLIP = 1 !< Twin-slip interaction types for fcc - integer, dimension(LATTICE_BCC_NSLIP,LATTICE_BCC_NTWIN), parameter :: & + integer, dimension(BCC_NSLIP,BCC_NTWIN), parameter :: & BCC_INTERACTIONTWINSLIP = 1 !< Twin-slip interaction types for bcc - integer, dimension(LATTICE_HEX_NSLIP,LATTICE_HEX_NTWIN), parameter :: & + integer, dimension(HEX_NSLIP,HEX_NTWIN), parameter :: & HEX_INTERACTIONTWINSLIP = reshape( [& 1, 1, 1, 5, 5, 5, 9, 9, 9, 13,13,13,13,13,13, 17,17,17,17,17,17,17,17,17,17,17,17, 21,21,21,21,21,21, & ! ----> slip (acting) 1, 1, 1, 5, 5, 5, 9, 9, 9, 13,13,13,13,13,13, 17,17,17,17,17,17,17,17,17,17,17,17, 21,21,21,21,21,21, & ! | @@ -1636,19 +1384,19 @@ function lattice_interaction_TwinBySlip(Ntwin,Nslip,interactionValues,structure) if (len_trim(structure) /= 3) & call IO_error(137,ext_msg='lattice_interaction_TwinBySlip: '//trim(structure)) - select case(structure(1:3)) + select case(structure) case('fcc') interactionTypes = FCC_INTERACTIONTWINSLIP - NtwinMax = LATTICE_FCC_NTWINSYSTEM - NslipMax = LATTICE_FCC_NSLIPSYSTEM + NtwinMax = FCC_NTWINSYSTEM + NslipMax = FCC_NSLIPSYSTEM case('bcc') interactionTypes = BCC_INTERACTIONTWINSLIP - NtwinMax = LATTICE_BCC_NTWINSYSTEM - NslipMax = LATTICE_BCC_NSLIPSYSTEM + NtwinMax = BCC_NTWINSYSTEM + NslipMax = BCC_NSLIPSYSTEM case('hex') interactionTypes = HEX_INTERACTIONTWINSLIP - NtwinMax = LATTICE_HEX_NTWINSYSTEM - NslipMax = LATTICE_HEX_NSLIPSYSTEM + NtwinMax = HEX_NTWINSYSTEM + NslipMax = HEX_NSLIPSYSTEM case default call IO_error(137,ext_msg='lattice_interaction_TwinBySlip: '//trim(structure)) end select @@ -1677,19 +1425,19 @@ function lattice_SchmidMatrix_slip(Nslip,structure,cOverA) result(SchmidMatrix) if (len_trim(structure) /= 3) & call IO_error(137,ext_msg='lattice_SchmidMatrix_slip: '//trim(structure)) - select case(structure(1:3)) + select case(structure) case('fcc') - NslipMax = LATTICE_FCC_NSLIPSYSTEM - slipSystems = LATTICE_FCC_SYSTEMSLIP + NslipMax = FCC_NSLIPSYSTEM + slipSystems = FCC_SYSTEMSLIP case('bcc') - NslipMax = LATTICE_BCC_NSLIPSYSTEM - slipSystems = LATTICE_BCC_SYSTEMSLIP + NslipMax = BCC_NSLIPSYSTEM + slipSystems = BCC_SYSTEMSLIP case('hex') - NslipMax = LATTICE_HEX_NSLIPSYSTEM - slipSystems = LATTICE_HEX_SYSTEMSLIP + NslipMax = HEX_NSLIPSYSTEM + slipSystems = HEX_SYSTEMSLIP case('bct') - NslipMax = LATTICE_BCT_NSLIPSYSTEM - slipSystems = LATTICE_BCT_SYSTEMSLIP + NslipMax = BCT_NSLIPSYSTEM + slipSystems = BCT_SYSTEMSLIP case default call IO_error(137,ext_msg='lattice_SchmidMatrix_slip: '//trim(structure)) end select @@ -1729,16 +1477,16 @@ function lattice_SchmidMatrix_twin(Ntwin,structure,cOverA) result(SchmidMatrix) if (len_trim(structure) /= 3) & call IO_error(137,ext_msg='lattice_SchmidMatrix_twin: '//trim(structure)) - select case(structure(1:3)) + select case(structure) case('fcc') - NtwinMax = LATTICE_FCC_NTWINSYSTEM - twinSystems = LATTICE_FCC_SYSTEMTWIN + NtwinMax = FCC_NTWINSYSTEM + twinSystems = FCC_SYSTEMTWIN case('bcc') - NtwinMax = LATTICE_BCC_NTWINSYSTEM - twinSystems = LATTICE_BCC_SYSTEMTWIN + NtwinMax = BCC_NTWINSYSTEM + twinSystems = BCC_SYSTEMTWIN case('hex') - NtwinMax = LATTICE_HEX_NTWINSYSTEM - twinSystems = LATTICE_HEX_SYSTEMTWIN + NtwinMax = HEX_NTWINSYSTEM + twinSystems = HEX_SYSTEMTWIN case default call IO_error(137,ext_msg='lattice_SchmidMatrix_twin: '//trim(structure)) end select @@ -1808,22 +1556,16 @@ function lattice_SchmidMatrix_cleavage(Ncleavage,structure,cOverA) result(Schmid if (len_trim(structure) /= 3) & call IO_error(137,ext_msg='lattice_SchmidMatrix_cleavage: '//trim(structure)) - select case(structure(1:3)) - case('iso') - NcleavageMax = LATTICE_ISO_NCLEAVAGESYSTEM - cleavageSystems = LATTICE_ISO_SYSTEMCLEAVAGE + select case(structure) case('ort') - NcleavageMax = LATTICE_ORT_NCLEAVAGESYSTEM - cleavageSystems = LATTICE_ORT_SYSTEMCLEAVAGE + NcleavageMax = ORT_NCLEAVAGESYSTEM + cleavageSystems = ORT_SYSTEMCLEAVAGE case('fcc') - NcleavageMax = LATTICE_FCC_NCLEAVAGESYSTEM - cleavageSystems = LATTICE_FCC_SYSTEMCLEAVAGE + NcleavageMax = FCC_NCLEAVAGESYSTEM + cleavageSystems = FCC_SYSTEMCLEAVAGE case('bcc') - NcleavageMax = LATTICE_BCC_NCLEAVAGESYSTEM - cleavageSystems = LATTICE_BCC_SYSTEMCLEAVAGE - case('hex') - NcleavageMax = LATTICE_HEX_NCLEAVAGESYSTEM - cleavageSystems = LATTICE_HEX_SYSTEMCLEAVAGE + NcleavageMax = BCC_NCLEAVAGESYSTEM + cleavageSystems = BCC_SYSTEMCLEAVAGE case default call IO_error(137,ext_msg='lattice_SchmidMatrix_cleavage: '//trim(structure)) end select @@ -1898,6 +1640,209 @@ function lattice_slip_transverse(Nslip,structure,cOverA) result(t) end function lattice_slip_transverse +!-------------------------------------------------------------------------------------------------- +!> @brief Labels for slip systems +!> details only active slip systems are considered +!-------------------------------------------------------------------------------------------------- +function lattice_labels_slip(Nslip,structure) result(labels) + + integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family + character(len=*), intent(in) :: structure !< lattice structure + + character(len=:), dimension(:), allocatable :: labels + + real(pReal), dimension(:,:), allocatable :: slipSystems + integer, dimension(:), allocatable :: NslipMax + + if (len_trim(structure) /= 3) & + call IO_error(137,ext_msg='lattice_labels_slip: '//trim(structure)) + + select case(structure) + case('fcc') + NslipMax = FCC_NSLIPSYSTEM + slipSystems = FCC_SYSTEMSLIP + case('bcc') + NslipMax = BCC_NSLIPSYSTEM + slipSystems = BCC_SYSTEMSLIP + case('hex') + NslipMax = HEX_NSLIPSYSTEM + slipSystems = HEX_SYSTEMSLIP + case('bct') + NslipMax = BCT_NSLIPSYSTEM + slipSystems = BCT_SYSTEMSLIP + case default + call IO_error(137,ext_msg='lattice_labels_slip: '//trim(structure)) + end select + + if (any(NslipMax(1:size(Nslip)) - Nslip < 0)) & + call IO_error(145,ext_msg='Nslip '//trim(structure)) + if (any(Nslip < 0)) & + call IO_error(144,ext_msg='Nslip '//trim(structure)) + + labels = getLabels(Nslip,NslipMax,slipSystems) + +end function lattice_labels_slip + + +!-------------------------------------------------------------------------------------------------- +!> @brief Return 3x3 tensor with symmetry according to given crystal structure +!-------------------------------------------------------------------------------------------------- +function lattice_applyLatticeSymmetry33(T,structure) result(T_sym) + + real(pReal), dimension(3,3) :: T_sym + + real(pReal), dimension(3,3), intent(in) :: T + character(len=*), intent(in) :: structure + + integer :: k + + T_sym = 0.0_pReal + + if (len_trim(structure) /= 3) & + call IO_error(137,ext_msg='lattice_applyLatticeSymmetry33: '//trim(structure)) + + select case(structure) + case('iso','fcc','bcc') + do k=1,3 + T_sym(k,k) = T(1,1) + enddo + case('hex') + T_sym(1,1) = T(1,1) + T_sym(2,2) = T(1,1) + T_sym(3,3) = T(3,3) + case('ort','bct') + T_sym(1,1) = T(1,1) + T_sym(2,2) = T(2,2) + T_sym(3,3) = T(3,3) + case default + call IO_error(137,ext_msg='lattice_applyLatticeSymmetry33: '//trim(structure)) + end select + +end function lattice_applyLatticeSymmetry33 + + +!-------------------------------------------------------------------------------------------------- +!> @brief Return stiffness matrix in 6x6 notation with symmetry according to given crystal structure +!> @details J. A. Rayne and B. S. Chandrasekhar Phys. Rev. 120, 1658 Erratum Phys. Rev. 122, 1962 +!-------------------------------------------------------------------------------------------------- +function applyLatticeSymmetryC66(C66,structure) result(C66_sym) + + real(pReal), dimension(6,6) :: C66_sym + + real(pReal), dimension(6,6), intent(in) :: C66 + character(len=*), intent(in) :: structure + + integer :: j,k + + C66_sym = 0.0_pReal + + if (len_trim(structure) /= 3) & + call IO_error(137,ext_msg='applyLatticeSymmetryC66: '//trim(structure)) + + select case(structure) + case ('iso') + do k=1,3 + do j=1,3 + C66_sym(k,j) = C66(1,2) + enddo + C66_sym(k,k) = C66(1,1) + C66_sym(k+3,k+3) = 0.5_pReal*(C66(1,1)-C66(1,2)) + enddo + case ('fcc','bcc') + do k=1,3 + do j=1,3 + C66_sym(k,j) = C66(1,2) + enddo + C66_sym(k,k) = C66(1,1) + C66_sym(k+3,k+3) = C66(4,4) + enddo + case ('hex') + C66_sym(1,1) = C66(1,1) + C66_sym(2,2) = C66(1,1) + C66_sym(3,3) = C66(3,3) + C66_sym(1,2) = C66(1,2) + C66_sym(2,1) = C66(1,2) + C66_sym(1,3) = C66(1,3) + C66_sym(3,1) = C66(1,3) + C66_sym(2,3) = C66(1,3) + C66_sym(3,2) = C66(1,3) + C66_sym(4,4) = C66(4,4) + C66_sym(5,5) = C66(4,4) + C66_sym(6,6) = 0.5_pReal*(C66(1,1)-C66(1,2)) + case ('ort') + C66_sym(1,1) = C66(1,1) + C66_sym(2,2) = C66(2,2) + C66_sym(3,3) = C66(3,3) + C66_sym(1,2) = C66(1,2) + C66_sym(2,1) = C66(1,2) + C66_sym(1,3) = C66(1,3) + C66_sym(3,1) = C66(1,3) + C66_sym(2,3) = C66(2,3) + C66_sym(3,2) = C66(2,3) + C66_sym(4,4) = C66(4,4) + C66_sym(5,5) = C66(5,5) + C66_sym(6,6) = C66(6,6) + case ('bct') + C66_sym(1,1) = C66(1,1) + C66_sym(2,2) = C66(1,1) + C66_sym(3,3) = C66(3,3) + C66_sym(1,2) = C66(1,2) + C66_sym(2,1) = C66(1,2) + C66_sym(1,3) = C66(1,3) + C66_sym(3,1) = C66(1,3) + C66_sym(2,3) = C66(1,3) + C66_sym(3,2) = C66(1,3) + C66_sym(4,4) = C66(4,4) + C66_sym(5,5) = C66(4,4) + C66_sym(6,6) = C66(6,6) + case default + call IO_error(137,ext_msg='applyLatticeSymmetryC66: '//trim(structure)) + end select + +end function applyLatticeSymmetryC66 + + +!-------------------------------------------------------------------------------------------------- +!> @brief Labels for twin systems +!> details only active twin systems are considered +!-------------------------------------------------------------------------------------------------- +function lattice_labels_twin(Ntwin,structure) result(labels) + + integer, dimension(:), intent(in) :: Ntwin !< number of active slip systems per family + character(len=*), intent(in) :: structure !< lattice structure + + character(len=:), dimension(:), allocatable :: labels + + real(pReal), dimension(:,:), allocatable :: twinSystems + integer, dimension(:), allocatable :: NtwinMax + + if (len_trim(structure) /= 3) & + call IO_error(137,ext_msg='lattice_labels_twin: '//trim(structure)) + + select case(structure) + case('fcc') + NtwinMax = FCC_NTWINSYSTEM + twinSystems = FCC_SYSTEMTWIN + case('bcc') + NtwinMax = BCC_NTWINSYSTEM + twinSystems = BCC_SYSTEMTWIN + case('hex') + NtwinMax = HEX_NTWINSYSTEM + twinSystems = HEX_SYSTEMTWIN + case default + call IO_error(137,ext_msg='lattice_labels_twin: '//trim(structure)) + end select + + if (any(NtwinMax(1:size(Ntwin)) - Ntwin < 0)) & + call IO_error(145,ext_msg='Ntwin '//trim(structure)) + if (any(Ntwin < 0)) & + call IO_error(144,ext_msg='Ntwin '//trim(structure)) + + labels = getLabels(Ntwin,NtwinMax,twinSystems) + +end function lattice_labels_twin + + !-------------------------------------------------------------------------------------------------- !> @brief Projection of the transverse direction onto the slip plane !> @details: This projection is used to calculate forest hardening for edge dislocations @@ -1922,91 +1867,6 @@ function slipProjection_transverse(Nslip,structure,cOverA) result(projection) end function slipProjection_transverse -!-------------------------------------------------------------------------------------------------- -!> @brief Labels for slip systems -!> details only active slip systems are considered -!-------------------------------------------------------------------------------------------------- -function lattice_labels_slip(Nslip,structure) result(labels) - - integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family - character(len=*), intent(in) :: structure !< lattice structure - - character(len=:), dimension(:), allocatable :: labels - - real(pReal), dimension(:,:), allocatable :: slipSystems - integer, dimension(:), allocatable :: NslipMax - - if (len_trim(structure) /= 3) & - call IO_error(137,ext_msg='lattice_labels_slip: '//trim(structure)) - - select case(structure(1:3)) - case('fcc') - NslipMax = LATTICE_FCC_NSLIPSYSTEM - slipSystems = LATTICE_FCC_SYSTEMSLIP - case('bcc') - NslipMax = LATTICE_BCC_NSLIPSYSTEM - slipSystems = LATTICE_BCC_SYSTEMSLIP - case('hex') - NslipMax = LATTICE_HEX_NSLIPSYSTEM - slipSystems = LATTICE_HEX_SYSTEMSLIP - case('bct') - NslipMax = LATTICE_BCT_NSLIPSYSTEM - slipSystems = LATTICE_BCT_SYSTEMSLIP - case default - call IO_error(137,ext_msg='lattice_labels_slip: '//trim(structure)) - end select - - if (any(NslipMax(1:size(Nslip)) - Nslip < 0)) & - call IO_error(145,ext_msg='Nslip '//trim(structure)) - if (any(Nslip < 0)) & - call IO_error(144,ext_msg='Nslip '//trim(structure)) - - labels = getLabels(Nslip,NslipMax,slipSystems) - -end function lattice_labels_slip - - -!-------------------------------------------------------------------------------------------------- -!> @brief Labels for twin systems -!> details only active twin systems are considered -!-------------------------------------------------------------------------------------------------- -function lattice_labels_twin(Ntwin,structure) result(labels) - - integer, dimension(:), intent(in) :: Ntwin !< number of active slip systems per family - character(len=*), intent(in) :: structure !< lattice structure - - character(len=:), dimension(:), allocatable :: labels - - real(pReal), dimension(:,:), allocatable :: twinSystems - integer, dimension(:), allocatable :: NtwinMax - - if (len_trim(structure) /= 3) & - call IO_error(137,ext_msg='lattice_labels_twin: '//trim(structure)) - - select case(structure(1:3)) - case('fcc') - NtwinMax = LATTICE_FCC_NTWINSYSTEM - twinSystems = LATTICE_FCC_SYSTEMTWIN - case('bcc') - NtwinMax = LATTICE_BCC_NTWINSYSTEM - twinSystems = LATTICE_BCC_SYSTEMTWIN - case('hex') - NtwinMax = LATTICE_HEX_NTWINSYSTEM - twinSystems = LATTICE_HEX_SYSTEMTWIN - case default - call IO_error(137,ext_msg='lattice_labels_twin: '//trim(structure)) - end select - - if (any(NtwinMax(1:size(Ntwin)) - Ntwin < 0)) & - call IO_error(145,ext_msg='Ntwin '//trim(structure)) - if (any(Ntwin < 0)) & - call IO_error(144,ext_msg='Ntwin '//trim(structure)) - - labels = getLabels(Ntwin,NtwinMax,twinSystems) - -end function lattice_labels_twin - - !-------------------------------------------------------------------------------------------------- !> @brief Projection of the slip direction onto the slip plane !> @details: This projection is used to calculate forest hardening for screw dislocations @@ -2048,19 +1908,19 @@ function coordinateSystem_slip(Nslip,structure,cOverA) result(coordinateSystem) if (len_trim(structure) /= 3) & call IO_error(137,ext_msg='coordinateSystem_slip: '//trim(structure)) - select case(structure(1:3)) + select case(structure) case('fcc') - NslipMax = LATTICE_FCC_NSLIPSYSTEM - slipSystems = LATTICE_FCC_SYSTEMSLIP + NslipMax = FCC_NSLIPSYSTEM + slipSystems = FCC_SYSTEMSLIP case('bcc') - NslipMax = LATTICE_BCC_NSLIPSYSTEM - slipSystems = LATTICE_BCC_SYSTEMSLIP + NslipMax = BCC_NSLIPSYSTEM + slipSystems = BCC_SYSTEMSLIP case('hex') - NslipMax = LATTICE_HEX_NSLIPSYSTEM - slipSystems = LATTICE_HEX_SYSTEMSLIP + NslipMax = HEX_NSLIPSYSTEM + slipSystems = HEX_SYSTEMSLIP case('bct') - NslipMax = LATTICE_BCT_NSLIPSYSTEM - slipSystems = LATTICE_BCT_SYSTEMSLIP + NslipMax = BCT_NSLIPSYSTEM + slipSystems = BCT_SYSTEMSLIP case default call IO_error(137,ext_msg='coordinateSystem_slip: '//trim(structure)) end select @@ -2080,40 +1940,40 @@ end function coordinateSystem_slip !-------------------------------------------------------------------------------------------------- function buildInteraction(reacting_used,acting_used,reacting_max,acting_max,values,matrix) - integer, dimension(:), intent(in) :: & - reacting_used, & !< # of reacting systems per family as specified in material.config - acting_used, & !< # of acting systems per family as specified in material.config - reacting_max, & !< max # of reacting systems per family for given lattice - acting_max !< max # of acting systems per family for given lattice - real(pReal), dimension(:), intent(in) :: values !< interaction values - integer, dimension(:,:), intent(in) :: matrix !< interaction types - real(pReal), dimension(sum(reacting_used),sum(acting_used)) :: buildInteraction + integer, dimension(:), intent(in) :: & + reacting_used, & !< # of reacting systems per family as specified in material.config + acting_used, & !< # of acting systems per family as specified in material.config + reacting_max, & !< max # of reacting systems per family for given lattice + acting_max !< max # of acting systems per family for given lattice + real(pReal), dimension(:), intent(in) :: values !< interaction values + integer, dimension(:,:), intent(in) :: matrix !< interaction types + real(pReal), dimension(sum(reacting_used),sum(acting_used)) :: buildInteraction - integer :: & - acting_family_index, acting_family, acting_system, & - reacting_family_index, reacting_family, reacting_system, & - i,j,k,l + integer :: & + acting_family_index, acting_family, acting_system, & + reacting_family_index, reacting_family, reacting_system, & + i,j,k,l - do acting_family = 1,size(acting_used,1) - acting_family_index = sum(acting_used(1:acting_family-1)) - do acting_system = 1,acting_used(acting_family) + do acting_family = 1,size(acting_used,1) + acting_family_index = sum(acting_used(1:acting_family-1)) + do acting_system = 1,acting_used(acting_family) - do reacting_family = 1,size(reacting_used,1) - reacting_family_index = sum(reacting_used(1:reacting_family-1)) - do reacting_system = 1,reacting_used(reacting_family) + do reacting_family = 1,size(reacting_used,1) + reacting_family_index = sum(reacting_used(1:reacting_family-1)) + do reacting_system = 1,reacting_used(reacting_family) - i = sum( acting_max(1: acting_family-1)) + acting_system - j = sum(reacting_max(1:reacting_family-1)) + reacting_system + i = sum( acting_max(1: acting_family-1)) + acting_system + j = sum(reacting_max(1:reacting_family-1)) + reacting_system - k = acting_family_index + acting_system - l = reacting_family_index + reacting_system + k = acting_family_index + acting_system + l = reacting_family_index + reacting_system - if (matrix(i,j) > size(values)) call IO_error(138,ext_msg='buildInteraction') + if (matrix(i,j) > size(values)) call IO_error(138,ext_msg='buildInteraction') - buildInteraction(l,k) = values(matrix(i,j)) + buildInteraction(l,k) = values(matrix(i,j)) - enddo; enddo - enddo; enddo + enddo; enddo + enddo; enddo end function buildInteraction @@ -2146,9 +2006,9 @@ function buildCoordinateSystem(active,potential,system,structure,cOverA) if (len_trim(structure) /= 3) & call IO_error(137,ext_msg='buildCoordinateSystem: '//trim(structure)) - if (trim(structure(1:3)) == 'bct' .and. cOverA > 2.0_pReal) & + if (trim(structure) == 'bct' .and. cOverA > 2.0_pReal) & call IO_error(131,ext_msg='buildCoordinateSystem:'//trim(structure)) - if (trim(structure(1:3)) == 'hex' .and. (cOverA < 1.0_pReal .or. cOverA > 2.0_pReal)) & + if (trim(structure) == 'hex' .and. (cOverA < 1.0_pReal .or. cOverA > 2.0_pReal)) & call IO_error(131,ext_msg='buildCoordinateSystem:'//trim(structure)) a = 0 @@ -2157,7 +2017,7 @@ function buildCoordinateSystem(active,potential,system,structure,cOverA) a = a + 1 p = sum(potential(1:f-1))+s - select case(trim(structure(1:3))) + select case(trim(structure)) case ('fcc','bcc','iso','ort','bct') direction = system(1:3,p) @@ -2215,8 +2075,8 @@ subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_fcc,a_bcc) x, y, z integer :: & i - real(pReal), dimension(3+3,LATTICE_FCC_NTRANS), parameter :: & - LATTICE_FCCTOHEX_SYSTEMTRANS = reshape(real( [& + real(pReal), dimension(3+3,FCC_NTRANS), parameter :: & + FCCTOHEX_SYSTEMTRANS = reshape(real( [& -2, 1, 1, 1, 1, 1, & 1,-2, 1, 1, 1, 1, & 1, 1,-2, 1, 1, 1, & @@ -2229,9 +2089,9 @@ subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_fcc,a_bcc) 2, 1,-1, -1, 1,-1, & -1,-2,-1, -1, 1,-1, & -1, 1, 2, -1, 1,-1 & - ],pReal),shape(LATTICE_FCCTOHEX_SYSTEMTRANS)) - real(pReal), dimension(4,LATTICE_fcc_Ntrans), parameter :: & - LATTICE_FCCTOBCC_SYSTEMTRANS = reshape([& + ],pReal),shape(FCCTOHEX_SYSTEMTRANS)) + real(pReal), dimension(4,fcc_Ntrans), parameter :: & + FCCTOBCC_SYSTEMTRANS = reshape([& 0.0, 1.0, 0.0, 10.26, & ! Pitsch OR (Ma & Hartmaier 2014, Table 3) 0.0,-1.0, 0.0, 10.26, & 0.0, 0.0, 1.0, 10.26, & @@ -2244,10 +2104,10 @@ subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_fcc,a_bcc) -1.0, 0.0, 0.0, 10.26, & 0.0, 1.0, 0.0, 10.26, & 0.0,-1.0, 0.0, 10.26 & - ],shape(LATTICE_FCCTOBCC_SYSTEMTRANS)) + ],shape(FCCTOBCC_SYSTEMTRANS)) - integer, dimension(9,LATTICE_fcc_Ntrans), parameter :: & - LATTICE_FCCTOBCC_BAINVARIANT = reshape( [& + integer, dimension(9,fcc_Ntrans), parameter :: & + FCCTOBCC_BAINVARIANT = reshape( [& 1, 0, 0, 0, 1, 0, 0, 0, 1, & ! Pitsch OR (Ma & Hartmaier 2014, Table 3) 1, 0, 0, 0, 1, 0, 0, 0, 1, & 1, 0, 0, 0, 1, 0, 0, 0, 1, & @@ -2260,10 +2120,10 @@ subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_fcc,a_bcc) 0, 0, 1, 1, 0, 0, 0, 1, 0, & 0, 0, 1, 1, 0, 0, 0, 1, 0, & 0, 0, 1, 1, 0, 0, 0, 1, 0 & - ],shape(LATTICE_FCCTOBCC_BAINVARIANT)) + ],shape(FCCTOBCC_BAINVARIANT)) - real(pReal), dimension(4,LATTICE_fcc_Ntrans), parameter :: & - LATTICE_FCCTOBCC_BAINROT = reshape([& + real(pReal), dimension(4,fcc_Ntrans), parameter :: & + FCCTOBCC_BAINROT = reshape([& 1.0, 0.0, 0.0, 45.0, & ! Rotate fcc austensite to bain variant 1.0, 0.0, 0.0, 45.0, & 1.0, 0.0, 0.0, 45.0, & @@ -2276,15 +2136,15 @@ subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_fcc,a_bcc) 0.0, 0.0, 1.0, 45.0, & 0.0, 0.0, 1.0, 45.0, & 0.0, 0.0, 1.0, 45.0 & - ],shape(LATTICE_FCCTOBCC_BAINROT)) + ],shape(FCCTOBCC_BAINROT)) if (a_bcc > 0.0_pReal .and. a_fcc > 0.0_pReal .and. dEq0(cOverA)) then ! fcc -> bcc transformation do i = 1,sum(Ntrans) - call R%fromAxisAngle(LATTICE_FCCTOBCC_SYSTEMTRANS(:,i),degrees=.true.,P=1) - call B%fromAxisAngle(LATTICE_FCCTOBCC_BAINROT(:,i), degrees=.true.,P=1) - x = real(LATTICE_FCCTOBCC_BAINVARIANT(1:3,i),pReal) - y = real(LATTICE_FCCTOBCC_BAINVARIANT(4:6,i),pReal) - z = real(LATTICE_FCCTOBCC_BAINVARIANT(7:9,i),pReal) + call R%fromAxisAngle(FCCTOBCC_SYSTEMTRANS(:,i),degrees=.true.,P=1) + call B%fromAxisAngle(FCCTOBCC_BAINROT(:,i), degrees=.true.,P=1) + x = real(FCCTOBCC_BAINVARIANT(1:3,i),pReal) + y = real(FCCTOBCC_BAINVARIANT(4:6,i),pReal) + z = real(FCCTOBCC_BAINVARIANT(7:9,i),pReal) U = (a_bcc/a_fcc)*math_outer(x,x) & + (a_bcc/a_fcc)*math_outer(y,y) * sqrt(2.0_pReal) & @@ -2299,8 +2159,8 @@ subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_fcc,a_bcc) sd(3,3) = cOverA/sqrt(8.0_pReal/3.0_pReal) do i = 1,sum(Ntrans) - x = LATTICE_FCCTOHEX_SYSTEMTRANS(1:3,i)/norm2(LATTICE_FCCTOHEX_SYSTEMTRANS(1:3,i)) - z = LATTICE_FCCTOHEX_SYSTEMTRANS(4:6,i)/norm2(LATTICE_FCCTOHEX_SYSTEMTRANS(4:6,i)) + x = FCCTOHEX_SYSTEMTRANS(1:3,i)/norm2(FCCTOHEX_SYSTEMTRANS(1:3,i)) + z = FCCTOHEX_SYSTEMTRANS(4:6,i)/norm2(FCCTOHEX_SYSTEMTRANS(4:6,i)) y = -math_cross(x,z) Q(1:3,1,i) = x Q(1:3,2,i) = y diff --git a/src/prec.f90 b/src/prec.f90 index 8d8e63e09..8ced9813a 100644 --- a/src/prec.f90 +++ b/src/prec.f90 @@ -11,30 +11,31 @@ module prec implicit none public + ! https://software.intel.com/en-us/blogs/2017/03/27/doctor-fortran-in-it-takes-all-kinds - integer, parameter, public :: pReal = IEEE_selected_real_kind(15,307) !< number with 15 significant digits, up to 1e+-307 (typically 64 bit) + integer, parameter :: pReal = IEEE_selected_real_kind(15,307) !< number with 15 significant digits, up to 1e+-307 (typically 64 bit) #if(INT==8) - integer, parameter, public :: pInt = selected_int_kind(18) !< number with at least up to +-1e18 (typically 64 bit) + integer, parameter :: pInt = selected_int_kind(18) !< number with at least up to +-1e18 (typically 64 bit) #else - integer, parameter, public :: pInt = selected_int_kind(9) !< number with at least up to +-1e9 (typically 32 bit) + integer, parameter :: pInt = selected_int_kind(9) !< number with at least up to +-1e9 (typically 32 bit) #endif - integer, parameter, public :: pLongInt = selected_int_kind(18) !< number with at least up to +-1e18 (typically 64 bit) - integer, parameter, public :: pStringLen = 256 !< default string length - integer, parameter, public :: pPathLen = 4096 !< maximum length of a path name on linux + integer, parameter :: pLongInt = selected_int_kind(18) !< number with at least up to +-1e18 (typically 64 bit) + integer, parameter :: pStringLen = 256 !< default string length + integer, parameter :: pPathLen = 4096 !< maximum length of a path name on linux - real(pReal), parameter, public :: tol_math_check = 1.0e-8_pReal !< tolerance for internal math self-checks (rotation) + real(pReal), parameter :: tol_math_check = 1.0e-8_pReal !< tolerance for internal math self-checks (rotation) - type, public :: group_float !< variable length datatype used for storage of state + type :: group_float !< variable length datatype used for storage of state real(pReal), dimension(:), pointer :: p end type group_float - type, public :: group_int + type :: group_int integer, dimension(:), pointer :: p end type group_int ! http://stackoverflow.com/questions/3948210/can-i-have-a-pointer-to-an-item-in-an-allocatable-array - type, public :: tState + type :: tState integer :: & sizeState = 0, & !< size of state sizeDotState = 0, & !< size of dot state, i.e. state(1:sizeDot) follows time evolution by dotState rates @@ -57,7 +58,7 @@ module prec RKCK45dotState end type - type, extends(tState), public :: tPlasticState + type, extends(tState) :: tPlasticState logical :: & nonlocal = .false. real(pReal), pointer, dimension(:,:) :: & @@ -65,22 +66,22 @@ module prec accumulatedSlip !< accumulated plastic slip end type - type, public :: tSourceState + type :: tSourceState type(tState), dimension(:), allocatable :: p !< tState for each active source mechanism in a phase end type - type, public :: tHomogMapping + type :: tHomogMapping integer, pointer, dimension(:,:) :: p end type real(pReal), private, parameter :: PREAL_EPSILON = epsilon(0.0_pReal) !< minimum positive number such that 1.0 + EPSILON /= 1.0. real(pReal), private, parameter :: PREAL_MIN = tiny(0.0_pReal) !< smallest normalized floating point number - integer, dimension(0), parameter, public :: & + integer, dimension(0), parameter :: & emptyIntArray = [integer::] - real(pReal), dimension(0), parameter, public :: & + real(pReal), dimension(0), parameter :: & emptyRealArray = [real(pReal)::] - character(len=pStringLen), dimension(0), parameter, public :: & + character(len=pStringLen), dimension(0), parameter :: & emptyStringArray = [character(len=pStringLen)::] private :: & diff --git a/src/quaternions.f90 b/src/quaternions.f90 index 37131800a..4e40ecc4e 100644 --- a/src/quaternions.f90 +++ b/src/quaternions.f90 @@ -10,7 +10,7 @@ module quaternions use IO implicit none - public + private real(pReal), parameter, public :: P = -1.0_pReal !< parameter for orientation conversion. @@ -95,9 +95,13 @@ module quaternions interface aimag module procedure aimag__ end interface aimag - - private :: & - unitTest + + public :: & + quaternions_init, & + assignment(=), & + conjg, aimag, & + log, exp, & + real contains @@ -511,6 +515,7 @@ subroutine unitTest q_2 = conjg(q_2) - inverse(q_2) if(any(dNeq0(q_2%asArray(),1.0e-15_pReal))) call IO_error(0,ext_msg='inverse/conjg') endif + if(dNeq(dot_product(qu,qu),dot_product(q,q))) call IO_error(0,ext_msg='dot_product') #if !(defined(__GFORTRAN__) && __GNUC__ < 9) if (norm2(aimag(q)) > 0.0_pReal) then diff --git a/src/source_damage_anisoBrittle.f90 b/src/source_damage_anisoBrittle.f90 index e5ed05799..91c90adc2 100644 --- a/src/source_damage_anisoBrittle.f90 +++ b/src/source_damage_anisoBrittle.f90 @@ -9,8 +9,8 @@ module source_damage_anisoBrittle use debug use IO use math - use material use discretization + use material use config use lattice use results @@ -21,17 +21,8 @@ module source_damage_anisoBrittle integer, dimension(:), allocatable :: & source_damage_anisoBrittle_offset, & !< which source is my current source mechanism? source_damage_anisoBrittle_instance !< instance of source mechanism - - integer, dimension(:,:), allocatable :: & - source_damage_anisoBrittle_Ncleavage !< number of cleavage systems per family - enum, bind(c) - enumerator :: undefined_ID, & - damage_drivingforce_ID - end enum - - - type :: tParameters !< container type for internal constitutive parameters + type :: tParameters !< container type for internal constitutive parameters real(pReal) :: & aTol, & sdot_0, & @@ -45,8 +36,8 @@ module source_damage_anisoBrittle totalNcleavage integer, dimension(:), allocatable :: & Ncleavage - integer(kind(undefined_ID)), allocatable, dimension(:) :: & - outputID !< ID of each post result output + character(len=pStringLen), allocatable, dimension(:) :: & + output end type tParameters type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance) @@ -67,101 +58,69 @@ contains !-------------------------------------------------------------------------------------------------- subroutine source_damage_anisoBrittle_init - integer :: Ninstance,phase,instance,source,sourceOffset - integer :: NofMyPhase,p ,i - integer(kind(undefined_ID)) :: & - outputID + integer :: Ninstance,sourceOffset,NofMyPhase,p + character(len=pStringLen) :: extmsg = '' - character(len=pStringLen) :: & - extmsg = '' - character(len=pStringLen), dimension(:), allocatable :: & - outputs + write(6,'(/,a)') ' <<<+- source_'//SOURCE_DAMAGE_ANISOBRITTLE_LABEL//' init -+>>>'; flush(6) - write(6,'(/,a)') ' <<<+- source_'//SOURCE_DAMAGE_ANISOBRITTLE_LABEL//' init -+>>>'; flush(6) - - Ninstance = count(phase_source == SOURCE_damage_anisoBrittle_ID) - if (Ninstance == 0) return - + Ninstance = count(phase_source == SOURCE_DAMAGE_ANISOBRITTLE_ID) if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) & write(6,'(a16,1x,i5,/)') '# instances:',Ninstance - - allocate(source_damage_anisoBrittle_offset(material_Nphase), source=0) - allocate(source_damage_anisoBrittle_instance(material_Nphase), source=0) - do phase = 1, material_Nphase - source_damage_anisoBrittle_instance(phase) = count(phase_source(:,1:phase) == source_damage_anisoBrittle_ID) - do source = 1, phase_Nsources(phase) - if (phase_source(source,phase) == source_damage_anisoBrittle_ID) & - source_damage_anisoBrittle_offset(phase) = source - enddo - enddo - - allocate(source_damage_anisoBrittle_Ncleavage(lattice_maxNcleavageFamily,Ninstance), source=0) + allocate(source_damage_anisoBrittle_offset (size(config_phase)), source=0) + allocate(source_damage_anisoBrittle_instance(size(config_phase)), source=0) allocate(param(Ninstance)) - - do p=1, size(config_phase) + + do p = 1, size(config_phase) + source_damage_anisoBrittle_instance(p) = count(phase_source(:,1:p) == SOURCE_DAMAGE_ANISOBRITTLE_ID) + do sourceOffset = 1, phase_Nsources(p) + if (phase_source(sourceOffset,p) == SOURCE_DAMAGE_ANISOBRITTLE_ID) then + source_damage_anisoBrittle_offset(p) = sourceOffset + exit + endif + enddo + if (all(phase_source(:,p) /= SOURCE_DAMAGE_ANISOBRITTLE_ID)) cycle associate(prm => param(source_damage_anisoBrittle_instance(p)), & config => config_phase(p)) - - prm%aTol = config%getFloat('anisobrittle_atol',defaultVal = 1.0e-3_pReal) + prm%Ncleavage = config%getInts('ncleavage',defaultVal=emptyIntArray) + prm%totalNcleavage = sum(prm%Ncleavage) + + prm%aTol = config%getFloat('anisobrittle_atol',defaultVal = 1.0e-3_pReal) prm%N = config%getFloat('anisobrittle_ratesensitivity') prm%sdot_0 = config%getFloat('anisobrittle_sdot0') - + + prm%critDisp = config%getFloats('anisobrittle_criticaldisplacement',requiredSize=size(prm%Ncleavage)) + prm%critLoad = config%getFloats('anisobrittle_criticalload', requiredSize=size(prm%Ncleavage)) + + prm%cleavage_systems = lattice_SchmidMatrix_cleavage(prm%Ncleavage,config%getString('lattice_structure'),& + config%getFloat('c/a',defaultVal=0.0_pReal)) + + ! expand: family => system + prm%critDisp = math_expand(prm%critDisp, prm%Ncleavage) + prm%critLoad = math_expand(prm%critLoad, prm%Ncleavage) + ! sanity checks - if (prm%aTol < 0.0_pReal) extmsg = trim(extmsg)//' anisobrittle_atol' - - if (prm%N <= 0.0_pReal) extmsg = trim(extmsg)//' anisobrittle_ratesensitivity' - if (prm%sdot_0 <= 0.0_pReal) extmsg = trim(extmsg)//' anisobrittle_sdot0' - - prm%Ncleavage = config%getInts('ncleavage',defaultVal=emptyIntArray) + if (prm%aTol < 0.0_pReal) extmsg = trim(extmsg)//' anisobrittle_atol' + if (prm%N <= 0.0_pReal) extmsg = trim(extmsg)//' anisobrittle_n' + if (prm%sdot_0 <= 0.0_pReal) extmsg = trim(extmsg)//' anisobrittle_sdot0' + if (any(prm%critLoad < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_critLoad' + if (any(prm%critDisp < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_critDisp' - prm%critDisp = config%getFloats('anisobrittle_criticaldisplacement',requiredSize=size(prm%Ncleavage)) - prm%critLoad = config%getFloats('anisobrittle_criticalload', requiredSize=size(prm%Ncleavage)) - - prm%cleavage_systems = lattice_SchmidMatrix_cleavage (prm%Ncleavage,config%getString('lattice_structure'),& - config%getFloat('c/a',defaultVal=0.0_pReal)) - - ! expand: family => system - prm%critDisp = math_expand(prm%critDisp, prm%Ncleavage) - prm%critLoad = math_expand(prm%critLoad, prm%Ncleavage) - - if (any(prm%critLoad < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_criticalload' - if (any(prm%critDisp < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_criticaldisplacement' !-------------------------------------------------------------------------------------------------- ! exit if any parameter is out of range - if (extmsg /= '') & - call IO_error(211,ext_msg=trim(extmsg)//'('//SOURCE_DAMAGE_ANISOBRITTLE_LABEL//')') + if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'('//SOURCE_DAMAGE_ANISOBRITTLE_LABEL//')') !-------------------------------------------------------------------------------------------------- ! output pararameters - outputs = config%getStrings('(output)',defaultVal=emptyStringArray) - allocate(prm%outputID(0)) - do i=1, size(outputs) - outputID = undefined_ID - select case(outputs(i)) - - case ('anisobrittle_drivingforce') - prm%outputID = [prm%outputID, damage_drivingforce_ID] + prm%output = config%getStrings('(output)',defaultVal=emptyStringArray) - end select - - enddo + NofMyPhase = count(material_phaseAt==p) * discretization_nIP + call material_allocateSourceState(p,sourceOffset,NofMyPhase,1,1,0) + sourceState(p)%p(sourceOffset)%aTolState=prm%aTol end associate - - phase = p - NofMyPhase=count(material_phaseAt==phase) * discretization_nIP - instance = source_damage_anisoBrittle_instance(phase) - sourceOffset = source_damage_anisoBrittle_offset(phase) - - - call material_allocateSourceState(phase,sourceOffset,NofMyPhase,1,1,0) - sourceState(phase)%p(sourceOffset)%aTolState=param(instance)%aTol - - - source_damage_anisoBrittle_Ncleavage(1:size(param(instance)%Ncleavage),instance) = param(instance)%Ncleavage enddo end subroutine source_damage_anisoBrittle_init @@ -178,49 +137,42 @@ subroutine source_damage_anisoBrittle_dotState(S, ipc, ip, el) el !< element real(pReal), intent(in), dimension(3,3) :: & S + integer :: & phase, & constituent, & - instance, & sourceOffset, & damageOffset, & homog, & - f, i, index_myFamily, index + i real(pReal) :: & traction_d, traction_t, traction_n, traction_crit phase = material_phaseAt(ipc,el) constituent = material_phasememberAt(ipc,ip,el) - instance = source_damage_anisoBrittle_instance(phase) sourceOffset = source_damage_anisoBrittle_offset(phase) homog = material_homogenizationAt(el) damageOffset = damageMapping(homog)%p(ip,el) - + + associate(prm => param(source_damage_anisoBrittle_instance(phase))) sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = 0.0_pReal - - index = 1 - do f = 1,lattice_maxNcleavageFamily - index_myFamily = sum(lattice_NcleavageSystem(1:f-1,phase)) ! at which index starts my family - do i = 1,source_damage_anisoBrittle_Ncleavage(f,instance) ! process each (active) cleavage system in family + do i = 1, prm%totalNcleavage - traction_d = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,1,index_myFamily+i,phase)) - traction_t = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,2,index_myFamily+i,phase)) - traction_n = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,3,index_myFamily+i,phase)) - - traction_crit = param(instance)%critLoad(index)* & - damage(homog)%p(damageOffset)*damage(homog)%p(damageOffset) + traction_d = math_mul33xx33(S,prm%cleavage_systems(1:3,1:3,1,i)) + traction_t = math_mul33xx33(S,prm%cleavage_systems(1:3,1:3,2,i)) + traction_n = math_mul33xx33(S,prm%cleavage_systems(1:3,1:3,3,i)) - sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = & - sourceState(phase)%p(sourceOffset)%dotState(1,constituent) + & - param(instance)%sdot_0* & - ((max(0.0_pReal, abs(traction_d) - traction_crit)/traction_crit)**param(instance)%N + & - (max(0.0_pReal, abs(traction_t) - traction_crit)/traction_crit)**param(instance)%N + & - (max(0.0_pReal, abs(traction_n) - traction_crit)/traction_crit)**param(instance)%N)/ & - param(instance)%critDisp(index) + traction_crit = prm%critLoad(i)*damage(homog)%p(damageOffset)**2.0_pReal + + sourceState(phase)%p(sourceOffset)%dotState(1,constituent) & + = sourceState(phase)%p(sourceOffset)%dotState(1,constituent) & + + prm%sdot_0 / prm%critDisp(i) & + * ((max(0.0_pReal, abs(traction_d) - traction_crit)/traction_crit)**prm%N + & + (max(0.0_pReal, abs(traction_t) - traction_crit)/traction_crit)**prm%N + & + (max(0.0_pReal, abs(traction_n) - traction_crit)/traction_crit)**prm%N) - index = index + 1 - enddo enddo + end associate end subroutine source_damage_anisoBrittle_dotState @@ -238,16 +190,17 @@ subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalph real(pReal), intent(out) :: & localphiDot, & dLocalphiDot_dPhi + integer :: & sourceOffset sourceOffset = source_damage_anisoBrittle_offset(phase) - - localphiDot = 1.0_pReal & - - sourceState(phase)%p(sourceOffset)%state(1,constituent)*phi - + dLocalphiDot_dPhi = -sourceState(phase)%p(sourceOffset)%state(1,constituent) - + + localphiDot = 1.0_pReal & + + dLocalphiDot_dPhi*phi + end subroutine source_damage_anisoBrittle_getRateAndItsTangent @@ -256,21 +209,20 @@ end subroutine source_damage_anisoBrittle_getRateAndItsTangent !-------------------------------------------------------------------------------------------------- subroutine source_damage_anisoBrittle_results(phase,group) - integer, intent(in) :: phase - character(len=*), intent(in) :: group - integer :: sourceOffset, o, instance - - instance = source_damage_anisoBrittle_instance(phase) - sourceOffset = source_damage_anisoBrittle_offset(phase) + integer, intent(in) :: phase + character(len=*), intent(in) :: group - associate(prm => param(instance), stt => sourceState(phase)%p(sourceOffset)%state) - outputsLoop: do o = 1,size(prm%outputID) - select case(prm%outputID(o)) - case (damage_drivingforce_ID) - call results_writeDataset(group,stt,'tbd','driving force','tbd') - end select - enddo outputsLoop - end associate + integer :: o + + associate(prm => param(source_damage_anisoBrittle_instance(phase)), & + stt => sourceState(phase)%p(source_damage_anisoBrittle_offset(phase))%state) + outputsLoop: do o = 1,size(prm%output) + select case(trim(prm%output(o))) + case ('anisobrittle_drivingforce') + call results_writeDataset(group,stt,'tbd','driving force','tbd') + end select + enddo outputsLoop + end associate end subroutine source_damage_anisoBrittle_results diff --git a/src/source_damage_anisoDuctile.f90 b/src/source_damage_anisoDuctile.f90 index fef897914..4dd18dac0 100644 --- a/src/source_damage_anisoDuctile.f90 +++ b/src/source_damage_anisoDuctile.f90 @@ -13,21 +13,15 @@ module source_damage_anisoDuctile use material use config use results - + implicit none private - + integer, dimension(:), allocatable :: & source_damage_anisoDuctile_offset, & !< which source is my current damage mechanism? source_damage_anisoDuctile_instance !< instance of damage source mechanism - - enum, bind(c) - enumerator :: undefined_ID, & - damage_drivingforce_ID - end enum - - - type, private :: tParameters !< container type for internal constitutive parameters + + type, private :: tParameters !< container type for internal constitutive parameters real(pReal) :: & aTol, & N @@ -37,13 +31,13 @@ module source_damage_anisoDuctile totalNslip integer, dimension(:), allocatable :: & Nslip - integer(kind(undefined_ID)), allocatable, dimension(:) :: & - outputID + character(len=pStringLen), allocatable, dimension(:) :: & + output end type tParameters - - type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance) - - + + type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance) + + public :: & source_damage_anisoDuctile_init, & source_damage_anisoDuctile_dotState, & @@ -58,93 +52,64 @@ contains !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- subroutine source_damage_anisoDuctile_init - - integer :: Ninstance,phase,instance,source,sourceOffset - integer :: NofMyPhase,p ,i - - integer(kind(undefined_ID)) :: & - outputID - - character(len=pStringLen) :: & - extmsg = '' - character(len=pStringLen), dimension(:), allocatable :: & - outputs - - write(6,'(/,a)') ' <<<+- source_'//SOURCE_DAMAGE_ANISODUCTILE_LABEL//' init -+>>>'; flush(6) - - Ninstance = count(phase_source == SOURCE_damage_anisoDuctile_ID) - if (Ninstance == 0) return - + + integer :: Ninstance,sourceOffset,NofMyPhase,p + character(len=pStringLen) :: extmsg = '' + + write(6,'(/,a)') ' <<<+- source_'//SOURCE_DAMAGE_ANISODUCTILE_LABEL//' init -+>>>'; flush(6) + + Ninstance = count(phase_source == SOURCE_DAMAGE_ANISODUCTILE_ID) if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) & write(6,'(a16,1x,i5,/)') '# instances:',Ninstance - - allocate(source_damage_anisoDuctile_offset(size(config_phase)), source=0) + + allocate(source_damage_anisoDuctile_offset (size(config_phase)), source=0) allocate(source_damage_anisoDuctile_instance(size(config_phase)), source=0) - do phase = 1, size(config_phase) - source_damage_anisoDuctile_instance(phase) = count(phase_source(:,1:phase) == source_damage_anisoDuctile_ID) - do source = 1, phase_Nsources(phase) - if (phase_source(source,phase) == source_damage_anisoDuctile_ID) & - source_damage_anisoDuctile_offset(phase) = source - enddo - enddo - allocate(param(Ninstance)) - - do p=1, size(config_phase) + + do p = 1, size(config_phase) + source_damage_anisoDuctile_instance(p) = count(phase_source(:,1:p) == SOURCE_DAMAGE_ANISODUCTILE_ID) + do sourceOffset = 1, phase_Nsources(p) + if (phase_source(sourceOffset,p) == SOURCE_DAMAGE_ANISODUCTILE_ID) then + source_damage_anisoDuctile_offset(p) = sourceOffset + exit + endif + enddo + if (all(phase_source(:,p) /= SOURCE_DAMAGE_ANISODUCTILE_ID)) cycle associate(prm => param(source_damage_anisoDuctile_instance(p)), & config => config_phase(p)) - - prm%aTol = config%getFloat('anisoductile_atol',defaultVal = 1.0e-3_pReal) - - prm%N = config%getFloat('anisoductile_ratesensitivity') + + prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray) prm%totalNslip = sum(prm%Nslip) - ! sanity checks - if (prm%aTol < 0.0_pReal) extmsg = trim(extmsg)//' anisoductile_atol' - - if (prm%N <= 0.0_pReal) extmsg = trim(extmsg)//' anisoductile_ratesensitivity' - - prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray) - + + prm%aTol = config%getFloat('anisoductile_atol',defaultVal = 1.0e-3_pReal) + prm%N = config%getFloat('anisoductile_ratesensitivity') + prm%critPlasticStrain = config%getFloats('anisoductile_criticalplasticstrain',requiredSize=size(prm%Nslip)) - - ! expand: family => system - prm%critPlasticStrain = math_expand(prm%critPlasticStrain, prm%Nslip) - - if (any(prm%critPlasticStrain < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_criticalplasticstrain' - + + ! expand: family => system + prm%critPlasticStrain = math_expand(prm%critPlasticStrain, prm%Nslip) + + ! sanity checks + if (prm%aTol < 0.0_pReal) extmsg = trim(extmsg)//' anisoductile_atol' + if (prm%N <= 0.0_pReal) extmsg = trim(extmsg)//' anisoductile_ratesensitivity' + if (any(prm%critPlasticStrain < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_criticalplasticstrain' + !-------------------------------------------------------------------------------------------------- ! exit if any parameter is out of range if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'('//SOURCE_DAMAGE_ANISODUCTILE_LABEL//')') !-------------------------------------------------------------------------------------------------- ! output pararameters - outputs = config%getStrings('(output)',defaultVal=emptyStringArray) - allocate(prm%outputID(0)) - do i=1, size(outputs) - outputID = undefined_ID - select case(outputs(i)) - - case ('anisoductile_drivingforce') - prm%outputID = [prm%outputID, damage_drivingforce_ID] - - end select - - enddo - + prm%output = config%getStrings('(output)',defaultVal=emptyStringArray) + + NofMyPhase=count(material_phaseAt==p) * discretization_nIP + call material_allocateSourceState(p,sourceOffset,NofMyPhase,1,1,0) + sourceState(p)%p(sourceOffset)%aTolState=prm%aTol + end associate - - phase = p - - NofMyPhase=count(material_phaseAt==phase) * discretization_nIP - instance = source_damage_anisoDuctile_instance(phase) - sourceOffset = source_damage_anisoDuctile_offset(phase) - - call material_allocateSourceState(phase,sourceOffset,NofMyPhase,1,1,0) - sourceState(phase)%p(sourceOffset)%aTolState=param(instance)%aTol - enddo - + end subroutine source_damage_anisoDuctile_init @@ -157,29 +122,29 @@ subroutine source_damage_anisoDuctile_dotState(ipc, ip, el) ipc, & !< component-ID of integration point ip, & !< integration point el !< element + integer :: & phase, & constituent, & sourceOffset, & - homog, damageOffset, & - instance, & + damageOffset, & + homog, & i - + phase = material_phaseAt(ipc,el) constituent = material_phasememberAt(ipc,ip,el) - instance = source_damage_anisoDuctile_instance(phase) sourceOffset = source_damage_anisoDuctile_offset(phase) homog = material_homogenizationAt(el) damageOffset = damageMapping(homog)%p(ip,el) - - - do i = 1, param(instance)%totalNslip - sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = & - sourceState(phase)%p(sourceOffset)%dotState(1,constituent) + & - plasticState(phase)%slipRate(i,constituent)/ & - ((damage(homog)%p(damageOffset))**param(instance)%N)/param(instance)%critPlasticStrain(i) + + associate(prm => param(source_damage_anisoDuctile_instance(phase))) + do i = 1, prm%totalNslip + sourceState(phase)%p(sourceOffset)%dotState(1,constituent) & + = sourceState(phase)%p(sourceOffset)%dotState(1,constituent) & + + plasticState(phase)%slipRate(i,constituent)/(damage(homog)%p(damageOffset)**prm%N)/prm%critPlasticStrain(i) enddo - + end associate + end subroutine source_damage_anisoDuctile_dotState @@ -196,16 +161,17 @@ subroutine source_damage_anisoDuctile_getRateAndItsTangent(localphiDot, dLocalph real(pReal), intent(out) :: & localphiDot, & dLocalphiDot_dPhi + integer :: & sourceOffset - + sourceOffset = source_damage_anisoDuctile_offset(phase) - - localphiDot = 1.0_pReal & - - sourceState(phase)%p(sourceOffset)%state(1,constituent) * phi - + dLocalphiDot_dPhi = -sourceState(phase)%p(sourceOffset)%state(1,constituent) - + + localphiDot = 1.0_pReal & + + dLocalphiDot_dPhi*phi + end subroutine source_damage_anisoDuctile_getRateAndItsTangent @@ -214,21 +180,20 @@ end subroutine source_damage_anisoDuctile_getRateAndItsTangent !-------------------------------------------------------------------------------------------------- subroutine source_damage_anisoDuctile_results(phase,group) - integer, intent(in) :: phase + integer, intent(in) :: phase character(len=*), intent(in) :: group - integer :: sourceOffset, o, instance - - instance = source_damage_anisoDuctile_instance(phase) - sourceOffset = source_damage_anisoDuctile_offset(phase) - associate(prm => param(instance), stt => sourceState(phase)%p(sourceOffset)%state) - outputsLoop: do o = 1,size(prm%outputID) - select case(prm%outputID(o)) - case (damage_drivingforce_ID) - call results_writeDataset(group,stt,'tbd','driving force','tbd') - end select - enddo outputsLoop - end associate + integer :: o + + associate(prm => param(source_damage_anisoDuctile_instance(phase)), & + stt => sourceState(phase)%p(source_damage_anisoDuctile_offset(phase))%state) + outputsLoop: do o = 1,size(prm%output) + select case(trim(prm%output(o))) + case ('anisoductile_drivingforce') + call results_writeDataset(group,stt,'tbd','driving force','tbd') + end select + enddo outputsLoop + end associate end subroutine source_damage_anisoDuctile_results diff --git a/src/source_damage_isoBrittle.f90 b/src/source_damage_isoBrittle.f90 index 53c0b77a7..8856dc4e0 100644 --- a/src/source_damage_isoBrittle.f90 +++ b/src/source_damage_isoBrittle.f90 @@ -16,34 +16,28 @@ module source_damage_isoBrittle implicit none private + integer, dimension(:), allocatable :: & source_damage_isoBrittle_offset, & source_damage_isoBrittle_instance - enum, bind(c) - enumerator :: & - undefined_ID, & - damage_drivingforce_ID - end enum - - type, private :: tParameters !< container type for internal constitutive parameters real(pReal) :: & critStrainEnergy, & N, & aTol - integer(kind(undefined_ID)), allocatable, dimension(:) :: & - outputID + character(len=pStringLen), allocatable, dimension(:) :: & + output end type tParameters - type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance) + 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 + source_damage_isoBrittle_results contains @@ -54,52 +48,41 @@ contains !-------------------------------------------------------------------------------------------------- subroutine source_damage_isoBrittle_init - integer :: Ninstance,phase,instance,source,sourceOffset - integer :: NofMyPhase,p,i - integer(kind(undefined_ID)) :: & - outputID - - character(len=pStringLen) :: & - extmsg = '' - character(len=pStringLen), dimension(:), allocatable :: & - outputs - - write(6,'(/,a)') ' <<<+- source_'//SOURCE_DAMAGE_ISOBRITTLE_LABEL//' init -+>>>'; flush(6) - - Ninstance = count(phase_source == SOURCE_damage_isoBrittle_ID) - if (Ninstance == 0) return - + integer :: Ninstance,sourceOffset,NofMyPhase,p + character(len=pStringLen) :: extmsg = '' + + write(6,'(/,a)') ' <<<+- source_'//SOURCE_DAMAGE_ISOBRITTLE_LABEL//' init -+>>>'; flush(6) + + Ninstance = count(phase_source == SOURCE_DAMAGE_ISOBRITTLE_ID) if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) & write(6,'(a16,1x,i5,/)') '# instances:',Ninstance - - allocate(source_damage_isoBrittle_offset(material_Nphase), source=0) - allocate(source_damage_isoBrittle_instance(material_Nphase), source=0) - do phase = 1, material_Nphase - source_damage_isoBrittle_instance(phase) = count(phase_source(:,1:phase) == source_damage_isoBrittle_ID) - do source = 1, phase_Nsources(phase) - if (phase_source(source,phase) == source_damage_isoBrittle_ID) & - source_damage_isoBrittle_offset(phase) = source - enddo - enddo - + + allocate(source_damage_isoBrittle_offset (size(config_phase)), source=0) + allocate(source_damage_isoBrittle_instance(size(config_phase)), source=0) allocate(param(Ninstance)) - - do p=1, size(config_phase) + + do p = 1, size(config_phase) + source_damage_isoBrittle_instance(p) = count(phase_source(:,1:p) == SOURCE_DAMAGE_ISOBRITTLE_ID) + do sourceOffset = 1, phase_Nsources(p) + if (phase_source(sourceOffset,p) == SOURCE_DAMAGE_ISOBRITTLE_ID) then + source_damage_isoBrittle_offset(p) = sourceOffset + exit + endif + enddo + if (all(phase_source(:,p) /= SOURCE_DAMAGE_ISOBRITTLE_ID)) cycle associate(prm => param(source_damage_isoBrittle_instance(p)), & config => config_phase(p)) - + prm%aTol = config%getFloat('isobrittle_atol',defaultVal = 1.0e-3_pReal) - prm%N = config%getFloat('isobrittle_n') prm%critStrainEnergy = config%getFloat('isobrittle_criticalstrainenergy') - + ! sanity checks - if (prm%aTol < 0.0_pReal) extmsg = trim(extmsg)//' isobrittle_atol' - - if (prm%N <= 0.0_pReal) extmsg = trim(extmsg)//' isobrittle_n' - if (prm%critStrainEnergy <= 0.0_pReal) extmsg = trim(extmsg)//' isobrittle_criticalstrainenergy' - + if (prm%aTol < 0.0_pReal) extmsg = trim(extmsg)//' isobrittle_atol' + if (prm%N <= 0.0_pReal) extmsg = trim(extmsg)//' isobrittle_n' + if (prm%critStrainEnergy <= 0.0_pReal) extmsg = trim(extmsg)//' isobrittle_criticalstrainenergy' + !-------------------------------------------------------------------------------------------------- ! exit if any parameter is out of range if (extmsg /= '') & @@ -107,34 +90,18 @@ subroutine source_damage_isoBrittle_init !-------------------------------------------------------------------------------------------------- ! output pararameters - outputs = config%getStrings('(output)',defaultVal=emptyStringArray) - allocate(prm%outputID(0)) - do i=1, size(outputs) - outputID = undefined_ID - select case(outputs(i)) - - case ('isobrittle_drivingforce') - prm%outputID = [prm%outputID, damage_drivingforce_ID] - - end select + prm%output = config%getStrings('(output)',defaultVal=emptyStringArray) - enddo + NofMyPhase = count(material_phaseAt==p) * discretization_nIP + call material_allocateSourceState(p,sourceOffset,NofMyPhase,1,1,1) + sourceState(p)%p(sourceOffset)%aTolState=prm%aTol end associate - - phase = p - - NofMyPhase = count(material_phaseAt==phase) * discretization_nIP - instance = source_damage_isoBrittle_instance(phase) - sourceOffset = source_damage_isoBrittle_offset(phase) - - call material_allocateSourceState(phase,sourceOffset,NofMyPhase,1,1,1) - sourceState(phase)%p(sourceOffset)%aTolState=param(instance)%aTol - enddo - + end subroutine source_damage_isoBrittle_init + !-------------------------------------------------------------------------------------------------- !> @brief calculates derived quantities from state !-------------------------------------------------------------------------------------------------- @@ -148,24 +115,26 @@ subroutine source_damage_isoBrittle_deltaState(C, Fe, ipc, ip, el) Fe real(pReal), intent(in), dimension(6,6) :: & C + integer :: & - phase, constituent, instance, sourceOffset + phase, & + constituent, & + sourceOffset + real(pReal), dimension(6) :: & + strain real(pReal) :: & - strain(6), & strainenergy - phase = material_phaseAt(ipc,el) !< phase ID at ipc,ip,el - constituent = material_phasememberAt(ipc,ip,el) !< state array offset for phase ID at ipc,ip,el - ! ToDo: capability for multiple instances of SAME source within given phase. Needs Ninstance loop from here on! - instance = source_damage_isoBrittle_instance(phase) !< instance of damage_isoBrittle source + phase = material_phaseAt(ipc,el) !< phase ID at ipc,ip,el + constituent = material_phasememberAt(ipc,ip,el) !< state array offset for phase ID at ipc,ip,el sourceOffset = source_damage_isoBrittle_offset(phase) - strain = 0.5_pReal*math_sym33to6(matmul(transpose(Fe),Fe)-math_I3) - strainenergy = 2.0_pReal*sum(strain*matmul(C,strain))/param(instance)%critStrainEnergy + associate(prm => param(source_damage_isoBrittle_instance(phase))) + strainenergy = 2.0_pReal*sum(strain*matmul(C,strain))/prm%critStrainEnergy ! ToDo: check strainenergy = 2.0_pReal*dot_product(strain,matmul(C,strain))/param(instance)%critStrainEnergy - + if (strainenergy > sourceState(phase)%p(sourceOffset)%subState0(1,constituent)) then sourceState(phase)%p(sourceOffset)%deltaState(1,constituent) = & strainenergy - sourceState(phase)%p(sourceOffset)%state(1,constituent) @@ -174,9 +143,11 @@ subroutine source_damage_isoBrittle_deltaState(C, Fe, ipc, ip, el) sourceState(phase)%p(sourceOffset)%subState0(1,constituent) - & sourceState(phase)%p(sourceOffset)%state(1,constituent) endif - + end associate + end subroutine source_damage_isoBrittle_deltaState - + + !-------------------------------------------------------------------------------------------------- !> @brief returns local part of nonlocal damage driving force !-------------------------------------------------------------------------------------------------- @@ -190,18 +161,19 @@ subroutine source_damage_isoBrittle_getRateAndItsTangent(localphiDot, dLocalphiD real(pReal), intent(out) :: & localphiDot, & dLocalphiDot_dPhi - integer :: & - instance, sourceOffset - instance = source_damage_isoBrittle_instance(phase) + integer :: & + sourceOffset + sourceOffset = source_damage_isoBrittle_offset(phase) - - localphiDot = (1.0_pReal - phi)**(param(instance)%N - 1.0_pReal) - & - phi*sourceState(phase)%p(sourceOffset)%state(1,constituent) - dLocalphiDot_dPhi = - (param(instance)%N - 1.0_pReal)* & - (1.0_pReal - phi)**max(0.0_pReal,param(instance)%N - 2.0_pReal) & + + associate(prm => param(source_damage_isoBrittle_instance(phase))) + localphiDot = (1.0_pReal - phi)**(prm%n - 1.0_pReal) & + - phi*sourceState(phase)%p(sourceOffset)%state(1,constituent) + dLocalphiDot_dPhi = - (prm%n - 1.0_pReal)* (1.0_pReal - phi)**max(0.0_pReal,prm%n - 2.0_pReal) & - sourceState(phase)%p(sourceOffset)%state(1,constituent) - + end associate + end subroutine source_damage_isoBrittle_getRateAndItsTangent @@ -210,21 +182,20 @@ end subroutine source_damage_isoBrittle_getRateAndItsTangent !-------------------------------------------------------------------------------------------------- subroutine source_damage_isoBrittle_results(phase,group) - integer, intent(in) :: phase - character(len=*), intent(in) :: group - integer :: sourceOffset, o, instance - - instance = source_damage_isoBrittle_instance(phase) - sourceOffset = source_damage_isoBrittle_offset(phase) + integer, intent(in) :: phase + character(len=*), intent(in) :: group - associate(prm => param(instance), stt => sourceState(phase)%p(sourceOffset)%state) - outputsLoop: do o = 1,size(prm%outputID) - select case(prm%outputID(o)) - case (damage_drivingforce_ID) - call results_writeDataset(group,stt,'tbd','driving force','tbd') - end select - enddo outputsLoop - end associate + integer :: o + + associate(prm => param(source_damage_isoBrittle_instance(phase)), & + stt => sourceState(phase)%p(source_damage_isoBrittle_offset(phase))%state) + outputsLoop: do o = 1,size(prm%output) + select case(trim(prm%output(o))) + case ('isobrittle_drivingforce') + call results_writeDataset(group,stt,'tbd','driving force','tbd') + end select + enddo outputsLoop + end associate end subroutine source_damage_isoBrittle_results diff --git a/src/source_damage_isoDuctile.f90 b/src/source_damage_isoDuctile.f90 index 6ee588d0c..2f7389592 100644 --- a/src/source_damage_isoDuctile.f90 +++ b/src/source_damage_isoDuctile.f90 @@ -5,42 +5,38 @@ !> @details to be done !-------------------------------------------------------------------------------------------------- module source_damage_isoDuctile - use prec - use debug - use IO - use discretization - use material - use config - use results + use prec + use debug + use IO + use discretization + use material + use config + use results - implicit none - private - integer, dimension(:), allocatable :: & - source_damage_isoDuctile_offset, & !< which source is my current damage mechanism? - source_damage_isoDuctile_instance !< instance of damage source mechanism + implicit none + private - enum, bind(c) - enumerator :: undefined_ID, & - damage_drivingforce_ID - end enum !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ToDo + 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 - real(pReal) :: & - critPlasticStrain, & - N, & - aTol - integer(kind(undefined_ID)), allocatable, dimension(:) :: & - outputID - end type tParameters + type, private :: tParameters !< container type for internal constitutive parameters + real(pReal) :: & + critPlasticStrain, & + N, & + aTol + character(len=pStringLen), allocatable, dimension(:) :: & + output + end type tParameters - type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance) + type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance) - public :: & - source_damage_isoDuctile_init, & - source_damage_isoDuctile_dotState, & - source_damage_isoDuctile_getRateAndItsTangent, & - source_damage_isoDuctile_Results + public :: & + source_damage_isoDuctile_init, & + source_damage_isoDuctile_dotState, & + source_damage_isoDuctile_getRateAndItsTangent, & + source_damage_isoDuctile_Results contains @@ -51,135 +47,115 @@ contains !-------------------------------------------------------------------------------------------------- subroutine source_damage_isoDuctile_init - integer :: Ninstance,phase,instance,source,sourceOffset - integer :: NofMyPhase,p,i - integer(kind(undefined_ID)) :: & - outputID + integer :: Ninstance,sourceOffset,NofMyPhase,p + character(len=pStringLen) :: extmsg = '' - character(len=pStringLen) :: & - extmsg = '' - character(len=pStringLen), dimension(:), allocatable :: & - outputs + write(6,'(/,a)') ' <<<+- source_'//SOURCE_DAMAGE_ISODUCTILE_LABEL//' init -+>>>'; flush(6) - write(6,'(/,a)') ' <<<+- source_'//SOURCE_DAMAGE_ISODUCTILE_LABEL//' init -+>>>' + Ninstance = count(phase_source == SOURCE_DAMAGE_ISODUCTILE_ID) + if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) & + write(6,'(a16,1x,i5,/)') '# instances:',Ninstance - Ninstance = count(phase_source == SOURCE_damage_isoDuctile_ID) - if (Ninstance == 0) return - - if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) & - write(6,'(a16,1x,i5,/)') '# instances:',Ninstance - - allocate(source_damage_isoDuctile_offset(material_Nphase), source=0) - allocate(source_damage_isoDuctile_instance(material_Nphase), source=0) - do phase = 1, material_Nphase - source_damage_isoDuctile_instance(phase) = count(phase_source(:,1:phase) == source_damage_isoDuctile_ID) - do source = 1, phase_Nsources(phase) - if (phase_source(source,phase) == source_damage_isoDuctile_ID) & - source_damage_isoDuctile_offset(phase) = source - enddo - enddo + allocate(source_damage_isoDuctile_offset (size(config_phase)), source=0) + allocate(source_damage_isoDuctile_instance(size(config_phase)), source=0) + allocate(param(Ninstance)) - allocate(param(Ninstance)) - - do p=1, size(config_phase) - if (all(phase_source(:,p) /= SOURCE_DAMAGE_ISODUCTILE_ID)) cycle - associate(prm => param(source_damage_isoDuctile_instance(p)), & - config => config_phase(p)) - - prm%aTol = config%getFloat('isoductile_atol',defaultVal = 1.0e-3_pReal) + do p = 1, size(config_phase) + source_damage_isoDuctile_instance(p) = count(phase_source(:,1:p) == SOURCE_DAMAGE_ISODUCTILE_ID) + do sourceOffset = 1, phase_Nsources(p) + if (phase_source(sourceOffset,p) == SOURCE_DAMAGE_ISODUCTILE_ID) then + source_damage_isoDuctile_offset(p) = sourceOffset + exit + endif + enddo + + if (all(phase_source(:,p) /= SOURCE_DAMAGE_ISODUCTILE_ID)) cycle + associate(prm => param(source_damage_isoDuctile_instance(p)), & + config => config_phase(p)) + + prm%aTol = config%getFloat('isoductile_atol',defaultVal = 1.0e-3_pReal) + prm%N = config%getFloat('isoductile_ratesensitivity') + prm%critPlasticStrain = config%getFloat('isoductile_criticalplasticstrain') + + ! sanity checks + if (prm%aTol < 0.0_pReal) extmsg = trim(extmsg)//' isoductile_atol' + if (prm%N <= 0.0_pReal) extmsg = trim(extmsg)//' isoductile_ratesensitivity' + if (prm%critPlasticStrain <= 0.0_pReal) extmsg = trim(extmsg)//' isoductile_criticalplasticstrain' - prm%N = config%getFloat('isoductile_ratesensitivity') - prm%critPlasticStrain = config%getFloat('isoductile_criticalplasticstrain') - - ! sanity checks - if (prm%aTol < 0.0_pReal) extmsg = trim(extmsg)//' isoductile_atol' - - if (prm%N <= 0.0_pReal) extmsg = trim(extmsg)//' isoductile_ratesensitivity' - if (prm%critPlasticStrain <= 0.0_pReal) extmsg = trim(extmsg)//' isoductile_criticalplasticstrain' - !-------------------------------------------------------------------------------------------------- ! exit if any parameter is out of range - if (extmsg /= '') & - call IO_error(211,ext_msg=trim(extmsg)//'('//SOURCE_DAMAGE_ISODUCTILE_LABEL//')') + if (extmsg /= '') & + call IO_error(211,ext_msg=trim(extmsg)//'('//SOURCE_DAMAGE_ISODUCTILE_LABEL//')') !-------------------------------------------------------------------------------------------------- ! output pararameters - outputs = config%getStrings('(output)',defaultVal=emptyStringArray) - allocate(prm%outputID(0)) - do i=1, size(outputs) - outputID = undefined_ID - select case(outputs(i)) - - case ('isoductile_drivingforce') - prm%outputID = [prm%outputID, damage_drivingforce_ID] + prm%output = config%getStrings('(output)',defaultVal=emptyStringArray) - end select + NofMyPhase=count(material_phaseAt==p) * discretization_nIP + call material_allocateSourceState(p,sourceOffset,NofMyPhase,1,1,0) + sourceState(p)%p(sourceOffset)%aTolState=prm%aTol - enddo + end associate + enddo - end associate - - phase = p - NofMyPhase=count(material_phaseAt==phase) * discretization_nIP - instance = source_damage_isoDuctile_instance(phase) - sourceOffset = source_damage_isoDuctile_offset(phase) - - call material_allocateSourceState(phase,sourceOffset,NofMyPhase,1,1,0) - sourceState(phase)%p(sourceOffset)%aTolState=param(instance)%aTol - - enddo - end subroutine source_damage_isoDuctile_init + !-------------------------------------------------------------------------------------------------- !> @brief calculates derived quantities from state !-------------------------------------------------------------------------------------------------- subroutine source_damage_isoDuctile_dotState(ipc, ip, el) - integer, intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element - integer :: & - phase, constituent, instance, homog, sourceOffset, damageOffset + integer, intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element - phase = material_phaseAt(ipc,el) - constituent = material_phasememberAt(ipc,ip,el) - instance = source_damage_isoDuctile_instance(phase) - sourceOffset = source_damage_isoDuctile_offset(phase) - homog = material_homogenizationAt(el) - damageOffset = damageMapping(homog)%p(ip,el) + integer :: & + phase, & + constituent, & + sourceOffset, & + damageOffset, & + homog - sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = & - sum(plasticState(phase)%slipRate(:,constituent))/ & - ((damage(homog)%p(damageOffset))**param(instance)%N)/ & - param(instance)%critPlasticStrain + phase = material_phaseAt(ipc,el) + constituent = material_phasememberAt(ipc,ip,el) + sourceOffset = source_damage_isoDuctile_offset(phase) + homog = material_homogenizationAt(el) + damageOffset = damageMapping(homog)%p(ip,el) + + associate(prm => param(source_damage_isoDuctile_instance(phase))) + sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = & + sum(plasticState(phase)%slipRate(:,constituent))/(damage(homog)%p(damageOffset)**prm%N)/prm%critPlasticStrain + end associate 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) - integer, intent(in) :: & - phase, & - constituent - real(pReal), intent(in) :: & - phi - real(pReal), intent(out) :: & - localphiDot, & - dLocalphiDot_dPhi - integer :: & - sourceOffset + integer, intent(in) :: & + phase, & + constituent + real(pReal), intent(in) :: & + phi + real(pReal), intent(out) :: & + localphiDot, & + dLocalphiDot_dPhi + + integer :: & + sourceOffset + + sourceOffset = source_damage_isoDuctile_offset(phase) + + dLocalphiDot_dPhi = -sourceState(phase)%p(sourceOffset)%state(1,constituent) + + localphiDot = 1.0_pReal & + + dLocalphiDot_dPhi*phi - sourceOffset = source_damage_isoDuctile_offset(phase) - - localphiDot = 1.0_pReal & - - sourceState(phase)%p(sourceOffset)%state(1,constituent) * phi - - dLocalphiDot_dPhi = -sourceState(phase)%p(sourceOffset)%state(1,constituent) - end subroutine source_damage_isoDuctile_getRateAndItsTangent @@ -188,23 +164,21 @@ end subroutine source_damage_isoDuctile_getRateAndItsTangent !-------------------------------------------------------------------------------------------------- subroutine source_damage_isoDuctile_results(phase,group) - integer, intent(in) :: phase - character(len=*), intent(in) :: group - integer :: sourceOffset, o, instance - - instance = source_damage_isoDuctile_instance(phase) - sourceOffset = source_damage_isoDuctile_offset(phase) + integer, intent(in) :: phase + character(len=*), intent(in) :: group - associate(prm => param(instance), stt => sourceState(phase)%p(sourceOffset)%state) - outputsLoop: do o = 1,size(prm%outputID) - select case(prm%outputID(o)) - case (damage_drivingforce_ID) - call results_writeDataset(group,stt,'tbd','driving force','tbd') - end select - enddo outputsLoop - end associate + integer :: o + + associate(prm => param(source_damage_isoDuctile_instance(phase)), & + stt => sourceState(phase)%p(source_damage_isoDuctile_offset(phase))%state) + outputsLoop: do o = 1,size(prm%output) + select case(trim(prm%output(o))) + case ('isoductile_drivingforce') + call results_writeDataset(group,stt,'tbd','driving force','tbd') + end select + enddo outputsLoop + end associate end subroutine source_damage_isoDuctile_results - end module source_damage_isoDuctile diff --git a/src/source_thermal_dissipation.f90 b/src/source_thermal_dissipation.f90 index e13742a90..c03d66c24 100644 --- a/src/source_thermal_dissipation.f90 +++ b/src/source_thermal_dissipation.f90 @@ -10,26 +10,26 @@ module source_thermal_dissipation use discretization use material use config - + implicit none private integer, dimension(:), allocatable :: & source_thermal_dissipation_offset, & !< which source is my current thermal dissipation mechanism? source_thermal_dissipation_instance !< instance of thermal dissipation source mechanism - + type :: tParameters !< container type for internal constitutive parameters real(pReal) :: & kappa end type tParameters - + type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance) - - + + public :: & source_thermal_dissipation_init, & source_thermal_dissipation_getRateAndItsTangent - + contains @@ -38,47 +38,47 @@ contains !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- subroutine source_thermal_dissipation_init - - integer :: Ninstance,instance,source,sourceOffset,NofMyPhase,p - - write(6,'(/,a)') ' <<<+- source_'//SOURCE_thermal_dissipation_label//' init -+>>>'; flush(6) - - - Ninstance = count(phase_source == SOURCE_thermal_dissipation_ID) - if (Ninstance == 0) return + + integer :: Ninstance,sourceOffset,NofMyPhase,p + + write(6,'(/,a)') ' <<<+- source_'//SOURCE_thermal_dissipation_label//' init -+>>>'; flush(6) + + Ninstance = count(phase_source == SOURCE_THERMAL_DISSIPATION_ID) if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) & write(6,'(a16,1x,i5,/)') '# instances:',Ninstance - - allocate(source_thermal_dissipation_offset(material_Nphase), source=0) - allocate(source_thermal_dissipation_instance(material_Nphase), source=0) + + allocate(source_thermal_dissipation_offset (size(config_phase)), source=0) + allocate(source_thermal_dissipation_instance(size(config_phase)), source=0) allocate(param(Ninstance)) - - do p = 1, material_Nphase - source_thermal_dissipation_instance(p) = count(phase_source(:,1:p) == SOURCE_thermal_dissipation_ID) - do source = 1, phase_Nsources(p) - if (phase_source(source,p) == SOURCE_thermal_dissipation_ID) & - source_thermal_dissipation_offset(p) = source - enddo - enddo - - do p=1, size(config_phase) + + do p = 1, size(config_phase) + source_thermal_dissipation_instance(p) = count(phase_source(:,1:p) == SOURCE_THERMAL_DISSIPATION_ID) + do sourceOffset = 1, phase_Nsources(p) + if (phase_source(sourceOffset,p) == SOURCE_THERMAL_DISSIPATION_ID) then + source_thermal_dissipation_offset(p) = sourceOffset + exit + endif + enddo + if (all(phase_source(:,p) /= SOURCE_THERMAL_DISSIPATION_ID)) cycle - instance = source_thermal_dissipation_instance(p) - param(instance)%kappa = config_phase(p)%getFloat('dissipation_coldworkcoeff') + associate(prm => param(source_thermal_dissipation_instance(p)), & + config => config_phase(p)) + + prm%kappa = config%getFloat('dissipation_coldworkcoeff') + NofMyPhase = count(material_phaseAt==p) * discretization_nIP - sourceOffset = source_thermal_dissipation_offset(p) - call material_allocateSourceState(p,sourceOffset,NofMyPhase,0,0,0) - + + end associate enddo - + end subroutine source_thermal_dissipation_init !-------------------------------------------------------------------------------------------------- -!> @brief returns dissipation rate +!> @brief Ninstances dissipation rate !-------------------------------------------------------------------------------------------------- -subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDOT_dT, Tstar, Lp, phase) +subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDot_dT, Tstar, Lp, phase) integer, intent(in) :: & phase @@ -86,17 +86,16 @@ subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDOT_dT, Tstar Tstar real(pReal), intent(in), dimension(3,3) :: & Lp + real(pReal), intent(out) :: & TDot, & - dTDOT_dT - integer :: & - instance - - instance = source_thermal_dissipation_instance(phase) - - TDot = param(instance)%kappa*sum(abs(Tstar*Lp)) - dTDOT_dT = 0.0_pReal - + dTDot_dT + + associate(prm => param(source_thermal_dissipation_instance(phase))) + TDot = prm%kappa*sum(abs(Tstar*Lp)) + dTDot_dT = 0.0_pReal + end associate + end subroutine source_thermal_dissipation_getRateAndItsTangent - + end module source_thermal_dissipation diff --git a/src/source_thermal_externalheat.f90 b/src/source_thermal_externalheat.f90 index ccc6ebc29..6cb8cdf6e 100644 --- a/src/source_thermal_externalheat.f90 +++ b/src/source_thermal_externalheat.f90 @@ -18,7 +18,7 @@ module source_thermal_externalheat source_thermal_externalheat_offset, & !< which source is my current thermal dissipation mechanism? source_thermal_externalheat_instance !< instance of thermal dissipation source mechanism - type :: tParameters !< container type for internal constitutive parameters + type :: tParameters !< container type for internal constitutive parameters real(pReal), dimension(:), allocatable :: & time, & heat_rate @@ -42,44 +42,41 @@ contains !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- subroutine source_thermal_externalheat_init - - integer :: maxNinstance,instance,source,sourceOffset,NofMyPhase,p - - write(6,'(/,a)') ' <<<+- source_'//SOURCE_thermal_externalheat_label//' init -+>>>'; flush(6) - - - maxNinstance = count(phase_source == SOURCE_thermal_externalheat_ID) - if (maxNinstance == 0) return + + integer :: Ninstance,sourceOffset,NofMyPhase,p + + write(6,'(/,a)') ' <<<+- source_'//SOURCE_thermal_externalheat_label//' init -+>>>'; flush(6) + + Ninstance = count(phase_source == SOURCE_thermal_externalheat_ID) if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) & - write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance - - allocate(source_thermal_externalheat_offset(material_Nphase), source=0) - allocate(source_thermal_externalheat_instance(material_Nphase), source=0) - - do p = 1, material_Nphase + write(6,'(a16,1x,i5,/)') '# instances:',Ninstance + + allocate(source_thermal_externalheat_offset (size(config_phase)), source=0) + allocate(source_thermal_externalheat_instance(size(config_phase)), source=0) + allocate(param(Ninstance)) + + do p = 1, size(config_phase) source_thermal_externalheat_instance(p) = count(phase_source(:,1:p) == SOURCE_thermal_externalheat_ID) - do source = 1, phase_Nsources(p) - if (phase_source(source,p) == SOURCE_thermal_externalheat_ID) & - source_thermal_externalheat_offset(p) = source - enddo - enddo - - allocate(param(maxNinstance)) - - do p=1, size(config_phase) + do sourceOffset = 1, phase_Nsources(p) + if (phase_source(sourceOffset,p) == SOURCE_thermal_externalheat_ID) then + source_thermal_externalheat_offset(p) = sourceOffset + exit + endif + enddo + if (all(phase_source(:,p) /= SOURCE_thermal_externalheat_ID)) cycle - instance = source_thermal_externalheat_instance(p) - sourceOffset = source_thermal_externalheat_offset(p) + associate(prm => param(source_thermal_externalheat_instance(p)), & + config => config_phase(p)) + + prm%time = config%getFloats('externalheat_time') + prm%nIntervals = size(prm%time) - 1 + + prm%heat_rate = config%getFloats('externalheat_rate',requiredSize = size(prm%time)) + NofMyPhase = count(material_phaseAt==p) * discretization_nIP - - param(instance)%time = config_phase(p)%getFloats('externalheat_time') - param(instance)%nIntervals = size(param(instance)%time) - 1 - - - param(instance)%heat_rate = config_phase(p)%getFloats('externalheat_rate',requiredSize = size(param(instance)%time)) - call material_allocateSourceState(p,sourceOffset,NofMyPhase,1,1,0) - + + end associate enddo end subroutine source_thermal_externalheat_init @@ -90,51 +87,54 @@ end subroutine source_thermal_externalheat_init !> @details state only contains current time to linearly interpolate given heat powers !-------------------------------------------------------------------------------------------------- subroutine source_thermal_externalheat_dotState(phase, of) - + integer, intent(in) :: & phase, & of + integer :: & sourceOffset - + sourceOffset = source_thermal_externalheat_offset(phase) - + sourceState(phase)%p(sourceOffset)%dotState(1,of) = 1.0_pReal ! state is current time end subroutine source_thermal_externalheat_dotState !-------------------------------------------------------------------------------------------------- -!> @brief returns local heat generation rate +!> @brief returns local heat generation rate !-------------------------------------------------------------------------------------------------- subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, dTDot_dT, phase, of) - + integer, intent(in) :: & phase, & of real(pReal), intent(out) :: & TDot, & dTDot_dT + integer :: & - instance, sourceOffset, interval + sourceOffset, interval real(pReal) :: & - frac_time - - instance = source_thermal_externalheat_instance(phase) + frac_time + sourceOffset = source_thermal_externalheat_offset(phase) - - do interval = 1, param(instance)%nIntervals ! scan through all rate segments - frac_time = (sourceState(phase)%p(sourceOffset)%state(1,of) - param(instance)%time(interval)) & - / (param(instance)%time(interval+1) - param(instance)%time(interval)) ! fractional time within segment + + associate(prm => param(source_thermal_externalheat_instance(phase))) + do interval = 1, prm%nIntervals ! scan through all rate segments + frac_time = (sourceState(phase)%p(sourceOffset)%state(1,of) - prm%time(interval)) & + / (prm%time(interval+1) - prm%time(interval)) ! fractional time within segment if ( (frac_time < 0.0_pReal .and. interval == 1) & - .or. (frac_time >= 1.0_pReal .and. interval == param(instance)%nIntervals) & + .or. (frac_time >= 1.0_pReal .and. interval == prm%nIntervals) & .or. (frac_time >= 0.0_pReal .and. frac_time < 1.0_pReal) ) & - TDot = param(instance)%heat_rate(interval ) * (1.0_pReal - frac_time) + & - param(instance)%heat_rate(interval+1) * frac_time ! interpolate heat rate between segment boundaries... + TDot = prm%heat_rate(interval ) * (1.0_pReal - frac_time) + & + prm%heat_rate(interval+1) * frac_time ! interpolate heat rate between segment boundaries... ! ...or extrapolate if outside of bounds enddo dTDot_dT = 0.0 - + end associate + end subroutine source_thermal_externalheat_getRateAndItsTangent - + end module source_thermal_externalheat diff --git a/src/thermal_conduction.f90 b/src/thermal_conduction.f90 index dad3fecd8..be72c07b6 100644 --- a/src/thermal_conduction.f90 +++ b/src/thermal_conduction.f90 @@ -11,7 +11,7 @@ module thermal_conduction use crystallite use source_thermal_dissipation use source_thermal_externalheat - + implicit none private @@ -19,14 +19,14 @@ module thermal_conduction character(len=pStringLen), allocatable, dimension(:) :: & output end type tParameters - + type(tparameters), dimension(:), allocatable :: & - param - + param + public :: & thermal_conduction_init, & thermal_conduction_getSourceAndItsTangent, & - thermal_conduction_getConductivity33, & + thermal_conduction_getConductivity, & thermal_conduction_getSpecificHeat, & thermal_conduction_getMassDensity, & thermal_conduction_putTemperatureAndItsRate, & @@ -41,37 +41,34 @@ contains !-------------------------------------------------------------------------------------------------- subroutine thermal_conduction_init - - integer :: maxNinstance,NofMyHomog,h - + integer :: Ninstance,NofMyHomog,h + write(6,'(/,a)') ' <<<+- thermal_'//THERMAL_CONDUCTION_label//' init -+>>>'; flush(6) - - maxNinstance = count(thermal_type == THERMAL_conduction_ID) - if (maxNinstance == 0) return - - allocate(param(maxNinstance)) - - do h = 1, size(thermal_type) + + Ninstance = count(thermal_type == THERMAL_conduction_ID) + allocate(param(Ninstance)) + + do h = 1, size(config_homogenization) if (thermal_type(h) /= THERMAL_conduction_ID) cycle associate(prm => param(thermal_typeInstance(h)),config => config_homogenization(h)) - + prm%output = config%getStrings('(output)',defaultVal=emptyStringArray) - + NofMyHomog=count(material_homogenizationAt==h) thermalState(h)%sizeState = 0 allocate(thermalState(h)%state0 (0,NofMyHomog)) allocate(thermalState(h)%subState0(0,NofMyHomog)) allocate(thermalState(h)%state (0,NofMyHomog)) - + thermalMapping(h)%p => material_homogenizationMemberAt deallocate(temperature (h)%p) allocate (temperature (h)%p(NofMyHomog), source=thermal_initialT(h)) deallocate(temperatureRate(h)%p) allocate (temperatureRate(h)%p(NofMyHomog), source=0.0_pReal) - + end associate enddo - + end subroutine thermal_conduction_init @@ -79,7 +76,7 @@ end subroutine thermal_conduction_init !> @brief returns heat generation rate !-------------------------------------------------------------------------------------------------- subroutine thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el) - + integer, intent(in) :: & ip, & !< integration point number el !< element number @@ -97,74 +94,74 @@ subroutine thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el) grain, & source, & constituent - + homog = material_homogenizationAt(el) offset = material_homogenizationMemberAt(ip,el) instance = thermal_typeInstance(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)) + 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 enddo - + Tdot = Tdot/real(homogenization_Ngrains(homog),pReal) dTdot_dT = dTdot_dT/real(homogenization_Ngrains(homog),pReal) - + end subroutine thermal_conduction_getSourceAndItsTangent - + !-------------------------------------------------------------------------------------------------- !> @brief returns homogenized thermal conductivity in reference configuration !-------------------------------------------------------------------------------------------------- -function thermal_conduction_getConductivity33(ip,el) - +function thermal_conduction_getConductivity(ip,el) + integer, intent(in) :: & ip, & !< integration point number el !< element number real(pReal), dimension(3,3) :: & - thermal_conduction_getConductivity33 + thermal_conduction_getConductivity integer :: & grain - - - thermal_conduction_getConductivity33 = 0.0_pReal + + + thermal_conduction_getConductivity = 0.0_pReal do grain = 1, homogenization_Ngrains(material_homogenizationAt(el)) - thermal_conduction_getConductivity33 = thermal_conduction_getConductivity33 + & - crystallite_push33ToRef(grain,ip,el,lattice_thermalConductivity33(:,:,material_phaseAt(grain,el))) + thermal_conduction_getConductivity = thermal_conduction_getConductivity + & + crystallite_push33ToRef(grain,ip,el,lattice_thermalConductivity(:,:,material_phaseAt(grain,el))) enddo - - thermal_conduction_getConductivity33 = thermal_conduction_getConductivity33 & - / real(homogenization_Ngrains(material_homogenizationAt(el)),pReal) - -end function thermal_conduction_getConductivity33 + + thermal_conduction_getConductivity = thermal_conduction_getConductivity & + / real(homogenization_Ngrains(material_homogenizationAt(el)),pReal) + +end function thermal_conduction_getConductivity !-------------------------------------------------------------------------------------------------- !> @brief returns homogenized specific heat capacity !-------------------------------------------------------------------------------------------------- function thermal_conduction_getSpecificHeat(ip,el) - + integer, intent(in) :: & ip, & !< integration point number el !< element number @@ -172,17 +169,17 @@ function thermal_conduction_getSpecificHeat(ip,el) thermal_conduction_getSpecificHeat integer :: & grain - + thermal_conduction_getSpecificHeat = 0.0_pReal - + do grain = 1, homogenization_Ngrains(material_homogenizationAt(el)) thermal_conduction_getSpecificHeat = thermal_conduction_getSpecificHeat & + lattice_specificHeat(material_phaseAt(grain,el)) enddo - + thermal_conduction_getSpecificHeat = thermal_conduction_getSpecificHeat & / real(homogenization_Ngrains(material_homogenizationAt(el)),pReal) - + end function thermal_conduction_getSpecificHeat @@ -198,18 +195,18 @@ function thermal_conduction_getMassDensity(ip,el) thermal_conduction_getMassDensity integer :: & grain - + thermal_conduction_getMassDensity = 0.0_pReal - - + + do grain = 1, homogenization_Ngrains(material_homogenizationAt(el)) thermal_conduction_getMassDensity = thermal_conduction_getMassDensity & + lattice_massDensity(material_phaseAt(grain,el)) enddo - + thermal_conduction_getMassDensity = thermal_conduction_getMassDensity & / real(homogenization_Ngrains(material_homogenizationAt(el)),pReal) - + end function thermal_conduction_getMassDensity @@ -226,15 +223,15 @@ subroutine thermal_conduction_putTemperatureAndItsRate(T,Tdot,ip,el) Tdot integer :: & homog, & - offset - + offset + homog = material_homogenizationAt(el) offset = thermalMapping(homog)%p(ip,el) temperature (homog)%p(offset) = T temperatureRate(homog)%p(offset) = Tdot end subroutine thermal_conduction_putTemperatureAndItsRate - + !-------------------------------------------------------------------------------------------------- !> @brief writes results to HDF5 output file @@ -245,7 +242,7 @@ subroutine thermal_conduction_results(homog,group) character(len=*), intent(in) :: group integer :: o - + associate(prm => param(damage_typeInstance(homog))) outputsLoop: do o = 1,size(prm%output) select case(trim(prm%output(o))) diff --git a/src/thermal_isothermal.f90 b/src/thermal_isothermal.f90 index 3cafc6402..630e8db6d 100644 --- a/src/thermal_isothermal.f90 +++ b/src/thermal_isothermal.f90 @@ -7,11 +7,8 @@ module thermal_isothermal use material implicit none - private + public - public :: & - thermal_isothermal_init - contains !--------------------------------------------------------------------------------------------------