From e555ce4827fe7e4315127ecf444a3c1bbca378c7 Mon Sep 17 00:00:00 2001 From: Luv Sharma Date: Thu, 21 Aug 2014 17:48:20 +0000 Subject: [PATCH] started introducing new state structure in homogenisation --- code/CPFEM.f90 | 29 +++- code/Makefile | 5 + code/homogenization.f90 | 85 +++++++++++- code/homogenization_RGC.f90 | 216 +++++++++++++++++++++++++++++- code/homogenization_isostrain.f90 | 31 ++++- code/homogenization_none.f90 | 35 ++++- code/material.f90 | 52 ++++++- code/prec.f90 | 9 ++ 8 files changed, 438 insertions(+), 24 deletions(-) diff --git a/code/CPFEM.f90 b/code/CPFEM.f90 index 281a643d9..771106a15 100644 --- a/code/CPFEM.f90 +++ b/code/CPFEM.f90 @@ -147,6 +147,10 @@ subroutine CPFEM_init use material, only: & homogenization_maxNgrains, & material_phase, & +#ifdef NEWSTATE + homogState, & + mappingHomogenization, & +#endif phase_plasticity, & plasticState use crystallite, only: & @@ -157,8 +161,10 @@ subroutine CPFEM_init crystallite_Tstar0_v, & crystallite_localPlasticity use homogenization, only: & - homogenization_sizeState, & - homogenization_state0 +#ifndef NEWSTATE + homogenization_state0, & +#endif + homogenization_sizeState implicit none integer(pInt) :: i,j,k,l,m,ph @@ -222,7 +228,11 @@ subroutine CPFEM_init do k = 1,mesh_NcpElems; do j = 1,mesh_maxNips do l = 1,homogenization_sizeState(j,k) m = m+1_pInt +#ifdef NEWSTATE + read(777,rec=m) homogState(mappingHomogenization(2,j,k))%state0(l,mappingHomogenization(1,j,k)) +#else read(777,rec=m) homogenization_state0(j,k)%p(l) +#endif enddo enddo; enddo @@ -304,6 +314,10 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature, dt, elFE, ip) microstructure_elemhomo, & plasticState, & damageState, & +#ifdef NEWSTATE + homogState, & + mappingHomogenization, & +#endif thermalState, & mappingConstitutive, & material_phase, & @@ -322,8 +336,10 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature, dt, elFE, ip) crystallite_temperature use homogenization, only: & homogenization_sizeState, & +#ifndef NEWSTATE homogenization_state, & homogenization_state0, & +#endif materialpoint_F, & materialpoint_F0, & materialpoint_P, & @@ -417,7 +433,12 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature, dt, elFE, ip) do k = 1,mesh_NcpElems do j = 1,mesh_maxNips if (homogenization_sizeState(j,k) > 0_pInt) & +#ifdef NEWSTATE + homogState(mappingHomogenization(2,j,k))%state0(:,mappingHomogenization(1,j,k)) = & + homogState(mappingHomogenization(2,j,k))%state(:,mappingHomogenization(1,j,k)) ! internal state of homogenization scheme +#else homogenization_state0(j,k)%p = homogenization_state(j,k)%p ! internal state of homogenization scheme +#endif enddo enddo !$OMP END PARALLEL DO @@ -468,7 +489,11 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature, dt, elFE, ip) do k = 1,mesh_NcpElems; do j = 1,mesh_maxNips do l = 1,homogenization_sizeState(j,k) m = m+1_pInt +#ifdef NEWSTATE + write(777,rec=m) homogState(mappingHomogenization(2,j,k))%state0(l,mappingHomogenization(1,j,k)) +#else write(777,rec=m) homogenization_state0(j,k)%p(l) +#endif enddo enddo; enddo close (777) diff --git a/code/Makefile b/code/Makefile index 2f09388f4..7c3079f74 100644 --- a/code/Makefile +++ b/code/Makefile @@ -117,6 +117,11 @@ RUN_PATH :=$(RUN_PATH),-rpath,$(HDF5_ROOT)/lib64,-rpath,$(HDF5_ROOT)/lib INCLUDE_DIRS +=-I$(HDF5_ROOT)/include -DHDF endif +#newstate +ifeq "$(STATE)" "NEWH" +INCLUDE_DIRS +=-DNEWSTATE +endif + ifdef STANDARD_CHECK STANDARD_CHECK_ifort =$(STANDARD_CHECK) STANDARD_CHECK_gfortran =$(STANDARD_CHECK) diff --git a/code/homogenization.f90 b/code/homogenization.f90 index d6223020b..e76d9071d 100644 --- a/code/homogenization.f90 +++ b/code/homogenization.f90 @@ -11,14 +11,16 @@ module homogenization pInt, & pReal, & p_vec - + !-------------------------------------------------------------------------------------------------- ! General variables for the homogenization at a material point implicit none private +#ifndef NEWSTATE type(p_vec), dimension(:,:), allocatable, public :: & homogenization_state0 !< pointer array to homogenization state at start of FE increment - real(pReal), dimension(:,:,:,:), allocatable, public :: & +#endif + real(pReal), dimension(:,:,:,:), allocatable, public :: & materialpoint_F0, & !< def grad of IP at start of FE increment materialpoint_F, & !< def grad of IP to be reached at end of FE increment materialpoint_P !< first P--K stress of IP @@ -26,10 +28,11 @@ module homogenization materialpoint_dPdF !< tangent of first P--K stress at IP real(pReal), dimension(:,:,:), allocatable, public :: & materialpoint_results !< results array of material point - +#ifndef NEWSTATE type(p_vec), dimension(:,:), allocatable, public, protected :: & homogenization_state !< pointer array to current homogenization state (end of converged time step) - integer(pInt), dimension(:,:), allocatable, public, protected :: & +#endif + integer(pInt), dimension(:,:), allocatable, public, protected :: & homogenization_sizeState !< size of state array per grain integer(pInt), public, protected :: & materialpoint_sizeResults, & @@ -175,9 +178,11 @@ subroutine homogenization_init() !-------------------------------------------------------------------------------------------------- ! allocate and initialize global variables +#ifndef NEWSTATE allocate(homogenization_state0(mesh_maxNips,mesh_NcpElems)) allocate(homogenization_subState0(mesh_maxNips,mesh_NcpElems)) allocate(homogenization_state(mesh_maxNips,mesh_NcpElems)) +#endif allocate(homogenization_sizeState(mesh_maxNips,mesh_NcpElems), source=0_pInt) allocate(homogenization_sizePostResults(mesh_maxNips,mesh_NcpElems), source=0_pInt) allocate(materialpoint_heat(mesh_maxNips,mesh_NcpElems), source=0.0_pReal) @@ -207,18 +212,34 @@ subroutine homogenization_init() #endif select case(homogenization_type(mesh_element(3,e))) case (HOMOGENIZATION_none_ID) +#ifdef NEWSTATE + homogenization_sizePostResults(i,e) = homogState(mappingHomogenization(2,i,e))%sizePostResults +#else homogenization_sizePostResults(i,e) = 0_pInt +#endif case (HOMOGENIZATION_ISOSTRAIN_ID) +#ifdef NEWSTATE + homogenization_sizePostResults(i,e) = homogState(mappingHomogenization(2,i,e))%sizePostResults +#else homogenization_sizePostResults(i,e) = homogenization_isostrain_sizePostResults(myInstance) +#endif case (HOMOGENIZATION_RGC_ID) if (homogenization_RGC_sizeState(myInstance) > 0_pInt) then +#ifdef NEWSTATE + homogenization_sizeState(i,e) = homogState(mappingHomogenization(2,i,e))%sizeState +#else allocate(homogenization_state0(i,e)%p(homogenization_RGC_sizeState(myInstance))) allocate(homogenization_subState0(i,e)%p(homogenization_RGC_sizeState(myInstance))) allocate(homogenization_state(i,e)%p(homogenization_RGC_sizeState(myInstance))) homogenization_state0(i,e)%p = 0.0_pReal homogenization_sizeState(i,e) = homogenization_RGC_sizeState(myInstance) +#endif endif +#ifdef NEWSTATE + homogenization_sizePostResults(i,e) = homogState(mappingHomogenization(2,i,e))%sizePostResults +#else homogenization_sizePostResults(i,e) = homogenization_RGC_sizePostResults(myInstance) +#endif end select enddo IpLooping enddo elementLooping @@ -247,9 +268,11 @@ subroutine homogenization_init() write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" if (iand(debug_level(debug_homogenization), debug_levelBasic) /= 0_pInt) then +#ifndef NEWSTATE write(6,'(a32,1x,7(i8,1x))') 'homogenization_state0: ', shape(homogenization_state0) write(6,'(a32,1x,7(i8,1x))') 'homogenization_subState0: ', shape(homogenization_subState0) write(6,'(a32,1x,7(i8,1x))') 'homogenization_state: ', shape(homogenization_state) +#endif write(6,'(a32,1x,7(i8,1x))') 'homogenization_sizeState: ', shape(homogenization_sizeState) write(6,'(a32,1x,7(i8,1x),/)') 'homogenization_sizePostResults: ', shape(homogenization_sizePostResults) write(6,'(a32,1x,7(i8,1x))') 'materialpoint_dPdF: ', shape(materialpoint_dPdF) @@ -301,6 +324,10 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt) plasticState, & damageState, & thermalState, & +#ifdef NEWSTATE + homogState, & + mappingHomogenization, & +#endif mappingConstitutive, & homogenization_Ngrains @@ -388,8 +415,14 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt) materialpoint_converged(i,e) = .false. ! pretend failed step of twice the required size materialpoint_requested(i,e) = .true. ! everybody requires calculation endforall +#ifdef NEWSTATE + forall(i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), homogenization_sizeState(i,e) > 0_pInt) & + homogState(mappingHomogenization(2,i,e))%subState0(:,mappingHomogenization(1,i,e)) = & + homogState(mappingHomogenization(2,i,e))%State0(:,mappingHomogenization(1,i,e)) ! ...internal homogenization state +#else forall(i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), homogenization_sizeState(i,e) > 0_pInt) & homogenization_subState0(i,e)%p = homogenization_state0(i,e)%p ! ...internal homogenization state +#endif enddo NiterationHomog = 0_pInt @@ -437,7 +470,12 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt) thermalState(mappingConstitutive(2,g,i,e))%state(:,mappingConstitutive(1,g,i,e)) end forall if (homogenization_sizeState(i,e) > 0_pInt) & - homogenization_subState0(i,e)%p = homogenization_state(i,e)%p ! ...internal state of homog scheme +#ifdef NEWSTATE + homogState(mappingHomogenization(2,i,e))%subState0(:,mappingHomogenization(1,i,e)) = & + homogState(mappingHomogenization(2,i,e))%state(:,mappingHomogenization(1,i,e)) +#else + homogenization_subState0(i,e)%p = homogenization_state(i,e)%p ! ...internal state of homog scheme +#endif materialpoint_subF0(1:3,1:3,i,e) = materialpoint_subF(1:3,1:3,i,e) ! ...def grad !$OMP FLUSH(materialpoint_subF0) elseif (materialpoint_requested(i,e)) then steppingNeeded ! already at final time (??) @@ -491,7 +529,12 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt) thermalState(mappingConstitutive(2,g,i,e))%partionedState0(:,mappingConstitutive(1,g,i,e)) end forall if (homogenization_sizeState(i,e) > 0_pInt) & +#ifdef NEWSTATE + homogState(mappingHomogenization(2,i,e))%state(:,mappingHomogenization(1,i,e)) = & + homogState(mappingHomogenization(2,i,e))%subState0(:,mappingHomogenization(1,i,e)) +#else homogenization_state(i,e)%p = homogenization_subState0(i,e)%p ! ...internal state of homog scheme +#endif endif endif converged @@ -703,12 +746,21 @@ subroutine homogenization_partitionDeformation(ip,el) materialpoint_subF(1:3,1:3,ip,el),& el) case (HOMOGENIZATION_RGC_ID) chosenHomogenization +#ifdef NEWSTATE + call homogenization_RGC_partitionDeformation(& + crystallite_partionedF(1:3,1:3,1:homogenization_maxNgrains,ip,el), & + materialpoint_subF(1:3,1:3,ip,el),& + ip, & + el) +#else call homogenization_RGC_partitionDeformation(& crystallite_partionedF(1:3,1:3,1:homogenization_maxNgrains,ip,el), & materialpoint_subF(1:3,1:3,ip,el),& homogenization_state(ip,el), & ip, & el) +#endif + end select chosenHomogenization end subroutine homogenization_partitionDeformation @@ -743,7 +795,17 @@ function homogenization_updateState(ip,el) case (HOMOGENIZATION_RGC_ID) chosenHomogenization homogenization_updateState = & - homogenization_RGC_updateState( homogenization_state(ip,el), & +#ifdef NEWSTATE + homogenization_RGC_updateState(crystallite_P(1:3,1:3,1:homogenization_maxNgrains,ip,el), & + crystallite_partionedF(1:3,1:3,1:homogenization_maxNgrains,ip,el), & + crystallite_partionedF0(1:3,1:3,1:homogenization_maxNgrains,ip,el),& + materialpoint_subF(1:3,1:3,ip,el),& + materialpoint_subdt(ip,el), & + crystallite_dPdF(1:3,1:3,1:3,1:3,1:homogenization_maxNgrains,ip,el), & + ip, & + el) +#else + homogenization_RGC_updateState(homogenization_state(ip,el), & homogenization_subState0(ip,el), & crystallite_P(1:3,1:3,1:homogenization_maxNgrains,ip,el), & crystallite_partionedF(1:3,1:3,1:homogenization_maxNgrains,ip,el), & @@ -753,6 +815,7 @@ function homogenization_updateState(ip,el) crystallite_dPdF(1:3,1:3,1:3,1:3,1:homogenization_maxNgrains,ip,el), & ip, & el) +#endif case default chosenHomogenization homogenization_updateState = .true. end select chosenHomogenization @@ -869,12 +932,22 @@ function homogenization_postResults(ip,el) materialpoint_P(1:3,1:3,ip,el), & materialpoint_F(1:3,1:3,ip,el)) case (HOMOGENIZATION_RGC_ID) chosenHomogenization +#ifdef NEWSTATE + homogenization_postResults = homogenization_RGC_postResults(& + ip, & + el, & + materialpoint_P(1:3,1:3,ip,el), & + materialpoint_F(1:3,1:3,ip,el)) + +#else homogenization_postResults = homogenization_RGC_postResults(& homogenization_state(ip,el),& ip, & el, & materialpoint_P(1:3,1:3,ip,el), & materialpoint_F(1:3,1:3,ip,el)) + +#endif end select chosenHomogenization end function homogenization_postResults diff --git a/code/homogenization_RGC.f90 b/code/homogenization_RGC.f90 index 5e2a3979b..1886a07c8 100644 --- a/code/homogenization_RGC.f90 +++ b/code/homogenization_RGC.f90 @@ -75,6 +75,9 @@ contains !-------------------------------------------------------------------------------------------------- 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 prec, only: & + pReal, & + pInt use debug, only: & debug_level, & debug_homogenization, & @@ -100,6 +103,13 @@ subroutine homogenization_RGC_init(fileUnit) 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 +#ifdef NEWSTATE + integer :: & + homog, & + NofMyHomog, & + instance, & + sizeHState +#endif integer(pInt) :: section=0_pInt, maxNinstance, i,j,e, output=-1_pInt, mySize, myInstance character(len=65536) :: & tag = '', & @@ -238,7 +248,7 @@ subroutine homogenization_RGC_init(fileUnit) write(6,'(a25,3(1x,e10.3))') 'cluster orientation: ',(homogenization_RGC_angles(j,i),j=1_pInt,3_pInt) enddo endif - + do i = 1_pInt,maxNinstance do j = 1_pInt,maxval(homogenization_Noutput) select case(homogenization_RGC_outputID(j,i)) @@ -267,14 +277,36 @@ subroutine homogenization_RGC_init(fileUnit) + 8_pInt ! (1) Average constitutive work, (2-4) Overall mismatch, (5) Average penalty energy, ! (6) Volume discrepancy, (7) Avg relaxation rate component, (8) Max relaxation rate component enddo +#ifdef NEWSTATE + initializeInstances: do homog = 1_pInt, material_Nhomogenization + + myhomog: if (homogenization_type(homog) == HOMOGENIZATION_RGC_ID) then + NofMyHomog = count(material_homog == homog) +! instance = phase_plasticityInstance(phase) +! allocate homogenization state arrays + sizeHState = homogenization_RGC_sizeState(homog) + homogState(homog)%sizeState = sizeHState + homogState(homog)%sizePostResults = homogenization_RGC_sizePostResults(homog) + allocate(homogState(homog)%state0 ( sizeHState,NofMyHomog), source=0.0_pReal) + allocate(homogState(homog)%subState0 ( sizeHState,NofMyHomog), source=0.0_pReal) + allocate(homogState(homog)%state ( sizeHState,NofMyHomog), source=0.0_pReal) + + endif myhomog + enddo initializeInstances +#endif end subroutine homogenization_RGC_init !-------------------------------------------------------------------------------------------------- !> @brief partitions the deformation gradient onto the constituents !-------------------------------------------------------------------------------------------------- +#ifdef NEWSTATE +subroutine homogenization_RGC_partitionDeformation(F,avgF,ip,el) +#else subroutine homogenization_RGC_partitionDeformation(F,avgF,state,ip,el) +#endif + use prec, only: & p_vec use debug, only: & @@ -285,7 +317,11 @@ subroutine homogenization_RGC_partitionDeformation(F,avgF,state,ip,el) mesh_element use material, only: & homogenization_maxNgrains, & - homogenization_Ngrains,& + homogenization_Ngrains,& +#ifdef NEWSTATE + homogState, & + mappingHomogenization, & +#endif homogenization_typeInstance use FEsolving, only: & theInc,& @@ -294,7 +330,9 @@ subroutine homogenization_RGC_partitionDeformation(F,avgF,state,ip,el) implicit none real(pReal), dimension (3,3,homogenization_maxNgrains), intent(out) :: F !< partioned F per grain real(pReal), dimension (3,3), intent(in) :: avgF !< averaged F +#ifndef NEWSTATE type(p_vec), intent(in) :: state +#endif integer(pInt), intent(in) :: & ip, & !< integration point number el !< element number @@ -326,7 +364,11 @@ subroutine homogenization_RGC_partitionDeformation(F,avgF,state,ip,el) iGrain3 = homogenization_RGC_grain1to3(iGrain,homID) do iFace = 1_pInt,nFace intFace = homogenization_RGC_getInterface(iFace,iGrain3) ! identifying 6 interfaces of each grain +#ifdef NEWSTATE + aVect = homogenization_RGC_relaxationVector(intFace,homID, ip, el) ! get the relaxation vectors for each interface from global relaxation vector array +#else aVect = homogenization_RGC_relaxationVector(intFace,state,homID) ! get the relaxation vectors for each interface from global relaxation vector array +#endif nVect = homogenization_RGC_interfaceNormal(intFace,ip,el) ! get the normal of each interface forall (i=1_pInt:3_pInt,j=1_pInt:3_pInt) & F(i,j,iGrain) = F(i,j,iGrain) + aVect(i)*nVect(j) ! calculating deformation relaxations due to interface relaxation @@ -355,7 +397,11 @@ end subroutine homogenization_RGC_partitionDeformation !> @brief update the internal state of the homogenization scheme and tell whether "done" and ! "happy" with result !-------------------------------------------------------------------------------------------------- +#ifdef NEWSTATE +function homogenization_RGC_updateState( P,F,F0,avgF,dt,dPdF,ip,el) +#else function homogenization_RGC_updateState( state, state0,P,F,F0,avgF,dt,dPdF,ip,el) +#endif use prec, only: & p_vec use debug, only: & @@ -371,6 +417,10 @@ function homogenization_RGC_updateState( state, state0,P,F,F0,avgF,dt,dPdF,ip,el use material, only: & homogenization_maxNgrains, & homogenization_typeInstance, & +#ifdef NEWSTATE + homogState, & + mappingHomogenization, & +#endif homogenization_Ngrains use numerics, only: & absTol_RGC, & @@ -384,8 +434,10 @@ function homogenization_RGC_updateState( state, state0,P,F,F0,avgF,dt,dPdF,ip,el refRelaxRate_RGC implicit none +#ifndef NEWSTATE type(p_vec), intent(inout) :: state !< current state type(p_vec), intent(in) :: state0 !< initial state +#endif real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: & P,& !< array of P F,& !< array of F @@ -429,16 +481,28 @@ function homogenization_RGC_updateState( state, state0,P,F,F0,avgF,dt,dPdF,ip,el ! allocate the size of the global relaxation arrays/jacobian matrices depending on the size of the cluster allocate(resid(3_pInt*nIntFaceTot), source=0.0_pReal) allocate(tract(nIntFaceTot,3), source=0.0_pReal) +#ifdef NEWSTATE + allocate(relax(3_pInt*nIntFaceTot)); relax= homogState(mappingHomogenization(2,ip,el))% & + state(1:3_pInt*nIntFaceTot,mappingHomogenization(1,ip,el)) + allocate(drelax(3_pInt*nIntFaceTot)); drelax= homogState(mappingHomogenization(2,ip,el))% & + state(1:3_pInt*nIntFaceTot,mappingHomogenization(1,ip,el)) - & + homogState(mappingHomogenization(2,ip,el))% & + state0(1:3_pInt*nIntFaceTot,mappingHomogenization(1,ip,el)) +#else allocate(relax(3_pInt*nIntFaceTot)); relax=state%p(1:3_pInt*nIntFaceTot) allocate(drelax(3_pInt*nIntFaceTot)); drelax=state%p(1:3_pInt*nIntFaceTot)-state0%p(1:3_pInt*nIntFaceTot) - +#endif !-------------------------------------------------------------------------------------------------- ! debugging the obtained state if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt) then !$OMP CRITICAL (write2out) write(6,'(1x,a30)')'Obtained state: ' do i = 1_pInt,3_pInt*nIntFaceTot +#ifdef NEWSTATE + write(6,'(1x,2(e15.8,1x))')homogState(mappingHomogenization(2,ip,el))%state(i,mappingHomogenization(1,ip,el)) +#else write(6,'(1x,2(e15.8,1x))')state%p(i) +#endif enddo write(6,*)' ' !$OMP END CRITICAL (write2out) @@ -552,8 +616,13 @@ function homogenization_RGC_updateState( state, state0,P,F,F0,avgF,dt,dPdF,ip,el !-------------------------------------------------------------------------------------------------- ! compute/update the state for postResult, i.e., all energy densities computed by time-integration +#ifdef NEWSTATE + constitutiveWork = homogState(mappingHomogenization(2,ip,el))%state(3*nIntFaceTot+1,mappingHomogenization(1,ip,el)) + penaltyEnergy = homogState(mappingHomogenization(2,ip,el))%state(3*nIntFaceTot+5,mappingHomogenization(1,ip,el)) +#else constitutiveWork = state%p(3*nIntFaceTot+1) penaltyEnergy = state%p(3*nIntFaceTot+5) +#endif do iGrain = 1_pInt,homogenization_Ngrains(mesh_element(3,el)) ! time-integration loop for the calculating the work and energy do i = 1_pInt,3_pInt do j = 1_pInt,3_pInt @@ -562,6 +631,26 @@ function homogenization_RGC_updateState( state, state0,P,F,F0,avgF,dt,dPdF,ip,el enddo enddo enddo +#ifdef NEWSTATE + homogState(mappingHomogenization(2,ip,el))% & + state(3*nIntFaceTot+1,mappingHomogenization(1,ip,el)) = constitutiveWork ! the bulk mechanical/constitutive work + homogState(mappingHomogenization(2,ip,el))% & + state(3*nIntFaceTot+2,mappingHomogenization(1,ip,el)) = sum(NN(1,:))/real(nGrain,pReal) ! the overall mismatch of all interface normal to e1-direction + homogState(mappingHomogenization(2,ip,el))% & + state(3*nIntFaceTot+3,mappingHomogenization(1,ip,el)) = sum(NN(2,:))/real(nGrain,pReal) ! the overall mismatch of all interface normal to e2-direction + homogState(mappingHomogenization(2,ip,el))% & + state(3*nIntFaceTot+4,mappingHomogenization(1,ip,el)) = sum(NN(3,:))/real(nGrain,pReal) ! the overall mismatch of all interface normal to e3-direction + homogState(mappingHomogenization(2,ip,el))% & + state(3*nIntFaceTot+5,mappingHomogenization(1,ip,el)) = penaltyEnergy ! the overall penalty energy + homogState(mappingHomogenization(2,ip,el))% & + state(3*nIntFaceTot+6,mappingHomogenization(1,ip,el)) = volDiscrep ! the overall volume discrepancy + homogState(mappingHomogenization(2,ip,el))% & + state(3*nIntFaceTot+7,mappingHomogenization(1,ip,el)) = & + sum(abs(drelax))/dt/real(3_pInt*nIntFaceTot,pReal) ! the average rate of relaxation vectors + homogState(mappingHomogenization(2,ip,el))% & + state(3*nIntFaceTot+8,mappingHomogenization(1,ip,el)) = maxval(abs(drelax))/dt ! the maximum rate of relaxation vectors + +#else state%p(3*nIntFaceTot+1) = constitutiveWork ! the bulk mechanical/constitutive work state%p(3*nIntFaceTot+2) = sum(NN(1,:))/real(nGrain,pReal) ! the overall mismatch of all interface normal to e1-direction state%p(3*nIntFaceTot+3) = sum(NN(2,:))/real(nGrain,pReal) ! the overall mismatch of all interface normal to e2-direction @@ -571,6 +660,7 @@ function homogenization_RGC_updateState( state, state0,P,F,F0,avgF,dt,dPdF,ip,el state%p(3*nIntFaceTot+7) = sum(abs(drelax))/dt/real(3_pInt*nIntFaceTot,pReal) ! the average rate of relaxation vectors state%p(3*nIntFaceTot+8) = maxval(abs(drelax))/dt ! the maximum rate of relaxation vectors +#endif if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt & .and. debug_e == el .and. debug_i == ip) then !$OMP CRITICAL (write2out) @@ -639,8 +729,8 @@ function homogenization_RGC_updateState( state, state0,P,F,F0,avgF,dt,dPdF,ip,el do i=1_pInt,3_pInt; do j=1_pInt,3_pInt; do k=1_pInt,3_pInt; do l=1_pInt,3_pInt smatrix(3*(iNum-1)+i,3*(iMun-1)+j) = smatrix(3*(iNum-1)+i,3*(iMun-1)+j) + dPdF(i,k,j,l,iGrN)*normN(k)*mornN(l) enddo;enddo;enddo;enddo - ! projecting the material tangent dPdF into the interface - ! to obtain the Jacobian matrix contribution of dPdF +! projecting the material tangent dPdF into the interface +! to obtain the Jacobian matrix contribution of dPdF endif enddo @@ -685,8 +775,13 @@ function homogenization_RGC_updateState( state, state0,P,F,F0,avgF,dt,dPdF,ip,el do ipert = 1_pInt,3_pInt*nIntFaceTot p_relax = relax p_relax(ipert) = relax(ipert) + pPert_RGC ! perturb the relaxation vector +#ifdef NEWSTATE + homogState(mappingHomogenization(2,ip,el))%state(1:3*nIntFaceTot,mappingHomogenization(1,ip,el)) = p_relax + call homogenization_RGC_grainDeformation(pF,avgF,ip,el) ! compute the grains deformation from perturbed state +#else state%p(1:3*nIntFaceTot) = p_relax call homogenization_RGC_grainDeformation(pF,avgF,state,ip,el) ! compute the grains deformation from perturbed state +#endif call homogenization_RGC_stressPenalty(pR,pNN,avgF,pF,ip,el,homID) ! compute stress penalty due to interface mismatch from perturbed state call homogenization_RGC_volumePenalty(pD,volDiscrep,pF,avgF,ip,el) ! compute stress penalty due to volume discrepancy from perturbed state @@ -801,7 +896,11 @@ function homogenization_RGC_updateState( state, state0,P,F,F0,avgF,dt,dPdF,ip,el enddo enddo relax = relax + drelax ! Updateing the state variable for the next iteration +#ifdef NEWSTATE + homogState(mappingHomogenization(2,ip,el))%state(1:3*nIntFaceTot,mappingHomogenization(1,ip,el)) = relax +#else state%p(1:3*nIntFaceTot) = relax +#endif if (any(abs(drelax) > maxdRelax_RGC)) then ! Forcing cutback when the incremental change of relaxation vector becomes too large homogenization_RGC_updateState = [.true.,.false.] !$OMP CRITICAL (write2out) @@ -817,7 +916,11 @@ function homogenization_RGC_updateState( state, state0,P,F,F0,avgF,dt,dPdF,ip,el !$OMP CRITICAL (write2out) write(6,'(1x,a30)')'Returned state: ' do i = 1_pInt,3_pInt*nIntFaceTot +#ifdef NEWSTATE + write(6,'(1x,2(e15.8,1x))')homogState(mappingHomogenization(2,ip,el))%state(i,mappingHomogenization(1,ip,el)) +#else write(6,'(1x,2(e15.8,1x))')state%p(i) +#endif enddo write(6,*)' ' flush(6) @@ -838,7 +941,14 @@ subroutine homogenization_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF, debug_homogenization,& debug_levelExtensive use mesh, only: mesh_element - use material, only: homogenization_maxNgrains,homogenization_Ngrains,homogenization_typeInstance + use material, only: & + homogenization_maxNgrains, & +#ifdef NEWSTATE + homogState, & + mappingHomogenization, & +#endif + homogenization_Ngrains, & + homogenization_typeInstance use math, only: math_Plain3333to99 implicit none @@ -881,7 +991,11 @@ end subroutine homogenization_RGC_averageStressAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief return array of homogenization results for post file inclusion !-------------------------------------------------------------------------------------------------- +#ifdef NEWSTATE +pure function homogenization_RGC_postResults(ip,el,avgP,avgF) +#else pure function homogenization_RGC_postResults(state,ip,el,avgP,avgF) +#endif use prec, only: & p_vec use mesh, only: & @@ -889,6 +1003,10 @@ pure function homogenization_RGC_postResults(state,ip,el,avgP,avgF) mesh_ipCoordinates use material, only: & homogenization_typeInstance,& +#ifdef NEWSTATE + homogState, & + mappingHomogenization, & +#endif homogenization_Noutput use crystallite, only: & crystallite_temperature @@ -900,8 +1018,10 @@ pure function homogenization_RGC_postResults(state,ip,el,avgP,avgF) real(pReal), dimension(3,3), intent(in) :: & avgP, & !< average stress at material point avgF !< average deformation gradient at material point +#ifndef NEWSTATE type(p_vec), intent(in) :: & state ! my State +#endif integer(pInt) homID,o,c,nIntFaceTot real(pReal), dimension(homogenization_RGC_sizePostResults(homogenization_typeInstance(mesh_element(3,el)))) :: & homogenization_RGC_postResults @@ -928,25 +1048,63 @@ pure function homogenization_RGC_postResults(state,ip,el,avgP,avgF) homogenization_RGC_postResults(c+1_pInt:c+3_pInt) = mesh_ipCoordinates(1:3,ip,el) ! current ip coordinates c = c + 3_pInt case (constitutivework_ID) +#ifdef NEWSTATE + homogenization_RGC_postResults(c+1) = homogState(mappingHomogenization(2,ip,el))% & + state(3*nIntFaceTot+1,mappingHomogenization(1,ip,el)) +#else homogenization_RGC_postResults(c+1) = state%p(3*nIntFaceTot+1) +#endif c = c + 1_pInt case (magnitudemismatch_ID) +#ifdef NEWSTATE + homogenization_RGC_postResults(c+1) = homogState(mappingHomogenization(2,ip,el))% & + state(3*nIntFaceTot+2,mappingHomogenization(1,ip,el)) + homogenization_RGC_postResults(c+2) = homogState(mappingHomogenization(2,ip,el))% & + state(3*nIntFaceTot+3,mappingHomogenization(1,ip,el)) + homogenization_RGC_postResults(c+3) = homogState(mappingHomogenization(2,ip,el))% & + state(3*nIntFaceTot+4,mappingHomogenization(1,ip,el)) +#else homogenization_RGC_postResults(c+1) = state%p(3*nIntFaceTot+2) homogenization_RGC_postResults(c+2) = state%p(3*nIntFaceTot+3) homogenization_RGC_postResults(c+3) = state%p(3*nIntFaceTot+4) +#endif c = c + 3_pInt case (penaltyenergy_ID) +#ifdef NEWSTATE + homogenization_RGC_postResults(c+1) = homogState(mappingHomogenization(2,ip,el))% & + state(3*nIntFaceTot+5,mappingHomogenization(1,ip,el)) + c = c + 1_pInt +#else homogenization_RGC_postResults(c+1) = state%p(3*nIntFaceTot+5) c = c + 1_pInt +#endif case (volumediscrepancy_ID) +#ifdef NEWSTATE + homogenization_RGC_postResults(c+1) = homogState(mappingHomogenization(2,ip,el))% & + state(3*nIntFaceTot+6,mappingHomogenization(1,ip,el)) + c = c + 1_pInt +#else homogenization_RGC_postResults(c+1) = state%p(3*nIntFaceTot+6) c = c + 1_pInt +#endif case (averagerelaxrate_ID) +#ifdef NEWSTATE + homogenization_RGC_postResults(c+1) = homogState(mappingHomogenization(2,ip,el))% & + state(3*nIntFaceTot+7,mappingHomogenization(1,ip,el)) + c = c + 1_pInt +#else homogenization_RGC_postResults(c+1) = state%p(3*nIntFaceTot+7) c = c + 1_pInt +#endif case (maximumrelaxrate_ID) +#ifdef NEWSTATE + homogenization_RGC_postResults(c+1) = homogState(mappingHomogenization(2,ip,el))% & + state(3*nIntFaceTot+8,mappingHomogenization(1,ip,el)) + c = c + 1_pInt +#else homogenization_RGC_postResults(c+1) = state%p(3*nIntFaceTot+8) c = c + 1_pInt +#endif end select enddo @@ -972,6 +1130,10 @@ subroutine homogenization_RGC_stressPenalty(rPen,nMis,avgF,fDef,ip,el,homID) math_invert33 use material, only: & homogenization_maxNgrains,& +#ifdef NEWSTATE + homogState, & + mappingHomogenization, & +#endif homogenization_Ngrains use numerics, only: & xSmoo_RGC @@ -1111,6 +1273,10 @@ subroutine homogenization_RGC_volumePenalty(vPen,vDiscrep,fDef,fAvg,ip,el) math_inv33 use material, only: & homogenization_maxNgrains,& +#ifdef NEWSTATE + homogState, & + mappingHomogenization, & +#endif homogenization_Ngrains use numerics, only: & maxVolDiscr_RGC,& @@ -1170,6 +1336,11 @@ function homogenization_RGC_surfaceCorrection(avgF,ip,el) use math, only: & math_invert33, & math_mul33x33 +#ifdef NEWSTATE + use material, only: & + homogState, & + mappingHomogenization +#endif implicit none real(pReal), dimension(3) :: homogenization_RGC_surfaceCorrection @@ -1239,14 +1410,28 @@ end function homogenization_RGC_equivalentModuli !-------------------------------------------------------------------------------------------------- !> @brief collect relaxation vectors of an interface !-------------------------------------------------------------------------------------------------- +#ifdef NEWSTATE +function homogenization_RGC_relaxationVector(intFace,homID, ip, el) +#else function homogenization_RGC_relaxationVector(intFace,state,homID) +#endif use prec, only: & p_vec +#ifdef NEWSTATE + use material, only: & + homogState, & + mappingHomogenization +#endif implicit none +#ifdef NEWSTATE + integer(pInt), intent(in) :: ip, el +#endif real(pReal), dimension (3) :: homogenization_RGC_relaxationVector integer(pInt), dimension (4), intent(in) :: intFace !< set of interface ID in 4D array (normal and position) +#ifndef NEWSTATE type(p_vec), intent(in) :: state !< set of global relaxation vectors +#endif integer(pInt), dimension (3) :: nGDim integer(pInt) :: & iNum, & @@ -1257,7 +1442,12 @@ function homogenization_RGC_relaxationVector(intFace,state,homID) homogenization_RGC_relaxationVector = 0.0_pReal nGDim = homogenization_RGC_Ngrains(1:3,homID) iNum = homogenization_RGC_interface4to1(intFace,homID) ! identify the position of the interface in global state array +#ifdef NEWSTATE + if (iNum > 0_pInt) homogenization_RGC_relaxationVector = homogState(mappingHomogenization(2,ip,el))% & + state((3*iNum-2):(3*iNum),mappingHomogenization(1,ip,el)) ! get the corresponding entries +#else if (iNum > 0_pInt) homogenization_RGC_relaxationVector = state%p((3*iNum-2):(3*iNum)) ! get the corresponding entries +#endif end function homogenization_RGC_relaxationVector @@ -1471,7 +1661,11 @@ end function homogenization_RGC_interface1to4 !> @brief calculating the grain deformation gradient (the same with ! homogenization_RGC_partionDeformation, but used only for perturbation scheme) !-------------------------------------------------------------------------------------------------- +#ifdef NEWSTATE +subroutine homogenization_RGC_grainDeformation(F, avgF, ip, el) +#else subroutine homogenization_RGC_grainDeformation(F, avgF, state, ip, el) +#endif use prec, only: & p_vec use mesh, only: & @@ -1479,12 +1673,18 @@ subroutine homogenization_RGC_grainDeformation(F, avgF, state, ip, el) use material, only: & homogenization_maxNgrains,& homogenization_Ngrains, & +#ifdef NEWSTATE + homogState, & + mappingHomogenization, & +#endif homogenization_typeInstance implicit none real(pReal), dimension (3,3,homogenization_maxNgrains), intent(out) :: F !< partioned F per grain real(pReal), dimension (3,3), intent(in) :: avgF !< +#ifndef NEWSTATE type(p_vec), intent(in) :: state +#endif integer(pInt), intent(in) :: & el, & !< element number ip !< integration point number @@ -1502,7 +1702,11 @@ subroutine homogenization_RGC_grainDeformation(F, avgF, state, ip, el) iGrain3 = homogenization_RGC_grain1to3(iGrain,homID) do iFace = 1_pInt,nFace intFace = homogenization_RGC_getInterface(iFace,iGrain3) +#ifdef NEWSTATE + aVect = homogenization_RGC_relaxationVector(intFace,homID, ip, el) +#else aVect = homogenization_RGC_relaxationVector(intFace,state,homID) +#endif nVect = homogenization_RGC_interfaceNormal(intFace,ip,el) forall (i=1_pInt:3_pInt,j=1_pInt:3_pInt) & F(i,j,iGrain) = F(i,j,iGrain) + aVect(i)*nVect(j) ! effective relaxations diff --git a/code/homogenization_isostrain.f90 b/code/homogenization_isostrain.f90 index 9842c9084..33fa9daac 100644 --- a/code/homogenization_isostrain.f90 +++ b/code/homogenization_isostrain.f90 @@ -51,6 +51,9 @@ contains !-------------------------------------------------------------------------------------------------- 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 prec, only: & + pReal, & + pInt use IO use material @@ -61,7 +64,15 @@ subroutine homogenization_isostrain_init(fileUnit) integer(pInt) :: & section = 0_pInt, i, j, output, mySize integer :: & - maxNinstance, k ! no pInt (stores a system dependen value from 'count' + maxNinstance, & +#ifdef NEWSTATE + homog, & + NofMyHomog, & + instance, & + sizeHState, & +#endif + k +! no pInt (stores a system dependen value from 'count' character(len=65536) :: & tag = '', & line = '' @@ -144,6 +155,24 @@ subroutine homogenization_isostrain_init(fileUnit) endif enddo +#ifdef NEWSTATE + initializeInstances: do homog = 1_pInt, material_Nhomogenization + + myhomog: if (homogenization_type(homog) == HOMOGENIZATION_ISOSTRAIN_ID) then + NofMyHomog = count(material_homog == homog) +! instance = phase_plasticityInstance(phase) + +! allocate homogenization state arrays + sizeHState = 0_pInt + homogState(homog)%sizeState = sizeHState + homogState(homog)%sizePostResults = homogenization_isostrain_sizePostResults(homog) + allocate(homogState(homog)%state0 ( sizeHState,NofMyHomog), source=0.0_pReal) + allocate(homogState(homog)%subState0 ( sizeHState,NofMyHomog), source=0.0_pReal) + allocate(homogState(homog)%state ( sizeHState,NofMyHomog), source=0.0_pReal) + + endif myhomog + enddo initializeInstances +#endif do i = 1,maxNinstance do j = 1_pInt,maxval(homogenization_Noutput) diff --git a/code/homogenization_none.f90 b/code/homogenization_none.f90 index bfb4a7d56..c58337719 100644 --- a/code/homogenization_none.f90 +++ b/code/homogenization_none.f90 @@ -21,18 +21,45 @@ contains !-------------------------------------------------------------------------------------------------- subroutine homogenization_none_init() use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) + use prec, only: & + pReal, & + pInt use IO, only: & IO_timeStamp - use material, only: & - HOMOGENIZATION_NONE_label - - implicit none + use material + implicit none +#ifdef NEWSTATE + integer :: & + homog, & + NofMyHomog, & + instance, & + sizeHState +#endif write(6,'(/,a)') ' <<<+- homogenization_'//HOMOGENIZATION_NONE_label//' init -+>>>' write(6,'(a)') ' $Id$' write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" +#ifdef NEWSTATE + initializeInstances: do homog = 1_pInt, material_Nhomogenization + + myhomog: if (homogenization_type(homog) == HOMOGENIZATION_none_ID) then + NofMyHomog = count(material_homog == homog) +! instance = phase_plasticityInstance(phase) + +! allocate homogenization state arrays + sizeHState = 0_pInt + homogState(homog)%sizeState = sizeHState + homogState(homog)%sizePostResults = 0_pInt + allocate(homogState(homog)%state0 ( sizeHState,NofMyHomog), source=0.0_pReal) + allocate(homogState(homog)%subState0 ( sizeHState,NofMyHomog), source=0.0_pReal) + allocate(homogState(homog)%state ( sizeHState,NofMyHomog), source=0.0_pReal) + + endif myhomog + enddo initializeInstances +#endif + end subroutine homogenization_none_init end module homogenization_none diff --git a/code/material.f90 b/code/material.f90 index 99b77d983..555d4b695 100644 --- a/code/material.f90 +++ b/code/material.f90 @@ -13,6 +13,9 @@ module material pReal, & pInt, & tState, & +#ifdef NEWSTATE + hState, & +#endif p_intvec implicit none @@ -116,12 +119,19 @@ module material integer(pInt), dimension(:,:,:), allocatable, public :: & material_phase !< phase (index) of each grain,IP,element - +#ifdef NEWSTATE + integer(pInt), dimension(:,:), allocatable, public :: & + material_homog !< homogenization (index) of each IP,element +#endif type(tState), allocatable, dimension(:), public :: & plasticState, & elasticState, & damageState, & thermalState +#ifdef NEWSTATE + type(hState), allocatable, dimension(:), public :: & + homogState +#endif integer(pInt), dimension(:,:,:), allocatable, public, protected :: & @@ -177,8 +187,14 @@ module material integer(pInt), dimension(:,:,:,:), allocatable, public, protected :: mappingConstitutive integer(pInt), dimension(:,:,:), allocatable, public, protected :: mappingCrystallite +#ifdef NEWSTATE + integer(pInt), dimension(:,:,:), allocatable, public, protected :: mappingHomogenization +#endif integer(pInt), dimension(:), allocatable :: ConstitutivePosition integer(pInt), dimension(:), allocatable :: CrystallitePosition +#ifdef NEWSTATE + integer(pInt), dimension(:), allocatable :: HomogenizationPosition +#endif public :: & @@ -242,7 +258,7 @@ subroutine material_init implicit none integer(pInt), parameter :: FILEUNIT = 200_pInt - integer(pInt) :: m,c,h, myDebug + integer(pInt) :: m,c,h, myDebug, myHomogInstance integer(pInt) :: & g, & !< grain number i, & !< integration point number @@ -273,7 +289,9 @@ subroutine material_init allocate(elasticState(material_Nphase)) allocate(damageState (material_Nphase)) allocate(thermalState(material_Nphase)) - +#ifdef NEWSTATE + allocate(homogState(material_Nhomogenization)) +#endif do m = 1_pInt,material_Nmicrostructure if(microstructure_crystallite(m) < 1_pInt .or. & microstructure_crystallite(m) > material_Ncrystallite) & @@ -313,11 +331,22 @@ subroutine material_init call material_populateGrains allocate(mappingConstitutive(2,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems),source=0_pInt) +#ifdef NEWSTATE + allocate(mappingHomogenization(2,mesh_maxNips,mesh_NcpElems),source=0_pInt) +#endif allocate(mappingCrystallite (2,homogenization_maxNgrains,mesh_NcpElems),source=0_pInt) allocate(ConstitutivePosition(material_Nphase),source=0_pInt) +#ifdef NEWSTATE + allocate(HomogenizationPosition(material_Nhomogenization),source=0_pInt) +#endif allocate(CrystallitePosition(material_Nphase),source=0_pInt) ElemLoop:do e = 1_pInt,mesh_NcpElems ! loop over elements + myHomogInstance = homogenization_typeInstance(mesh_element(3,e)) IPloop:do i = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,e))) ! loop over IPs +#ifdef NEWSTATE + HomogenizationPosition(myHomogInstance) = HomogenizationPosition(myHomogInstance)+1_pInt + mappingHomogenization(1:2,i,e) = [HomogenizationPosition(myHomogInstance),myHomogInstance] +#endif GrainLoop:do g = 1_pInt,homogenization_Ngrains(mesh_element(3,e)) ! loop over grains phase = material_phase(g,i,e) ConstitutivePosition(phase) = ConstitutivePosition(phase)+1_pInt ! not distinguishing between instances of same phase @@ -941,9 +970,9 @@ subroutine material_populateGrains real(pReal), dimension (3) :: orientation real(pReal), dimension (3,3) :: symOrientation integer(pInt), dimension (:), allocatable :: phaseOfGrain, textureOfGrain - integer(pInt) :: t,e,i,g,j,m,c,r,homog,micro,sgn,hme, myDebug, & + integer(pInt) :: t,e,i,ii,g,j,m,c,r,homog,micro,sgn,hme, myDebug, & phaseID,textureID,dGrains,myNgrains,myNorientations,myNconstituents, & - grain,constituentGrain,ipGrain,symExtension, ip + grain,constituentGrain,ipGrain,symExtension, ip, HomogInstType real(pReal) :: extreme,rnd integer(pInt), dimension (:,:), allocatable :: Nelems ! counts number of elements in homog, micro array type(p_intvec), dimension (:,:), allocatable :: elemsOfHomogMicro ! lists element number in homog, micro array @@ -952,12 +981,25 @@ subroutine material_populateGrains allocate(material_volume(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), source=0.0_pReal) allocate(material_phase(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), source=0_pInt) +#ifdef NEWSTATE + allocate(material_homog(mesh_maxNips,mesh_NcpElems), source=0_pInt) +#endif allocate(material_texture(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), source=0_pInt) allocate(material_EulerAngles(3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems),source=0.0_pReal) allocate(Ngrains(material_Nhomogenization,material_Nmicrostructure), source=0_pInt) allocate(Nelems(material_Nhomogenization,material_Nmicrostructure), source=0_pInt) +#ifdef NEWSTATE +! populating homogenization schemes in each +!-------------------------------------------------------------------------------------------------- + do e = 1_pInt, mesh_NcpElems +! do i = 1_pInt:FE_Nips(FE_geomtype(mesh_element(2,e))) + material_homog(1_pInt:FE_Nips(FE_geomtype(mesh_element(2,e))),e)=homogenization_typeInstance(mesh_element(3,e)) +! enddo + enddo + +#endif !-------------------------------------------------------------------------------------------------- ! precounting of elements for each homog/micro pair do e = 1_pInt, mesh_NcpElems diff --git a/code/prec.f90 b/code/prec.f90 index b602114f2..4ce28177a 100644 --- a/code/prec.f90 +++ b/code/prec.f90 @@ -77,6 +77,15 @@ module prec RK4dotState real(pReal), allocatable, dimension(:,:,:) :: RKCK45dotState end type +#ifdef NEWSTATE + type, public :: hState + integer(pInt) :: sizeState = 0_pInt , & + sizePostResults = 0_pInt + real(pReal), allocatable, dimension(:,:) :: state, & ! material points, state size + state0, & + subState0 + end type +#endif public :: & prec_init