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

added functions math_mul33x3_complex math_Plain99to3333 to math.f90. 

mpie_spectral.f90 uses both functions
math_Plain99to3333 is used for inversion of c0
math_mul33x3_complex is used for equilibrium check in fourier space
also did some cleaning up on mpie_spectral.f90
This commit is contained in:
Martin Diehl 2010-09-22 12:04:43 +00:00
parent 80618c9814
commit 6ea8623f65
2 changed files with 140 additions and 117 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

34
code/math.f90
View File

@ -509,6 +509,24 @@ real(pReal), dimension(4,36), parameter :: math_symOperations = &
endfunction endfunction
!**************************************************************************
! matrix multiplication complex(33) x real(3) = complex(3)
!**************************************************************************
pure function math_mul33x3_complex(A,B)
use prec, only: pReal, pInt
implicit none
integer(pInt) i
complex(pReal), dimension(3,3), intent(in) :: A
real(pReal), dimension(3), intent(in) :: B
complex(pReal), dimension(3) :: math_mul33x3_complex
forall (i=1:3) math_mul33x3_complex(i) = A(i,1)*B(1) + A(i,2)*B(2) + A(i,3)*B(3)
return
endfunction
!************************************************************************** !**************************************************************************
! matrix multiplication 66x6 = 6 ! matrix multiplication 66x6 = 6
@ -1165,7 +1183,23 @@ pure function math_transpose3x3(A)
endfunction endfunction
!********************************************************************
! plain matrix 9x9 into 3x3x3x3 tensor
!********************************************************************
pure function math_Plain99to3333(m99)
use prec, only: pReal,pInt
implicit none
real(pReal), dimension(9,9), intent(in) :: m99
real(pReal), dimension(3,3,3,3) :: math_Plain99to3333
integer(pInt) i,j
forall (i=1:9,j=1:9) math_Plain99to3333(mapPlain(1,i),mapPlain(2,i),&
mapPlain(1,j),mapPlain(2,j)) = m99(i,j)
return
endfunction
!******************************************************************** !********************************************************************
! convert symmetric 3x3x3x3 tensor into Mandel matrix 6x6 ! convert symmetric 3x3x3x3 tensor into Mandel matrix 6x6

159
code/mpie_spectral.f90
View File

