diff --git a/code/IO.f90 b/code/IO.f90 index 59632cd0a..0038d767d 100644 --- a/code/IO.f90 +++ b/code/IO.f90 @@ -926,6 +926,24 @@ endfunction character(len=80) msg select case (ID) + case (40) + msg = 'path rectification error' + case (41) + msg = 'path too long' + case (42) + msg = 'missing descriptive information in spectral mesh' + case (43) + msg = 'resolution error in spectral mesh' + case (44) + msg = 'dimension error in spectral mesh' + case (45) + msg = 'error opening spectral loadcase' + case (46) + msg = 'missing parameter in spectral loadcase' + case (47) + msg = 'mask consistency violated in spectral loadcase' + case (48) + msg = 'stress symmetry violated in spectral loadcase' case (50) msg = 'Error writing constitutive output description' case (100) @@ -933,11 +951,7 @@ endfunction case (101) msg = 'Cannot open input file' case (102) - msg = 'missing descriptive information in spectral config file' - case (103) - msg = 'resolution error in spectral config file' - case (104) - msg = 'dimension error in spectral config file' + msg = 'argument count error (mesh and loadcase) for mpie_spectral' case (105) msg = 'Error reading from ODF file' case (110) diff --git a/code/mesh.f90 b/code/mesh.f90 index 44d87b14f..339a2f2bb 100644 --- a/code/mesh.f90 +++ b/code/mesh.f90 @@ -2137,9 +2137,9 @@ subroutine mesh_marc_count_cpSizes (unit) ! --- sanity checks --- - if (.not. gotDimension .or. .not. gotResolution) call IO_error(102) - if (a < 2 .or. b < 2 .or. c < 2) call IO_error(103) - if (x <= 0.0_pReal .or. y <= 0.0_pReal .or. z <= 0.0_pReal) call IO_error(104) + if (.not. gotDimension .or. .not. gotResolution) call IO_error(42) + if (a < 2 .or. b < 2 .or. c < 2) call IO_error(43) + if (x <= 0.0_pReal .or. y <= 0.0_pReal .or. z <= 0.0_pReal) call IO_error(44) forall (n = 0:mesh_Nnodes-1) mesh_node(1,n+1) = x * dble(mod(n,a) / (a-1.0_pReal)) diff --git a/code/mpie_spectral.f90 b/code/mpie_spectral.f90 index 60562ed81..054ac2bfb 100644 --- a/code/mpie_spectral.f90 +++ b/code/mpie_spectral.f90 @@ -26,290 +26,331 @@ include "prec.f90" ! uses nothing else MODULE mpie_interface - -character(len=64), parameter :: FEsolver = 'Spectral' -character(len=5), parameter :: InputFileExtension = '.mesh' + use prec, only: pInt, pReal + character(len=64), parameter :: FEsolver = 'Spectral' + character(len=5), parameter :: InputFileExtension = '.mesh' CONTAINS -subroutine mpie_interface_init +!******************************************************************** +! initialize interface module +! +!******************************************************************** +subroutine mpie_interface_init() + write(6,*) write(6,*) '<<<+- mpie_spectral init -+>>>' write(6,*) '$Id$' write(6,*) + return end subroutine +!******************************************************************** +! extract working directory from loadcase file +! possibly based on current working dir +!******************************************************************** function getSolverWorkingDirectoryName() - implicit none - character(1024) cwd,outname,getSolverWorkingDirectoryName - character(len=*), parameter :: pathSep = achar(47)//achar(92) ! /, \ - print *, 'start of func' + implicit none + + character(len=1024) cwd,outname,getSolverWorkingDirectoryName + character(len=*), parameter :: pathSep = achar(47)//achar(92) ! /, \ call getarg(2,outname) ! path to loadFile if (scan(outname,pathSep) == 1) then ! absolute path given as command line argument getSolverWorkingDirectoryName = outname(1:scan(outname,pathSep,back=.true.)) - print *, 'abs',scan(outname,pathSep,back=.true.) else call getcwd(cwd) getSolverWorkingDirectoryName = trim(cwd)//'/'//outname(1:scan(outname,pathSep,back=.true.)) - print *, 'rel',scan(outname,pathSep,back=.true.),getSolverWorkingDirectoryName endif + getSolverWorkingDirectoryName = rectifyPath(getSolverWorkingDirectoryName) + return -! getSolverWorkingDirectoryName = rectifyPath(getSolverWorkingdirectory) + end function +!******************************************************************** +! basename of meshfile from command line arguments +! +!******************************************************************** function getSolverJobName() - use prec + + use prec, only: pInt + implicit none - character(1024) getSolverJobName, outName + character(1024) getSolverJobName, outName, cwd character(len=*), parameter :: pathSep = achar(47)//achar(92) ! /, \ - integer(pInt) extPos + integer(pInt) posExt,posSep + + getSolverJobName = '' - getSolverJobName='' call getarg(1,outName) - extPos = len_trim(outName)-4 - getSolverJobName = outName(scan(outName,pathSep,back=.true.)+1:extPos) -! write(6,*) 'getSolverJobName', 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 + call getcwd(cwd) + getSolverJobName = makeRelativePath(getSolverWorkingDirectoryName(),& + rectifyPath(trim(cwd)//'/'//getSolverJobName)) + endif + + return end function -!function removes ../ and ./ in Path -function rectifyPath(Path) -implicit none -character(len=1024) Path, rectifyPath -integer i,j,k,l -!remove ./ from path -l = min(1024,len_trim(Path)) -rectifyPath = path -do i=l,2,-1 +!******************************************************************** +! realive path of loadcase from command line arguments +! +!******************************************************************** +function getLoadcaseName() + + use prec, only: pInt + + implicit none + + character(len=1024) getLoadcaseName, outName, cwd + character(len=*), parameter :: pathSep = achar(47)//achar(92) ! /, \ + integer(pInt) posExt,posSep + + call getarg(2,getLoadcaseName) + + if (scan(getLoadcaseName,pathSep) /= 1) then ! relative path given as command line argument + call getcwd(cwd) + getLoadcaseName = rectifyPath(trim(cwd)//'/'//getLoadcaseName) + endif + + getLoadcaseName = makeRelativePath(getSolverWorkingDirectoryName(),& + getLoadcaseName) + + return +end function + + +!******************************************************************** +! remove ../ and ./ from path +! +!******************************************************************** +function rectifyPath(path) + + use prec, only: pInt + + implicit none + + character(len=*) path + character(len=len_trim(path)) rectifyPath + integer(pInt) i,j,k,l + + !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) /= '.' ) & rectifyPath(i-1:l) = rectifyPath(i+1:l)//' ' -end do + end do -!remove ../ from rectifyPath and change directory -i=0 -do while (i<=len_trim(rectifyPath)) - k=0 - j=0 - if(rectifyPath(len_trim(rectifyPath)-2-i:len_trim(rectifyPath)-i)=='../') then !search for ../ (directory down) - j=i - k=i - if(len_trim(rectifyPath)-5-i<=0) print*, 'enter error message here' !invalid rectifyPath (below root) - do while(rectifyPath(len_trim(rectifyPath)-j-5:len_trim(rectifyPath)-j-3)=='../') !search for more ../ in front of first hit - j=j+3 - end do - j=((j-i)/3+1)*2 !calculate number of ../ - do while(j/=0) !find position to break rectifyPath - i=i+1 - if(rectifyPath(len_trim(rectifyPath)-i:len_trim(rectifyPath)-i)=='/') j=j-1 - end do - if(i>len_trim(rectifyPath)) print*, 'enter error message here' !invalid path (below root) - rectifyPath(len_trim(rectifyPath)-i:len_trim(rectifyPath))=rectifyPath(len_trim(rectifyPath)-k:len_trim(rectifyPath)+i-k) - i=k-1 - end if - i=i+1 -end do -end function rectifyPath + !remove ../ and corresponding directory from rectifyPath + l = len_trim(rectifyPath) + + i = index(rectifyPath(i:l),'../') + j = 0_pInt + + do while (i > j) + 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 + return + + end function rectifyPath - - -! make out of two absolute Paths (a,b) relative Path from a to b +!******************************************************************** +! relative path from absolute a to absolute b +! +!******************************************************************** function makeRelativePath(a,b) + + use prec, only: pInt + implicit none - character (len=1024) :: makeRelativePath,a,b - integer i,posLastCommonSlash,remainingSlashes + + character (len=*) :: a,b + character (len=1024) :: makeRelativePath + integer(pInt) i,posLastCommonSlash,remainingSlashes posLastCommonSlash = 0 remainingSlashes = 0 - do i = 1,min(len_trim(a),len_trim(b)) + + 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 enddo do i = posLastCommonSlash+1,len_trim(a) if (a(i:i) == '/') remainingSlashes = remainingSlashes + 1 enddo - makeRelativePath=repeat('../',remainingSlashes)//b(posLastCommonSlash+1:len_trim(b)) + makeRelativePath = repeat('../',remainingSlashes)//b(posLastCommonSlash+1:len_trim(b)) + + return end function makeRelativePath END MODULE - include "IO.f90" ! uses prec - include "numerics.f90" ! uses prec, IO - include "math.f90" ! uses prec, numerics - include "debug.f90" ! uses prec, numerics - include "FEsolving.f90" ! uses prec, IO - include "mesh.f90" ! uses prec, math, IO, FEsolving - include "material.f90" ! uses prec, math, IO, mesh - include "lattice.f90" ! uses prec, math, IO, material - include "constitutive_phenopowerlaw.f90" ! uses prec, math, IO, lattice, material, debug - include "constitutive_j2.f90" ! uses prec, math, IO, lattice, material, debug - include "constitutive_dislotwin.f90" ! uses prec, math, IO, lattice, material, debug - include "constitutive_nonlocal.f90" ! uses prec, math, IO, lattice, material, debug - include "constitutive.f90" ! uses prec, IO, math, lattice, mesh, debug - include "crystallite.f90" ! uses prec, math, IO, numerics - include "homogenization_isostrain.f90" ! uses prec, math, IO, - include "homogenization_RGC.f90" ! uses prec, math, IO, numerics, mesh: added <<>> - include "homogenization.f90" ! uses prec, math, IO, numerics - include "CPFEM.f90" ! uses prec, math, IO, numerics, debug, FEsolving, mesh, lattice, constitutive, crystallite + + +include "IO.f90" ! uses prec +include "numerics.f90" ! uses prec, IO +include "math.f90" ! uses prec, numerics +include "debug.f90" ! uses prec, numerics +include "FEsolving.f90" ! uses prec, IO +include "mesh.f90" ! uses prec, math, IO, FEsolving +include "material.f90" ! uses prec, math, IO, mesh +include "lattice.f90" ! uses prec, math, IO, material +include "constitutive_phenopowerlaw.f90" ! uses prec, math, IO, lattice, material, debug +include "constitutive_j2.f90" ! uses prec, math, IO, lattice, material, debug +include "constitutive_dislotwin.f90" ! uses prec, math, IO, lattice, material, debug +include "constitutive_nonlocal.f90" ! uses prec, math, IO, lattice, material, debug +include "constitutive.f90" ! uses prec, IO, math, lattice, mesh, debug +include "crystallite.f90" ! uses prec, math, IO, numerics +include "homogenization_isostrain.f90" ! uses prec, math, IO, +include "homogenization_RGC.f90" ! uses prec, math, IO, numerics, mesh: added <<>> +include "homogenization.f90" ! uses prec, math, IO, numerics +include "CPFEM.f90" ! uses prec, math, IO, numerics, debug, FEsolving, mesh, lattice, constitutive, crystallite +!******************************************************************** program mpie_spectral +!******************************************************************** - use prec use mpie_interface - + use prec, only: pInt, pReal + use IO + use math, only: math_transpose3x3 + implicit none - character(len=1024) path, line - integer(pInt), parameter :: maxNchunks = 24 ! 4 identifiers, 18 values for the matrices and 2 scalars + + real(pReal), dimension (:,:,:), allocatable :: velocityGrad, & + stress ! velocity gradient and stress BC + real(pReal), dimension(:), allocatable :: timeIncrement ! length of increment + integer(pInt), dimension(:), allocatable :: steps ! number of steps + logical, dimension(:,:,:,:), allocatable :: mask ! BC mask + + character(len=1024) path,line + logical, dimension(9) :: maskvector + integer(pInt), parameter :: maxNchunks = 24 ! 4 identifiers, 18 values for the matrices and 2 scalars integer(pInt), dimension (1+maxNchunks*2) :: pos - real(pReal), dimension (:,:), allocatable :: l, s ! velocity gradient and stress BC - real(pReal), dimension(:), allocatable :: t, n ! length of time and step number - character, dimension(:,:), allocatable :: mask ! BC mask - integer(pInt) unit, N_l, N_s, N_t, N_n, i, j, k ! numbers of identifiers, loop variables - + real(pReal), dimension(9) :: valuevector + integer(pInt) unit, N_l, N_s, N_t, N_n, N, i,j,k,l ! numbers of identifiers, loop variables + + if (IargC() < 2) call IO_error(102) + call mpie_interface_init() - if (IargC() < 2) then - print *,'buh' - else - path = rectifyPath(getSolverWorkingDirectoryName()) - endif - -! initialize variables -unit=2 -N_l=0 -N_s=0 -N_t=0 -N_n=0 -if (IO_open_file(unit,path)) rewind(unit) ! open file -do - read(unit,'(a1024)',END=101) line - if (IO_isBlank(line)) cycle ! skip empty lines - pos = IO_stringPos(line,maxNchunks) - do i = 1,maxNchunks,1 - select case (IO_lc(IO_stringValue(line,pos,i))) - case('l') - N_l=N_l+1 - case('s') - N_s=N_s+1 - case('t') - N_t=N_t+1 - case('n') - N_n=N_n+1 - 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)) then ! sanity check - print*, 'insert error message code here' ! error message for incomplete input file - end if -enddo + path = getLoadcaseName() + maskvector = '' + unit = 234_pInt + N_l = 0_pInt + N_s = 0_pInt + N_t = 0_pInt + N_n = 0_pInt + + if (.not. IO_open_file(unit,path)) call IO_error(45,ext_msg=path) + + rewind(unit) + do + read(unit,'(a1024)',END=101) line + if (IO_isBlank(line)) cycle ! skip empty lines + pos = IO_stringPos(line,maxNchunks) + do i = 1,maxNchunks,1 + select case (IO_lc(IO_stringValue(line,pos,i))) + case('l','velocitygrad') + N_l = N_l+1 + case('s','stress') + N_s = N_s+1 + case('t','time','delta') + N_t = N_t+1 + case('n','incs','increments','steps') + N_n = N_n+1 + 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 + + enddo ! allocate memory depending on lines in input file -101 allocate (l(9,N_l)) -allocate (s(9,N_s)) -allocate (mask(9,N_s)) -allocate (t(N_t)) -allocate (n(N_n)) +101 N = N_l + allocate (velocityGrad(3,3,N)); velocityGrad = 0.0_pReal + allocate (stress(3,3,N)); stress = 0.0_pReal + allocate (mask(3,3,2,N)); mask = .false. + allocate (timeIncrement(N)); timeIncrement = 0.0_pReal + allocate (steps(N)); steps = 0_pInt -! initialize variables -do i=1,9 - do j=1,N_l - mask(i,j)='x' - end do -end do -do i=1,9 - do j=1,N_l - l(i,j)=0 - end do -end do -do i=1,9 - do j=1,N_s - s(i,j)=0 - end do -end do -i=0 -j=0 - -rewind(unit) -do - read(unit,'(a1024)',END=100) line - if (IO_isBlank(line)) cycle ! build BC mask from input file - j=j+1 - pos = IO_stringPos(line,maxNchunks) - do i = 1,maxNchunks,2 - select case (IO_lc(IO_stringValue(line,pos,i))) - case('l') - do k=1,9 - if((IO_lc(IO_stringValue(line,pos,i+k)))/='-') mask(k,j)='l' - end do - if(((mask(2,j)/='x')).and.((mask(4,j)=='x'))& ! if one non-diagonal element is defined, the - &.or.((mask(4,j)/='x')).and.((mask(2,j)=='x'))& ! correspondig one should not be empty - &.or.((mask(3,j)/='x')).and.((mask(7,j)=='x'))& - &.or.((mask(7,j)/='x')).and.((mask(3,j)=='x'))& - &.or.((mask(6,j)/='x')).and.((mask(8,j)=='x'))& - &.or.((mask(8,j)/='x')).and.((mask(6,j)=='x'))) print*, 'enter error message here' - case('s') - do k=1,9 - if((IO_lc(IO_stringValue(line,pos,i+k)))/='-') then - if(mask(k,j)=='l') then - print*, 'enter error message here' ! stress and velocity gradient bc at the same place - else - mask(k,j)='s' - end if - end if - end do - end select - enddo - enddo -100 rewind(unit) - -do i=1,9 - do j=1,N_l - if(mask(i,j)=='x') print*,'enter error message here' !check if sufficient Nr. of BCs are found - end do -end do - -j=0 -i=0 -do + rewind(unit) + j = 0_pInt + do read(unit,'(a1024)',END=200) line if (IO_isBlank(line)) cycle ! skip empty lines j=j+1 pos = IO_stringPos(line,maxNchunks) - do i = 1,maxNchunks,2 - select case (IO_lc(IO_stringValue(line,pos,i))) - case('l') - do k=1,9 - if(mask(k,j)=='l') L(k,j) = IO_floatValue(line,pos,i+k) ! assign values for the velocity gradient matrix - end do - case('s') - do k=1,9 - if(mask(k,j)=='s') then - s(k,j) = IO_floatValue(line,pos,i+k) ! assign values for the stress BC - select case(k) - case(4) - if(s(4,j)/=s(2,j)) print*, 'enter error code here' !non-symmetric stress BC - end select - else - end if - end do - case('t') ! assign the scalars - t(j) = IO_floatValue(line,pos,i+1) - case('n') - n(j) = IO_floatValue(line,pos,i+1) - end select - enddo + do i = 1,maxNchunks,2 + select case (IO_lc(IO_stringValue(line,pos,i))) + case('l','velocitygrad') + valuevector = 0.0_pReal + forall (k = 1:9) maskvector(k) = IO_stringValue(line,pos,i+k) /= '#' + do k = 1,9 + if (maskvector(k)) valuevector(k) = IO_floatValue(line,pos,i+k) ! assign values for the velocity gradient matrix + enddo + mask(:,:,1,j) = reshape(maskvector,(/3,3/)) + velocityGrad(:,:,j) = reshape(valuevector,(/3,3/)) + case('s','stress') + valuevector = 0.0_pReal + forall (k = 1:9) maskvector(k) = IO_stringValue(line,pos,i+k) /= '#' + do k = 1,9 + if (maskvector(k)) valuevector(k) = IO_floatValue(line,pos,i+k) ! assign values for the stress matrix + enddo + mask(:,:,2,j) = reshape(maskvector,(/3,3/)) + stress(:,:,j) = reshape(valuevector,(/3,3/)) + case('t','time','delta') ! increment time + timeIncrement(j) = IO_floatValue(line,pos,i+1) + case('n','incs','increments','steps') ! steps + steps(j) = IO_intValue(line,pos,i+1) + end select enddo -200 end program mpie_spectral + enddo +200 close(unit) + ! consistency checks + do j = 1,N + if (any(mask(:,:,1,j) == mask(:,:,2,j))) & + call IO_error(47,j) ! mask consistency + if (any(math_transpose3x3(stress(:,:,j)) + stress(:,:,j) /= 2.0_pReal * stress(:,:,j))) & + call IO_error(48,j) ! stress symmetry + + print '(a,/,3(3(f12.6,x)/))','L',velocityGrad(:,:,j) + print '(a,/,3(3(f12.6,x)/))','stress',stress(:,:,j) + print '(a,/,3(3(l,x)/))','mask',mask(:,:,1,j) + print *,'time',timeIncrement(j) + print *,'incs',steps(j) + print *, '' + enddo + +end program mpie_spectral + + +!******************************************************************** +! quit subroutine to satisfy IO_error +! +!******************************************************************** subroutine quit(id) use prec - implicit none integer(pInt) id