223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

mpie_spectral and mpie_interface: changed workingdir from pathToLoadFile to pathToGeomFile 

mpie_spectral and numerics: added switch to prevent pre calculation of gamma_hat. slower, but saves memory
3Dvisualize: started to add support for gmsh (not fully working yet)
reconstruct: new version of f2py/Fortran subroutines for output of results from spectral method
This commit is contained in:
Martin Diehl 2011-02-07 14:35:42 +00:00
parent faba13f7fd
commit 7a7ca1aab7
7 changed files with 428 additions and 546 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

261
code/mpie_spectral.f90
View File

@ -5,10 +5,10 @@
! written by P. Eisenlohr, ! written by P. Eisenlohr,
! F. Roters, ! F. Roters,
! L. Hantcherli, ! L. Hantcherli,
! W.A. Counts ! W.A. Counts,
! D.D. Tjahjanto ! D.D. Tjahjanto,
! C. Kords ! C. Kords,
! M. Diehl ! M. Diehl,
! R. Lebensohn ! R. Lebensohn
! !
! MPI fuer Eisenforschung, Duesseldorf ! MPI fuer Eisenforschung, Duesseldorf
@ -17,9 +17,9 @@
! Usage: ! Usage:
! - start program with mpie_spectral PathToGeomFile/NameOfGeom.geom ! - start program with mpie_spectral PathToGeomFile/NameOfGeom.geom
! PathToLoadFile/NameOfLoadFile.load ! PathToLoadFile/NameOfLoadFile.load
! - PathToLoadFile will be the working directory ! - PathToGeomFile will be the working directory
! - make sure the file "material.config" exists in the working ! - make sure the file "material.config" exists in the working
! directory ! directory. For further configuration use "numerics.config"
!******************************************************************** !********************************************************************
program mpie_spectral program mpie_spectral
!******************************************************************** !********************************************************************
@ -29,7 +29,8 @@ program mpie_spectral
use IO use IO
use math use math
use CPFEM, only: CPFEM_general, CPFEM_initAll use CPFEM, only: CPFEM_general, CPFEM_initAll
use numerics, only: mpieNumThreadsInt, err_div_tol, err_stress_tol, err_stress_tolrel, err_defgrad_tol, itmax use numerics, only: err_div_tol, err_stress_tol, err_stress_tolrel, err_defgrad_tol,&
itmax, fast_execution, mpieNumThreadsInt
use homogenization, only: materialpoint_sizeResults, materialpoint_results use homogenization, only: materialpoint_sizeResults, materialpoint_results
!$ use OMP_LIB ! the openMP function library !$ use OMP_LIB ! the openMP function library
@ -70,7 +71,7 @@ program mpie_spectral
real(pReal), dimension(3,3,3) :: temp333_Real real(pReal), dimension(3,3,3) :: temp333_Real
real(pReal), dimension(3,3,3,3) :: dPdF, c0, s0 real(pReal), dimension(3,3,3,3) :: dPdF, c0, s0
real(pReal), dimension(6) :: cstress ! cauchy stress in Mandel notation real(pReal), dimension(6) :: cstress ! cauchy stress in Mandel notation
real(pReal), dimension(6,6) :: dsde, c066, s066 real(pReal), dimension(6,6) :: dsde, c066, s066 ! Mandel notation of 4th order tensors
real(pReal), dimension(:,:,:), allocatable :: ddefgrad real(pReal), dimension(:,:,:), allocatable :: ddefgrad
real(pReal), dimension(:,:,:,:,:), allocatable :: pstress_field, defgrad, defgradold real(pReal), dimension(:,:,:,:,:), allocatable :: pstress_field, defgrad, defgradold
@ -79,27 +80,23 @@ program mpie_spectral
complex(pReal), dimension(3,3) :: temp33_Complex complex(pReal), dimension(3,3) :: temp33_Complex
real(pReal), dimension(3,3) :: xinormdyad real(pReal), dimension(3,3) :: xinormdyad
real(pReal), dimension(:,:,:,:,:,:,:), allocatable :: gamma_hat real(pReal), dimension(:,:,:,:,:,:,:), allocatable :: gamma_hat
real(pReal), dimension(:,:,:,:), allocatable :: xi
integer(pInt), dimension(3) :: k_s integer(pInt), dimension(3) :: k_s
real(pReal), dimension(3) :: xi, xi_middle
integer*8, dimension(2,3,3) :: plan_fft integer*8, dimension(2,3,3) :: plan_fft
! convergence etc. ! loop variables, convergence etc.
real(pReal) err_div, err_stress, err_defgrad, sigma0 real(pReal) guessmode, err_div, err_stress, err_defgrad, sigma0
integer(pInt) ierr
logical errmatinv
! loop variables etc.
real(pReal) guessmode ! flip-flop to guess defgrad fluctuation field evolution
integer(pInt) i, j, k, l, m, n, p integer(pInt) i, j, k, l, m, n, p
integer(pInt) loadcase, ielem, iter, calcmode, CPFEM_mode integer(pInt) loadcase, ielem, iter, calcmode, CPFEM_mode, ierr
logical errmatinv
real(pReal) temperature ! not used, but needed for call to CPFEM_general real(pReal) temperature ! not used, but needed for call to CPFEM_general
!Initializing !Initializing
!$ call omp_set_num_threads(mpieNumThreadsInt) ! set number of threads for parallel execution set by MPIE_NUM_THREADS !$ call omp_set_num_threads(mpieNumThreadsInt) ! set number of threads for parallel execution set by MPIE_NUM_THREADS
bc_maskvector = '' bc_maskvector = ''
unit = 234_pInt unit = 234_pInt
ones = 1.0_pReal; zeroes = 0.0_pReal ones = 1.0_pReal; zeroes = 0.0_pReal
N_l = 0_pInt; N_s = 0_pInt N_l = 0_pInt; N_s = 0_pInt
@ -184,10 +181,10 @@ program mpie_spectral
do i = 1, N_Loadcases do i = 1, N_Loadcases
if (any(bc_mask(:,:,1,i) == bc_mask(:,:,2,i))) call IO_error(47,i) ! bc_mask consistency if (any(bc_mask(:,:,1,i) == bc_mask(:,:,2,i))) call IO_error(47,i) ! bc_mask consistency
print '(a,/,3(3(f12.6,x)/))','L',math_transpose3x3(bc_velocityGrad(:,:,i)) print '(a,/,3(3(f12.6,x)/))','L' ,math_transpose3x3(bc_velocityGrad(:,:,i))
print '(a,/,3(3(f12.6,x)/))','bc_stress',math_transpose3x3(bc_stress(:,:,i)) print '(a,/,3(3(f12.6,x)/))','bc_stress',math_transpose3x3(bc_stress(:,:,i))
print '(a,/,3(3(l,x)/))','bc_mask for velocitygrad',transpose(bc_mask(:,:,1,i)) print '(a,/,3(3(l,x)/))', 'bc_mask for velocitygrad',transpose(bc_mask(:,:,1,i))
print '(a,/,3(3(l,x)/))','bc_mask for stress',transpose(bc_mask(:,:,2,i)) print '(a,/,3(3(l,x)/))', 'bc_mask for stress' ,transpose(bc_mask(:,:,2,i))
print *,'time',bc_timeIncrement(i) print *,'time',bc_timeIncrement(i)
print *,'incs',bc_steps(i) print *,'incs',bc_steps(i)
print *, '' print *, ''
@ -241,33 +238,15 @@ program mpie_spectral
print '(a,/,f6.1,f6.1,f6.1)','dimension x y z', geomdimension print '(a,/,f6.1,f6.1,f6.1)','dimension x y z', geomdimension
print *,'homogenization',homog print *,'homogenization',homog
allocate (workfft(resolution(1)/2+1,resolution(2),resolution(3),3,3)); workfft = 0.0_pReal allocate (defgrad(resolution(1),resolution(2),resolution(3),3,3)); defgrad = 0.0_pReal
allocate (gamma_hat(resolution(1)/2+1,resolution(2),resolution(3),3,3,3,3)); gamma_hat = 0.0_pReal allocate (defgradold(resolution(1),resolution(2),resolution(3),3,3)); defgradold = 0.0_pReal
allocate (xi(resolution(1)/2+1,resolution(2),resolution(3),3)); xi = 0.0_pReal
allocate (pstress_field(resolution(1),resolution(2),resolution(3),3,3)); pstress_field = 0.0_pReal
allocate (defgrad(resolution(1),resolution(2),resolution(3),3,3)); defgrad = 0.0_pReal
allocate (defgradold(resolution(1),resolution(2),resolution(3),3,3)); defgradold = 0.0_pReal
allocate (ddefgrad(resolution(1),resolution(2),resolution(3))); ddefgrad = 0.0_pReal
! Initialization of fftw (see manual on fftw.org for more details)
call dfftw_init_threads(ierr) !toDo: add error code
call dfftw_plan_with_nthreads(mpieNumThreadsInt)
do m = 1,3; do n = 1,3
call dfftw_plan_dft_r2c_3d(plan_fft(1,m,n),resolution(1),resolution(2),resolution(3),&
pstress_field(:,:,:,m,n), workfft(:,:,:,m,n), FFTW_PATIENT)
call dfftw_plan_dft_c2r_3d(plan_fft(2,m,n),resolution(1),resolution(2),resolution(3),&
workfft(:,:,:,m,n), ddefgrad(:,:,:), FFTW_PATIENT)
enddo; enddo
!try to make just one call instead of 9 for the r2c transform. Not working yet
! call dfftw_plan_many_dft_r2c_(plan_fft(1,1,1),3,(/resolution(1),resolution(2),resolution(3)/),9,&
!pstress_field(:,:,:,:,:),NULL,(/resolution(1),resolution(2),resolution(3)/),1, workfft(:,:,:,m,n),NULL, FFTW_PATIENT)
!
prodnn = resolution(1)*resolution(2)*resolution(3) prodnn = resolution(1)*resolution(2)*resolution(3)
wgt = 1_pReal/real(prodnn, pReal) wgt = 1.0_pReal/real(prodnn, pReal)
defgradAim = math_I3 defgradAim = math_I3
defgradAimOld = math_I3 defgradAimOld = math_I3
defgrad_av = math_I3 defgrad_av = math_I3
! Initialization of CPFEM_general (= constitutive law) and of deformation gradient field ! Initialization of CPFEM_general (= constitutive law) and of deformation gradient field
call CPFEM_initAll(temperature,1_pInt,1_pInt) call CPFEM_initAll(temperature,1_pInt,1_pInt)
ielem = 0_pInt ielem = 0_pInt
@ -279,41 +258,69 @@ program mpie_spectral
call CPFEM_general(2,math_I3,math_I3,temperature,0.0_pReal,ielem,1_pInt,cstress,dsde,pstress,dPdF) call CPFEM_general(2,math_I3,math_I3,temperature,0.0_pReal,ielem,1_pInt,cstress,dsde,pstress,dPdF)
c066 = c066 + dsde c066 = c066 + dsde
enddo; enddo; enddo enddo; enddo; enddo
c066 = c066 * wgt c066 = c066 * wgt
c0 = math_mandel66to3333(c066) c0 = math_mandel66to3333(c066)
call math_invert(6, c066, s066,i, errmatinv) call math_invert(6, c066, s066,i, errmatinv)
s0 = math_mandel66to3333(s066) s0 = math_mandel66to3333(s066)
!calculation of xinormdyad (to calculate gamma_hat) and xi (waves, for proof of equilibrium) !calculation of calculate gamma_hat field in the case of fast execution (needs a lot of memory)
do k = 1, resolution(3) if(fast_execution) then
k_s(3) = k-1 allocate (gamma_hat(resolution(1)/2+1,resolution(2),resolution(3),3,3,3,3)); gamma_hat = 0.0_pReal
if(k > resolution(3)/2+1) k_s(3) = k_s(3)-resolution(3) do k = 1, resolution(3)
do j = 1, resolution(2) k_s(3) = k-1
k_s(2) = j-1 if(k > resolution(3)/2+1) k_s(3) = k_s(3)-resolution(3)
if(j > resolution(2)/2+1) k_s(2) = k_s(2)-resolution(2) do j = 1, resolution(2)
do i = 1, resolution(1)/2+1 k_s(2) = j-1
k_s(1) = i-1 if(j > resolution(2)/2+1) k_s(2) = k_s(2)-resolution(2)
xi(i,j,k,3) = 0.0_pReal do i = 1, resolution(1)/2+1
if(resolution(3) > 1) xi(i,j,k,3) = real(k_s(3), pReal)/geomdimension(3) k_s(1) = i-1
xi(i,j,k,2) = real(k_s(2), pReal)/geomdimension(2) xi(3) = 0.0_pReal !for the 2D case
xi(i,j,k,1) = real(k_s(1), pReal)/geomdimension(1) if(resolution(3) > 1) xi(3) = real(k_s(3), pReal)/geomdimension(3) !3D case
if (any(xi(i,j,k,:) /= 0.0_pReal)) then xi(2) = real(k_s(2), pReal)/geomdimension(2)
do l = 1,3; do m = 1,3 xi(1) = real(k_s(1), pReal)/geomdimension(1)
xinormdyad(l,m) = xi(i,j,k, l)*xi(i,j,k, m)/sum(xi(i,j,k,:)**2) if (any(xi /= 0.0_pReal)) then
enddo; enddo do l = 1,3; do m = 1,3
else xinormdyad(l,m) = xi(l)*xi(m)/sum(xi**2)
xinormdyad = 0.0_pReal enddo; enddo
endif else
temp33_Real = math_mul3333xx33(c0, xinormdyad) xinormdyad = 0.0_pReal
temp33_Real = math_inv3x3(temp33_Real) endif
do l=1,3; do m=1,3; do n=1,3; do p=1,3 temp33_Real = math_mul3333xx33(c0, xinormdyad)
gamma_hat(i,j,k, l,m,n,p) = - (0.5*temp33_Real(l,n)+0.5*temp33_Real(n,l)) *& temp33_Real = math_inv3x3(temp33_Real)
(0.5*xinormdyad(m,p)+0.5*xinormdyad(p,m)) do l=1,3; do m=1,3; do n=1,3; do p=1,3
enddo; enddo; enddo; enddo gamma_hat(i,j,k, l,m,n,p) = - (0.5*temp33_Real(l,n)+0.5*temp33_Real(n,l)) *&
enddo; enddo; enddo (0.5*xinormdyad(m,p)+0.5*xinormdyad(p,m))
enddo; enddo; enddo; enddo
enddo; enddo; enddo
else ! or allocate just one fourth order tensor
allocate (gamma_hat(1,1,1,3,3,3,3)); gamma_hat = 0.0_pReal
endif
! calculate xi for the calculation of divergence in Fourier space (middle frequency)
xi_middle(3) = 0.0_pReal
if(resolution(3) > 1) xi_middle(3) = real(resolution(3)/2, pReal)/geomdimension(3) !3D case
xi_middle(2) = real(resolution(2)/2, pReal)/geomdimension(2)
xi_middle(1) = real(resolution(1)/2, pReal)/geomdimension(1)
! Initialization of fftw (see manual on fftw.org for more details)
allocate (workfft(resolution(1)/2+1,resolution(2),resolution(3),3,3)); workfft = 0.0_pReal
allocate (ddefgrad(resolution(1),resolution(2),resolution(3))); ddefgrad = 0.0_pReal
allocate (pstress_field(resolution(1),resolution(2),resolution(3),3,3)); pstress_field = 0.0_pReal
call dfftw_init_threads(ierr) !toDo: add error code
call dfftw_plan_with_nthreads(mpieNumThreadsInt)
! Do r2c Transform r2c in one step
call dfftw_plan_many_dft_r2c(plan_fft(1,1,1),3,(/resolution(1),resolution(2),resolution(3)/),9,&
pstress_field(:,:,:,:,:),(/resolution(1),resolution(2),resolution(3)/),1,prodnn,workfft(:,:,:,:,:),&
(/resolution(1)/2+1,resolution(2),resolution(3)/),1,(resolution(1)/2+1)*resolution(2)*resolution(3),FFTW_PATIENT)
do m = 1,3; do n = 1,3 ! do the back transform for each single component (saves memory)
call dfftw_plan_dft_c2r(plan_fft(2,m,n),3,(/resolution(1),resolution(2),resolution(3)/),&
workfft(:,:,:,m,n), ddefgrad(:,:,:), FFTW_PATIENT)
enddo; enddo
! write header of output file ! write header of output file
open(538,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//'.spectralOut',form='UNFORMATTED') open(538,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())&
//'_'//trim(getLoadcase())&
//'.spectralOut',form='UNFORMATTED')
write(538), 'load',trim(getLoadcaseName()) write(538), 'load',trim(getLoadcaseName())
write(538), 'workingdir',trim(getSolverWorkingDirectoryName()) write(538), 'workingdir',trim(getSolverWorkingDirectoryName())
write(538), 'geometry',trim(getSolverJobName())//InputFileExtension write(538), 'geometry',trim(getSolverJobName())//InputFileExtension
@ -373,7 +380,8 @@ program mpie_spectral
(err_div > err_div_tol .or. & (err_div > err_div_tol .or. &
err_stress > err_stress_tol .or. & err_stress > err_stress_tol .or. &
err_defgrad > err_defgrad_tol)) err_defgrad > err_defgrad_tol))
iter = iter + 1 iter = iter + 1_pInt
print*, ' '
print '(3(A,I5.5,tr2))', ' Loadcase = ',loadcase, ' Step = ',steps,'Iteration = ',iter print '(3(A,I5.5,tr2))', ' Loadcase = ',loadcase, ' Step = ',steps,'Iteration = ',iter
cstress_av = 0.0_pReal cstress_av = 0.0_pReal
!************************************************************* !*************************************************************
@ -407,7 +415,7 @@ program mpie_spectral
pstress_av(m,n) = sum(pstress_field(:,:,:,m,n)) * wgt pstress_av(m,n) = sum(pstress_field(:,:,:,m,n)) * wgt
defgrad_av(m,n) = sum(defgrad(:,:,:,m,n)) * wgt defgrad_av(m,n) = sum(defgrad(:,:,:,m,n)) * wgt
enddo; enddo enddo; enddo
err_stress = maxval(abs(mask_stress * (pstress_av - bc_stress(:,:,loadcase)))) err_stress = maxval(abs(mask_stress * (pstress_av - bc_stress(:,:,loadcase))))
err_stress_tol = maxval(abs(pstress_av))*err_stress_tolrel err_stress_tol = maxval(abs(pstress_av))*err_stress_tolrel
@ -431,19 +439,19 @@ program mpie_spectral
err_div = 2 * err_div_tol err_div = 2 * err_div_tol
err_defgrad = maxval(abs(mask_defgrad * (defgrad_av - defgradAim))) err_defgrad = maxval(abs(mask_defgrad * (defgrad_av - defgradAim)))
print '(a,/,3(3(f12.7,x)/))', ' Deformation Gradient: ',math_transpose3x3(defgrad_av) print '(a,/,3(3(f12.7,x)/))', ' Deformation Gradient: ',math_transpose3x3(defgrad_av)
print '(a,/,3(3(f10.4,x)/))', ' Cauchy Stress [MPa]: ',math_transpose3x3(cstress_av)/1.e6 print '(a,/,3(3(f10.4,x)/))', ' Cauchy Stress [MPa]: ' ,math_transpose3x3(cstress_av)/1.e6
print '(2(a,E8.2))', ' error stress ',err_stress,' Tol. = ', err_stress_tol print '(2(a,E8.2))', ' error stress ',err_stress, ' Tol. = ', err_stress_tol
print '(2(a,E8.2))', ' error deformation gradient ',err_defgrad,' Tol. = ', err_defgrad_tol*0.8 print '(2(a,E8.2))', ' error deformation gradient ',err_defgrad,' Tol. = ', err_defgrad_tol*0.8
if(err_stress < err_stress_tol*0.8) then if(err_stress < err_stress_tol*0.8) then
calcmode = 1_pInt calcmode = 1_pInt
endif endif
! Using the spectral method to calculate the change of deformation gradient, check divergence of stress field in fourier space ! Using the spectral method to calculate the change of deformation gradient, check divergence of stress field in fourier space
case (1) case (1)
print *, 'Update Stress Field (constitutive evaluation P(F))' print *, 'Update Stress Field (constitutive evaluation P(F))'
ielem = 0_pInt ielem = 0_pInt
do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1) do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1)
ielem = ielem + 1 ielem = ielem + 1_pInt
call CPFEM_general(3, defgradold(i,j,k,:,:), defgrad(i,j,k,:,:),& call CPFEM_general(3, defgradold(i,j,k,:,:), defgrad(i,j,k,:,:),&
temperature,timeinc,ielem,1_pInt,& temperature,timeinc,ielem,1_pInt,&
cstress,dsde, pstress, dPdF) cstress,dsde, pstress, dPdF)
@ -451,7 +459,7 @@ program mpie_spectral
ielem = 0_pInt ielem = 0_pInt
do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1) do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1)
ielem = ielem + 1 ielem = ielem + 1_pInt
call CPFEM_general(2,& call CPFEM_general(2,&
defgradold(i,j,k,:,:), defgrad(i,j,k,:,:),& defgradold(i,j,k,:,:), defgrad(i,j,k,:,:),&
temperature,timeinc,ielem,1_pInt,& temperature,timeinc,ielem,1_pInt,&
@ -459,46 +467,83 @@ program mpie_spectral
pstress_field(i,j,k,:,:) = pstress pstress_field(i,j,k,:,:) = pstress
cstress_av = cstress_av + math_mandel6to33(cstress) cstress_av = cstress_av + math_mandel6to33(cstress)
enddo; enddo; enddo enddo; enddo; enddo
cstress_av = cstress_av * wgt
do m = 1,3; do n = 1,3
pstress_av(m,n) = sum(pstress_field(:,:,:,m,n))*wgt
enddo; enddo
print *, 'Calculating equilibrium using spectral method' print *, 'Calculating equilibrium using spectral method'
err_div = 0.0_pReal; sigma0 = 0.0_pReal err_div = 0.0_pReal; sigma0 = 0.0_pReal
do m = 1,3; do n = 1,3 call dfftw_execute_dft_r2c(plan_fft(1,1,1), pstress_field,workfft) ! FFT of pstress
call dfftw_execute_dft_r2c(plan_fft(1,m,n), pstress_field(:,:,:,m,n),workfft(:,:,:,m,n))
if(n==3) sigma0 = max(sigma0, sum(abs(workfft(1,1,1,m,:)))) ! L infinity Norm of stress tensor do m = 1,3
enddo; enddo sigma0 = max(sigma0, sum(abs(workfft(1,1,1,m,:)))) ! L infinity Norm of stress tensor
err_div = (maxval(abs(math_mul33x3_complex(workfft(resolution(1)/2+1,resolution(2)/2+1,resolution(3)/2+1,:,:),& enddo
xi(resolution(1)/2+1,resolution(2)/2+1,resolution(3)/2+1,:))))) ! L infinity Norm of div(stress) err_div = (maxval(abs(math_mul33x3_complex(workfft(resolution(1)/2+1,resolution(2)/2+1,resolution(3)/2+1,:,:),xi_middle)))) ! L infinity Norm of div(stress)
err_div = err_div/sigma0 !weighting of error
if(fast_execution) then ! fast execution with stored gamma_hat
do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1)/2+1 do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1)/2+1
temp33_Complex = 0.0_pReal temp33_Complex = 0.0_pReal
do m = 1,3; do n = 1,3 do m = 1,3; do n = 1,3
temp33_Complex(m,n) = sum(gamma_hat(i,j,k,m,n,:,:) * workfft(i,j,k,:,:)) temp33_Complex(m,n) = sum(gamma_hat(i,j,k,m,n,:,:) * workfft(i,j,k,:,:))
enddo; enddo enddo; enddo
workfft(i,j,k,:,:) = temp33_Complex(:,:) workfft(i,j,k,:,:) = temp33_Complex(:,:)
enddo; enddo; enddo enddo; enddo; enddo
workfft(1,1,1,:,:) = defgrad_av - math_I3
else ! memory saving version, in-time calculation of gamma_hat
do k = 1, resolution(3)
k_s(3) = k-1
if(k > resolution(3)/2+1) k_s(3) = k_s(3)-resolution(3)
do j = 1, resolution(2)
k_s(2) = j-1
if(j > resolution(2)/2+1) k_s(2) = k_s(2)-resolution(2)
do i = 1, resolution(1)/2+1
k_s(1) = i-1
xi(3) = 0.0_pReal !for the 2D case
if(resolution(3) > 1) xi(3) = real(k_s(3), pReal)/geomdimension(3) !3D case
xi(2) = real(k_s(2), pReal)/geomdimension(2)
xi(1) = real(k_s(1), pReal)/geomdimension(1)
if (any(xi(:) /= 0.0_pReal)) then
do l = 1,3; do m = 1,3
xinormdyad(l,m) = xi(l)*xi( m)/sum(xi**2)
enddo; enddo
else
xinormdyad = 0.0_pReal
endif
temp33_Real = math_mul3333xx33(c0, xinormdyad)
temp33_Real = math_inv3x3(temp33_Real)
do l=1,3; do m=1,3; do n=1,3; do p=1,3
gamma_hat(1,1,1, l,m,n,p) = - (0.5*temp33_Real(l,n)+0.5*temp33_Real(n,l)) *&
(0.5*xinormdyad(m,p)+0.5*xinormdyad(p,m))
enddo; enddo; enddo; enddo
temp33_Complex = 0.0_pReal
do m = 1,3; do n = 1,3
temp33_Complex(m,n) = sum(gamma_hat(1,1,1, m,n,:,:) * workfft(i,j,k,:,:))
enddo; enddo
workfft(i,j,k,:,:) = temp33_Complex(:,:)
enddo
enddo
enddo
endif
err_div = err_div/sigma0 !weighting of error workfft(1,1,1,:,:) = defgrad_av - math_I3 !zero frequency
cstress_av = cstress_av * wgt
do m = 1,3; do n = 1,3 do m = 1,3; do n = 1,3
call dfftw_execute_dft_c2r(plan_fft(2,m,n), workfft(:,:,:,m,n),ddefgrad(:,:,:)) call dfftw_execute_dft_c2r(plan_fft(2,m,n), workfft(:,:,:,m,n),ddefgrad(:,:,:))
defgrad(:,:,:,m,n) = defgrad(:,:,:,m,n) + ddefgrad * wgt defgrad(:,:,:,m,n) = defgrad(:,:,:,m,n) + ddefgrad * wgt
pstress_av(m,n) = sum(pstress_field(:,:,:,m,n))*wgt defgrad_av(m,n) = sum(defgrad(:,:,:,m,n))*wgt
defgrad_av(m,n) = sum(defgrad(:,:,:,m,n))*wgt defgrad(:,:,:,m,n) = defgrad(:,:,:,m,n) + (defgradAim(m,n) - defgrad_av(m,n)) !anticipated target minus current state
defgrad(:,:,:,m,n) = defgrad(:,:,:,m,n) + (defgradAim(m,n) - defgrad_av(m,n)) !anticipated target minus current state
enddo; enddo enddo; enddo
err_stress = maxval(abs(mask_stress * (pstress_av - bc_stress(:,:,loadcase)))) err_stress = maxval(abs(mask_stress * (pstress_av - bc_stress(:,:,loadcase))))
err_stress_tol = maxval(abs(pstress_av))*err_stress_tolrel !accecpt relativ error specified err_stress_tol = maxval(abs(pstress_av))*err_stress_tolrel !accecpt relativ error specified
err_defgrad = maxval(abs(mask_defgrad * (defgrad_av - defgradAim))) err_defgrad = maxval(abs(mask_defgrad * (defgrad_av - defgradAim)))
print '(2(a,E8.2))', ' error divergence ',err_div,' Tol. = ', err_div_tol print '(2(a,E8.2))', ' error divergence ',err_div, ' Tol. = ', err_div_tol
print '(2(a,E8.2))', ' error stress ',err_stress,' Tol. = ', err_stress_tol print '(2(a,E8.2))', ' error stress ',err_stress, ' Tol. = ', err_stress_tol
print '(2(a,E8.2))', ' error deformation gradient ',err_defgrad,' Tol. = ', err_defgrad_tol print '(2(a,E8.2))', ' error deformation gradient ',err_defgrad,' Tol. = ', err_defgrad_tol
if((err_stress > err_stress_tol .or. err_defgrad > err_defgrad_tol) .and. err_div < err_div_tol) then ! change to calculation of BCs, reset damper etc. if((err_stress > err_stress_tol .or. err_defgrad > err_defgrad_tol) .and. err_div < err_div_tol) then ! change to calculation of BCs, reset damper etc.
calcmode = 0 calcmode = 0_pInt
defgradAimCorr = 0.0_pReal defgradAimCorr = 0.0_pReal
damper = damper * 0.9_pReal damper = damper * 0.9_pReal
endif endif
@ -508,9 +553,9 @@ program mpie_spectral
write(538) materialpoint_results(:,1,:) !write to output file write(538) materialpoint_results(:,1,:) !write to output file
print '(a,x,f12.7)' , ' Determinant of Deformation Aim:', math_det3x3(defgradAim) print '(a,x,f12.7)' , ' Determinant of Deformation Aim:', math_det3x3(defgradAim)
print '(a,/,3(3(f12.7,x)/))', ' Deformation Aim: ',defgradAim(1:3,:) print '(a,/,3(3(f12.7,x)/))', ' Deformation Aim: ',math_transpose3x3(defgradAim)
print '(a,/,3(3(f12.7,x)/))', ' Deformation Gradient: ',defgrad_av(1:3,:) print '(a,/,3(3(f12.7,x)/))', ' Deformation Gradient: ',math_transpose3x3(defgrad_av)
print '(a,/,3(3(f10.4,x)/))', ' Cauchy Stress [MPa]: ',cstress_av(1:3,:)/1.e6 print '(a,/,3(3(f10.4,x)/))', ' Cauchy Stress [MPa]: ',math_transpose3x3(cstress_av)/1.e6
print '(A)', '************************************************************' print '(A)', '************************************************************'
enddo ! end looping over steps in current loadcase enddo ! end looping over steps in current loadcase
enddo ! end looping over loadcases enddo ! end looping over loadcases

