From 438bf95105e9b149c00c1df85d39239e9406c204 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 14 Feb 2011 17:21:31 +0000 Subject: [PATCH] added ft-based reconstruction of deformed configuration to postprocessingMath.f90 and postprocessingMath.pyf also added function to calculate average of tensor some polishing in mpie_spectral.f90, added sanity check to see im resolution is power of 2 --- code/mpie_spectral.f90 | 7 +- processing/post/make_postprocessingMath | 5 +- processing/post/postprocessingMath.f90 | 129 ++++++++++++++++++++++++ processing/post/postprocessingMath.pyf | 64 ++++-------- processing/post/spectral_post.py | 8 +- 5 files changed, 161 insertions(+), 52 deletions(-) diff --git a/code/mpie_spectral.f90 b/code/mpie_spectral.f90 index 087f79b58..513eeb9ff 100644 --- a/code/mpie_spectral.f90 +++ b/code/mpie_spectral.f90 @@ -36,7 +36,7 @@ program mpie_spectral implicit none include 'fftw3.f' !header file for fftw3 (declaring variables). Library files are also needed - +! compile FFTW 3.2.2 with ./configure --enable-threads ! variables to read from loadcase and geom file real(pReal), dimension(9) :: valuevector ! stores information temporarily from loadcase file integer(pInt), parameter :: maxNchunksInput = 24 ! 4 identifiers, 18 values for the matrices and 2 scalars @@ -232,7 +232,9 @@ program mpie_spectral if (gotDimension .and. gotHomogenization .and. gotResolution) exit enddo 100 close(unit) - + + if(mod(resolution(1),2)/=0 .or. mod(resolution(2),2)/=0 .or. mod(resolution(3),2)/=0) call IO_error(102) !!ToDo: add correct error to IO + print '(a,/,i4,i4,i4)','resolution a b c', resolution print '(a,/,f6.1,f6.1,f6.1)','dimension x y z', geomdimension print *,'homogenization',homog @@ -305,7 +307,6 @@ program mpie_spectral ! 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 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) diff --git a/processing/post/make_postprocessingMath b/processing/post/make_postprocessingMath index cebc877b1..534fc4631 100644 --- a/processing/post/make_postprocessingMath +++ b/processing/post/make_postprocessingMath @@ -5,10 +5,11 @@ # 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. +#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 +# use ./configure --enable-portable-binary --enable-shared for the compilation of fftw3.2.2 f2py -c --f90flags="-heap-arrays 500000000" \ -postprocessingMath.pyf postprocessingMath.f90 +postprocessingMath.pyf postprocessingMath.f90 -L ./libfftw3.a diff --git a/processing/post/postprocessingMath.f90 b/processing/post/postprocessingMath.f90 index 18d4b0f01..05f239166 100644 --- a/processing/post/postprocessingMath.f90 +++ b/processing/post/postprocessingMath.f90 @@ -268,6 +268,8 @@ end function math_mul33x33 END SUBROUTINE math_spectral1 end module math + +!subroutines below are for postprocessing with python !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ subroutine mesh(res_x,res_y,res_z,geomdim,defgrad_av,centroids,nodes) !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ @@ -426,6 +428,133 @@ subroutine deformed(res_x,res_y,res_z,geomdim,defgrad,defgrad_av,coord_avgCorner enddo; enddo; enddo end subroutine deformed +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +subroutine deformed_fft(res_x,res_y,res_z,geomdim,defgrad,defgrad_av,scaling,coords) +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + + implicit none + integer res_x, res_y, res_z + real*8 geomdim(3) + real*8 defgrad(res_x,res_y,res_z,3,3) + complex*16 defgrad_fft(res_x,res_y,res_z,3,3) + real*8 defgrad_av(3,3) + real*8 scaling + real*8 coords(res_x, res_y,res_z,3) + complex*16 coords_fft(res_x/2+1,res_y,res_z,3) + real*8 waves(res_x/2+1,res_y,res_z,3) + include 'fftw3.f' !header file for fftw3 (declaring variables). Library files are also needed + integer*8 :: plan_fft(2) + real*8, parameter :: pi = 3.14159265358979323846264338327950288419716939937510 + + real*8 zero + real*8 temp33_Real(3,3) + integer i, j, k, ierr + integer k_s(3) + real*8 step(3) + complex*16 img + + img = cmplx(0.0,1.0) + + do k = 1, res_z + k_s(3) = k-1 + if(k > res_z/2+1) k_s(3) = k_s(3)-res_z + do j = 1, res_y + k_s(2) = j-1 + if(j > res_y/2+1) k_s(2) = k_s(2)-res_y + do i = 1, res_x/2+1 + k_s(1) = i-1 + waves(i,j,k,:) = real(k_s)/geomdim + enddo; enddo; enddo + + + call dfftw_plan_many_dft(plan_fft(1),3,(/res_x,res_y,res_z/),9,& + defgrad_fft,(/res_x,res_y,res_z/),1,res_x*res_y*res_z,& + defgrad_fft,(/res_x,res_y,res_z/),1,res_x*res_y*res_z,FFTW_FORWARD,FFTW_PATIENT) + + call dfftw_plan_many_dft_c2r(plan_fft(2),3,(/res_x,res_y,res_z/),3,& + coords_fft,(/res_x/2+1,res_y,res_z/),1,(res_x/2+1)*res_y*res_z,& + coords, (/res_x, res_y,res_z/),1, res_x* res_y*res_z,FFTW_PATIENT) + + coords_fft=0.0 + defgrad_fft = defgrad + + call dfftw_execute_dft(plan_fft(1), defgrad_fft, defgrad_fft) + + do k = 1, res_z; do j = 1, res_y; do i = 1, res_x/2+1 + if(i/=1) then + coords_fft(i,j,k,1) = defgrad_fft(i,j,k,1,1)/(waves(i,j,k,1)*img) + coords_fft(i,j,k,2) = defgrad_fft(i,j,k,2,1)/(waves(i,j,k,1)*img) + coords_fft(i,j,k,3) = defgrad_fft(i,j,k,3,1)/(waves(i,j,k,1)*img) + endif + if(j/=1) then + coords_fft(i,j,k,1) = coords_fft(i,j,k,1) + defgrad_fft(i,j,k,1,2)/(waves(i,j,k,2)*img) + coords_fft(i,j,k,2) = coords_fft(i,j,k,2) + defgrad_fft(i,j,k,2,2)/(waves(i,j,k,2)*img) + coords_fft(i,j,k,3) = coords_fft(i,j,k,3) + defgrad_fft(i,j,k,3,2)/(waves(i,j,k,2)*img) + endif + if(k/=1) then + coords_fft(i,j,k,1) = coords_fft(i,j,k,1) + defgrad_fft(i,j,k,1,3)/(waves(i,j,k,3)*img) + coords_fft(i,j,k,2) = coords_fft(i,j,k,2) + defgrad_fft(i,j,k,2,3)/(waves(i,j,k,3)*img) + coords_fft(i,j,k,3) = coords_fft(i,j,k,3) + defgrad_fft(i,j,k,3,3)/(waves(i,j,k,3)*img) + endif + enddo; enddo; enddo + + call dfftw_execute_dft_c2r(plan_fft(2), coords_fft, coords) + + coords = coords/real(res_x*res_y*res_z) + + do k = 1, res_z; do j = 1, res_y; do i = 1, res_x + coords(i,j,k,:) = coords(i,j,k,:)/(geomdim*2.0*pi) + enddo; enddo; enddo + + step(1) = geomdim(1)/real(res_x) + step(2) = geomdim(2)/real(res_y) + step(3) = geomdim(3)/real(res_z) + + + temp33_Real(1,:) = matmul(defgrad_av,step/2.0) & + - matmul(defgrad_av,(/zero,zero,step(3)/)) ! start below origin + + do k = 1, res_z; do j = 1, res_y; do i = 1, res_x + if((j==1).and.(i==1)) then + temp33_Real(1,:) = temp33_Real(1,:) + matmul(defgrad_av,& + (/zero,zero,step(3)/)) + temp33_Real(2,:) = temp33_Real(1,:) + temp33_Real(3,:) = temp33_Real(1,:) + coords(i,j,k,:) = coords(i,j,k,:) + temp33_Real(1,:) + else + if(i==1) then + temp33_Real(2,:) = temp33_Real(2,:) + matmul(defgrad_av,& + (/zero,step(2),zero/)) + temp33_Real(3,:) = temp33_Real(2,:) + coords(i,j,k,:) = coords(i,j,k,:) + temp33_Real(2,:) + else + temp33_Real(3,:) = temp33_Real(3,:) + matmul(defgrad_av,& + (/step(1),zero,zero/)) + coords(i,j,k,:) = coords(i,j,k,:) + temp33_Real(3,:) + endif + endif + enddo; enddo; enddo + +end subroutine deformed_fft + +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +subroutine tensor_avg(res_x,res_y,res_z,tensor,avg) +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + + implicit none + integer res_x, res_y, res_z + real*8 tensor(res_x,res_y,res_z,3,3) + real*8 avg(3,3) + real*8 wgt + integer m,n + + wgt = 1/real(res_x*res_y*res_z) + + do m = 1,3; do n = 1,3 + avg(m,n) = sum(tensor(:,:,:,m,n)) * wgt + enddo; enddo +end subroutine tensor_avg + !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ subroutine logstrain_spat(res_x,res_y,res_z,defgrad,logstrain_field) !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ diff --git a/processing/post/postprocessingMath.pyf b/processing/post/postprocessingMath.pyf index a983d30a0..442fa46a8 100644 --- a/processing/post/postprocessingMath.pyf +++ b/processing/post/postprocessingMath.pyf @@ -3,50 +3,6 @@ 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 @@ -68,6 +24,26 @@ python module postprocessingMath ! in 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 deformed_fft(res_x,res_y,res_z,geomdim,defgrad,defgrad_av,scaling,coords) ! in :postprocessingMath:postprocessingMath.f90 + integer intent(in) :: res_x + integer intent(in) :: res_y + integer intent(in) :: res_z + real*8 intent(in) :: scaling + 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) :: coords + complex*16 dimension(res_x/2+1,res_y,res_z,3),depend(res_x,res_y,res_z) :: coords_fft + real*8 dimension(res_x,res_y,res_z,3,3),intent(in),depend(res_x,res_y,res_z) :: defgrad + complex*16 dimension(res_x,res_y,res_z,3,3),depend(res_x,res_y,res_z) :: defgrad_fft + real*8 dimension(res_x/2+1,res_y,res_z,3),depend(res_x,res_y,res_z) :: waves + end subroutine deformed_fft + subroutine tensor_avg(res_x,res_y,res_z,tensor,avg) ! 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(3,3),intent(out) :: avg + end subroutine tensor_avg 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 diff --git a/processing/post/spectral_post.py b/processing/post/spectral_post.py index fe7f66c96..692d09407 100644 --- a/processing/post/spectral_post.py +++ b/processing/post/spectral_post.py @@ -341,11 +341,13 @@ for i in range(249,250): p_stress = readTensor(p.resolution,p.file,p.N_element_scalars,c_pos,58) 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) + defgrad_av = postprocessingMath.tensor_avg(res_x,res_y,res_z,defgrad) centroids_coord = postprocessingMath.deformed(res_x,res_y,res_z,p.dimension,defgrad,defgrad_av) + centroids_coord2 = postprocessingMath.deformed_fft(res_x,res_y,res_z,p.dimension,defgrad,defgrad_av,1.0) 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) + ms2 = postprocessingMath.mesh(p.resolution[0],p.resolution[1],p.resolution[2],p.dimension,defgrad_av,centroids_coord2) + writeVtkAscii(name+'-mesh-usual-%s.vtk'%i,ms,logstrain[:,:,:,1,2],p.resolution) + writeVtkAscii(name+'-mesh-fft-%s.vtk'%i,ms2,logstrain[:,:,:,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()