@ -33,60 +33,59 @@ program mpie_spectral
implicit none implicit none
include 'fftw3.f' !header file for fftw3 (declaring variables). Library file is also needed include 'fftw3.f' !header file for fftw3 (declaring variables). Library file is also needed
!variables to read in from loadcase and mesh file ! variables to read from loadcase and mesh file
real(pReal), dimension(9) :: valuevector ! stores information temporarily from loadcase 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 integer(pInt), parameter :: maxNchunksInput = 24 ! 4 identifiers, 18 values for the matrices and 2 scalars
integer(pInt), dimension (1+maxNchunksInput*2) :: posInput integer(pInt), dimension (1+maxNchunksInput*2) :: posInput
integer(pInt), parameter :: maxNchunksMesh = 7 ! 4 identifiers, 3 values integer(pInt), parameter :: maxNchunksMesh = 7 ! 4 identifiers, 3 values
integer(pInt), dimension (1+2*maxNchunksMesh) :: posMesh integer(pInt), dimension (1+2*maxNchunksMesh) :: posMesh
integer(pInt) unit, N_l, N_s, N_t, N_n ! numbers of identifiers integer(pInt) unit, N_l, N_s, N_t, N_n ! numbers of identifiers
character(len=1024) path, line
logical gotResolution,gotDimension,gotHomogenization logical gotResolution,gotDimension,gotHomogenization
logical, dimension(9) :: bc_maskvector logical, dimension(9) :: bc_maskvector
character(len=1024) path, line
! variables storing information from loadcase file ! variables storing information from loadcase file
integer(pInt) N_Loadcases, steps
integer(pInt), dimension(:), allocatable :: bc_steps ! number of steps
real(pReal) timeinc real(pReal) timeinc
real(pReal), dimension (:,:,:), allocatable :: bc_velocityGrad, & real(pReal), dimension (:,:,:), allocatable :: bc_velocityGrad, &
bc_stress ! velocity gradient and stress BC bc_stress ! velocity gradient and stress BC
real(pReal), dimension(:), allocatable :: bc_timeIncrement ! length of increment real(pReal), dimension(:), allocatable :: bc_timeIncrement ! length of increment
integer(pInt) N_Loadcases, steps
integer(pInt), dimension(:), allocatable :: bc_steps ! number of steps
logical, dimension(:,:,:,:), allocatable :: bc_mask ! mask of boundary conditions logical, dimension(:,:,:,:), allocatable :: bc_mask ! mask of boundary conditions
! variables storing information from mesh file ! variables storing information from mesh file
integer(pInt) homog, prodnn
real(pReal) wgt real(pReal) wgt
integer(pInt), dimension(3) :: resolution
real(pReal), dimension(3) :: meshdimension real(pReal), dimension(3) :: meshdimension
integer(pInt) homog, prodnn
integer(pInt), dimension(3) :: resolution
! stress etc. ! stress etc.
real(pReal), dimension(3,3) :: pstress, defgradmacro, pstress_av, defgrad_av, temp33_Real
real(pReal), dimension(3,3,3,3) :: dPdF, c0, s0
real(pReal), dimension(6) :: cstress ! cauchy stress in Mandel notation (not needed) real(pReal), dimension(6) :: cstress ! cauchy stress in Mandel notation (not needed)
real(pReal), dimension(3,3) :: pstress ! Piola-Kirchhoff stress in Matrix notation real(pReal), dimension(6,6) :: dsde
real(pReal), dimension(3,3,3,3) :: dPdF, c0, s0 ! ??, reference stiffnes, compliance real(pReal), dimension(9,9) :: s099
real(pReal), dimension(6,6) :: dsde, s066
real(pReal), dimension(3,3) :: defgradmacro
real(pReal), dimension(3,3) :: pstress_av, defgrad_av, temp33_Real
real(pReal), dimension(:,:,:,:,:), allocatable :: pstress_field, defgrad, defgradold
real(pReal), dimension(:,:,:), allocatable :: ddefgrad real(pReal), dimension(:,:,:), allocatable :: ddefgrad
real(pReal), dimension(:,:,:,:,:), allocatable :: pstress_field, defgrad, defgradold
! variables storing information for spectral method ! variables storing information for spectral method
complex(pReal), dimension(:,:,:,:,:), allocatable :: workfft complex(pReal), dimension(:,:,:,:,:), allocatable :: workfft
complex(pReal), dimension(3,3) :: temp33_Complex complex(pReal), dimension(3,3) :: temp33_Complex
real(pReal), dimension(:,:,:,:,:,:,:), allocatable :: gamma_hat real(pReal), dimension(:,:,:,:,:,:,:), allocatable :: gamma_hat
real(pReal), dimension(:,:,:,:,:), allocatable :: xknormdyad real(pReal), dimension(:,:,:,:,:), allocatable :: xinormdyad
real(pReal), dimension(:,:,:,:), allocatable :: xi real(pReal), dimension(:,:,:,:), allocatable :: xi
integer(pInt), dimension(3) :: k_s integer(pInt), dimension(3) :: k_s
integer*8, dimension(2,3,3) :: plan_fft integer*8, dimension(2,3,3) :: plan_fft
! convergency etc. ! convergency etc.
logical errmatinv
integer(pInt) itmax, ierr
real(pReal) error, err_div, sigma0 real(pReal) error, err_div, sigma0
integer(pInt) itmax, ierr
logical errmatinv
! loop variables etc. ! loop variables etc.
real(pReal) guessmode ! flip-flop to guess defgrad fluctuation field evolution
integer(pInt) i, j, k, l, m, n, p integer(pInt) i, j, k, l, m, n, p
integer(pInt) loadcase, ielem, iter, calcmode integer(pInt) loadcase, ielem, iter, calcmode
real(pReal) guessmode ! flip-flop to guess defgrad fluctuation field evolution
real(pReal) temperature ! not used, but needed real(pReal) temperature ! not used, but needed
@ -115,7 +114,7 @@ program mpie_spectral
gotResolution =.false.; gotDimension =.false.; gotHomogenization = .false. gotResolution =.false.; gotDimension =.false.; gotHomogenization = .false.
if (IargC() < 2) call IO_error(102) ! check for correct number of arguments given if (IargC() /= 2) call IO_error(102) ! check for correct number of given arguments
! Reading the loadcase file and assign variables ! Reading the loadcase file and assign variables
path = getLoadcaseName() path = getLoadcaseName()
@ -251,7 +250,7 @@ program mpie_spectral
allocate (workfft(resolution(1)/2+1,resolution(2),resolution(3),3,3)); workfft = 0.0_pReal allocate (workfft(resolution(1)/2+1,resolution(2),resolution(3),3,3)); workfft = 0.0_pReal
allocate (gamma_hat(resolution(1)/2+1,resolution(2),resolution(3),3,3,3,3)); gamma_hat = 0.0_pReal allocate (gamma_hat(resolution(1)/2+1,resolution(2),resolution(3),3,3,3,3)); gamma_hat = 0.0_pReal
allocate (xknormdyad(resolution(1)/2+1,resolution(2),resolution(3),3,3)); xknormdyad = 0.0_pReal allocate (xinormdyad(resolution(1)/2+1,resolution(2),resolution(3),3,3)); xinormdyad = 0.0_pReal
allocate (xi(resolution(1)/2+1,resolution(2),resolution(3),3)); xi = 0.0_pReal allocate (xi(resolution(1)/2+1,resolution(2),resolution(3),3)); xi = 0.0_pReal
allocate (pstress_field(resolution(1),resolution(2),resolution(3),3,3)); pstress_field = 0.0_pReal allocate (pstress_field(resolution(1),resolution(2),resolution(3),3,3)); pstress_field = 0.0_pReal
allocate (displacement(resolution(1),resolution(2),resolution(3),3)); displacement = 0.0_pReal allocate (displacement(resolution(1),resolution(2),resolution(3),3)); displacement = 0.0_pReal
@ -259,6 +258,7 @@ program mpie_spectral
allocate (defgradold(resolution(1),resolution(2),resolution(3),3,3)); defgradold = 0.0_pReal allocate (defgradold(resolution(1),resolution(2),resolution(3),3,3)); defgradold = 0.0_pReal
allocate (ddefgrad(resolution(1),resolution(2),resolution(3))); ddefgrad = 0.0_pReal allocate (ddefgrad(resolution(1),resolution(2),resolution(3))); ddefgrad = 0.0_pReal
! Initialization of fftw (see manual on fftw.org for more details)
call dfftw_init_threads(ierr) call dfftw_init_threads(ierr)
call dfftw_plan_with_nthreads(4) call dfftw_plan_with_nthreads(4)
do m = 1,3; do n = 1,3 do m = 1,3; do n = 1,3
@ -272,6 +272,7 @@ program mpie_spectral
wgt = 1_pReal/real(prodnn, pReal) wgt = 1_pReal/real(prodnn, pReal)
defgradmacro = math_I3 defgradmacro = math_I3
! Initialization of CPFEM_general (= constitutive law) and of deformation gradient field
ielem = 0_pInt ielem = 0_pInt
do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1) do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1)
defgradold(i,j,k,:,:) = math_I3 !no deformation at the beginning defgradold(i,j,k,:,:) = math_I3 !no deformation at the beginning
@ -280,7 +281,7 @@ program mpie_spectral
call CPFEM_general(2,math_I3,math_I3,temperature,0.0_pReal,ielem,1_pInt,cstress,dsde,pstress,dPdF) call CPFEM_general(2,math_I3,math_I3,temperature,0.0_pReal,ielem,1_pInt,cstress,dsde,pstress,dPdF)
enddo; enddo; enddo enddo; enddo; enddo
!calculation of xknormdyad (needed to calculate gamma_hat) and xi (waves, needed for proof of equilibrium) !calculation of xinormdyad (needed to calculate gamma_hat) and xi (waves, needed for proof of equilibrium)
do k = 1, resolution(3) do k = 1, resolution(3)
k_s(3) = k-1 k_s(3) = k-1
if(k > resolution(3)/2+1) k_s(3) = k_s(3)-resolution(3) if(k > resolution(3)/2+1) k_s(3) = k_s(3)-resolution(3)
@ -289,14 +290,14 @@ program mpie_spectral
if(j > resolution(2)/2+1) k_s(2) = k_s(2)-resolution(2) if(j > resolution(2)/2+1) k_s(2) = k_s(2)-resolution(2)
do i = 1, resolution(1)/2+1 do i = 1, resolution(1)/2+1
k_s(1) = i-1 k_s(1) = i-1
xi(i,j,k,3) = .0_pReal xi(i,j,k,3) = 0.0_pReal
if(resolution(3) > 1) xi(i,j,k,3) = real(k_s(3), pReal)/meshdimension(3) if(resolution(3) > 1) xi(i,j,k,3) = real(k_s(3), pReal)/meshdimension(3)
xi(i,j,k,2) = real(k_s(2), pReal)/meshdimension(2) xi(i,j,k,2) = real(k_s(2), pReal)/meshdimension(2)
xi(i,j,k,1) = real(k_s(1), pReal)/meshdimension(1) xi(i,j,k,1) = real(k_s(1), pReal)/meshdimension(1)
if (any(xi(i,j,k,:) /= .0_pReal)) then if (any(xi(i,j,k,:) /= 0.0_pReal)) then
do l = 1,3; do m = 1,3 do l = 1,3; do m = 1,3
xknormdyad(i,j,k, l,m) = xi(i,j,k, l)*xi(i,j,k, m)/sum(xi(i,j,k,:)**2) xinormdyad(i,j,k, l,m) = xi(i,j,k, l)*xi(i,j,k, m)/sum(xi(i,j,k,:)**2)
enddo; enddo enddo; enddo
endif endif
@ -312,23 +313,23 @@ program mpie_spectral
!************************************************************* !*************************************************************
timeinc = bc_timeIncrement(loadcase)/bc_steps(loadcase) timeinc = bc_timeIncrement(loadcase)/bc_steps(loadcase)
guessmode = 0.0_pReal ! change of load case guessmode = 0.0_pReal ! change of load case, homogeneous guess for the first step
!************************************************************* !*************************************************************
! loop oper steps defined in input file for current loadcase ! loop oper steps defined in input file for current loadcase
do steps = 1, bc_steps(loadcase) do steps = 1, bc_steps(loadcase)
!************************************************************* !*************************************************************
defgradmacro = defgradmacro& defgradmacro = defgradmacro& ! update macroscopic displacement gradient (BC of defgrad)
+ math_mul33x33(bc_velocityGrad(:,:,loadcase), defgradmacro)*timeinc !update macroscopic displacement gradient (stores the desired BCs of defgrad) + math_mul33x33(bc_velocityGrad(:,:,loadcase), defgradmacro)*timeinc
do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1) do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1)
temp33_Real = defgrad(i,j,k,:,:) temp33_Real = defgrad(i,j,k,:,:)
defgrad(i,j,k,:,:) = defgrad(i,j,k,:,:)& defgrad(i,j,k,:,:) = defgrad(i,j,k,:,:)& ! old fluctuations as guess for new step, no fluctuations for new loadcase
+ guessmode * (defgrad(i,j,k,:,:) - defgradold(i,j,k,:,:))& ! old fluctuations as guess for new step + guessmode * (defgrad(i,j,k,:,:) - defgradold(i,j,k,:,:))&
+ (1.0_pReal-guessmode) * math_mul33x33(bc_velocityGrad(:,:,loadcase),defgradold(i,j,k,:,:))*timeinc ! no fluctuations for new loadcase + (1.0_pReal-guessmode) * math_mul33x33(bc_velocityGrad(:,:,loadcase),defgradold(i,j,k,:,:))*timeinc
defgradold(i,j,k,:,:) = temp33_Real defgradold(i,j,k,:,:) = temp33_Real
enddo; enddo; enddo enddo; enddo; enddo
guessmode = 1_pReal ! keep guessing along former trajectory guessmode = 1_pReal ! keep guessing along former trajectory during same loadcase
calcmode = 1_pInt calcmode = 1_pInt
iter = 0_pInt iter = 0_pInt
err_div= 2_pInt * error err_div= 2_pInt * error
@ -339,9 +340,10 @@ program mpie_spectral
iter = iter + 1 iter = iter + 1
print '(A,I5.5,tr2,A,I5.5)', ' Step = ',steps,'Iteration = ',iter print '(A,I5.5,tr2,A,I5.5)', ' Step = ',steps,'Iteration = ',iter
!************************************************************* !*************************************************************
err_div = .0_pReal; sigma0 = .0_pReal err_div = 0.0_pReal; sigma0 = 0.0_pReal
pstress_av = .0_pReal; defgrad_av=.0_pReal pstress_av = 0.0_pReal; defgrad_av = 0.0_pReal
! Calculate stress field for current deformation gradient using CPFEM_general
print *, 'Update Stress Field' print *, 'Update Stress Field'
ielem = 0_pInt ielem = 0_pInt
do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1) do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1)
@ -351,7 +353,7 @@ program mpie_spectral
cstress,dsde, pstress, dPdF) cstress,dsde, pstress, dPdF)
enddo; enddo; enddo enddo; enddo; enddo
c0 = .0_pReal c0 = 0.0_pReal
ielem = 0_pInt ielem = 0_pInt
do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1) do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1)
ielem = ielem + 1 ielem = ielem + 1
@ -364,34 +366,35 @@ program mpie_spectral
pstress_field(i,j,k,:,:) = pstress pstress_field(i,j,k,:,:) = pstress
pstress_av = pstress_av + pstress ! average stress pstress_av = pstress_av + pstress ! average stress
enddo; enddo; enddo enddo; enddo; enddo
pstress_av = pstress_av*wgt ! do the weighting of average stress pstress_av = pstress_av*wgt ! do the weighting of average stress
if(iter==1) then !update gamma_hat with new reference stiffness ! Update gamma_hat with new reference stiffness and calculate new average compliance (for stress BC)
call math_invert(6,math_mandel3333to66(c0),s066,i,errmatinv) !i is just a dummy variable if(iter==1) then
call math_invert(9, math_plain3333to99(c0),s099,i, errmatinv)
if(errmatinv) call IO_error(45,ext_msg = "problem in c0 inversion") ! todo: change number and add message to io.f90 (and remove No. 48) if(errmatinv) call IO_error(45,ext_msg = "problem in c0 inversion") ! todo: change number and add message to io.f90 (and remove No. 48)
s0 = math_mandel66to3333(s066)*real(prodnn, pReal) s0 = math_plain99to3333(s099) * real(prodnn, pReal)
c0 = c0 * wgt c0 = c0 * wgt
do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1)/2+1 do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1)/2+1
temp33_Real = .0_pReal temp33_Real = 0.0_pReal
do l = 1,3; do m = 1,3; do n = 1,3; do p = 1,3 do l = 1,3; do m = 1,3; do n = 1,3; do p = 1,3
temp33_Real(l,m) = temp33_Real(l,m)+c0(l,n,m,p)*xknormdyad(i,j,k, n,p) temp33_Real(l,m) = temp33_Real(l,m) + c0(l,n,m,p) * xinormdyad(i,j,k, n,p)
enddo; enddo; enddo; enddo enddo; enddo; enddo; enddo
temp33_Real = math_inv3x3(temp33_Real) temp33_Real = math_inv3x3(temp33_Real)
do l=1,3; do m=1,3; do n=1,3; do p=1,3 do l=1,3; do m=1,3; do n=1,3; do p=1,3
gamma_hat(i,j,k, l,m,n,p) = -temp33_Real(l,n)*xknormdyad(i,j,k, m,p) gamma_hat(i,j,k, l,m,n,p) = - temp33_Real(l,n) * xinormdyad(i,j,k, m,p)
enddo; enddo; enddo; enddo enddo; enddo; enddo; enddo
enddo; enddo; enddo enddo; enddo; enddo
endif endif
! Using the spectral method to calculate the change of deformation gradient, check divergence of stress field in fourier space
print *, 'Update Deformation Gradient Field' print *, 'Update Deformation Gradient Field'
do m = 1,3; do n = 1,3 do m = 1,3; do n = 1,3
call dfftw_execute_dft_r2c(plan_fft(1,m,n), pstress_field(:,:,:,m,n),workfft(:,:,:,m,n)) call dfftw_execute_dft_r2c(plan_fft(1,m,n), pstress_field(:,:,:,m,n),workfft(:,:,:,m,n))
if(n == 3) sigma0 = max(sigma0, sum(abs(real(workfft(1,1,1,m,:))))) ! L infinity Norm of stress tensor in Fourier space if(n == 3) sigma0 = max(sigma0, sum(abs(real(workfft(1,1,1,m,:))))) ! L infinity Norm of stress tensor
enddo; enddo enddo; enddo
do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1)/2+1 do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1)/2+1
err_div = err_div + (maxval(abs(math_mul33x3c(workfft(i,j,k,:,:),xi(i,j,k,:))))) ! L infinity Norm of div(stress tensor) in Fourier space err_div = err_div + (maxval(abs(math_mul33x3_complex(workfft(i,j,k,:,:),xi(i,j,k,:))))) ! L infinity Norm of div(stress)
temp33_Complex = .0_pReal temp33_Complex = .0_pReal
do m = 1,3; do n = 1,3 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,j,k,:,:))
@ -399,7 +402,7 @@ program mpie_spectral
workfft(i,j,k,:,:) = temp33_Complex(:,:) workfft(i,j,k,:,:) = temp33_Complex(:,:)
enddo; enddo; enddo enddo; enddo; enddo
err_div = err_div/(real(prodnn/resolution(1)*(resolution(1)/2+1)))/sigma0 !calculate error (divergence of stress field) err_div = err_div/(real(prodnn/resolution(1)*(resolution(1)/2+1)))/sigma0 !weighting of error
do m = 1,3; do n = 1,3 do m = 1,3; do n = 1,3
call dfftw_execute_dft_c2r(plan_fft(2,m,n), workfft(:,:,:,m,n),ddefgrad(:,:,:)) call dfftw_execute_dft_c2r(plan_fft(2,m,n), workfft(:,:,:,m,n),ddefgrad(:,:,:))
@ -422,7 +425,7 @@ program mpie_spectral
enddo; enddo enddo; enddo
print '(2(a,E8.2))', ' Error = ',err_div,' Criteria = ', error print '(2(a,E8.2))', ' Error = ',err_div,' Criteria = ', error
print '(A)', '----------------------------------' print '(A)', '---------------------------------------'
enddo ! end looping when convergency is achieved enddo ! end looping when convergency is achieved
@ -436,9 +439,9 @@ program mpie_spectral
print '(A,3(E10.4,tr2))', ' ', pstress_av(3,:) print '(A,3(E10.4,tr2))', ' ', pstress_av(3,:)
print '(A)', '************************************************************' print '(A)', '************************************************************'
!gsmh output ! Postprocessing (gsmh output)
temp33_Real(1,:) = 0.0_pReal temp33_Real(1,:) = 0.0_pReal; temp33_Real(1,3) = -(real(resolution(3))/meshdimension(3)) ! start just below origin
temp33_Real(1,3) = -1.0_pReal
do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1) do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1)
if((j==1).and.(i==1)) then if((j==1).and.(i==1)) then
temp33_Real(1,:) = temp33_Real(1,:) + math_mul33x3(defgrad(i,j,k,:,:),(/0.0_pReal,0.0_pReal,(real(resolution(3))/meshdimension(3))/)) temp33_Real(1,:) = temp33_Real(1,:) + math_mul33x3(defgrad(i,j,k,:,:),(/0.0_pReal,0.0_pReal,(real(resolution(3))/meshdimension(3))/))
@ -456,59 +459,45 @@ program mpie_spectral
endif endif
endif endif
enddo; enddo; enddo enddo; enddo; enddo
write(nriter, *) iter
write(nrstep, *) steps write(nriter, *) iter; write(nrstep, *) steps
nrstep = 'stress'//trim(adjustl(nrstep))//'-'//trim(adjustl(nriter))//'_cpfem.msh' open(589,file = 'stress' //trim(adjustl(nrstep))//'-'//trim(adjustl(nriter))//'_cpfem.msh')
open(589,file = nrstep) open(588,file = 'disgrad'//trim(adjustl(nrstep))//'-'//trim(adjustl(nriter))//'_cpfem.msh')
write(589, '(A, /, A, /, A, /, A, /, I10)'), '$MeshFormat', '2.1 0 8', '$EndMeshFormat', '$Nodes', prodnn write(589, '(4(A, /), I10)'), '$MeshFormat', '2.1 0 8', '$EndMeshFormat', '$Nodes', prodnn
write(588, '(4(A, /), I10)'), '$MeshFormat', '2.1 0 8', '$EndMeshFormat', '$Nodes', prodnn
ielem = 0_pInt ielem = 0_pInt
do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1) do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1)
ielem = ielem + 1 ielem = ielem + 1
write(589, '(I10,tr2,E12.6,tr2,E12.6,tr2,E12.6)'), ielem, displacement(i,j,k,:) !real(i), real(j), real(k) write(589, '(I10, 3(tr2, E12.6))'), ielem, displacement(i,j,k,:)
write(588, '(I10, 3(tr2, E12.6))'), ielem, displacement(i,j,k,:)
enddo; enddo; enddo enddo; enddo; enddo
write(589, '(A, /, A, /, I10)'), '$EndNodes', '$Elements', prodnn
write(589, '(2(A, /), I10)'), '$EndNodes', '$Elements', prodnn
write(588, '(2(A, /), I10)'), '$EndNodes', '$Elements', prodnn
do i = 1, prodnn do i = 1, prodnn
write(589, '(I10, A, I10)'), i, ' 15 2 1 2', i write(589, '(I10, A, I10)'), i, ' 15 2 1 2', i
write(588, '(I10, A, I10)'), i, ' 15 2 1 2', i
enddo enddo
write(589, '(A)'), '$EndElements'
write(589, '(A, /, A, /, A, /, A, /, A, /, A, /, A, /, A, /, I10)'), '$NodeData', '1',&
'"'//trim(adjustl(nrstep))//'"', '1','0.0', '3', '0', '9', prodnn
ielem = 0_pInt
do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1)
ielem = ielem + 1
write(589, '(i10,tr2,E12.6,tr2,E12.6,tr2,E12.6,tr2,E12.6,tr2,E12.6,tr2,E12.6,&
tr2,E12.6,tr2,E12.6,tr2,E12.6,tr2)'), ielem, pstress_field(i,j,k,:,:)
enddo; enddo; enddo
write(nriter, *) iter
write(nrstep, *) steps
write(589, *), '$EndNodeData'
close(589)
nrstep = 'defgrad'//trim(adjustl(nrstep))//'-'//trim(adjustl(nriter))//'_cpfem.msh'
open(589,file = nrstep)
write(589, '(A, /, A, /, A, /, A, /, I10)'), '$MeshFormat', '2.1 0 8', '$EndMeshFormat', '$Nodes', prodnn
ielem = 0_pInt
do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1)
ielem = ielem + 1
write(589, '(I10,tr2,E12.6,tr2,E12.6,tr2,E12.6)'), ielem, displacement(i,j,k,:) !real(i), real(j), real(k)
enddo; enddo; enddo
write(589, '(A, /, A, /, I10)'), '$EndNodes', '$Elements', prodnn
do i = 1, prodnn
write(589, '(I10, A, I10)'), i, ' 15 2 1 2', i
enddo
write(589, '(A)'), '$EndElements' write(589, '(A)'), '$EndElements'
write(589, '(A, /, A, /, A, /, A, /, A, /, A, /, A, /, A, /, I10)'), '$NodeData', '1',& write(588, '(A)'), '$EndElements'
'"'//trim(adjustl(nrstep))//'"', '1','0.0', '3', '0', '9', prodnn write(589, '(8(A, /), I10)'), '$NodeData', '1','"'//trim(adjustl('stress'//trim(adjustl(nrstep))//&
'-'//trim(adjustl(nriter))//'_cpfem.msh'))//'"','1','0.0', '3', '0', '9', prodnn
write(588, '(8(A, /), I10)'), '$NodeData', '1','"'//trim(adjustl('disgrad'//trim(adjustl(nrstep))//&
'-'//trim(adjustl(nriter))//'_cpfem.msh'))//'"','1','0.0', '3', '0', '9', prodnn
ielem = 0_pInt ielem = 0_pInt
do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1) do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1)
ielem = ielem + 1 ielem = ielem + 1
write(589, '(i10,tr2,E12.6,tr2,E12.6,tr2,E12.6,tr2,E12.6,tr2,E12.6,tr2,E12.6,& write(589, '(i10, 9(tr2, E12.6))'), ielem, pstress_field(i,j,k,:,:)
tr2,E12.6,tr2,E12.6,tr2,E12.6,tr2)'), ielem, defgrad(i,j,k,:,:) - math_I3 write(588, '(i10, 9(tr2, E12.6))'), ielem, defgrad(i,j,k,:,:) - math_I3
enddo; enddo; enddo enddo; enddo; enddo
write(589, *), '$EndNodeData' write(589, *), '$EndNodeData'
close(589) write(588, *), '$EndNodeData'
!end gmsh close(589); close(588)
enddo ! end looping over steps in current loadcase enddo ! end looping over steps in current loadcase
enddo ! end looping over loadcases enddo ! end looping over loadcases