29
code/mpie_spectral_interface.f90
View File

@ -34,7 +34,7 @@ function getSolverWorkingDirectoryName()
character(len=1024) cwd,outname,getSolverWorkingDirectoryName character(len=1024) cwd,outname,getSolverWorkingDirectoryName
character(len=*), parameter :: pathSep = achar(47)//achar(92) ! forwardslash, backwardslash character(len=*), parameter :: pathSep = achar(47)//achar(92) ! forwardslash, backwardslash
call getarg(2,outname) ! path to loadFile call getarg(1,outname) ! path to loadFile
if (scan(outname,pathSep) == 1) then ! absolute path given as command line argument if (scan(outname,pathSep) == 1) then ! absolute path given as command line argument
getSolverWorkingDirectoryName = outname(1:scan(outname,pathSep,back=.true.)) getSolverWorkingDirectoryName = outname(1:scan(outname,pathSep,back=.true.))
@ -84,6 +84,33 @@ function getSolverJobName()
return return
endfunction endfunction
!********************************************************************
! name of load case file exluding extension
!
!********************************************************************
function getLoadCase()
use prec, only: pInt
implicit none
character(1024) getLoadCase, outName
character(len=*), parameter :: pathSep = achar(47)//achar(92) ! /, \
integer(pInt) posExt,posSep
getLoadCase = ''
posSep=1
call getarg(2,outName)
posExt = scan(outName,'.',back=.true.)
posSep = scan(outName,pathSep,back=.true.)
if (posExt <= posSep) posExt = len_trim(outName)+1 ! no extension present
getLoadCase = outName(posSep+1:posExt-1) ! name of load case file exluding extension
return
endfunction
!******************************************************************** !********************************************************************
! relative path of loadcase from command line arguments ! relative path of loadcase from command line arguments

