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

changed declaration of two variables from real to int, spectral method is now working. 

started to commend Ricardos code, layoutet loops and removed some redundant variables.
until now, no error calculation is done. at the moment calculations are in an infinite loop
This commit is contained in:
Martin Diehl 2010-07-02 14:10:36 +00:00
parent 661bb97800
commit d6ba9d54b6
1 changed files with 395 additions and 489 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

652
code/mpie_spectral.f90
View File

@ -44,7 +44,7 @@ subroutine mpie_interface_init()
write(6,*) write(6,*)
return return
end subroutine endsubroutine
!******************************************************************** !********************************************************************
! extract working directory from loadcase file ! extract working directory from loadcase file
@ -70,7 +70,7 @@ function getSolverWorkingDirectoryName()
return return
end function endfunction
!******************************************************************** !********************************************************************
! basename of meshfile from command line arguments ! basename of meshfile from command line arguments
@ -105,7 +105,7 @@ function getSolverJobName()
getSolverJobName = makeRelativePath(getSolverWorkingDirectoryName(),& getSolverJobName = makeRelativePath(getSolverWorkingDirectoryName(),&
getSolverJobName) getSolverJobName)
return return
end function endfunction
!******************************************************************** !********************************************************************
@ -121,7 +121,7 @@ function getLoadcaseName()
character(len=1024) getLoadcaseName, outName, cwd character(len=1024) getLoadcaseName, outName, cwd
character(len=*), parameter :: pathSep = achar(47)//achar(92) ! /, \ character(len=*), parameter :: pathSep = achar(47)//achar(92) ! /, \
integer(pInt) posExt,posSep integer(pInt) posExt,posSep
posExt=0 !not sure if its needed posExt = 0 !not sure if its needed
call getarg(2,getLoadcaseName) call getarg(2,getLoadcaseName)
posExt = scan(getLoadcaseName,'.',back=.true.) posExt = scan(getLoadcaseName,'.',back=.true.)
@ -138,7 +138,7 @@ function getLoadcaseName()
getLoadcaseName = makeRelativePath(getSolverWorkingDirectoryName(),& getLoadcaseName = makeRelativePath(getSolverWorkingDirectoryName(),&
getLoadcaseName) getLoadcaseName)
return return
end function endfunction
!******************************************************************** !********************************************************************
@ -158,10 +158,10 @@ function rectifyPath(path)
!remove ./ from path !remove ./ from path
l = len_trim(path) l = len_trim(path)
rectifyPath = path rectifyPath = path
do i=l,2,-1 do i = l,2,-1
if ( rectifyPath(i-1:i)=='./' .and. rectifyPath(i-2:i-2) /= '.' ) & if ( rectifyPath(i-1:i) == './' .and. rectifyPath(i-2:i-2) /= '.' ) &
rectifyPath(i-1:l) = rectifyPath(i+1:l)//' ' rectifyPath(i-1:l) = rectifyPath(i+1:l)//' '
end do enddo
!remove ../ and corresponding directory from rectifyPath !remove ../ and corresponding directory from rectifyPath
l = len_trim(rectifyPath) l = len_trim(rectifyPath)
@ -173,11 +173,11 @@ function rectifyPath(path)
j = scan(rectifyPath(:i-2),'/',back=.true.) j = scan(rectifyPath(:i-2),'/',back=.true.)
rectifyPath(j+1:l) = rectifyPath(i+3:l)//repeat(' ',2+i-j) rectifyPath(j+1:l) = rectifyPath(i+3:l)//repeat(' ',2+i-j)
i = j+index(rectifyPath(j+1:l),'../') i = j+index(rectifyPath(j+1:l),'../')
end do enddo
if(len_trim(rectifyPath)==0) rectifyPath='/' if(len_trim(rectifyPath) == 0) rectifyPath = '/'
return return
end function rectifyPath endfunction rectifyPath
!******************************************************************** !********************************************************************
@ -205,7 +205,7 @@ function makeRelativePath(a,b)
enddo enddo
makeRelativePath = repeat('../',remainingSlashes)//b(posLastCommonSlash+1:len_trim(b)) makeRelativePath = repeat('../',remainingSlashes)//b(posLastCommonSlash+1:len_trim(b))
return return
end function makeRelativePath endfunction makeRelativePath
END MODULE END MODULE
@ -269,18 +269,11 @@ program mpie_spectral
integer(pInt) unit, N_l, N_s, N_t, N_n, N, i, j, k, l ! numbers of identifiers, loop variables integer(pInt) unit, N_l, N_s, N_t, N_n, N, i, j, k, l ! numbers of identifiers, loop variables
integer(pInt) a, b, c, e, homog integer(pInt) a, b, c, e, homog
real(pReal) x, y, z real(pReal) x, y, z
!------------------------- !-------------------------
!begin RL !begin RL
!------------------------- !-------------------------
! begin change by m.diehl
!integer(pInt), parameter :: npts1 = 32
!integer(pInt), parameter :: npts2 = 32
!integer(pInt), parameter :: npts3 = 32
integer(pInt) npts1, npts2, npts3 ! will be allocated later
! end change by m.diehl
real(pReal), dimension(:), allocatable :: datafft real(pReal), dimension(:), allocatable :: datafft
real(pReal), dimension(:,:,:,:,:), allocatable :: workfft,workfftim,sg,disgrad,defgradold real(pReal), dimension(:,:,:,:,:), allocatable :: workfft,workfftim,sg,disgrad,defgradold
@ -291,13 +284,15 @@ program mpie_spectral
real(pReal), dimension(3,3) :: defgrad0,defgrad real(pReal), dimension(3,3) :: defgrad0,defgrad
real(pReal), dimension(3,3) :: udot,scauchy,scauav,aux33 real(pReal), dimension(3,3) :: udot,scauchy,scauav,aux33
real(pReal), dimension(3) :: delt,nn,xk,xk2 integer(pInt), dimension(3) :: nn
real(pReal), dimension(2) :: nn2 integer(pInt), dimension(2) :: nn2
real(pReal), dimension(3) :: delt,xk,xk2
real(pReal), dimension(6) :: aux6 real(pReal), dimension(6) :: aux6
real(pReal), dimension(3,3,3,3) :: c0,s0,g1 real(pReal), dimension(3,3,3,3) :: c0,s0,g1
real(pReal), dimension(6,6) :: c066,s066 real(pReal), dimension(6,6) :: c066,s066
integer(pInt) itmax,nsteps,jload,ielem,ii,jj,k1,kxx,kyy,kzz,kx,ky,kz,idum,iter,imicro,m1,n1,p,q integer(pInt) itmax, jload, ielem, ii, jj, k1, kxx, kyy, kzz, kx, ky, kz, idum, iter, imicro, m1, n1, p, q
real(pReal) prodnn,wgt,error,tdot,erre,errs,evm,svm,det,xknorm real(pReal) prodnn,wgt,error,tdot,erre,errs,evm,svm,det,xknorm
@ -320,11 +315,11 @@ program mpie_spectral
print*,'Loadcase: ',trim(path) print*,'Loadcase: ',trim(path)
print*,'Workingdir: ',trim(getSolverWorkingDirectoryName()) print*,'Workingdir: ',trim(getSolverWorkingDirectoryName())
if (.not. IO_open_file(unit,path)) call IO_error(45,ext_msg=path) if (.not. IO_open_file(unit,path)) call IO_error(45,ext_msg = path)
rewind(unit) rewind(unit)
do do
read(unit,'(a1024)',END=101) line read(unit,'(a1024)',END = 101) line
if (IO_isBlank(line)) cycle ! skip empty lines if (IO_isBlank(line)) cycle ! skip empty lines
posInput = IO_stringPos(line,maxNchunksInput) posInput = IO_stringPos(line,maxNchunksInput)
do i = 1,maxNchunksInput,1 do i = 1,maxNchunksInput,1
@ -340,7 +335,7 @@ program mpie_spectral
end select end select
enddo ! count all identifiers to allocate memory and do sanity check enddo ! count all identifiers to allocate memory and do sanity check
if ((N_l /= N_s).or.(N_s /= N_t).or.(N_t /= N_n)) & ! sanity check if ((N_l /= N_s).or.(N_s /= N_t).or.(N_t /= N_n)) & ! sanity check
call IO_error(46,ext_msg=path) ! error message for incomplete input file call IO_error(46,ext_msg = path) ! error message for incomplete input file
enddo enddo
@ -355,9 +350,9 @@ program mpie_spectral
rewind(unit) rewind(unit)
j = 0_pInt j = 0_pInt
do do
read(unit,'(a1024)',END=200) line read(unit,'(a1024)',END = 200) line
if (IO_isBlank(line)) cycle ! skip empty lines if (IO_isBlank(line)) cycle ! skip empty lines
j=j+1 j = j+1
posInput = IO_stringPos(line,maxNchunksInput) posInput = IO_stringPos(line,maxNchunksInput)
do i = 1,maxNchunksInput,2 do i = 1,maxNchunksInput,2
select case (IO_lc(IO_stringValue(line,posInput,i))) select case (IO_lc(IO_stringValue(line,posInput,i)))
@ -414,11 +409,11 @@ program mpie_spectral
gotHomogenization = .false. gotHomogenization = .false.
path = getSolverJobName() path = getSolverJobName()
print*,'JobName: ',trim(path) print*,'JobName: ',trim(path)
if (.not. IO_open_file(unit,trim(path)//InputFileExtension)) call IO_error(101,ext_msg=path) if (.not. IO_open_file(unit,trim(path)//InputFileExtension)) call IO_error(101,ext_msg = path)
rewind(unit) rewind(unit)
do do
read(unit,'(a1024)',END=100) line read(unit,'(a1024)',END = 100) line
if (IO_isBlank(line)) cycle ! skip empty lines if (IO_isBlank(line)) cycle ! skip empty lines
posMesh = IO_stringPos(line,maxNchunksMesh) posMesh = IO_stringPos(line,maxNchunksMesh)
@ -461,267 +456,212 @@ program mpie_spectral
print *, '' print *, ''
temperature = 300.0_pReal temperature = 300.0_pReal
!-------------------------
!RL
!------------------------- !-------------------------
!begin change m.diehl !begin RL
npts1=a !-------------------------
npts2=b allocate (datafft(2*a*b*c))
npts3=c allocate (workfft(3,3,a,b,c))
!end change m.diehl allocate (workfftim(3,3,a,b,c))
allocate (sg(3,3,a,b,c))
allocate (disgrad(3,3,a,b,c))
allocate (defgradold(3,3,a,b,c))
allocate (datafft(2*npts1*npts2*npts3)) error = 0.00001
allocate (workfft(3,3,npts1,npts2,npts3)) itmax = 100
allocate (workfftim(3,3,npts1,npts2,npts3))
allocate (sg(3,3,npts1,npts2,npts3))
allocate (disgrad(3,3,npts1,npts2,npts3))
allocate (defgradold(3,3,npts1,npts2,npts3))
error=0.00001 delt(1) = 1.
itmax=100 delt(2) = 1.
delt(3) = 1.
delt(1)=1. nn(1) = a
delt(2)=1. nn(2) = b
delt(3)=1. nn(3) = c
nn(1)=npts1 nn2(1) = a
nn(2)=npts2 nn2(2) = b
nn(3)=npts3
nn2(1)=npts1 prodnn = nn(1)*nn(2)*nn(3)
nn2(2)=npts2 wgt = 1./prodnn
prodnn=nn(1)*nn(2)*nn(3)
wgt=1./prodnn
! C_0 and S_0 CALCULATION ! C_0 and S_0 CALCULATION
!! !!
!! PHILIP: FE_exec_elem? !! PHILIP: FE_exec_elem?
!! !!
stress = 0. !should be initialized somewhere
dsde = 0. !should be initialized somewhere
c0 = 0. !stiffness of reference material
c066 = 0. ! other way of notating c0
ielem=0 do ielem = 1, prodnn !call each element with identity (math_i3) to initialize with high stress
c0=0.
c066=0.
do k=1,npts3
do j=1,npts2
do i=1,npts1
ielem=ielem+1
call CPFEM_general(3,math_i3,math_i3,temperature,0.0_pReal,ielem,1_pInt,stress,dsde) call CPFEM_general(3,math_i3,math_i3,temperature,0.0_pReal,ielem,1_pInt,stress,dsde)
c066 = c066+dsde
c066=c066+dsde c0 = c0+math_Mandel66to3333(dsde)
c0=c0+math_Mandel66to3333(dsde)
enddo
enddo
enddo enddo
c066=c066*wgt c066 = c066*wgt
c0=c0*wgt c0 = c0*wgt
call math_invert(6,c066,s066,idum,errmatinv) call math_invert(6,c066,s066,idum,errmatinv)
if(errmatinv) then if(errmatinv) then
write(*,*) 'ERROR IN C0 INVERSION' write(*,*) 'ERROR IN C0 INVERSION'
stop stop
endif endif
s0=math_Mandel66to3333(s066) s0 = math_Mandel66to3333(s066)
! INITIALIZATION BEFORE STARTING WITH LOADINGS ! INITIALIZATION BEFORE STARTING WITH LOADINGS
disgradmacro=0. disgradmacro = 0.
do k=1,npts3
do j=1,npts2
do i=1,npts1
defgradold(:,:,i,j,k)=math_i3(:,:)
do k = 1, c
do j = 1, b
do i = 1, a
defgradold(:,:,i,j,k) = math_i3(:,:)
enddo enddo
enddo enddo
enddo enddo
do jload = 1,N do jload = 1,N !Loop over loadcases defined in the loadcase file
udot(:,:) = bc_velocityGrad(:,:,jload)
scauchy(:,:) = bc_stress(:,:,jload)
iudot = 0
iscau = 0
udot(:,:)=bc_velocityGrad(:,:,jload) do i = 1, 3 !convert information about rb's from bc_mask in corresponding arrays
scauchy(:,:)=bc_stress(:,:,jload) do j = 1, 3
if(bc_mask(i,j,1,jload)) iudot(i,j) = 1
iudot=0 if(bc_mask(i,j,2,jload)) iscau(i,j) = 1
iscau=0
do i=1,3
do j=1,3
if(bc_mask(i,j,1,jload)) iudot(i,j)=1
if(bc_mask(i,j,2,jload)) iscau(i,j)=1
enddo enddo
enddo enddo
nsteps=bc_steps(jload) tdot = bc_timeIncrement(jload)/bc_steps(jload)
tdot=bc_timeIncrement(jload)/bc_steps(jload) ! evm = dvm*tdot ?
do imicro = 1, bc_steps(jload) ! loop oper steps defined in input file for current loadcase
! evm=dvm*tdot ?
do imicro=1,nsteps
write(*,*) '***************************************************' write(*,*) '***************************************************'
write(*,*) 'STEP = ',imicro write(*,*) 'STEP = ',imicro
! write(21,*) 'STEP = ',imicro
! INITIALIZATION BEFORE NEW TIME STEP ! INITIALIZATION BEFORE NEW TIME STEP
disgradmacro = disgradmacro+udot*tdot !update macroscopic displacementgradient
disgradmacro=disgradmacro+udot*tdot ddisgradmacro = 0.
ddisgradmacro=0.
ielem=0 ielem=0
do k=1,npts3 do k = 1, c !loop over FPs
do j=1,npts2 do j = 1, b
do i=1,npts1 do i = 1, a
ielem = ielem+1
ielem=ielem+1
disgrad(:,:,i,j,k)=disgradmacro(:,:) disgrad(:,:,i,j,k)=disgradmacro(:,:)
defgrad0(:,:)=defgradold(:,:,i,j,k) defgrad0(:,:)=defgradold(:,:,i,j,k)
defgrad(:,:)=math_i3(:,:)+disgrad(:,:,i,j,k) defgrad(:,:)=math_i3(:,:)+disgrad(:,:,i,j,k) ! calculate new defgrad
call CPFEM_general(3,defgrad0,defgrad,temperature,0.0_pReal,ielem,1_pInt,stress,dsde) call CPFEM_general(3,defgrad0,defgrad,temperature,0.0_pReal,ielem,1_pInt,stress,dsde)
! sg(:,:,i,j,k)=math_Mandel6to33(stress) ? ! sg(:,:,i,j,k)=math_Mandel6to33(stress) ?
enddo enddo
enddo enddo
enddo enddo
ielem=0 ielem = 0
do k=1,npts3 do k = 1, c !loop over FPs
do j=1,npts2 do j = 1, b
do i=1,npts1 do i = 1, a
ielem = ielem+1
ielem=ielem+1 disgrad(:,:,i,j,k) = disgradmacro(:,:)
disgrad(:,:,i,j,k)=disgradmacro(:,:)
defgrad0(:,:)=defgradold(:,:,i,j,k)
defgrad(:,:)=math_i3(:,:)+disgrad(:,:,i,j,k)
defgrad0(:,:) = defgradold(:,:,i,j,k)
defgrad(:,:) = math_i3(:,:)+disgrad(:,:,i,j,k)
call CPFEM_general(1,defgrad0,defgrad,temperature,0.0_pReal,ielem,1_pInt,stress,dsde) call CPFEM_general(1,defgrad0,defgrad,temperature,0.0_pReal,ielem,1_pInt,stress,dsde)
sg(:,:,i,j,k) = math_Mandel6to33(stress)
sg(:,:,i,j,k)=math_Mandel6to33(stress)
enddo enddo
enddo enddo
enddo enddo
ddisgradmacroacum=0. ddisgradmacroacum = 0.
iter=0 iter = 0
erre=2.*error erre = 2.*error
errs=2.*error errs = 2.*error
do while(iter.lt.itmax.and.(errs.gt.error.or.erre.gt.error))
iter=iter+1
do while(iter <= itmax.and.(errs.gt.error.or.erre.gt.error))
iter = iter+1
write(*,*)'ITER = ',iter write(*,*)'ITER = ',iter
write(*,*) 'DIRECT FFT OF STRESS FIELD' write(*,*) 'DIRECT FFT OF STRESS FIELD'
do ii = 1,3
do ii=1,3 do jj = 1,3
do jj=1,3 k1 = 0
do k = 1,c
k1=0
do k=1,npts3
! write(*,'(1H+,a,i2,2(a,i4))') ! write(*,'(1H+,a,i2,2(a,i4))')
! # 'STRESS - COMPONENT',ii,jj,' - Z = ',k,' OUT OF ',npts ! # 'STRESS - COMPONENT',ii,jj,' - Z = ',k,' OUT OF ',npts
do j = 1,b
do j=1,npts2 do i = 1,a
do i=1,npts1 k1 = k1+1
datafft(k1) = sg(ii,jj,i,j,k)
k1=k1+1 k1 = k1+1
datafft(k1)=sg(ii,jj,i,j,k) datafft(k1) = 0.
k1=k1+1
datafft(k1)=0.
enddo enddo
enddo enddo
enddo enddo
if(c.gt.1) then
if(npts3.gt.1) then
call fourn(datafft,nn,3,1) call fourn(datafft,nn,3,1)
else else
call fourn(datafft,nn2,2,1) call fourn(datafft,nn2,2,1)
endif endif
k1 = 0
k1=0 do k = 1, c
do j = 1, b
do k=1,npts3 do i = 1, a
do j=1,npts2 k1 = k1+1
do i=1,npts1 workfft(ii,jj,i,j,k) = datafft(k1)
k1 = k1+1
k1=k1+1 workfftim(ii,jj,i,j,k) = datafft(k1)
workfft(ii,jj,i,j,k)=datafft(k1)
k1=k1+1
workfftim(ii,jj,i,j,k)=datafft(k1)
enddo enddo
enddo enddo
enddo enddo
enddo enddo
enddo enddo
write(*,*) 'CALCULATING G^pq,ij : SG^ij ...' write(*,*) 'CALCULATING G^pq,ij : SG^ij ...'
do kzz = 1,c
do kyy = 1,b
do kxx = 1,a
if(kxx.le.a/2) kx = kxx-1
if(kxx.gt.a/2) kx = kxx-a-1
do kzz=1,npts3 if(kyy.le.b/2) ky = kyy-1
do kyy=1,npts2 if(kyy.gt.b/2) ky = kyy-b-1
do kxx=1,npts1
if(kxx.le.npts1/2) kx=kxx-1 if(kzz.le.c/2) kz = kzz-1
if(kxx.gt.npts1/2) kx=kxx-npts1-1 if(kzz.gt.c/2) kz = kzz-c-1
if(kyy.le.npts2/2) ky=kyy-1 xk(1) = kx/(delt(1)*nn(1))
if(kyy.gt.npts2/2) ky=kyy-npts2-1 xk(2) = ky/(delt(2)*nn(2))
if(kzz.le.npts3/2) kz=kzz-1 if(c.gt.1) then
if(kzz.gt.npts3/2) kz=kzz-npts3-1 xk(3) = kz/(delt(3)*nn(3))
xk(1)=kx/(delt(1)*nn(1))
xk(2)=ky/(delt(2)*nn(2))
if(npts3.gt.1) then
xk(3)=kz/(delt(3)*nn(3))
else else
xk(3)=0. xk(3) = 0.
endif endif
xknorm=sqrt(xk(1)**2+xk(2)**2+xk(3)**2) xknorm = sqrt(xk(1)**2+xk(2)**2+xk(3)**2)
if (xknorm.ne.0.) then if (xknorm.ne.0.) then
do i=1,3 do i = 1,3
xk2(i)=xk(i)/(xknorm*xknorm*2.*pi) xk2(i) = xk(i)/(xknorm*xknorm*2.*pi)
xk(i)=xk(i)/xknorm xk(i) = xk(i)/xknorm
enddo enddo
endif endif
do i=1,3 do i = 1,3
do k=1,3 do k = 1,3
aux33(i,k)=0. aux33(i,k) = 0.
do j=1,3 do j = 1,3
do l=1,3 do l = 1,3
aux33(i,k)=aux33(i,k)+c0(i,j,k,l)*xk(j)*xk(l) aux33(i,k) = aux33(i,k)+c0(i,j,k,l)*xk(j)*xk(l)
enddo enddo
enddo enddo
enddo enddo
enddo enddo
aux33 = math_inv3x3(aux33)
aux33=math_inv3x3(aux33)
! call minv(aux33,3,det,minv1,minv2) ! call minv(aux33,3,det,minv1,minv2)
! if(det.eq.0) then ! if(det.eq.0) then
@ -730,41 +670,34 @@ program mpie_spectral
! pause ! pause
! endif ! endif
do p=1,3 do p = 1,3
do q=1,3 do q = 1,3
do i=1,3 do i = 1,3
do j=1,3 do j = 1,3
g1(p,q,i,j)=-aux33(p,i)*xk(q)*xk(j) g1(p,q,i,j) = -aux33(p,i)*xk(q)*xk(j)
enddo enddo
enddo enddo
enddo enddo
enddo enddo
do i=1,3 do i = 1,3
do j=1,3 do j = 1,3
ddisgrad(i,j) = 0.
ddisgrad(i,j)=0. ddisgradim(i,j) = 0.
ddisgradim(i,j)=0.
! if(kx.eq.0.and.ky.eq.0.and.kz.eq.0.) goto 4 ! if(kx.eq.0.and.ky.eq.0.and.kz.eq.0.) goto 4
if(.not.(kx.eq.0.and.ky.eq.0.and.kz.eq.0.)) then if(.not.(kx.eq.0.and.ky.eq.0.and.kz.eq.0.)) then
do k = 1,3
do k=1,3 do l = 1,3
do l=1,3 ddisgrad(i,j) = ddisgrad(i,j)+g1(i,j,k,l)*workfft(k,l,kxx,kyy,kzz)
ddisgradim(i,j) = ddisgradim(i,j)+g1(i,j,k,l)*workfftim(k,l,kxx,kyy,kzz)
ddisgrad(i,j)=ddisgrad(i,j)+g1(i,j,k,l)*workfft(k,l,kxx,kyy,kzz)
ddisgradim(i,j)=ddisgradim(i,j)+g1(i,j,k,l)*workfftim(k,l,kxx,kyy,kzz)
enddo enddo
enddo enddo
endif endif
enddo enddo
enddo enddo
workfft(:,:,kxx,kyy,kzz)=ddisgrad(:,:) workfft(:,:,kxx,kyy,kzz) = ddisgrad(:,:)
workfftim(:,:,kxx,kyy,kzz)=ddisgradim(:,:) workfftim(:,:,kxx,kyy,kzz) = ddisgradim(:,:)
enddo enddo
enddo enddo
@ -772,42 +705,39 @@ program mpie_spectral
write(*,*) 'INVERSE FFT TO GET DISPLACEMENT GRADIENT FIELD' write(*,*) 'INVERSE FFT TO GET DISPLACEMENT GRADIENT FIELD'
do m1=1,3 do m1 = 1,3
do n1=1,3 do n1 = 1,3
k1=0 k1 = 0
do k=1,npts3 do k = 1,c
do j=1,npts2 do j = 1,b
do i=1,npts1 do i = 1,a
k1=k1+1 k1 = k1+1
datafft(k1)=workfft(m1,n1,i,j,k) datafft(k1) = workfft(m1,n1,i,j,k)
k1=k1+1
datafft(k1)=workfftim(m1,n1,i,j,k)
k1 = k1+1
datafft(k1) = workfftim(m1,n1,i,j,k)
enddo enddo
enddo enddo
enddo enddo
if(npts3.gt.1) then if(c.gt.1) then
call fourn(datafft,nn,3,-1) call fourn(datafft,nn,3,-1)
else else
call fourn(datafft,nn2,2,-1) call fourn(datafft,nn2,2,-1)
endif endif
datafft = datafft*wgt
datafft=datafft*wgt k1 = 0
do k = 1,c
k1=0 do j = 1,b
do k=1,npts3 do i = 1,a
do j=1,npts2
do i=1,npts1
k1=k1+1
disgrad(m1,n1,i,j,k)=disgrad(m1,n1,i,j,k)+ddisgradmacro(m1,n1)+datafft(k1)
k1=k1+1
k1 = k1+1
disgrad(m1,n1,i,j,k) = disgrad(m1,n1,i,j,k)+ddisgradmacro(m1,n1)+datafft(k1)
k1 = k1+1
enddo enddo
enddo enddo
enddo enddo
@ -816,103 +746,81 @@ program mpie_spectral
enddo enddo
write(*,*) 'UPDATE STRESS FIELD' write(*,*) 'UPDATE STRESS FIELD'
!!! call evpal !!! call evpal
ielem = 0
do k = 1,c
do j = 1,b
do i = 1,a
ielem=0 ielem = ielem+1
do k=1,npts3 defgrad0(:,:) = defgradold(:,:,i,j,k)
do j=1,npts2 defgrad(:,:) = math_i3(:,:)+disgrad(:,:,i,j,k)
do i=1,npts1
ielem=ielem+1
defgrad0(:,:)=defgradold(:,:,i,j,k)
defgrad(:,:)=math_i3(:,:)+disgrad(:,:,i,j,k)
call CPFEM_general(3,defgrad0,defgrad,temperature,0.0_pReal,ielem,1_pInt,stress,dsde) call CPFEM_general(3,defgrad0,defgrad,temperature,0.0_pReal,ielem,1_pInt,stress,dsde)
sg(:,:,i,j,k)=math_Mandel6to33(stress) sg(:,:,i,j,k) = math_Mandel6to33(stress)
enddo enddo
enddo enddo
enddo enddo
ielem=0 ielem = 0
do k=1,npts3 do k = 1,c
do j=1,npts2 do j = 1,b
do i=1,npts1 do i = 1,a
ielem=ielem+1 ielem = ielem+1
defgrad0(:,:)=defgradold(:,:,i,j,k) defgrad0(:,:) = defgradold(:,:,i,j,k)
defgrad(:,:)=math_i3(:,:)+disgrad(:,:,i,j,k) defgrad(:,:) = math_i3(:,:)+disgrad(:,:,i,j,k)
call CPFEM_general(2,defgrad0,defgrad,temperature,0.0_pReal,ielem,1_pInt,stress,dsde) call CPFEM_general(2,defgrad0,defgrad,temperature,0.0_pReal,ielem,1_pInt,stress,dsde)
sg(:,:,i,j,k)=math_Mandel6to33(stress) sg(:,:,i,j,k) = math_Mandel6to33(stress)
enddo enddo
enddo enddo
enddo enddo
!!! call get_smacro !!! call get_smacro
! AVERAGE STRESS ! AVERAGE STRESS
scauav=0. scauav = 0.
do k=1,npts3
do j=1,npts2
do i=1,npts1
scauav(:,:)=scauav(:,:)+sg(:,:,i,j,k)*wgt
do k = 1,c
do j = 1,b
do i = 1,a
scauav(:,:) = scauav(:,:)+sg(:,:,i,j,k)*wgt
enddo enddo
enddo enddo
enddo enddo
! MIXED BC ! MIXED BC
do i=1,3 do i = 1,3
do j=1,3 do j = 1,3
ddisgradmacro(i,j) = 0.
ddisgradmacro(i,j)=0.
if(iudot(i,j).eq.0) then if(iudot(i,j).eq.0) then
do k = 1,3
do k=1,3 do l = 1,3
do l=1,3 ddisgradmacro(i,j) = ddisgradmacro(i,j)+s0(i,j,k,l)*iscau(k,l)*(scauchy(k,l)-scauav(k,l))
ddisgradmacro(i,j)=ddisgradmacro(i,j)+s0(i,j,k,l)*iscau(k,l)*(scauchy(k,l)-scauav(k,l))
enddo enddo
enddo enddo
endif endif
enddo enddo
enddo enddo
ddisgradmacroacum = ddisgradmacroacum+ddisgradmacro
ddisgradmacroacum=ddisgradmacroacum+ddisgradmacro ! write(*,*) 'DDEFGRADMACRO(1,1),(2,2) = ',ddefgradmacro(1,1),ddefgradmacro(2,2)
!! svm = ?
! write(*,*) 'DDEFGRADMACRO(1,1),(2,2)=',ddefgradmacro(1,1),ddefgradmacro(2,2) ! erre = erre/evm
! errs = errs/svm
!! svm= ? write(*,*) 'STRESS FIELD ERROR = ',errs
write(*,*) 'STRAIN FIELD ERROR = ',erre
! erre=erre/evm
! errs=errs/svm
write(*,*) 'STRESS FIELD ERROR =',errs
write(*,*) 'STRAIN FIELD ERROR =',erre
! write(21,101) iter,erre,errs,svm ! write(21,101) iter,erre,errs,svm
!101 format(i3,4(1x,e10.4),10(1x,F7.4)) !101 format(i3,4(1x,e10.4),10(1x,F7.4))
enddo ! WHILE ENDDO enddo ! WHILE ENDDO
disgradmacroactual=disgradmacro+ddisgradmacroacum disgradmacroactual = disgradmacro+ddisgradmacroacum
!! write(*,*) 'defgradmacro(1,1),defgradmacro(2,2),defgradmacro(3,3)' !! write(*,*) 'defgradmacro(1,1),defgradmacro(2,2),defgradmacro(3,3)'
!! write(*,*) defgradmacroactual(1,1),defgradmacroactual(2,2),defgradmacroactual(3,3) !! write(*,*) defgradmacroactual(1,1),defgradmacroactual(2,2),defgradmacroactual(3,3)
@ -921,24 +829,19 @@ program mpie_spectral
!! write(*,*) 'scauav(1,1),scauav(2,2),scauav(3,3)' !! write(*,*) 'scauav(1,1),scauav(2,2),scauav(3,3)'
!! write(*,*) scauav(1,1),scauav(2,2),scauav(3,3) !! write(*,*) scauav(1,1),scauav(2,2),scauav(3,3)
do k=1,npts3 do k = 1,c
do j=1,npts2 do j = 1,b
do i=1,npts1 do i = 1,a
defgradold(:,:,i,j,k) = math_i3(:,:)+disgrad(:,:,i,j,k)
defgradold(:,:,i,j,k)=math_i3(:,:)+disgrad(:,:,i,j,k)
enddo enddo
enddo enddo
enddo enddo
enddo ! IMICRO ENDDO enddo ! IMICRO ENDDO
enddo ! JLOAD ENDDO Ende looping over loadcases
enddo ! JLOAD ENDDO
!------------------------- !-------------------------
!RL !end RL
!------------------------- !-------------------------
end program mpie_spectral end program mpie_spectral
@ -955,11 +858,15 @@ subroutine quit(id)
stop stop
end subroutine end subroutine
SUBROUTINE fourn(data,nn,ndim,isign)
!
!********************************************************************
! fourn subroutine (fourier transform)
! FROM NUMERICAL RECIPES IN F77 (FIXED FORMAT), ! FROM NUMERICAL RECIPES IN F77 (FIXED FORMAT),
! CONVERTED INTO FREE FORMAT (RL @ MPIE, JUNE 2010) ! CONVERTED INTO FREE FORMAT (RL @ MPIE, JUNE 2010)
! !********************************************************************
subroutine fourn(data,nn,ndim,isign)
INTEGER isign,ndim,nn(ndim) INTEGER isign,ndim,nn(ndim)
REAL data(*) REAL data(*)
! INTEGER i1,i2,i2rev,i3,i3rev,ibit,idim,ifp1,ifp2,ip1,ip2,ip3,k1, ! INTEGER i1,i2,i2rev,i3,i3rev,ibit,idim,ifp1,ifp2,ip1,ip2,ip3,k1,
@ -967,96 +874,95 @@ end subroutine
INTEGER i1,i2,i2rev,i3,i3rev,ibit,idim,ifp1,ifp2,ip1,ip2,ip3,k1,k2,n,nprev,nrem,ntot INTEGER i1,i2,i2rev,i3,i3rev,ibit,idim,ifp1,ifp2,ip1,ip2,ip3,k1,k2,n,nprev,nrem,ntot
REAL tempi,tempr REAL tempi,tempr
DOUBLE PRECISION theta,wi,wpi,wpr,wr,wtemp DOUBLE PRECISION theta,wi,wpi,wpr,wr,wtemp
ntot=1 ntot = 1
! do 11 idim=1,ndim ! do 11 idim = 1,ndim
do idim=1,ndim ! 11 do idim = 1,ndim ! 11
! ntot = ntot*nn(idim)
ntot=ntot*nn(idim)
!11 continue !11 continue
enddo ! 11 enddo ! 11
! !
nprev=1 nprev = 1
! do 18 idim=1,ndim ! do 18 idim = 1,ndim
do idim=1,ndim ! 18 do idim = 1,ndim ! 18
n=nn(idim) n = nn(idim)
nrem=ntot/(n*nprev) nrem = ntot/(n*nprev)
ip1=2*nprev ip1 = 2*nprev
ip2=ip1*n ip2 = ip1*n
ip3=ip2*nrem ip3 = ip2*nrem
i2rev=1 i2rev = 1
! do 14 i2=1,ip2,ip1 ! do 14 i2 = 1,ip2,ip1
do i2=1,ip2,ip1 ! 14 do i2 = 1,ip2,ip1 ! 14
if(i2.lt.i2rev)then if(i2.lt.i2rev)then
! do 13 i1=i2,i2+ip1-2,2 ! do 13 i1 = i2,i2+ip1-2,2
! do 12 i3=i1,ip3,ip2 ! do 12 i3 = i1,ip3,ip2
do i1=i2,i2+ip1-2,2 ! 13 do i1 = i2,i2+ip1-2,2 ! 13
do i3=i1,ip3,ip2 ! 12 do i3 = i1,ip3,ip2 ! 12
i3rev=i2rev+i3-i2 i3rev = i2rev+i3-i2
tempr=data(i3) tempr = data(i3)
tempi=data(i3+1) tempi = data(i3+1)
data(i3)=data(i3rev) data(i3) = data(i3rev)
data(i3+1)=data(i3rev+1) data(i3+1) = data(i3rev+1)
data(i3rev)=tempr data(i3rev) = tempr
data(i3rev+1)=tempi data(i3rev+1) = tempi
!12 continue
!13 continue !13 continue
enddo ! 12 enddo ! 12
enddo ! 13 enddo ! 13
endif endif
ibit=ip2/2 ibit = ip2/2
!1 if ((ibit.ge.ip1).and.(i2rev.gt.ibit)) then !1 if ((ibit.ge.ip1).and.(i2rev.gt.ibit)) then
do while ((ibit.ge.ip1).and.(i2rev.gt.ibit)) ! if 1 do while ((ibit.ge.ip1).and.(i2rev.gt.ibit)) ! if 1
i2rev=i2rev-ibit i2rev = i2rev-ibit
ibit=ibit/2 ibit = ibit/2
! goto 1 ! goto 1
! endif ! endif
enddo ! do while (if 1) enddo ! do while (if 1)
i2rev=i2rev+ibit i2rev = i2rev+ibit
!14 continue !14 continue
enddo ! 14 enddo ! 14
ifp1=ip1 ifp1 = ip1
!2 if(ifp1.lt.ip2)then !2 if(ifp1.lt.ip2)then
do while (ifp1.lt.ip2) ! if 2 do while (ifp1.lt.ip2) ! if 2
ifp2=2*ifp1 ifp2 = 2*ifp1
theta=isign*6.28318530717959d0/(ifp2/ip1) theta = isign*6.28318530717959d0/(ifp2/ip1)
wpr=-2.d0*sin(0.5d0*theta)**2 wpr = -2.d0*sin(0.5d0*theta)**2
wpi=sin(theta) wpi = sin(theta)
wr=1.d0 wr = 1.d0
wi=0.d0 wi = 0.d0
! do 17 i3=1,ifp1,ip1 ! do 17 i3 = 1,ifp1,ip1
! do 16 i1=i3,i3+ip1-2,2 ! do 16 i1 = i3,i3+ip1-2,2
! do 15 i2=i1,ip3,ifp2 ! do 15 i2 = i1,ip3,ifp2
do i3=1,ifp1,ip1 ! 17 do i3 = 1,ifp1,ip1 ! 17
do i1=i3,i3+ip1-2,2 ! 16 do i1 = i3,i3+ip1-2,2 ! 16
do i2=i1,ip3,ifp2 ! 15 do i2 = i1,ip3,ifp2 ! 15
k1=i2 k1 = i2
k2=k1+ifp1 k2 = k1+ifp1
tempr=sngl(wr)*data(k2)-sngl(wi)*data(k2+1) tempr = sngl(wr)*data(k2)-sngl(wi)*data(k2+1)
tempi=sngl(wr)*data(k2+1)+sngl(wi)*data(k2) tempi = sngl(wr)*data(k2+1)+sngl(wi)*data(k2)
data(k2)=data(k1)-tempr data(k2) = data(k1)-tempr
data(k2+1)=data(k1+1)-tempi data(k2+1) = data(k1+1)-tempi
data(k1)=data(k1)+tempr data(k1) = data(k1)+tempr
data(k1+1)=data(k1+1)+tempi data(k1+1) = data(k1+1)+tempi
!15 continue !15 continue
!16 continue !16 continue
enddo ! 15 enddo ! 15
enddo ! 16 enddo ! 16
wtemp=wr wtemp = wr
wr=wr*wpr-wi*wpi+wr wr = wr*wpr-wi*wpi+wr
wi=wi*wpr+wtemp*wpi+wi wi = wi*wpr+wtemp*wpi+wi
!17 continue !17 continue
enddo ! 17 enddo ! 17
ifp1=ifp2 ifp1 = ifp2
! goto 2 ! goto 2
! endif ! endif
enddo ! do while (if 2) enddo ! do while (if 2)
nprev=n*nprev nprev = n*nprev
!18 continue !18 continue
enddo ! 18 enddo ! 18
return return