From 65ae9799207924137fdfb605214a48b7145a1c47 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 12 Dec 2013 17:09:59 +0000 Subject: [PATCH] indroduced sourced allocation, enums where applicable (some parts still missing). fixed bug when having recursive file input --- code/DAMASK_spectral_driver.f90 | 14 +- code/IO.f90 | 4 +- code/constitutive.f90 | 38 ++-- code/constitutive_dislotwin.f90 | 10 +- code/constitutive_j2.f90 | 5 +- code/constitutive_none.f90 | 2 +- code/constitutive_nonlocal.f90 | 4 +- code/constitutive_phenopowerlaw.f90 | 25 +-- code/constitutive_titanmod.f90 | 21 +- code/crystallite.f90 | 309 +++++++++++++++++----------- code/homogenization.f90 | 24 +-- code/homogenization_RGC.f90 | 25 ++- code/homogenization_isostrain.f90 | 109 +++++----- code/material.f90 | 244 ++++++++++++---------- 14 files changed, 471 insertions(+), 363 deletions(-) diff --git a/code/DAMASK_spectral_driver.f90 b/code/DAMASK_spectral_driver.f90 index ab83105be..357ef812b 100644 --- a/code/DAMASK_spectral_driver.f90 +++ b/code/DAMASK_spectral_driver.f90 @@ -110,7 +110,7 @@ program DAMASK_spectral_Driver integer(pInt), parameter :: maxNchunks = (1_pInt + 9_pInt)*3_pInt + & ! deformation, rotation, and stress (1_pInt + 1_pInt)*5_pInt + & ! time, (log)incs, temp, restartfrequency, and outputfrequency 1_pInt, & ! dropguessing - myUnit = 234_pInt !< file unit, DAMASK IO does not support newunit feature + FILEUNIT = 234_pInt !< file unit, DAMASK IO does not support newunit feature integer(pInt), dimension(1_pInt + maxNchunks*2_pInt) :: positions ! this is longer than needed for geometry parsing integer(pInt) :: & @@ -165,10 +165,10 @@ program DAMASK_spectral_Driver !-------------------------------------------------------------------------------------------------- ! reading basic information from load case file and allocate data structure containing load cases - call IO_open_file(myUnit,trim(loadCaseFile)) - rewind(myUnit) + call IO_open_file(FILEUNIT,trim(loadCaseFile)) + rewind(FILEUNIT) do - line = IO_read(myUnit) + line = IO_read(FILEUNIT) if (trim(line) == '#EOF#') exit if (IO_isBlank(line)) cycle ! skip empty lines positions = IO_stringPos(line,maxNchunks) @@ -191,9 +191,9 @@ program DAMASK_spectral_Driver !-------------------------------------------------------------------------------------------------- ! reading the load case and assign values to the allocated data structure - rewind(myUnit) + rewind(FILEUNIT) do - line = IO_read(myUnit) + line = IO_read(FILEUNIT) if (trim(line) == '#EOF#') exit if (IO_isBlank(line)) cycle ! skip empty lines currentLoadCase = currentLoadCase + 1_pInt @@ -289,7 +289,7 @@ program DAMASK_spectral_Driver loadCases(currentLoadCase)%rotation = math_plain9to33(temp_valueVector) end select enddo; enddo - close(myUnit) + close(FILEUNIT) ! reorder phase field data to remove redundant non-active fields nActivePhaseFields = 0_pInt do i = 1, maxPhaseFields diff --git a/code/IO.f90 b/code/IO.f90 index 6156e09cf..95cb07b75 100644 --- a/code/IO.f90 +++ b/code/IO.f90 @@ -1452,6 +1452,8 @@ subroutine IO_error(error_ID,el,ip,g,ext_msg) msg = 'could not assemble input files' case (104_pInt) msg = '{input} recursion limit reached' + case (105_pInt) + msg = 'unknown output:' !-------------------------------------------------------------------------------------------------- ! material error messages and related messages in mesh @@ -1493,8 +1495,6 @@ subroutine IO_error(error_ID,el,ip,g,ext_msg) msg = 'unknown material parameter:' case (211_pInt) msg = 'material parameter out of bounds:' - case (212_pInt) - msg = 'unknown plasticity output:' case (214_pInt) msg = 'stiffness parameter close to zero:' diff --git a/code/constitutive.f90 b/code/constitutive.f90 index 6d1471ce6..ef08a8dd4 100644 --- a/code/constitutive.f90 +++ b/code/constitutive.f90 @@ -130,7 +130,7 @@ subroutine constitutive_init use constitutive_nonlocal implicit none - integer(pInt), parameter :: fileunit = 200_pInt + integer(pInt), parameter :: FILEUNIT = 200_pInt integer(pInt) :: & g, & !< grain number i, & !< integration point number @@ -152,15 +152,15 @@ subroutine constitutive_init !-------------------------------------------------------------------------------------------------- ! parse plasticities from config file - if (.not. IO_open_jobFile_stat(fileunit,material_localFileExt)) & ! no local material configuration present... - call IO_open_file(fileunit,material_configFile) ! ... open material.config file - call constitutive_none_init(fileunit) - call constitutive_j2_init(fileunit) - call constitutive_phenopowerlaw_init(fileunit) - call constitutive_titanmod_init(fileunit) - call constitutive_dislotwin_init(fileunit) - call constitutive_nonlocal_init(fileunit) - close(fileunit) + if (.not. IO_open_jobFile_stat(FILEUNIT,material_localFileExt)) & ! no local material configuration present... + call IO_open_file(FILEUNIT,material_configFile) ! ... open material.config file + call constitutive_none_init(FILEUNIT) + call constitutive_j2_init(FILEUNIT) + call constitutive_phenopowerlaw_init(FILEUNIT) + call constitutive_titanmod_init(FILEUNIT) + call constitutive_dislotwin_init(FILEUNIT) + call constitutive_nonlocal_init(FILEUNIT) + close(FILEUNIT) write(6,'(/,a)') ' <<<+- constitutive init -+>>>' write(6,'(a)') ' $Id$' @@ -169,7 +169,7 @@ subroutine constitutive_init !-------------------------------------------------------------------------------------------------- ! write description file for constitutive phase output - call IO_write_jobFile(fileunit,'outputConstitutive') + call IO_write_jobFile(FILEUNIT,'outputConstitutive') do p = 1_pInt,material_Nphase i = phase_plasticityInstance(p) ! which instance of a plasticity is present phase knownPlasticity = .true. ! assume valid @@ -201,15 +201,15 @@ subroutine constitutive_init case default knownPlasticity = .false. end select - write(fileunit,'(/,a,/)') '['//trim(phase_name(p))//']' + write(FILEUNIT,'(/,a,/)') '['//trim(phase_name(p))//']' if (knownPlasticity) then - write(fileunit,'(a)') '(plasticity)'//char(9)//trim(outputName) + write(FILEUNIT,'(a)') '(plasticity)'//char(9)//trim(outputName) do e = 1_pInt,phase_Noutput(p) - write(fileunit,'(a,i4)') trim(thisOutput(e,i))//char(9),thisSize(e,i) + write(FILEUNIT,'(a,i4)') trim(thisOutput(e,i))//char(9),thisSize(e,i) enddo endif enddo - close(fileunit) + close(FILEUNIT) !-------------------------------------------------------------------------------------------------- ! allocation of states @@ -217,7 +217,7 @@ subroutine constitutive_init iMax = mesh_maxNips eMax = mesh_NcpElems - allocate(constitutive_state0(cMax,iMax,eMax)) + allocate(constitutive_state0(cMax,iMax,eMax)) allocate(constitutive_partionedState0(cMax,iMax,eMax)) allocate(constitutive_subState0(cMax,iMax,eMax)) allocate(constitutive_state(cMax,iMax,eMax)) @@ -226,9 +226,9 @@ subroutine constitutive_init allocate(constitutive_deltaState(cMax,iMax,eMax)) allocate(constitutive_dotState_backup(cMax,iMax,eMax)) allocate(constitutive_aTolState(cMax,iMax,eMax)) - allocate(constitutive_sizeDotState(cMax,iMax,eMax)) ; constitutive_sizeDotState = 0_pInt - allocate(constitutive_sizeState(cMax,iMax,eMax)) ; constitutive_sizeState = 0_pInt - allocate(constitutive_sizePostResults(cMax,iMax,eMax)); constitutive_sizePostResults = 0_pInt + allocate(constitutive_sizeDotState(cMax,iMax,eMax), source=0_pInt) + allocate(constitutive_sizeState(cMax,iMax,eMax), source=0_pInt) + allocate(constitutive_sizePostResults(cMax,iMax,eMax), source=0_pInt) if (any(numerics_integrator == 1_pInt)) then allocate(constitutive_previousDotState(cMax,iMax,eMax)) allocate(constitutive_previousDotState2(cMax,iMax,eMax)) diff --git a/code/constitutive_dislotwin.f90 b/code/constitutive_dislotwin.f90 index bf15d0405..ab97fd6d4 100644 --- a/code/constitutive_dislotwin.f90 +++ b/code/constitutive_dislotwin.f90 @@ -318,16 +318,16 @@ subroutine constitutive_dislotwin_init(fileUnit) line = IO_read(fileUnit) enddo - do while (trim(line) /= IO_EOF) ! read thru sections of phase part + do while (trim(line) /= IO_EOF) ! read through sections of phase part line = IO_read(fileUnit) if (IO_isBlank(line)) cycle ! skip empty lines if (IO_getTag(line,'<','>') /= '') then ! stop at next part - line = IO_read(fileUnit, .true.) ! reset IO_read + line = IO_read(fileUnit, .true.) ! reset IO_read exit endif if (IO_getTag(line,'[',']') /= '') then ! next section section = section + 1_pInt ! advance section counter - cycle + cycle ! skip to next line endif if (section > 0_pInt ) then ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran if (phase_plasticity(section) == PLASTICITY_DISLOTWIN_ID) then ! one of my sections @@ -378,6 +378,8 @@ subroutine constitutive_dislotwin_init(fileUnit) constitutive_dislotwin_outputID(constitutive_dislotwin_Noutput(i),i) = sb_eigenvalues_ID case ('sb_eigenvectors') constitutive_dislotwin_outputID(constitutive_dislotwin_Noutput(i),i) = sb_eigenvectors_ID + case default + call IO_error(105_pInt,ext_msg=IO_stringValue(line,positions,2_pInt)//' ('//PLASTICITY_DISLOTWIN_label//')') end select case ('lattice_structure') structure = IO_lc(IO_stringValue(line,positions,2_pInt)) @@ -671,8 +673,6 @@ subroutine constitutive_dislotwin_init(fileUnit) mySize = 3_pInt case(sb_eigenvectors_ID) mySize = 9_pInt - case default - call IO_error(212_pInt,ext_msg=constitutive_dislotwin_output(o,i)//' ('//PLASTICITY_DISLOTWIN_label//')') end select if (mySize > 0_pInt) then ! any meaningful output found diff --git a/code/constitutive_j2.f90 b/code/constitutive_j2.f90 index 36317e1bd..2e57867ad 100644 --- a/code/constitutive_j2.f90 +++ b/code/constitutive_j2.f90 @@ -131,7 +131,7 @@ subroutine constitutive_j2_init(fileUnit) structure = '' character(len=65536) :: & tag = '', & - line = '' ! to start initialized + line = '' write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_J2_label//' init -+>>>' write(6,'(a)') ' $Id$' @@ -201,6 +201,8 @@ subroutine constitutive_j2_init(fileUnit) constitutive_j2_outputID(constitutive_j2_Noutput(i),i) = flowstress_ID case ('strainrate') constitutive_j2_outputID(constitutive_j2_Noutput(i),i) = strainrate_ID + case default + call IO_error(105_pInt,ext_msg=IO_stringValue(line,positions,2_pInt)//' ('//PLASTICITY_J2_label//')') end select case ('lattice_structure') structure = IO_lc(IO_stringValue(line,positions,2_pInt)) @@ -290,7 +292,6 @@ subroutine constitutive_j2_init(fileUnit) case(flowstress_ID,strainrate_ID) mySize = 1_pInt case default - call IO_error(212_pInt,ext_msg=constitutive_j2_output(o,i)//' ('//PLASTICITY_J2_label//')') end select if (mySize > 0_pInt) then ! any meaningful output found diff --git a/code/constitutive_none.f90 b/code/constitutive_none.f90 index 0f3b87f30..5e39881dd 100644 --- a/code/constitutive_none.f90 +++ b/code/constitutive_none.f90 @@ -124,7 +124,7 @@ subroutine constitutive_none_init(fileUnit) endif if (IO_getTag(line,'[',']') /= '') then ! next section section = section + 1_pInt ! advance section counter - cycle + cycle ! skip to next line endif if (section > 0_pInt ) then ! do not short-circuit here (.and. with next if-statement). It's not safe in Fortran if (phase_plasticity(section) == PLASTICITY_NONE_ID) then ! one of my sections diff --git a/code/constitutive_nonlocal.f90 b/code/constitutive_nonlocal.f90 index 157e6a937..9350e02b9 100644 --- a/code/constitutive_nonlocal.f90 +++ b/code/constitutive_nonlocal.f90 @@ -630,6 +630,8 @@ do while (trim(line) /= IO_EOF) constitutive_nonlocal_outputID(Noutput(i),i) = accumulatedshear_ID case('dislocationstress') constitutive_nonlocal_outputID(Noutput(i),i) = dislocationstress_ID + case default + call IO_error(105_pInt,ext_msg=IO_stringValue(line,positions,2_pInt)//' ('//PLASTICITY_NONLOCAL_label//')') end select case ('lattice_structure') structure = IO_lc(IO_stringValue(line,positions,2_pInt)) @@ -1117,8 +1119,6 @@ instancesLoop: do i = 1,maxNmatIDs case(dislocationstress_ID) mySize = 6_pInt case default - call IO_error(212_pInt,ext_msg=constitutive_nonlocal_output(o,i)//& - '('//PLASTICITY_NONLOCAL_label//')') end select if (mySize > 0_pInt) then ! any meaningful output found diff --git a/code/constitutive_phenopowerlaw.f90 b/code/constitutive_phenopowerlaw.f90 index ddc7f0cbb..3c5c7c167 100644 --- a/code/constitutive_phenopowerlaw.f90 +++ b/code/constitutive_phenopowerlaw.f90 @@ -232,13 +232,17 @@ subroutine constitutive_phenopowerlaw_init(fileUnit) source=0.0_pReal) rewind(fileUnit) - do while (trim(line) /= '#EOF#' .and. IO_lc(IO_getTag(line,'<','>')) /= 'phase') ! wind forward to + do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= 'phase') ! wind forward to line = IO_read(fileUnit) enddo - do while (trim(line) /= '#EOF#') ! read through sections of phase part + + do while (trim(line) /= IO_EOF) ! read through sections of phase part line = IO_read(fileUnit) if (IO_isBlank(line)) cycle ! skip empty lines - if (IO_getTag(line,'<','>') /= '') exit ! stop at next part + 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 ! advance section counter cycle ! skip to next line @@ -276,6 +280,8 @@ subroutine constitutive_phenopowerlaw_init(fileUnit) constitutive_phenopowerlaw_outputID(constitutive_phenopowerlaw_Noutput(i),i) = resolvedstress_twin_ID case ('totalvolfrac') constitutive_phenopowerlaw_outputID(constitutive_phenopowerlaw_Noutput(i),i) = totalvolfrac_ID + case default + call IO_error(105_pInt,ext_msg=IO_stringValue(line,positions,2_pInt)//' ('//PLASTICITY_PHENOPOWERLAW_label//')') end select case ('lattice_structure') structure = IO_lc(IO_stringValue(line,positions,2_pInt)) @@ -484,20 +490,16 @@ subroutine constitutive_phenopowerlaw_init(fileUnit) ! allocation of variables whose size depends on the total number of active slip systems allocate(constitutive_phenopowerlaw_hardeningMatrix_SlipSlip(maxval(constitutive_phenopowerlaw_totalNslip),& ! slip resistance from slip activity maxval(constitutive_phenopowerlaw_totalNslip),& - maxNinstance)) + maxNinstance), source=0.0_pReal) allocate(constitutive_phenopowerlaw_hardeningMatrix_SlipTwin(maxval(constitutive_phenopowerlaw_totalNslip),& ! slip resistance from twin activity maxval(constitutive_phenopowerlaw_totalNtwin),& - maxNinstance)) + maxNinstance), source=0.0_pReal) allocate(constitutive_phenopowerlaw_hardeningMatrix_TwinSlip(maxval(constitutive_phenopowerlaw_totalNtwin),& ! twin resistance from slip activity maxval(constitutive_phenopowerlaw_totalNslip),& - maxNinstance)) + maxNinstance), source=0.0_pReal) allocate(constitutive_phenopowerlaw_hardeningMatrix_TwinTwin(maxval(constitutive_phenopowerlaw_totalNtwin),& ! twin resistance from twin activity maxval(constitutive_phenopowerlaw_totalNtwin),& - maxNinstance)) - constitutive_phenopowerlaw_hardeningMatrix_SlipSlip = 0.0_pReal - constitutive_phenopowerlaw_hardeningMatrix_SlipTwin = 0.0_pReal - constitutive_phenopowerlaw_hardeningMatrix_TwinSlip = 0.0_pReal - constitutive_phenopowerlaw_hardeningMatrix_TwinTwin = 0.0_pReal + maxNinstance), source=0.0_pReal) instancesLoop: do i = 1_pInt,maxNinstance outputsLoop: do o = 1_pInt,constitutive_phenopowerlaw_Noutput(i) @@ -519,7 +521,6 @@ subroutine constitutive_phenopowerlaw_init(fileUnit) ) mySize = 1_pInt case default - call IO_error(212_pInt,ext_msg=constitutive_phenopowerlaw_output(o,i)//' ('//PLASTICITY_PHENOPOWERLAW_label//')') end select outputFound: if (mySize > 0_pInt) then diff --git a/code/constitutive_titanmod.f90 b/code/constitutive_titanmod.f90 index e10d388e6..5effdd7a1 100644 --- a/code/constitutive_titanmod.f90 +++ b/code/constitutive_titanmod.f90 @@ -200,7 +200,7 @@ module constitutive_titanmod !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -subroutine constitutive_titanmod_init(myFile) +subroutine constitutive_titanmod_init(fileUnit) 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,& @@ -215,7 +215,7 @@ subroutine constitutive_titanmod_init(myFile) use lattice implicit none - integer(pInt), intent(in) :: myFile + integer(pInt), intent(in) :: fileUnit integer(pInt), parameter :: MAXNCHUNKS = LATTICE_maxNinteraction + 1_pInt integer(pInt), dimension(1_pInt+2_pInt*MAXNCHUNKS) :: positions @@ -392,15 +392,18 @@ subroutine constitutive_titanmod_init(myFile) allocate(constitutive_titanmod_interactionTwinTwin(lattice_maxNinteraction,maxNinstance)) constitutive_titanmod_interactionTwinTwin = 0.0_pReal - rewind(myFile) - do while (trim(line) /= '#EOF#' .and. IO_lc(IO_getTag(line,'<','>')) /= 'phase') ! wind forward to - line = IO_read(myFile) + rewind(fileUnit) + do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= 'phase') ! wind forward to + line = IO_read(fileUnit) enddo - do while (trim(line) /= '#EOF#') ! read through sections of phase part - line = IO_read(myFile) + do while (trim(line) /= IO_EOF) ! read through sections of phase part + line = IO_read(fileUnit) if (IO_isBlank(line)) cycle ! skip empty lines - if (IO_getTag(line,'<','>') /= '') exit ! stop at next part + 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 ! advance section counter cycle ! skip to next line @@ -843,7 +846,7 @@ subroutine constitutive_titanmod_init(myFile) 'shear_total') mySize = 1_pInt case default - call IO_error(212_pInt,ext_msg=constitutive_titanmod_output(o,i)// & + call IO_error(105_pInt,ext_msg=constitutive_titanmod_output(o,i)// & ' ('//PLASTICITY_TITANMOD_label//')') end select diff --git a/code/crystallite.f90 b/code/crystallite.f90 index bcf2389bf..49fb51015 100644 --- a/code/crystallite.f90 +++ b/code/crystallite.f90 @@ -32,37 +32,37 @@ module crystallite implicit none private - character(len=64), dimension(:,:), allocatable, private :: & + character(len=64), dimension(:,:), allocatable, private :: & crystallite_output !< name of each post result output - integer(pInt), public, protected :: & + integer(pInt), public, protected :: & crystallite_maxSizePostResults !< description not available - integer(pInt), dimension(:), allocatable, public, protected :: & + integer(pInt), dimension(:), allocatable, public, protected :: & crystallite_sizePostResults !< description not available - integer(pInt), dimension(:,:), allocatable, private :: & + integer(pInt), dimension(:,:), allocatable, private :: & crystallite_sizePostResult !< description not available - integer(pInt), dimension(:,:,:), allocatable, private :: & + integer(pInt), dimension(:,:,:), allocatable, private :: & crystallite_symmetryID !< crystallographic symmetry 1=cubic 2=hexagonal, needed in all orientation calcs - real(pReal), dimension(:,:), allocatable, public :: & + real(pReal), dimension(:,:), allocatable, public :: & crystallite_temperature !< temperature (same on all components on one IP) - real(pReal), dimension(:,:,:), allocatable, public, protected :: & + real(pReal), dimension(:,:,:), allocatable, public, protected :: & crystallite_heat !< heat source - real(pReal), dimension(:,:,:), allocatable, public :: & + real(pReal), dimension(:,:,:), allocatable, public :: & crystallite_dt !< requested time increment of each grain - real(pReal), dimension(:,:,:), allocatable, private :: & + real(pReal), dimension(:,:,:), allocatable, private :: & crystallite_subdt, & !< substepped time increment of each grain crystallite_subFrac, & !< already calculated fraction of increment crystallite_subStep !< size of next integration step - real(pReal), dimension(:,:,:,:), allocatable, public :: & + real(pReal), dimension(:,:,:,:), allocatable, public :: & crystallite_Tstar_v, & !< current 2nd Piola-Kirchhoff stress vector (end of converged time step) crystallite_Tstar0_v, & !< 2nd Piola-Kirchhoff stress vector at start of FE inc crystallite_partionedTstar0_v !< 2nd Piola-Kirchhoff stress vector at start of homog inc - real(pReal), dimension(:,:,:,:), allocatable, private :: & + real(pReal), dimension(:,:,:,:), allocatable, private :: & crystallite_subTstar0_v, & !< 2nd Piola-Kirchhoff stress vector at start of crystallite inc crystallite_orientation, & !< orientation as quaternion crystallite_orientation0, & !< initial orientation as quaternion crystallite_rotation !< grain rotation away from initial orientation as axis-angle (in degrees) in crystal reference frame - real(pReal), dimension(:,:,:,:,:), allocatable, public :: & + real(pReal), dimension(:,:,:,:,:), allocatable, public :: & crystallite_Fp, & !< current plastic def grad (end of converged time step) crystallite_Fp0, & !< plastic def grad at start of FE inc crystallite_partionedFp0,& !< plastic def grad at start of homog inc @@ -73,7 +73,7 @@ module crystallite crystallite_Lp0, & !< plastic velocitiy grad at start of FE inc crystallite_partionedLp0,& !< plastic velocity grad at start of homog inc crystallite_P !< 1st Piola-Kirchhoff stress per grain - real(pReal), dimension(:,:,:,:,:), allocatable, private :: & + real(pReal), dimension(:,:,:,:,:), allocatable, private :: & crystallite_Fe, & !< current "elastic" def grad (end of converged time step) crystallite_subFe0,& !< "elastic" def grad at start of crystallite inc crystallite_invFp, & !< inverse of current plastic def grad (end of converged time step) @@ -82,26 +82,55 @@ module crystallite crystallite_subF0, & !< def grad at start of crystallite inc crystallite_subLp0,& !< plastic velocity grad at start of crystallite inc crystallite_disorientation !< disorientation between two neighboring ips (only calculated for single grain IPs) - real(pReal), dimension(:,:,:,:,:,:,:), allocatable, public :: & + real(pReal), dimension(:,:,:,:,:,:,:), allocatable, public :: & crystallite_dPdF, & !< current individual dPdF per grain (end of converged time step) crystallite_dPdF0, & !< individual dPdF per grain at start of FE inc crystallite_partioneddPdF0 !< individual dPdF per grain at start of homog inc - real(pReal), dimension(:,:,:,:,:,:,:), allocatable, private :: & + real(pReal), dimension(:,:,:,:,:,:,:), allocatable, private :: & crystallite_fallbackdPdF !< dPdF fallback for non-converged grains (elastic prediction) - logical, dimension(:,:,:), allocatable, public :: & + logical, dimension(:,:,:), allocatable, public :: & crystallite_requested !< flag to request crystallite calculation - logical, dimension(:,:,:), allocatable, public, protected :: & + logical, dimension(:,:,:), allocatable, public, protected :: & crystallite_converged, & !< convergence flag crystallite_localPlasticity !< indicates this grain to have purely local constitutive law - logical, dimension(:,:,:), allocatable, private :: & + logical, dimension(:,:,:), allocatable, private :: & crystallite_todo !< flag to indicate need for further computation - logical, dimension(:,:), allocatable, private :: & + logical, dimension(:,:), allocatable, private :: & crystallite_clearToWindForward, & !< description not available crystallite_clearToCutback, & !< description not available crystallite_syncSubFrac, & !< description not available crystallite_syncSubFracCompleted, & !< description not available crystallite_neighborEnforcedCutback !< description not available + real(pReal), dimension(:,:,:), allocatable, private :: & + constitutive_j2_Cslip_66 + enum, bind(c) + enumerator :: undefined_ID, & + phase_ID, & + texture_ID, & + volume_ID, & + grainrotationx_ID, & + grainrotationy_ID, & + grainrotationz_ID, & + heat_ID, & + orientation_ID, & + grainrotation_ID, & + eulerangles_ID, & + defgrad_ID, & + fe_ID, & + fp_ID, & + lp_ID, & + e_ID, & + ee_ID, & + p_ID, & + s_ID, & + elasmatrix_ID, & + neighboringip_ID, & + neighboringelement_ID + end enum + integer(kind(undefined_ID)),dimension(:,:), allocatable, private :: & + crystallite_outputID !< ID of each post result output + public :: & crystallite_init, & crystallite_stressAndItsTangent, & @@ -158,7 +187,8 @@ subroutine crystallite_init(temperature) IO_stringPos, & IO_stringValue, & IO_write_jobFile, & - IO_error + IO_error, & + IO_EOF use material use lattice, only: & lattice_symmetryType @@ -176,7 +206,7 @@ subroutine crystallite_init(temperature) implicit none real(pReal), intent(in) :: temperature integer(pInt), parameter :: & - MYUNIT = 200_pInt, & + FILEUNIT = 200_pInt, & MAXNCHUNKS = 2_pInt integer(pInt), dimension(1+2*MAXNCHUNKS) :: positions @@ -189,7 +219,7 @@ subroutine crystallite_init(temperature) eMax, & !< maximum number of elements nMax, & !< maximum number of ip neighbors myNgrains, & !< number of grains in current IP - section, & + section = 0_pInt, & j, & p, & output, & @@ -197,83 +227,86 @@ subroutine crystallite_init(temperature) myPhase, & myMat character(len=65536) :: & - tag, & - line + tag = '', & + line= '' write(6,'(/,a)') ' <<<+- crystallite init -+>>>' write(6,'(a)') ' $Id$' write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" - - line = '' - section = 0_pInt gMax = homogenization_maxNgrains iMax = mesh_maxNips eMax = mesh_NcpElems nMax = mesh_maxNipNeighbors - allocate(crystallite_temperature(iMax,eMax)); crystallite_temperature = temperature - allocate(crystallite_heat(gMax,iMax,eMax)); crystallite_heat = 0.0_pReal - allocate(crystallite_Tstar0_v(6,gMax,iMax,eMax)); crystallite_Tstar0_v = 0.0_pReal - allocate(crystallite_partionedTstar0_v(6,gMax,iMax,eMax)); crystallite_partionedTstar0_v = 0.0_pReal - allocate(crystallite_subTstar0_v(6,gMax,iMax,eMax)); crystallite_subTstar0_v = 0.0_pReal - allocate(crystallite_Tstar_v(6,gMax,iMax,eMax)); crystallite_Tstar_v = 0.0_pReal - allocate(crystallite_P(3,3,gMax,iMax,eMax)); crystallite_P = 0.0_pReal - allocate(crystallite_F0(3,3,gMax,iMax,eMax)); crystallite_F0 = 0.0_pReal - allocate(crystallite_partionedF0(3,3,gMax,iMax,eMax)); crystallite_partionedF0 = 0.0_pReal - allocate(crystallite_partionedF(3,3,gMax,iMax,eMax)); crystallite_partionedF = 0.0_pReal - allocate(crystallite_subF0(3,3,gMax,iMax,eMax)); crystallite_subF0 = 0.0_pReal - allocate(crystallite_subF(3,3,gMax,iMax,eMax)); crystallite_subF = 0.0_pReal - allocate(crystallite_Fp0(3,3,gMax,iMax,eMax)); crystallite_Fp0 = 0.0_pReal - allocate(crystallite_partionedFp0(3,3,gMax,iMax,eMax)); crystallite_partionedFp0 = 0.0_pReal - allocate(crystallite_subFp0(3,3,gMax,iMax,eMax)); crystallite_subFp0 = 0.0_pReal - allocate(crystallite_Fp(3,3,gMax,iMax,eMax)); crystallite_Fp = 0.0_pReal - allocate(crystallite_invFp(3,3,gMax,iMax,eMax)); crystallite_invFp = 0.0_pReal - allocate(crystallite_Fe(3,3,gMax,iMax,eMax)); crystallite_Fe = 0.0_pReal - allocate(crystallite_subFe0(3,3,gMax,iMax,eMax)); crystallite_subFe0 = 0.0_pReal - allocate(crystallite_Lp0(3,3,gMax,iMax,eMax)); crystallite_Lp0 = 0.0_pReal - allocate(crystallite_partionedLp0(3,3,gMax,iMax,eMax)); crystallite_partionedLp0 = 0.0_pReal - allocate(crystallite_subLp0(3,3,gMax,iMax,eMax)); crystallite_subLp0 = 0.0_pReal - allocate(crystallite_Lp(3,3,gMax,iMax,eMax)); crystallite_Lp = 0.0_pReal - allocate(crystallite_dPdF(3,3,3,3,gMax,iMax,eMax)); crystallite_dPdF = 0.0_pReal - allocate(crystallite_dPdF0(3,3,3,3,gMax,iMax,eMax)); crystallite_dPdF0 = 0.0_pReal - allocate(crystallite_partioneddPdF0(3,3,3,3,gMax,iMax,eMax)); crystallite_partioneddPdF0 = 0.0_pReal - allocate(crystallite_fallbackdPdF(3,3,3,3,gMax,iMax,eMax)); crystallite_fallbackdPdF = 0.0_pReal - allocate(crystallite_dt(gMax,iMax,eMax)); crystallite_dt = 0.0_pReal - allocate(crystallite_subdt(gMax,iMax,eMax)); crystallite_subdt = 0.0_pReal - allocate(crystallite_subFrac(gMax,iMax,eMax)); crystallite_subFrac = 0.0_pReal - allocate(crystallite_subStep(gMax,iMax,eMax)); crystallite_subStep = 0.0_pReal - allocate(crystallite_orientation(4,gMax,iMax,eMax)); crystallite_orientation = 0.0_pReal - allocate(crystallite_orientation0(4,gMax,iMax,eMax)); crystallite_orientation0 = 0.0_pReal - allocate(crystallite_rotation(4,gMax,iMax,eMax)); crystallite_rotation = 0.0_pReal - allocate(crystallite_disorientation(4,nMax,gMax,iMax,eMax)); crystallite_disorientation = 0.0_pReal - allocate(crystallite_symmetryID(gMax,iMax,eMax)); crystallite_symmetryID = 0_pInt - allocate(crystallite_localPlasticity(gMax,iMax,eMax)); crystallite_localPlasticity = .true. - allocate(crystallite_requested(gMax,iMax,eMax)); crystallite_requested = .false. - allocate(crystallite_todo(gMax,iMax,eMax)); crystallite_todo = .false. - allocate(crystallite_converged(gMax,iMax,eMax)); crystallite_converged = .true. - allocate(crystallite_clearToWindForward(iMax,eMax)); crystallite_clearToWindForward = .true. - allocate(crystallite_syncSubFrac(iMax,eMax)); crystallite_syncSubFrac = .false. - allocate(crystallite_syncSubFracCompleted(iMax,eMax)); crystallite_syncSubFracCompleted = .false. - allocate(crystallite_clearToCutback(iMax,eMax)); crystallite_clearToCutback = .true. - allocate(crystallite_neighborEnforcedCutback(iMax,eMax)); crystallite_neighborEnforcedCutback = .false. + allocate(crystallite_temperature(iMax,eMax), source=temperature) + allocate(crystallite_heat(gMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_Tstar0_v(6,gMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_partionedTstar0_v(6,gMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_subTstar0_v(6,gMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_Tstar_v(6,gMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_P(3,3,gMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_F0(3,3,gMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_partionedF0(3,3,gMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_partionedF(3,3,gMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_subF0(3,3,gMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_subF(3,3,gMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_Fp0(3,3,gMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_partionedFp0(3,3,gMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_subFp0(3,3,gMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_Fp(3,3,gMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_invFp(3,3,gMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_Fe(3,3,gMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_subFe0(3,3,gMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_Lp0(3,3,gMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_partionedLp0(3,3,gMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_subLp0(3,3,gMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_Lp(3,3,gMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_dPdF(3,3,3,3,gMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_dPdF0(3,3,3,3,gMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_partioneddPdF0(3,3,3,3,gMax,iMax,eMax),source=0.0_pReal) + allocate(crystallite_fallbackdPdF(3,3,3,3,gMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_dt(gMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_subdt(gMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_subFrac(gMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_subStep(gMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_orientation(4,gMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_orientation0(4,gMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_rotation(4,gMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_disorientation(4,nMax,gMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_symmetryID(gMax,iMax,eMax), source=0_pInt) + allocate(crystallite_localPlasticity(gMax,iMax,eMax), source=.true.) + allocate(crystallite_requested(gMax,iMax,eMax), source=.false.) + allocate(crystallite_todo(gMax,iMax,eMax), source=.false.) + allocate(crystallite_converged(gMax,iMax,eMax), source=.true.) + allocate(crystallite_clearToWindForward(iMax,eMax), source=.true.) + allocate(crystallite_syncSubFrac(iMax,eMax), source=.false.) + allocate(crystallite_syncSubFracCompleted(iMax,eMax), source=.false.) + allocate(crystallite_clearToCutback(iMax,eMax), source=.true.) + allocate(crystallite_neighborEnforcedCutback(iMax,eMax), source=.false.) allocate(crystallite_output(maxval(crystallite_Noutput), & - material_Ncrystallite)) ; crystallite_output = '' - allocate(crystallite_sizePostResults(material_Ncrystallite)) ; crystallite_sizePostResults = 0_pInt + material_Ncrystallite)) ; crystallite_output = '' + allocate(crystallite_outputID(maxval(crystallite_Noutput), & + material_Ncrystallite), source=undefined_ID) + allocate(crystallite_sizePostResults(material_Ncrystallite),source=0_pInt) allocate(crystallite_sizePostResult(maxval(crystallite_Noutput), & - material_Ncrystallite)) ; crystallite_sizePostResult = 0_pInt + material_Ncrystallite), source=0_pInt) - if (.not. IO_open_jobFile_stat(myUnit,material_localFileExt)) & ! no local material configuration present... - call IO_open_file(myUnit,material_configFile) ! ...open material.config file - do while (trim(line) /= '#EOF#' .and. IO_lc(IO_getTag(line,'<','>')) /= material_partCrystallite) ! wind forward to - line = IO_read(myUnit) + if (.not. IO_open_jobFile_stat(FILEUNIT,material_localFileExt)) & ! no local material configuration present... + call IO_open_file(FILEUNIT,material_configFile) ! ...open material.config file + rewind(FILEUNIT) + do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= material_partCrystallite) ! wind forward to + line = IO_read(FILEUNIT) enddo - do while (trim(line) /= '#EOF#') ! read thru sections of phase part - line = IO_read(myUnit) + do while (trim(line) /= IO_EOF) ! read through sections of crystallite part + line = IO_read(FILEUNIT) if (IO_isBlank(line)) cycle ! skip empty lines - if (IO_getTag(line,'<','>') /= '') exit ! stop at next part + 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 output = 0_pInt ! reset output counter @@ -285,26 +318,70 @@ subroutine crystallite_init(temperature) case ('(output)') output = output + 1_pInt crystallite_output(output,section) = IO_lc(IO_stringValue(line,positions,2_pInt)) + select case(IO_lc(IO_stringValue(line,positions,2_pInt))) + case ('phase') + crystallite_outputID = phase_ID + case ('texture') + crystallite_outputID = texture_ID + case ('volume') + crystallite_outputID = volume_ID + case ('grainrotationx') + crystallite_outputID = grainrotationx_ID + case ('grainrotationy') + crystallite_outputID = grainrotationy_ID + case ('grainrotationz') + crystallite_outputID = grainrotationx_ID + case ('heat') + crystallite_outputID = heat_ID + case ('orientation') + crystallite_outputID = orientation_ID + case ('grainrotation') + crystallite_outputID = grainrotation_ID + case ('eulerangles') + crystallite_outputID = eulerangles_ID + case ('defgrad','f') + crystallite_outputID = defgrad_ID + case ('fe') + crystallite_outputID = fe_ID + case ('fp') + crystallite_outputID = fp_ID + case ('e') + crystallite_outputID = e_ID + case ('ee') + crystallite_outputID = ee_ID + case ('p','firstpiola','1piola') + crystallite_outputID = p_ID + case ('s','tstar','secondpiola','2ndpiola') + crystallite_outputID = s_ID + case ('elasmatrix') + crystallite_outputID = elasmatrix_ID + case ('neighboringip') + crystallite_outputID = neighboringip_ID + case ('neighboringelement') + crystallite_outputID = neighboringelement_ID + case default + call IO_error(105_pInt,ext_msg=IO_stringValue(line,positions,2_pInt)//' (Crystallite)') + end select end select endif enddo - close(myUnit) + close(FILEUNIT) do i = 1_pInt,material_Ncrystallite do j = 1_pInt,crystallite_Noutput(i) - select case(crystallite_output(j,i)) - case('phase','texture','volume','grainrotationx','grainrotationy','grainrotationz','heat') + select case(crystallite_outputID(j,i)) + case(phase_ID,texture_ID,volume_ID,grainrotationx_ID,grainrotationy_ID,grainrotationz_ID,heat_ID) mySize = 1_pInt - case('orientation','grainrotation') ! orientation as quaternion, or deviation from initial grain orientation in axis-angle form (angle in degrees) + case(orientation_ID,grainrotation_ID) ! orientation as quaternion, or deviation from initial grain orientation in axis-angle form (angle in degrees) mySize = 4_pInt - case('eulerangles') + case(eulerangles_ID) mySize = 3_pInt - case('defgrad','f','fe','fp','lp','e','ee','p','firstpiola','1stpiola','s','tstar','secondpiola','2ndpiola') + case(defgrad_ID,fe_ID,fp_ID,lp_ID,e_ID,ee_ID,p_ID,s_ID) mySize = 9_pInt - case('elasmatrix') + case(elasmatrix_ID) mySize = 36_pInt - case('neighboringip','neighboringelement') + case(neighboringip_ID,neighboringelement_ID) mySize = mesh_maxNipNeighbors case default mySize = 0_pInt @@ -326,16 +403,16 @@ subroutine crystallite_init(temperature) !-------------------------------------------------------------------------------------------------- ! write description file for crystallite output - call IO_write_jobFile(myUnit,'outputCrystallite') + call IO_write_jobFile(FILEUNIT,'outputCrystallite') do p = 1_pInt,material_Ncrystallite - write(myUnit,'(/,a,/)') '['//trim(crystallite_name(p))//']' + write(FILEUNIT,'(/,a,/)') '['//trim(crystallite_name(p))//']' do e = 1_pInt,crystallite_Noutput(p) - write(myUnit,'(a,i4)') trim(crystallite_output(e,p))//char(9),crystallite_sizePostResult(e,p) + write(FILEUNIT,'(a,i4)') trim(crystallite_output(e,p))//char(9),crystallite_sizePostResult(e,p) enddo enddo - close(myUnit) + close(FILEUNIT) !-------------------------------------------------------------------------------------------------- ! initialize @@ -3422,42 +3499,42 @@ function crystallite_postResults(ipc, ip, el) do o = 1_pInt,crystallite_Noutput(crystID) mySize = 0_pInt - select case(crystallite_output(o,crystID)) - case ('phase') + select case(crystallite_outputID(o,crystID)) + case (phase_ID) mySize = 1_pInt crystallite_postResults(c+1) = real(material_phase(ipc,ip,el),pReal) ! phaseID of grain - case ('texture') + case (texture_ID) mySize = 1_pInt crystallite_postResults(c+1) = real(material_texture(ipc,ip,el),pReal) ! textureID of grain - case ('volume') + case (volume_ID) mySize = 1_pInt detF = math_det33(crystallite_partionedF(1:3,1:3,ipc,ip,el)) ! V_current = det(F) * V_reference crystallite_postResults(c+1) = detF * mesh_ipVolume(ip,el) & / homogenization_Ngrains(mesh_element(3,el)) ! grain volume (not fraction but absolute) - case ('heat') + case (heat_ID) mySize = 1_pInt crystallite_postResults(c+1) = crystallite_heat(ipc,ip,el) ! heat production - case ('orientation') + case (orientation_ID) mySize = 4_pInt crystallite_postResults(c+1:c+mySize) = crystallite_orientation(1:4,ipc,ip,el) ! grain orientation as quaternion - case ('eulerangles') + case (eulerangles_ID) mySize = 3_pInt crystallite_postResults(c+1:c+mySize) = inDeg & * math_qToEuler(crystallite_orientation(1:4,ipc,ip,el)) ! grain orientation as Euler angles in degree - case ('grainrotation') + case (grainrotation_ID) mySize = 4_pInt crystallite_postResults(c+1:c+mySize) = & math_qToEulerAxisAngle(crystallite_rotation(1:4,ipc,ip,el)) ! grain rotation away from initial orientation as axis-angle in sample reference coordinates crystallite_postResults(c+4) = inDeg * crystallite_postResults(c+4) ! angle in degree - case ('grainrotationx') + case (grainrotationx_ID) mySize = 1_pInt rotation = math_qToEulerAxisAngle(crystallite_rotation(1:4,ipc,ip,el)) ! grain rotation away from initial orientation as axis-angle in sample reference coordinates crystallite_postResults(c+1) = inDeg * rotation(1) * rotation(4) ! angle in degree - case ('grainrotationy') + case (grainrotationy_ID) mySize = 1_pInt rotation = math_qToEulerAxisAngle(crystallite_rotation(1:4,ipc,ip,el)) ! grain rotation away from initial orientation as axis-angle in sample reference coordinates crystallite_postResults(c+1) = inDeg * rotation(2) * rotation(4) ! angle in degree - case ('grainrotationz') + case (grainrotationz_ID) mySize = 1_pInt rotation = math_qToEulerAxisAngle(crystallite_rotation(1:4,ipc,ip,el)) ! grain rotation away from initial orientation as axis-angle in sample reference coordinates crystallite_postResults(c+1) = inDeg * rotation(3) * rotation(4) ! angle in degree @@ -3465,49 +3542,49 @@ function crystallite_postResults(ipc, ip, el) ! remark: tensor output is of the form 11,12,13, 21,22,23, 31,32,33 ! thus row index i is slow, while column index j is fast. reminder: "row is slow" - case ('defgrad','f') + case (defgrad_ID) mySize = 9_pInt crystallite_postResults(c+1:c+mySize) = & reshape(math_transpose33(crystallite_partionedF(1:3,1:3,ipc,ip,el)),[mySize]) - case ('e') + case (e_ID) mySize = 9_pInt crystallite_postResults(c+1:c+mySize) = 0.5_pReal * reshape((math_mul33x33( & math_transpose33(crystallite_partionedF(1:3,1:3,ipc,ip,el)), & crystallite_partionedF(1:3,1:3,ipc,ip,el)) - math_I3),[mySize]) - case ('fe') + case (fe_ID) mySize = 9_pInt crystallite_postResults(c+1:c+mySize) = & reshape(math_transpose33(crystallite_Fe(1:3,1:3,ipc,ip,el)),[mySize]) - case ('ee') + case (ee_ID) Ee = 0.5_pReal *(math_mul33x33(math_transpose33(crystallite_Fe(1:3,1:3,ipc,ip,el)), & crystallite_Fe(1:3,1:3,ipc,ip,el)) - math_I3) mySize = 9_pInt crystallite_postResults(c+1:c+mySize) = reshape(Ee,[mySize]) - case ('fp') + case (fp_ID) mySize = 9_pInt crystallite_postResults(c+1:c+mySize) = & reshape(math_transpose33(crystallite_Fp(1:3,1:3,ipc,ip,el)),[mySize]) - case ('lp') + case (lp_ID) mySize = 9_pInt crystallite_postResults(c+1:c+mySize) = & reshape(math_transpose33(crystallite_Lp(1:3,1:3,ipc,ip,el)),[mySize]) - case ('p','firstpiola','1stpiola') + case (p_ID) mySize = 9_pInt crystallite_postResults(c+1:c+mySize) = & reshape(math_transpose33(crystallite_P(1:3,1:3,ipc,ip,el)),[mySize]) - case ('s','tstar','secondpiola','2ndpiola') + case (s_ID) mySize = 9_pInt crystallite_postResults(c+1:c+mySize) = & reshape(math_Mandel6to33(crystallite_Tstar_v(1:6,ipc,ip,el)),[mySize]) - case ('elasmatrix') + case (elasmatrix_ID) mySize = 36_pInt crystallite_postResults(c+1:c+mySize) = reshape(constitutive_homogenizedC(ipc,ip,el),[mySize]) - case('neighboringelement') + case(neighboringelement_ID) mySize = mesh_maxNipNeighbors crystallite_postResults(c+1:c+mySize) = 0.0_pReal forall (n = 1_pInt:FE_NipNeighbors(FE_celltype(FE_geomtype(mesh_element(2,el))))) & crystallite_postResults(c+n) = real(mesh_ipNeighborhood(1,n,ip,el),pReal) - case('neighboringip') + case(neighboringip_ID) mySize = mesh_maxNipNeighbors crystallite_postResults(c+1:c+mySize) = 0.0_pReal forall (n = 1_pInt:FE_NipNeighbors(FE_celltype(FE_geomtype(mesh_element(2,el))))) & diff --git a/code/homogenization.f90 b/code/homogenization.f90 index de0075a28..e86cf27bb 100644 --- a/code/homogenization.f90 +++ b/code/homogenization.f90 @@ -123,7 +123,7 @@ subroutine homogenization_init() use homogenization_RGC implicit none - integer(pInt), parameter :: fileunit = 200_pInt + integer(pInt), parameter :: FILEUNIT = 200_pInt integer(pInt) :: e,i,p,myInstance integer(pInt), dimension(:,:), pointer :: thisSize character(len=64), dimension(:,:), pointer :: thisOutput @@ -132,15 +132,15 @@ subroutine homogenization_init() !-------------------------------------------------------------------------------------------------- ! parse homogenization from config file - if (.not. IO_open_jobFile_stat(fileunit,material_localFileExt)) & ! no local material configuration present... - call IO_open_file(fileunit,material_configFile) ! ... open material.config file - call homogenization_isostrain_init(fileunit) - call homogenization_RGC_init(fileunit) - close(fileunit) + if (.not. IO_open_jobFile_stat(FILEUNIT,material_localFileExt)) & ! no local material configuration present... + call IO_open_file(FILEUNIT,material_configFile) ! ... open material.config file + call homogenization_isostrain_init(FILEUNIT) + call homogenization_RGC_init(FILEUNIT) + close(FILEUNIT) !-------------------------------------------------------------------------------------------------- ! write description file for homogenization output - call IO_write_jobFile(fileunit,'outputHomogenization') + call IO_write_jobFile(FILEUNIT,'outputHomogenization') do p = 1,material_Nhomogenization i = homogenization_typeInstance(p) ! which instance of this homogenization type knownHomogenization = .true. ! assume valid @@ -156,16 +156,16 @@ subroutine homogenization_init() case default knownHomogenization = .false. end select - write(fileunit,'(/,a,/)') '['//trim(homogenization_name(p))//']' + write(FILEUNIT,'(/,a,/)') '['//trim(homogenization_name(p))//']' if (knownHomogenization) then - write(fileunit,'(a)') '(type)'//char(9)//trim(outputName) - write(fileunit,'(a,i4)') '(ngrains)'//char(9),homogenization_Ngrains(p) + write(FILEUNIT,'(a)') '(type)'//char(9)//trim(outputName) + write(FILEUNIT,'(a,i4)') '(ngrains)'//char(9),homogenization_Ngrains(p) do e = 1,homogenization_Noutput(p) - write(fileunit,'(a,i4)') trim(thisOutput(e,i))//char(9),thisSize(e,i) + write(FILEUNIT,'(a,i4)') trim(thisOutput(e,i))//char(9),thisSize(e,i) enddo endif enddo - close(fileunit) + close(FILEUNIT) !-------------------------------------------------------------------------------------------------- ! allocate and initialize global variables diff --git a/code/homogenization_RGC.f90 b/code/homogenization_RGC.f90 index e462cd593..ffbc89f2c 100644 --- a/code/homogenization_RGC.f90 +++ b/code/homogenization_RGC.f90 @@ -75,7 +75,7 @@ contains !-------------------------------------------------------------------------------------------------- !> @brief allocates all neccessary fields, reads information from material configuration file !-------------------------------------------------------------------------------------------------- -subroutine homogenization_RGC_init(myUnit) +subroutine homogenization_RGC_init(fileUnit) use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) use debug, only: & debug_level, & @@ -99,12 +99,13 @@ subroutine homogenization_RGC_init(myUnit) use material implicit none - integer(pInt), intent(in) :: myUnit !< file pointer to material configuration + integer(pInt), intent(in) :: fileUnit !< file pointer to material configuration 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 = '' + character(len=65536) :: & + tag = '', & + line = '' write(6,'(/,a)') ' <<<+- homogenization_'//HOMOGENIZATION_RGC_label//' init -+>>>' write(6,'(a)') ' $Id$' @@ -128,16 +129,18 @@ subroutine homogenization_RGC_init(myUnit) allocate(homogenization_RGC_orientation(3,3,mesh_maxNips,mesh_NcpElems)) homogenization_RGC_orientation = spread(spread(math_I3,3,mesh_maxNips),4,mesh_NcpElems) ! initialize to identity - rewind(myUnit) - - do while (trim(line) /= '#EOF#' .and. IO_lc(IO_getTag(line,'<','>')) /= material_partHomogenization) ! wind forward to - line = IO_read(myUnit) + rewind(fileUnit) + do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= material_partHomogenization) ! wind forward to + line = IO_read(fileUnit) enddo - do while (trim(line) /= '#EOF#') - line = IO_read(myUnit) + 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,'<','>') /= '') exit ! stop at next part + 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 output = 0_pInt ! reset output counter diff --git a/code/homogenization_isostrain.f90 b/code/homogenization_isostrain.f90 index 41727883c..e9e19c13a 100644 --- a/code/homogenization_isostrain.f90 +++ b/code/homogenization_isostrain.f90 @@ -29,28 +29,31 @@ module homogenization_isostrain implicit none private - integer(pInt), dimension(:), allocatable, public, protected :: & + 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 - integer(pInt), dimension(:), allocatable, private :: & + integer(pInt), dimension(:), allocatable, private :: & homogenization_isostrain_Ngrains enum, bind(c) - enumerator :: nconstituents, & - temperature, & - ipcoords, & - avgdefgrad, & - avgfirstpiola - enumerator :: parallel, & - average + enumerator :: undefined_ID, & + nconstituents_ID, & + temperature_ID, & + ipcoords_ID, & + avgdefgrad_ID, & + avgfirstpiola_ID end enum - integer(kind(nconstituents)), dimension(:,:), allocatable, private :: & + 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)), dimension(:), allocatable, private :: & + integer(kind(average_ID)), dimension(:), allocatable, private :: & homogenization_isostrain_mapping !< ID of each post result output @@ -65,7 +68,7 @@ contains !-------------------------------------------------------------------------------------------------- !> @brief allocates all neccessary fields, reads information from material configuration file !-------------------------------------------------------------------------------------------------- -subroutine homogenization_isostrain_init(myUnit) +subroutine homogenization_isostrain_init(fileUnit) 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, & @@ -74,16 +77,16 @@ subroutine homogenization_isostrain_init(myUnit) use material implicit none - integer(pInt), intent(in) :: myUnit + integer(pInt), intent(in) :: fileUnit integer(pInt), parameter :: MAXNCHUNKS = 2_pInt integer(pInt), dimension(1_pInt+2_pInt*MAXNCHUNKS) :: positions integer(pInt) :: & - section, i, j, output, mySize + section = 0_pInt, i, j, output, mySize integer :: & maxNinstance, k ! no pInt (stores a system dependen value from 'count' character(len=65536) :: & tag = '', & - line = '' ! to start initialized + line = '' write(6,'(/,a)') ' <<<+- homogenization_'//HOMOGENIZATION_ISOSTRAIN_label//' init -+>>>' write(6,'(a)') ' $Id$' @@ -93,32 +96,29 @@ subroutine homogenization_isostrain_init(myUnit) maxNinstance = count(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID) 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 = average + allocate(homogenization_isostrain_sizeState(maxNinstance), source=0_pInt) + allocate(homogenization_isostrain_sizePostResults(maxNinstance), source=0_pInt) + allocate(homogenization_isostrain_sizePostResult(maxval(homogenization_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)) - homogenization_isostrain_outputID = -1 + allocate(homogenization_isostrain_outputID(maxval(homogenization_Noutput),maxNinstance), & + source=undefined_ID) - rewind(myUnit) - section = 0_pInt - - do while (trim(line) /= '#EOF#' .and. IO_lc(IO_getTag(line,'<','>')) /= material_partHomogenization) ! wind forward to - line = IO_read(myUnit) + rewind(fileUnit) + do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= material_partHomogenization) ! wind forward to + line = IO_read(fileUnit) enddo - do while (trim(line) /= '#EOF#') - line = IO_read(myUnit) + 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,'<','>') /= '') exit ! stop at next part + 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 output = 0_pInt ! reset output counter @@ -134,15 +134,15 @@ subroutine homogenization_isostrain_init(myUnit) homogenization_isostrain_output(output,i) = IO_lc(IO_stringValue(line,positions,2_pInt)) select case(homogenization_isostrain_output(output,i)) case('nconstituents','ngrains') - homogenization_isostrain_outputID(output,i) = nconstituents + homogenization_isostrain_outputID(output,i) = nconstituents_ID case('temperature') - homogenization_isostrain_outputID(output,i) = temperature + homogenization_isostrain_outputID(output,i) = nconstituents_ID case('ipcoords') - homogenization_isostrain_outputID(output,i) = ipcoords + homogenization_isostrain_outputID(output,i) = nconstituents_ID case('avgdefgrad','avgf') - homogenization_isostrain_outputID(output,i) = avgdefgrad + homogenization_isostrain_outputID(output,i) = nconstituents_ID case('avgp','avgfirstpiola','avg1stpiola') - homogenization_isostrain_outputID(output,i) = avgfirstpiola + homogenization_isostrain_outputID(output,i) = nconstituents_ID case default mySize = 0_pInt end select @@ -151,11 +151,10 @@ subroutine homogenization_isostrain_init(myUnit) case ('mapping') select case(IO_lc(IO_stringValue(line,positions,2_pInt))) case ('parallel','sum') - homogenization_isostrain_mapping(i) = parallel + homogenization_isostrain_mapping(i) = nconstituents_ID case ('average','mean','avg') - homogenization_isostrain_mapping(i) = average + homogenization_isostrain_mapping(i) = nconstituents_ID case default - print*, 'There should be an error here' end select end select endif @@ -167,11 +166,11 @@ subroutine homogenization_isostrain_init(myUnit) do j = 1_pInt,maxval(homogenization_Noutput) select case(homogenization_isostrain_outputID(j,i)) - case(nconstituents, temperature) + case(nconstituents_ID, temperature_ID) mySize = 1_pInt - case(ipcoords) + case(ipcoords_ID) mySize = 3_pInt - case(avgdefgrad, avgfirstpiola) + case(avgdefgrad_ID, avgfirstpiola_ID) mySize = 9_pInt case default mySize = 0_pInt @@ -238,10 +237,10 @@ subroutine homogenization_isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P Ngrains = homogenization_Ngrains(mesh_element(3,el)) select case (homogenization_isostrain_mapping(homID)) - case (parallel) + case (parallel_ID) avgP = sum(P,3) dAvgPdAvgF = sum(dPdF,5) - case (average) + case (average_ID) avgP = sum(P,3) /real(Ngrains,pReal) dAvgPdAvgF = sum(dPdF,5)/real(Ngrains,pReal) end select @@ -285,19 +284,19 @@ pure function homogenization_isostrain_postResults(ip,el,avgP,avgF) do o = 1_pInt,homogenization_Noutput(mesh_element(3,el)) select case(homogenization_isostrain_outputID(o,homID)) - case (nconstituents) + case (nconstituents_ID) homogenization_isostrain_postResults(c+1_pInt) = real(homogenization_isostrain_Ngrains(homID),pReal) c = c + 1_pInt - case (temperature) + case (temperature_ID) homogenization_isostrain_postResults(c+1_pInt) = crystallite_temperature(ip,el) c = c + 1_pInt - case (avgdefgrad) + case (avgdefgrad_ID) homogenization_isostrain_postResults(c+1_pInt:c+9_pInt) = reshape(avgF,[9]) c = c + 9_pInt - case (avgfirstpiola) + case (avgfirstpiola_ID) homogenization_isostrain_postResults(c+1_pInt:c+9_pInt) = reshape(avgP,[9]) c = c + 9_pInt - case (ipcoords) + 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 diff --git a/code/material.f90 b/code/material.f90 index 4a078e659..ad297d66a 100644 --- a/code/material.f90 +++ b/code/material.f90 @@ -46,14 +46,21 @@ module material HOMOGENIZATION_RGC_label = 'rgc' enum, bind(c) - enumerator :: ELASTICITY_hooke_ID - enumerator :: PLASTICITY_none_ID, & + enumerator :: ELASTICITY_undefined_ID, & + ELASTICITY_hooke_ID + end enum + enum, bind(c) + enumerator :: PLASTICITY_undefined_ID, & + PLASTICITY_none_ID, & PLASTICITY_J2_ID, & PLASTICITY_phenopowerlaw_ID, & PLASTICITY_dislotwin_ID, & PLASTICITY_titanmod_ID, & PLASTICITY_nonlocal_ID - enumerator :: HOMOGENIZATION_isostrain_ID, & + end enum + enum, bind(c) + enumerator :: HOMOGENIZATION_undefined_ID, & + HOMOGENIZATION_isostrain_ID, & HOMOGENIZATION_RGC_ID end enum @@ -66,7 +73,7 @@ module material MATERIAL_partCrystallite = 'crystallite', & !< keyword for crystallite part MATERIAL_partPhase = 'phase' !< keyword for phase part - integer(kind(ELASTICITY_hooke_ID)), dimension(:), allocatable, public, protected :: & + integer(kind(ELASTICITY_undefined_ID)), dimension(:), allocatable, public, protected :: & phase_elasticity, & !< elasticity of each phase phase_plasticity, & !< plasticity of each phase homogenization_type !< type of each homogenization @@ -149,13 +156,16 @@ module material public :: & material_init, & - ELASTICITY_hooke_ID ,& + ELASTICITY_undefined_ID, & + ELASTICITY_hooke_ID, & + PLASTICITY_undefined_ID, & PLASTICITY_none_ID, & PLASTICITY_J2_ID, & PLASTICITY_phenopowerlaw_ID, & PLASTICITY_dislotwin_ID, & PLASTICITY_titanmod_ID, & PLASTICITY_nonlocal_ID, & + HOMOGENIZATION_undefined_ID, & HOMOGENIZATION_isostrain_ID, & HOMOGENIZATION_RGC_ID @@ -189,7 +199,7 @@ subroutine material_init debug_levelExtensive implicit none - integer(pInt), parameter :: fileunit = 200_pInt + integer(pInt), parameter :: FILEUNIT = 200_pInt integer(pInt) :: m,c,h, myDebug myDebug = debug_level(debug_material) @@ -198,19 +208,19 @@ subroutine material_init write(6,'(a16,a)') ' Current time : ',IO_timeStamp() #include "compilation_info.f90" - if (.not. IO_open_jobFile_stat(fileunit,material_localFileExt)) & ! no local material configuration present... - call IO_open_file(fileunit,material_configFile) ! ...open material.config file - call material_parseHomogenization(fileunit,material_partHomogenization) + if (.not. IO_open_jobFile_stat(FILEUNIT,material_localFileExt)) & ! no local material configuration present... + call IO_open_file(FILEUNIT,material_configFile) ! ...open material.config file + call material_parseHomogenization(FILEUNIT,material_partHomogenization) if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Homogenization parsed' - call material_parseMicrostructure(fileunit,material_partMicrostructure) + call material_parseMicrostructure(FILEUNIT,material_partMicrostructure) if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Microstructure parsed' - call material_parseCrystallite(fileunit,material_partCrystallite) + call material_parseCrystallite(FILEUNIT,material_partCrystallite) if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Crystallite parsed' - call material_parseTexture(fileunit,material_partTexture) + call material_parseTexture(FILEUNIT,material_partTexture) if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Texture parsed' - call material_parsePhase(fileunit,material_partPhase) + call material_parsePhase(FILEUNIT,material_partPhase) if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Phase parsed' - close(fileunit) + close(FILEUNIT) do m = 1_pInt,material_Nmicrostructure if(microstructure_crystallite(m) < 1_pInt .or. & @@ -262,7 +272,7 @@ end subroutine material_init !-------------------------------------------------------------------------------------------------- !> @brief parses the homogenization part in the material configuration file !-------------------------------------------------------------------------------------------------- -subroutine material_parseHomogenization(myFile,myPart) +subroutine material_parseHomogenization(fileUnit,myPart) use IO, only: & IO_read, & IO_globalTagInPart, & @@ -274,25 +284,25 @@ subroutine material_parseHomogenization(myFile,myPart) IO_isBlank, & IO_stringValue, & IO_intValue, & - IO_stringPos + IO_stringPos, & + IO_EOF use mesh, only: & mesh_element implicit none character(len=*), intent(in) :: myPart - integer(pInt), intent(in) :: myFile + integer(pInt), intent(in) :: fileUnit - integer(pInt), parameter :: maxNchunks = 2_pInt + integer(pInt), parameter :: MAXNCHUNKS = 2_pInt - integer(pInt), dimension(1+2*maxNchunks) :: positions + integer(pInt), dimension(1+2*MAXNCHUNKS) :: positions integer(pInt) :: Nsections, section, s - character(len=65536) :: tag - character(len=65536) :: line + character(len=65536) :: & + tag, line logical :: echo - - echo = IO_globalTagInPart(myFile,myPart,'/echo/') - - Nsections = IO_countSections(myFile,myPart) + + echo = IO_globalTagInPart(fileUnit,myPart,'/echo/') + Nsections = IO_countSections(fileUnit,myPart) material_Nhomogenization = Nsections if (Nsections < 1_pInt) call IO_error(160_pInt,ext_msg=myPart) @@ -304,28 +314,31 @@ subroutine material_parseHomogenization(myFile,myPart) allocate(homogenization_active(Nsections)); homogenization_active = .false. forall (s = 1_pInt:Nsections) homogenization_active(s) = any(mesh_element(3,:) == s) ! current homogenization used in model? Homogenization view, maximum operations depend on maximum number of homog schemes - homogenization_Noutput = IO_countTagInPart(myFile,myPart,'(output)',Nsections) + homogenization_Noutput = IO_countTagInPart(fileUnit,myPart,'(output)',Nsections) - rewind(myFile) - line = '' - section = 0_pInt - - do while (trim(line) /= '#EOF#' .and. IO_lc(IO_getTag(line,'<','>')) /= myPart) ! wind forward to myPart - line = IO_read(myFile) + rewind(fileUnit) + line = '' ! to have it initialized + section = 0_pInt ! - " - + do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= myPart) ! wind forward to myPart + line = IO_read(fileUnit) enddo if (echo) write(6,'(/,1x,a)') trim(line) ! echo part header - do while (trim(line) /= '#EOF#') - line = IO_read(myFile) + + do while (trim(line) /= IO_EOF) ! read through sections of material part + line = IO_read(fileUnit) if (IO_isBlank(line)) cycle ! skip empty lines - if (IO_getTag(line,'<','>') /= '') exit ! stop at next part + if (IO_getTag(line,'<','>') /= '') then ! stop at next part + line = IO_read(fileUnit, .true.) ! reset IO_read + exit + endif if (echo) write(6,'(2x,a)') trim(line) ! echo back read lines if (IO_getTag(line,'[',']') /= '') then ! next section section = section + 1_pInt - homogenization_name(section) = IO_getTag(line,'[',']') + homogenization_name(section) = IO_getTag(line,'[',']') endif if (section > 0_pInt) then - 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 ('type') @@ -353,7 +366,7 @@ end subroutine material_parseHomogenization !-------------------------------------------------------------------------------------------------- !> @brief parses the microstructure part in the material configuration file !-------------------------------------------------------------------------------------------------- -subroutine material_parseMicrostructure(myFile,myPart) +subroutine material_parseMicrostructure(fileUnit,myPart) use IO use mesh, only: & mesh_element, & @@ -361,19 +374,19 @@ subroutine material_parseMicrostructure(myFile,myPart) implicit none character(len=*), intent(in) :: myPart - integer(pInt), intent(in) :: myFile + integer(pInt), intent(in) :: fileUnit - integer(pInt), parameter :: maxNchunks = 7_pInt + integer(pInt), parameter :: MAXNCHUNKS = 7_pInt - integer(pInt), dimension(1_pInt+2_pInt*maxNchunks) :: positions + integer(pInt), dimension(1_pInt+2_pInt*MAXNCHUNKS) :: positions integer(pInt) :: Nsections, section, constituent, e, i - character(len=65536) :: tag - character(len=65536) :: line + character(len=65536) :: & + tag, line logical :: echo - echo = IO_globalTagInPart(myFile,myPart,'/echo/') + echo = IO_globalTagInPart(fileUnit,myPart,'/echo/') - Nsections = IO_countSections(myFile,myPart) + Nsections = IO_countSections(fileUnit,myPart) material_Nmicrostructure = Nsections if (Nsections < 1_pInt) call IO_error(160_pInt,ext_msg=myPart) @@ -385,9 +398,9 @@ subroutine material_parseMicrostructure(myFile,myPart) forall (e = 1_pInt:mesh_NcpElems) microstructure_active(mesh_element(4,e)) = .true. ! current microstructure used in model? Elementwise view, maximum N operations for N elements - microstructure_Nconstituents = IO_countTagInPart(myFile,myPart,'(constituent)',Nsections) + microstructure_Nconstituents = IO_countTagInPart(fileUnit,myPart,'(constituent)',Nsections) microstructure_maxNconstituents = maxval(microstructure_Nconstituents) - microstructure_elemhomo = IO_spotTagInPart(myFile,myPart,'/elementhomogeneous/',Nsections) + microstructure_elemhomo = IO_spotTagInPart(fileUnit,myPart,'/elementhomogeneous/',Nsections) allocate(microstructure_phase (microstructure_maxNconstituents,Nsections)) microstructure_phase = 0_pInt @@ -396,20 +409,22 @@ subroutine material_parseMicrostructure(myFile,myPart) allocate(microstructure_fraction(microstructure_maxNconstituents,Nsections)) microstructure_fraction = 0.0_pReal - rewind(myFile) + rewind(fileUnit) line = '' ! to have it initialized section = 0_pInt ! - " - constituent = 0_pInt ! - " - - - do while (trim(line) /= '#EOF#' .and. IO_lc(IO_getTag(line,'<','>')) /= myPart) ! wind forward to myPart - line = IO_read(myFile) + do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= myPart) ! wind forward to myPart + line = IO_read(fileUnit) enddo if (echo) write(6,'(/,1x,a)') trim(line) ! echo part header - do while (trim(line) /= '#EOF#') - line = IO_read(myFile) + do while (trim(line) /= IO_EOF) ! read through sections of material part + line = IO_read(fileUnit) if (IO_isBlank(line)) cycle ! skip empty lines - if (IO_getTag(line,'<','>') /= '') exit ! stop at next part + if (IO_getTag(line,'<','>') /= '') then ! stop at next part + line = IO_read(fileUnit, .true.) ! reset IO_read + exit + endif if (echo) write(6,'(2x,a)') trim(line) ! echo back read lines if (IO_getTag(line,'[',']') /= '') then ! next section section = section + 1_pInt @@ -417,7 +432,7 @@ subroutine material_parseMicrostructure(myFile,myPart) microstructure_name(section) = IO_getTag(line,'[',']') endif if (section > 0_pInt) then - 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 ('crystallite') @@ -445,7 +460,7 @@ end subroutine material_parseMicrostructure !-------------------------------------------------------------------------------------------------- !> @brief parses the crystallite part in the material configuration file !-------------------------------------------------------------------------------------------------- -subroutine material_parseCrystallite(myFile,myPart) +subroutine material_parseCrystallite(fileUnit,myPart) use IO, only: & IO_read, & IO_countSections, & @@ -454,41 +469,44 @@ subroutine material_parseCrystallite(myFile,myPart) IO_globalTagInPart, & IO_getTag, & IO_lc, & - IO_isBlank + IO_isBlank, & + IO_EOF implicit none character(len=*), intent(in) :: myPart - integer(pInt), intent(in) :: myFile + integer(pInt), intent(in) :: fileUnit integer(pInt) :: Nsections, & section character(len=65536) :: line logical :: echo - echo = IO_globalTagInPart(myFile,myPart,'/echo/') + echo = IO_globalTagInPart(fileUnit,myPart,'/echo/') - Nsections = IO_countSections(myFile,myPart) + Nsections = IO_countSections(fileUnit,myPart) material_Ncrystallite = Nsections if (Nsections < 1_pInt) call IO_error(160_pInt,ext_msg=myPart) allocate(crystallite_name(Nsections)); crystallite_name = '' allocate(crystallite_Noutput(Nsections)); crystallite_Noutput = 0_pInt - crystallite_Noutput = IO_countTagInPart(myFile,myPart,'(output)',Nsections) + crystallite_Noutput = IO_countTagInPart(fileUnit,myPart,'(output)',Nsections) - rewind(myFile) - line = '' - section = 0_pInt - - do while (trim(line) /= '#EOF#' .and. IO_lc(IO_getTag(line,'<','>')) /= myPart) ! wind forward to myPart - line = IO_read(myFile) + rewind(fileUnit) + line = '' ! to have it initialized + section = 0_pInt ! - " - + do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= myPart) ! wind forward to myPart + line = IO_read(fileUnit) enddo if (echo) write(6,'(/,1x,a)') trim(line) ! echo part header - do while (trim(line) /= '#EOF#') - line = IO_read(myFile) + do while (trim(line) /= IO_EOF) ! read through sections of material part + line = IO_read(fileUnit) if (IO_isBlank(line)) cycle ! skip empty lines - if (IO_getTag(line,'<','>') /= '') exit ! stop at next part + if (IO_getTag(line,'<','>') /= '') then ! stop at next part + line = IO_read(fileUnit, .true.) ! reset IO_read + exit + endif if (echo) write(6,'(2x,a)') trim(line) ! echo back read lines if (IO_getTag(line,'[',']') /= '') then ! next section section = section + 1_pInt @@ -502,7 +520,7 @@ end subroutine material_parseCrystallite !-------------------------------------------------------------------------------------------------- !> @brief parses the phase part in the material configuration file !-------------------------------------------------------------------------------------------------- -subroutine material_parsePhase(myFile,myPart) +subroutine material_parsePhase(fileUnit,myPart) use IO, only: & IO_read, & IO_globalTagInPart, & @@ -514,23 +532,24 @@ subroutine material_parsePhase(myFile,myPart) IO_lc, & IO_isBlank, & IO_stringValue, & - IO_stringPos + IO_stringPos, & + IO_EOF implicit none character(len=*), intent(in) :: myPart - integer(pInt), intent(in) :: myFile + integer(pInt), intent(in) :: fileUnit - integer(pInt), parameter :: maxNchunks = 2_pInt + integer(pInt), parameter :: MAXNCHUNKS = 2_pInt - integer(pInt), dimension(1+2*maxNchunks) :: positions - integer(pInt) Nsections, section - character(len=65536) :: tag - character(len=65536) :: line + integer(pInt), dimension(1+2*MAXNCHUNKS) :: positions + integer(pInt) :: Nsections, section + character(len=65536) :: & + tag,line logical :: echo - echo = IO_globalTagInPart(myFile,myPart,'/echo/') + echo = IO_globalTagInPart(fileUnit,myPart,'/echo/') - Nsections = IO_countSections(myFile,myPart) + Nsections = IO_countSections(fileUnit,myPart) material_Nphase = Nsections if (Nsections < 1_pInt) call IO_error(160_pInt,ext_msg=myPart) @@ -542,29 +561,31 @@ subroutine material_parsePhase(myFile,myPart) allocate(phase_Noutput(Nsections)); phase_Noutput = 0_pInt allocate(phase_localPlasticity(Nsections)); phase_localPlasticity = .false. - phase_Noutput = IO_countTagInPart(myFile,myPart,'(output)',Nsections) - phase_localPlasticity = .not. IO_spotTagInPart(myFile,myPart,'/nonlocal/',Nsections) + phase_Noutput = IO_countTagInPart(fileUnit,myPart,'(output)',Nsections) + phase_localPlasticity = .not. IO_spotTagInPart(fileUnit,myPart,'/nonlocal/',Nsections) - rewind(myFile) - line = '' - section = 0_pInt - - do while (trim(line) /= '#EOF#' .and. IO_lc(IO_getTag(line,'<','>')) /= myPart) ! wind forward to myPart - line = IO_read(myFile) + rewind(fileUnit) + line = '' ! to have it initialized + section = 0_pInt ! - " - + do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= myPart) ! wind forward to myPart + line = IO_read(fileUnit) enddo if (echo) write(6,'(/,1x,a)') trim(line) ! echo part header - do while (trim(line) /= '#EOF#') - line = IO_read(myFile) + do while (trim(line) /= IO_EOF) ! read through sections of material part + line = IO_read(fileUnit) if (IO_isBlank(line)) cycle ! skip empty lines - if (IO_getTag(line,'<','>') /= '') exit ! stop at next part + if (IO_getTag(line,'<','>') /= '') then ! stop at next part + line = IO_read(fileUnit, .true.) ! reset IO_read + exit + endif if (echo) write(6,'(2x,a)') trim(line) ! echo back read lines if (IO_getTag(line,'[',']') /= '') then ! next section section = section + 1_pInt phase_name(section) = IO_getTag(line,'[',']') endif if (section > 0_pInt) then - 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 ('elasticity') @@ -603,7 +624,7 @@ end subroutine material_parsePhase !-------------------------------------------------------------------------------------------------- !> @brief parses the texture part in the material configuration file !-------------------------------------------------------------------------------------------------- -subroutine material_parseTexture(myFile,myPart) +subroutine material_parseTexture(fileUnit,myPart) use IO, only: & IO_read, & IO_globalTagInPart, & @@ -616,7 +637,8 @@ subroutine material_parseTexture(myFile,myPart) IO_isBlank, & IO_floatValue, & IO_stringValue, & - IO_stringPos + IO_stringPos, & + IO_EOF use math, only: & inRad, & math_sampleRandomOri, & @@ -625,19 +647,19 @@ subroutine material_parseTexture(myFile,myPart) implicit none character(len=*), intent(in) :: myPart - integer(pInt), intent(in) :: myFile + integer(pInt), intent(in) :: fileUnit - integer(pInt), parameter :: maxNchunks = 13_pInt + integer(pInt), parameter :: MAXNCHUNKS = 13_pInt - integer(pInt), dimension(1+2*maxNchunks) :: positions + integer(pInt), dimension(1+2*MAXNCHUNKS) :: positions integer(pInt) :: Nsections, section, gauss, fiber, j character(len=65536) :: tag character(len=65536) :: line logical :: echo - echo = IO_globalTagInPart(myFile,myPart,'/echo/') + echo = IO_globalTagInPart(fileUnit,myPart,'/echo/') - Nsections = IO_countSections(myFile,myPart) + Nsections = IO_countSections(fileUnit,myPart) material_Ntexture = Nsections if (Nsections < 1_pInt) call IO_error(160_pInt,ext_msg=myPart) @@ -647,9 +669,9 @@ subroutine material_parseTexture(myFile,myPart) allocate(texture_Ngauss(Nsections)); texture_Ngauss = 0_pInt allocate(texture_Nfiber(Nsections)); texture_Nfiber = 0_pInt - texture_Ngauss = IO_countTagInPart(myFile,myPart,'(gauss)', Nsections) + & - IO_countTagInPart(myFile,myPart,'(random)',Nsections) - texture_Nfiber = IO_countTagInPart(myFile,myPart,'(fiber)', Nsections) + texture_Ngauss = IO_countTagInPart(fileUnit,myPart,'(gauss)', Nsections) + & + IO_countTagInPart(fileUnit,myPart,'(random)',Nsections) + texture_Nfiber = IO_countTagInPart(fileUnit,myPart,'(fiber)', Nsections) texture_maxNgauss = maxval(texture_Ngauss) texture_maxNfiber = maxval(texture_Nfiber) allocate(texture_Gauss (5,texture_maxNgauss,Nsections)); texture_Gauss = 0.0_pReal @@ -659,21 +681,23 @@ subroutine material_parseTexture(myFile,myPart) texture_transformation(1:3,1:3,j) = math_I3 enddo - rewind(myFile) + rewind(fileUnit) line = '' ! to have in initialized section = 0_pInt ! - " - gauss = 0_pInt ! - " - fiber = 0_pInt ! - " - - - do while (trim(line) /= '#EOF#' .and. IO_lc(IO_getTag(line,'<','>')) /= myPart) ! wind forward to myPart - line = IO_read(myFile) + do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= myPart) ! wind forward to myPart + line = IO_read(fileUnit) enddo if (echo) write(6,'(/,1x,a)') trim(line) ! echo part header - do while (trim(line) /= '#EOF#') - line = IO_read(myFile) + do while (trim(line) /= IO_EOF) + line = IO_read(fileUnit) if (IO_isBlank(line)) cycle ! skip empty lines - if (IO_getTag(line,'<','>') /= '') exit ! stop at next part + if (IO_getTag(line,'<','>') /= '') then ! stop at next part + line = IO_read(fileUnit, .true.) ! reset IO_read + exit + endif if (echo) write(6,'(2x,a)') trim(line) ! echo back read lines if (IO_getTag(line,'[',']') /= '') then ! next section section = section + 1_pInt @@ -682,7 +706,7 @@ subroutine material_parseTexture(myFile,myPart) texture_name(section) = IO_getTag(line,'[',']') endif if (section > 0_pInt) then - positions = IO_stringPos(line,maxNchunks) + positions = IO_stringPos(line,MAXNCHUNKS) tag = IO_lc(IO_stringValue(line,positions,1_pInt)) ! extract key textureType: select case(tag)