6
code/numerics.f90
View File

@ -50,8 +50,10 @@ real(pReal) relevantStrain, & ! strain
err_stress_tol, & ! absolut stress error, will be computed from err_stress_tolrel (dont prescribe a value) err_stress_tol, & ! absolut stress error, will be computed from err_stress_tolrel (dont prescribe a value)
err_stress_tolrel, & ! factor to multiply with highest stress to get err_stress_tol err_stress_tolrel, & ! factor to multiply with highest stress to get err_stress_tol
err_defgrad_tol ! tolerance for error of defgrad compared to prescribed defgrad err_defgrad_tol ! tolerance for error of defgrad compared to prescribed defgrad
logical fast_execution ! for fast execution (pre calculation of gamma_hat)
integer(pInt) itmax , & ! maximum number of iterations integer(pInt) itmax , & ! maximum number of iterations
!* Random seeding parameters !* Random seeding parameters
fixedSeed ! fixed seeding for pseudo-random number generator fixedSeed ! fixed seeding for pseudo-random number generator
!* OpenMP variable !* OpenMP variable
@ -144,6 +146,7 @@ subroutine numerics_init()
err_defgrad_tol = 1.0e-3 err_defgrad_tol = 1.0e-3
err_stress_tolrel = 0.01 err_stress_tolrel = 0.01
itmax = 20_pInt itmax = 20_pInt
fast_execution = .false.
!* Random seeding parameters: added <<<updated 27.08.2009>>> !* Random seeding parameters: added <<<updated 27.08.2009>>>
fixedSeed = 0_pInt fixedSeed = 0_pInt
@ -253,6 +256,8 @@ subroutine numerics_init()
err_stress_tolrel = IO_floatValue(line,positions,2) err_stress_tolrel = IO_floatValue(line,positions,2)
case ('itmax') case ('itmax')
itmax = IO_intValue(line,positions,2) itmax = IO_intValue(line,positions,2)
case ('fast_execution')
fast_execution = IO_intValue(line,positions,2) > 0_pInt
!* Random seeding parameters !* Random seeding parameters
case ('fixed_seed') case ('fixed_seed')
@ -317,6 +322,7 @@ subroutine numerics_init()
write(6,'(a24,x,e8.1)') 'err_defgrad_tol: ',err_defgrad_tol write(6,'(a24,x,e8.1)') 'err_defgrad_tol: ',err_defgrad_tol
write(6,'(a24,x,e8.1)') 'err_stress_tolrel: ',err_stress_tolrel write(6,'(a24,x,e8.1)') 'err_stress_tolrel: ',err_stress_tolrel
write(6,'(a24,x,i8)') 'itmax: ',itmax write(6,'(a24,x,i8)') 'itmax: ',itmax
write(6,'(a24,x,L8)') 'fast_execution: ',fast_execution
write(6,*) write(6,*)
!* Random seeding parameters !* Random seeding parameters

