From 59d22d47b2361aee2637065701e19258fe118208 Mon Sep 17 00:00:00 2001 From: Denny Tjahjanto Date: Tue, 24 Nov 2009 15:00:25 +0000 Subject: [PATCH] introduction of a possibility to have homogeneous element (all ips in one element are identical, sort of reduced integration) and bugs-fixing in crystallite.f90, homogenization_RGC.f90, numerics.f90. --- code/crystallite.f90 | 1 - code/homogenization_RGC.f90 | 4 +-- code/material.config | 2 ++ code/material.f90 | 51 +++++++++++++++++++++++++++---------- code/numerics.f90 | 8 +++--- 5 files changed, 45 insertions(+), 21 deletions(-) diff --git a/code/crystallite.f90 b/code/crystallite.f90 index 9cb706d51..5cc2869bd 100644 --- a/code/crystallite.f90 +++ b/code/crystallite.f90 @@ -300,7 +300,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco) real(pReal), dimension(3,3,3,3) :: dPdF_pos,dPdF_neg ! ------ initialize to starting condition ------ - centralDifference = .true. !$OMP CRITICAL (write2out) ! write (6,*) diff --git a/code/homogenization_RGC.f90 b/code/homogenization_RGC.f90 index 2cdfe2652..ff5a12c51 100644 --- a/code/homogenization_RGC.f90 +++ b/code/homogenization_RGC.f90 @@ -303,7 +303,7 @@ function homogenization_RGC_updateState(& allocate(resid(3*nIntFaceTot)); resid = 0.0_pReal allocate(tract(nIntFaceTot,3)); tract = 0.0_pReal allocate(relax(3*nIntFaceTot)); relax = state%p(1:3*nIntFaceTot) - allocate(drelax(3*nIntFaceTot)); & + allocate(drelax(3*nIntFaceTot)) drelax = state%p(1:3*nIntFaceTot) - state0%p(1:3*nIntFaceTot) !* Debugging the obtained state if (RGCdebug) then @@ -525,7 +525,7 @@ function homogenization_RGC_updateState(& !* Construct the Jacobian matrix of the numerical viscosity tangent allocate(rmatrix(3*nIntFaceTot,3*nIntFaceTot)); rmatrix = 0.0_pReal - forall (i=1,3*nIntFaceTot) & + forall (i=1:3*nIntFaceTot) & rmatrix(i,i) = viscModus_RGC*ratePower_RGC/dt*(abs(drelax(i)/dt))**(ratePower_RGC - 1.0_pReal) !* Debugging the global Jacobian matrix of numerical viscosity tangent if (RGCdebugJacobi) then diff --git a/code/material.config b/code/material.config index cdb13052a..a40256d7b 100644 --- a/code/material.config +++ b/code/material.config @@ -27,6 +27,7 @@ Ngrains 2 ##################### [Aluminum_Poly] +/elementhomogeneous/ # put this flag to set ips identical in one element (something like reduced integration) (constituent) phase 3 texture 1 fraction 1.0 [Aluminum_001] @@ -36,6 +37,7 @@ Ngrains 2 (constituent) phase 1 texture 1 fraction 1.0 [DP_Steel] +/elementhomogeneous/ (constituent) phase 1 texture 1 fraction 0.82 (constituent) phase 2 texture 1 fraction 0.18 diff --git a/code/material.f90 b/code/material.f90 index 862aef56e..ede6d8a72 100644 --- a/code/material.f90 +++ b/code/material.f90 @@ -49,6 +49,7 @@ integer(pInt), dimension(:), allocatable :: homogenization_Ngrains, & texture_Nfiber ! number of Fiber components per texture logical, dimension(:), allocatable :: homogenization_active, & ! microstructure_active, & ! + microstructure_elemhomo, & ! flag to indicate homogeneous microstructure distribution over element's IPs phase_localConstitution ! flags phases with local constitutive law integer(pInt), dimension(:,:), allocatable :: microstructure_phase, & ! phase IDs of each microstructure microstructure_texture ! texture IDs of each microstructure @@ -102,9 +103,9 @@ subroutine material_init() write (6,'(x,a32,x,a16,x,i4)') homogenization_name(i),homogenization_type(i),homogenization_Ngrains(i) enddo write (6,*) - write (6,'(a32,x,a12)') 'microstructure ','constituents' + write (6,'(a32,x,a12,x,a13)') 'microstructure ','constituents','homogeneous' do i = 1,material_Nmicrostructure - write (6,'(a32,x,i4)') microstructure_name(i),microstructure_Nconstituents(i) + write (6,'(a32,4x,i4,8x,l)') microstructure_name(i),microstructure_Nconstituents(i),microstructure_elemhomo(i) if (microstructure_Nconstituents(i) > 0_pInt) then do j = 1,microstructure_Nconstituents(i) write (6,'(a1,x,a32,x,a32,x,f6.4)') '>',phase_name(microstructure_phase(j,i)),& @@ -207,14 +208,17 @@ subroutine material_parseMicrostructure(file,myPart) Nsections = IO_countSections(file,myPart) material_Nmicrostructure = Nsections - allocate(microstructure_name(Nsections)); microstructure_name = '' + allocate(microstructure_name(Nsections)); microstructure_name = '' allocate(microstructure_Nconstituents(Nsections)) allocate(microstructure_active(Nsections)) + allocate(microstructure_elemhomo(Nsections)) forall (i = 1:Nsections) microstructure_active(i) = any(mesh_element(4,:) == i) ! current microstructure used in model? microstructure_Nconstituents = IO_countTagInPart(file,myPart,'(constituent)',Nsections) microstructure_maxNconstituents = maxval(microstructure_Nconstituents) + microstructure_elemhomo = IO_spotTagInPart(file,myPart,'/elementhomogeneous/',Nsections) + allocate(microstructure_phase (microstructure_maxNconstituents,Nsections)); microstructure_phase = 0_pInt allocate(microstructure_texture (microstructure_maxNconstituents,Nsections)); microstructure_texture = 0_pInt allocate(microstructure_fraction(microstructure_maxNconstituents,Nsections)); microstructure_fraction = 0.0_pReal @@ -475,7 +479,12 @@ subroutine material_populateGrains() call IO_error(130,e,0,0) if (micro < 1 .or. micro > material_Nmicrostructure) & ! out of bounds call IO_error(140,e,0,0) - Ngrains(homog,micro) = Ngrains(homog,micro) + homogenization_Ngrains(homog) * FE_Nips(mesh_element(2,e)) + if (microstructure_elemhomo(micro)) then + dGrains = homogenization_Ngrains(homog) + else + dGrains = homogenization_Ngrains(homog) * FE_Nips(mesh_element(2,e)) + endif + Ngrains(homog,micro) = Ngrains(homog,micro) + dGrains enddo allocate(volumeOfGrain(maxval(Ngrains))) ! reserve memory for maximum case @@ -499,10 +508,15 @@ subroutine material_populateGrains() grain = 0_pInt ! microstructure grain index do e = 1,mesh_NcpElems ! check each element if (mesh_element(3,e) == homog .and. mesh_element(4,e) == micro) then ! my combination of homog and micro - forall (i = 1:FE_Nips(mesh_element(2,e))) & ! loop over IPs - volumeOfGrain(grain+(i-1)*dGrains+1:grain+i*dGrains) = & - mesh_ipVolume(i,e)/dGrains ! assign IPvolume/Ngrains to all grains of IP - grain = grain + FE_Nips(mesh_element(2,e)) * dGrains ! wind forward by Nips*NgrainsPerIP + if (microstructure_elemhomo(micro)) then ! homogeneous distribution of grains over each element's IPs + volumeOfGrain(grain+1:grain+dGrains) = sum(mesh_ipVolume(1:FE_Nips(mesh_element(2,e)),e))/dGrains + grain = grain + dGrains ! wind forward by NgrainsPerIP + else + forall (i = 1:FE_Nips(mesh_element(2,e))) & ! loop over IPs + volumeOfGrain(grain+(i-1)*dGrains+1:grain+i*dGrains) = & + mesh_ipVolume(i,e)/dGrains ! assign IPvolume/Ngrains to all grains of IP + grain = grain + FE_Nips(mesh_element(2,e)) * dGrains ! wind forward by Nips*NgrainsPerIP + endif endif enddo ! ---------------------------------------------------------------------------- divide myNgrains as best over constituents @@ -604,15 +618,24 @@ subroutine material_populateGrains() grain = 0_pInt ! microstructure grain index do e = 1,mesh_NcpElems ! check each element if (mesh_element(3,e) == homog .and. mesh_element(4,e) == micro) then ! my combination of homog and micro - forall (i = 1:FE_Nips(mesh_element(2,e)), g = 1:dGrains) ! loop over IPs and grains - material_volume(g,i,e) = volumeOfGrain(grain+(i-1)*dGrains+g) - material_phase(g,i,e) = phaseOfGrain(grain+(i-1)*dGrains+g) - material_EulerAngles(:,g,i,e) = orientationOfGrain(:,grain+(i-1)*dGrains+g) - end forall + if (microstructure_elemhomo(micro)) then ! homogeneous distribution of grains over each element's IPs + forall (i = 1:FE_Nips(mesh_element(2,e)), g = 1:dGrains) ! loop over IPs and grains + material_volume(g,i,e) = volumeOfGrain(grain+g) + material_phase(g,i,e) = phaseOfGrain(grain+g) + material_EulerAngles(:,g,i,e) = orientationOfGrain(:,grain+g) + end forall + grain = grain + dGrains ! wind forward by NgrainsPerIP + else + forall (i = 1:FE_Nips(mesh_element(2,e)), g = 1:dGrains) ! loop over IPs and grains + material_volume(g,i,e) = volumeOfGrain(grain+(i-1)*dGrains+g) + material_phase(g,i,e) = phaseOfGrain(grain+(i-1)*dGrains+g) + material_EulerAngles(:,g,i,e) = orientationOfGrain(:,grain+(i-1)*dGrains+g) + end forall + grain = grain + FE_Nips(mesh_element(2,e)) * dGrains ! wind forward by Nips*NgrainsPerIP + endif ! write (6,*) e ! write (6,*) material_phase(:,:,e) ! write (6,*) material_EulerAngles(:,:,:,e) - grain = grain + FE_Nips(mesh_element(2,e)) * dGrains ! wind forward by Nips*NgrainsPerIP endif enddo diff --git a/code/numerics.f90 b/code/numerics.f90 index 9689817f2..fb3303478 100644 --- a/code/numerics.f90 +++ b/code/numerics.f90 @@ -109,7 +109,7 @@ subroutine numerics_init() pPert_RGC = 1.0e-7 xSmoo_RGC = 1.0e-5 ratePower_RGC = 1.0e+0 ! Newton viscosity (linear model) - viscModul_RGC = 0.0e+0 ! No viscosity is applied + viscModus_RGC = 0.0e+0 ! No viscosity is applied maxdRelax_RGC = 1.0e+0 !* Random seeding parameters: added <<>> @@ -186,7 +186,7 @@ subroutine numerics_init() case ('viscosityrate_rgc') ratePower_RGC = IO_floatValue(line,positions,2) case ('viscositymodulus_rgc') - viscModul_RGC = IO_floatValue(line,positions,2) + viscModus_RGC = IO_floatValue(line,positions,2) case ('maxrelaxation_rgc') maxdRelax_RGC = IO_floatValue(line,positions,2) @@ -238,7 +238,7 @@ subroutine numerics_init() write(6,'(a24,x,e8.1)') 'perturbPenalty_RGC: ',pPert_RGC write(6,'(a24,x,e8.1)') 'relevantMismatch_RGC: ',xSmoo_RGC write(6,'(a24,x,e8.1)') 'viscosityrate_RGC: ',ratePower_RGC - write(6,'(a24,x,e8.1)') 'viscositymodulus_RGC: ',viscModul_RGC + write(6,'(a24,x,e8.1)') 'viscositymodulus_RGC: ',viscModus_RGC write(6,'(a24,x,e8.1)') 'maxrelaxation_RGC: ',maxdRelax_RGC write(6,*) @@ -277,7 +277,7 @@ subroutine numerics_init() if (pPert_RGC <= 0.0_pReal) call IO_error(276) !! oops !! if (xSmoo_RGC <= 0.0_pReal) call IO_error(277) if (ratePower_RGC < 0.0_pReal) call IO_error(278) - if (viscModul_RGC < 0.0_pReal) call IO_error(278) + if (viscModus_RGC < 0.0_pReal) call IO_error(278) if (maxdRelax_RGC <= 0.0_pReal) call IO_error(288) if (fixedSeed <= 0_pInt) write(6,'(a)') 'Random is random!'