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

mpie_spectral.f90: changed fourier transormation, now using the advanced interface to do the tranformation both ways with just one call. now also using the same variable for pk-stress in both domains and the change of deformation gradient in both domains. 

postprocessing: renamed name of python/f2py modul from "reconstruct" to "postprocessingMath", added some numerical operations to use for postprocessing.
This commit is contained in:
Martin Diehl 2011-02-09 17:47:28 +00:00
parent 803e1a8c05
commit 3d7fad6ba9
7 changed files with 711 additions and 378 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

99
code/mpie_spectral.f90
View File

@ -60,7 +60,7 @@ 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.
@ -68,21 +68,19 @@ program mpie_spectral
pstress, pstress_av, cstress_av, defgrad_av,&
defgradAim, defgradAimOld, defgradAimCorr, defgradAimCorrPrev,&
mask_stress, mask_defgrad
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 ! 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
@ -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
@ -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
! 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
@ -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
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
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,26 +511,28 @@ 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)))
@ -560,9 +559,7 @@ 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

14
processing/post/make_postprocessingMath Normal file
View File

@ -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

15
processing/post/make_reconstruct.py
View File

@ -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

525
processing/post/postprocessingMath.f90 Normal file
View File

@ -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

105
processing/post/postprocessingMath.pyf Normal file
View File

@ -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

157
processing/post/reconstruct.f90
View File

@ -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

156
processing/post/spectral_post.py
View File

@ -6,7 +6,7 @@
# computed using the FEM solvers. Until now, its capable to handle elements with one IP in a regular order
# 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
@ -475,13 +336,16 @@ 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)
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()