62
processing/post/3Dvisualize
View File

@ -277,6 +277,68 @@ def vtk_writeASCII_mesh(mesh,data,res):
[[['\t'.join(map(str,data[type][item][:,j,k]))] for j in range(res[1])] for k in range(res[2])], [[['\t'.join(map(str,data[type][item][:,j,k]))] for j in range(res[1])] for k in range(res[2])],
] ]
# vtk = open(filename, 'w')
# output(cmd,{'filepointer':vtk},'File')
# vtk.close()
return cmds
# ++++++++++++++++++++++++++++++++++++++++++++++++++++
def gmsh_writeASCII_mesh(mesh,data,res):
# ++++++++++++++++++++++++++++++++++++++++++++++++++++
""" function writes data array defined on a hexahedral mesh (geometry) """
N1 = (res[0]+1)*(res[1]+1)*(res[2]+1)
N = res[0]*res[1]*res[2]
cmds = [\
'$MeshFormat',
'2.1 0 8',
'$EndMeshFormat',
'$Nodes',
'%i float'%N1,
[[['\t'.join(map(str,l,mesh[i,j,k])) for l in range(1,N1+1) for i in range(res[0]+1)] for j in range(res[1]+1)] for k in range(res[2]+1)],
'$EndNodes',
'$Elements',
'%i'%N,
]
n_elem = 0
for i in range (res[2]):
for j in range (res[1]):
for k in range (res[0]):
base = i*(res[1]+1)*(res[2]+1)+j*(res[1]+1)+k
n_elem +=1
cmds.append('\t'.join(map(str,[ \
n_elem,
'5',
base,
base+1,
base+res[1]+2,
base+res[1]+1,
base+(res[1]+1)*(res[2]+1),
base+(res[1]+1)*(res[2]+1)+1,
base+(res[1]+1)*(res[2]+1)+res[1]+2,
base+(res[1]+1)*(res[2]+1)+res[1]+1,
])))
cmds += [\
'ElementData',
'1',
'%s'%item, # name of the view
'0.0', # thats the time value
'3',
'0', # time step
'1',
'%i'%N
]
for type in data:
for item in data[type]:
cmds += [\
'%s %s float'%(type.upper(),item),
'LOOKUP_TABLE default',
[[['\t'.join(map(str,data[type][item][:,j,k]))] for j in range(res[1])] for k in range(res[2])],
]
# vtk = open(filename, 'w') # vtk = open(filename, 'w')
# output(cmd,{'filepointer':vtk},'File') # output(cmd,{'filepointer':vtk},'File')
# vtk.close() # vtk.close()

