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

@@ -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
@@ -233,6 +233,8 @@ program mpie_spectral
  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)   

processing/post/make_postprocessingMath

@@ -9,6 +9,7 @@
 #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
 

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)
 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

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

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)
+  ms2 = postprocessingMath.mesh(p.resolution[0],p.resolution[1],p.resolution[2],p.dimension,defgrad_av,centroids_coord2)
-  writeVtkAscii(name+'-mesh-2-%s.vtk'%i,ms,logstrain2[:,:,:,1,2],p.resolution)
+  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()