diff --git a/code/mpie_spectral.f90 b/code/mpie_spectral.f90 index c7349d86f..087f79b58 100644 --- a/code/mpie_spectral.f90 +++ b/code/mpie_spectral.f90 @@ -60,30 +60,28 @@ program mpie_spectral ! variables storing information from geom file real(pReal) wgt real(pReal), dimension(3) :: geomdimension - integer(pInt) homog, prodnn + integer(pInt) homog integer(pInt), dimension(3) :: resolution ! stress etc. real(pReal), dimension(3,3) :: ones, zeroes, temp33_Real, damper,& pstress, pstress_av, cstress_av, defgrad_av,& defgradAim, defgradAimOld, defgradAimCorr, defgradAimCorrPrev,& - mask_stress, mask_defgrad - real(pReal), dimension(3,3,3) :: temp333_Real + mask_stress, mask_defgrad 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 ! Mandel notation of 4th order tensors - real(pReal), dimension(:,:,:), allocatable :: ddefgrad - real(pReal), dimension(:,:,:,:,:), allocatable :: pstress_field, defgrad, defgradold + real(pReal), dimension(:,:,:,:,:), allocatable :: workfft, defgrad, defgradold ! variables storing information for spectral method - complex(pReal), dimension(:,:,:,:,:), allocatable :: workfft + complex(pReal) :: img complex(pReal), dimension(3,3) :: temp33_Complex real(pReal), dimension(3,3) :: xinormdyad real(pReal), dimension(:,:,:,:,:,:,:), allocatable :: gamma_hat - integer(pInt), dimension(3) :: k_s real(pReal), dimension(3) :: xi, xi_middle - integer*8, dimension(2,3,3) :: plan_fft - + integer(pInt), dimension(3) :: k_s + integer*8, dimension(2) :: plan_fft + ! loop variables, convergence etc. real(pReal) guessmode, err_div, err_stress, err_defgrad, sigma0 integer(pInt) i, j, k, l, m, n, p @@ -98,6 +96,7 @@ program mpie_spectral bc_maskvector = '' unit = 234_pInt ones = 1.0_pReal; zeroes = 0.0_pReal + img = cmplx(0.0,1.0) N_l = 0_pInt; N_s = 0_pInt N_t = 0_pInt; N_n = 0_pInt @@ -114,7 +113,7 @@ program mpie_spectral print*,'Workingdir: ',trim(getSolverWorkingDirectoryName()) if (.not. IO_open_file(unit,path)) call IO_error(45,ext_msg = path) - + rewind(unit) do read(unit,'(a1024)',END = 101) line @@ -238,11 +237,10 @@ program mpie_spectral print '(a,/,f6.1,f6.1,f6.1)','dimension x y z', geomdimension print *,'homogenization',homog - allocate (defgrad(resolution(1),resolution(2),resolution(3),3,3)); defgrad = 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 - prodnn = resolution(1)*resolution(2)*resolution(3) - wgt = 1.0_pReal/real(prodnn, pReal) + wgt = 1.0_pReal/real(resolution(1)*resolution(2)*resolution(3), pReal) defgradAim = math_I3 defgradAimOld = math_I3 defgrad_av = math_I3 @@ -297,25 +295,23 @@ program mpie_spectral endif ! calculate xi for the calculation of divergence in Fourier space (middle frequency) - xi_middle(3) = 0.0_pReal + xi_middle(3) = 0.0_pReal !2D case 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) + allocate (workfft(resolution(1)+2,resolution(2),resolution(3),3,3)); workfft = 0.0_pReal + ! 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 + call dfftw_plan_with_nthreads(mpieNumThreadsInt) +! Do r2c/c2r Transform in one step + call dfftw_plan_many_dft_r2c(plan_fft(1),3,(/resolution(1),resolution(2),resolution(3)/),9,& + workfft,(/resolution(1) +2,resolution(2),resolution(3)/),1,(resolution(1) +2)*resolution(2)*resolution(3),& + workfft,(/resolution(1)/2+1,resolution(2),resolution(3)/),1,(resolution(1)/2+1)*resolution(2)*resolution(3),FFTW_PATIENT) + call dfftw_plan_many_dft_c2r(plan_fft(2),3,(/resolution(1),resolution(2),resolution(3)/),9,& + workfft,(/resolution(1)/2+1,resolution(2),resolution(3)/),1,(resolution(1)/2+1)*resolution(2)*resolution(3),& + workfft,(/resolution(1) +2,resolution(2),resolution(3)/),1,(resolution(1) +2)*resolution(2)*resolution(3),FFTW_PATIENT) ! write header of output file open(538,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())& @@ -384,6 +380,7 @@ program mpie_spectral print*, ' ' print '(3(A,I5.5,tr2))', ' Loadcase = ',loadcase, ' Step = ',steps,'Iteration = ',iter cstress_av = 0.0_pReal + workfft = 0.0_pReal !needed because of the padding for FFTW !************************************************************* ! adjust defgrad to fulfill BCs @@ -397,7 +394,7 @@ program mpie_spectral temperature,timeinc,ielem,1_pInt,& cstress,dsde, pstress, dPdF) enddo; enddo; enddo - + ielem = 0_pInt do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1) ielem = ielem + 1_pInt @@ -406,13 +403,13 @@ program mpie_spectral temperature,timeinc,ielem,1_pInt,& cstress,dsde, pstress, dPdF) CPFEM_mode = 2_pInt - pstress_field(i,j,k,:,:) = pstress + workfft(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 + pstress_av(m,n) = sum(workfft(1:resolution(1),:,:,m,n)) * wgt defgrad_av(m,n) = sum(defgrad(:,:,:,m,n)) * wgt enddo; enddo @@ -456,7 +453,6 @@ program mpie_spectral temperature,timeinc,ielem,1_pInt,& cstress,dsde, pstress, dPdF) enddo; enddo; enddo - ielem = 0_pInt do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1) ielem = ielem + 1_pInt @@ -464,31 +460,32 @@ program mpie_spectral defgradold(i,j,k,:,:), defgrad(i,j,k,:,:),& temperature,timeinc,ielem,1_pInt,& cstress,dsde, pstress, dPdF) - pstress_field(i,j,k,:,:) = pstress + workfft(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 + pstress_av(m,n) = sum(workfft(1:resolution(1),:,:,m,n))*wgt enddo; enddo print *, 'Calculating equilibrium using spectral method' err_div = 0.0_pReal; sigma0 = 0.0_pReal - 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 + call dfftw_execute_dft_r2c(plan_fft(1),workfft,workfft) ! FFT of pstress + do m = 1,3 ! L infinity Norm of stress tensor + sigma0 = max(sigma0, sum(abs(workfft(1,1,1,m,:) + (workfft(2,1,1,m,:))*img))) 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 + err_div = (maxval(abs(math_mul33x3_complex(workfft(resolution(1)+1,resolution(2)/2+1,resolution(3)/2+1,:,:)+& ! L infinity Norm of div(stress) + workfft(resolution(1)+2,resolution(2)/2+1,resolution(3)/2+1,:,:)*img,xi_middle)))) + 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,:,:)) + temp33_Complex(m,n) = sum(gamma_hat(i,j,k, m,n,:,:) *(workfft(i*2-1,j,k,:,:)& + + workfft(i*2 ,j,k,:,:)*img)) enddo; enddo - workfft(i,j,k,:,:) = temp33_Complex(:,:) + workfft(i*2-1,j,k,:,:) = real (temp33_Complex) + workfft(i*2 ,j,k,:,:) = aimag(temp33_Complex) enddo; enddo; enddo else ! memory saving version, in-time calculation of gamma_hat @@ -514,30 +511,32 @@ program mpie_spectral 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)) + 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,:,:)) + temp33_Complex(m,n) = sum(gamma_hat(1,1,1,m,n,:,:) *(workfft(i*2-1,j,k,:,:)& + +workfft(i*2 ,j,k,:,:)*img)) enddo; enddo - workfft(i,j,k,:,:) = temp33_Complex(:,:) - enddo - enddo - enddo + workfft(i*2-1,j,k,:,:) = real (temp33_Complex) + workfft(i*2 ,j,k,:,:) = aimag(temp33_Complex) + enddo; enddo; enddo endif - workfft(1,1,1,:,:) = defgrad_av - math_I3 !zero frequency + workfft(1,1,1,:,:) = defgrad_av - math_I3 !zero frequency (real part) + workfft(2,1,1,:,:) = 0.0_pReal !zero frequency (imaginary part) + + call dfftw_execute_dft_c2r(plan_fft(2),workfft,workfft) + defgrad = defgrad + workfft(1:resolution(1),:,:,:,:)*wgt 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 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 deformation gradient ',err_defgrad,' Tol. = ', err_defgrad_tol @@ -560,10 +559,8 @@ program mpie_spectral enddo ! end looping over steps in current loadcase enddo ! end looping over loadcases close(538) - do i=1,2; do m = 1,3; do n = 1,3 - call dfftw_destroy_plan(plan_fft(i,m,n)) - enddo; enddo; enddo - + call dfftw_destroy_plan(plan_fft(1)); call dfftw_destroy_plan(plan_fft(2)) + end program mpie_spectral !******************************************************************** diff --git a/processing/post/make_postprocessingMath b/processing/post/make_postprocessingMath new file mode 100644 index 000000000..cebc877b1 --- /dev/null +++ b/processing/post/make_postprocessingMath @@ -0,0 +1,14 @@ +#!/bin/bash + +# This script is used to compile the python module used for geometry reconstruction. +# It uses the fortran wrapper f2py that is included in the numpy package to construct the +# module postprocessingMath.so out of the fortran code postprocessingMath.f90 +# written by M. Diehl, m.diehl@mpie.de + +#for the generation of the pyf file. +#f2py -m postprocessingMath -h postprocessingMath.pyf --overwrite-signature postprocessingMath.f90 + +# --f90flags="-heap-arrays 500000000" prevents segmentation fault for large arrays +f2py -c --f90flags="-heap-arrays 500000000" \ +postprocessingMath.pyf postprocessingMath.f90 + diff --git a/processing/post/make_reconstruct.py b/processing/post/make_reconstruct.py deleted file mode 100644 index fbbdc808b..000000000 --- a/processing/post/make_reconstruct.py +++ /dev/null @@ -1,15 +0,0 @@ -#!/bin/bash - -# This script is used to compile the python module used for geometry reconstruction. -# It uses the fortran wrapper f2py that is included in the numpy package to construct the -# module reconstruct.so out of the fortran code reconstruct.f90 -# written by M. Diehl, m.diehl@mpie.de - - -#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/postprocessingMath.f90 b/processing/post/postprocessingMath.f90 new file mode 100644 index 000000000..18d4b0f01 --- /dev/null +++ b/processing/post/postprocessingMath.f90 @@ -0,0 +1,525 @@ +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!all function below are taken from math.f90 +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +module math + +real*8, parameter :: pi = 3.14159265358979323846264338327950288419716939937510 + +! *** 3x3 Identity *** + real*8, dimension(3,3), parameter :: math_I3 = & + reshape( (/ & + 1.0,0.0,0.0, & + 0.0,1.0,0.0, & + 0.0,0.0,1.0 /),(/3,3/)) + +contains +!************************************************************************** +! matrix multiplication 33x33 = 3x3 +!************************************************************************** +pure function math_mul33x33(A,B) + + implicit none + + integer i,j + real*8, dimension(3,3), intent(in) :: A,B + real*8, dimension(3,3) :: math_mul33x33 + + forall (i=1:3,j=1:3) math_mul33x33(i,j) = & + A(i,1)*B(1,j) + A(i,2)*B(2,j) + A(i,3)*B(3,j) + return + +end function math_mul33x33 + +!************************************************************************** +! Cramer inversion of 3x3 matrix (subroutine) +!************************************************************************** + PURE SUBROUTINE math_invert3x3(A, InvA, DetA, error) + +! Bestimmung der Determinanten und Inversen einer 3x3-Matrix +! A = Matrix A +! InvA = Inverse of A +! DetA = Determinant of A +! error = logical + + implicit none + + logical, intent(out) :: error + + real*8,dimension(3,3),intent(in) :: A + real*8,dimension(3,3),intent(out) :: InvA + real*8, intent(out) :: DetA + + DetA = A(1,1) * ( A(2,2) * A(3,3) - A(2,3) * A(3,2) )& + - A(1,2) * ( A(2,1) * A(3,3) - A(2,3) * A(3,1) )& + + A(1,3) * ( A(2,1) * A(3,2) - A(2,2) * A(3,1) ) + + if (DetA <= tiny(DetA)) then + error = .true. + else + InvA(1,1) = ( A(2,2) * A(3,3) - A(2,3) * A(3,2) ) / DetA + InvA(2,1) = ( -A(2,1) * A(3,3) + A(2,3) * A(3,1) ) / DetA + InvA(3,1) = ( A(2,1) * A(3,2) - A(2,2) * A(3,1) ) / DetA + + InvA(1,2) = ( -A(1,2) * A(3,3) + A(1,3) * A(3,2) ) / DetA + InvA(2,2) = ( A(1,1) * A(3,3) - A(1,3) * A(3,1) ) / DetA + InvA(3,2) = ( -A(1,1) * A(3,2) + A(1,2) * A(3,1) ) / DetA + + InvA(1,3) = ( A(1,2) * A(2,3) - A(1,3) * A(2,2) ) / DetA + InvA(2,3) = ( -A(1,1) * A(2,3) + A(1,3) * A(2,1) ) / DetA + InvA(3,3) = ( A(1,1) * A(2,2) - A(1,2) * A(2,1) ) / DetA + + error = .false. + endif + return + + END SUBROUTINE math_invert3x3 + +!******************************************************************** +! determinant of a 3x3 matrix +!******************************************************************** + pure function math_det3x3(m) + + implicit none + + real*8, dimension(3,3), intent(in) :: m + real*8 math_det3x3 + + math_det3x3 = m(1,1)*(m(2,2)*m(3,3)-m(2,3)*m(3,2)) & + -m(1,2)*(m(2,1)*m(3,3)-m(2,3)*m(3,1)) & + +m(1,3)*(m(2,1)*m(3,2)-m(2,2)*m(3,1)) + return + + end function math_det3x3 + +!**************************************************************** + pure subroutine math_pDecomposition(FE,U,R,error) +!-----FE = R.U +!**************************************************************** + implicit none + + real*8, intent(in) :: FE(3,3) + real*8, intent(out) :: R(3,3), U(3,3) + logical, intent(out) :: error + real*8 CE(3,3),EW1,EW2,EW3,EB1(3,3),EB2(3,3),EB3(3,3),UI(3,3),det + + error = .false. + ce = math_mul33x33(transpose(FE),FE) + + CALL math_spectral1(CE,EW1,EW2,EW3,EB1,EB2,EB3) + U=DSQRT(EW1)*EB1+DSQRT(EW2)*EB2+DSQRT(EW3)*EB3 + call math_invert3x3(U,UI,det,error) + if (.not. error) R = math_mul33x33(FE,UI) + + return + + end subroutine math_pDecomposition + +!************************************************************************** +! Cramer inversion of 3x3 matrix (function) +!************************************************************************** + pure function math_inv3x3(A) + +! direct Cramer inversion of matrix A. +! returns all zeroes if not possible, i.e. if det close to zero + + implicit none + + real*8,dimension(3,3),intent(in) :: A + real*8 DetA + + real*8,dimension(3,3) :: math_inv3x3 + + math_inv3x3 = 0.0 + + DetA = A(1,1) * ( A(2,2) * A(3,3) - A(2,3) * A(3,2) )& + - A(1,2) * ( A(2,1) * A(3,3) - A(2,3) * A(3,1) )& + + A(1,3) * ( A(2,1) * A(3,2) - A(2,2) * A(3,1) ) + + if (DetA > tiny(DetA)) then + math_inv3x3(1,1) = ( A(2,2) * A(3,3) - A(2,3) * A(3,2) ) / DetA + math_inv3x3(2,1) = ( -A(2,1) * A(3,3) + A(2,3) * A(3,1) ) / DetA + math_inv3x3(3,1) = ( A(2,1) * A(3,2) - A(2,2) * A(3,1) ) / DetA + + math_inv3x3(1,2) = ( -A(1,2) * A(3,3) + A(1,3) * A(3,2) ) / DetA + math_inv3x3(2,2) = ( A(1,1) * A(3,3) - A(1,3) * A(3,1) ) / DetA + math_inv3x3(3,2) = ( -A(1,1) * A(3,2) + A(1,2) * A(3,1) ) / DetA + + math_inv3x3(1,3) = ( A(1,2) * A(2,3) - A(1,3) * A(2,2) ) / DetA + math_inv3x3(2,3) = ( -A(1,1) * A(2,3) + A(1,3) * A(2,1) ) / DetA + math_inv3x3(3,3) = ( A(1,1) * A(2,2) - A(1,2) * A(2,1) ) / DetA + endif + return + + end function math_inv3x3 + +!********************************************************************** +! HAUPTINVARIANTEN HI1M, HI2M, HI3M DER 3X3 MATRIX M +!********************************************************************** + + PURE SUBROUTINE math_hi(M,HI1M,HI2M,HI3M) + implicit none + + real*8, intent(in) :: M(3,3) + real*8, intent(out) :: HI1M, HI2M, HI3M + + HI1M=M(1,1)+M(2,2)+M(3,3) + HI2M=HI1M**2/2.0-(M(1,1)**2+M(2,2)**2+M(3,3)**2)/2.0-M(1,2)*M(2,1)-M(1,3)*M(3,1)-M(2,3)*M(3,2) + HI3M=math_det3x3(M) +! QUESTION: is 3rd equiv det(M) ?? if yes, use function math_det !agreed on YES + return + + END SUBROUTINE math_hi + +!********************************************************************** + pure subroutine math_spectral1(M,EW1,EW2,EW3,EB1,EB2,EB3) +!**** EIGENWERTE UND EIGENWERTBASIS DER SYMMETRISCHEN 3X3 MATRIX M + + implicit none + + real*8, intent(in) :: M(3,3) + real*8, intent(out) :: EB1(3,3),EB2(3,3),EB3(3,3),EW1,EW2,EW3 + real*8 HI1M,HI2M,HI3M,TOL,R,S,T,P,Q,RHO,PHI,Y1,Y2,Y3,D1,D2,D3 + real*8 C1,C2,C3,M1(3,3),M2(3,3),M3(3,3),arg + TOL=1.e-14 + CALL math_hi(M,HI1M,HI2M,HI3M) + R=-HI1M + S= HI2M + T=-HI3M + P=S-R**2.0/3.0 + Q=2.0/27.0*R**3.0-R*S/3.0+T + EB1=0.0 + EB2=0.0 + EB3=0.0 + IF((ABS(P).LT.TOL).AND.(ABS(Q).LT.TOL))THEN +! DREI GLEICHE EIGENWERTE + EW1=HI1M/3.0 + EW2=EW1 + EW3=EW1 +! this is not really correct, but this way U is calculated +! correctly in PDECOMPOSITION (correct is EB?=I) + EB1(1,1)=1.0 + EB2(2,2)=1.0 + EB3(3,3)=1.0 + ELSE + RHO=DSQRT(-3.0*P**3.0)/9.0 + arg=-Q/RHO/2.0 + if(arg.GT.1) arg=1 + if(arg.LT.-1) arg=-1 + PHI=DACOS(arg) + Y1=2*RHO**(1.0/3.0)*DCOS(PHI/3.0) + Y2=2*RHO**(1.0/3.0)*DCOS(PHI/3.0+2.0/3.0*PI) + Y3=2*RHO**(1.0/3.0)*DCOS(PHI/3.0+4.0/3.0*PI) + EW1=Y1-R/3.0 + EW2=Y2-R/3.0 + EW3=Y3-R/3.0 + C1=ABS(EW1-EW2) + C2=ABS(EW2-EW3) + C3=ABS(EW3-EW1) + + IF(C1.LT.TOL) THEN +! EW1 is equal to EW2 + D3=1.0/(EW3-EW1)/(EW3-EW2) + M1=M-EW1*math_I3 + M2=M-EW2*math_I3 + EB3=math_mul33x33(M1,M2)*D3 + + EB1=math_I3-EB3 +! both EB2 and EW2 are set to zero so that they do not +! contribute to U in PDECOMPOSITION + EW2=0.0 + ELSE IF(C2.LT.TOL) THEN +! EW2 is equal to EW3 + D1=1.0/(EW1-EW2)/(EW1-EW3) + M2=M-math_I3*EW2 + M3=M-math_I3*EW3 + EB1=math_mul33x33(M2,M3)*D1 + EB2=math_I3-EB1 +! both EB3 and EW3 are set to zero so that they do not +! contribute to U in PDECOMPOSITION + EW3=0.0 + ELSE IF(C3.LT.TOL) THEN +! EW1 is equal to EW3 + D2=1.0/(EW2-EW1)/(EW2-EW3) + M1=M-math_I3*EW1 + M3=M-math_I3*EW3 + EB2=math_mul33x33(M1,M3)*D2 + EB1=math_I3-EB2 +! both EB3 and EW3 are set to zero so that they do not +! contribute to U in PDECOMPOSITION + EW3=0.0 + ELSE +! all three eigenvectors are different + D1=1.0/(EW1-EW2)/(EW1-EW3) + D2=1.0/(EW2-EW1)/(EW2-EW3) + D3=1.0/(EW3-EW1)/(EW3-EW2) + M1=M-EW1*math_I3 + M2=M-EW2*math_I3 + M3=M-EW3*math_I3 + EB1=math_mul33x33(M2,M3)*D1 + EB2=math_mul33x33(M1,M3)*D2 + EB3=math_mul33x33(M1,M2)*D3 + + END IF + END IF + RETURN + END SUBROUTINE math_spectral1 + + end module math +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +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 + +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +subroutine logstrain_spat(res_x,res_y,res_z,defgrad,logstrain_field) +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + use math + implicit none + integer res_x, res_y, res_z + integer i, j, k + real*8 defgrad(res_x,res_y,res_z,3,3) + real*8 logstrain_field(res_x,res_y,res_z,3,3) + real*8 temp33_Real(3,3), temp33_Real2(3,3) + real*8 eigenvectorbasis(3,3,3) + real*8 eigenvalue(3) + logical errmatinv + + do k = 1, res_z; do j = 1, res_y; do i = 1, res_x + call math_pDecomposition(defgrad(i,j,k,:,:),temp33_Real2,temp33_Real,errmatinv) !store R in temp33_Real + temp33_Real2 = math_inv3x3(temp33_Real) + temp33_Real = math_mul33x33(defgrad(i,j,k,:,:),temp33_Real2) ! v = F o inv(R), store in temp33_Real2 + call math_spectral1(temp33_Real, eigenvalue(1), eigenvalue(2), eigenvalue(3),& + eigenvectorbasis(1,:,:), eigenvectorbasis(2,:,:), eigenvectorbasis(3,:,:)) + eigenvalue = log(sqrt(eigenvalue)) + logstrain_field(i,j,k,:,:) = eigenvalue(1)*eigenvectorbasis(1,:,:)+& + eigenvalue(2)*eigenvectorbasis(2,:,:)+& + eigenvalue(3)*eigenvectorbasis(3,:,:) + enddo; enddo; enddo + end subroutine logstrain_spat + + !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +subroutine logstrain_mat(res_x,res_y,res_z,defgrad,logstrain_field) +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + use math + implicit none + integer res_x, res_y, res_z + integer i, j, k + real*8 defgrad(res_x,res_y,res_z,3,3) + real*8 logstrain_field(res_x,res_y,res_z,3,3) + real*8 temp33_Real(3,3), temp33_Real2(3,3) + real*8 eigenvectorbasis(3,3,3) + real*8 eigenvalue(3) + logical errmatinv + + do k = 1, res_z; do j = 1, res_y; do i = 1, res_x + call math_pDecomposition(defgrad(i,j,k,:,:),temp33_Real,temp33_Real2,errmatinv) !store U in temp33_Real + call math_spectral1(temp33_Real, eigenvalue(1), eigenvalue(2), eigenvalue(3),& + eigenvectorbasis(1,:,:), eigenvectorbasis(2,:,:), eigenvectorbasis(3,:,:)) + eigenvalue = log(sqrt(eigenvalue)) + logstrain_field(i,j,k,:,:) = eigenvalue(1)*eigenvectorbasis(1,:,:)+& + eigenvalue(2)*eigenvectorbasis(2,:,:)+& + eigenvalue(3)*eigenvectorbasis(3,:,:) + enddo; enddo; enddo + end subroutine logstrain_mat + +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +subroutine calculate_cauchy(res_x,res_y,res_z,defgrad,p_stress,c_stress) +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + use math + implicit none + integer res_x, res_y, res_z + integer i, j, k + real*8 defgrad(res_x,res_y,res_z,3,3) + real*8 p_stress(res_x,res_y,res_z,3,3) + real*8 c_stress(res_x,res_y,res_z,3,3) + real*8 jacobi + c_stress = 0.0 + do k = 1, res_z; do j = 1, res_y; do i = 1, res_x + jacobi = math_det3x3(defgrad(i,j,k,:,:)) + c_stress(i,j,k,:,:) = matmul(p_stress(i,j,k,:,:),transpose(defgrad(i,j,k,:,:)))/jacobi + enddo; enddo; enddo +end subroutine calculate_cauchy + +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +subroutine calculate_mises(res_x,res_y,res_z,tensor,vm) +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + implicit none + integer res_x, res_y, res_z + integer i, j, k + real*8 tensor(res_x,res_y,res_z,3,3) + real*8 vm(res_x,res_y,res_z,1) + real*8 deviator(3,3) + real*8 delta(3,3) + real*8 J_2 + + delta =0.0 + delta(1,1) = 1.0 + delta(2,2) = 1.0 + delta(3,3) = 1.0 + do k = 1, res_z; do j = 1, res_y; do i = 1, res_x + deviator = tensor(i,j,k,:,:) - 1.0/3.0*tensor(i,j,k,1,1)*tensor(i,j,k,2,2)*tensor(i,j,k,3,3)*delta + J_2 = deviator(1,1)*deviator(2,2)& + + deviator(2,2)*deviator(3,3)& + + deviator(1,1)*deviator(3,3)& + - (deviator(1,2))**2& + - (deviator(2,3))**2& + - (deviator(1,3))**2 + vm(i,j,k,:) = sqrt(3*J_2) + enddo; enddo; enddo +end subroutine calculate_mises diff --git a/processing/post/postprocessingMath.pyf b/processing/post/postprocessingMath.pyf new file mode 100644 index 000000000..a983d30a0 --- /dev/null +++ b/processing/post/postprocessingMath.pyf @@ -0,0 +1,105 @@ +! -*- f90 -*- +! Note: the context of this file is case sensitive. + +python module postprocessingMath ! in + interface ! in :postprocessingMath + module math ! in :postprocessingMath:postprocessingMath.f90 + real*8 parameter,optional,dimension(3,3) :: math_i3=reshape((/1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0/),(/3,3/)) + real*8 parameter,optional :: pi=3.14159265359 + function math_mul33x33(a,b) ! in :postprocessingMath:postprocessingMath.f90:math + real*8 dimension(3,3),intent(in) :: a + real*8 dimension(3,3),intent(in) :: b + real*8 dimension(3,3) :: math_mul33x33 + end function math_mul33x33 + subroutine math_invert3x3(a,inva,deta,error) ! in :postprocessingMath:postprocessingMath.f90:math + real*8 dimension(3,3),intent(in) :: a + real*8 dimension(3,3),intent(out) :: inva + real*8 intent(out) :: deta + logical intent(out) :: error + end subroutine math_invert3x3 + function math_det3x3(m) ! in :postprocessingMath:postprocessingMath.f90:math + real*8 dimension(3,3),intent(in) :: m + real*8 :: math_det3x3 + end function math_det3x3 + subroutine math_pdecomposition(fe,u,r,error) ! in :postprocessingMath:postprocessingMath.f90:math + real*8 dimension(3,3),intent(in) :: fe + real*8 dimension(3,3),intent(out) :: u + real*8 dimension(3,3),intent(out) :: r + logical intent(out) :: error + end subroutine math_pdecomposition + function math_inv3x3(a) ! in :postprocessingMath:postprocessingMath.f90:math + real*8 dimension(3,3),intent(in) :: a + real*8 dimension(3,3) :: math_inv3x3 + end function math_inv3x3 + subroutine math_hi(m,hi1m,hi2m,hi3m) ! in :postprocessingMath:postprocessingMath.f90:math + real*8 dimension(3,3),intent(in) :: m + real*8 intent(out) :: hi1m + real*8 intent(out) :: hi2m + real*8 intent(out) :: hi3m + end subroutine math_hi + subroutine math_spectral1(m,ew1,ew2,ew3,eb1,eb2,eb3) ! in :postprocessingMath:postprocessingMath.f90:math + real*8 dimension(3,3),intent(in) :: m + real*8 intent(out) :: ew1 + real*8 intent(out) :: ew2 + real*8 intent(out) :: ew3 + real*8 dimension(3,3),intent(out) :: eb1 + real*8 dimension(3,3),intent(out) :: eb2 + real*8 dimension(3,3),intent(out) :: eb3 + end subroutine math_spectral1 + end module math + subroutine mesh(res_x,res_y,res_z,geomdim,defgrad_av,centroids,nodes) ! in :postprocessingMath:postprocessingMath.f90 + integer intent(in) :: res_x + integer intent(in) :: res_y + integer intent(in) :: res_z + real*8 dimension(3),intent(in) :: geomdim + real*8 dimension(3,3),intent(in) :: defgrad_av + real*8 dimension(res_x,res_y,res_z,3),intent(in),depend(res_x,res_y,res_z) :: centroids + real*8 dimension(res_x + 2,res_y + 2,res_z +2 ,3),depend(res_x,res_y,res_z) :: centroids + real*8 dimension(res_x + 1,res_y + 1,res_z + 1,3),intent(out),depend(res_x,res_y,res_z) :: nodes + end subroutine mesh + subroutine deformed(res_x,res_y,res_z,geomdim,defgrad,defgrad_av,coord_avgCorner) ! in :postprocessingMath:postprocessingMath.f90 + integer intent(in) :: res_x + integer intent(in) :: res_y + integer intent(in) :: res_z + real*8 dimension(3),intent(in) :: geomdim + real*8 dimension(3,3),intent(in) :: defgrad_av + real*8 dimension(res_x,res_y,res_z,3),intent(out),depend(res_x,res_y,res_z) :: coord_avgCorner + real*8 dimension(8,6,res_x,res_y,res_z,3),depend(res_x,res_y,res_z) :: coord + real*8 dimension(8,res_x,res_y,res_z,3),depend(res_x,res_y,res_z) :: coord_avgOrder + real*8 dimension(res_x,res_y,res_z,3,3),intent(in),depend(res_x,res_y,res_z) :: defgrad + end subroutine deformed + subroutine logstrain_mat(res_x,res_y,res_z,defgrad,logstrain_field) ! in :postprocessingMath:postprocessingMath.f90 + integer intent(in) :: res_x + integer intent(in) :: res_y + integer intent(in) :: res_z + real*8 dimension(res_x,res_y,res_z,3,3),intent(in),depend(res_x,res_y,res_z) :: defgrad + real*8 dimension(res_x,res_y,res_z,3,3),intent(out),depend(res_x,res_y,res_z) :: logstrain_field + end subroutine logstrain_mat + subroutine logstrain_spat(res_x,res_y,res_z,defgrad,logstrain_field) ! in :postprocessingMath:postprocessingMath.f90 + integer intent(in) :: res_x + integer intent(in) :: res_y + integer intent(in) :: res_z + real*8 dimension(res_x,res_y,res_z,3,3),intent(in),depend(res_x,res_y,res_z) :: defgrad + real*8 dimension(res_x,res_y,res_z,3,3),intent(out),depend(res_x,res_y,res_z) :: logstrain_field + end subroutine logstrain_spat + subroutine calculate_cauchy(res_x,res_y,res_z,defgrad,p_stress,c_stress) ! in :postprocessingMath:postprocessingMath.f90 + integer intent(in) :: res_x + integer intent(in) :: res_y + integer intent(in) :: res_z + real*8 dimension(res_x,res_y,res_z,3,3),intent(in),depend(res_x,res_y,res_z) :: defgrad + real*8 dimension(res_x,res_y,res_z,3,3),intent(in),depend(res_x,res_y,res_z) :: p_stress + real*8 dimension(res_x,res_y,res_z,3,3),intent(out),depend(res_x,res_y,res_z) :: c_stress + end subroutine calculate_cauchy + subroutine calculate_mises(res_x,res_y,res_z,tensor,vm) ! in :postprocessingMath:postprocessingMath.f90 + integer intent(in) :: res_x + integer intent(in) :: res_y + integer intent(in) :: res_z + real*8 dimension(res_x,res_y,res_z,3,3),intent(in),depend(res_x,res_y,res_z) :: tensor + real*8 dimension(res_x,res_y,res_z,1),intent(out),depend(res_x,res_y,res_z) :: vm + end subroutine calculate_mises + end interface +end python module postprocessingMath + +! This file was auto-generated with f2py (version:2_5972). +! See http://cens.ioc.ee/projects/f2py2e/ +! modified by m.diehl diff --git a/processing/post/reconstruct.f90 b/processing/post/reconstruct.f90 deleted file mode 100644 index c1d825bdc..000000000 --- a/processing/post/reconstruct.f90 +++ /dev/null @@ -1,157 +0,0 @@ -!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -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 add25f670..fe7f66c96 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, reconstruct +import os,sys,re,array,struct,numpy, time, postprocessingMath class vector: x,y,z = [None,None,None] @@ -159,145 +159,6 @@ def calculateCauchyStress(p_stress,defgrad,res): c_stress[x,y,z] = numpy.dot(p_stress[x,y,z],numpy.transpose(defgrad[x,y,z]))/jacobi return c_stress -#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -def calculateVonMises(tensor,res): -#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - vonMises = numpy.zeros([res[0],res[1],res[2]],'d') - deviator = numpy.zeros([3,3],'d') - delta = numpy.zeros([3,3],'d') - delta[0,0] = 1.0 - delta[1,1] = 1.0 - delta[2,2] = 1.0 - for z in range(res[2]): - for y in range(res[1]): - for x in range(res[0]): - deviator = tensor[x,y,z] - 1.0/3.0*tensor[x,y,z,0,0]*tensor[x,y,z,1,1]*tensor[x,y,z,2,2]*delta - J_2 = deviator[0,0]*deviator[1,1]\ - + deviator[1,1]*deviator[2,2]\ - + deviator[0,0]*deviator[2,2]\ - - (deviator[0,1])**2\ - - (deviator[1,2])**2\ - - (deviator[0,2])**2 - vonMises[x,y,z] = numpy.sqrt(3*J_2) - return vonMises - -#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -def mesh(res,geomdim,defgrad_av,centroids): -#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - - neighbor = numpy.array([[0, 0, 0], - [1, 0, 0], - [1, 1, 0], - [0, 1, 0], - [0, 0, 1], - [1, 0, 1], - [1, 1, 1], - [0, 1, 1]]) - - wrappedCentroids = numpy.zeros([res[0]+2,res[1]+2,res[2]+2,3],'d') - nodes = numpy.zeros([res[0]+1,res[1]+1,res[2]+1,3],'d') - wrappedCentroids[1:-1,1:-1,1:-1] = centroids - diag = numpy.ones(3,'i') - shift = numpy.zeros(3,'i') - lookup = numpy.zeros(3,'i') - - for k in range(res[2]+2): - for j in range(res[1]+2): - for i in range(res[0]+2): - if (k==0 or k==res[2]+1 or \ - j==0 or j==res[1]+1 or \ - i==0 or i==res[0]+1 ): - me = numpy.array([i,j,k],'i') - shift = numpy.sign(res+diag-2*me)*(numpy.abs(res+diag-2*me)/(res+diag)) - lookup = me-diag+shift*res - wrappedCentroids[i,j,k] = centroids[lookup[0],lookup[1],lookup[2]]- \ - numpy.dot(defgrad_av, shift*geomdim) - for k in range(res[2]+1): - for j in range(res[1]+1): - 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 - - return nodes - -# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -def centroids(res,geomdimension,defgrad): -# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - corner = numpy.array([[0, 0, 0], - [1, 0, 0], - [1, 1, 0], - [0, 1, 0], - [1, 1, 1], - [0, 1, 1], - [0, 0, 1], - [1, 0, 1]]) - step = numpy.array([[ 1, 1, 1], - [-1, 1, 1], - [-1,-1, 1], - [ 1,-1, 1], - [-1,-1,-1], - [ 1,-1,-1], - [ 1, 1,-1], - [-1, 1,-1]]) - - order = numpy.array([[0, 1, 2], - [0, 2, 1], - [1, 0, 2], - [1, 2, 0], - [2, 0, 1], - [2, 1, 0]]) - - cornerCoords = numpy.zeros([8,res[0],res[1],res[2],3], 'd') - coord = numpy.zeros([8,6,res[0],res[1],res[2],3], 'd') - centroids = numpy.zeros([res[0],res[1],res[2],3], 'd') - myStep = numpy.zeros(3,'d') - rear = numpy.zeros(3, 'i') - init = numpy.zeros(3, 'i') - oppo = numpy.zeros(3, 'i') - me = numpy.zeros(3, 'i') - 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) - for o in range(6): # orders - for k in range(init[order[o,2]],oppo[order[o,2]]+step[s,order[o,2]],step[s,order[o,2]]): - rear[order[o,1]] = init[order[o,1]] - for j in range(init[order[o,1]],oppo[order[o,1]]+step[s,order[o,1]],step[s,order[o,1]]): - rear[order[o,0]] = init[order[o,0]] - for i in range(init[order[o,0]],oppo[order[o,0]]+step[s,order[o,0]],step[s,order[o,0]]): - me[order[o,0]] = i - me[order[o,1]] = j - me[order[o,2]] = k - if (numpy.all(me == init)): - coord[s,o,me[0],me[1],me[2]] = geomdimension * (numpy.dot(defgrad_av,corner[s]) + \ - numpy.dot(defgrad[me[0],me[1],me[2]],0.5*step[s]/res)) - else: - myStep = (me-rear)*geomdimension/res - coord[s,o,me[0],me[1],me[2]] = coord[s,o,rear[0],rear[1],rear[2]]+ \ - 0.5*numpy.dot(defgrad[me[0],me[1],me[2]] + \ - defgrad[rear[0],rear[1],rear[2]],myStep) - - rear[:] = me[:] - cornerCoords[s] = numpy.average(coord[s],0) - for k in range(res[2]): - for j in range(res[1]): - for i in range(res[0]): - parameter_coords=(2.0*numpy.array([i,j,k])-res+fones)/(res-fones) - pos = (fones + parameter_coords) - neg = (fones - parameter_coords) - centroids[i,j,k] = ( cornerCoords[0,i,j,k] *neg[0]*neg[1]*neg[2]\ - + cornerCoords[1,i,j,k] *pos[0]*neg[1]*neg[2]\ - + cornerCoords[2,i,j,k] *pos[0]*pos[1]*neg[2]\ - + cornerCoords[3,i,j,k] *neg[0]*pos[1]*neg[2]\ - + cornerCoords[4,i,j,k] *pos[0]*pos[1]*pos[2]\ - + cornerCoords[5,i,j,k] *neg[0]*pos[1]*pos[2]\ - + cornerCoords[6,i,j,k] *neg[0]*neg[1]*pos[2]\ - + cornerCoords[7,i,j,k] *pos[0]*neg[1]*pos[2])*0.125 - return centroids, defgrad_av # function writes scalar values to a mesh (geometry) def writeVtkAscii(filename,geometry,scalar,resolution): @@ -445,7 +306,7 @@ print 'Post Processing for Material subroutine for BVP solution using spectral m print '*********************************************************************************\n' #reading in the header of the results file -name = 'dipl32' +name = 'dipl10' p = MPIEspectral_result(name+'.spectralOut') p.extrapolation('') print p @@ -474,14 +335,17 @@ res_z=p.resolution[2] 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) + defgrad = readTensor(p.resolution,p.file,p.N_element_scalars,c_pos,16) + logstrain = postprocessingMath.logstrain_mat(res_x,res_y,res_z,defgrad) + logstrain2 = postprocessingMath.logstrain_spat(res_x,res_y,res_z,defgrad) 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) + c_stress = postprocessingMath.calculate_cauchy(res_x,res_y,res_z,defgrad,p_stress) + vm = postprocessingMath.calculate_mises(res_x,res_y,res_z,c_stress) 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) + centroids_coord = postprocessingMath.deformed(res_x,res_y,res_z,p.dimension,defgrad,defgrad_av) + ms = postprocessingMath.mesh(p.resolution[0],p.resolution[1],p.resolution[2],p.dimension,defgrad_av,centroids_coord) + writeVtkAscii(name+'-mesh-1-%s.vtk'%i,ms,logstrain[:,:,:,1,2],p.resolution) + writeVtkAscii(name+'-mesh-2-%s.vtk'%i,ms,logstrain2[:,:,:,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()