9
processing/post/make_reconstruct.py
View File

@ -5,4 +5,11 @@
# module reconstruct.so out of the fortran code reconstruct.f90 # module reconstruct.so out of the fortran code reconstruct.f90
# written by M. Diehl, m.diehl@mpie.de # written by M. Diehl, m.diehl@mpie.de
f2py -c -m reconstruct reconstruct.f90
#f2py -m reconstruct2 -h reconstruct2.pyf --overwrite-signature reconstruct_new.f90
#f2py -m reconstruct -h reconstruct.pyf reconstruct.f90
# f2py -c \
# --f90flags="-heap-arrays 500000000" \ # preventing segmentation fault for large arrays
# reconstruct.pyf \
# reconstruct.f90
f2py -c --f90flags="-heap-arrays 500000000" reconstruct.pyf reconstruct.f90

578
processing/post/reconstruct.f90
View File

@ -1,421 +1,157 @@
! -*- f90 -*- !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
function coordinates2(res_x,res_y,res_z,geomdimension,defgrad) subroutine mesh(res_x,res_y,res_z,geomdim,defgrad_av,centroids,nodes)
implicit none !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
integer, parameter :: pDouble = selected_real_kind(15,50) implicit none
integer i,j,k, l,m, s,o, loop, res_x, res_y, res_z real*8 geomdim(3)
integer, dimension(3) :: res, init, oppo, me, rear integer res_x, res_y, res_z
real*8, dimension(3,3) :: defgrad_av real*8 wrappedCentroids(res_x+2,res_y+2,res_z+2,3)
integer, dimension(3) :: resolution real*8 nodes(res_x+1,res_y+1,res_z+1,3)
real*8, dimension(3) :: geomdimension, myStep real*8 centroids(res_x ,res_y ,res_z ,3)
real*8, dimension(3,3) :: temp33_Real
real*8, dimension(3,res_x, res_y, res_z) :: coordinates2 integer, dimension(3,8) :: neighbor = reshape((/ &
real*8, dimension(3,3,res_x, res_y, res_z) :: defgrad 0, 0, 0,&
real*8, dimension(3,res_x,res_y,res_z,8) :: cornerCoords 1, 0, 0,&
real*8, dimension(3,2+res_x,2+res_y,2+res_z,6,8) :: coord 1, 1, 0,&
!f2py intent(in) res_x 0, 1, 0,&
!f2py intent(in) res_y 0, 0, 1,&
!f2py intent(in) res_z 1, 0, 1,&
!f2py intent(in) geomdimension 1, 1, 1,&
!f2py intent(in) defgrad 0, 1, 1 &
!f2py intent(out) coordinates2 /), &
!f2py depend(res_x, res_y, res_z) coordinates2 (/3,8/))
!f2py depend(res_x, res_y, res_z) defgrad
integer i,j,k,n
! integer, dimension(3,8) :: corner = reshape((/ & real*8, dimension(3,3) :: defgrad_av
! 0, 0, 0,& integer, dimension(3) :: diag, shift, lookup, me, res
! 1, 0, 0,&
! 1, 1, 0,& diag = 1
! 0, 1, 0,& shift = 0
! 1, 1, 1,& lookup = 0
! 0, 1, 1,&
! 0, 0, 1,& res = (/res_x,res_y,res_z/)
! 1, 0, 1 &
! /), & wrappedCentroids=0.0
! (/3,8/)) wrappedCentroids(2:res_x+1,2:res_y+1,2:res_z+1,:) = centroids
! integer, dimension(3,8) :: step = reshape((/ & do k=0, res_z+1
! 1, 1, 1,& do j=0, res_y+1
! -1, 1, 1,& do i=0, res_x+1
! -1,-1, 1,& if (k==0 .or. k==res_z+1 .or. &
! 1,-1, 1,& j==0 .or. j==res_y+1 .or. &
! -1,-1,-1,& i==0 .or. i==res_x+1 ) then
! 1,-1,-1,& me = (/i,j,k/)
! 1, 1,-1,& shift = sign(abs(res+diag-2*me)/(res+diag),res+diag-2*me)
! -1, 1,-1 & lookup = me-diag+shift*res
! /), & wrappedCentroids(i+1,j+1,k+1,:) = centroids(lookup(1)+1,lookup(2)+1,lookup(3)+1,:)- &
! (/3,8/)) matmul(defgrad_av, shift*geomdim)
endif
! integer, dimension(3,6) :: order = reshape((/ & enddo; enddo; enddo
! 1, 2, 3,& do k=0, res_z
! 1, 3, 2,& do j=0, res_y
! 2, 1, 3,& do i=0, res_x
! 2, 3, 1,& do n=1,8
! 3, 1, 2,& nodes(i+1,j+1,k+1,:) = nodes(i+1,j+1,k+1,:) + wrappedCentroids(i+1+neighbor(n,1),j+1+neighbor(n,2),k+1+neighbor(n,3),:)
! 3, 2, 1 & enddo; enddo; enddo; enddo
! /), & nodes = nodes/8.0
! (/3,6/))
end subroutine mesh
resolution = (/res_x,res_y,res_z/)
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
write(6,*) 'defgrad', defgrad subroutine deformed(res_x,res_y,res_z,geomdim,defgrad,defgrad_av,coord_avgCorner)
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
do i=1, 3; do j=1,3 implicit none
defgrad_av(i,j) = sum(defgrad(i,j,:,:,:)) /real(res_x*res_y*res_z) real*8 geomdim(3)
enddo; enddo integer res_x, res_y, res_z
real*8 coord(8,6,res_x,res_y,res_z,3)
do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1) real*8 coord_avgOrder(8,res_x,res_y,res_z,3)
if((k==1).and.(j==1).and.(i==1)) then real*8 coord_avgCorner(res_x,res_y,res_z,3)
temp33_Real = real(0.0) real*8 defgrad(res_x,res_y,res_z,3,3)
else integer, dimension(3,8) :: corner = reshape((/ &
if((j==1).and.(i==1)) then 0, 0, 0,&
temp33_Real(1,:) = temp33_Real(1,:) + matmul(defgrad(:,:,i,j,k),& 1, 0, 0,&
(/real(0.0),real(0.0),real(geomdimension(3))/(real(resolution(3)))/)) 1, 1, 0,&
temp33_Real(2,:) = temp33_Real(1,:) 0, 1, 0,&
temp33_Real(3,:) = temp33_Real(1,:) 1, 1, 1,&
coordinates2(:,i,j,k) = temp33_Real(1,:) 0, 1, 1,&
else 0, 0, 1,&
if(i==1) then 1, 0, 1 &
temp33_Real(2,:) = temp33_Real(2,:) + matmul(defgrad(:,:,i,j,k),& /), &
(/real(0.0),real(geomdimension(2))/(real(resolution(2))),real(0.0)/)) (/3,8/))
temp33_Real(3,:) = temp33_Real(2,:) integer, dimension(3,8) :: step = reshape((/ &
coordinates2(:,i,j,k) = temp33_Real(2,:) 1, 1, 1,&
else -1, 1, 1,&
temp33_Real(3,:) = temp33_Real(3,:) + matmul(defgrad(:,:,i,j,k),& -1,-1, 1,&
(/real(geomdimension(1))/(real(resolution(1))),real(0.0),real(0.0)/)) 1,-1, 1,&
coordinates2(:,i,j,k) = temp33_Real(3,:) -1,-1,-1,&
endif 1,-1,-1,&
endif 1, 1,-1,&
endif -1, 1,-1 &
enddo; enddo; enddo /), &
do i=1, res_x; do j = 1, res_y; do k = 1, res_z (/3,8/))
coordinates2(:,i,j,k) = coordinates2(:,i,j,k)+ matmul(defgrad_av,(/geomdimension(1)/real(res_x),geomdimension(2)/real(res_y),geomdimension(3)/real(res_z)/)) integer, dimension(3,6) :: order = reshape((/ &
enddo; enddo; enddo 1, 2, 3,&
1, 3, 2,&
res = (/res_x,res_y,res_z/) 2, 1, 3,&
do i=1,3; do j=1,3 2, 3, 1,&
defgrad_av(i,j) = sum(defgrad(i,j,:,:,:)) /real(res(1)*res(2)*res(3)) 3, 1, 2,&
enddo; enddo 3, 2, 1 &
/), &
! do s = 1,8 (/3,6/))
! init = corner(:,s)*(res-(/1,1,1/))
! oppo = corner(:,1+mod(s-1+4,8))*(res-(/1,1,1/)) real*8 myStep(3), fones(3), parameter_coords(3)
! do o = 1,6 real*8 defgrad_av(3,3)
! do k = init(order(3,o)),oppo(order(3,o)),step(order(3,o),s) real*8 negative(3), positive(3)
! rear(order(2,o)) = init(order(2,o)) integer rear(3), init(3), ones(3), oppo(3), me(3), res(3)
! do j = init(order(2,o)),oppo(order(2,o)),step(order(2,o),s) integer i, j, k, s, o
! rear(order(1,o)) = init(order(1,o))
! do i = init(order(1,o)),oppo(order(1,o)),step(order(1,o),s) ones = 1
! ! print*, order(:,o) fones = 1.0
! me(order(1,o)) = i coord_avgOrder=0.0
! me(order(2,o)) = j
! me(order(3,o)) = k res = (/res_x,res_y,res_z/)
! ! print*, me
! ! if (all(me == init)) then do s = 0, 7 ! corners (from 0 to 7)
! ! coord(:,1+me(1),1+me(2),1+me(3),o,s) = 0.0 !& init = corner(:,s+1)*(res-ones) +ones
! ! geomdimension*(matmul(defgrad_av,real(corner(:,s),pDouble)) + & oppo = corner(:,mod((s+4),8)+1)*(res-ones) +ones
! ! matmul(defGrad(:,:,1+me(1),1+me(2),1+me(3)),step(:,s)/res/2.0_pDouble)) do o=1,6 ! orders ! from 1 to 6)
! ! else do k = init(order(3,o)), oppo(order(3,o)), step(order(3,o),s+1)
! ! myStep = (me-rear)*geomdimension/res rear(order(2,o)) = init(order(2,o))
! ! coord(:,1+me(1),1+me(2),1+me(3),o,s) = coord(:,1+rear(1),1+rear(2),1+rear(3),o,s) + & do j = init(order(2,o)), oppo(order(2,o)), step(order(2,o),s+1)
! ! 0.5_pDouble*matmul(defGrad(:,:,1+me(1),1+me(2),1+me(3)) + & rear(order(1,o)) = init(order(1,o))
! ! defGrad(:,:,1+rear(1),1+rear(2),1+rear(3)),& do i = init(order(1,o)), oppo(order(1,o)), step(order(1,o),s+1)
! ! myStep) me(order(1,o)) = i
! ! endif me(order(2,o)) = j
! ! rear = me me(order(3,o)) = k
! enddo if ( (me(1)==init(1)).and.(me(2)==init(2)).and. (me(3)==init(3)) ) then
! enddo coord(s+1,o,me(1),me(2),me(3),:) = geomdim * (matmul(defgrad_av,corner(:,s+1)) + &
! enddo matmul(defgrad(me(1),me(2),me(3),:,:),0.5*step(:,s+1)/res))
! enddo ! orders
! ! cornerCoords(:,:,:,:,s) = sum(coord(:,:,:,:,:,s),5)/6.0_pDouble else
! enddo ! corners myStep = (me-rear)*geomdim/res
coord(s+1,o,me(1),me(2),me(3),:) = coord(s+1,o,rear(1),rear(2),rear(3),:) + &
! coordinates = sum(cornerCoords(:,:,:,:,:),5)/8.0_pDouble ! plain average no shape functions... 0.5*matmul(defgrad(me(1),me(2),me(3),:,:) + &
end function defgrad(rear(1),rear(2),rear(3),:,:),myStep)
endif
rear = me
enddo; enddo; enddo; enddo
! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ do i=1,6
subroutine coordinates5(res_x,res_y,res_z,geomdimension,defgrad) coord_avgOrder(s+1,:,:,:,:) = coord_avgOrder(s+1,:,:,:,:) + coord(s+1,i,:,:,:,:)/6.0
! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ enddo
implicit none enddo
integer, parameter :: pDouble = selected_real_kind(15,50)
integer i,j,k, l,m, s,o, loop, res_x, res_y, res_z do k=0, res_z-1
integer, dimension(3) :: res, init, oppo, me, rear do j=0, res_y-1
real*8, dimension(3) :: geomdimension do i=0, res_x-1
real*8, dimension(3) :: myStep parameter_coords = (2.0*(/i+0.0,j+0.0,k+0.0/)-real(res)+fones)/(real(res)-fones)
real*8, dimension(3,3) :: defGrad_av positive = fones + parameter_coords
integer, dimension(3,8) :: corner = reshape((/ & negative = fones - parameter_coords
0, 0, 0,& coord_avgCorner(i+1,j+1,k+1,:) = ( coord_avgOrder(1,i+1,j+1,k+1,:) *negative(1)*negative(2)*negative(3)&
1, 0, 0,& + coord_avgOrder(2,i+1,j+1,k+1,:) *positive(1)*negative(2)*negative(3)&
1, 1, 0,& + coord_avgOrder(3,i+1,j+1,k+1,:) *positive(1)*positive(2)*negative(3)&
0, 1, 0,& + coord_avgOrder(4,i+1,j+1,k+1,:) *negative(1)*positive(2)*negative(3)&
1, 1, 1,& + coord_avgOrder(5,i+1,j+1,k+1,:) *positive(1)*positive(2)*positive(3)&
0, 1, 1,& + coord_avgOrder(6,i+1,j+1,k+1,:) *negative(1)*positive(2)*positive(3)&
0, 0, 1,& + coord_avgOrder(7,i+1,j+1,k+1,:) *negative(1)*negative(2)*positive(3)&
1, 0, 1 & + coord_avgOrder(8,i+1,j+1,k+1,:) *positive(1)*negative(2)*positive(3))*0.125
/), & enddo; enddo; enddo
(/3,8/)) end subroutine deformed
integer, dimension(3,8) :: step = reshape((/ &
1, 1, 1,&
-1, 1, 1,&
-1,-1, 1,&
1,-1, 1,&
-1,-1,-1,&
1,-1,-1,&
1, 1,-1,&
-1, 1,-1 &
/), &
(/3,8/))
integer, dimension(3,6) :: order = reshape((/ &
1, 2, 3,&
1, 3, 2,&
2, 1, 3,&
2, 3, 1,&
3, 1, 2,&
3, 2, 1 &
/), &
(/3,6/))
real*8 defGrad(3,3,res_x,res_y,res_z)
real*8 coordinates(3,res_x,res_y,res_z)
real*8, dimension(3,res_x,res_y,res_z,8) :: cornerCoords
real*8, dimension(3,2+res_x,2+res_y,2+res_z,6,8) :: coord
!f2py intent(in) res_x, res_y, res_z
!f2py intent(in) geomdimension
!f2py intent(out) coordinates
!f2py intent(in) defgrad
!f2py depend(res_x, res_y, res_z) coordinates
!f2py depend(res_x, res_y, res_z) defgrad
!f2py depend(res_x, res_y, res_z) cornerCoords
!f2py depend(res_x, res_y, res_z) coord
res = (/res_x,res_y,res_z/)
do i=1,3; do j=1,3
defgrad_av(i,j) = sum(defgrad(i,j,:,:,:)) /real(res(1)*res(2)*res(3))
enddo; enddo
print*, 'defgra', defgrad
do s = 1,8
init = corner(:,s)*(res-(/1,1,1/))
oppo = corner(:,1+mod(s-1+4,8))*(res-(/1,1,1/))
do o = 1,6
do k = init(order(3,o)),oppo(order(3,o)),step(order(3,o),s)
rear(order(2,o)) = init(order(2,o))
do j = init(order(2,o)),oppo(order(2,o)),step(order(2,o),s)
rear(order(1,o)) = init(order(1,o))
do i = init(order(1,o)),oppo(order(1,o)),step(order(1,o),s)
me(order(1,o)) = i
me(order(2,o)) = j
me(order(3,o)) = k
if (all(me == init)) then
coord(:,1+me(1),1+me(2),1+me(3),o,s) = &
geomdimension*(matmul(defgrad_av,real(corner(:,s),pDouble)) + &
matmul(defGrad(:,:,1+me(1),1+me(2),1+me(3)),step(:,s)/res/2.0_pDouble))
else
myStep = (me-rear)*geomdimension/res
coord(:,1+me(1),1+me(2),1+me(3),o,s) = coord(:,1+rear(1),1+rear(2),1+rear(3),o,s) + &
0.5_pDouble*matmul(defGrad(:,:,1+me(1),1+me(2),1+me(3)) + &
defGrad(:,:,1+rear(1),1+rear(2),1+rear(3)),&
myStep)
endif
rear = me
enddo
enddo
enddo
enddo ! orders
cornerCoords(:,:,:,:,s) = sum(coord(:,:,:,:,:,s),5)/6.0_pDouble
enddo ! corners
coordinates = sum(cornerCoords(:,:,:,:,:),5)/8.0_pDouble ! plain average no shape functions...
end subroutine
! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
subroutine mesh(res_x,res_y,res_z,geomdim,defgrad_av,centroids,nodes)
! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
implicit none
integer, parameter :: pDouble = selected_real_kind(15,50)
integer i,j,k, n
integer res_x, res_y, res_z
integer, dimension(3) :: res, shift, lookup, me
real*8, dimension(3) :: geomdim, diag = (/1,1,1/)
real*8, dimension(3,3) :: defgrad_av
real*8, dimension(3,res_x, res_y, res_z) :: centroids
real*8, dimension(3,res_x+2, res_y+2, res_z+2) :: wrappedCentroids
real*8, dimension(3,res_x+1, res_y+1, res_z+1) :: nodes
integer, dimension(3,8) :: neighbor = reshape((/ &
0, 0, 0,&
1, 0, 0,&
1, 1, 0,&
0, 1, 0,&
0, 0, 1,&
1, 0, 1,&
1, 1, 1,&
0, 1, 1 &
/), &
(/3,8/))
!f2py intent(in) res_x, res_y, res_z
!f2py intent(in) centroids
!f2py intent(in) defgrad_av
!f2py intent(in) geomdim
!f2py intent(out) nodes
!f2py depend(res_x, res_y, res_z) centroids
!f2py depend(res_x, res_y, res_z) nodes
res = (/res_x,res_y,res_z/)
wrappedCentroids(:,2:res(1)+1,2:res(2)+1,2:res(3)+1) = centroids
do k = 0,res(3)+1
do j = 0,res(2)+1
do i = 0,res(1)+1
if (&
k==0 .or. k==res(3)+1 .or. &
j==0 .or. j==res(2)+1 .or. &
i==0 .or. i==res(1)+1 &
) then
me = (/i,j,k/)
shift = (res+diag-2*me)/(res+diag)
lookup = me+shift*res
wrappedCentroids(:,1+i,1+j,1+k) = centroids(:,lookup(1),lookup(2),lookup(3)) - &
matmul(defgrad_av,shift * geomdim)
endif
enddo
enddo
enddo
nodes = 0.0_pDouble
do k = 0,res(3)
do j = 0,res(2)
do i = 0,res(1)
do n = 1,8
nodes(:,1+i,1+j,1+k) = nodes(:,1+i,1+j,1+k) + wrappedCentroids(:,1+i+neighbor(1,n),&
1+j+neighbor(2,n),&
1+k+neighbor(3,n))
enddo
nodes(:,1+i,1+j,1+k) = nodes(:,1+i,1+j,1+k) / 8.0_pDouble
enddo
enddo
enddo
end subroutine mesh
!below some code I used for gmsh postprocessing. Might be helpful
!!gmsh output
! character(len=1024) :: nriter
! character(len=1024) :: nrstep
! character(len=1024) :: nrloadcase
! real(pReal), dimension(:,:,:,:), allocatable :: displacement
! real(pReal), dimension(3,3) :: temp33_Real2
! real(pReal), dimension(3,3,3) :: Eigenvectorbasis
! real(pReal), dimension(3) :: Eigenvalue
! real(pReal) determinant
!!gmsh output
!!Postprocessing (gsmh output)
! do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1)
! if((k==1).and.(j==1).and.(i==1)) then
! temp33_Real =0.0_pReal
! else
! if((j==1).and.(i==1)) then
! temp33_Real(1,:) = temp33_Real(1,:) + math_mul33x3(defgrad(i,j,k,:,:),&
! (/0.0_pReal,0.0_pReal,(real(resolution(3))/meshdimension(3))/))
! temp33_Real(2,:) = temp33_Real(1,:)
! temp33_Real(3,:) = temp33_Real(1,:)
! displacement(i,j,k,:) = temp33_Real(1,:)
! else
! if(i==1) then
! temp33_Real(2,:) = temp33_Real(2,:) + math_mul33x3(defgrad(i,j,k,:,:),&
! (/0.0_pReal,(real(resolution(2))/meshdimension(2)),0.0_pReal/))
! temp33_Real(3,:) = temp33_Real(2,:)
! displacement(i,j,k,:) = temp33_Real(2,:)
! else
! temp33_Real(3,:) = temp33_Real(3,:) + math_mul33x3(defgrad(i,j,k,:,:),&
! (/(real(resolution(1))/meshdimension(1)),0.0_pReal,0.0_pReal/))
! displacement(i,j,k,:) = temp33_Real(3,:)
! endif
! endif
! endif
! enddo; enddo; enddo
! write(nrloadcase, *) loadcase; write(nriter, *) iter; write(nrstep, *) steps
! open(589,file = 'stress' //trim(adjustl(nrloadcase))//'-'//trim(adjustl(nrstep))//'-'//trim(adjustl(nriter))//'_cpfem.msh')
! open(588,file = 'logstrain'//trim(adjustl(nrloadcase))//'-'//trim(adjustl(nrstep))//'-'//trim(adjustl(nriter))//'_cpfem.msh')
! write(589, '(4(A, /), I10)'), '$MeshFormat', '2.1 0 8', '$EndMeshFormat', '$Nodes', prodnn
! write(588, '(4(A, /), I10)'), '$MeshFormat', '2.1 0 8', '$EndMeshFormat', '$Nodes', prodnn
! ielem = 0_pInt
! do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1)
! ielem = ielem + 1
! write(589, '(I10, 3(tr2, E12.6))'), ielem, displacement(i,j,k,:) !for deformed configuration
! write(588, '(I10, 3(tr2, E12.6))'), ielem, displacement(i,j,k,:)
!! write(589, '(4(I10,tr2))'), ielem, i-1,j-1,k-1 !for undeformed configuration
!!write(588, '(4(I10,tr2))'), ielem, i-1,j-1,k-1
! enddo; enddo; enddo
! write(589, '(2(A, /), I10)'), '$EndNodes', '$Elements', prodnn
! write(588, '(2(A, /), I10)'), '$EndNodes', '$Elements', prodnn
! do i = 1, prodnn
! write(589, '(I10, A, I10)'), i, ' 15 2 1 2', i
! write(588, '(I10, A, I10)'), i, ' 15 2 1 2', i
! enddo
! write(589, '(A)'), '$EndElements'
! write(588, '(A)'), '$EndElements'
! write(589, '(8(A, /), I10)'), '$NodeData', '1','"'//trim(adjustl('stress'//trim(adjustl(nrloadcase))//'-'//&
! trim(adjustl(nrstep))//'-'//trim(adjustl(nriter))//'_cpfem.msh'))//'"','1','0.0', '3', '0', '9', prodnn
! write(588, '(8(A, /), I10)'), '$NodeData', '1','"'//trim(adjustl('logstrain'//trim(adjustl(nrloadcase))//'-'//&
! trim(adjustl(nrstep))//'-'//trim(adjustl(nriter))//'_cpfem.msh'))//'"','1','0.0', '3', '0', '9', prodnn
! ielem = 0_pInt
! do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1)
! ielem = ielem + 1
! write(589, '(i10, 9(tr2, E14.8))'), ielem, cstress_field(i,j,k,:,:)
! call math_pDecomposition(defgrad(i,j,k,:,:),temp33_Real2,temp33_Real,errmatinv) !store R in temp33_Real
! call math_invert3x3(temp33_Real, temp33_Real2, determinant, errmatinv) !inverse of R in temp33_Real2
! temp33_Real = math_mul33x33(defgrad(i,j,k,:,:), temp33_Real2) ! v = F o inv(R), store in temp33_Real
! call math_spectral1(temp33_Real,Eigenvalue(1), Eigenvalue(2), Eigenvalue(3),&
! Eigenvectorbasis(1,:,:),Eigenvectorbasis(2,:,:),Eigenvectorbasis(3,:,:))
! eigenvalue = log(sqrt(eigenvalue))
! temp33_Real = eigenvalue(1)*Eigenvectorbasis(1,:,:)+eigenvalue(2)*Eigenvectorbasis(2,:,:)+eigenvalue(3)*Eigenvectorbasis(3,:,:)
! write(588, '(i10, 9(tr2, E14.8))'), ielem, temp33_Real
! enddo; enddo; enddo
! write(589, *), '$EndNodeData'
! write(588, *), '$EndNodeData'
! close(589); close(588); close(540)
! enddo ! end looping over steps in current loadcase
! enddo ! end looping over loadcases
! close(539); close(538)
! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
function coordinates3(res_x,res_y,res_z,geomdimension,defgrad)
! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
implicit none
integer i,j,k, l,m, s,o, loop, res_x, res_y, res_z
real*8 defgrad_av(3,3)
real*8 geomdimension(3)
integer res(3)
real*8 defgrad(3,3,res_x,res_y,res_z)
real*8 coordinates3(3,res_x,res_y,res_z)
!f2py intent(in) res_x, res_y, res_z
!f2py depend(res_x, res_y, res_z) coordinates3
!f2py depend(res_x, res_y, res_z) defgrad
!f2py intent(in) geomdimension
!f2py intent(out) coordinates3
!f2py intent(in) defgrad
res = (/res_x,res_y,res_z/)
do i=1,3; do j=1,3
defgrad_av(i,j) = sum(defgrad(i,j,:,:,:)) /real(res(1)*res(2)*res(3))
enddo; enddo
print*, 'defgra', defgrad
coordinates3 = 0.0
end function

