From b651f334fe3afb875db3fa989d708b7fdadecebf Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 11 Oct 2013 16:01:53 +0000 Subject: [PATCH] set useful default values to remove dummy functions in individual homogenization schemes (stateInit, averageTemperature, and updateState (isostrain only)) --- code/homogenization.f90 | 161 +++++++++++------------- code/homogenization_RGC.f90 | 65 +++------- code/homogenization_isostrain.f90 | 200 ++++++++++++------------------ 3 files changed, 164 insertions(+), 262 deletions(-) diff --git a/code/homogenization.f90 b/code/homogenization.f90 index 64abf2690..0acfe630d 100644 --- a/code/homogenization.f90 +++ b/code/homogenization.f90 @@ -214,7 +214,7 @@ subroutine homogenization_init(Temperature) allocate(homogenization_state0(i,e)%p(homogenization_isostrain_sizeState(myInstance))) allocate(homogenization_subState0(i,e)%p(homogenization_isostrain_sizeState(myInstance))) allocate(homogenization_state(i,e)%p(homogenization_isostrain_sizeState(myInstance))) - homogenization_state0(i,e)%p = homogenization_isostrain_stateInit(myInstance) + homogenization_state0(i,e)%p = 0.0_pReal homogenization_sizeState(i,e) = homogenization_isostrain_sizeState(myInstance) endif homogenization_sizePostResults(i,e) = homogenization_isostrain_sizePostResults(myInstance) @@ -223,7 +223,7 @@ subroutine homogenization_init(Temperature) allocate(homogenization_state0(i,e)%p(homogenization_RGC_sizeState(myInstance))) allocate(homogenization_subState0(i,e)%p(homogenization_RGC_sizeState(myInstance))) allocate(homogenization_state(i,e)%p(homogenization_RGC_sizeState(myInstance))) - homogenization_state0(i,e)%p = homogenization_RGC_stateInit(myInstance) + homogenization_state0(i,e)%p = 0.0_pReal homogenization_sizeState(i,e) = homogenization_RGC_sizeState(myInstance) endif homogenization_sizePostResults(i,e) = homogenization_RGC_sizePostResults(myInstance) @@ -575,7 +575,7 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt) 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 homogenization_averageTemperature(i,e) + materialpoint_Temperature(i,e) = homogenization_averageTemperature(i,e) enddo IpLooping4 enddo elementLooping4 !$OMP END PARALLEL DO @@ -652,7 +652,7 @@ end subroutine materialpoint_postResults !-------------------------------------------------------------------------------------------------- !> @brief partition material point def grad onto constituents !-------------------------------------------------------------------------------------------------- -subroutine homogenization_partitionDeformation(i,e) +subroutine homogenization_partitionDeformation(ip,el) use mesh, only: & mesh_element use material, only: & @@ -670,25 +670,25 @@ subroutine homogenization_partitionDeformation(i,e) implicit none integer(pInt), intent(in) :: & - i, & !< integration point - e !< element number + ip, & !< integration point + el !< element number - chosenHomogenization: select case(homogenization_type(mesh_element(3,e))) + chosenHomogenization: select case(homogenization_type(mesh_element(3,el))) case (homogenization_isostrain_label) chosenHomogenization call homogenization_isostrain_partitionDeformation(& - crystallite_partionedF(1:3,1:3,1:homogenization_maxNgrains,i,e), & - crystallite_partionedF0(1:3,1:3,1:homogenization_maxNgrains,i,e),& - materialpoint_subF(1:3,1:3,i,e),& - homogenization_state(i,e), & - i, & - e) + crystallite_partionedF(1:3,1:3,1:homogenization_maxNgrains,ip,el), & + crystallite_partionedF0(1:3,1:3,1:homogenization_maxNgrains,ip,el),& + materialpoint_subF(1:3,1:3,ip,el),& + homogenization_state(ip,el), & + ip, & + el) case (homogenization_RGC_label) chosenHomogenization - call homogenization_RGC_partitionDeformation(crystallite_partionedF(1:3,1:3,1:homogenization_maxNgrains,i,e), & - crystallite_partionedF0(1:3,1:3,1:homogenization_maxNgrains,i,e),& - materialpoint_subF(1:3,1:3,i,e),& - homogenization_state(i,e), & - i, & - e) + call homogenization_RGC_partitionDeformation(crystallite_partionedF(1:3,1:3,1:homogenization_maxNgrains,ip,el), & + crystallite_partionedF0(1:3,1:3,1:homogenization_maxNgrains,ip,el),& + materialpoint_subF(1:3,1:3,ip,el),& + homogenization_state(ip,el), & + ip, & + el) end select chosenHomogenization end subroutine homogenization_partitionDeformation @@ -698,7 +698,7 @@ end subroutine homogenization_partitionDeformation !> @brief update the internal state of the homogenization scheme and tell whether "done" and !> "happy" with result !-------------------------------------------------------------------------------------------------- -function homogenization_updateState(i,e) +function homogenization_updateState(ip,el) use mesh, only: & mesh_element use material, only: & @@ -709,39 +709,32 @@ function homogenization_updateState(i,e) crystallite_dPdF, & crystallite_partionedF,& crystallite_partionedF0 - use homogenization_isostrain, only: & - homogenization_isostrain_updateState, & - homogenization_isostrain_label use homogenization_RGC, only: & homogenization_RGC_updateState, & homogenization_RGC_label implicit none integer(pInt), intent(in) :: & - i, & !< integration point - e !< element number + ip, & !< integration point + el !< element number logical, dimension(2) :: homogenization_updateState - chosenHomogenization: select case(homogenization_type(mesh_element(3,e))) - case (homogenization_isostrain_label) chosenHomogenization - homogenization_updateState = & - homogenization_isostrain_updateState( homogenization_state(i,e), & - crystallite_P(1:3,1:3,1:homogenization_maxNgrains,i,e), & - crystallite_dPdF(1:3,1:3,1:3,1:3,1:homogenization_maxNgrains,i,e), & - i, & - e) + chosenHomogenization: select case(homogenization_type(mesh_element(3,el))) + case (homogenization_RGC_label) chosenHomogenization homogenization_updateState = & - homogenization_RGC_updateState( homogenization_state(i,e), & - homogenization_subState0(i,e), & - crystallite_P(1:3,1:3,1:homogenization_maxNgrains,i,e), & - crystallite_partionedF(1:3,1:3,1:homogenization_maxNgrains,i,e), & - crystallite_partionedF0(1:3,1:3,1:homogenization_maxNgrains,i,e),& - materialpoint_subF(1:3,1:3,i,e),& - materialpoint_subdt(i,e), & - crystallite_dPdF(1:3,1:3,1:3,1:3,1:homogenization_maxNgrains,i,e), & - i, & - e) + homogenization_RGC_updateState( homogenization_state(ip,el), & + homogenization_subState0(ip,el), & + 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),& + 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), & + ip, & + e)l + case default chosenHomogenization + homogenization_updateState = .true. end select chosenHomogenization end function homogenization_updateState @@ -750,7 +743,7 @@ end function homogenization_updateState !-------------------------------------------------------------------------------------------------- !> @brief derive average stress and stiffness from constituent quantities !-------------------------------------------------------------------------------------------------- -subroutine homogenization_averageStressAndItsTangent(i,e) +subroutine homogenization_averageStressAndItsTangent(ip,el) use mesh, only: & mesh_element use material, only: & @@ -767,69 +760,61 @@ subroutine homogenization_averageStressAndItsTangent(i,e) implicit none integer(pInt), intent(in) :: & - i, & !< integration point - e !< element number + ip, & !< integration point + el !< element number - chosenHomogenization: select case(homogenization_type(mesh_element(3,e))) + chosenHomogenization: select case(homogenization_type(mesh_element(3,el))) case (homogenization_isostrain_label) chosenHomogenization - call homogenization_isostrain_averageStressAndItsTangent(materialpoint_P(1:3,1:3,i,e), & - materialpoint_dPdF(1:3,1:3,1:3,1:3,i,e),& - crystallite_P(1:3,1:3,1:homogenization_maxNgrains,i,e), & - crystallite_dPdF(1:3,1:3,1:3,1:3,1:homogenization_maxNgrains,i,e), & - i, & - e) + 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), & + ip, & + el) case (homogenization_RGC_label) chosenHomogenization - call homogenization_RGC_averageStressAndItsTangent( materialpoint_P(1:3,1:3,i,e), & - materialpoint_dPdF(1:3,1:3,1:3,1:3,i,e),& - crystallite_P(1:3,1:3,1:homogenization_maxNgrains,i,e), & - crystallite_dPdF(1:3,1:3,1:3,1:3,1:homogenization_maxNgrains,i,e), & - i, & - e) + 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), & + ip, & + el) end select chosenHomogenization end subroutine homogenization_averageStressAndItsTangent !-------------------------------------------------------------------------------------------------- -!> @brief derive average temperature from constituent quantities +!> @brief derive average temperature from constituent quantities (does not depend on choosen +!! homogenization scheme) !-------------------------------------------------------------------------------------------------- -subroutine homogenization_averageTemperature(i,e) +real(pReal) function homogenization_averageTemperature(ip,el) use mesh, only: & mesh_element use material, only: & - homogenization_type, & - homogenization_maxNgrains + homogenization_Ngrains use crystallite, only: & crystallite_Temperature - use homogenization_isostrain, only: & - homogenization_isostrain_averageTemperature, & - homogenization_isostrain_label - use homogenization_RGC, only: & - homogenization_RGC_averageTemperature, & - homogenization_RGC_label implicit none integer(pInt), intent(in) :: & - i, & !< integration point - e !< element number - - chosenHomogenization: select case(homogenization_type(mesh_element(3,e))) - case (homogenization_isostrain_label) chosenHomogenization - materialpoint_Temperature(i,e) = & - homogenization_isostrain_averageTemperature(crystallite_Temperature(1:homogenization_maxNgrains,i,e), i, e) - case (homogenization_RGC_label) chosenHomogenization - materialpoint_Temperature(i,e) = & - homogenization_RGC_averageTemperature(crystallite_Temperature(1:homogenization_maxNgrains,i,e), i, e) - end select chosenHomogenization + ip, & !< integration point number + el !< element number + integer(pInt) :: & + Ngrains -end subroutine homogenization_averageTemperature +!-------------------------------------------------------------------------------------------------- +! computing the average temperature + Ngrains = homogenization_Ngrains(mesh_element(3,el)) + homogenization_averageTemperature= sum(crystallite_Temperature(1:Ngrains,ip,el))/real(Ngrains,pReal) + +end function homogenization_averageTemperature !-------------------------------------------------------------------------------------------------- !> @brief return array of homogenization results for post file inclusion. call only, !> if homogenization_sizePostResults(i,e) > 0 !! !-------------------------------------------------------------------------------------------------- -function homogenization_postResults(i,e) +function homogenization_postResults(ip,el) use mesh, only: & mesh_element use material, only: & @@ -843,16 +828,16 @@ function homogenization_postResults(i,e) implicit none integer(pInt), intent(in) :: & - i, & !< integration point - e !< element number - real(pReal), dimension(homogenization_sizePostResults(i,e)) :: homogenization_postResults + ip, & !< integration point + el !< element number + real(pReal), dimension(homogenization_sizePostResults(ip,el)) :: homogenization_postResults homogenization_postResults = 0.0_pReal - chosenHomogenization: select case (homogenization_type(mesh_element(3,e))) + chosenHomogenization: select case (homogenization_type(mesh_element(3,el))) case (homogenization_isostrain_label) chosenHomogenization - homogenization_postResults = homogenization_isostrain_postResults(homogenization_state(i,e),i,e) + homogenization_postResults = homogenization_isostrain_postResults(homogenization_state(ip,el),ip,el) case (homogenization_RGC_label) chosenHomogenization - homogenization_postResults = homogenization_RGC_postResults(homogenization_state(i,e),i,e) + homogenization_postResults = homogenization_RGC_postResults(homogenization_state(ip,el),ip,el) end select chosenHomogenization end function homogenization_postResults diff --git a/code/homogenization_RGC.f90 b/code/homogenization_RGC.f90 index e08429602..28fbc820d 100644 --- a/code/homogenization_RGC.f90 +++ b/code/homogenization_RGC.f90 @@ -52,14 +52,11 @@ module homogenization_RGC real(pReal), dimension(:), allocatable, private :: & homogenization_RGC_xiAlpha, & homogenization_RGC_ciAlpha - public :: & homogenization_RGC_init, & - homogenization_RGC_stateInit, & homogenization_RGC_partitionDeformation, & homogenization_RGC_averageStressAndItsTangent, & homogenization_RGC_updateState, & - homogenization_RGC_averageTemperature, & homogenization_RGC_postResults private :: & homogenization_RGC_stressPenalty, & @@ -105,8 +102,8 @@ subroutine homogenization_RGC_init(myFile) implicit none integer(pInt), intent(in) :: myFile !< file pointer to material configuration - integer(pInt), parameter :: maxNchunks = 4_pInt - integer(pInt), dimension(1_pInt+2_pInt*maxNchunks) :: positions + integer(pInt), parameter :: MAXNCHUNKS = 4_pInt + integer(pInt), dimension(1_pInt+2_pInt*MAXNCHUNKS) :: positions integer(pInt) ::section=0_pInt, maxNinstance, i,j,e, output=-1_pInt, mySize, myInstance character(len=65536) :: tag character(len=65536) :: line = '' @@ -150,7 +147,7 @@ subroutine homogenization_RGC_init(myFile) if (section > 0_pInt ) then ! do not short-circuit here (.and. with next if-statement). It's not safe in Fortran if (trim(homogenization_type(section)) == HOMOGENIZATION_RGC_label) then ! one of my sections i = homogenization_typeInstance(section) ! which instance of my type is present homogenization - positions = IO_stringPos(line,maxNchunks) + positions = IO_stringPos(line,MAXNCHUNKS) tag = IO_lc(IO_stringValue(line,positions,1_pInt)) ! extract key select case(tag) case ('(output)') @@ -249,20 +246,6 @@ subroutine homogenization_RGC_init(myFile) end subroutine homogenization_RGC_init -!-------------------------------------------------------------------------------------------------- -!> @brief sets the initial homogenization state -!-------------------------------------------------------------------------------------------------- -function homogenization_RGC_stateInit(myInstance) - - implicit none - integer(pInt), intent(in) :: myInstance - real(pReal), dimension(homogenization_RGC_sizeState(myInstance)) :: homogenization_RGC_stateInit - - homogenization_RGC_stateInit = 0.0_pReal - -end function homogenization_RGC_stateInit - - !-------------------------------------------------------------------------------------------------- !> @brief partitions the deformation gradient onto the constituents !-------------------------------------------------------------------------------------------------- @@ -877,27 +860,7 @@ subroutine homogenization_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF, avgP = sum(P,3)/real(Ngrains,pReal) dAvgPdAvgF = sum(dPdF,5)/real(Ngrains,pReal) -end subroutine - - -!-------------------------------------------------------------------------------------------------- -!> @brief derive average temperature from constituent quantities -!-------------------------------------------------------------------------------------------------- -real(pReal) pure function homogenization_RGC_averageTemperature(Temperature,ip,el) - use mesh, only: mesh_element - use material, only: homogenization_maxNgrains, homogenization_Ngrains - - implicit none - real(pReal), dimension (homogenization_maxNgrains), intent(in) :: Temperature - integer(pInt), intent(in) :: ip,el - integer(pInt) :: Ngrains - -!-------------------------------------------------------------------------------------------------- -! computing the average temperature - Ngrains = homogenization_Ngrains(mesh_element(3,el)) - homogenization_RGC_averageTemperature = sum(Temperature(1:Ngrains))/real(Ngrains,pReal) - -end function homogenization_RGC_averageTemperature +end subroutine homogenization_RGC_averageStressAndItsTangent !-------------------------------------------------------------------------------------------------- @@ -990,9 +953,9 @@ subroutine homogenization_RGC_stressPenalty(rPen,nMis,avgF,fDef,ip,el,homID) 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) :: 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 @@ -1016,7 +979,7 @@ subroutine homogenization_RGC_stressPenalty(rPen,nMis,avgF,fDef,ip,el,homID) !$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) @@ -1129,8 +1092,8 @@ subroutine homogenization_RGC_volumePenalty(vPen,vDiscrep,fDef,fAvg,ip,el, homID integer(pInt), intent(in) :: ip,& ! integration point el real(pReal), dimension (homogenization_maxNgrains) :: gVol - integer(pInt) homID,iGrain,nGrain,i,j -! + integer(pInt) :: homID,iGrain,nGrain,i,j + nGrain = homogenization_Ngrains(mesh_element(3,el)) !-------------------------------------------------------------------------------------------------- @@ -1182,10 +1145,10 @@ function homogenization_RGC_surfaceCorrection(avgF,ip,el) el !< element number real(pReal), dimension(3,3) :: invC,avgC real(pReal), dimension(3) :: nVect - real(pReal) detF + real(pReal) :: detF integer(pInt), dimension(4) :: intFace - integer(pInt) i,j,iBase - logical error + integer(pInt) :: i,j,iBase + logical :: error avgC = math_mul33x33(transpose(avgF),avgF) call math_invert33(avgC,invC,detF,error) @@ -1219,7 +1182,7 @@ function homogenization_RGC_equivalentModuli(grainID,ip,el) el !< element number real(pReal), dimension (6,6) :: elasTens real(pReal), dimension(2) :: homogenization_RGC_equivalentModuli - real(pReal) cEquiv_11,cEquiv_12,cEquiv_44 + real(pReal) :: cEquiv_11,cEquiv_12,cEquiv_44 elasTens = constitutive_homogenizedC(grainID,ip,el) diff --git a/code/homogenization_isostrain.f90 b/code/homogenization_isostrain.f90 index 914e2e3dc..f863a9f20 100644 --- a/code/homogenization_isostrain.f90 +++ b/code/homogenization_isostrain.f90 @@ -29,15 +29,15 @@ module homogenization_isostrain implicit none private - character (len=*), parameter, public :: & + character (len=*), parameter, public :: & HOMOGENIZATION_ISOSTRAIN_label = 'isostrain' - integer(pInt), dimension(:), allocatable, public :: & + integer(pInt), dimension(:), allocatable, public, protected :: & homogenization_isostrain_sizeState, & homogenization_isostrain_sizePostResults - integer(pInt), dimension(:,:), allocatable, target, public :: & + integer(pInt), dimension(:,:), allocatable, target, public :: & homogenization_isostrain_sizePostResult - character(len=64), dimension(:,:), allocatable, target, public :: & + character(len=64), dimension(:,:), allocatable, target, public :: & homogenization_isostrain_output !< name of each post result output character(len=64), dimension(:), allocatable, private :: & homogenization_isostrain_mapping @@ -46,11 +46,8 @@ module homogenization_isostrain public :: & homogenization_isostrain_init, & - homogenization_isostrain_stateInit, & homogenization_isostrain_partitionDeformation, & - homogenization_isostrain_updateState, & homogenization_isostrain_averageStressAndItsTangent, & - homogenization_isostrain_averageTemperature, & homogenization_isostrain_postResults contains @@ -60,14 +57,20 @@ contains !-------------------------------------------------------------------------------------------------- subroutine homogenization_isostrain_init(myFile) use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) - use math, only: math_Mandel3333to66, math_Voigt66to3333 + use math, only: & + math_Mandel3333to66, & + math_Voigt66to3333 use IO use material - integer(pInt), intent(in) :: myFile - integer(pInt), parameter :: maxNchunks = 2_pInt - integer(pInt), dimension(1_pInt+2_pInt*maxNchunks) :: positions - integer(pInt) section, i, j, output, mySize - integer :: maxNinstance, k ! no pInt (stores a system dependen value from 'count' + + implicit none + integer(pInt), intent(in) :: myFile + integer(pInt), parameter :: MAXNCHUNKS = 2_pInt + integer(pInt), dimension(1_pInt+2_pInt*MAXNCHUNKS) :: positions + integer(pInt) :: & + section, i, j, output, mySize + integer :: & + maxNinstance, k ! no pInt (stores a system dependen value from 'count' character(len=65536) :: & tag = '', & line = '' ! to start initialized @@ -80,14 +83,18 @@ subroutine homogenization_isostrain_init(myFile) maxNinstance = count(homogenization_type == HOMOGENIZATION_ISOSTRAIN_label) if (maxNinstance == 0) return - allocate(homogenization_isostrain_sizeState(maxNinstance)) ; homogenization_isostrain_sizeState = 0_pInt - allocate(homogenization_isostrain_sizePostResults(maxNinstance)); homogenization_isostrain_sizePostResults = 0_pInt - allocate(homogenization_isostrain_sizePostResult(maxval(homogenization_Noutput), & - maxNinstance)); homogenization_isostrain_sizePostResult = 0_pInt - allocate(homogenization_isostrain_Ngrains(maxNinstance)); homogenization_isostrain_Ngrains = 0_pInt - allocate(homogenization_isostrain_mapping(maxNinstance)); homogenization_isostrain_mapping = 'avg' - allocate(homogenization_isostrain_output(maxval(homogenization_Noutput), & - maxNinstance)) ; homogenization_isostrain_output = '' + allocate(homogenization_isostrain_sizeState(maxNinstance)) + homogenization_isostrain_sizeState = 0_pInt + allocate(homogenization_isostrain_sizePostResults(maxNinstance)) + homogenization_isostrain_sizePostResults = 0_pInt + allocate(homogenization_isostrain_sizePostResult(maxval(homogenization_Noutput),maxNinstance)) + homogenization_isostrain_sizePostResult = 0_pInt + allocate(homogenization_isostrain_Ngrains(maxNinstance)) + homogenization_isostrain_Ngrains = 0_pInt + allocate(homogenization_isostrain_mapping(maxNinstance)) + homogenization_isostrain_mapping = 'avg' + allocate(homogenization_isostrain_output(maxval(homogenization_Noutput),maxNinstance)) + homogenization_isostrain_output = '' rewind(myFile) section = 0_pInt @@ -107,7 +114,7 @@ subroutine homogenization_isostrain_init(myFile) if (section > 0_pInt ) then ! do not short-circuit here (.and. with next if-statement). It's not safe in Fortran if (trim(homogenization_type(section)) == HOMOGENIZATION_ISOSTRAIN_label) then ! one of my sections i = homogenization_typeInstance(section) ! which instance of my type is present homogenization - positions = IO_stringPos(line,maxNchunks) + positions = IO_stringPos(line,MAXNCHUNKS) tag = IO_lc(IO_stringValue(line,positions,1_pInt)) ! extract key select case(tag) case ('(output)') @@ -144,92 +151,60 @@ subroutine homogenization_isostrain_init(myFile) end subroutine homogenization_isostrain_init -!-------------------------------------------------------------------------------------------------- -!> @brief sets the initial homogenization stated -!-------------------------------------------------------------------------------------------------- -function homogenization_isostrain_stateInit(myInstance) - use prec, only: & - pReal - - implicit none - integer(pInt), intent(in) :: myInstance - real(pReal), dimension(homogenization_isostrain_sizeState(myInstance)) :: & - homogenization_isostrain_stateInit - - homogenization_isostrain_stateInit = 0.0_pReal - -end function homogenization_isostrain_stateInit - - !-------------------------------------------------------------------------------------------------- !> @brief partitions the deformation gradient onto the constituents !-------------------------------------------------------------------------------------------------- -subroutine homogenization_isostrain_partitionDeformation(F,F0,avgF,state,i,e) - use prec, only: pReal,p_vec - use mesh, only: mesh_element - use material, only: homogenization_maxNgrains,homogenization_Ngrains +subroutine homogenization_isostrain_partitionDeformation(F,F0,avgF,state,ip,el) + use prec, only: & + pReal, & + p_vec + 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,homogenization_maxNgrains), intent(in) :: F0 ! initial partioned def grad per grain - real(pReal), dimension (3,3), intent(in) :: avgF ! my average def grad - type(p_vec), intent(in) :: state ! my state - integer(pInt), intent(in) :: & - i, & !< integration point number - e !< element number + real(pReal), dimension (3,3,homogenization_maxNgrains), intent(out) :: F ! partioned def grad per grain + real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: F0 ! initial partioned def grad per grain + real(pReal), dimension (3,3), intent(in) :: avgF ! my average def grad + type(p_vec), intent(in) :: state ! my state + integer(pInt), intent(in) :: & + ip, & !< integration point number + el !< element number - F = spread(avgF,3,homogenization_Ngrains(mesh_element(3,e))) + F = spread(avgF,3,homogenization_Ngrains(mesh_element(3,el))) end subroutine homogenization_isostrain_partitionDeformation -!-------------------------------------------------------------------------------------------------- -!> @brief update the internal state of the homogenization scheme and tell whether "done" and -! "happy" with result -!-------------------------------------------------------------------------------------------------- -function homogenization_isostrain_updateState(state,P,dPdF,i,e) - use prec, only: & - pReal,& - p_vec - use material, only: & - homogenization_maxNgrains - - implicit none - type(p_vec), intent(inout) :: state !< my state - 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) :: & - i, & !< integration point number - e !< element number - logical, dimension(2) :: homogenization_isostrain_updateState - - homogenization_isostrain_updateState = .true. ! homogenization at material point converged (done and happy) - -end function homogenization_isostrain_updateState - - !-------------------------------------------------------------------------------------------------- !> @brief derive average stress and stiffness from constituent quantities !-------------------------------------------------------------------------------------------------- -subroutine homogenization_isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,i,e) +subroutine homogenization_isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,ip,el) use prec, only: & pReal use mesh, only: & mesh_element - use material, only: homogenization_maxNgrains, homogenization_Ngrains, homogenization_typeInstance + 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) :: & - i, & !< integration point number - e !< 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 + 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) :: & + ip, & !< integration point number + el !< element number + integer(pInt) :: & + homID, & + Ngrains - homID = homogenization_typeInstance(mesh_element(3,e)) - Ngrains = homogenization_Ngrains(mesh_element(3,e)) + homID = homogenization_typeInstance(mesh_element(3,el)) + Ngrains = homogenization_Ngrains(mesh_element(3,el)) select case (homogenization_isostrain_mapping(homID)) case ('parallel','sum') @@ -246,35 +221,10 @@ subroutine homogenization_isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P end subroutine homogenization_isostrain_averageStressAndItsTangent -!-------------------------------------------------------------------------------------------------- -!> @brief derive average temperature from constituent quantities -!-------------------------------------------------------------------------------------------------- -real(pReal) pure function homogenization_isostrain_averageTemperature(Temperature,i,e) - use prec, only: & - pReal - use mesh, only: & - mesh_element - use material, only: & - homogenization_maxNgrains, & - homogenization_Ngrains - - implicit none - real(pReal), dimension (homogenization_maxNgrains), intent(in) :: Temperature - integer(pInt), intent(in) :: & - i, & !< integration point number - e !< element number - integer(pInt) :: Ngrains - - Ngrains = homogenization_Ngrains(mesh_element(3,e)) - homogenization_isostrain_averageTemperature = sum(Temperature(1:Ngrains))/real(Ngrains,pReal) - -end function homogenization_isostrain_averageTemperature - - !-------------------------------------------------------------------------------------------------- !> @brief return array of homogenization results for post file inclusion !-------------------------------------------------------------------------------------------------- -pure function homogenization_isostrain_postResults(state,i,e) +pure function homogenization_isostrain_postResults(state,ip,el) use prec, only: & pReal,& p_vec @@ -285,19 +235,23 @@ pure function homogenization_isostrain_postResults(state,i,e) homogenization_Noutput implicit none - type(p_vec), intent(in) :: state + type(p_vec), intent(in) :: state integer(pInt), intent(in) :: & - i, & !< integration point number - e !< element number - integer(pInt) :: homID,o,c - real(pReal), dimension(homogenization_isostrain_sizePostResults& - (homogenization_typeInstance(mesh_element(3,e)))) :: homogenization_isostrain_postResults - + ip, & !< integration point number + el !< element number + 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,e)) + homID = homogenization_typeInstance(mesh_element(3,el)) homogenization_isostrain_postResults = 0.0_pReal - do o = 1_pInt,homogenization_Noutput(mesh_element(3,e)) + do o = 1_pInt,homogenization_Noutput(mesh_element(3,el)) select case(homogenization_isostrain_output(o,homID)) case ('ngrains') homogenization_isostrain_postResults(c+1_pInt) = real(homogenization_isostrain_Ngrains(homID),pReal)