From d6ba9d54b658e7771618ae3e30aa8c57be92a2ef Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 2 Jul 2010 14:10:36 +0000 Subject: [PATCH] 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 --- code/mpie_spectral.f90 | 884 ++++++++++++++++++----------------------- 1 file changed, 395 insertions(+), 489 deletions(-) diff --git a/code/mpie_spectral.f90 b/code/mpie_spectral.f90 index 75e899fce..be3f04b2b 100644 --- a/code/mpie_spectral.f90 +++ b/code/mpie_spectral.f90 @@ -44,7 +44,7 @@ subroutine mpie_interface_init() write(6,*) return -end subroutine +endsubroutine !******************************************************************** ! extract working directory from loadcase file @@ -70,7 +70,7 @@ function getSolverWorkingDirectoryName() return -end function +endfunction !******************************************************************** ! basename of meshfile from command line arguments @@ -105,7 +105,7 @@ function getSolverJobName() getSolverJobName = makeRelativePath(getSolverWorkingDirectoryName(),& getSolverJobName) return -end function +endfunction !******************************************************************** @@ -121,7 +121,7 @@ function getLoadcaseName() character(len=1024) getLoadcaseName, outName, cwd character(len=*), parameter :: pathSep = achar(47)//achar(92) ! /, \ integer(pInt) posExt,posSep - posExt=0 !not sure if its needed + posExt = 0 !not sure if its needed call getarg(2,getLoadcaseName) posExt = scan(getLoadcaseName,'.',back=.true.) @@ -138,7 +138,7 @@ function getLoadcaseName() getLoadcaseName = makeRelativePath(getSolverWorkingDirectoryName(),& getLoadcaseName) return -end function +endfunction !******************************************************************** @@ -158,10 +158,10 @@ function rectifyPath(path) !remove ./ from path l = len_trim(path) rectifyPath = path - do i=l,2,-1 - if ( rectifyPath(i-1:i)=='./' .and. rectifyPath(i-2:i-2) /= '.' ) & + do i = l,2,-1 + if ( rectifyPath(i-1:i) == './' .and. rectifyPath(i-2:i-2) /= '.' ) & rectifyPath(i-1:l) = rectifyPath(i+1:l)//' ' - end do + enddo !remove ../ and corresponding directory from rectifyPath l = len_trim(rectifyPath) @@ -173,11 +173,11 @@ function rectifyPath(path) j = scan(rectifyPath(:i-2),'/',back=.true.) rectifyPath(j+1:l) = rectifyPath(i+3:l)//repeat(' ',2+i-j) i = j+index(rectifyPath(j+1:l),'../') - end do - if(len_trim(rectifyPath)==0) rectifyPath='/' + enddo + if(len_trim(rectifyPath) == 0) rectifyPath = '/' return - end function rectifyPath + endfunction rectifyPath !******************************************************************** @@ -205,7 +205,7 @@ function makeRelativePath(a,b) enddo makeRelativePath = repeat('../',remainingSlashes)//b(posLastCommonSlash+1:len_trim(b)) return -end function makeRelativePath +endfunction makeRelativePath 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) a, b, c, e, homog real(pReal) x, y, z + !------------------------- !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 :: workfft,workfftim,sg,disgrad,defgradold @@ -291,13 +284,15 @@ program mpie_spectral real(pReal), dimension(3,3) :: defgrad0,defgrad real(pReal), dimension(3,3) :: udot,scauchy,scauav,aux33 - real(pReal), dimension(3) :: delt,nn,xk,xk2 - real(pReal), dimension(2) :: nn2 + integer(pInt), dimension(3) :: nn + integer(pInt), dimension(2) :: nn2 + + real(pReal), dimension(3) :: delt,xk,xk2 real(pReal), dimension(6) :: aux6 real(pReal), dimension(3,3,3,3) :: c0,s0,g1 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 @@ -320,11 +315,11 @@ program mpie_spectral print*,'Loadcase: ',trim(path) 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) do - read(unit,'(a1024)',END=101) line + read(unit,'(a1024)',END = 101) line if (IO_isBlank(line)) cycle ! skip empty lines posInput = IO_stringPos(line,maxNchunksInput) do i = 1,maxNchunksInput,1 @@ -340,7 +335,7 @@ program mpie_spectral end select 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 - 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 @@ -355,9 +350,9 @@ program mpie_spectral rewind(unit) j = 0_pInt do - read(unit,'(a1024)',END=200) line + read(unit,'(a1024)',END = 200) line if (IO_isBlank(line)) cycle ! skip empty lines - j=j+1 + j = j+1 posInput = IO_stringPos(line,maxNchunksInput) do i = 1,maxNchunksInput,2 select case (IO_lc(IO_stringValue(line,posInput,i))) @@ -409,16 +404,16 @@ program mpie_spectral y = 1_pReal z = 1_pReal - gotResolution = .false. - gotDimension = .false. - gotHomogenization = .false. + gotResolution = .false. + gotDimension = .false. + gotHomogenization = .false. path = getSolverJobName() 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) do - read(unit,'(a1024)',END=100) line + read(unit,'(a1024)',END = 100) line if (IO_isBlank(line)) cycle ! skip empty lines posMesh = IO_stringPos(line,maxNchunksMesh) @@ -455,274 +450,219 @@ program mpie_spectral enddo 100 close(unit) - print '(a,/,i3,i3,i3)','resolution a b c',a,b,c - print '(a,/,f6.2,f6.2,f6.2)','dimension x y z',x, y, z - print *,'homogenization',homog - print *, '' + print '(a,/,i3,i3,i3)','resolution a b c',a,b,c + print '(a,/,f6.2,f6.2,f6.2)','dimension x y z',x, y, z + print *,'homogenization',homog + print *, '' temperature = 300.0_pReal - !------------------------- -!RL + !------------------------- -!begin change m.diehl - npts1=a - npts2=b - npts3=c -!end change m.diehl +!begin RL +!------------------------- + allocate (datafft(2*a*b*c)) + allocate (workfft(3,3,a,b,c)) + 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)) - allocate (workfft(3,3,npts1,npts2,npts3)) - 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 + itmax = 100 - error=0.00001 - itmax=100 + delt(1) = 1. + delt(2) = 1. + delt(3) = 1. - delt(1)=1. - delt(2)=1. - delt(3)=1. + nn(1) = a + nn(2) = b + nn(3) = c + + nn2(1) = a + nn2(2) = b - nn(1)=npts1 - nn(2)=npts2 - nn(3)=npts3 - - nn2(1)=npts1 - nn2(2)=npts2 - - prodnn=nn(1)*nn(2)*nn(3) - wgt=1./prodnn + prodnn = nn(1)*nn(2)*nn(3) + wgt = 1./prodnn ! C_0 and S_0 CALCULATION - !! !! 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 - 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) - - c066=c066+dsde - c0=c0+math_Mandel66to3333(dsde) - - enddo - enddo + do ielem = 1, prodnn !call each element with identity (math_i3) to initialize with high stress + call CPFEM_general(3,math_i3,math_i3,temperature,0.0_pReal,ielem,1_pInt,stress,dsde) + c066 = c066+dsde + c0 = c0+math_Mandel66to3333(dsde) enddo - c066=c066*wgt - c0=c0*wgt + c066 = c066*wgt + c0 = c0*wgt + call math_invert(6,c066,s066,idum,errmatinv) - if(errmatinv) then - write(*,*) 'ERROR IN C0 INVERSION' - stop + write(*,*) 'ERROR IN C0 INVERSION' + stop endif - s0=math_Mandel66to3333(s066) + s0 = math_Mandel66to3333(s066) ! 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(:,:) - - enddo - enddo - enddo - - do jload = 1,N - - udot(:,:)=bc_velocityGrad(:,:,jload) - scauchy(:,:)=bc_stress(:,:,jload) - - iudot=0 - 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 + do k = 1, c + do j = 1, b + do i = 1, a + defgradold(:,:,i,j,k) = math_i3(:,:) + enddo + enddo enddo - enddo - - nsteps=bc_steps(jload) - tdot=bc_timeIncrement(jload)/bc_steps(jload) -! evm=dvm*tdot ? - - do imicro=1,nsteps - - write(*,*) '***************************************************' - write(*,*) 'STEP = ',imicro -! write(21,*) 'STEP = ',imicro - -! INITIALIZATION BEFORE NEW TIME STEP - - disgradmacro=disgradmacro+udot*tdot - ddisgradmacro=0. - - ielem=0 - - do k=1,npts3 - do j=1,npts2 - do i=1,npts1 - - ielem=ielem+1 - - disgrad(:,:,i,j,k)=disgradmacro(:,:) - - 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) - -! sg(:,:,i,j,k)=math_Mandel6to33(stress) ? - - enddo - enddo - enddo - - ielem=0 - - do k=1,npts3 - do j=1,npts2 - do i=1,npts1 - - ielem=ielem+1 - - disgrad(:,:,i,j,k)=disgradmacro(:,:) - - 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) - - sg(:,:,i,j,k)=math_Mandel6to33(stress) - - enddo - enddo + do jload = 1,N !Loop over loadcases defined in the loadcase file + udot(:,:) = bc_velocityGrad(:,:,jload) + scauchy(:,:) = bc_stress(:,:,jload) + iudot = 0 + iscau = 0 + + do i = 1, 3 !convert information about rb's from bc_mask in corresponding arrays + 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 - ddisgradmacroacum=0. + 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 + write(*,*) '***************************************************' + write(*,*) 'STEP = ',imicro + +! INITIALIZATION BEFORE NEW TIME STEP + disgradmacro = disgradmacro+udot*tdot !update macroscopic displacementgradient + ddisgradmacro = 0. + ielem=0 - iter=0 - erre=2.*error - errs=2.*error + do k = 1, c !loop over FPs + do j = 1, b + do i = 1, a + ielem = ielem+1 + disgrad(:,:,i,j,k)=disgradmacro(:,:) - do while(iter.lt.itmax.and.(errs.gt.error.or.erre.gt.error)) + defgrad0(:,:)=defgradold(:,:,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) +! sg(:,:,i,j,k)=math_Mandel6to33(stress) ? + enddo + enddo + enddo + + ielem = 0 + + do k = 1, c !loop over FPs + do j = 1, b + do i = 1, a + ielem = ielem+1 + disgrad(:,:,i,j,k) = disgradmacro(:,:) + + 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) + sg(:,:,i,j,k) = math_Mandel6to33(stress) + enddo + enddo + enddo + + ddisgradmacroacum = 0. - iter=iter+1 - - write(*,*)'ITER = ',iter - write(*,*) 'DIRECT FFT OF STRESS FIELD' - - do ii=1,3 - do jj=1,3 - - k1=0 - - do k=1,npts3 + iter = 0 + erre = 2.*error + errs = 2.*error + do while(iter <= itmax.and.(errs.gt.error.or.erre.gt.error)) + iter = iter+1 + write(*,*)'ITER = ',iter + write(*,*) 'DIRECT FFT OF STRESS FIELD' + do ii = 1,3 + do jj = 1,3 + k1 = 0 + do k = 1,c ! write(*,'(1H+,a,i2,2(a,i4))') ! # 'STRESS - COMPONENT',ii,jj,' - Z = ',k,' OUT OF ',npts + do j = 1,b + do i = 1,a + k1 = k1+1 + datafft(k1) = sg(ii,jj,i,j,k) + k1 = k1+1 + datafft(k1) = 0. + enddo + enddo + enddo + if(c.gt.1) then + call fourn(datafft,nn,3,1) + else + call fourn(datafft,nn2,2,1) + endif + k1 = 0 + do k = 1, c + do j = 1, b + do i = 1, a + k1 = k1+1 + workfft(ii,jj,i,j,k) = datafft(k1) + k1 = k1+1 + workfftim(ii,jj,i,j,k) = datafft(k1) + enddo + enddo + enddo + enddo + enddo + 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 + + if(kyy.le.b/2) ky = kyy-1 + if(kyy.gt.b/2) ky = kyy-b-1 + + if(kzz.le.c/2) kz = kzz-1 + if(kzz.gt.c/2) kz = kzz-c-1 - do j=1,npts2 - do i=1,npts1 + xk(1) = kx/(delt(1)*nn(1)) + xk(2) = ky/(delt(2)*nn(2)) - k1=k1+1 - datafft(k1)=sg(ii,jj,i,j,k) - k1=k1+1 - datafft(k1)=0. - - enddo - enddo - enddo - - if(npts3.gt.1) then - call fourn(datafft,nn,3,1) - else - call fourn(datafft,nn2,2,1) - endif - - k1=0 - - do k=1,npts3 - do j=1,npts2 - do i=1,npts1 - - k1=k1+1 - workfft(ii,jj,i,j,k)=datafft(k1) - - k1=k1+1 - workfftim(ii,jj,i,j,k)=datafft(k1) - - enddo - enddo - enddo - - enddo - enddo - - write(*,*) 'CALCULATING G^pq,ij : SG^ij ...' - - do kzz=1,npts3 - do kyy=1,npts2 - do kxx=1,npts1 - - if(kxx.le.npts1/2) kx=kxx-1 - if(kxx.gt.npts1/2) kx=kxx-npts1-1 - - if(kyy.le.npts2/2) ky=kyy-1 - if(kyy.gt.npts2/2) ky=kyy-npts2-1 - - if(kzz.le.npts3/2) kz=kzz-1 - if(kzz.gt.npts3/2) kz=kzz-npts3-1 - - 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 - xk(3)=0. - endif - - xknorm=sqrt(xk(1)**2+xk(2)**2+xk(3)**2) - - if (xknorm.ne.0.) then - do i=1,3 - xk2(i)=xk(i)/(xknorm*xknorm*2.*pi) - xk(i)=xk(i)/xknorm - enddo - endif - - do i=1,3 - do k=1,3 - aux33(i,k)=0. - do j=1,3 - do l=1,3 - aux33(i,k)=aux33(i,k)+c0(i,j,k,l)*xk(j)*xk(l) - enddo - enddo - enddo - enddo - - aux33=math_inv3x3(aux33) + if(c.gt.1) then + xk(3) = kz/(delt(3)*nn(3)) + else + xk(3) = 0. + endif + + xknorm = sqrt(xk(1)**2+xk(2)**2+xk(3)**2) + if (xknorm.ne.0.) then + do i = 1,3 + xk2(i) = xk(i)/(xknorm*xknorm*2.*pi) + xk(i) = xk(i)/xknorm + enddo + endif + + do i = 1,3 + do k = 1,3 + aux33(i,k) = 0. + do j = 1,3 + do l = 1,3 + aux33(i,k) = aux33(i,k)+c0(i,j,k,l)*xk(j)*xk(l) + enddo + enddo + enddo + enddo + aux33 = math_inv3x3(aux33) + ! call minv(aux33,3,det,minv1,minv2) ! if(det.eq.0) then ! write(*,*) kx,ky,kz,' --> SINGULAR SYSTEM' @@ -730,189 +670,157 @@ program mpie_spectral ! pause ! endif - do p=1,3 - do q=1,3 - do i=1,3 - do j=1,3 - g1(p,q,i,j)=-aux33(p,i)*xk(q)*xk(j) - enddo - enddo - enddo - enddo - - do i=1,3 - do j=1,3 - - ddisgrad(i,j)=0. - ddisgradim(i,j)=0. + do p = 1,3 + do q = 1,3 + do i = 1,3 + do j = 1,3 + g1(p,q,i,j) = -aux33(p,i)*xk(q)*xk(j) + enddo + enddo + enddo + enddo + do i = 1,3 + do j = 1,3 + ddisgrad(i,j) = 0. + ddisgradim(i,j) = 0. ! 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 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) + enddo + enddo + endif + enddo + enddo - do k=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) - - enddo + workfft(:,:,kxx,kyy,kzz) = ddisgrad(:,:) + workfftim(:,:,kxx,kyy,kzz) = ddisgradim(:,:) + + enddo + enddo enddo - endif + write(*,*) 'INVERSE FFT TO GET DISPLACEMENT GRADIENT FIELD' - enddo - enddo + do m1 = 1,3 + do n1 = 1,3 - workfft(:,:,kxx,kyy,kzz)=ddisgrad(:,:) - workfftim(:,:,kxx,kyy,kzz)=ddisgradim(:,:) + k1 = 0 - enddo - enddo - enddo + do k = 1,c + do j = 1,b + do i = 1,a + + k1 = k1+1 + datafft(k1) = workfft(m1,n1,i,j,k) + + k1 = k1+1 + datafft(k1) = workfftim(m1,n1,i,j,k) + enddo + enddo + enddo - write(*,*) 'INVERSE FFT TO GET DISPLACEMENT GRADIENT FIELD' + if(c.gt.1) then + call fourn(datafft,nn,3,-1) + else + call fourn(datafft,nn2,2,-1) + endif + datafft = datafft*wgt - do m1=1,3 - do n1=1,3 + k1 = 0 + do k = 1,c + do j = 1,b + do i = 1,a - k1=0 + 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 - do k=1,npts3 - do j=1,npts2 - do i=1,npts1 - - k1=k1+1 - datafft(k1)=workfft(m1,n1,i,j,k) - - k1=k1+1 - datafft(k1)=workfftim(m1,n1,i,j,k) - - enddo - enddo - enddo - - if(npts3.gt.1) then - call fourn(datafft,nn,3,-1) - else - call fourn(datafft,nn2,2,-1) - endif - - datafft=datafft*wgt - - k1=0 - do k=1,npts3 - 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 - - enddo - enddo - enddo - - enddo - enddo - - write(*,*) 'UPDATE STRESS FIELD' + enddo + enddo + write(*,*) 'UPDATE STRESS FIELD' !!! 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 - do j=1,npts2 - do i=1,npts1 + defgrad0(:,:) = defgradold(:,:,i,j,k) + defgrad(:,:) = math_i3(:,:)+disgrad(:,:,i,j,k) - ielem=ielem+1 + call CPFEM_general(3,defgrad0,defgrad,temperature,0.0_pReal,ielem,1_pInt,stress,dsde) - defgrad0(:,:)=defgradold(:,:,i,j,k) - defgrad(:,:)=math_i3(:,:)+disgrad(:,:,i,j,k) + sg(:,:,i,j,k) = math_Mandel6to33(stress) + enddo + enddo + enddo - call CPFEM_general(3,defgrad0,defgrad,temperature,0.0_pReal,ielem,1_pInt,stress,dsde) + ielem = 0 - sg(:,:,i,j,k)=math_Mandel6to33(stress) + do k = 1,c + do j = 1,b + do i = 1,a - enddo - enddo - enddo + ielem = ielem+1 - ielem=0 + defgrad0(:,:) = defgradold(:,:,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) - do k=1,npts3 - do j=1,npts2 - do i=1,npts1 - - ielem=ielem+1 - - defgrad0(:,:)=defgradold(:,:,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) - - sg(:,:,i,j,k)=math_Mandel6to33(stress) - - enddo - enddo - enddo + sg(:,:,i,j,k) = math_Mandel6to33(stress) + enddo + enddo + enddo !!! call get_smacro - ! 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 - - enddo - enddo - enddo - + do k = 1,c + do j = 1,b + do i = 1,a + scauav(:,:) = scauav(:,:)+sg(:,:,i,j,k)*wgt + enddo + enddo + enddo ! MIXED BC - do i=1,3 - do j=1,3 - - ddisgradmacro(i,j)=0. - - if(iudot(i,j).eq.0) then - - do k=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)) - - enddo - enddo - - endif - - enddo - enddo - - ddisgradmacroacum=ddisgradmacroacum+ddisgradmacro - -! write(*,*) 'DDEFGRADMACRO(1,1),(2,2)=',ddefgradmacro(1,1),ddefgradmacro(2,2) - -!! svm= ? - -! erre=erre/evm -! errs=errs/svm - - write(*,*) 'STRESS FIELD ERROR =',errs - write(*,*) 'STRAIN FIELD ERROR =',erre - + do i = 1,3 + do j = 1,3 + ddisgradmacro(i,j) = 0. + if(iudot(i,j).eq.0) then + do k = 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)) + enddo + enddo + endif + enddo + enddo + ddisgradmacroacum = ddisgradmacroacum+ddisgradmacro +! write(*,*) 'DDEFGRADMACRO(1,1),(2,2) = ',ddefgradmacro(1,1),ddefgradmacro(2,2) +!! svm = ? +! erre = erre/evm +! errs = errs/svm + write(*,*) 'STRESS FIELD ERROR = ',errs + write(*,*) 'STRAIN FIELD ERROR = ',erre ! write(21,101) iter,erre,errs,svm !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(*,*) 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) - do k=1,npts3 - do j=1,npts2 - do i=1,npts1 - - defgradold(:,:,i,j,k)=math_i3(:,:)+disgrad(:,:,i,j,k) - - enddo - enddo - enddo - - enddo ! IMICRO ENDDO - - enddo ! JLOAD ENDDO + do k = 1,c + do j = 1,b + do i = 1,a + defgradold(:,:,i,j,k) = math_i3(:,:)+disgrad(:,:,i,j,k) + enddo + enddo + enddo + enddo ! IMICRO ENDDO + enddo ! JLOAD ENDDO Ende looping over loadcases !------------------------- -!RL +!end RL !------------------------- - end program mpie_spectral @@ -955,11 +858,15 @@ subroutine quit(id) stop end subroutine - SUBROUTINE fourn(data,nn,ndim,isign) -! + + +!******************************************************************** +! fourn subroutine (fourier transform) ! FROM NUMERICAL RECIPES IN F77 (FIXED FORMAT), ! CONVERTED INTO FREE FORMAT (RL @ MPIE, JUNE 2010) -! +!******************************************************************** +subroutine fourn(data,nn,ndim,isign) + INTEGER isign,ndim,nn(ndim) REAL data(*) ! 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 REAL tempi,tempr DOUBLE PRECISION theta,wi,wpi,wpr,wr,wtemp - ntot=1 -! do 11 idim=1,ndim - do idim=1,ndim ! 11 -! - ntot=ntot*nn(idim) + ntot = 1 +! do 11 idim = 1,ndim + do idim = 1,ndim ! 11 + ntot = ntot*nn(idim) !11 continue enddo ! 11 ! - nprev=1 -! do 18 idim=1,ndim - do idim=1,ndim ! 18 - n=nn(idim) - nrem=ntot/(n*nprev) - ip1=2*nprev - ip2=ip1*n - ip3=ip2*nrem - i2rev=1 -! do 14 i2=1,ip2,ip1 - do i2=1,ip2,ip1 ! 14 + nprev = 1 +! do 18 idim = 1,ndim + do idim = 1,ndim ! 18 + n = nn(idim) + nrem = ntot/(n*nprev) + ip1 = 2*nprev + ip2 = ip1*n + ip3 = ip2*nrem + i2rev = 1 +! do 14 i2 = 1,ip2,ip1 + do i2 = 1,ip2,ip1 ! 14 if(i2.lt.i2rev)then -! do 13 i1=i2,i2+ip1-2,2 -! do 12 i3=i1,ip3,ip2 - do i1=i2,i2+ip1-2,2 ! 13 - do i3=i1,ip3,ip2 ! 12 - i3rev=i2rev+i3-i2 - tempr=data(i3) - tempi=data(i3+1) - data(i3)=data(i3rev) - data(i3+1)=data(i3rev+1) - data(i3rev)=tempr - data(i3rev+1)=tempi -!12 continue +! do 13 i1 = i2,i2+ip1-2,2 +! do 12 i3 = i1,ip3,ip2 + do i1 = i2,i2+ip1-2,2 ! 13 + do i3 = i1,ip3,ip2 ! 12 + i3rev = i2rev+i3-i2 + tempr = data(i3) + tempi = data(i3+1) + data(i3) = data(i3rev) + data(i3+1) = data(i3rev+1) + data(i3rev) = tempr + data(i3rev+1) = tempi + !13 continue enddo ! 12 enddo ! 13 endif - ibit=ip2/2 + ibit = ip2/2 !1 if ((ibit.ge.ip1).and.(i2rev.gt.ibit)) then do while ((ibit.ge.ip1).and.(i2rev.gt.ibit)) ! if 1 - i2rev=i2rev-ibit - ibit=ibit/2 + i2rev = i2rev-ibit + ibit = ibit/2 ! goto 1 ! endif enddo ! do while (if 1) - i2rev=i2rev+ibit + i2rev = i2rev+ibit !14 continue enddo ! 14 - ifp1=ip1 + ifp1 = ip1 !2 if(ifp1.lt.ip2)then do while (ifp1.lt.ip2) ! if 2 - ifp2=2*ifp1 - theta=isign*6.28318530717959d0/(ifp2/ip1) - wpr=-2.d0*sin(0.5d0*theta)**2 - wpi=sin(theta) - wr=1.d0 - wi=0.d0 -! do 17 i3=1,ifp1,ip1 -! do 16 i1=i3,i3+ip1-2,2 -! do 15 i2=i1,ip3,ifp2 - do i3=1,ifp1,ip1 ! 17 - do i1=i3,i3+ip1-2,2 ! 16 - do i2=i1,ip3,ifp2 ! 15 - k1=i2 - k2=k1+ifp1 - tempr=sngl(wr)*data(k2)-sngl(wi)*data(k2+1) - tempi=sngl(wr)*data(k2+1)+sngl(wi)*data(k2) - data(k2)=data(k1)-tempr + ifp2 = 2*ifp1 + theta = isign*6.28318530717959d0/(ifp2/ip1) + wpr = -2.d0*sin(0.5d0*theta)**2 + wpi = sin(theta) + wr = 1.d0 + wi = 0.d0 +! do 17 i3 = 1,ifp1,ip1 +! do 16 i1 = i3,i3+ip1-2,2 +! do 15 i2 = i1,ip3,ifp2 + do i3 = 1,ifp1,ip1 ! 17 + do i1 = i3,i3+ip1-2,2 ! 16 + do i2 = i1,ip3,ifp2 ! 15 + k1 = i2 + k2 = k1+ifp1 + tempr = sngl(wr)*data(k2)-sngl(wi)*data(k2+1) + tempi = sngl(wr)*data(k2+1)+sngl(wi)*data(k2) + data(k2) = data(k1)-tempr - data(k2+1)=data(k1+1)-tempi - data(k1)=data(k1)+tempr - data(k1+1)=data(k1+1)+tempi + data(k2+1) = data(k1+1)-tempi + data(k1) = data(k1)+tempr + data(k1+1) = data(k1+1)+tempi !15 continue !16 continue enddo ! 15 enddo ! 16 - wtemp=wr - wr=wr*wpr-wi*wpi+wr - wi=wi*wpr+wtemp*wpi+wi + wtemp = wr + wr = wr*wpr-wi*wpi+wr + wi = wi*wpr+wtemp*wpi+wi !17 continue enddo ! 17 - ifp1=ifp2 + ifp1 = ifp2 ! goto 2 ! endif enddo ! do while (if 2) - nprev=n*nprev + nprev = n*nprev !18 continue enddo ! 18 return