diff --git a/code/mpie_spectral.f90 b/code/mpie_spectral.f90 index 6b5a09df4..75e899fce 100644 --- a/code/mpie_spectral.f90 +++ b/code/mpie_spectral.f90 @@ -92,22 +92,24 @@ function getSolverJobName() posExt = scan(outName,'.',back=.true.) posSep = scan(outName,pathSep,back=.true.) - if (posExt <= posSep) posExt = len_trim(outName) ! no extension present - - getSolverJobName = outName(posSep+1:posExt-1) ! path to mesh file (excl. extension) - - if (scan(getSolverJobName,pathSep) /= 1) then ! relative path given as command line argument + if (posExt <= posSep) posExt = len_trim(outName)+1 ! no extension present + getSolverJobName = outName(1:posExt-1) ! path to mesh file (excl. extension) + + if (scan(getSolverJobName,pathSep) /= 1) then ! relative path given as command line argument call getcwd(cwd) - getSolverJobName = makeRelativePath(getSolverWorkingDirectoryName(),& - rectifyPath(trim(cwd)//'/'//getSolverJobName)) + getSolverJobName = rectifyPath(trim(cwd)//'/'//getSolverJobName) + else + getSolverJobName = rectifyPath(getSolverJobName) endif - + + getSolverJobName = makeRelativePath(getSolverWorkingDirectoryName(),& + getSolverJobName) return end function !******************************************************************** -! realive path of loadcase from command line arguments +! relative path of loadcase from command line arguments ! !******************************************************************** function getLoadcaseName() @@ -119,17 +121,22 @@ 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 call getarg(2,getLoadcaseName) - + posExt = scan(getLoadcaseName,'.',back=.true.) + posSep = scan(getLoadcaseName,pathSep,back=.true.) + + if (posExt <= posSep) getLoadcaseName = trim(getLoadcaseName)//('.load') ! no extension present if (scan(getLoadcaseName,pathSep) /= 1) then ! relative path given as command line argument call getcwd(cwd) getLoadcaseName = rectifyPath(trim(cwd)//'/'//getLoadcaseName) + else + getLoadcaseName = rectifyPath(getLoadcaseName) endif getLoadcaseName = makeRelativePath(getSolverWorkingDirectoryName(),& getLoadcaseName) - return end function @@ -167,6 +174,7 @@ function rectifyPath(path) 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='/' return end function rectifyPath @@ -188,7 +196,6 @@ function makeRelativePath(a,b) posLastCommonSlash = 0 remainingSlashes = 0 - do i = 1,min(1024,len_trim(a),len_trim(b)) if (a(i:i) /= b(i:i)) exit if (a(i:i) == '/') posLastCommonSlash = i @@ -197,7 +204,6 @@ function makeRelativePath(a,b) if (a(i:i) == '/') remainingSlashes = remainingSlashes + 1 enddo makeRelativePath = repeat('../',remainingSlashes)//b(posLastCommonSlash+1:len_trim(b)) - return end function makeRelativePath @@ -233,7 +239,8 @@ program mpie_spectral use mpie_interface use prec, only: pInt, pReal use IO - use math, only: math_I3,math_transpose3x3 +! use math, only: math_I3,math_transpose3x3,math_Mandel66to3333 + use math use CPFEM, only: CPFEM_general implicit none @@ -262,6 +269,43 @@ 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 + + integer(pInt), dimension (3,3) :: iudot,iscau + + real(pReal), dimension(3,3) :: disgradmacro,disgradmacroactual + real(pReal), dimension(3,3) :: ddisgradmacro,ddisgradmacroacum,ddisgrad,ddisgradim + 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 + 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 + + real(pReal) prodnn,wgt,error,tdot,erre,errs,evm,svm,det,xknorm + + logical errmatinv + +!------------------------- +!end RL +!------------------------- if (IargC() < 2) call IO_error(102) @@ -272,6 +316,9 @@ program mpie_spectral N_s = 0_pInt N_t = 0_pInt N_n = 0_pInt + + print*,'Loadcase: ',trim(path) + print*,'Workingdir: ',trim(getSolverWorkingDirectoryName()) if (.not. IO_open_file(unit,path)) call IO_error(45,ext_msg=path) @@ -366,7 +413,8 @@ program mpie_spectral gotDimension = .false. gotHomogenization = .false. path = getSolverJobName() - if (.not. IO_open_file(unit,path)) call IO_error(101,ext_msg=path) + print*,'JobName: ',trim(path) + if (.not. IO_open_file(unit,trim(path)//InputFileExtension)) call IO_error(101,ext_msg=path) rewind(unit) do @@ -405,9 +453,491 @@ program mpie_spectral end select if (gotDimension .and. gotHomogenization .and. gotResolution) exit 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 *, '' temperature = 300.0_pReal - call CPFEM_general(2,math_i3,math_i3,temperature,0.0_pReal,1_pInt,1_pInt,stress,dsde) + !------------------------- +!RL +!------------------------- +!begin change m.diehl + npts1=a + npts2=b + npts3=c +!end change m.diehl + + 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 + + delt(1)=1. + delt(2)=1. + delt(3)=1. + + 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 + +! C_0 and S_0 CALCULATION + +!! +!! PHILIP: FE_exec_elem? +!! + + 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 + enddo + + c066=c066*wgt + c0=c0*wgt + + call math_invert(6,c066,s066,idum,errmatinv) + + if(errmatinv) then + write(*,*) 'ERROR IN C0 INVERSION' + stop + endif + + s0=math_Mandel66to3333(s066) + +! INITIALIZATION BEFORE STARTING WITH LOADINGS + + 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 + 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 + enddo + + ddisgradmacroacum=0. + + iter=0 + erre=2.*error + errs=2.*error + + do while(iter.lt.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,npts3 + +! write(*,'(1H+,a,i2,2(a,i4))') +! # 'STRESS - COMPONENT',ii,jj,' - Z = ',k,' OUT OF ',npts + + do j=1,npts2 + do i=1,npts1 + + 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) + +! call minv(aux33,3,det,minv1,minv2) +! if(det.eq.0) then +! write(*,*) kx,ky,kz,' --> SINGULAR SYSTEM' +! stop +! 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. + +! 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 + + 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 + + workfft(:,:,kxx,kyy,kzz)=ddisgrad(:,:) + workfftim(:,:,kxx,kyy,kzz)=ddisgradim(:,:) + + enddo + enddo + enddo + + write(*,*) 'INVERSE FFT TO GET DISPLACEMENT GRADIENT FIELD' + + do m1=1,3 + do n1=1,3 + + k1=0 + + 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' + +!!! call evpal + + ielem=0 + + 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(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 + + 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 + +!!! call get_smacro + +! AVERAGE STRESS + + scauav=0. + + do k=1,npts3 + do j=1,npts2 + do i=1,npts1 + + 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 + +! write(21,101) iter,erre,errs,svm +!101 format(i3,4(1x,e10.4),10(1x,F7.4)) + + enddo ! WHILE ENDDO + + disgradmacroactual=disgradmacro+ddisgradmacroacum + +!! write(*,*) 'defgradmacro(1,1),defgradmacro(2,2),defgradmacro(3,3)' +!! write(*,*) defgradmacroactual(1,1),defgradmacroactual(2,2),defgradmacroactual(3,3) +!! write(*,*) 'defgradmacro(1,1)/defgradmacro(3,3)' +!! write(*,*) defgradmacroactual(1,1)/defgradmacroactual(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 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 + +!------------------------- +!RL +!------------------------- end program mpie_spectral @@ -424,3 +954,110 @@ subroutine quit(id) stop end subroutine + + SUBROUTINE fourn(data,nn,ndim,isign) +! +! FROM NUMERICAL RECIPES IN F77 (FIXED FORMAT), +! CONVERTED INTO FREE FORMAT (RL @ MPIE, JUNE 2010) +! + INTEGER isign,ndim,nn(ndim) + REAL data(*) +! 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 + DOUBLE PRECISION theta,wi,wpi,wpr,wr,wtemp + 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 + 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 +!13 continue + enddo ! 12 + enddo ! 13 + endif + 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 +! goto 1 +! endif + enddo ! do while (if 1) + + i2rev=i2rev+ibit +!14 continue + enddo ! 14 + 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 + + 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 +!17 continue + enddo ! 17 + ifp1=ifp2 +! goto 2 +! endif + enddo ! do while (if 2) + + nprev=n*nprev +!18 continue + enddo ! 18 + return + END subroutine