29
processing/post/spectral_post.py
View File

@ -6,7 +6,7 @@
# computed using the FEM solvers. Until now, its capable to handle elements with one IP in a regular order # computed using the FEM solvers. Until now, its capable to handle elements with one IP in a regular order
# written by M. Diehl, m.diehl@mpie.de # written by M. Diehl, m.diehl@mpie.de
import os,sys,re,array,struct,numpy, time import os,sys,re,array,struct,numpy, time, reconstruct
class vector: class vector:
x,y,z = [None,None,None] x,y,z = [None,None,None]
@ -63,8 +63,6 @@ class MPIEspectral_result:
self.N_increments = self._keyedInt('increments') self.N_increments = self._keyedInt('increments')
self.N_element_scalars = self._keyedInt('materialpoint_sizeResults') self.N_element_scalars = self._keyedInt('materialpoint_sizeResults')
self.resolution = self._keyedPackedArray('resolution',3,'i') self.resolution = self._keyedPackedArray('resolution',3,'i')
#print self.resolution
#self.resolution = numpy.array([10,10,10],'i')
self.N_nodes = (self.resolution[0]+1)*(self.resolution[1]+1)*(self.resolution[2]+1) self.N_nodes = (self.resolution[0]+1)*(self.resolution[1]+1)*(self.resolution[2]+1)
self.N_elements = self.resolution[0]*self.resolution[1]*self.resolution[2] self.N_elements = self.resolution[0]*self.resolution[1]*self.resolution[2]
self.dimension = self._keyedPackedArray('dimension',3,'d') self.dimension = self._keyedPackedArray('dimension',3,'d')
@ -95,7 +93,7 @@ class MPIEspectral_result:
def _keyedPackedArray(self,identifier,length = 3,type = 'd'): def _keyedPackedArray(self,identifier,length = 3,type = 'd'):
match = {'d': 8,'i': 4} match = {'d': 8,'i': 4}
self.file.seek(0) self.file.seek(0)
m = re.search('%s%s'%(identifier,'(.{%i})'%(match[type])*length),self.file.read(2048)) m = re.search('%s%s'%(identifier,'(.{%i})'%(match[type])*length),self.file.read(2048),re.DOTALL)
values = [] values = []
if m: if m:
for i in m.groups(): for i in m.groups():
@ -105,7 +103,7 @@ class MPIEspectral_result:
def _keyedInt(self,identifier): def _keyedInt(self,identifier):
value = None value = None
self.file.seek(0) self.file.seek(0)
m = re.search('%s%s'%(identifier,'(.{4})'),self.file.read(2048)) m = re.search('%s%s'%(identifier,'(.{4})'),self.file.read(2048),re.DOTALL)
if m: if m:
value = struct.unpack('i',m.group(1))[0] value = struct.unpack('i',m.group(1))[0]
return value return value
@ -113,7 +111,7 @@ class MPIEspectral_result:
def _keyedString(self,identifier): def _keyedString(self,identifier):
value = None value = None
self.file.seek(0) self.file.seek(0)
m = re.search(r'(.{4})%s(.*?)\1'%identifier,self.file.read(2048)) m = re.search(r'(.{4})%s(.*?)\1'%identifier,self.file.read(2048),re.DOTALL)
if m: if m:
value = m.group(2) value = m.group(2)
return value return value
@ -219,7 +217,7 @@ def mesh(res,geomdim,defgrad_av,centroids):
for i in range(res[0]+1): for i in range(res[0]+1):
for n in range(8): for n in range(8):
nodes[i,j,k] += wrappedCentroids[i+neighbor[n,0],j+neighbor[n,1],k+neighbor[n,2]] nodes[i,j,k] += wrappedCentroids[i+neighbor[n,0],j+neighbor[n,1],k+neighbor[n,2]]
nodes[:,:,:] /=8.0 nodes[:,:,:] /= 8.0
return nodes return nodes
@ -261,7 +259,7 @@ def centroids(res,geomdimension,defgrad):
ones = numpy.ones(3, 'i') ones = numpy.ones(3, 'i')
fones = numpy.ones(3, 'd') fones = numpy.ones(3, 'd')
defgrad_av=numpy.average(numpy.average(numpy.average(defgrad,0),0),0) defgrad_av=numpy.average(numpy.average(numpy.average(defgrad,0),0),0)
print defgrad_av
for s in range(8):# corners for s in range(8):# corners
init = corner[s]*(res-ones) init = corner[s]*(res-ones)
oppo = corner[(s+4)%8]*(res-ones) oppo = corner[(s+4)%8]*(res-ones)
@ -473,17 +471,18 @@ res_z=p.resolution[2]
# print(struct.unpack('d',p.file.read(8))) # print(struct.unpack('d',p.file.read(8)))
ms=numpy.zeros([res_x,res_y,res_z,3], 'd')
for i in range(40,46): for i in range(249,250):
c_pos = p.dataOffset + i*(p.N_element_scalars*8*p.N_elements + 8) #8 accounts for header&footer c_pos = p.dataOffset + i*(p.N_element_scalars*8*p.N_elements + 8) #8 accounts for header&footer
defgrad = readTensor(p.resolution,p.file,p.N_element_scalars,c_pos,16) defgrad = readTensor(p.resolution,p.file,p.N_element_scalars,c_pos,16)
p_stress = readTensor(p.resolution,p.file,p.N_element_scalars,c_pos,52) p_stress = readTensor(p.resolution,p.file,p.N_element_scalars,c_pos,58)
#c_stress = calculateCauchyStress(p_stress,defgrad,p.resolution) #c_stress = calculateCauchyStress(p_stress,defgrad,p.resolution)
#grain = calculateVonMises(c_stress,p.resolution) #grain = calculateVonMises(c_stress,p.resolution)
centroids_coord, defgrad_av = centroids(p.resolution,p.dimension,defgrad) defgrad_av=numpy.average(numpy.average(numpy.average(defgrad,0),0),0)
ms = mesh(p.resolution,p.dimension,defgrad_av,centroids_coord) centroids_coord = reconstruct.deformed(res_x,res_y,res_z,p.dimension,defgrad,defgrad_av)
writeVtkAscii(name+'-mesh-%i.vtk'%i,ms,p_stress[:,:,:,1,2],p.resolution) ms = reconstruct.mesh(p.resolution[0],p.resolution[1],p.resolution[2],p.dimension,defgrad_av,centroids_coord)
writeVtkAsciidefgrad_av(name+'-box-%i.vtk'%i,p.dimension,defgrad_av) writeVtkAscii(name+'-mesh-fortran-%s.vtk'%i,ms,p_stress[:,:,:,1,2],p.resolution)
#writeVtkAsciidefgrad_av(name+'-box-%i.vtk'%i,p.dimension,defgrad_av)
#writeVtkAsciiDots(name+'-points-%i.vtk'%i,centroids_coord,grain,p.resolution) #writeVtkAsciiDots(name+'-points-%i.vtk'%i,centroids_coord,grain,p.resolution)
sys.stdout.flush() sys.stdout.flush()