diff --git a/code/mpie_spectral.f90 b/code/mpie_spectral.f90 index 1b11e6cec..c7349d86f 100644 --- a/code/mpie_spectral.f90 +++ b/code/mpie_spectral.f90 @@ -5,10 +5,10 @@ ! written by P. Eisenlohr, ! F. Roters, ! L. Hantcherli, -! W.A. Counts -! D.D. Tjahjanto -! C. Kords -! M. Diehl +! W.A. Counts, +! D.D. Tjahjanto, +! C. Kords, +! M. Diehl, ! R. Lebensohn ! ! MPI fuer Eisenforschung, Duesseldorf @@ -17,9 +17,9 @@ ! Usage: ! - start program with mpie_spectral PathToGeomFile/NameOfGeom.geom ! 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 -! directory +! directory. For further configuration use "numerics.config" !******************************************************************** program mpie_spectral !******************************************************************** @@ -29,7 +29,8 @@ program mpie_spectral use IO use math 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 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,3) :: dPdF, c0, s0 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 :: pstress_field, defgrad, defgradold @@ -79,27 +80,23 @@ program mpie_spectral complex(pReal), dimension(3,3) :: temp33_Complex real(pReal), dimension(3,3) :: xinormdyad real(pReal), dimension(:,:,:,:,:,:,:), allocatable :: gamma_hat - real(pReal), dimension(:,:,:,:), allocatable :: xi integer(pInt), dimension(3) :: k_s + real(pReal), dimension(3) :: xi, xi_middle integer*8, dimension(2,3,3) :: plan_fft -! convergence etc. - real(pReal) 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 +! loop variables, convergence etc. + real(pReal) guessmode, err_div, err_stress, err_defgrad, sigma0 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 + !Initializing !$ call omp_set_num_threads(mpieNumThreadsInt) ! set number of threads for parallel execution set by MPIE_NUM_THREADS bc_maskvector = '' unit = 234_pInt - ones = 1.0_pReal; zeroes = 0.0_pReal N_l = 0_pInt; N_s = 0_pInt @@ -184,10 +181,10 @@ program mpie_spectral do i = 1, N_Loadcases 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(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 velocitygrad',transpose(bc_mask(:,:,1,i)) + print '(a,/,3(3(l,x)/))', 'bc_mask for stress' ,transpose(bc_mask(:,:,2,i)) print *,'time',bc_timeIncrement(i) print *,'incs',bc_steps(i) print *, '' @@ -241,33 +238,15 @@ program mpie_spectral print '(a,/,f6.1,f6.1,f6.1)','dimension x y z', geomdimension print *,'homogenization',homog - allocate (workfft(resolution(1)/2+1,resolution(2),resolution(3),3,3)); workfft = 0.0_pReal - allocate (gamma_hat(resolution(1)/2+1,resolution(2),resolution(3),3,3,3,3)); gamma_hat = 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) - ! + 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 + prodnn = resolution(1)*resolution(2)*resolution(3) - wgt = 1_pReal/real(prodnn, pReal) + wgt = 1.0_pReal/real(prodnn, pReal) defgradAim = math_I3 defgradAimOld = math_I3 defgrad_av = math_I3 + ! Initialization of CPFEM_general (= constitutive law) and of deformation gradient field call CPFEM_initAll(temperature,1_pInt,1_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) c066 = c066 + dsde enddo; enddo; enddo - c066 = c066 * wgt - c0 = math_mandel66to3333(c066) - call math_invert(6, c066, s066,i, errmatinv) - s0 = math_mandel66to3333(s066) + c066 = c066 * wgt + c0 = math_mandel66to3333(c066) + call math_invert(6, c066, s066,i, errmatinv) + s0 = math_mandel66to3333(s066) -!calculation of xinormdyad (to calculate gamma_hat) and xi (waves, for proof of equilibrium) - 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(i,j,k,3) = 0.0_pReal - if(resolution(3) > 1) xi(i,j,k,3) = real(k_s(3), pReal)/geomdimension(3) - xi(i,j,k,2) = real(k_s(2), pReal)/geomdimension(2) - xi(i,j,k,1) = real(k_s(1), pReal)/geomdimension(1) - if (any(xi(i,j,k,:) /= 0.0_pReal)) then - do l = 1,3; do m = 1,3 - xinormdyad(l,m) = xi(i,j,k, l)*xi(i,j,k, m)/sum(xi(i,j,k,:)**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(i,j,k, 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 - enddo; enddo; enddo +!calculation of calculate gamma_hat field in the case of fast execution (needs a lot of memory) + if(fast_execution) then + allocate (gamma_hat(resolution(1)/2+1,resolution(2),resolution(3),3,3,3,3)); gamma_hat = 0.0_pReal + 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(i,j,k, 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 + 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 - 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), 'workingdir',trim(getSolverWorkingDirectoryName()) write(538), 'geometry',trim(getSolverJobName())//InputFileExtension @@ -373,7 +380,8 @@ program mpie_spectral (err_div > err_div_tol .or. & err_stress > err_stress_tol .or. & 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 cstress_av = 0.0_pReal !************************************************************* @@ -407,7 +415,7 @@ program mpie_spectral pstress_av(m,n) = sum(pstress_field(:,:,:,m,n)) * wgt defgrad_av(m,n) = sum(defgrad(:,:,:,m,n)) * wgt enddo; enddo - + err_stress = maxval(abs(mask_stress * (pstress_av - bc_stress(:,:,loadcase)))) err_stress_tol = maxval(abs(pstress_av))*err_stress_tolrel @@ -431,19 +439,19 @@ program mpie_spectral err_div = 2 * err_div_tol 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(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 '(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 deformation gradient ',err_defgrad,' Tol. = ', err_defgrad_tol*0.8 if(err_stress < err_stress_tol*0.8) then calcmode = 1_pInt endif ! 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))' ielem = 0_pInt 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,:,:),& temperature,timeinc,ielem,1_pInt,& cstress,dsde, pstress, dPdF) @@ -451,7 +459,7 @@ program mpie_spectral ielem = 0_pInt 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,& defgradold(i,j,k,:,:), defgrad(i,j,k,:,:),& temperature,timeinc,ielem,1_pInt,& @@ -459,46 +467,83 @@ program mpie_spectral pstress_field(i,j,k,:,:) = pstress cstress_av = cstress_av + math_mandel6to33(cstress) 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' 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,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 - enddo; enddo - err_div = (maxval(abs(math_mul33x3_complex(workfft(resolution(1)/2+1,resolution(2)/2+1,resolution(3)/2+1,:,:),& - xi(resolution(1)/2+1,resolution(2)/2+1,resolution(3)/2+1,:))))) ! L infinity Norm of div(stress) - + call dfftw_execute_dft_r2c(plan_fft(1,1,1), pstress_field,workfft) ! FFT of pstress + + do m = 1,3 + sigma0 = max(sigma0, sum(abs(workfft(1,1,1,m,:)))) ! L infinity Norm of stress tensor + enddo + 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 temp33_Complex = 0.0_pReal do m = 1,3; do n = 1,3 temp33_Complex(m,n) = sum(gamma_hat(i,j,k,m,n,:,:) * workfft(i,j,k,:,:)) enddo; enddo workfft(i,j,k,:,:) = temp33_Complex(:,:) - enddo; enddo; enddo - workfft(1,1,1,:,:) = defgrad_av - math_I3 + enddo; enddo; enddo + + 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 - - cstress_av = cstress_av * wgt + workfft(1,1,1,:,:) = defgrad_av - math_I3 !zero frequency do m = 1,3; do n = 1,3 - call dfftw_execute_dft_c2r(plan_fft(2,m,n), workfft(:,:,:,m,n),ddefgrad(:,:,:)) - 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(:,:,:,m,n) = defgrad(:,:,:,m,n) + (defgradAim(m,n) - defgrad_av(m,n)) !anticipated target minus current state + call dfftw_execute_dft_c2r(plan_fft(2,m,n), workfft(:,:,:,m,n),ddefgrad(:,:,:)) + defgrad(:,:,:,m,n) = defgrad(:,:,:,m,n) + ddefgrad * 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 enddo; enddo - 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_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 stress ',err_stress,' Tol. = ', err_stress_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 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. - calcmode = 0 + calcmode = 0_pInt defgradAimCorr = 0.0_pReal damper = damper * 0.9_pReal endif @@ -508,9 +553,9 @@ program mpie_spectral write(538) materialpoint_results(:,1,:) !write to output file 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 Gradient: ',defgrad_av(1:3,:) - print '(a,/,3(3(f10.4,x)/))', ' Cauchy Stress [MPa]: ',cstress_av(1:3,:)/1.e6 + print '(a,/,3(3(f12.7,x)/))', ' Deformation Aim: ',math_transpose3x3(defgradAim) + 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)', '************************************************************' enddo ! end looping over steps in current loadcase enddo ! end looping over loadcases diff --git a/code/mpie_spectral_interface.f90 b/code/mpie_spectral_interface.f90 index 5f83ec8c0..4a4c57801 100644 --- a/code/mpie_spectral_interface.f90 +++ b/code/mpie_spectral_interface.f90 @@ -34,7 +34,7 @@ function getSolverWorkingDirectoryName() character(len=1024) cwd,outname,getSolverWorkingDirectoryName 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 getSolverWorkingDirectoryName = outname(1:scan(outname,pathSep,back=.true.)) @@ -84,6 +84,33 @@ function getSolverJobName() return 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 diff --git a/code/numerics.f90 b/code/numerics.f90 index 74b542ffa..0c43dacab 100644 --- a/code/numerics.f90 +++ b/code/numerics.f90 @@ -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_tolrel, & ! factor to multiply with highest stress to get err_stress_tol 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 + !* Random seeding parameters fixedSeed ! fixed seeding for pseudo-random number generator !* OpenMP variable @@ -144,6 +146,7 @@ subroutine numerics_init() err_defgrad_tol = 1.0e-3 err_stress_tolrel = 0.01 itmax = 20_pInt + fast_execution = .false. !* Random seeding parameters: added <<>> fixedSeed = 0_pInt @@ -253,6 +256,8 @@ subroutine numerics_init() err_stress_tolrel = IO_floatValue(line,positions,2) case ('itmax') itmax = IO_intValue(line,positions,2) + case ('fast_execution') + fast_execution = IO_intValue(line,positions,2) > 0_pInt !* Random seeding parameters 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_stress_tolrel: ',err_stress_tolrel write(6,'(a24,x,i8)') 'itmax: ',itmax + write(6,'(a24,x,L8)') 'fast_execution: ',fast_execution write(6,*) !* Random seeding parameters diff --git a/processing/post/3Dvisualize b/processing/post/3Dvisualize index 8b205423a..36baf5119 100755 --- a/processing/post/3Dvisualize +++ b/processing/post/3Dvisualize @@ -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])], ] +# 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') # output(cmd,{'filepointer':vtk},'File') # vtk.close() diff --git a/processing/post/make_reconstruct.py b/processing/post/make_reconstruct.py index 3cbafd3f7..fbbdc808b 100644 --- a/processing/post/make_reconstruct.py +++ b/processing/post/make_reconstruct.py @@ -5,4 +5,11 @@ # module reconstruct.so out of the fortran code reconstruct.f90 # 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 \ No newline at end of file diff --git a/processing/post/reconstruct.f90 b/processing/post/reconstruct.f90 index 8b116858b..c1d825bdc 100644 --- a/processing/post/reconstruct.f90 +++ b/processing/post/reconstruct.f90 @@ -1,421 +1,157 @@ -! -*- f90 -*- -function coordinates2(res_x,res_y,res_z,geomdimension,defgrad) - implicit none - integer, parameter :: pDouble = selected_real_kind(15,50) - integer i,j,k, l,m, s,o, loop, res_x, res_y, res_z - integer, dimension(3) :: res, init, oppo, me, rear - real*8, dimension(3,3) :: defgrad_av - integer, dimension(3) :: resolution - real*8, dimension(3) :: geomdimension, myStep - real*8, dimension(3,3) :: temp33_Real - real*8, dimension(3,res_x, res_y, res_z) :: coordinates2 - real*8, dimension(3,3,res_x, res_y, res_z) :: defgrad - 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 - !f2py intent(in) res_y - !f2py intent(in) res_z - !f2py intent(in) geomdimension - !f2py intent(in) defgrad - !f2py intent(out) coordinates2 - !f2py depend(res_x, res_y, res_z) coordinates2 - !f2py depend(res_x, res_y, res_z) defgrad - - ! integer, dimension(3,8) :: corner = reshape((/ & - ! 0, 0, 0,& - ! 1, 0, 0,& - ! 1, 1, 0,& - ! 0, 1, 0,& - ! 1, 1, 1,& - ! 0, 1, 1,& - ! 0, 0, 1,& - ! 1, 0, 1 & - ! /), & - ! (/3,8/)) - - ! 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/)) - - resolution = (/res_x,res_y,res_z/) - - write(6,*) 'defgrad', defgrad - - do i=1, 3; do j=1,3 - defgrad_av(i,j) = sum(defgrad(i,j,:,:,:)) /real(res_x*res_y*res_z) - enddo; enddo - - 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 = real(0.0) - else - if((j==1).and.(i==1)) then - temp33_Real(1,:) = temp33_Real(1,:) + matmul(defgrad(:,:,i,j,k),& - (/real(0.0),real(0.0),real(geomdimension(3))/(real(resolution(3)))/)) - temp33_Real(2,:) = temp33_Real(1,:) - temp33_Real(3,:) = temp33_Real(1,:) - coordinates2(:,i,j,k) = temp33_Real(1,:) - else - if(i==1) then - temp33_Real(2,:) = temp33_Real(2,:) + matmul(defgrad(:,:,i,j,k),& - (/real(0.0),real(geomdimension(2))/(real(resolution(2))),real(0.0)/)) - temp33_Real(3,:) = temp33_Real(2,:) - coordinates2(:,i,j,k) = temp33_Real(2,:) - else - temp33_Real(3,:) = temp33_Real(3,:) + matmul(defgrad(:,:,i,j,k),& - (/real(geomdimension(1))/(real(resolution(1))),real(0.0),real(0.0)/)) - coordinates2(:,i,j,k) = temp33_Real(3,:) - endif - endif - endif - enddo; enddo; enddo - do i=1, res_x; do j = 1, res_y; do k = 1, res_z - 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)/)) - enddo; enddo; enddo - - 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 - - ! 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) - ! ! print*, order(:,o) - ! me(order(1,o)) = i - ! me(order(2,o)) = j - ! me(order(3,o)) = k - ! ! print*, me -! ! if (all(me == init)) then -! ! coord(:,1+me(1),1+me(2),1+me(3),o,s) = 0.0 !& - ! ! 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 function - - - -! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -subroutine coordinates5(res_x,res_y,res_z,geomdimension,defgrad) -! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - implicit none - integer, parameter :: pDouble = selected_real_kind(15,50) - integer i,j,k, l,m, s,o, loop, res_x, res_y, res_z - integer, dimension(3) :: res, init, oppo, me, rear - real*8, dimension(3) :: geomdimension - real*8, dimension(3) :: myStep - real*8, dimension(3,3) :: defGrad_av - integer, dimension(3,8) :: corner = reshape((/ & - 0, 0, 0,& - 1, 0, 0,& - 1, 1, 0,& - 0, 1, 0,& - 1, 1, 1,& - 0, 1, 1,& - 0, 0, 1,& - 1, 0, 1 & - /), & - (/3,8/)) - - 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 \ No newline at end of file +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +subroutine mesh(res_x,res_y,res_z,geomdim,defgrad_av,centroids,nodes) +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + implicit none + real*8 geomdim(3) + integer res_x, res_y, res_z + real*8 wrappedCentroids(res_x+2,res_y+2,res_z+2,3) + real*8 nodes(res_x+1,res_y+1,res_z+1,3) + real*8 centroids(res_x ,res_y ,res_z ,3) + + 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/)) + + integer i,j,k,n + real*8, dimension(3,3) :: defgrad_av + integer, dimension(3) :: diag, shift, lookup, me, res + + diag = 1 + shift = 0 + lookup = 0 + + res = (/res_x,res_y,res_z/) + + wrappedCentroids=0.0 + wrappedCentroids(2:res_x+1,2:res_y+1,2:res_z+1,:) = centroids + + do k=0, res_z+1 + do j=0, res_y+1 + do i=0, res_x+1 + if (k==0 .or. k==res_z+1 .or. & + j==0 .or. j==res_y+1 .or. & + i==0 .or. i==res_x+1 ) then + me = (/i,j,k/) + shift = sign(abs(res+diag-2*me)/(res+diag),res+diag-2*me) + lookup = me-diag+shift*res + wrappedCentroids(i+1,j+1,k+1,:) = centroids(lookup(1)+1,lookup(2)+1,lookup(3)+1,:)- & + matmul(defgrad_av, shift*geomdim) + endif + enddo; enddo; enddo + do k=0, res_z + do j=0, res_y + do i=0, res_x + do n=1,8 + 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),:) + enddo; enddo; enddo; enddo + nodes = nodes/8.0 + + end subroutine mesh + +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +subroutine deformed(res_x,res_y,res_z,geomdim,defgrad,defgrad_av,coord_avgCorner) +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + implicit none + real*8 geomdim(3) + integer res_x, res_y, res_z + real*8 coord(8,6,res_x,res_y,res_z,3) + real*8 coord_avgOrder(8,res_x,res_y,res_z,3) + real*8 coord_avgCorner(res_x,res_y,res_z,3) + real*8 defgrad(res_x,res_y,res_z,3,3) + integer, dimension(3,8) :: corner = reshape((/ & + 0, 0, 0,& + 1, 0, 0,& + 1, 1, 0,& + 0, 1, 0,& + 1, 1, 1,& + 0, 1, 1,& + 0, 0, 1,& + 1, 0, 1 & + /), & + (/3,8/)) + 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 myStep(3), fones(3), parameter_coords(3) + real*8 defgrad_av(3,3) + real*8 negative(3), positive(3) + integer rear(3), init(3), ones(3), oppo(3), me(3), res(3) + integer i, j, k, s, o + + ones = 1 + fones = 1.0 + coord_avgOrder=0.0 + + res = (/res_x,res_y,res_z/) + + do s = 0, 7 ! corners (from 0 to 7) + init = corner(:,s+1)*(res-ones) +ones + oppo = corner(:,mod((s+4),8)+1)*(res-ones) +ones + do o=1,6 ! orders ! from 1 to 6) + do k = init(order(3,o)), oppo(order(3,o)), step(order(3,o),s+1) + rear(order(2,o)) = init(order(2,o)) + do j = init(order(2,o)), oppo(order(2,o)), step(order(2,o),s+1) + rear(order(1,o)) = init(order(1,o)) + do i = init(order(1,o)), oppo(order(1,o)), step(order(1,o),s+1) + me(order(1,o)) = i + me(order(2,o)) = j + me(order(3,o)) = k + if ( (me(1)==init(1)).and.(me(2)==init(2)).and. (me(3)==init(3)) ) then + coord(s+1,o,me(1),me(2),me(3),:) = geomdim * (matmul(defgrad_av,corner(:,s+1)) + & + matmul(defgrad(me(1),me(2),me(3),:,:),0.5*step(:,s+1)/res)) + + else + 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),:) + & + 0.5*matmul(defgrad(me(1),me(2),me(3),:,:) + & + defgrad(rear(1),rear(2),rear(3),:,:),myStep) + endif + rear = me + enddo; enddo; enddo; enddo + do i=1,6 + coord_avgOrder(s+1,:,:,:,:) = coord_avgOrder(s+1,:,:,:,:) + coord(s+1,i,:,:,:,:)/6.0 + enddo + enddo + + do k=0, res_z-1 + do j=0, res_y-1 + do i=0, res_x-1 + parameter_coords = (2.0*(/i+0.0,j+0.0,k+0.0/)-real(res)+fones)/(real(res)-fones) + positive = fones + parameter_coords + negative = fones - parameter_coords + coord_avgCorner(i+1,j+1,k+1,:) = ( coord_avgOrder(1,i+1,j+1,k+1,:) *negative(1)*negative(2)*negative(3)& + + coord_avgOrder(2,i+1,j+1,k+1,:) *positive(1)*negative(2)*negative(3)& + + coord_avgOrder(3,i+1,j+1,k+1,:) *positive(1)*positive(2)*negative(3)& + + coord_avgOrder(4,i+1,j+1,k+1,:) *negative(1)*positive(2)*negative(3)& + + coord_avgOrder(5,i+1,j+1,k+1,:) *positive(1)*positive(2)*positive(3)& + + coord_avgOrder(6,i+1,j+1,k+1,:) *negative(1)*positive(2)*positive(3)& + + coord_avgOrder(7,i+1,j+1,k+1,:) *negative(1)*negative(2)*positive(3)& + + coord_avgOrder(8,i+1,j+1,k+1,:) *positive(1)*negative(2)*positive(3))*0.125 + enddo; enddo; enddo +end subroutine deformed \ No newline at end of file diff --git a/processing/post/spectral_post.py b/processing/post/spectral_post.py index 4372c6bee..add25f670 100644 --- a/processing/post/spectral_post.py +++ b/processing/post/spectral_post.py @@ -6,7 +6,7 @@ # 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 -import os,sys,re,array,struct,numpy, time +import os,sys,re,array,struct,numpy, time, reconstruct class vector: x,y,z = [None,None,None] @@ -63,8 +63,6 @@ class MPIEspectral_result: self.N_increments = self._keyedInt('increments') self.N_element_scalars = self._keyedInt('materialpoint_sizeResults') 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_elements = self.resolution[0]*self.resolution[1]*self.resolution[2] self.dimension = self._keyedPackedArray('dimension',3,'d') @@ -95,7 +93,7 @@ class MPIEspectral_result: def _keyedPackedArray(self,identifier,length = 3,type = 'd'): match = {'d': 8,'i': 4} 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 = [] if m: for i in m.groups(): @@ -105,7 +103,7 @@ class MPIEspectral_result: def _keyedInt(self,identifier): value = None 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: value = struct.unpack('i',m.group(1))[0] return value @@ -113,7 +111,7 @@ class MPIEspectral_result: def _keyedString(self,identifier): value = None 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: value = m.group(2) return value @@ -219,7 +217,7 @@ def mesh(res,geomdim,defgrad_av,centroids): for i in range(res[0]+1): for n in range(8): 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 @@ -261,7 +259,7 @@ def centroids(res,geomdimension,defgrad): ones = numpy.ones(3, 'i') fones = numpy.ones(3, 'd') defgrad_av=numpy.average(numpy.average(numpy.average(defgrad,0),0),0) - + print defgrad_av for s in range(8):# corners init = corner[s]*(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))) - -for i in range(40,46): +ms=numpy.zeros([res_x,res_y,res_z,3], 'd') +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 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) #grain = calculateVonMises(c_stress,p.resolution) - centroids_coord, defgrad_av = centroids(p.resolution,p.dimension,defgrad) - ms = mesh(p.resolution,p.dimension,defgrad_av,centroids_coord) - writeVtkAscii(name+'-mesh-%i.vtk'%i,ms,p_stress[:,:,:,1,2],p.resolution) - writeVtkAsciidefgrad_av(name+'-box-%i.vtk'%i,p.dimension,defgrad_av) + defgrad_av=numpy.average(numpy.average(numpy.average(defgrad,0),0),0) + centroids_coord = reconstruct.deformed(res_x,res_y,res_z,p.dimension,defgrad,defgrad_av) + ms = reconstruct.mesh(p.resolution[0],p.resolution[1],p.resolution[2],p.dimension,defgrad_av,centroids_coord) + 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) sys.stdout.flush()