diff --git a/code/homogenization_isostrain.f90 b/code/homogenization_isostrain.f90 index 5671c504e..a377c1d2f 100644 --- a/code/homogenization_isostrain.f90 +++ b/code/homogenization_isostrain.f90 @@ -16,56 +16,49 @@ ! You should have received a copy of the GNU General Public License ! along with DAMASK. If not, see . ! -!############################################################## -!* $Id$ -!***************************************************** -!* Module: HOMOGENIZATION_ISOSTRAIN * -!***************************************************** -!* contains: * -!***************************************************** - -! [isostrain] -! type isostrain -! Ngrains 6 -! (output) Ngrains - +!-------------------------------------------------------------------------------------------------- +! $Id$ +!-------------------------------------------------------------------------------------------------- +!> @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 +!-------------------------------------------------------------------------------------------------- module homogenization_isostrain - - use prec, only: pInt + use prec, only: & + pInt implicit none - character (len=*), parameter :: & + private + character (len=*), parameter, public :: & homogenization_isostrain_label = 'isostrain' - integer(pInt),dimension(:), allocatable :: & + integer(pInt), dimension(:), allocatable, public :: & homogenization_isostrain_sizeState, & - homogenization_isostrain_Ngrains, & homogenization_isostrain_sizePostResults - - integer(pInt), dimension(:,:), allocatable, target :: & + 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 - character(len=64), dimension(:,:), allocatable, target :: & - homogenization_isostrain_output ! name of each post result output + integer(pInt), dimension(:), allocatable, private :: & + homogenization_isostrain_Ngrains + public :: & + homogenization_isostrain_init, & + homogenization_isostrain_stateInit, & + homogenization_isostrain_partitionDeformation, & + homogenization_isostrain_updateState, & + homogenization_isostrain_averageStressAndItsTangent, & + homogenization_isostrain_averageTemperature, & + homogenization_isostrain_postResults contains -!**************************************** -!* - homogenization_isostrain_init -!* - homogenization_isostrain_stateInit -!* - homogenization_isostrain_deformationPartititon -!* - homogenization_isostrain_stateUpdate -!* - homogenization_isostrain_averageStressAndItsTangent -!* - homogenization_isostrain_postResults -!**************************************** - -!************************************** -!* Module initialization * -!************************************** -subroutine homogenization_isostrain_init(myFile) ! file pointer to material configuration - use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) - use prec, only: pInt +!-------------------------------------------------------------------------------------------------- +!> @brief allocates all neccessary fields, reads information from material configuration file +!-------------------------------------------------------------------------------------------------- +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 IO use material @@ -73,7 +66,7 @@ subroutine homogenization_isostrain_init(myFile) ! fil 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' + integer :: maxNinstance, k ! no pInt (stores a system dependen value from 'count' character(len=64) :: tag character(len=1024) :: line = '' ! to start initialized @@ -98,22 +91,22 @@ subroutine homogenization_isostrain_init(myFile) ! fil rewind(myFile) section = 0_pInt - do while (IO_lc(IO_getTag(line,'<','>')) /= material_partHomogenization) ! wind forward to + do while (IO_lc(IO_getTag(line,'<','>')) /= material_partHomogenization) ! wind forward to read(myFile,'(a1024)',END=100) line enddo - do ! read thru sections of phase part + do ! read thru sections of phase part read(myFile,'(a1024)',END=100) line - if (IO_isBlank(line)) cycle ! skip empty lines - if (IO_getTag(line,'<','>') /= '') exit ! stop at next part - if (IO_getTag(line,'[',']') /= '') then ! next section + if (IO_isBlank(line)) cycle ! skip empty lines + if (IO_getTag(line,'<','>') /= '') exit ! stop at next part + if (IO_getTag(line,'[',']') /= '') then ! next section section = section + 1_pInt - output = 0_pInt ! reset output counter + output = 0_pInt ! reset output counter endif - if (section > 0 .and. homogenization_type(section) == homogenization_isostrain_label) then ! one of my sections - i = homogenization_typeInstance(section) ! which instance of my type is present homogenization + if (section > 0 .and. 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) - tag = IO_lc(IO_stringValue(line,positions,1_pInt)) ! extract key + tag = IO_lc(IO_stringValue(line,positions,1_pInt)) ! extract key select case(tag) case ('(output)') output = output + 1_pInt @@ -124,10 +117,7 @@ subroutine homogenization_isostrain_init(myFile) ! fil endif enddo -100 do k = 1,maxNinstance ! sanity checks - enddo - - do k = 1,maxNinstance +100 do k = 1,maxNinstance homogenization_isostrain_sizeState(i) = 0_pInt do j = 1_pInt,maxval(homogenization_Noutput) @@ -138,7 +128,7 @@ subroutine homogenization_isostrain_init(myFile) ! fil mySize = 0_pInt end select - if (mySize > 0_pInt) then ! any meaningful output found + if (mySize > 0_pInt) then ! any meaningful output found homogenization_isostrain_sizePostResult(j,i) = mySize homogenization_isostrain_sizePostResults(i) = & homogenization_isostrain_sizePostResults(i) + mySize @@ -149,11 +139,12 @@ subroutine homogenization_isostrain_init(myFile) ! fil end subroutine homogenization_isostrain_init -!********************************************************************* -!* initial homogenization state * -!********************************************************************* +!-------------------------------------------------------------------------------------------------- +!> @brief sets the initial homogenization stated +!-------------------------------------------------------------------------------------------------- function homogenization_isostrain_stateInit(myInstance) - use prec, only: pReal,pInt + use prec, only: & + pReal implicit none integer(pInt), intent(in) :: myInstance @@ -165,151 +156,132 @@ function homogenization_isostrain_stateInit(myInstance) end function homogenization_isostrain_stateInit -!******************************************************************** -! partition material point def grad onto constituents -!******************************************************************** -subroutine homogenization_isostrain_partitionDeformation(& - F, & ! partioned def grad per grain -! - F0, & ! initial partioned def grad per grain - avgF, & ! my average def grad - state, & ! my state - ip, & ! my integration point - el & ! my element - ) - use prec, only: pReal,pInt,p_vec +!-------------------------------------------------------------------------------------------------- +!> @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 implicit none - real(pReal), dimension (3,3,homogenization_maxNgrains), intent(out) :: F - real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: F0 - real(pReal), dimension (3,3), intent(in) :: avgF - type(p_vec), intent(in) :: state - integer(pInt), intent(in) :: ip,el - integer(pInt) i - -! homID = homogenization_typeInstance(mesh_element(3,el)) - forall (i = 1_pInt:homogenization_Ngrains(mesh_element(3,el))) & - F(1:3,1:3,i) = avgF + 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 + + F = spread(avgF,3,homogenization_Ngrains(mesh_element(3,e))) end subroutine homogenization_isostrain_partitionDeformation -!******************************************************************** -! update the internal state of the homogenization scheme -! and tell whether "done" and "happy" with result -!******************************************************************** -function homogenization_isostrain_updateState(& - state, & ! my state -! - P, & ! array of current grain stresses - dPdF, & ! array of current grain stiffnesses - ip, & ! my integration point - el & ! my element - ) - - use prec, only: pReal,pInt,p_vec - use material, only: homogenization_maxNgrains +!-------------------------------------------------------------------------------------------------- +!> @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 - -!* Definition of variables - type(p_vec), intent(inout) :: state - real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: P - real(pReal), dimension (3,3,3,3,homogenization_maxNgrains), intent(in) :: dPdF - integer(pInt), intent(in) :: ip,el -! integer(pInt) homID + 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 -! homID = homogenization_typeInstance(mesh_element(3,el)) - homogenization_isostrain_updateState = .true. ! homogenization at material point converged (done and happy) + homogenization_isostrain_updateState = .true. ! homogenization at material point converged (done and happy) end function homogenization_isostrain_updateState -!******************************************************************** -! derive average stress and stiffness from constituent quantities -!******************************************************************** -subroutine homogenization_isostrain_averageStressAndItsTangent(& - avgP, & ! average stress at material point - dAvgPdAvgF, & ! average stiffness at material point -! - P, & ! array of current grain stresses - dPdF, & ! array of current grain stiffnesses - ip, & ! my integration point - el & ! my element - ) - - use prec, only: pReal,pInt,p_vec - use mesh, only: mesh_element +!-------------------------------------------------------------------------------------------------- +!> @brief derive average stress and stiffness from constituent quantities +!-------------------------------------------------------------------------------------------------- +subroutine homogenization_isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,i,e) + use prec, only: & + pReal + use mesh, only: & + mesh_element use material, only: homogenization_maxNgrains, homogenization_Ngrains implicit none - real(pReal), dimension (3,3), intent(out) :: avgP - real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF - real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: P - real(pReal), dimension (3,3,3,3,homogenization_maxNgrains), intent(in) :: dPdF - integer(pInt), intent(in) :: ip,el - integer(pInt) 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) :: & + i, & !< integration point number + e !< element number + integer(pInt) :: Ngrains -! homID = homogenization_typeInstance(mesh_element(3,el)) - Ngrains = homogenization_Ngrains(mesh_element(3,el)) + Ngrains = homogenization_Ngrains(mesh_element(3,e)) avgP = sum(P,3)/real(Ngrains,pReal) dAvgPdAvgF = sum(dPdF,5)/real(Ngrains,pReal) end subroutine homogenization_isostrain_averageStressAndItsTangent -!******************************************************************** -! derive average stress and stiffness from constituent quantities -!******************************************************************** -pure function homogenization_isostrain_averageTemperature(& - Temperature, & ! temperature - ip, & ! my integration point - el & ! my element - ) - - use prec, only: pReal,pInt,p_vec - use mesh, only: mesh_element - use material, only: homogenization_maxNgrains, homogenization_Ngrains +!-------------------------------------------------------------------------------------------------- +!> @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) :: ip,el - real(pReal) homogenization_isostrain_averageTemperature - integer(pInt) Ngrains + integer(pInt), intent(in) :: & + i, & !< integration point number + e !< element number + integer(pInt) :: Ngrains -! homID = homogenization_typeInstance(mesh_element(3,el)) - Ngrains = homogenization_Ngrains(mesh_element(3,el)) + Ngrains = homogenization_Ngrains(mesh_element(3,e)) homogenization_isostrain_averageTemperature = sum(Temperature(1:Ngrains))/real(Ngrains,pReal) end function homogenization_isostrain_averageTemperature -!******************************************************************** -! return array of homogenization results for post file inclusion -!******************************************************************** -pure function homogenization_isostrain_postResults(& - state, & ! my state - ip, & ! my integration point - el & ! my element - ) - - use prec, only: pReal,pInt,p_vec - use mesh, only: mesh_element - use material, only: homogenization_typeInstance,homogenization_Noutput +!-------------------------------------------------------------------------------------------------- +!> @brief return array of homogenization results for post file inclusion +!-------------------------------------------------------------------------------------------------- +pure function homogenization_isostrain_postResults(state,i,e) + use prec, only: & + pReal,& + p_vec + use mesh, only: & + mesh_element + use material, only: & + homogenization_typeInstance, & + homogenization_Noutput implicit none type(p_vec), intent(in) :: state - integer(pInt), intent(in) :: ip,el + 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,el)))) :: homogenization_isostrain_postResults + (homogenization_typeInstance(mesh_element(3,e)))) :: homogenization_isostrain_postResults + c = 0_pInt - homID = homogenization_typeInstance(mesh_element(3,el)) + homID = homogenization_typeInstance(mesh_element(3,e)) homogenization_isostrain_postResults = 0.0_pReal - do o = 1_pInt,homogenization_Noutput(mesh_element(3,el)) + do o = 1_pInt,homogenization_Noutput(mesh_element(3,e)) select case(homogenization_isostrain_output(o,homID)) case ('ngrains') homogenization_isostrain_postResults(c+1_pInt) = real(homogenization_isostrain_Ngrains(homID),pReal)