diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 2a141da56..ac41158a1 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -45,10 +45,10 @@ module homogenization materialpoint_stressAndItsTangent, & materialpoint_postResults private :: & - homogenization_partitionDeformation, & - homogenization_updateState, & - homogenization_averageStressAndItsTangent, & - homogenization_postResults + partitionDeformation, & + updateState, & + averageStressAndItsTangent, & + postResults contains @@ -118,12 +118,9 @@ subroutine homogenization_init !-------------------------------------------------------------------------------------------------- ! parse homogenization from config file - if (any(homogenization_type == HOMOGENIZATION_NONE_ID)) & - call homogenization_none_init() - if (any(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID)) & - call homogenization_isostrain_init(FILEUNIT) - if (any(homogenization_type == HOMOGENIZATION_RGC_ID)) & - call homogenization_RGC_init(FILEUNIT) + if (any(homogenization_type == HOMOGENIZATION_NONE_ID)) call homogenization_none_init + if (any(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID)) call homogenization_isostrain_init + if (any(homogenization_type == HOMOGENIZATION_RGC_ID)) call homogenization_RGC_init !-------------------------------------------------------------------------------------------------- ! parse thermal from config file @@ -156,17 +153,14 @@ subroutine homogenization_init select case(homogenization_type(p)) ! split per homogenization type case (HOMOGENIZATION_NONE_ID) outputName = HOMOGENIZATION_NONE_label - thisNoutput => null() thisOutput => null() thisSize => null() case (HOMOGENIZATION_ISOSTRAIN_ID) outputName = HOMOGENIZATION_ISOSTRAIN_label - thisNoutput => homogenization_isostrain_Noutput - thisOutput => homogenization_isostrain_output - thisSize => homogenization_isostrain_sizePostResult + thisOutput => null() + thisSize => null() case (HOMOGENIZATION_RGC_ID) outputName = HOMOGENIZATION_RGC_label - thisNoutput => homogenization_RGC_Noutput thisOutput => homogenization_RGC_output thisSize => homogenization_RGC_sizePostResult case default @@ -176,8 +170,9 @@ subroutine homogenization_init if (valid) then write(FILEUNIT,'(a)') '(type)'//char(9)//trim(outputName) write(FILEUNIT,'(a,i4)') '(ngrains)'//char(9),homogenization_Ngrains(p) - if (homogenization_type(p) /= HOMOGENIZATION_NONE_ID) then - do e = 1,thisNoutput(i) + if (homogenization_type(p) /= HOMOGENIZATION_NONE_ID .and. & + homogenization_type(p) /= HOMOGENIZATION_ISOSTRAIN_ID) then + do e = 1,size(thisOutput(:,i)) write(FILEUNIT,'(a,i4)') trim(thisOutput(e,i))//char(9),thisSize(e,i) enddo endif @@ -605,7 +600,7 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt) IpLooping2: do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) if ( materialpoint_requested(i,e) .and. & ! process requested but... .not. materialpoint_doneAndHappy(1,i,e)) then ! ...not yet done material points - call homogenization_partitionDeformation(i,e) ! partition deformation onto constituents + call partitionDeformation(i,e) ! partition deformation onto constituents crystallite_dt(1:myNgrains,i,e) = materialpoint_subdt(i,e) ! propagate materialpoint dt to grains crystallite_requested(1:myNgrains,i,e) = .true. ! request calculation for constituents else @@ -631,7 +626,7 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt) if (.not. materialpoint_converged(i,e)) then materialpoint_doneAndHappy(1:2,i,e) = [.true.,.false.] else - materialpoint_doneAndHappy(1:2,i,e) = homogenization_updateState(i,e) + materialpoint_doneAndHappy(1:2,i,e) = updateState(i,e) materialpoint_converged(i,e) = all(materialpoint_doneAndHappy(1:2,i,e)) ! converged if done and happy endif endif @@ -652,7 +647,7 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt) !$OMP PARALLEL DO elementLooping4: do e = FEsolving_execElem(1),FEsolving_execElem(2) IpLooping4: do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) - call homogenization_averageStressAndItsTangent(i,e) + call averageStressAndItsTangent(i,e) enddo IpLooping4 enddo elementLooping4 !$OMP END PARALLEL DO @@ -710,7 +705,7 @@ subroutine materialpoint_postResults thePos = thePos + 1_pInt if (theSize > 0_pInt) then ! any homogenization results to mention? - materialpoint_results(thePos+1:thePos+theSize,i,e) = homogenization_postResults(i,e) ! tell homogenization results + materialpoint_results(thePos+1:thePos+theSize,i,e) = postResults(i,e) ! tell homogenization results thePos = thePos + theSize endif @@ -734,12 +729,12 @@ end subroutine materialpoint_postResults !-------------------------------------------------------------------------------------------------- !> @brief partition material point def grad onto constituents !-------------------------------------------------------------------------------------------------- -subroutine homogenization_partitionDeformation(ip,el) +subroutine partitionDeformation(ip,el) use mesh, only: & mesh_element use material, only: & homogenization_type, & - homogenization_maxNgrains, & + homogenization_Ngrains, & HOMOGENIZATION_NONE_ID, & HOMOGENIZATION_ISOSTRAIN_ID, & HOMOGENIZATION_RGC_ID @@ -758,38 +753,36 @@ subroutine homogenization_partitionDeformation(ip,el) chosenHomogenization: select case(homogenization_type(mesh_element(3,el))) case (HOMOGENIZATION_NONE_ID) chosenHomogenization - crystallite_partionedF(1:3,1:3,1:homogenization_maxNgrains,ip,el) = 0.0_pReal - crystallite_partionedF(1:3,1:3,1:1,ip,el) = & - spread(materialpoint_subF(1:3,1:3,ip,el),3,1) + crystallite_partionedF(1:3,1:3,1,ip,el) = materialpoint_subF(1:3,1:3,ip,el) case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization call homogenization_isostrain_partitionDeformation(& - crystallite_partionedF(1:3,1:3,1:homogenization_maxNgrains,ip,el), & - materialpoint_subF(1:3,1:3,ip,el),& - el) + crystallite_partionedF(1:3,1:3,1:homogenization_Ngrains(mesh_element(3,el)),ip,el), & + materialpoint_subF(1:3,1:3,ip,el)) + case (HOMOGENIZATION_RGC_ID) chosenHomogenization call homogenization_RGC_partitionDeformation(& - crystallite_partionedF(1:3,1:3,1:homogenization_maxNgrains,ip,el), & + crystallite_partionedF(1:3,1:3,1:homogenization_Ngrains(mesh_element(3,el)),ip,el), & materialpoint_subF(1:3,1:3,ip,el),& ip, & el) end select chosenHomogenization -end subroutine homogenization_partitionDeformation +end subroutine partitionDeformation !-------------------------------------------------------------------------------------------------- !> @brief update the internal state of the homogenization scheme and tell whether "done" and !> "happy" with result !-------------------------------------------------------------------------------------------------- -function homogenization_updateState(ip,el) +function updateState(ip,el) use mesh, only: & mesh_element use material, only: & homogenization_type, & thermal_type, & damage_type, & - homogenization_maxNgrains, & + homogenization_Ngrains, & HOMOGENIZATION_RGC_ID, & THERMAL_adiabatic_ID, & DAMAGE_local_ID @@ -809,27 +802,27 @@ function homogenization_updateState(ip,el) integer(pInt), intent(in) :: & ip, & !< integration point el !< element number - logical, dimension(2) :: homogenization_updateState + logical, dimension(2) :: updateState - homogenization_updateState = .true. + updateState = .true. chosenHomogenization: select case(homogenization_type(mesh_element(3,el))) case (HOMOGENIZATION_RGC_ID) chosenHomogenization - homogenization_updateState = & - homogenization_updateState .and. & - homogenization_RGC_updateState(crystallite_P(1:3,1:3,1:homogenization_maxNgrains,ip,el), & - crystallite_partionedF(1:3,1:3,1:homogenization_maxNgrains,ip,el), & - crystallite_partionedF0(1:3,1:3,1:homogenization_maxNgrains,ip,el),& + updateState = & + updateState .and. & + homogenization_RGC_updateState(crystallite_P(1:3,1:3,1:homogenization_Ngrains(mesh_element(3,el)),ip,el), & + crystallite_partionedF(1:3,1:3,1:homogenization_Ngrains(mesh_element(3,el)),ip,el), & + crystallite_partionedF0(1:3,1:3,1:homogenization_Ngrains(mesh_element(3,el)),ip,el),& materialpoint_subF(1:3,1:3,ip,el),& materialpoint_subdt(ip,el), & - crystallite_dPdF(1:3,1:3,1:3,1:3,1:homogenization_maxNgrains,ip,el), & + crystallite_dPdF(1:3,1:3,1:3,1:3,1:homogenization_Ngrains(mesh_element(3,el)),ip,el), & ip, & el) end select chosenHomogenization chosenThermal: select case (thermal_type(mesh_element(3,el))) case (THERMAL_adiabatic_ID) chosenThermal - homogenization_updateState = & - homogenization_updateState .and. & + updateState = & + updateState .and. & thermal_adiabatic_updateState(materialpoint_subdt(ip,el), & ip, & el) @@ -837,25 +830,26 @@ function homogenization_updateState(ip,el) chosenDamage: select case (damage_type(mesh_element(3,el))) case (DAMAGE_local_ID) chosenDamage - homogenization_updateState = & - homogenization_updateState .and. & + updateState = & + updateState .and. & damage_local_updateState(materialpoint_subdt(ip,el), & ip, & el) end select chosenDamage -end function homogenization_updateState +end function updateState !-------------------------------------------------------------------------------------------------- !> @brief derive average stress and stiffness from constituent quantities !-------------------------------------------------------------------------------------------------- -subroutine homogenization_averageStressAndItsTangent(ip,el) +subroutine averageStressAndItsTangent(ip,el) use mesh, only: & mesh_element use material, only: & homogenization_type, & - homogenization_maxNgrains, & + homogenization_typeInstance, & + homogenization_Ngrains, & HOMOGENIZATION_NONE_ID, & HOMOGENIZATION_ISOSTRAIN_ID, & HOMOGENIZATION_RGC_ID @@ -873,36 +867,39 @@ subroutine homogenization_averageStressAndItsTangent(ip,el) chosenHomogenization: select case(homogenization_type(mesh_element(3,el))) case (HOMOGENIZATION_NONE_ID) chosenHomogenization - materialpoint_P(1:3,1:3,ip,el) = sum(crystallite_P(1:3,1:3,1:1,ip,el),3) - materialpoint_dPdF(1:3,1:3,1:3,1:3,ip,el) & - = sum(crystallite_dPdF(1:3,1:3,1:3,1:3,1:1,ip,el),5) + materialpoint_P(1:3,1:3,ip,el) = crystallite_P(1:3,1:3,1,ip,el) + materialpoint_dPdF(1:3,1:3,1:3,1:3,ip,el) = crystallite_dPdF(1:3,1:3,1:3,1:3,1,ip,el) case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization call homogenization_isostrain_averageStressAndItsTangent(& materialpoint_P(1:3,1:3,ip,el), & materialpoint_dPdF(1:3,1:3,1:3,1:3,ip,el),& - crystallite_P(1:3,1:3,1:homogenization_maxNgrains,ip,el), & - crystallite_dPdF(1:3,1:3,1:3,1:3,1:homogenization_maxNgrains,ip,el), & - el) + crystallite_P(1:3,1:3,1:homogenization_Ngrains(mesh_element(3,el)),ip,el), & + crystallite_dPdF(1:3,1:3,1:3,1:3,1:homogenization_Ngrains(mesh_element(3,el)),ip,el), & + homogenization_typeInstance(mesh_element(3,el))) + case (HOMOGENIZATION_RGC_ID) chosenHomogenization call homogenization_RGC_averageStressAndItsTangent(& materialpoint_P(1:3,1:3,ip,el), & materialpoint_dPdF(1:3,1:3,1:3,1:3,ip,el),& - crystallite_P(1:3,1:3,1:homogenization_maxNgrains,ip,el), & - crystallite_dPdF(1:3,1:3,1:3,1:3,1:homogenization_maxNgrains,ip,el), & - el) + crystallite_P(1:3,1:3,1:homogenization_Ngrains(mesh_element(3,el)),ip,el), & + crystallite_dPdF(1:3,1:3,1:3,1:3,1:homogenization_Ngrains(mesh_element(3,el)),ip,el), & + homogenization_typeInstance(mesh_element(3,el))) end select chosenHomogenization -end subroutine homogenization_averageStressAndItsTangent +end subroutine averageStressAndItsTangent + !-------------------------------------------------------------------------------------------------- !> @brief return array of homogenization results for post file inclusion. call only, !> if homogenization_sizePostResults(i,e) > 0 !! !-------------------------------------------------------------------------------------------------- -function homogenization_postResults(ip,el) +function postResults(ip,el) use mesh, only: & mesh_element use material, only: & + material_homogenizationAt, & + homogenization_typeInstance,& mappingHomogenization, & homogState, & thermalState, & @@ -919,8 +916,6 @@ function homogenization_postResults(ip,el) DAMAGE_none_ID, & DAMAGE_local_ID, & DAMAGE_nonlocal_ID - use homogenization_isostrain, only: & - homogenization_isostrain_postResults use homogenization_RGC, only: & homogenization_RGC_postResults use thermal_adiabatic, only: & @@ -939,60 +934,46 @@ function homogenization_postResults(ip,el) real(pReal), dimension( homogState (mappingHomogenization(2,ip,el))%sizePostResults & + thermalState (mappingHomogenization(2,ip,el))%sizePostResults & + damageState (mappingHomogenization(2,ip,el))%sizePostResults) :: & - homogenization_postResults + postResults integer(pInt) :: & - startPos, endPos + startPos, endPos ,& + of, instance - homogenization_postResults = 0.0_pReal + postResults = 0.0_pReal startPos = 1_pInt endPos = homogState(mappingHomogenization(2,ip,el))%sizePostResults chosenHomogenization: select case (homogenization_type(mesh_element(3,el))) - case (HOMOGENIZATION_NONE_ID) chosenHomogenization - case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization - homogenization_postResults(startPos:endPos) = & - homogenization_isostrain_postResults(& - ip, & - el, & - materialpoint_P(1:3,1:3,ip,el), & - materialpoint_F(1:3,1:3,ip,el)) case (HOMOGENIZATION_RGC_ID) chosenHomogenization - homogenization_postResults(startPos:endPos) = & - homogenization_RGC_postResults(& - ip, & - el, & - materialpoint_P(1:3,1:3,ip,el), & - materialpoint_F(1:3,1:3,ip,el)) + instance = homogenization_typeInstance(material_homogenizationAt(el)) + of = mappingHomogenization(1,ip,el) + postResults(startPos:endPos) = homogenization_RGC_postResults(instance,of) + end select chosenHomogenization startPos = endPos + 1_pInt endPos = endPos + thermalState(mappingHomogenization(2,ip,el))%sizePostResults chosenThermal: select case (thermal_type(mesh_element(3,el))) - case (THERMAL_isothermal_ID) chosenThermal case (THERMAL_adiabatic_ID) chosenThermal - homogenization_postResults(startPos:endPos) = & - thermal_adiabatic_postResults(ip, el) + postResults(startPos:endPos) = thermal_adiabatic_postResults(ip, el) case (THERMAL_conduction_ID) chosenThermal - homogenization_postResults(startPos:endPos) = & - thermal_conduction_postResults(ip, el) + postResults(startPos:endPos) = thermal_conduction_postResults(ip, el) + end select chosenThermal startPos = endPos + 1_pInt endPos = endPos + damageState(mappingHomogenization(2,ip,el))%sizePostResults chosenDamage: select case (damage_type(mesh_element(3,el))) - case (DAMAGE_none_ID) chosenDamage case (DAMAGE_local_ID) chosenDamage - homogenization_postResults(startPos:endPos) = & - damage_local_postResults(ip, el) - + postResults(startPos:endPos) = damage_local_postResults(ip, el) case (DAMAGE_nonlocal_ID) chosenDamage - homogenization_postResults(startPos:endPos) = & - damage_nonlocal_postResults(ip, el) + postResults(startPos:endPos) = damage_nonlocal_postResults(ip, el) + end select chosenDamage -end function homogenization_postResults +end function postResults end module homogenization diff --git a/src/homogenization_RGC.f90 b/src/homogenization_RGC.f90 index 1d7bc6f86..ef81043eb 100644 --- a/src/homogenization_RGC.f90 +++ b/src/homogenization_RGC.f90 @@ -1,9 +1,10 @@ !-------------------------------------------------------------------------------------------------- +!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @author Denny Tjahjanto, Max-Planck-Institut für Eisenforschung GmbH !> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH !> @brief Relaxed grain cluster (RGC) homogenization scheme -!> Ngrains is defined as p x q x r (cluster) +!> Nconstituents is defined as p x q x r (cluster) !-------------------------------------------------------------------------------------------------- module homogenization_RGC use prec, only: & @@ -12,39 +13,63 @@ module homogenization_RGC implicit none private - integer(pInt), dimension(:), allocatable, public :: & - homogenization_RGC_sizeState, & - homogenization_RGC_sizePostResults integer(pInt), dimension(:,:), allocatable,target, public :: & homogenization_RGC_sizePostResult character(len=64), dimension(:,:), allocatable,target, public :: & homogenization_RGC_output ! name of each post result output - integer(pInt), dimension(:), allocatable,target, public :: & - homogenization_RGC_Noutput !< number of outputs per homog instance - integer(pInt), dimension(:,:), allocatable, private :: & - homogenization_RGC_Ngrains - real(pReal), dimension(:,:), allocatable, private :: & - homogenization_RGC_dAlpha, & - homogenization_RGC_angles - real(pReal), dimension(:,:,:,:), allocatable, private :: & - homogenization_RGC_orientation - real(pReal), dimension(:), allocatable, private :: & - homogenization_RGC_xiAlpha, & - homogenization_RGC_ciAlpha + enum, bind(c) - enumerator :: undefined_ID, & - constitutivework_ID, & - penaltyenergy_ID, & - volumediscrepancy_ID, & - averagerelaxrate_ID,& - maximumrelaxrate_ID,& - ipcoords_ID,& - magnitudemismatch_ID,& - avgdefgrad_ID,& - avgfirstpiola_ID + enumerator :: & + undefined_ID, & + constitutivework_ID, & + penaltyenergy_ID, & + volumediscrepancy_ID, & + averagerelaxrate_ID,& + maximumrelaxrate_ID,& + magnitudemismatch_ID end enum - integer(kind(undefined_ID)), dimension(:,:), allocatable, private :: & - homogenization_RGC_outputID !< ID of each post result output + + type, private :: tParameters + integer(pInt), dimension(:), allocatable :: & + Nconstituents + real(pReal) :: & + xiAlpha, & + ciAlpha + real(pReal), dimension(:), allocatable :: & + dAlpha, & + angles + integer(pInt) :: & + of_debug = 0_pInt + integer(kind(undefined_ID)), dimension(:), allocatable :: & + outputID + end type + + type, private :: tRGCstate + real(pReal), pointer, dimension(:) :: & + work, & + penaltyEnergy + real(pReal), pointer, dimension(:,:) :: & + relaxationVector + end type tRGCstate + + type, private :: tRGCdependentState + real(pReal), allocatable, dimension(:) :: & + volumeDiscrepancy, & + relaxationRate_avg, & + relaxationRate_max + real(pReal), allocatable, dimension(:,:) :: & + mismatch + real(pReal), allocatable, dimension(:,:,:) :: & + orientation + end type tRGCdependentState + + type(tparameters), dimension(:), allocatable, private :: & + param + type(tRGCstate), dimension(:), allocatable, private :: & + state, & + state0 + type(tRGCdependentState), dimension(:), allocatable, private :: & + dependentState public :: & homogenization_RGC_init, & @@ -53,313 +78,239 @@ module homogenization_RGC homogenization_RGC_updateState, & homogenization_RGC_postResults private :: & - homogenization_RGC_stressPenalty, & - homogenization_RGC_volumePenalty, & - homogenization_RGC_grainDeformation, & - homogenization_RGC_surfaceCorrection, & - homogenization_RGC_equivalentModuli, & - homogenization_RGC_relaxationVector, & - homogenization_RGC_interfaceNormal, & - homogenization_RGC_getInterface, & - homogenization_RGC_grain1to3, & - homogenization_RGC_grain3to1, & - homogenization_RGC_interface4to1, & - homogenization_RGC_interface1to4 + relaxationVector, & + interfaceNormal, & + getInterface, & + grain1to3, & + grain3to1, & + interface4to1, & + interface1to4 contains !-------------------------------------------------------------------------------------------------- !> @brief allocates all necessary fields, reads information from material configuration file !-------------------------------------------------------------------------------------------------- -subroutine homogenization_RGC_init(fileUnit) +subroutine homogenization_RGC_init() #if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & compiler_options #endif - use prec, only: & - pReal, & - pInt use debug, only: & - debug_level, & - debug_homogenization, & - debug_levelBasic, & - debug_levelExtensive +#ifdef DEBUG + debug_i, & + debug_e, & +#endif + debug_level, & + debug_homogenization, & + debug_levelBasic use math, only: & - math_Mandel3333to66,& - math_Voigt66to3333, & - math_I3, & - math_sampleRandomOri,& - math_EulerToR,& + math_EulerToR, & INRAD - use mesh, only: & - mesh_maxNips, & - mesh_NcpElems,& - mesh_element, & - FE_Nips, & - FE_geomtype - use IO - use material - use config + use IO, only: & + IO_error, & + IO_timeStamp + use material, only: & +#ifdef DEBUG + material_homogenizationAt, & + mappingHomogenization, & +#endif + homogenization_type, & + material_homog, & + homogState, & + HOMOGENIZATION_RGC_ID, & + HOMOGENIZATION_RGC_LABEL, & + homogenization_typeInstance, & + homogenization_Noutput, & + homogenization_Ngrains + use config, only: & + config_homogenization implicit none - integer(pInt), intent(in) :: fileUnit !< file pointer to material configuration - integer(pInt), allocatable, dimension(:) :: chunkPos - integer :: & - homog, & - NofMyHomog, & - o, & - instance, & - sizeHState - integer(pInt) :: section=0_pInt, maxNinstance, i,j,e, mySize, myInstance - character(len=65536) :: & - tag = '', & - line = '' + integer(pInt) :: & + Ninstance, & + h, i, & + NofMyHomog, outputSize, & + sizeState, nIntFaceTot + + character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::] + + integer(kind(undefined_ID)) :: & + outputID + + character(len=65536), dimension(:), allocatable :: & + outputs write(6,'(/,a)') ' <<<+- homogenization_'//HOMOGENIZATION_RGC_label//' init -+>>>' write(6,'(/,a)') ' Tjahjanto et al., International Journal of Material Forming, 2(1):939–942, 2009' - write(6,'(/,a)') ' https://doi.org/10.1007/s12289-009-0619-1' + write(6,'(a)') ' https://doi.org/10.1007/s12289-009-0619-1' write(6,'(/,a)') ' Tjahjanto et al., Modelling and Simulation in Materials Science and Engineering, 18:015006, 2010' - write(6,'(/,a)') ' https://doi.org/10.1088/0965-0393/18/1/015006' + write(6,'(a)') ' https://doi.org/10.1088/0965-0393/18/1/015006' write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" - maxNinstance = int(count(homogenization_type == HOMOGENIZATION_RGC_ID),pInt) - if (maxNinstance == 0_pInt) return - if (iand(debug_level(debug_HOMOGENIZATION),debug_levelBasic) /= 0_pInt) & - write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance - allocate(homogenization_RGC_sizeState(maxNinstance), source=0_pInt) - allocate(homogenization_RGC_sizePostResults(maxNinstance), source=0_pInt) - allocate(homogenization_RGC_Noutput(maxNinstance), source=0_pInt) - allocate(homogenization_RGC_Ngrains(3,maxNinstance), source=0_pInt) - allocate(homogenization_RGC_ciAlpha(maxNinstance), source=0.0_pReal) - allocate(homogenization_RGC_xiAlpha(maxNinstance), source=0.0_pReal) - allocate(homogenization_RGC_dAlpha(3,maxNinstance), source=0.0_pReal) - allocate(homogenization_RGC_angles(3,maxNinstance), source=400.0_pReal) - allocate(homogenization_RGC_output(maxval(homogenization_Noutput),maxNinstance)) - homogenization_RGC_output='' - allocate(homogenization_RGC_outputID(maxval(homogenization_Noutput),maxNinstance),source=undefined_ID) - allocate(homogenization_RGC_sizePostResult(maxval(homogenization_Noutput),maxNinstance),& - source=0_pInt) - allocate(homogenization_RGC_orientation(3,3,mesh_maxNips,mesh_NcpElems), source=0.0_pReal) - homogenization_RGC_orientation = spread(spread(math_I3,3,mesh_maxNips),4,mesh_NcpElems) ! initialize to identity + Ninstance = int(count(homogenization_type == HOMOGENIZATION_RGC_ID),pInt) + if (iand(debug_level(debug_HOMOGENIZATION),debug_levelBasic) /= 0_pInt) & + write(6,'(a16,1x,i5,/)') '# instances:',Ninstance + + allocate(param(Ninstance)) + allocate(state(Ninstance)) + allocate(state0(Ninstance)) + allocate(dependentState(Ninstance)) - rewind(fileUnit) - do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>'))/=material_partHomogenization) ! wind forward to - line = IO_read(fileUnit) - enddo + allocate(homogenization_RGC_sizePostResult(maxval(homogenization_Noutput),Ninstance),source=0_pInt) + allocate(homogenization_RGC_output(maxval(homogenization_Noutput),Ninstance)) + homogenization_RGC_output='' - parsingFile: do while (trim(line) /= IO_EOF) ! read through sections of homogenization part - line = IO_read(fileUnit) - if (IO_isBlank(line)) cycle ! skip empty lines - if (IO_getTag(line,'<','>') /= '') then ! stop at next part - line = IO_read(fileUnit, .true.) ! reset IO_read - exit + do h = 1_pInt, size(homogenization_type) + if (homogenization_type(h) /= HOMOGENIZATION_RGC_ID) cycle + associate(prm => param(homogenization_typeInstance(h)), & + stt => state(homogenization_typeInstance(h)), & + st0 => state0(homogenization_typeInstance(h)), & + dst => dependentState(homogenization_typeInstance(h)), & + config => config_homogenization(h)) + +#ifdef DEBUG + if (h==material_homogenizationAt(debug_e)) then + prm%of_debug = mappingHomogenization(1,debug_i,debug_e) endif - if (IO_getTag(line,'[',']') /= '') then ! next section - section = section + 1_pInt - cycle - endif - if (section > 0_pInt ) then ! do not short-circuit here (.and. with next if-statement). It's not safe in Fortran - if (homogenization_type(section) == HOMOGENIZATION_RGC_ID) then ! one of my sections - i = homogenization_typeInstance(section) ! which instance of my type is present homogenization - chunkPos = IO_stringPos(line) - tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key - select case(tag) - case ('(output)') - homogenization_RGC_Noutput(i) = homogenization_RGC_Noutput(i) + 1_pInt - homogenization_RGC_output(homogenization_RGC_Noutput(i),i) = & - IO_lc(IO_stringValue(line,chunkPos,2_pInt)) - select case(IO_lc(IO_stringValue(line,chunkPos,2_pInt))) - case('constitutivework') - homogenization_RGC_outputID(homogenization_RGC_Noutput(i),i) = constitutivework_ID - case('penaltyenergy') - homogenization_RGC_outputID(homogenization_RGC_Noutput(i),i) = penaltyenergy_ID - case('volumediscrepancy') - homogenization_RGC_outputID(homogenization_RGC_Noutput(i),i) = volumediscrepancy_ID - case('averagerelaxrate') - homogenization_RGC_outputID(homogenization_RGC_Noutput(i),i) = averagerelaxrate_ID - case('maximumrelaxrate') - homogenization_RGC_outputID(homogenization_RGC_Noutput(i),i) = maximumrelaxrate_ID - case('magnitudemismatch') - homogenization_RGC_outputID(homogenization_RGC_Noutput(i),i) = magnitudemismatch_ID - case('ipcoords') - homogenization_RGC_outputID(homogenization_RGC_Noutput(i),i) = ipcoords_ID - case('avgdefgrad','avgf') - homogenization_RGC_outputID(homogenization_RGC_Noutput(i),i) = avgdefgrad_ID - case('avgp','avgfirstpiola','avg1stpiola') - homogenization_RGC_outputID(homogenization_RGC_Noutput(i),i) = avgfirstpiola_ID - case default - homogenization_RGC_Noutput(i) = homogenization_RGC_Noutput(i) -1_pInt ! correct for invalid +#endif - end select - case ('clustersize') - homogenization_RGC_Ngrains(1,i) = IO_intValue(line,chunkPos,2_pInt) - homogenization_RGC_Ngrains(2,i) = IO_intValue(line,chunkPos,3_pInt) - homogenization_RGC_Ngrains(3,i) = IO_intValue(line,chunkPos,4_pInt) - if (homogenization_Ngrains(section) /= product(homogenization_RGC_Ngrains(1:3,i))) & - call IO_error(211_pInt,ext_msg=trim(tag)//' ('//HOMOGENIZATION_RGC_label//')') - case ('scalingparameter') - homogenization_RGC_xiAlpha(i) = IO_floatValue(line,chunkPos,2_pInt) - case ('overproportionality') - homogenization_RGC_ciAlpha(i) = IO_floatValue(line,chunkPos,2_pInt) - case ('grainsize') - homogenization_RGC_dAlpha(1,i) = IO_floatValue(line,chunkPos,2_pInt) - homogenization_RGC_dAlpha(2,i) = IO_floatValue(line,chunkPos,3_pInt) - homogenization_RGC_dAlpha(3,i) = IO_floatValue(line,chunkPos,4_pInt) - case ('clusterorientation') - homogenization_RGC_angles(1,i) = IO_floatValue(line,chunkPos,2_pInt) - homogenization_RGC_angles(2,i) = IO_floatValue(line,chunkPos,3_pInt) - homogenization_RGC_angles(3,i) = IO_floatValue(line,chunkPos,4_pInt) + prm%Nconstituents = config%getInts('clustersize',requiredShape=[3]) + if (homogenization_Ngrains(h) /= product(prm%Nconstituents)) & + call IO_error(211_pInt,ext_msg='clustersize ('//HOMOGENIZATION_RGC_label//')') - end select - endif - endif - enddo parsingFile + prm%xiAlpha = config%getFloat('scalingparameter') + prm%ciAlpha = config%getFloat('overproportionality') -!-------------------------------------------------------------------------------------------------- -! * assigning cluster orientations - elementLooping: do e = 1_pInt,mesh_NcpElems - if (homogenization_type(mesh_element(3,e)) == HOMOGENIZATION_RGC_ID) then - myInstance = homogenization_typeInstance(mesh_element(3,e)) - if (all (homogenization_RGC_angles(1:3,myInstance) >= 399.9_pReal)) then - homogenization_RGC_orientation(1:3,1:3,1,e) = math_EulerToR(math_sampleRandomOri()) - do i = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,e))) - if (microstructure_elemhomo(mesh_element(4,e))) then - homogenization_RGC_orientation(1:3,1:3,i,e) = homogenization_RGC_orientation(1:3,1:3,1,e) - else - homogenization_RGC_orientation(1:3,1:3,i,e) = math_EulerToR(math_sampleRandomOri()) - endif - enddo - else - do i = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,e))) - homogenization_RGC_orientation(1:3,1:3,i,e) = & - math_EulerToR(homogenization_RGC_angles(1:3,myInstance)*inRad) - enddo - endif - endif - enddo elementLooping + prm%dAlpha = config%getFloats('grainsize', requiredShape=[3]) + prm%angles = config%getFloats('clusterorientation',requiredShape=[3]) - if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt) then - do i = 1_pInt,maxNinstance - write(6,'(a15,1x,i4,/)') 'instance: ', i - write(6,'(a25,3(1x,i8))') 'cluster size: ',(homogenization_RGC_Ngrains(j,i),j=1_pInt,3_pInt) - write(6,'(a25,1x,e10.3)') 'scaling parameter: ', homogenization_RGC_xiAlpha(i) - write(6,'(a25,1x,e10.3)') 'over-proportionality: ', homogenization_RGC_ciAlpha(i) - write(6,'(a25,3(1x,e10.3))') 'grain size: ',(homogenization_RGC_dAlpha(j,i),j=1_pInt,3_pInt) - write(6,'(a25,3(1x,e10.3))') 'cluster orientation: ',(homogenization_RGC_angles(j,i),j=1_pInt,3_pInt) - enddo - endif -!-------------------------------------------------------------------------------------------------- - initializeInstances: do homog = 1_pInt, material_Nhomogenization - myHomog: if (homogenization_type(homog) == HOMOGENIZATION_RGC_ID) then - NofMyHomog = count(material_homog == homog) - instance = homogenization_typeInstance(homog) + outputs = config%getStrings('(output)',defaultVal=emptyStringArray) + allocate(prm%outputID(0)) -! * Determine size of postResults array - outputsLoop: do o = 1_pInt, homogenization_RGC_Noutput(instance) - select case(homogenization_RGC_outputID(o,instance)) - case(constitutivework_ID,penaltyenergy_ID,volumediscrepancy_ID, & - averagerelaxrate_ID,maximumrelaxrate_ID) - mySize = 1_pInt - case(ipcoords_ID,magnitudemismatch_ID) - mySize = 3_pInt - case(avgdefgrad_ID,avgfirstpiola_ID) - mySize = 9_pInt - case default - mySize = 0_pInt - end select - - outputFound: if (mySize > 0_pInt) then - homogenization_RGC_sizePostResult(o,instance) = mySize - homogenization_RGC_sizePostResults(instance) = & - homogenization_RGC_sizePostResults(instance) + mySize - endif outputFound - enddo outputsLoop + do i=1_pInt, size(outputs) + outputID = undefined_ID + select case(outputs(i)) - sizeHState = & - 3_pInt*(homogenization_RGC_Ngrains(1,instance)-1_pInt)* & - homogenization_RGC_Ngrains(2,instance)*homogenization_RGC_Ngrains(3,instance) & - + 3_pInt*homogenization_RGC_Ngrains(1,instance)*(homogenization_RGC_Ngrains(2,instance)-1_pInt)* & - homogenization_RGC_Ngrains(3,instance) & - + 3_pInt*homogenization_RGC_Ngrains(1,instance)*homogenization_RGC_Ngrains(2,instance)* & - (homogenization_RGC_Ngrains(3,instance)-1_pInt) & - + 8_pInt ! (1) Average constitutive work, (2-4) Overall mismatch, (5) Average penalty energy, - ! (6) Volume discrepancy, (7) Avg relaxation rate component, (8) Max relaxation rate component - -! allocate state arrays - homogState(homog)%sizeState = sizeHState - homogState(homog)%sizePostResults = homogenization_RGC_sizePostResults(instance) - allocate(homogState(homog)%state0 (sizeHState,NofMyHomog), source=0.0_pReal) - allocate(homogState(homog)%subState0(sizeHState,NofMyHomog), source=0.0_pReal) - allocate(homogState(homog)%state (sizeHState,NofMyHomog), source=0.0_pReal) + case('constitutivework') + outputID = constitutivework_ID + outputSize = 1_pInt + case('penaltyenergy') + outputID = penaltyenergy_ID + outputSize = 1_pInt + case('volumediscrepancy') + outputID = volumediscrepancy_ID + outputSize = 1_pInt + case('averagerelaxrate') + outputID = averagerelaxrate_ID + outputSize = 1_pInt + case('maximumrelaxrate') + outputID = maximumrelaxrate_ID + outputSize = 1_pInt + case('magnitudemismatch') + outputID = magnitudemismatch_ID + outputSize = 3_pInt - endif myHomog - enddo initializeInstances + end select + + if (outputID /= undefined_ID) then + homogenization_RGC_output(i,homogenization_typeInstance(h)) = outputs(i) + homogenization_RGC_sizePostResult(i,homogenization_typeInstance(h)) = outputSize + prm%outputID = [prm%outputID , outputID] + endif + + enddo + + NofMyHomog = count(material_homog == h) + nIntFaceTot = 3_pInt*( (prm%Nconstituents(1)-1_pInt)*prm%Nconstituents(2)*prm%Nconstituents(3) & + + prm%Nconstituents(1)*(prm%Nconstituents(2)-1_pInt)*prm%Nconstituents(3) & + + prm%Nconstituents(1)*prm%Nconstituents(2)*(prm%Nconstituents(3)-1_pInt)) + sizeState = nIntFaceTot & + + size(['avg constitutive work ','average penalty energy']) + + homogState(h)%sizeState = sizeState + homogState(h)%sizePostResults = sum(homogenization_RGC_sizePostResult(:,homogenization_typeInstance(h))) + allocate(homogState(h)%state0 (sizeState,NofMyHomog), source=0.0_pReal) + allocate(homogState(h)%subState0(sizeState,NofMyHomog), source=0.0_pReal) + allocate(homogState(h)%state (sizeState,NofMyHomog), source=0.0_pReal) + + stt%relaxationVector => homogState(h)%state(1:nIntFaceTot,:) + st0%relaxationVector => homogState(h)%state0(1:nIntFaceTot,:) + stt%work => homogState(h)%state(nIntFaceTot+1,:) + stt%penaltyEnergy => homogState(h)%state(nIntFaceTot+2,:) + + allocate(dst%volumeDiscrepancy( NofMyHomog)) + allocate(dst%relaxationRate_avg( NofMyHomog)) + allocate(dst%relaxationRate_max( NofMyHomog)) + allocate(dst%mismatch( 3,NofMyHomog)) + +!-------------------------------------------------------------------------------------------------- +! assigning cluster orientations + dependentState(homogenization_typeInstance(h))%orientation = spread(math_EulerToR(prm%angles*inRad),3,NofMyHomog) + !dst%orientation = spread(math_EulerToR(prm%angles*inRad),3,NofMyHomog) ifort version 18.0.1 crashes (for whatever reason) + + end associate + + enddo - - end subroutine homogenization_RGC_init !-------------------------------------------------------------------------------------------------- !> @brief partitions the deformation gradient onto the constituents !-------------------------------------------------------------------------------------------------- -subroutine homogenization_RGC_partitionDeformation(F,avgF,ip,el) +subroutine homogenization_RGC_partitionDeformation(F,avgF,instance,of) +#ifdef DEBUG use debug, only: & debug_level, & debug_homogenization, & debug_levelExtensive - use mesh, only: & - mesh_element - use material, only: & - homogenization_maxNgrains, & - homogenization_Ngrains,& - homogenization_typeInstance +#endif implicit none - real(pReal), dimension (3,3,homogenization_maxNgrains), intent(out) :: F !< partioned F per grain - real(pReal), dimension (3,3), intent(in) :: avgF !< averaged F - integer(pInt), intent(in) :: & - ip, & !< integration point number - el !< element number + real(pReal), dimension (:,:,:), intent(out) :: F !< partioned F per grain + + real(pReal), dimension (:,:), intent(in) :: avgF !< averaged F + integer(pInt), intent(in) :: & + instance, & + of + real(pReal), dimension (3) :: aVect,nVect integer(pInt), dimension (4) :: intFace integer(pInt), dimension (3) :: iGrain3 - integer(pInt) :: homID, iGrain,iFace,i,j - integer(pInt), parameter :: nFace = 6_pInt + integer(pInt) :: iGrain,iFace,i,j + associate(prm => param(instance)) + !-------------------------------------------------------------------------------------------------- ! compute the deformation gradient of individual grains due to relaxations - homID = homogenization_typeInstance(mesh_element(3,el)) F = 0.0_pReal - do iGrain = 1_pInt,homogenization_Ngrains(mesh_element(3,el)) - iGrain3 = homogenization_RGC_grain1to3(iGrain,homID) - do iFace = 1_pInt,nFace - intFace = homogenization_RGC_getInterface(iFace,iGrain3) ! identifying 6 interfaces of each grain - - aVect = homogenization_RGC_relaxationVector(intFace,homID, ip, el) ! get the relaxation vectors for each interface from global relaxation vector array - - nVect = homogenization_RGC_interfaceNormal(intFace,ip,el) ! get the normal of each interface + do iGrain = 1_pInt,product(prm%Nconstituents) + iGrain3 = grain1to3(iGrain,prm%Nconstituents) + do iFace = 1_pInt,6_pInt + intFace = getInterface(iFace,iGrain3) ! identifying 6 interfaces of each grain + aVect = relaxationVector(intFace,instance,of) ! get the relaxation vectors for each interface from global relaxation vector array + nVect = interfaceNormal(intFace,instance,of) forall (i=1_pInt:3_pInt,j=1_pInt:3_pInt) & - F(i,j,iGrain) = F(i,j,iGrain) + aVect(i)*nVect(j) ! calculating deformation relaxations due to interface relaxation + F(i,j,iGrain) = F(i,j,iGrain) + aVect(i)*nVect(j) ! calculating deformation relaxations due to interface relaxation enddo F(1:3,1:3,iGrain) = F(1:3,1:3,iGrain) + avgF ! resulting relaxed deformation gradient -!-------------------------------------------------------------------------------------------------- -! debugging the grain deformation gradients +#ifdef DEBUG if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt) then - !$OMP CRITICAL (write2out) write(6,'(1x,a32,1x,i3)')'Deformation gradient of grain: ',iGrain do i = 1_pInt,3_pInt write(6,'(1x,3(e15.8,1x))')(F(i,j,iGrain), j = 1_pInt,3_pInt) enddo write(6,*)' ' flush(6) - !$OMP END CRITICAL (write2out) endif - +#endif enddo + + end associate end subroutine homogenization_RGC_partitionDeformation @@ -371,22 +322,18 @@ end subroutine homogenization_RGC_partitionDeformation function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) use prec, only: & dEq0 +#ifdef DEBUG use debug, only: & debug_level, & debug_homogenization,& - debug_levelExtensive, & - debug_e, & - debug_i + debug_levelExtensive +#endif use math, only: & - math_invert - use mesh, only: & - mesh_element + math_invert2 use material, only: & - homogenization_maxNgrains, & + material_homogenizationAt, & homogenization_typeInstance, & - homogState, & - mappingHomogenization, & - homogenization_Ngrains + mappingHomogenization use numerics, only: & absTol_RGC, & relTol_RGC, & @@ -400,112 +347,112 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) implicit none - real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: & + real(pReal), dimension (:,:,:), intent(in) :: & P,& !< array of P F,& !< array of F F0 !< array of initial F - real(pReal), dimension (3,3,3,3,homogenization_maxNgrains), intent(in) :: dPdF !< array of current grain stiffness - real(pReal), dimension (3,3), intent(in) :: avgF !< average F - real(pReal), intent(in) :: dt !< time increment - integer(pInt), intent(in) :: & + real(pReal), dimension (:,:,:,:,:), intent(in) :: dPdF !< array of current grain stiffness + real(pReal), dimension (3,3), intent(in) :: avgF !< average F + real(pReal), intent(in) :: dt !< time increment + integer(pInt), intent(in) :: & ip, & !< integration point number el !< element number + logical, dimension(2) :: homogenization_RGC_updateState + integer(pInt), dimension (4) :: intFaceN,intFaceP,faceID - integer(pInt), dimension (3) :: nGDim,iGr3N,iGr3P,stresLoc - integer(pInt), dimension (2) :: residLoc - integer(pInt) homID,iNum,i,j,nIntFaceTot,iGrN,iGrP,iMun,iFace,k,l,ipert,iGrain,nGrain - real(pReal), dimension (3,3,homogenization_maxNgrains) :: R,pF,pR,D,pD - real(pReal), dimension (3,homogenization_maxNgrains) :: NN,pNN - real(pReal), dimension (3) :: normP,normN,mornP,mornN - real(pReal) :: residMax,stresMax,constitutiveWork,penaltyEnergy,volDiscrep - logical error - - integer(pInt), parameter :: nFace = 6_pInt - + integer(pInt), dimension (3) :: nGDim,iGr3N,iGr3P + integer(pInt) :: instance,iNum,i,j,nIntFaceTot,iGrN,iGrP,iMun,iFace,k,l,ipert,iGrain,nGrain, of + real(pReal), dimension (3,3,size(P,3)) :: R,pF,pR,D,pD + real(pReal), dimension (3,size(P,3)) :: NN,devNull + real(pReal), dimension (3) :: normP,normN,mornP,mornN + real(pReal) :: residMax,stresMax + logical :: error real(pReal), dimension(:,:), allocatable :: tract,jmatrix,jnverse,smatrix,pmatrix,rmatrix real(pReal), dimension(:), allocatable :: resid,relax,p_relax,p_resid,drelax +#ifdef DEBUG + integer(pInt), dimension (3) :: stresLoc + integer(pInt), dimension (2) :: residLoc +#endif zeroTimeStep: if(dEq0(dt)) then homogenization_RGC_updateState = .true. ! pretend everything is fine and return return endif zeroTimeStep + instance = homogenization_typeInstance(material_homogenizationAt(el)) + of = mappingHomogenization(1,ip,el) + + associate(stt => state(instance), st0 => state0(instance), dst => dependentState(instance), prm => param(instance)) + !-------------------------------------------------------------------------------------------------- ! get the dimension of the cluster (grains and interfaces) - homID = homogenization_typeInstance(mesh_element(3,el)) - nGDim = homogenization_RGC_Ngrains(1:3,homID) - nGrain = homogenization_Ngrains(mesh_element(3,el)) - nIntFaceTot = (nGDim(1)-1_pInt)*nGDim(2)*nGDim(3) + nGDim(1)*(nGDim(2)-1_pInt)*nGDim(3) & - + nGDim(1)*nGDim(2)*(nGDim(3)-1_pInt) + nGDim = prm%Nconstituents + nGrain = product(nGDim) + nIntFaceTot = (nGDim(1)-1_pInt)*nGDim(2)*nGDim(3) & + + nGDim(1)*(nGDim(2)-1_pInt)*nGDim(3) & + + nGDim(1)*nGDim(2)*(nGDim(3)-1_pInt) !-------------------------------------------------------------------------------------------------- ! allocate the size of the global relaxation arrays/jacobian matrices depending on the size of the cluster allocate(resid(3_pInt*nIntFaceTot), source=0.0_pReal) allocate(tract(nIntFaceTot,3), source=0.0_pReal) - allocate(relax(3_pInt*nIntFaceTot)); relax= homogState(mappingHomogenization(2,ip,el))% & - state(1:3_pInt*nIntFaceTot,mappingHomogenization(1,ip,el)) - allocate(drelax(3_pInt*nIntFaceTot)); drelax= homogState(mappingHomogenization(2,ip,el))% & - state(1:3_pInt*nIntFaceTot,mappingHomogenization(1,ip,el)) - & - homogState(mappingHomogenization(2,ip,el))% & - state0(1:3_pInt*nIntFaceTot,mappingHomogenization(1,ip,el)) -!-------------------------------------------------------------------------------------------------- -! debugging the obtained state + relax = stt%relaxationVector(:,of) + drelax = stt%relaxationVector(:,of) - st0%relaxationVector(:,of) + +#ifdef DEBUG if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt) then - !$OMP CRITICAL (write2out) write(6,'(1x,a30)')'Obtained state: ' - do i = 1_pInt,3_pInt*nIntFaceTot - write(6,'(1x,2(e15.8,1x))')homogState(mappingHomogenization(2,ip,el))%state(i,mappingHomogenization(1,ip,el)) + do i = 1_pInt,size(stt%relaxationVector(:,of)) + write(6,'(1x,2(e15.8,1x))') stt%relaxationVector(i,of) enddo write(6,*)' ' - !$OMP END CRITICAL (write2out) endif +#endif !-------------------------------------------------------------------------------------------------- ! computing interface mismatch and stress penalty tensor for all interfaces of all grains - call homogenization_RGC_stressPenalty(R,NN,avgF,F,ip,el,homID) + call stressPenalty(R,NN,avgF,F,ip,el,instance,of) !-------------------------------------------------------------------------------------------------- ! calculating volume discrepancy and stress penalty related to overall volume discrepancy - call homogenization_RGC_volumePenalty(D,volDiscrep,F,avgF,ip,el) + call volumePenalty(D,dst%volumeDiscrepancy(of),avgF,F,nGrain,instance,of) -!-------------------------------------------------------------------------------------------------- -! debugging the mismatch, stress and penalties of grains +#ifdef DEBUG if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt) then - !$OMP CRITICAL (write2out) do iGrain = 1_pInt,nGrain write(6,'(1x,a30,1x,i3,1x,a4,3(1x,e15.8))')'Mismatch magnitude of grain(',iGrain,') :',& NN(1,iGrain),NN(2,iGrain),NN(3,iGrain) write(6,'(/,1x,a30,1x,i3)')'Stress and penalties of grain: ',iGrain do i = 1_pInt,3_pInt write(6,'(1x,3(e15.8,1x),1x,3(e15.8,1x),1x,3(e15.8,1x))')(P(i,j,iGrain), j = 1_pInt,3_pInt), & - (R(i,j,iGrain), j = 1_pInt,3_pInt), & - (D(i,j,iGrain), j = 1_pInt,3_pInt) + (R(i,j,iGrain), j = 1_pInt,3_pInt), & + (D(i,j,iGrain), j = 1_pInt,3_pInt) enddo write(6,*)' ' enddo - !$OMP END CRITICAL (write2out) endif +#endif !------------------------------------------------------------------------------------------------ ! computing the residual stress from the balance of traction at all (interior) interfaces do iNum = 1_pInt,nIntFaceTot - faceID = homogenization_RGC_interface1to4(iNum,homID) ! identifying the interface ID in local coordinate system (4-dimensional index) + faceID = interface1to4(iNum,param(instance)%Nconstituents) ! identifying the interface ID in local coordinate system (4-dimensional index) !-------------------------------------------------------------------------------------------------- ! identify the left/bottom/back grain (-|N) iGr3N = faceID(2:4) ! identifying the grain ID in local coordinate system (3-dimensional index) - iGrN = homogenization_RGC_grain3to1(iGr3N,homID) ! translate the local grain ID into global coordinate system (1-dimensional index) - intFaceN = homogenization_RGC_getInterface(2_pInt*faceID(1),iGr3N) - normN = homogenization_RGC_interfaceNormal(intFaceN,ip,el) ! get the interface normal + iGrN = grain3to1(iGr3N,param(instance)%Nconstituents) ! translate the local grain ID into global coordinate system (1-dimensional index) + intFaceN = getInterface(2_pInt*faceID(1),iGr3N) + normN = interfaceNormal(intFaceN,instance,of) !-------------------------------------------------------------------------------------------------- ! identify the right/up/front grain (+|P) iGr3P = iGr3N iGr3P(faceID(1)) = iGr3N(faceID(1))+1_pInt ! identifying the grain ID in local coordinate system (3-dimensional index) - iGrP = homogenization_RGC_grain3to1(iGr3P,homID) ! translate the local grain ID into global coordinate system (1-dimensional index) - intFaceP = homogenization_RGC_getInterface(2_pInt*faceID(1)-1_pInt,iGr3P) - normP = homogenization_RGC_interfaceNormal(intFaceP,ip,el) ! get the interface normal + iGrP = grain3to1(iGr3P,param(instance)%Nconstituents) ! translate the local grain ID into global coordinate system (1-dimensional index) + intFaceP = getInterface(2_pInt*faceID(1)-1_pInt,iGr3P) + normP = interfaceNormal(intFaceP,instance,of) !-------------------------------------------------------------------------------------------------- ! compute the residual of traction at the interface (in local system, 4-dimensional index) @@ -519,29 +466,25 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) enddo enddo -!-------------------------------------------------------------------------------------------------- -! debugging the residual stress +#ifdef DEBUG if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt) then - !$OMP CRITICAL (write2out) write(6,'(1x,a30,1x,i3)')'Traction at interface: ',iNum write(6,'(1x,3(e15.8,1x))')(tract(iNum,j), j = 1_pInt,3_pInt) write(6,*)' ' - !$OMP END CRITICAL (write2out) endif +#endif enddo !-------------------------------------------------------------------------------------------------- ! convergence check for stress residual stresMax = maxval(abs(P)) ! get the maximum of first Piola-Kirchhoff (material) stress - stresLoc = int(maxloc(abs(P)),pInt) ! get the location of the maximum stress residMax = maxval(abs(tract)) ! get the maximum of the residual - residLoc = int(maxloc(abs(tract)),pInt) ! get the position of the maximum residual -!-------------------------------------------------------------------------------------------------- -! Debugging the convergent criteria +#ifdef DEBUG if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt & - .and. debug_e == el .and. debug_i == ip) then - !$OMP CRITICAL (write2out) + .and. prm%of_debug == of) then + stresLoc = int(maxloc(abs(P)),pInt) ! get the location of the maximum stress + residLoc = int(maxloc(abs(tract)),pInt) ! get the position of the maximum residual write(6,'(1x,a)')' ' write(6,'(1x,a,1x,i2,1x,i4)')'RGC residual check ...',ip,el write(6,'(1x,a15,1x,e15.8,1x,a7,i3,1x,a12,i2,i2)')'Max stress: ',stresMax, & @@ -549,8 +492,8 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) write(6,'(1x,a15,1x,e15.8,1x,a7,i3,1x,a12,i2)')'Max residual: ',residMax, & '@ iface',residLoc(1),'in direction',residLoc(2) flush(6) - !$OMP END CRITICAL (write2out) endif +#endif homogenization_RGC_updateState = .false. @@ -558,86 +501,63 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) ! If convergence reached => done and happy if (residMax < relTol_RGC*stresMax .or. residMax < absTol_RGC) then homogenization_RGC_updateState = .true. - +#ifdef DEBUG if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt & - .and. debug_e == el .and. debug_i == ip) then - !$OMP CRITICAL (write2out) - write(6,'(1x,a55,/)')'... done and happy' - flush(6) - !$OMP END CRITICAL (write2out) - endif + .and. prm%of_debug == of) write(6,'(1x,a55,/)')'... done and happy' + flush(6) +#endif !-------------------------------------------------------------------------------------------------- ! compute/update the state for postResult, i.e., all energy densities computed by time-integration - constitutiveWork = homogState(mappingHomogenization(2,ip,el))%state(3*nIntFaceTot+1,mappingHomogenization(1,ip,el)) - penaltyEnergy = homogState(mappingHomogenization(2,ip,el))%state(3*nIntFaceTot+5,mappingHomogenization(1,ip,el)) - do iGrain = 1_pInt,homogenization_Ngrains(mesh_element(3,el)) ! time-integration loop for the calculating the work and energy - do i = 1_pInt,3_pInt - do j = 1_pInt,3_pInt - constitutiveWork = constitutiveWork + P(i,j,iGrain)*(F(i,j,iGrain) - F0(i,j,iGrain))/real(nGrain,pReal) - penaltyEnergy = penaltyEnergy + R(i,j,iGrain)*(F(i,j,iGrain) - F0(i,j,iGrain))/real(nGrain,pReal) - enddo - enddo + do iGrain = 1_pInt,product(prm%Nconstituents) + do i = 1_pInt,3_pInt;do j = 1_pInt,3_pInt + stt%work(of) = stt%work(of) & + + P(i,j,iGrain)*(F(i,j,iGrain) - F0(i,j,iGrain))/real(nGrain,pReal) + stt%penaltyEnergy(of) = stt%penaltyEnergy(of) & + + R(i,j,iGrain)*(F(i,j,iGrain) - F0(i,j,iGrain))/real(nGrain,pReal) + enddo; enddo enddo - homogState(mappingHomogenization(2,ip,el))% & - state(3*nIntFaceTot+1,mappingHomogenization(1,ip,el)) = constitutiveWork ! the bulk mechanical/constitutive work - homogState(mappingHomogenization(2,ip,el))% & - state(3*nIntFaceTot+2,mappingHomogenization(1,ip,el)) = sum(NN(1,:))/real(nGrain,pReal) ! the overall mismatch of all interface normal to e1-direction - homogState(mappingHomogenization(2,ip,el))% & - state(3*nIntFaceTot+3,mappingHomogenization(1,ip,el)) = sum(NN(2,:))/real(nGrain,pReal) ! the overall mismatch of all interface normal to e2-direction - homogState(mappingHomogenization(2,ip,el))% & - state(3*nIntFaceTot+4,mappingHomogenization(1,ip,el)) = sum(NN(3,:))/real(nGrain,pReal) ! the overall mismatch of all interface normal to e3-direction - homogState(mappingHomogenization(2,ip,el))% & - state(3*nIntFaceTot+5,mappingHomogenization(1,ip,el)) = penaltyEnergy ! the overall penalty energy - homogState(mappingHomogenization(2,ip,el))% & - state(3*nIntFaceTot+6,mappingHomogenization(1,ip,el)) = volDiscrep ! the overall volume discrepancy - homogState(mappingHomogenization(2,ip,el))% & - state(3*nIntFaceTot+7,mappingHomogenization(1,ip,el)) = & - sum(abs(drelax))/dt/real(3_pInt*nIntFaceTot,pReal) ! the average rate of relaxation vectors - homogState(mappingHomogenization(2,ip,el))% & - state(3*nIntFaceTot+8,mappingHomogenization(1,ip,el)) = maxval(abs(drelax))/dt ! the maximum rate of relaxation vectors - if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt & - .and. debug_e == el .and. debug_i == ip) then - !$OMP CRITICAL (write2out) - write(6,'(1x,a30,1x,e15.8)') 'Constitutive work: ',constitutiveWork - write(6,'(1x,a30,3(1x,e15.8))')'Magnitude mismatch: ',sum(NN(1,:))/real(nGrain,pReal), & - sum(NN(2,:))/real(nGrain,pReal), & - sum(NN(3,:))/real(nGrain,pReal) - write(6,'(1x,a30,1x,e15.8)') 'Penalty energy: ',penaltyEnergy - write(6,'(1x,a30,1x,e15.8,/)') 'Volume discrepancy: ',volDiscrep - write(6,'(1x,a30,1x,e15.8)') 'Maximum relaxation rate: ',maxval(abs(drelax))/dt - write(6,'(1x,a30,1x,e15.8,/)') 'Average relaxation rate: ',sum(abs(drelax))/dt/real(3_pInt*nIntFaceTot,pReal) - flush(6) - !$OMP END CRITICAL (write2out) - endif + dst%mismatch(1:3,of) = sum(NN,2)/real(nGrain,pReal) + dst%relaxationRate_avg(of) = sum(abs(drelax))/dt/real(3_pInt*nIntFaceTot,pReal) + dst%relaxationRate_max(of) = maxval(abs(drelax))/dt + +#ifdef DEBUG + if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt & + .and. prm%of_debug == of) then + write(6,'(1x,a30,1x,e15.8)') 'Constitutive work: ',stt%work(of) + write(6,'(1x,a30,3(1x,e15.8))')'Magnitude mismatch: ',dst%mismatch(1,of), & + dst%mismatch(2,of), & + dst%mismatch(3,of) + write(6,'(1x,a30,1x,e15.8)') 'Penalty energy: ', stt%penaltyEnergy(of) + write(6,'(1x,a30,1x,e15.8,/)') 'Volume discrepancy: ', dst%volumeDiscrepancy(of) + write(6,'(1x,a30,1x,e15.8)') 'Maximum relaxation rate: ', dst%relaxationRate_max(of) + write(6,'(1x,a30,1x,e15.8,/)') 'Average relaxation rate: ', dst%relaxationRate_avg(of) + flush(6) + endif +#endif - deallocate(tract,resid,relax,drelax) return !-------------------------------------------------------------------------------------------------- ! if residual blows-up => done but unhappy elseif (residMax > relMax_RGC*stresMax .or. residMax > absMax_RGC) then ! try to restart when residual blows up exceeding maximum bound homogenization_RGC_updateState = [.true.,.false.] ! with direct cut-back - - if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt & - .and. debug_e == el .and. debug_i == ip) then - !$OMP CRITICAL (write2out) - write(6,'(1x,a55,/)')'... broken' - flush(6) - !$OMP END CRITICAL (write2out) - endif - deallocate(tract,resid,relax,drelax) - return - else ! proceed with computing the Jacobian and state update +#ifdef DEBUG if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt & - .and. debug_e == el .and. debug_i == ip) then - !$OMP CRITICAL (write2out) - write(6,'(1x,a55,/)')'... not yet done' - flush(6) - !$OMP END CRITICAL (write2out) - endif + .and. prm%of_debug == of) write(6,'(1x,a,/)') '... broken' + flush(6) +#endif + + return + + else ! proceed with computing the Jacobian and state update +#ifdef DEBUG + if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt & + .and. prm%of_debug == of) write(6,'(1x,a,/)') '... not yet done' + flush(6) +#endif endif @@ -649,21 +569,22 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) ! ... of the constitutive stress tangent, assembled from dPdF or material constitutive model "smatrix" allocate(smatrix(3*nIntFaceTot,3*nIntFaceTot), source=0.0_pReal) do iNum = 1_pInt,nIntFaceTot - faceID = homogenization_RGC_interface1to4(iNum,homID) ! assembling of local dPdF into global Jacobian matrix + faceID = interface1to4(iNum,param(instance)%Nconstituents) ! assembling of local dPdF into global Jacobian matrix !-------------------------------------------------------------------------------------------------- ! identify the left/bottom/back grain (-|N) iGr3N = faceID(2:4) ! identifying the grain ID in local coordinate sytem - iGrN = homogenization_RGC_grain3to1(iGr3N,homID) ! translate into global grain ID - intFaceN = homogenization_RGC_getInterface(2_pInt*faceID(1),iGr3N) ! identifying the connecting interface in local coordinate system - normN = homogenization_RGC_interfaceNormal(intFaceN,ip,el) ! get the interface normal - do iFace = 1_pInt,nFace - intFaceN = homogenization_RGC_getInterface(iFace,iGr3N) ! identifying all interfaces that influence relaxation of the above interface - mornN = homogenization_RGC_interfaceNormal(intFaceN,ip,el) ! get normal of the interfaces - iMun = homogenization_RGC_interface4to1(intFaceN,homID) ! translate the interfaces ID into local 4-dimensional index + iGrN = grain3to1(iGr3N,param(instance)%Nconstituents) ! translate into global grain ID + intFaceN = getInterface(2_pInt*faceID(1),iGr3N) ! identifying the connecting interface in local coordinate system + normN = interfaceNormal(intFaceN,instance,of) + do iFace = 1_pInt,6_pInt + intFaceN = getInterface(iFace,iGr3N) ! identifying all interfaces that influence relaxation of the above interface + mornN = interfaceNormal(intFaceN,instance,of) + iMun = interface4to1(intFaceN,param(instance)%Nconstituents) ! translate the interfaces ID into local 4-dimensional index if (iMun > 0) then ! get the corresponding tangent do i=1_pInt,3_pInt; do j=1_pInt,3_pInt; do k=1_pInt,3_pInt; do l=1_pInt,3_pInt - smatrix(3*(iNum-1)+i,3*(iMun-1)+j) = smatrix(3*(iNum-1)+i,3*(iMun-1)+j) + dPdF(i,k,j,l,iGrN)*normN(k)*mornN(l) + smatrix(3*(iNum-1)+i,3*(iMun-1)+j) = smatrix(3*(iNum-1)+i,3*(iMun-1)+j) & + + dPdF(i,k,j,l,iGrN)*normN(k)*mornN(l) enddo;enddo;enddo;enddo ! projecting the material tangent dPdF into the interface ! to obtain the Jacobian matrix contribution of dPdF @@ -674,33 +595,32 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) ! identify the right/up/front grain (+|P) iGr3P = iGr3N iGr3P(faceID(1)) = iGr3N(faceID(1))+1_pInt ! identifying the grain ID in local coordinate sytem - iGrP = homogenization_RGC_grain3to1(iGr3P,homID) ! translate into global grain ID - intFaceP = homogenization_RGC_getInterface(2_pInt*faceID(1)-1_pInt,iGr3P) ! identifying the connecting interface in local coordinate system - normP = homogenization_RGC_interfaceNormal(intFaceP,ip,el) ! get the interface normal - do iFace = 1_pInt,nFace - intFaceP = homogenization_RGC_getInterface(iFace,iGr3P) ! identifying all interfaces that influence relaxation of the above interface - mornP = homogenization_RGC_interfaceNormal(intFaceP,ip,el) ! get normal of the interfaces - iMun = homogenization_RGC_interface4to1(intFaceP,homID) ! translate the interfaces ID into local 4-dimensional index + iGrP = grain3to1(iGr3P,param(instance)%Nconstituents) ! translate into global grain ID + intFaceP = getInterface(2_pInt*faceID(1)-1_pInt,iGr3P) ! identifying the connecting interface in local coordinate system + normP = interfaceNormal(intFaceP,instance,of) + do iFace = 1_pInt,6_pInt + intFaceP = getInterface(iFace,iGr3P) ! identifying all interfaces that influence relaxation of the above interface + mornP = interfaceNormal(intFaceP,instance,of) + iMun = interface4to1(intFaceP,param(instance)%Nconstituents) ! translate the interfaces ID into local 4-dimensional index if (iMun > 0_pInt) then ! get the corresponding tangent do i=1_pInt,3_pInt; do j=1_pInt,3_pInt; do k=1_pInt,3_pInt; do l=1_pInt,3_pInt - smatrix(3*(iNum-1)+i,3*(iMun-1)+j) = smatrix(3*(iNum-1)+i,3*(iMun-1)+j) + dPdF(i,k,j,l,iGrP)*normP(k)*mornP(l) + smatrix(3*(iNum-1)+i,3*(iMun-1)+j) = smatrix(3*(iNum-1)+i,3*(iMun-1)+j) & + + dPdF(i,k,j,l,iGrP)*normP(k)*mornP(l) enddo;enddo;enddo;enddo endif enddo enddo -!-------------------------------------------------------------------------------------------------- -! debugging the global Jacobian matrix of stress tangent +#ifdef DEBUG if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt) then - !$OMP CRITICAL (write2out) write(6,'(1x,a30)')'Jacobian matrix of stress' do i = 1_pInt,3_pInt*nIntFaceTot write(6,'(1x,100(e11.4,1x))')(smatrix(i,j), j = 1_pInt,3_pInt*nIntFaceTot) enddo write(6,*)' ' flush(6) - !$OMP END CRITICAL (write2out) endif +#endif !-------------------------------------------------------------------------------------------------- ! ... of the stress penalty tangent (mismatch penalty and volume penalty, computed using numerical @@ -708,34 +628,35 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) allocate(pmatrix(3*nIntFaceTot,3*nIntFaceTot), source=0.0_pReal) allocate(p_relax(3*nIntFaceTot), source=0.0_pReal) allocate(p_resid(3*nIntFaceTot), source=0.0_pReal) + do ipert = 1_pInt,3_pInt*nIntFaceTot p_relax = relax p_relax(ipert) = relax(ipert) + pPert_RGC ! perturb the relaxation vector - homogState(mappingHomogenization(2,ip,el))%state(1:3*nIntFaceTot,mappingHomogenization(1,ip,el)) = p_relax - call homogenization_RGC_grainDeformation(pF,avgF,ip,el) ! compute the grains deformation from perturbed state - call homogenization_RGC_stressPenalty(pR,pNN,avgF,pF,ip,el,homID) ! compute stress penalty due to interface mismatch from perturbed state - call homogenization_RGC_volumePenalty(pD,volDiscrep,pF,avgF,ip,el) ! compute stress penalty due to volume discrepancy from perturbed state + stt%relaxationVector(:,of) = p_relax + call grainDeformation(pF,avgF,instance,of) ! rain deformation from perturbed state + call stressPenalty(pR,DevNull, avgF,pF,ip,el,instance,of) ! stress penalty due to interface mismatch from perturbed state + call volumePenalty(pD,devNull(1,1), avgF,pF,nGrain,instance,of) ! stress penalty due to volume discrepancy from perturbed state !-------------------------------------------------------------------------------------------------- ! computing the global stress residual array from the perturbed state p_resid = 0.0_pReal do iNum = 1_pInt,nIntFaceTot - faceID = homogenization_RGC_interface1to4(iNum,homID) ! identifying the interface ID in local coordinate system (4-dimensional index) + faceID = interface1to4(iNum,param(instance)%Nconstituents) ! identifying the interface ID in local coordinate system (4-dimensional index) !-------------------------------------------------------------------------------------------------- ! identify the left/bottom/back grain (-|N) - iGr3N = faceID(2:4) ! identifying the grain ID in local coordinate system (3-dimensional index) - iGrN = homogenization_RGC_grain3to1(iGr3N,homID) ! translate the local grain ID into global coordinate system (1-dimensional index) - intFaceN = homogenization_RGC_getInterface(2_pInt*faceID(1),iGr3N) ! identifying the interface ID of the grain - normN = homogenization_RGC_interfaceNormal(intFaceN,ip,el) ! get the corresponding interface normal + iGr3N = faceID(2:4) ! identify the grain ID in local coordinate system (3-dimensional index) + iGrN = grain3to1(iGr3N,param(instance)%Nconstituents) ! translate the local grain ID into global coordinate system (1-dimensional index) + intFaceN = getInterface(2_pInt*faceID(1),iGr3N) ! identify the interface ID of the grain + normN = interfaceNormal(intFaceN,instance,of) !-------------------------------------------------------------------------------------------------- ! identify the right/up/front grain (+|P) iGr3P = iGr3N - iGr3P(faceID(1)) = iGr3N(faceID(1))+1_pInt ! identifying the grain ID in local coordinate system (3-dimensional index) - iGrP = homogenization_RGC_grain3to1(iGr3P,homID) ! translate the local grain ID into global coordinate system (1-dimensional index) - intFaceP = homogenization_RGC_getInterface(2_pInt*faceID(1)-1_pInt,iGr3P) ! identifying the interface ID of the grain - normP = homogenization_RGC_interfaceNormal(intFaceP,ip,el) ! get the corresponding normal + iGr3P(faceID(1)) = iGr3N(faceID(1))+1_pInt ! identify the grain ID in local coordinate system (3-dimensional index) + iGrP = grain3to1(iGr3P,param(instance)%Nconstituents) ! translate the local grain ID into global coordinate system (1-dimensional index) + intFaceP = getInterface(2_pInt*faceID(1)-1_pInt,iGr3P) ! identify the interface ID of the grain + normP = interfaceNormal(intFaceP,instance,of) !-------------------------------------------------------------------------------------------------- ! compute the residual stress (contribution of mismatch and volume penalties) from perturbed state @@ -750,18 +671,16 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) pmatrix(:,ipert) = p_resid/pPert_RGC enddo -!-------------------------------------------------------------------------------------------------- -! debugging the global Jacobian matrix of penalty tangent +#ifdef DEBUG if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0_pInt) then - !$OMP CRITICAL (write2out) write(6,'(1x,a30)')'Jacobian matrix of penalty' do i = 1_pInt,3_pInt*nIntFaceTot write(6,'(1x,100(e11.4,1x))')(pmatrix(i,j), j = 1_pInt,3_pInt*nIntFaceTot) enddo write(6,*)' ' flush(6) - !$OMP END CRITICAL (write2out) endif +#endif !-------------------------------------------------------------------------------------------------- ! ... of the numerical viscosity traction "rmatrix" @@ -769,65 +688,56 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) forall (i=1_pInt:3_pInt*nIntFaceTot) & rmatrix(i,i) = viscModus_RGC*viscPower_RGC/(refRelaxRate_RGC*dt)* & ! tangent due to numerical viscosity traction appears (abs(drelax(i))/(refRelaxRate_RGC*dt))**(viscPower_RGC - 1.0_pReal) ! only in the main diagonal term - - - -!-------------------------------------------------------------------------------------------------- -! debugging the global Jacobian matrix of numerical viscosity tangent + +#ifdef DEBUG if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0_pInt) then - !$OMP CRITICAL (write2out) write(6,'(1x,a30)')'Jacobian matrix of penalty' do i = 1_pInt,3_pInt*nIntFaceTot write(6,'(1x,100(e11.4,1x))')(rmatrix(i,j), j = 1_pInt,3_pInt*nIntFaceTot) enddo write(6,*)' ' flush(6) - !$OMP END CRITICAL (write2out) endif +#endif !-------------------------------------------------------------------------------------------------- ! The overall Jacobian matrix summarizing contributions of smatrix, pmatrix, rmatrix allocate(jmatrix(3*nIntFaceTot,3*nIntFaceTot)); jmatrix = smatrix + pmatrix + rmatrix - + +#ifdef DEBUG if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0_pInt) then - !$OMP CRITICAL (write2out) write(6,'(1x,a30)')'Jacobian matrix (total)' do i = 1_pInt,3_pInt*nIntFaceTot write(6,'(1x,100(e11.4,1x))')(jmatrix(i,j), j = 1_pInt,3_pInt*nIntFaceTot) enddo write(6,*)' ' flush(6) - !$OMP END CRITICAL (write2out) endif +#endif !-------------------------------------------------------------------------------------------------- ! computing the update of the state variable (relaxation vectors) using the Jacobian matrix allocate(jnverse(3_pInt*nIntFaceTot,3_pInt*nIntFaceTot),source=0.0_pReal) - call math_invert(size(jmatrix,1),jmatrix,jnverse,error) ! Compute the inverse of the overall Jacobian matrix + call math_invert2(jnverse,error,jmatrix) -!-------------------------------------------------------------------------------------------------- -! debugging the inverse Jacobian matrix +#ifdef DEBUG if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0_pInt) then - !$OMP CRITICAL (write2out) write(6,'(1x,a30)')'Jacobian inverse' do i = 1_pInt,3_pInt*nIntFaceTot write(6,'(1x,100(e11.4,1x))')(jnverse(i,j), j = 1_pInt,3_pInt*nIntFaceTot) enddo write(6,*)' ' flush(6) - !$OMP END CRITICAL (write2out) endif +#endif !-------------------------------------------------------------------------------------------------- ! calculate the state update (global relaxation vectors) for the next Newton-Raphson iteration drelax = 0.0_pReal - do i = 1_pInt,3_pInt*nIntFaceTot - do j = 1_pInt,3_pInt*nIntFaceTot - drelax(i) = drelax(i) - jnverse(i,j)*resid(j) ! Calculate the correction for the state variable - enddo - enddo - relax = relax + drelax ! Updateing the state variable for the next iteration - homogState(mappingHomogenization(2,ip,el))%state(1:3*nIntFaceTot,mappingHomogenization(1,ip,el)) = relax + do i = 1_pInt,3_pInt*nIntFaceTot;do j = 1_pInt,3_pInt*nIntFaceTot + drelax(i) = drelax(i) - jnverse(i,j)*resid(j) ! Calculate the correction for the state variable + enddo; enddo + stt%relaxationVector(:,of) = relax + drelax ! Updateing the state variable for the next iteration if (any(abs(drelax) > maxdRelax_RGC)) then ! Forcing cutback when the incremental change of relaxation vector becomes too large homogenization_RGC_updateState = [.true.,.false.] !$OMP CRITICAL (write2out) @@ -837,20 +747,303 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) !$OMP END CRITICAL (write2out) endif -!-------------------------------------------------------------------------------------------------- -! debugging the return state +#ifdef DEBUG if (iand(debug_homogenization, debug_levelExtensive) > 0_pInt) then - !$OMP CRITICAL (write2out) write(6,'(1x,a30)')'Returned state: ' - do i = 1_pInt,3_pInt*nIntFaceTot - write(6,'(1x,2(e15.8,1x))')homogState(mappingHomogenization(2,ip,el))%state(i,mappingHomogenization(1,ip,el)) + do i = 1_pInt,size(stt%relaxationVector(:,of)) + write(6,'(1x,2(e15.8,1x))') stt%relaxationVector(i,of) enddo write(6,*)' ' flush(6) - !$OMP END CRITICAL (write2out) endif +#endif - deallocate(tract,resid,jmatrix,jnverse,relax,drelax,pmatrix,smatrix,p_relax,p_resid) + end associate + + contains + !-------------------------------------------------------------------------------------------------- + !> @brief calculate stress-like penalty due to deformation mismatch + !-------------------------------------------------------------------------------------------------- + subroutine stressPenalty(rPen,nMis,avgF,fDef,ip,el,instance,of) + use math, only: & + math_civita + use numerics, only: & + xSmoo_RGC + + implicit none + real(pReal), dimension (:,:,:), intent(out) :: rPen !< stress-like penalty + real(pReal), dimension (:,:), intent(out) :: nMis !< total amount of mismatch + + real(pReal), dimension (:,:,:), intent(in) :: fDef !< deformation gradients + real(pReal), dimension (3,3), intent(in) :: avgF !< initial effective stretch tensor + integer(pInt), intent(in) :: ip,el,instance,of + + integer(pInt), dimension (4) :: intFace + integer(pInt), dimension (3) :: iGrain3,iGNghb3,nGDim + real(pReal), dimension (3,3) :: gDef,nDef + real(pReal), dimension (3) :: nVect,surfCorr + real(pReal), dimension (2) :: Gmoduli + integer(pInt) :: iGrain,iGNghb,iFace,i,j,k,l + real(pReal) :: muGrain,muGNghb,nDefNorm,bgGrain,bgGNghb + real(pReal), parameter :: nDefToler = 1.0e-10_pReal +#ifdef DEBUG + logical :: debugActive +#endif + + nGDim = param(instance)%Nconstituents + rPen = 0.0_pReal + nMis = 0.0_pReal + + !-------------------------------------------------------------------------------------------------- + ! get the correction factor the modulus of penalty stress representing the evolution of area of + ! the interfaces due to deformations + + surfCorr = surfaceCorrection(avgF,instance,of) + + associate(prm => param(instance)) + +#ifdef DEBUG + debugActive = iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt & + .and. prm%of_debug == of + + if (debugActive) then + write(6,'(1x,a20,2(1x,i3))')'Correction factor: ',ip,el + write(6,*) surfCorr + endif +#endif + + !-------------------------------------------------------------------------------------------------- + ! computing the mismatch and penalty stress tensor of all grains + grainLoop: do iGrain = 1_pInt,product(prm%Nconstituents) + Gmoduli = equivalentModuli(iGrain,ip,el) + muGrain = Gmoduli(1) ! collecting the equivalent shear modulus of grain + bgGrain = Gmoduli(2) ! and the lengthh of Burgers vector + iGrain3 = grain1to3(iGrain,prm%Nconstituents) ! get the grain ID in local 3-dimensional index (x,y,z)-position + + interfaceLoop: do iFace = 1_pInt,6_pInt + intFace = getInterface(iFace,iGrain3) ! get the 4-dimensional index of the interface in local numbering system of the grain + nVect = interfaceNormal(intFace,instance,of) + iGNghb3 = iGrain3 ! identify the neighboring grain across the interface + iGNghb3(abs(intFace(1))) = iGNghb3(abs(intFace(1))) & + + int(real(intFace(1),pReal)/real(abs(intFace(1)),pReal),pInt) + where(iGNghb3 < 1) iGNghb3 = nGDim + where(iGNghb3 >nGDim) iGNghb3 = 1_pInt + iGNghb = grain3to1(iGNghb3,prm%Nconstituents) ! get the ID of the neighboring grain + Gmoduli = equivalentModuli(iGNghb,ip,el) ! collect the shear modulus and Burgers vector of the neighbor + muGNghb = Gmoduli(1) + bgGNghb = Gmoduli(2) + gDef = 0.5_pReal*(fDef(1:3,1:3,iGNghb) - fDef(1:3,1:3,iGrain)) ! difference/jump in deformation gradeint across the neighbor + + !-------------------------------------------------------------------------------------------------- + ! compute the mismatch tensor of all interfaces + nDefNorm = 0.0_pReal + nDef = 0.0_pReal + do i = 1_pInt,3_pInt; do j = 1_pInt,3_pInt + do k = 1_pInt,3_pInt; do l = 1_pInt,3_pInt + nDef(i,j) = nDef(i,j) - nVect(k)*gDef(i,l)*math_civita(j,k,l) ! compute the interface mismatch tensor from the jump of deformation gradient + enddo; enddo + nDefNorm = nDefNorm + nDef(i,j)**2.0_pReal ! compute the norm of the mismatch tensor + enddo; enddo + nDefNorm = max(nDefToler,sqrt(nDefNorm)) ! approximation to zero mismatch if mismatch is zero (singularity) + nMis(abs(intFace(1)),iGrain) = nMis(abs(intFace(1)),iGrain) + nDefNorm ! total amount of mismatch experienced by the grain (at all six interfaces) +#ifdef DEBUG + if (debugActive) then + write(6,'(1x,a20,i2,1x,a20,1x,i3)')'Mismatch to face: ',intFace(1),'neighbor grain: ',iGNghb + write(6,*) transpose(nDef) + write(6,'(1x,a20,e11.4)')'with magnitude: ',nDefNorm + endif +#endif + + !-------------------------------------------------------------------------------------------------- + ! compute the stress penalty of all interfaces + do i = 1_pInt,3_pInt; do j = 1_pInt,3_pInt; do k = 1_pInt,3_pInt; do l = 1_pInt,3_pInt + rPen(i,j,iGrain) = rPen(i,j,iGrain) + 0.5_pReal*(muGrain*bgGrain + muGNghb*bgGNghb)*prm%xiAlpha & + *surfCorr(abs(intFace(1)))/prm%dAlpha(abs(intFace(1))) & + *cosh(prm%ciAlpha*nDefNorm) & + *0.5_pReal*nVect(l)*nDef(i,k)/nDefNorm*math_civita(k,l,j) & + *tanh(nDefNorm/xSmoo_RGC) + enddo; enddo;enddo; enddo + enddo interfaceLoop +#ifdef DEBUG + if (debugActive) then + write(6,'(1x,a20,i2)')'Penalty of grain: ',iGrain + write(6,*) transpose(rPen(1:3,1:3,iGrain)) + endif +#endif + + enddo grainLoop + + end associate + + end subroutine stressPenalty + + + !-------------------------------------------------------------------------------------------------- + !> @brief calculate stress-like penalty due to volume discrepancy + !-------------------------------------------------------------------------------------------------- + subroutine volumePenalty(vPen,vDiscrep,fAvg,fDef,nGrain,instance,of) + use math, only: & + math_det33, & + math_inv33 + use numerics, only: & + maxVolDiscr_RGC,& + volDiscrMod_RGC,& + volDiscrPow_RGC + + implicit none + real(pReal), dimension (:,:,:), intent(out) :: vPen ! stress-like penalty due to volume + real(pReal), intent(out) :: vDiscrep ! total volume discrepancy + + real(pReal), dimension (:,:,:), intent(in) :: fDef ! deformation gradients + real(pReal), dimension (3,3), intent(in) :: fAvg ! overall deformation gradient + integer(pInt), intent(in) :: & + Ngrain, & + instance, & + of + + real(pReal), dimension(size(vPen,3)) :: gVol + integer(pInt) :: i + + !-------------------------------------------------------------------------------------------------- + ! compute the volumes of grains and of cluster + vDiscrep = math_det33(fAvg) ! compute the volume of the cluster + do i = 1_pInt,nGrain + gVol(i) = math_det33(fDef(1:3,1:3,i)) ! compute the volume of individual grains + vDiscrep = vDiscrep - gVol(i)/real(nGrain,pReal) ! calculate the difference/dicrepancy between + ! the volume of the cluster and the the total volume of grains + enddo + + !-------------------------------------------------------------------------------------------------- + ! calculate the stress and penalty due to volume discrepancy + vPen = 0.0_pReal + do i = 1_pInt,nGrain + vPen(:,:,i) = -1.0_pReal/real(nGrain,pReal)*volDiscrMod_RGC*volDiscrPow_RGC/maxVolDiscr_RGC* & + sign((abs(vDiscrep)/maxVolDiscr_RGC)**(volDiscrPow_RGC - 1.0),vDiscrep)* & + gVol(i)*transpose(math_inv33(fDef(:,:,i))) + +#ifdef DEBUG + if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt & + .and. param(instance)%of_debug == of) then + write(6,'(1x,a30,i2)')'Volume penalty of grain: ',i + write(6,*) transpose(vPen(:,:,i)) + endif +#endif + enddo + + end subroutine volumePenalty + + + !-------------------------------------------------------------------------------------------------- + !> @brief compute the correction factor accouted for surface evolution (area change) due to + ! deformation + !-------------------------------------------------------------------------------------------------- + function surfaceCorrection(avgF,instance,of) + use math, only: & + math_invert33, & + math_mul33x33 + + implicit none + real(pReal), dimension(3) :: surfaceCorrection + real(pReal), dimension(3,3), intent(in) :: avgF !< average F + integer(pInt), intent(in) :: & + instance, & + of + real(pReal), dimension(3,3) :: invC + real(pReal), dimension(3) :: nVect + real(pReal) :: detF + integer(pInt) :: i,j,iBase + logical :: error + + call math_invert33(math_mul33x33(transpose(avgF),avgF),invC,detF,error) + + surfaceCorrection = 0.0_pReal + do iBase = 1_pInt,3_pInt + nVect = interfaceNormal([iBase,1_pInt,1_pInt,1_pInt],instance,of) + do i = 1_pInt,3_pInt; do j = 1_pInt,3_pInt + surfaceCorrection(iBase) = surfaceCorrection(iBase) + invC(i,j)*nVect(i)*nVect(j) ! compute the component of (the inverse of) the stretch in the direction of the normal + enddo; enddo + surfaceCorrection(iBase) = sqrt(surfaceCorrection(iBase))*detF ! get the surface correction factor (area contraction/enlargement) + enddo + + end function surfaceCorrection + + + !-------------------------------------------------------------------------------------------------- + !> @brief compute the equivalent shear and bulk moduli from the elasticity tensor + !-------------------------------------------------------------------------------------------------- + function equivalentModuli(grainID,ip,el) + use constitutive, only: & + constitutive_homogenizedC + + implicit none + real(pReal), dimension(2) :: equivalentModuli + integer(pInt), intent(in) :: & + grainID,& + ip, & !< integration point number + el !< element number + real(pReal), dimension(6,6) :: elasTens + real(pReal) :: & + cEquiv_11, & + cEquiv_12, & + cEquiv_44 + + elasTens = constitutive_homogenizedC(grainID,ip,el) + + !-------------------------------------------------------------------------------------------------- + ! compute the equivalent shear modulus after Turterltaub and Suiker, JMPS (2005) + cEquiv_11 = (elasTens(1,1) + elasTens(2,2) + elasTens(3,3))/3.0_pReal + cEquiv_12 = (elasTens(1,2) + elasTens(2,3) + elasTens(3,1) + & + elasTens(1,3) + elasTens(2,1) + elasTens(3,2))/6.0_pReal + cEquiv_44 = (elasTens(4,4) + elasTens(5,5) + elasTens(6,6))/3.0_pReal + equivalentModuli(1) = 0.2_pReal*(cEquiv_11 - cEquiv_12) + 0.6_pReal*cEquiv_44 + + !-------------------------------------------------------------------------------------------------- + ! obtain the length of Burgers vector (could be model dependend) + equivalentModuli(2) = 2.5e-10_pReal + + end function equivalentModuli + + + !-------------------------------------------------------------------------------------------------- + !> @brief calculating the grain deformation gradient (the same with + ! homogenization_RGC_partitionDeformation, but used only for perturbation scheme) + !-------------------------------------------------------------------------------------------------- + subroutine grainDeformation(F, avgF, instance, of) + + implicit none + real(pReal), dimension (:,:,:), intent(out) :: F !< partioned F per grain + + real(pReal), dimension (:,:), intent(in) :: avgF !< averaged F + integer(pInt), intent(in) :: & + instance, & + of + + real(pReal), dimension (3) :: aVect,nVect + integer(pInt), dimension (4) :: intFace + integer(pInt), dimension (3) :: iGrain3 + integer(pInt) :: iGrain,iFace,i,j + + !------------------------------------------------------------------------------------------------- + ! compute the deformation gradient of individual grains due to relaxations + + associate(prm => param(instance)) + + F = 0.0_pReal + do iGrain = 1_pInt,product(prm%Nconstituents) + iGrain3 = grain1to3(iGrain,prm%Nconstituents) + do iFace = 1_pInt,6_pInt + intFace = getInterface(iFace,iGrain3) + aVect = relaxationVector(intFace,instance,of) + nVect = interfaceNormal(intFace,instance,of) + forall (i=1_pInt:3_pInt,j=1_pInt:3_pInt) & + F(i,j,iGrain) = F(i,j,iGrain) + aVect(i)*nVect(j) ! effective relaxations + enddo + F(1:3,1:3,iGrain) = F(1:3,1:3,iGrain) + avgF ! relaxed deformation gradient + enddo + + end associate + + end subroutine grainDeformation end function homogenization_RGC_updateState @@ -858,51 +1051,18 @@ end function homogenization_RGC_updateState !-------------------------------------------------------------------------------------------------- !> @brief derive average stress and stiffness from constituent quantities !-------------------------------------------------------------------------------------------------- -subroutine homogenization_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,el) - use debug, only: & - debug_level, & - debug_homogenization,& - debug_levelExtensive - use mesh, only: mesh_element - use material, only: & - homogenization_maxNgrains, & - homogenization_Ngrains, & - homogenization_typeInstance - use math, only: math_Plain3333to99 - +subroutine homogenization_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,instance) + implicit none real(pReal), dimension (3,3), intent(out) :: avgP !< average stress at material point real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF !< average stiffness at material point - real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: P !< array of current grain stresses - real(pReal), dimension (3,3,3,3,homogenization_maxNgrains), intent(in) :: dPdF !< array of current grain stiffnesses - integer(pInt), intent(in) :: el !< element number - real(pReal), dimension (9,9) :: dPdF99 - integer(pInt) :: homID, i, j, Ngrains, iGrain + real(pReal), dimension (:,:,:), intent(in) :: P !< partitioned stresses + real(pReal), dimension (:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses + integer(pInt), intent(in) :: instance - homID = homogenization_typeInstance(mesh_element(3,el)) - Ngrains = homogenization_Ngrains(mesh_element(3,el)) - -!-------------------------------------------------------------------------------------------------- -! debugging the grain tangent - if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0_pInt) then - !$OMP CRITICAL (write2out) - do iGrain = 1_pInt,Ngrains - dPdF99 = math_Plain3333to99(dPdF(1:3,1:3,1:3,1:3,iGrain)) - write(6,'(1x,a30,1x,i3)')'Stress tangent of grain: ',iGrain - do i = 1_pInt,9_pInt - write(6,'(1x,(e15.8,1x))') (dPdF99(i,j), j = 1_pInt,9_pInt) - enddo - write(6,*)' ' - enddo - flush(6) - !$OMP END CRITICAL (write2out) - endif - -!-------------------------------------------------------------------------------------------------- -! computing the average first Piola-Kirchhoff stress P and the average tangent dPdF - avgP = sum(P,3)/real(Ngrains,pReal) - dAvgPdAvgF = sum(dPdF,5)/real(Ngrains,pReal) + avgP = sum(P,3) /real(product(param(instance)%Nconstituents),pReal) + dAvgPdAvgF = sum(dPdF,5)/real(product(param(instance)%Nconstituents),pReal) end subroutine homogenization_RGC_averageStressAndItsTangent @@ -910,632 +1070,271 @@ end subroutine homogenization_RGC_averageStressAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief return array of homogenization results for post file inclusion !-------------------------------------------------------------------------------------------------- -pure function homogenization_RGC_postResults(ip,el,avgP,avgF) - use mesh, only: & - mesh_element, & - mesh_ipCoordinates - use material, only: & - homogenization_typeInstance,& - homogState, & - mappingHomogenization, & - homogenization_Noutput - +pure function homogenization_RGC_postResults(instance,of) result(postResults) + implicit none integer(pInt), intent(in) :: & - ip, & !< integration point number - el !< element number - real(pReal), dimension(3,3), intent(in) :: & - avgP, & !< average stress at material point - avgF !< average deformation gradient at material point + instance, & + of + + integer(pInt) :: & + o,c + real(pReal), dimension(sum(homogenization_RGC_sizePostResult(:,instance))) :: & + postResults - integer(pInt) homID,o,c,nIntFaceTot - real(pReal), dimension(homogenization_RGC_sizePostResults(homogenization_typeInstance(mesh_element(3,el)))) :: & - homogenization_RGC_postResults - - homID = homogenization_typeInstance(mesh_element(3,el)) - nIntFaceTot=(homogenization_RGC_Ngrains(1,homID)-1_pInt)*homogenization_RGC_Ngrains(2,homID)*homogenization_RGC_Ngrains(3,homID)& - + homogenization_RGC_Ngrains(1,homID)*(homogenization_RGC_Ngrains(2,homID)-1_pInt)*homogenization_RGC_Ngrains(3,homID)& - + homogenization_RGC_Ngrains(1,homID)*homogenization_RGC_Ngrains(2,homID)*(homogenization_RGC_Ngrains(3,homID)-1_pInt) + associate(stt => state(instance), dst => dependentState(instance), prm => param(instance)) c = 0_pInt - homogenization_RGC_postResults = 0.0_pReal - do o = 1_pInt,homogenization_Noutput(mesh_element(3,el)) - select case(homogenization_RGC_outputID(o,homID)) - case (avgdefgrad_ID) - homogenization_RGC_postResults(c+1_pInt:c+9_pInt) = reshape(transpose(avgF),[9]) - c = c + 9_pInt - case (avgfirstpiola_ID) - homogenization_RGC_postResults(c+1_pInt:c+9_pInt) = reshape(transpose(avgP),[9]) - c = c + 9_pInt - case (ipcoords_ID) - homogenization_RGC_postResults(c+1_pInt:c+3_pInt) = mesh_ipCoordinates(1:3,ip,el) ! current ip coordinates - c = c + 3_pInt + + outputsLoop: do o = 1_pInt,size(prm%outputID) + select case(prm%outputID(o)) + case (constitutivework_ID) - homogenization_RGC_postResults(c+1) = homogState(mappingHomogenization(2,ip,el))% & - state(3*nIntFaceTot+1,mappingHomogenization(1,ip,el)) + postResults(c+1) = stt%work(of) c = c + 1_pInt case (magnitudemismatch_ID) - homogenization_RGC_postResults(c+1) = homogState(mappingHomogenization(2,ip,el))% & - state(3*nIntFaceTot+2,mappingHomogenization(1,ip,el)) - homogenization_RGC_postResults(c+2) = homogState(mappingHomogenization(2,ip,el))% & - state(3*nIntFaceTot+3,mappingHomogenization(1,ip,el)) - homogenization_RGC_postResults(c+3) = homogState(mappingHomogenization(2,ip,el))% & - state(3*nIntFaceTot+4,mappingHomogenization(1,ip,el)) + postResults(c+1:c+3) = dst%mismatch(1:3,of) c = c + 3_pInt case (penaltyenergy_ID) - homogenization_RGC_postResults(c+1) = homogState(mappingHomogenization(2,ip,el))% & - state(3*nIntFaceTot+5,mappingHomogenization(1,ip,el)) + postResults(c+1) = stt%penaltyEnergy(of) c = c + 1_pInt case (volumediscrepancy_ID) - homogenization_RGC_postResults(c+1) = homogState(mappingHomogenization(2,ip,el))% & - state(3*nIntFaceTot+6,mappingHomogenization(1,ip,el)) + postResults(c+1) = dst%volumeDiscrepancy(of) c = c + 1_pInt case (averagerelaxrate_ID) - homogenization_RGC_postResults(c+1) = homogState(mappingHomogenization(2,ip,el))% & - state(3*nIntFaceTot+7,mappingHomogenization(1,ip,el)) + postResults(c+1) = dst%relaxationrate_avg(of) c = c + 1_pInt case (maximumrelaxrate_ID) - homogenization_RGC_postResults(c+1) = homogState(mappingHomogenization(2,ip,el))% & - state(3*nIntFaceTot+8,mappingHomogenization(1,ip,el)) + postResults(c+1) = dst%relaxationrate_max(of) c = c + 1_pInt end select - enddo - -end function homogenization_RGC_postResults - - -!-------------------------------------------------------------------------------------------------- -!> @brief calculate stress-like penalty due to deformation mismatch -!-------------------------------------------------------------------------------------------------- -subroutine homogenization_RGC_stressPenalty(rPen,nMis,avgF,fDef,ip,el,homID) - use debug, only: & - debug_level, & - debug_homogenization,& - debug_levelExtensive, & - debug_e, & - debug_i - use mesh, only: & - mesh_element - use constitutive, only: & - constitutive_homogenizedC - use math, only: & - math_civita - use material, only: & - homogenization_maxNgrains,& - homogenization_Ngrains - use numerics, only: & - xSmoo_RGC - - implicit none - real(pReal), dimension (3,3,homogenization_maxNgrains), intent(out) :: rPen !< stress-like penalty - real(pReal), dimension (3,homogenization_maxNgrains), intent(out) :: nMis !< total amount of mismatch - real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: fDef !< deformation gradients - real(pReal), dimension (3,3), intent(in) :: avgF !< initial effective stretch tensor - integer(pInt), intent(in) :: ip,el - integer(pInt), dimension (4) :: intFace - integer(pInt), dimension (3) :: iGrain3,iGNghb3,nGDim - real(pReal), dimension (3,3) :: gDef,nDef - real(pReal), dimension (3) :: nVect,surfCorr - real(pReal), dimension (2) :: Gmoduli - integer(pInt) :: homID,iGrain,iGNghb,iFace,i,j,k,l - real(pReal) :: muGrain,muGNghb,nDefNorm,bgGrain,bgGNghb - - integer(pInt), parameter :: nFace = 6_pInt - real(pReal), parameter :: nDefToler = 1.0e-10_pReal - - nGDim = homogenization_RGC_Ngrains(1:3,homID) - rPen = 0.0_pReal - nMis = 0.0_pReal - -!-------------------------------------------------------------------------------------------------- -! get the correction factor the modulus of penalty stress representing the evolution of area of -! the interfaces due to deformations - surfCorr = homogenization_RGC_surfaceCorrection(avgF,ip,el) - -!-------------------------------------------------------------------------------------------------- -! debugging the surface correction factor - if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt & - .and. debug_e == el .and. debug_i == ip) then - !$OMP CRITICAL (write2out) - write(6,'(1x,a20,2(1x,i3))')'Correction factor: ',ip,el - write(6,'(1x,3(e11.4,1x))')(surfCorr(i), i = 1,3) - !$OMP END CRITICAL (write2out) - endif - -!-------------------------------------------------------------------------------------------------- -! computing the mismatch and penalty stress tensor of all grains - do iGrain = 1_pInt,homogenization_Ngrains(mesh_element(3,el)) - Gmoduli = homogenization_RGC_equivalentModuli(iGrain,ip,el) - muGrain = Gmoduli(1) ! collecting the equivalent shear modulus of grain - bgGrain = Gmoduli(2) ! and the lengthh of Burgers vector - iGrain3 = homogenization_RGC_grain1to3(iGrain,homID) ! get the grain ID in local 3-dimensional index (x,y,z)-position - -!* Looping over all six interfaces of each grain - do iFace = 1_pInt,nFace - intFace = homogenization_RGC_getInterface(iFace,iGrain3) ! get the 4-dimensional index of the interface in local numbering system of the grain - nVect = homogenization_RGC_interfaceNormal(intFace,ip,el) ! get the interface normal - iGNghb3 = iGrain3 ! identify the neighboring grain across the interface - iGNghb3(abs(intFace(1))) = iGNghb3(abs(intFace(1))) + int(real(intFace(1),pReal)/real(abs(intFace(1)),pReal),pInt) - if (iGNghb3(1) < 1) iGNghb3(1) = nGDim(1) ! with periodicity along e1 direction - if (iGNghb3(1) > nGDim(1)) iGNghb3(1) = 1_pInt - if (iGNghb3(2) < 1) iGNghb3(2) = nGDim(2) ! with periodicity along e2 direction - if (iGNghb3(2) > nGDim(2)) iGNghb3(2) = 1_pInt - if (iGNghb3(3) < 1) iGNghb3(3) = nGDim(3) ! with periodicity along e3 direction - if (iGNghb3(3) > nGDim(3)) iGNghb3(3) = 1_pInt - iGNghb = homogenization_RGC_grain3to1(iGNghb3,homID) ! get the ID of the neighboring grain - Gmoduli = homogenization_RGC_equivalentModuli(iGNghb,ip,el) ! collecting the shear modulus and Burgers vector of the neighbor - muGNghb = Gmoduli(1) - bgGNghb = Gmoduli(2) - gDef = 0.5_pReal*(fDef(1:3,1:3,iGNghb) - fDef(1:3,1:3,iGrain)) ! compute the difference/jump in deformation gradeint across the neighbor - -!-------------------------------------------------------------------------------------------------- -! compute the mismatch tensor of all interfaces - nDefNorm = 0.0_pReal - nDef = 0.0_pReal - do i = 1_pInt,3_pInt; do j = 1_pInt,3_pInt - do k = 1_pInt,3_pInt; do l = 1_pInt,3_pInt - nDef(i,j) = nDef(i,j) - nVect(k)*gDef(i,l)*math_civita(j,k,l) ! compute the interface mismatch tensor from the jump of deformation gradient - enddo; enddo - nDefNorm = nDefNorm + nDef(i,j)*nDef(i,j) ! compute the norm of the mismatch tensor - enddo; enddo - nDefNorm = max(nDefToler,sqrt(nDefNorm)) ! approximation to zero mismatch if mismatch is zero (singularity) - nMis(abs(intFace(1)),iGrain) = nMis(abs(intFace(1)),iGrain) + nDefNorm ! total amount of mismatch experienced by the grain (at all six interfaces) - -!-------------------------------------------------------------------------------------------------- -! debuggin the mismatch tensor - if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt & - .and. debug_e == el .and. debug_i == ip) then - !$OMP CRITICAL (write2out) - write(6,'(1x,a20,i2,1x,a20,1x,i3)')'Mismatch to face: ',intFace(1),'neighbor grain: ',iGNghb - do i = 1,3 - write(6,'(1x,3(e11.4,1x))')(nDef(i,j), j = 1,3) - enddo - write(6,'(1x,a20,e11.4)')'with magnitude: ',nDefNorm - !$OMP END CRITICAL (write2out) - endif - -!-------------------------------------------------------------------------------------------------- -! compute the stress penalty of all interfaces - do i = 1_pInt,3_pInt; do j = 1_pInt,3_pInt - do k = 1_pInt,3_pInt; do l = 1_pInt,3_pInt - rPen(i,j,iGrain) = rPen(i,j,iGrain) + 0.5_pReal*(muGrain*bgGrain + muGNghb*bgGNghb)*homogenization_RGC_xiAlpha(homID) & - *surfCorr(abs(intFace(1)))/homogenization_RGC_dAlpha(abs(intFace(1)),homID) & - *cosh(homogenization_RGC_ciAlpha(homID)*nDefNorm) & - *0.5_pReal*nVect(l)*nDef(i,k)/nDefNorm*math_civita(k,l,j) & - *tanh(nDefNorm/xSmoo_RGC) - enddo; enddo - enddo; enddo - enddo -!-------------------------------------------------------------------------------------------------- -! debugging the stress-like penalty - if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt & - .and. debug_e == el .and. debug_i == ip) then - !$OMP CRITICAL (write2out) - write(6,'(1x,a20,i2)')'Penalty of grain: ',iGrain - do i = 1,3 - write(6,'(1x,3(e11.4,1x))')(rPen(i,j,iGrain), j = 1,3) - enddo - !$OMP END CRITICAL (write2out) - endif - - enddo - -end subroutine homogenization_RGC_stressPenalty - - -!-------------------------------------------------------------------------------------------------- -!> @brief calculate stress-like penalty due to volume discrepancy -!-------------------------------------------------------------------------------------------------- -subroutine homogenization_RGC_volumePenalty(vPen,vDiscrep,fDef,fAvg,ip,el) - use debug, only: & - debug_level, & - debug_homogenization,& - debug_levelExtensive, & - debug_e, & - debug_i - use mesh, only: & - mesh_element - use math, only: & - math_det33, & - math_inv33 - use material, only: & - homogenization_maxNgrains,& - homogenization_Ngrains - use numerics, only: & - maxVolDiscr_RGC,& - volDiscrMod_RGC,& - volDiscrPow_RGC - - implicit none - real(pReal), dimension (3,3,homogenization_maxNgrains), intent(out) :: vPen ! stress-like penalty due to volume - real(pReal), intent(out) :: vDiscrep ! total volume discrepancy - real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: fDef ! deformation gradients - real(pReal), dimension (3,3), intent(in) :: fAvg ! overall deformation gradient - integer(pInt), intent(in) :: ip,& ! integration point - el - real(pReal), dimension (homogenization_maxNgrains) :: gVol - integer(pInt) :: iGrain,nGrain,i,j - - nGrain = homogenization_Ngrains(mesh_element(3,el)) - -!-------------------------------------------------------------------------------------------------- -! compute the volumes of grains and of cluster - vDiscrep = math_det33(fAvg) ! compute the volume of the cluster - do iGrain = 1_pInt,nGrain - gVol(iGrain) = math_det33(fDef(1:3,1:3,iGrain)) ! compute the volume of individual grains - vDiscrep = vDiscrep - gVol(iGrain)/real(nGrain,pReal) ! calculate the difference/dicrepancy between - ! the volume of the cluster and the the total volume of grains - enddo - -!-------------------------------------------------------------------------------------------------- -! calculate the stress and penalty due to volume discrepancy - vPen = 0.0_pReal - do iGrain = 1_pInt,nGrain - vPen(:,:,iGrain) = -1.0_pReal/real(nGrain,pReal)*volDiscrMod_RGC*volDiscrPow_RGC/maxVolDiscr_RGC* & - sign((abs(vDiscrep)/maxVolDiscr_RGC)**(volDiscrPow_RGC - 1.0),vDiscrep)* & - gVol(iGrain)*transpose(math_inv33(fDef(:,:,iGrain))) - -!-------------------------------------------------------------------------------------------------- -! debugging the stress-like penalty - if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt & - .and. debug_e == el .and. debug_i == ip) then - !$OMP CRITICAL (write2out) - write(6,'(1x,a30,i2)')'Volume penalty of grain: ',iGrain - do i = 1,3 - write(6,'(1x,3(e11.4,1x))')(vPen(i,j,iGrain), j = 1,3) - enddo - !$OMP END CRITICAL (write2out) - endif - enddo - -end subroutine homogenization_RGC_volumePenalty - - -!-------------------------------------------------------------------------------------------------- -!> @brief compute the correction factor accouted for surface evolution (area change) due to -! deformation -!-------------------------------------------------------------------------------------------------- -function homogenization_RGC_surfaceCorrection(avgF,ip,el) - use math, only: & - math_invert33, & - math_mul33x33 + enddo outputsLoop - implicit none - real(pReal), dimension(3) :: homogenization_RGC_surfaceCorrection - real(pReal), dimension(3,3), intent(in) :: avgF !< average F - integer(pInt), intent(in) :: ip,& !< integration point number - el !< element number - real(pReal), dimension(3,3) :: invC,avgC - real(pReal), dimension(3) :: nVect - real(pReal) :: detF - integer(pInt), dimension(4) :: intFace - integer(pInt) :: i,j,iBase - logical :: error - - avgC = math_mul33x33(transpose(avgF),avgF) - call math_invert33(avgC,invC,detF,error) - homogenization_RGC_surfaceCorrection = 0.0_pReal - do iBase = 1_pInt,3_pInt - intFace = [iBase,1_pInt,1_pInt,1_pInt] - nVect = homogenization_RGC_interfaceNormal(intFace,ip,el) ! get the normal of the interface - do i = 1_pInt,3_pInt; do j = 1_pInt,3_pInt - homogenization_RGC_surfaceCorrection(iBase) = & ! compute the component of (the inverse of) the stretch in the direction of the normal - homogenization_RGC_surfaceCorrection(iBase) + invC(i,j)*nVect(i)*nVect(j) - enddo; enddo - homogenization_RGC_surfaceCorrection(iBase) = & ! get the surface correction factor (area contraction/enlargement) - sqrt(homogenization_RGC_surfaceCorrection(iBase))*detF - enddo - -end function homogenization_RGC_surfaceCorrection - - -!-------------------------------------------------------------------------------------------------- -!> @brief compute the equivalent shear and bulk moduli from the elasticity tensor -!-------------------------------------------------------------------------------------------------- -function homogenization_RGC_equivalentModuli(grainID,ip,el) - use constitutive, only: & - constitutive_homogenizedC - - implicit none - integer(pInt), intent(in) :: & - grainID,& - ip, & !< integration point number - el !< element number - real(pReal), dimension (6,6) :: elasTens - real(pReal), dimension(2) :: homogenization_RGC_equivalentModuli - real(pReal) :: & - cEquiv_11, & - cEquiv_12, & - cEquiv_44 - - elasTens = constitutive_homogenizedC(grainID,ip,el) - -!-------------------------------------------------------------------------------------------------- -! compute the equivalent shear modulus after Turterltaub and Suiker, JMPS (2005) - cEquiv_11 = (elasTens(1,1) + elasTens(2,2) + elasTens(3,3))/3.0_pReal - cEquiv_12 = (elasTens(1,2) + elasTens(2,3) + elasTens(3,1) + & - elasTens(1,3) + elasTens(2,1) + elasTens(3,2))/6.0_pReal - cEquiv_44 = (elasTens(4,4) + elasTens(5,5) + elasTens(6,6))/3.0_pReal - homogenization_RGC_equivalentModuli(1) = 0.2_pReal*(cEquiv_11 - cEquiv_12) + 0.6_pReal*cEquiv_44 - -!-------------------------------------------------------------------------------------------------- -! obtain the length of Burgers vector (could be model dependend) - homogenization_RGC_equivalentModuli(2) = 2.5e-10_pReal - -end function homogenization_RGC_equivalentModuli + end associate + +end function homogenization_RGC_postResults !-------------------------------------------------------------------------------------------------- !> @brief collect relaxation vectors of an interface !-------------------------------------------------------------------------------------------------- -function homogenization_RGC_relaxationVector(intFace,homID, ip, el) - use material, only: & - homogState, & - mappingHomogenization +pure function relaxationVector(intFace,instance,of) implicit none - integer(pInt), intent(in) :: ip, el - real(pReal), dimension (3) :: homogenization_RGC_relaxationVector + integer(pInt), intent(in) :: instance,of + + real(pReal), dimension (3) :: relaxationVector integer(pInt), dimension (4), intent(in) :: intFace !< set of interface ID in 4D array (normal and position) - integer(pInt), dimension (3) :: nGDim integer(pInt) :: & - iNum, & - homID !< homogenization ID + iNum !-------------------------------------------------------------------------------------------------- ! collect the interface relaxation vector from the global state array - homogenization_RGC_relaxationVector = 0.0_pReal - nGDim = homogenization_RGC_Ngrains(1:3,homID) - iNum = homogenization_RGC_interface4to1(intFace,homID) ! identify the position of the interface in global state array - if (iNum > 0_pInt) homogenization_RGC_relaxationVector = homogState(mappingHomogenization(2,ip,el))% & - state((3*iNum-2):(3*iNum),mappingHomogenization(1,ip,el)) ! get the corresponding entries -end function homogenization_RGC_relaxationVector + iNum = interface4to1(intFace,param(instance)%Nconstituents) ! identify the position of the interface in global state array + if (iNum > 0_pInt) then + relaxationVector = state(instance)%relaxationVector((3*iNum-2):(3*iNum),of) + else + relaxationVector = 0.0_pReal + endif + +end function relaxationVector !-------------------------------------------------------------------------------------------------- !> @brief identify the normal of an interface !-------------------------------------------------------------------------------------------------- -function homogenization_RGC_interfaceNormal(intFace,ip,el) - use debug, only: & - debug_homogenization,& - debug_levelExtensive +pure function interfaceNormal(intFace,instance,of) use math, only: & math_mul33x3 implicit none - real(pReal), dimension (3) :: homogenization_RGC_interfaceNormal + real(pReal), dimension (3) :: interfaceNormal integer(pInt), dimension (4), intent(in) :: intFace !< interface ID in 4D array (normal and position) integer(pInt), intent(in) :: & - ip, & !< integration point number - el !< element number + instance, & + of integer(pInt) :: nPos !-------------------------------------------------------------------------------------------------- ! get the normal of the interface, identified from the value of intFace(1) - homogenization_RGC_interfaceNormal = 0.0_pReal + interfaceNormal = 0.0_pReal nPos = abs(intFace(1)) ! identify the position of the interface in global state array - homogenization_RGC_interfaceNormal(nPos) = real(intFace(1)/abs(intFace(1)),pReal) ! get the normal vector w.r.t. cluster axis + interfaceNormal(nPos) = real(intFace(1)/abs(intFace(1)),pReal) ! get the normal vector w.r.t. cluster axis - homogenization_RGC_interfaceNormal = & - math_mul33x3(homogenization_RGC_orientation(1:3,1:3,ip,el),homogenization_RGC_interfaceNormal) - ! map the normal vector into sample coordinate system (basis) + interfaceNormal = math_mul33x3(dependentState(instance)%orientation(1:3,1:3,of),interfaceNormal) ! map the normal vector into sample coordinate system (basis) -end function homogenization_RGC_interfaceNormal +end function interfaceNormal !-------------------------------------------------------------------------------------------------- !> @brief collect six faces of a grain in 4D (normal and position) !-------------------------------------------------------------------------------------------------- -function homogenization_RGC_getInterface(iFace,iGrain3) +pure function getInterface(iFace,iGrain3) implicit none - integer(pInt), dimension (4) :: homogenization_RGC_getInterface + integer(pInt), dimension (4) :: getInterface integer(pInt), dimension (3), intent(in) :: iGrain3 !< grain ID in 3D array integer(pInt), intent(in) :: iFace !< face index (1..6) mapped like (-e1,-e2,-e3,+e1,+e2,+e3) or iDir = (-1,-2,-3,1,2,3) integer(pInt) :: iDir !* Direction of interface normal iDir = (int(real(iFace-1_pInt,pReal)/2.0_pReal,pInt)+1_pInt)*(-1_pInt)**iFace - homogenization_RGC_getInterface(1) = iDir + getInterface(1) = iDir !-------------------------------------------------------------------------------------------------- ! identify the interface position by the direction of its normal - homogenization_RGC_getInterface(2:4) = iGrain3 - if (iDir < 0_pInt) & ! to have a correlation with coordinate/position in real space - homogenization_RGC_getInterface(1_pInt-iDir) = homogenization_RGC_getInterface(1_pInt-iDir)-1_pInt + getInterface(2:4) = iGrain3 + if (iDir < 0_pInt) getInterface(1_pInt-iDir) = getInterface(1_pInt-iDir)-1_pInt ! to have a correlation with coordinate/position in real space + +end function getInterface -end function homogenization_RGC_getInterface !-------------------------------------------------------------------------------------------------- !> @brief map grain ID from in 1D (global array) to in 3D (local position) !-------------------------------------------------------------------------------------------------- -function homogenization_RGC_grain1to3(grain1,homID) +pure function grain1to3(grain1,nGDim) implicit none - integer(pInt), dimension (3) :: homogenization_RGC_grain1to3 - integer(pInt), intent(in) :: & - grain1,& !< grain ID in 1D array - homID !< homogenization ID - integer(pInt), dimension (3) :: nGDim + integer(pInt), dimension(3) :: grain1to3 + integer(pInt), intent(in) :: grain1 !< grain ID in 1D array + integer(pInt), dimension(3), intent(in) :: nGDim -!-------------------------------------------------------------------------------------------------- -! get the grain position - nGDim = homogenization_RGC_Ngrains(1:3,homID) - homogenization_RGC_grain1to3(3) = 1_pInt+(grain1-1_pInt)/(nGDim(1)*nGDim(2)) - homogenization_RGC_grain1to3(2) = 1_pInt+mod((grain1-1_pInt)/nGDim(1),nGDim(2)) - homogenization_RGC_grain1to3(1) = 1_pInt+mod((grain1-1_pInt),nGDim(1)) + grain1to3 = 1_pInt + [mod((grain1-1_pInt),nGDim(1)), & + mod((grain1-1_pInt)/nGDim(1),nGDim(2)), & + (grain1-1_pInt)/(nGDim(1)*nGDim(2))] -end function homogenization_RGC_grain1to3 +end function grain1to3 !-------------------------------------------------------------------------------------------------- !> @brief map grain ID from in 3D (local position) to in 1D (global array) !-------------------------------------------------------------------------------------------------- -pure function homogenization_RGC_grain3to1(grain3,homID) +integer(pInt) pure function grain3to1(grain3,nGDim) implicit none - integer(pInt), dimension (3), intent(in) :: grain3 !< grain ID in 3D array (pos.x,pos.y,pos.z) - integer(pInt) :: homogenization_RGC_grain3to1 - integer(pInt), dimension (3) :: nGDim - integer(pInt), intent(in) :: homID ! homogenization ID + integer(pInt), dimension(3), intent(in) :: grain3 !< grain ID in 3D array (pos.x,pos.y,pos.z) + integer(pInt), dimension(3), intent(in) :: nGDim -!-------------------------------------------------------------------------------------------------- -! get the grain ID - nGDim = homogenization_RGC_Ngrains(1:3,homID) - homogenization_RGC_grain3to1 = grain3(1) + nGDim(1)*(grain3(2)-1_pInt) + nGDim(1)*nGDim(2)*(grain3(3)-1_pInt) + grain3to1 = grain3(1) & + + nGDim(1)*(grain3(2)-1_pInt) & + + nGDim(1)*nGDim(2)*(grain3(3)-1_pInt) -end function homogenization_RGC_grain3to1 +end function grain3to1 !-------------------------------------------------------------------------------------------------- !> @brief maps interface ID from 4D (normal and local position) into 1D (global array) !-------------------------------------------------------------------------------------------------- -integer(pInt) pure function homogenization_RGC_interface4to1(iFace4D, homID) +integer(pInt) pure function interface4to1(iFace4D, nGDim) implicit none - integer(pInt), dimension (4), intent(in) :: iFace4D !< interface ID in 4D array (n.dir,pos.x,pos.y,pos.z) - integer(pInt), dimension (3) :: nGDim,nIntFace - integer(pInt), intent(in) :: homID !< homogenization ID + integer(pInt), dimension(4), intent(in) :: iFace4D !< interface ID in 4D array (n.dir,pos.x,pos.y,pos.z) + integer(pInt), dimension(3), intent(in) :: nGDim - nGDim = homogenization_RGC_Ngrains(1:3,homID) -!-------------------------------------------------------------------------------------------------- -! compute the total number of interfaces, which ... - nIntFace(1) = (nGDim(1)-1_pInt)*nGDim(2)*nGDim(3) ! ... normal //e1 - nIntFace(2) = nGDim(1)*(nGDim(2)-1_pInt)*nGDim(3) ! ... normal //e2 - nIntFace(3) = nGDim(1)*nGDim(2)*(nGDim(3)-1_pInt) ! ... normal //e3 + select case(abs(iFace4D(1))) - homogenization_RGC_interface4to1 = -1_pInt - -!-------------------------------------------------------------------------------------------------- -! get the corresponding interface ID in 1D global array - if (abs(iFace4D(1)) == 1_pInt) then ! interface with normal //e1 - homogenization_RGC_interface4to1 = iFace4D(3) + nGDim(2)*(iFace4D(4)-1_pInt) & - + nGDim(2)*nGDim(3)*(iFace4D(2)-1_pInt) - if ((iFace4D(2) == 0_pInt) .or. (iFace4D(2) == nGDim(1))) homogenization_RGC_interface4to1 = 0_pInt - elseif (abs(iFace4D(1)) == 2_pInt) then ! interface with normal //e2 - homogenization_RGC_interface4to1 = iFace4D(4) + nGDim(3)*(iFace4D(2)-1_pInt) & - + nGDim(3)*nGDim(1)*(iFace4D(3)-1_pInt) + nIntFace(1) - if ((iFace4D(3) == 0_pInt) .or. (iFace4D(3) == nGDim(2))) homogenization_RGC_interface4to1 = 0_pInt - elseif (abs(iFace4D(1)) == 3_pInt) then ! interface with normal //e3 - homogenization_RGC_interface4to1 = iFace4D(2) + nGDim(1)*(iFace4D(3)-1_pInt) & - + nGDim(1)*nGDim(2)*(iFace4D(4)-1_pInt) + nIntFace(1) + nIntFace(2) - if ((iFace4D(4) == 0_pInt) .or. (iFace4D(4) == nGDim(3))) homogenization_RGC_interface4to1 = 0_pInt - endif + case(1_pInt) + if ((iFace4D(2) == 0_pInt) .or. (iFace4D(2) == nGDim(1))) then + interface4to1 = 0_pInt + else + interface4to1 = iFace4D(3) + nGDim(2)*(iFace4D(4)-1_pInt) & + + nGDim(2)*nGDim(3)*(iFace4D(2)-1_pInt) + endif -end function homogenization_RGC_interface4to1 + case(2_pInt) + if ((iFace4D(3) == 0_pInt) .or. (iFace4D(3) == nGDim(2))) then + interface4to1 = 0_pInt + else + interface4to1 = iFace4D(4) + nGDim(3)*(iFace4D(2)-1_pInt) & + + nGDim(3)*nGDim(1)*(iFace4D(3)-1_pInt) & + + (nGDim(1)-1_pInt)*nGDim(2)*nGDim(3) ! total number of interfaces normal //e1 + endif + + case(3_pInt) + if ((iFace4D(4) == 0_pInt) .or. (iFace4D(4) == nGDim(3))) then + interface4to1 = 0_pInt + else + interface4to1 = iFace4D(2) + nGDim(1)*(iFace4D(3)-1_pInt) & + + nGDim(1)*nGDim(2)*(iFace4D(4)-1_pInt) & + + (nGDim(1)-1_pInt)*nGDim(2)*nGDim(3) & ! total number of interfaces normal //e1 + + nGDim(1)*(nGDim(2)-1_pInt)*nGDim(3) ! total number of interfaces normal //e2 + endif + + case default + interface4to1 = -1_pInt + + end select + +end function interface4to1 !-------------------------------------------------------------------------------------------------- !> @brief maps interface ID from 1D (global array) into 4D (normal and local position) !-------------------------------------------------------------------------------------------------- -function homogenization_RGC_interface1to4(iFace1D, homID) +pure function interface1to4(iFace1D, nGDim) implicit none - integer(pInt), dimension (4) :: homogenization_RGC_interface1to4 - integer(pInt), intent(in) :: iFace1D !< interface ID in 1D array - integer(pInt), dimension (3) :: nGDim,nIntFace - integer(pInt), intent(in) :: homID !< homogenization ID - - nGDim = homogenization_RGC_Ngrains(:,homID) + integer(pInt), dimension(4) :: interface1to4 + integer(pInt), intent(in) :: iFace1D !< interface ID in 1D array + integer(pInt), dimension(3), intent(in) :: nGDim + integer(pInt), dimension(3) :: nIntFace !-------------------------------------------------------------------------------------------------- ! compute the total number of interfaces, which ... - nIntFace(1) = (nGDim(1)-1_pInt)*nGDim(2)*nGDim(3) ! ... normal //e1 - nIntFace(2) = nGDim(1)*(nGDim(2)-1_pInt)*nGDim(3) ! ... normal //e2 - nIntFace(3) = nGDim(1)*nGDim(2)*(nGDim(3)-1_pInt) ! ... normal //e3 + nIntFace = [(nGDim(1)-1_pInt)*nGDim(2)*nGDim(3), & ! ... normal //e1 + nGDim(1)*(nGDim(2)-1_pInt)*nGDim(3), & ! ... normal //e2 + nGDim(1)*nGDim(2)*(nGDim(3)-1_pInt)] ! ... normal //e3 !-------------------------------------------------------------------------------------------------- ! get the corresponding interface ID in 4D (normal and local position) if (iFace1D > 0 .and. iFace1D <= nIntFace(1)) then ! interface with normal //e1 - homogenization_RGC_interface1to4(1) = 1_pInt - homogenization_RGC_interface1to4(3) = mod((iFace1D-1_pInt),nGDim(2))+1_pInt - homogenization_RGC_interface1to4(4) = mod(& + interface1to4(1) = 1_pInt + interface1to4(3) = mod((iFace1D-1_pInt),nGDim(2))+1_pInt + interface1to4(4) = mod(& int(& real(iFace1D-1_pInt,pReal)/& real(nGDim(2),pReal)& ,pInt)& ,nGDim(3))+1_pInt - homogenization_RGC_interface1to4(2) = int(& + interface1to4(2) = int(& real(iFace1D-1_pInt,pReal)/& real(nGDim(2),pReal)/& real(nGDim(3),pReal)& ,pInt)+1_pInt elseif (iFace1D > nIntFace(1) .and. iFace1D <= (nIntFace(2) + nIntFace(1))) then ! interface with normal //e2 - homogenization_RGC_interface1to4(1) = 2_pInt - homogenization_RGC_interface1to4(4) = mod((iFace1D-nIntFace(1)-1_pInt),nGDim(3))+1_pInt - homogenization_RGC_interface1to4(2) = mod(& + interface1to4(1) = 2_pInt + interface1to4(4) = mod((iFace1D-nIntFace(1)-1_pInt),nGDim(3))+1_pInt + interface1to4(2) = mod(& int(& real(iFace1D-nIntFace(1)-1_pInt,pReal)/& real(nGDim(3),pReal)& ,pInt)& ,nGDim(1))+1_pInt - homogenization_RGC_interface1to4(3) = int(& + interface1to4(3) = int(& real(iFace1D-nIntFace(1)-1_pInt,pReal)/& real(nGDim(3),pReal)/& real(nGDim(1),pReal)& ,pInt)+1_pInt elseif (iFace1D > nIntFace(2) + nIntFace(1) .and. iFace1D <= (nIntFace(3) + nIntFace(2) + nIntFace(1))) then ! interface with normal //e3 - homogenization_RGC_interface1to4(1) = 3_pInt - homogenization_RGC_interface1to4(2) = mod((iFace1D-nIntFace(2)-nIntFace(1)-1_pInt),nGDim(1))+1_pInt - homogenization_RGC_interface1to4(3) = mod(& + interface1to4(1) = 3_pInt + interface1to4(2) = mod((iFace1D-nIntFace(2)-nIntFace(1)-1_pInt),nGDim(1))+1_pInt + interface1to4(3) = mod(& int(& real(iFace1D-nIntFace(2)-nIntFace(1)-1_pInt,pReal)/& real(nGDim(1),pReal)& ,pInt)& ,nGDim(2))+1_pInt - homogenization_RGC_interface1to4(4) = int(& + interface1to4(4) = int(& real(iFace1D-nIntFace(2)-nIntFace(1)-1_pInt,pReal)/& real(nGDim(1),pReal)/& real(nGDim(2),pReal)& ,pInt)+1_pInt endif -end function homogenization_RGC_interface1to4 +end function interface1to4 -!-------------------------------------------------------------------------------------------------- -!> @brief calculating the grain deformation gradient (the same with -! homogenization_RGC_partionDeformation, but used only for perturbation scheme) -!-------------------------------------------------------------------------------------------------- -subroutine homogenization_RGC_grainDeformation(F, avgF, ip, el) - use mesh, only: & - mesh_element - use material, only: & - homogenization_maxNgrains,& - homogenization_Ngrains, & - homogenization_typeInstance - - implicit none - real(pReal), dimension (3,3,homogenization_maxNgrains), intent(out) :: F !< partioned F per grain - real(pReal), dimension (3,3), intent(in) :: avgF !< - integer(pInt), intent(in) :: & - el, & !< element number - ip !< integration point number - real(pReal), dimension (3) :: aVect,nVect - integer(pInt), dimension (4) :: intFace - integer(pInt), dimension (3) :: iGrain3 - integer(pInt) :: homID, iGrain,iFace,i,j - integer(pInt), parameter :: nFace = 6_pInt - -!-------------------------------------------------------------------------------------------------- -! compute the deformation gradient of individual grains due to relaxations - homID = homogenization_typeInstance(mesh_element(3,el)) - F = 0.0_pReal - do iGrain = 1_pInt,homogenization_Ngrains(mesh_element(3,el)) - iGrain3 = homogenization_RGC_grain1to3(iGrain,homID) - do iFace = 1_pInt,nFace - intFace = homogenization_RGC_getInterface(iFace,iGrain3) - aVect = homogenization_RGC_relaxationVector(intFace,homID, ip, el) - nVect = homogenization_RGC_interfaceNormal(intFace,ip,el) - forall (i=1_pInt:3_pInt,j=1_pInt:3_pInt) & - F(i,j,iGrain) = F(i,j,iGrain) + aVect(i)*nVect(j) ! effective relaxations - enddo - F(1:3,1:3,iGrain) = F(1:3,1:3,iGrain) + avgF ! relaxed deformation gradient - enddo - -end subroutine homogenization_RGC_grainDeformation - end module homogenization_RGC diff --git a/src/homogenization_isostrain.f90 b/src/homogenization_isostrain.f90 index 24aedf75f..42c0c9287 100644 --- a/src/homogenization_isostrain.f90 +++ b/src/homogenization_isostrain.f90 @@ -1,4 +1,5 @@ !-------------------------------------------------------------------------------------------------- +!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH !> @brief Isostrain (full constraint Taylor assuption) homogenization scheme @@ -6,220 +7,119 @@ module homogenization_isostrain use prec, only: & pInt - + implicit none private - integer(pInt), dimension(:), allocatable, public, protected :: & - homogenization_isostrain_sizePostResults - integer(pInt), dimension(:,:), allocatable, target, public :: & - homogenization_isostrain_sizePostResult - - character(len=64), dimension(:,:), allocatable, target, public :: & - homogenization_isostrain_output !< name of each post result output - integer(pInt), dimension(:), allocatable, target, public :: & - homogenization_isostrain_Noutput !< number of outputs per homog instance - integer(pInt), dimension(:), allocatable, private :: & - homogenization_isostrain_Ngrains enum, bind(c) - enumerator :: undefined_ID, & - nconstituents_ID, & - ipcoords_ID, & - avgdefgrad_ID, & - avgfirstpiola_ID + enumerator :: & + parallel_ID, & + average_ID end enum - enum, bind(c) - enumerator :: parallel_ID, & - average_ID - end enum - integer(kind(undefined_ID)), dimension(:,:), allocatable, private :: & - homogenization_isostrain_outputID !< ID of each post result output - integer(kind(average_ID)), dimension(:), allocatable, private :: & - homogenization_isostrain_mapping !< mapping type + type, private :: tParameters !< container type for internal constitutive parameters + integer(pInt) :: & + Nconstituents + integer(kind(average_ID)) :: & + mapping + end type + + type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance) public :: & homogenization_isostrain_init, & homogenization_isostrain_partitionDeformation, & - homogenization_isostrain_averageStressAndItsTangent, & - homogenization_isostrain_postResults + homogenization_isostrain_averageStressAndItsTangent contains !-------------------------------------------------------------------------------------------------- !> @brief allocates all neccessary fields, reads information from material configuration file !-------------------------------------------------------------------------------------------------- -subroutine homogenization_isostrain_init(fileUnit) +subroutine homogenization_isostrain_init() #if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & compiler_options #endif - use prec, only: & - pReal use debug, only: & debug_HOMOGENIZATION, & debug_level, & debug_levelBasic - use IO - use material - use config + use IO, only: & + IO_timeStamp, & + IO_error + use material, only: & + homogenization_type, & + material_homog, & + homogState, & + HOMOGENIZATION_ISOSTRAIN_ID, & + HOMOGENIZATION_ISOSTRAIN_LABEL, & + homogenization_typeInstance + use config, only: & + config_homogenization implicit none - integer(pInt), intent(in) :: fileUnit - integer(pInt), allocatable, dimension(:) :: chunkPos integer(pInt) :: & - section = 0_pInt, i, mySize, o - integer :: & - maxNinstance, & - homog, & - instance - integer :: & - NofMyHomog ! no pInt (stores a system dependen value from 'count' + Ninstance, & + h, & + NofMyHomog character(len=65536) :: & - tag = '', & - line = '' - + tag = '' + write(6,'(/,a)') ' <<<+- homogenization_'//HOMOGENIZATION_ISOSTRAIN_label//' init -+>>>' write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" - maxNinstance = count(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID) - if (maxNinstance == 0) return - + Ninstance = int(count(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID),pInt) if (iand(debug_level(debug_HOMOGENIZATION),debug_levelBasic) /= 0_pInt) & - write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance - allocate(homogenization_isostrain_sizePostResults(maxNinstance), source=0_pInt) - allocate(homogenization_isostrain_sizePostResult(maxval(homogenization_Noutput),maxNinstance), & - source=0_pInt) - allocate(homogenization_isostrain_Noutput(maxNinstance), source=0_pInt) - allocate(homogenization_isostrain_Ngrains(maxNinstance), source=0_pInt) - allocate(homogenization_isostrain_mapping(maxNinstance), source=average_ID) - allocate(homogenization_isostrain_output(maxval(homogenization_Noutput),maxNinstance)) - homogenization_isostrain_output = '' - allocate(homogenization_isostrain_outputID(maxval(homogenization_Noutput),maxNinstance), & - source=undefined_ID) + write(6,'(a16,1x,i5,/)') '# instances:',Ninstance + + allocate(param(Ninstance)) ! one container of parameters per instance + + do h = 1_pInt, size(homogenization_type) + if (homogenization_type(h) /= HOMOGENIZATION_ISOSTRAIN_ID) cycle + + associate(prm => param(homogenization_typeInstance(h)),& + config => config_homogenization(h)) + + prm%Nconstituents = config_homogenization(h)%getInt('nconstituents') + tag = 'sum' + select case(trim(config%getString('mapping',defaultVal = tag))) + case ('sum') + prm%mapping = parallel_ID + case ('avg') + prm%mapping = average_ID + case default + call IO_error(211_pInt,ext_msg=trim(tag)//' ('//HOMOGENIZATION_isostrain_label//')') + end select + + NofMyHomog = count(material_homog == h) + homogState(h)%sizeState = 0_pInt + homogState(h)%sizePostResults = 0_pInt + allocate(homogState(h)%state0 (0_pInt,NofMyHomog)) + allocate(homogState(h)%subState0(0_pInt,NofMyHomog)) + allocate(homogState(h)%state (0_pInt,NofMyHomog)) + + end associate - rewind(fileUnit) - do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= material_partHomogenization)! wind forward to - line = IO_read(fileUnit) enddo - parsingFile: do while (trim(line) /= IO_EOF) ! read through sections of homogenization part - line = IO_read(fileUnit) - if (IO_isBlank(line)) cycle ! skip empty lines - if (IO_getTag(line,'<','>') /= '') then ! stop at next part - line = IO_read(fileUnit, .true.) ! reset IO_read - exit - endif - if (IO_getTag(line,'[',']') /= '') then ! next section - section = section + 1_pInt - cycle - endif - if (section > 0_pInt ) then ! do not short-circuit here (.and. with next if-statement). It's not safe in Fortran - if (homogenization_type(section) == HOMOGENIZATION_ISOSTRAIN_ID) then ! one of my sections - i = homogenization_typeInstance(section) ! which instance of my type is present homogenization - chunkPos = IO_stringPos(line) - tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key - select case(tag) - case ('(output)') - select case(IO_lc(IO_stringValue(line,chunkPos,2_pInt))) - case('nconstituents','ngrains') - homogenization_isostrain_Noutput(i) = homogenization_isostrain_Noutput(i) + 1_pInt - homogenization_isostrain_outputID(homogenization_isostrain_Noutput(i),i) = nconstituents_ID - homogenization_isostrain_output(homogenization_isostrain_Noutput(i),i) = & - IO_lc(IO_stringValue(line,chunkPos,2_pInt)) - case('ipcoords') - homogenization_isostrain_Noutput(i) = homogenization_isostrain_Noutput(i) + 1_pInt - homogenization_isostrain_outputID(homogenization_isostrain_Noutput(i),i) = ipcoords_ID - homogenization_isostrain_output(homogenization_isostrain_Noutput(i),i) = & - IO_lc(IO_stringValue(line,chunkPos,2_pInt)) - case('avgdefgrad','avgf') - homogenization_isostrain_Noutput(i) = homogenization_isostrain_Noutput(i) + 1_pInt - homogenization_isostrain_outputID(homogenization_isostrain_Noutput(i),i) = avgdefgrad_ID - homogenization_isostrain_output(homogenization_isostrain_Noutput(i),i) = & - IO_lc(IO_stringValue(line,chunkPos,2_pInt)) - case('avgp','avgfirstpiola','avg1stpiola') - homogenization_isostrain_Noutput(i) = homogenization_isostrain_Noutput(i) + 1_pInt - homogenization_isostrain_outputID(homogenization_isostrain_Noutput(i),i) = avgfirstpiola_ID - homogenization_isostrain_output(homogenization_isostrain_Noutput(i),i) = & - IO_lc(IO_stringValue(line,chunkPos,2_pInt)) - - end select - case ('nconstituents','ngrains') - homogenization_isostrain_Ngrains(i) = IO_intValue(line,chunkPos,2_pInt) - case ('mapping') - select case(IO_lc(IO_stringValue(line,chunkPos,2_pInt))) - case ('parallel','sum') - homogenization_isostrain_mapping(i) = parallel_ID - case ('average','mean','avg') - homogenization_isostrain_mapping(i) = average_ID - case default - call IO_error(211_pInt,ext_msg=trim(tag)//' ('//HOMOGENIZATION_isostrain_label//')') - end select - - end select - endif - endif - enddo parsingFile - - initializeInstances: do homog = 1_pInt, material_Nhomogenization - myHomog: if (homogenization_type(homog) == HOMOGENIZATION_ISOSTRAIN_ID) then - NofMyHomog = count(material_homog == homog) - instance = homogenization_typeInstance(homog) - -! * Determine size of postResults array - outputsLoop: do o = 1_pInt, homogenization_isostrain_Noutput(instance) - select case(homogenization_isostrain_outputID(o,instance)) - case(nconstituents_ID) - mySize = 1_pInt - case(ipcoords_ID) - mySize = 3_pInt - case(avgdefgrad_ID, avgfirstpiola_ID) - mySize = 9_pInt - case default - mySize = 0_pInt - end select - - outputFound: if (mySize > 0_pInt) then - homogenization_isostrain_sizePostResult(o,instance) = mySize - homogenization_isostrain_sizePostResults(instance) = & - homogenization_isostrain_sizePostResults(instance) + mySize - endif outputFound - enddo outputsLoop - -! allocate state arrays - homogState(homog)%sizeState = 0_pInt - homogState(homog)%sizePostResults = homogenization_isostrain_sizePostResults(instance) - allocate(homogState(homog)%state0 (0_pInt,NofMyHomog), source=0.0_pReal) - allocate(homogState(homog)%subState0(0_pInt,NofMyHomog), source=0.0_pReal) - allocate(homogState(homog)%state (0_pInt,NofMyHomog), source=0.0_pReal) - - endif myHomog - enddo initializeInstances - end subroutine homogenization_isostrain_init !-------------------------------------------------------------------------------------------------- !> @brief partitions the deformation gradient onto the constituents !-------------------------------------------------------------------------------------------------- -subroutine homogenization_isostrain_partitionDeformation(F,avgF,el) +subroutine homogenization_isostrain_partitionDeformation(F,avgF) use prec, only: & pReal - use mesh, only: & - mesh_element - use material, only: & - homogenization_maxNgrains, & - homogenization_Ngrains implicit none - real(pReal), dimension (3,3,homogenization_maxNgrains), intent(out) :: F !< partioned def grad per grain - real(pReal), dimension (3,3), intent(in) :: avgF !< my average def grad - integer(pInt), intent(in) :: & - el !< element number - F=0.0_pReal - F(1:3,1:3,1:homogenization_Ngrains(mesh_element(3,el)))= & - spread(avgF,3,homogenization_Ngrains(mesh_element(3,el))) + real(pReal), dimension (:,:,:), intent(out) :: F !< partitioned deformation gradient + + real(pReal), dimension (3,3), intent(in) :: avgF !< average deformation gradient at material point + + F = spread(avgF,3,size(F,3)) end subroutine homogenization_isostrain_partitionDeformation @@ -227,90 +127,31 @@ end subroutine homogenization_isostrain_partitionDeformation !-------------------------------------------------------------------------------------------------- !> @brief derive average stress and stiffness from constituent quantities !-------------------------------------------------------------------------------------------------- -subroutine homogenization_isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,el) +subroutine homogenization_isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,instance) use prec, only: & pReal - use mesh, only: & - mesh_element - use material, only: & - homogenization_maxNgrains, & - homogenization_Ngrains, & - homogenization_typeInstance implicit none - real(pReal), dimension (3,3), intent(out) :: avgP !< average stress at material point - real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF !< average stiffness at material point - real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: P !< array of current grain stresses - real(pReal), dimension (3,3,3,3,homogenization_maxNgrains), intent(in) :: dPdF !< array of current grain stiffnesses - integer(pInt), intent(in) :: el !< element number - integer(pInt) :: & - homID, & - Ngrains + real(pReal), dimension (3,3), intent(out) :: avgP !< average stress at material point + real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF !< average stiffness at material point - homID = homogenization_typeInstance(mesh_element(3,el)) - Ngrains = homogenization_Ngrains(mesh_element(3,el)) + real(pReal), dimension (:,:,:), intent(in) :: P !< partitioned stresses + real(pReal), dimension (:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses + integer(pInt), intent(in) :: instance - select case (homogenization_isostrain_mapping(homID)) + associate(prm => param(instance)) + + select case (prm%mapping) case (parallel_ID) avgP = sum(P,3) dAvgPdAvgF = sum(dPdF,5) case (average_ID) - avgP = sum(P,3) /real(Ngrains,pReal) - dAvgPdAvgF = sum(dPdF,5)/real(Ngrains,pReal) + avgP = sum(P,3) /real(prm%Nconstituents,pReal) + dAvgPdAvgF = sum(dPdF,5)/real(prm%Nconstituents,pReal) end select + + end associate end subroutine homogenization_isostrain_averageStressAndItsTangent - -!-------------------------------------------------------------------------------------------------- -!> @brief return array of homogenization results for post file inclusion -!-------------------------------------------------------------------------------------------------- -pure function homogenization_isostrain_postResults(ip,el,avgP,avgF) - use prec, only: & - pReal - use mesh, only: & - mesh_element, & - mesh_ipCoordinates - use material, only: & - homogenization_typeInstance, & - homogenization_Noutput - - implicit none - integer(pInt), intent(in) :: & - ip, & !< integration point number - el !< element number - real(pReal), dimension(3,3), intent(in) :: & - avgP, & !< average stress at material point - avgF !< average deformation gradient at material point - real(pReal), dimension(homogenization_isostrain_sizePostResults & - (homogenization_typeInstance(mesh_element(3,el)))) :: & - homogenization_isostrain_postResults - - integer(pInt) :: & - homID, & - o, c - - c = 0_pInt - homID = homogenization_typeInstance(mesh_element(3,el)) - homogenization_isostrain_postResults = 0.0_pReal - - do o = 1_pInt,homogenization_Noutput(mesh_element(3,el)) - select case(homogenization_isostrain_outputID(o,homID)) - case (nconstituents_ID) - homogenization_isostrain_postResults(c+1_pInt) = real(homogenization_isostrain_Ngrains(homID),pReal) - c = c + 1_pInt - case (avgdefgrad_ID) - homogenization_isostrain_postResults(c+1_pInt:c+9_pInt) = reshape(transpose(avgF),[9]) - c = c + 9_pInt - case (avgfirstpiola_ID) - homogenization_isostrain_postResults(c+1_pInt:c+9_pInt) = reshape(transpose(avgP),[9]) - c = c + 9_pInt - case (ipcoords_ID) - homogenization_isostrain_postResults(c+1_pInt:c+3_pInt) = mesh_ipCoordinates(1:3,ip,el) ! current ip coordinates - c = c + 3_pInt - end select - enddo - -end function homogenization_isostrain_postResults - end module homogenization_isostrain diff --git a/src/homogenization_none.f90 b/src/homogenization_none.f90 index c33aabe89..04ea55abe 100644 --- a/src/homogenization_none.f90 +++ b/src/homogenization_none.f90 @@ -2,7 +2,7 @@ !> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH -!> @brief dummy homogenization homogenization scheme +!> @brief dummy homogenization homogenization scheme for 1 constituent per material point !-------------------------------------------------------------------------------------------------- module homogenization_none @@ -24,35 +24,46 @@ subroutine homogenization_none_init() compiler_options #endif use prec, only: & - pReal, & - pInt + pInt + use debug, only: & + debug_HOMOGENIZATION, & + debug_level, & + debug_levelBasic use IO, only: & IO_timeStamp - use material - use config - + + use material, only: & + homogenization_type, & + material_homog, & + homogState, & + HOMOGENIZATION_NONE_LABEL, & + HOMOGENIZATION_NONE_ID + implicit none integer(pInt) :: & - homog, & + Ninstance, & + h, & NofMyHomog write(6,'(/,a)') ' <<<+- homogenization_'//HOMOGENIZATION_NONE_label//' init -+>>>' write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" - initializeInstances: do homog = 1_pInt, material_Nhomogenization + Ninstance = int(count(homogenization_type == HOMOGENIZATION_NONE_ID),pInt) + if (iand(debug_level(debug_HOMOGENIZATION),debug_levelBasic) /= 0_pInt) & + write(6,'(a16,1x,i5,/)') '# instances:',Ninstance + + do h = 1_pInt, size(homogenization_type) + if (homogenization_type(h) /= HOMOGENIZATION_NONE_ID) cycle - myhomog: if (homogenization_type(homog) == HOMOGENIZATION_none_ID) then - NofMyHomog = count(material_homog == homog) - homogState(homog)%sizeState = 0_pInt - homogState(homog)%sizePostResults = 0_pInt - allocate(homogState(homog)%state0 (0_pInt,NofMyHomog), source=0.0_pReal) - allocate(homogState(homog)%subState0(0_pInt,NofMyHomog), source=0.0_pReal) - allocate(homogState(homog)%state (0_pInt,NofMyHomog), source=0.0_pReal) - - endif myhomog - enddo initializeInstances + NofMyHomog = count(material_homog == h) + homogState(h)%sizeState = 0_pInt + homogState(h)%sizePostResults = 0_pInt + allocate(homogState(h)%state0 (0_pInt,NofMyHomog)) + allocate(homogState(h)%subState0(0_pInt,NofMyHomog)) + allocate(homogState(h)%state (0_pInt,NofMyHomog)) + enddo end subroutine homogenization_none_init diff --git a/src/math.f90 b/src/math.f90 index 12f56707a..28c7175e3 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -1998,6 +1998,7 @@ end function math_symmetricEulers !-------------------------------------------------------------------------------------------------- !> @brief eigenvalues and eigenvectors of symmetric matrix m +! ToDo: has wrong oder of arguments !-------------------------------------------------------------------------------------------------- subroutine math_eigenValuesVectorsSym(m,values,vectors,error) @@ -2021,9 +2022,10 @@ end subroutine math_eigenValuesVectorsSym !-------------------------------------------------------------------------------------------------- !> @brief eigenvalues and eigenvectors of symmetric 33 matrix m using an analytical expression !> and the general LAPACK powered version for arbritrary sized matrices as fallback -!> @author Joachim Kopp, Max–Planck–Institut für Kernphysik, Heidelberg (Copyright (C) 2006) +!> @author Joachim Kopp, Max-Planck-Institut für Kernphysik, Heidelberg (Copyright (C) 2006) !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @details See http://arxiv.org/abs/physics/0610206 (DSYEVH3) +! ToDo: has wrong oder of arguments !-------------------------------------------------------------------------------------------------- subroutine math_eigenValuesVectorsSym33(m,values,vectors)