From 06be437bc9302406466e2ef86a6b7a97bd3af8f1 Mon Sep 17 00:00:00 2001 From: Krishna Komerla Date: Tue, 19 Jun 2012 13:31:15 +0000 Subject: [PATCH] added minRes to regridding function and writing out of new geometry file updated f2py wrapper to enable the use of init functions. added 2 new error messages to io --- code/DAMASK_run.py | 64 +++++++++++++++++ code/IO.f90 | 4 ++ code/damask.core.pyf | 48 ++++++++++--- code/mesh.f90 | 159 +++++++++++++++++++++++++++++++------------ 4 files changed, 219 insertions(+), 56 deletions(-) create mode 100644 code/DAMASK_run.py diff --git a/code/DAMASK_run.py b/code/DAMASK_run.py new file mode 100644 index 000000000..4c529f6d1 --- /dev/null +++ b/code/DAMASK_run.py @@ -0,0 +1,64 @@ +#!/usr/bin/env python + +import numpy,os,damask,string,sys,subprocess,re +from optparse import OptionParser, Option + +# ----------------------------- +class extendableOption(Option): +# ----------------------------- +# used for definition of new option parser action 'extend', which enables to take multiple option arguments +# taken from online tutorial http://docs.python.org/library/optparse.html + + ACTIONS = Option.ACTIONS + ("extend",) + STORE_ACTIONS = Option.STORE_ACTIONS + ("extend",) + TYPED_ACTIONS = Option.TYPED_ACTIONS + ("extend",) + ALWAYS_TYPED_ACTIONS = Option.ALWAYS_TYPED_ACTIONS + ("extend",) + + def take_action(self, action, dest, opt, value, values, parser): + if action == "extend": + lvalue = value.split(",") + values.ensure_value(dest, []).extend(lvalue) + else: + Option.take_action(self, action, dest, opt, value, values, parser) + +# -------------------------------------------------------------------- +# MAIN +# -------------------------------------------------------------------- + +parser = OptionParser(option_class=extendableOption, usage='%prog options [file[s]]', description = """ +Add column(s) with derived values according to user defined arithmetic operation between column(s). +Columns can be specified either by label or index. + +Example: distance to IP coordinates -- "math.sqrt( #ip.x#**2 + #ip.y#**2 + #ip.z#**2 )" +""" + string.replace('$Id: addCalculation.py 1355 2012-02-23 13:53:12Z MPIE\p.eisenlohr $','\n','\\n') +) + +parser.add_option('-l','--load', '--loadcase', dest='loadcase', type='string', \ + help='PathToLoadFile/NameOfLoadFile.load. "PathToLoadFile" will be the working directory.') +parser.add_option('-g','--geom', '--geometry', dest='geometry', type='string', \ + help='PathToGeomFile/NameOfGeomFile.load.') + +parser.set_defaults(loadcase= '') +parser.set_defaults(geometry= '') + +(options,filenames) = parser.parse_args() + +out=subprocess.Popen(['DAMASK_spectral.exe', '-l', '%s'%options.loadcase, '-g', '%s'%options.geometry],stderr=subprocess.PIPE) +stderr = out.communicate() +stderrLines = string.split(stderr[1],'\n') +exitCode = int(stderrLines[-2]) +print 'exit code', exitCode +if exitCode==2: + start = int(re.search('\d',stderrLines[0]).group(0)) + print 'restart at ', start +#------------regridding---------------------------------------------- +#-------------------------------------------------------------------- +damask.core.prec.prec_init() +damask.core.damask_interface.damask_interface_init(options.loadcase,options.geometry) +damask.core.io.io_init() +damask.core.numerics.numerics_init() +damask.core.debug.debug_init() +damask.core.math.math_init() +damask.core.fesolving.fe_init() +damask.core.mesh.mesh_init(1,1) +damask.core.mesh.mesh_regrid() diff --git a/code/IO.f90 b/code/IO.f90 index ac6fb349a..308bb0258 100644 --- a/code/IO.f90 +++ b/code/IO.f90 @@ -1339,6 +1339,10 @@ subroutine IO_error(error_ID,e,i,g,ext_msg) msg = 'updating of gamma operator not possible if it is pre calculated' case (880_pInt) msg = 'mismatch of microstructure count and a*b*c in geom file' + case (890_pInt) + msg = 'invalid input for regridding a 2D model' + case (891_pInt) + msg = 'invalid input for regridding a 3D model' !* Error messages related to parsing of Abaqus input file diff --git a/code/damask.core.pyf b/code/damask.core.pyf index fc1744729..91c94ae78 100644 --- a/code/damask.core.pyf +++ b/code/damask.core.pyf @@ -14,18 +14,36 @@ python module core ! in interface ! in :core + module prec + subroutine prec_init + end subroutine prec_init + end module prec module damask_interface ! in :damask_interface:DAMASK_spectral_interface.f90 - subroutine damask_interface_init(loadcaseParameterIn,geometryParameterIn) ! in :damask_interface:DAMASK_spectral_interface.f90 character(len=1024), intent(in) :: loadcaseParameterIn character(len=1024), intent(in) :: geometryParameterIn end subroutine damask_interface_init + end module damask_interface - end module damask_interface + module io + subroutine io_init + end subroutine io_init + end module io + + module numerics + subroutine numerics_init + end subroutine numerics_init + end module numerics + module debug + subroutine debug_init + end subroutine debug_init + end module debug - module math ! in :math:math.f90 + module math ! in :math:math.f90 + subroutine math_init + end subroutine math_init subroutine volume_compare(res,geomdim,defgrad,nodes,volume_mismatch) ! in :math:math.f90 ! input variables @@ -167,24 +185,32 @@ python module core ! in real*8, dimension(3), intent(in) :: geomdim integer, intent(in) :: queryPoints integer, intent(in) :: domainPoints - real*8, dimension(spatialDim,queryPoints), intent(in), depend(spatialDim,queryPoints) :: querySet - real*8, dimension(spatialDim,domainPoints), intent(in), depend(spatialDim,domainPoints) :: domainSet + real*8, dimension(spatialDim,queryPoints), intent(in), depend(spatialDim,queryPoints) :: querySet + real*8, dimension(spatialDim,domainPoints), intent(in), depend(spatialDim,domainPoints) :: domainSet ! output variable - integer, dimension(queryPoints), intent(out), depend(queryPoints) :: indices + integer, dimension(queryPoints), intent(out), depend(queryPoints) :: indices ! depending on input - real*8, dimension(spatialDim,(3**spatialDim)*domainPoints), depend(spatialDim,domainPoints) :: domainSetLarge + real*8, dimension(spatialDim,(3**spatialDim)*domainPoints), depend(spatialDim,domainPoints) :: domainSetLarge end subroutine math_nearestNeighborSearch end module math - module mesh ! in :mesh:mesh.f90 + module fesolving + subroutine fe_init + end subroutine fe_init + end module fesolving - function mesh_regrid(resNewInput) ! in :mesh:mesh.f90 + module mesh ! in :mesh:mesh.f90 + subroutine mesh_init(ip,element) + integer, parameter :: ip = 1 + integer, parameter :: element = 1 + end subroutine mesh_init + + function mesh_regrid(resNewInput,minRes) ! in :mesh:mesh.f90 integer, dimension(3) :: mesh_regrid integer, dimension(3), intent(in), optional :: resNewInput + integer, dimension(3), intent(in), optional :: minRes end function mesh_regrid - end module mesh - end interface end python module core diff --git a/code/mesh.f90 b/code/mesh.f90 index a7a0a3d53..0f79e89f5 100644 --- a/code/mesh.f90 +++ b/code/mesh.f90 @@ -777,7 +777,7 @@ function mesh_spectral_getHomogenization(fileUnit) integer(pInt) :: mesh_spectral_getHomogenization character(len=1024) :: line, & keyword - integer(pInt) :: i, j + integer(pInt) :: i logical :: gotHomogenization = .false. integer(pInt) :: myUnit @@ -822,7 +822,7 @@ end function mesh_spectral_getHomogenization !-------------------------------------------------------------------------------------------------- !> @brief Performes a regridding from saved restart information !-------------------------------------------------------------------------------------------------- -function mesh_regrid(resNewInput) +function mesh_regrid(resNewInput,minRes) use prec, only: & pInt, & pReal @@ -830,18 +830,22 @@ function mesh_regrid(resNewInput) getSolverJobName use IO, only: & IO_read_jobBinaryFile ,& - IO_write_jobBinaryFile + IO_write_jobBinaryFile, & + IO_write_jobFile, & + IO_error use math, only: & math_nearestNeighborSearch, & deformed_FFT, & math_mul33x3 - - integer(pInt), dimension(3) :: mesh_regrid + character(len=1024):: formatString, N_Digits integer(pInt), dimension(3), optional, intent(in) :: resNewInput + integer(pInt), dimension(3), optional, intent(in) :: minRes + integer(pInt), dimension(3) :: mesh_regrid, ratio + integer(pInt), dimension(3,2) :: possibleResNew integer(pInt):: maxsize, i, j, k, ielem, Npoints, NpointsNew, spatialDim integer(pInt), dimension(3) :: res, resNew integer(pInt), dimension(:), allocatable :: indices - real(pReal), dimension(3) :: geomdim + real(pReal), dimension(3) :: geomdim,geomdimNew real(pReal), dimension(3,3) :: Favg real(pReal), dimension(:,:), allocatable :: & coordinatesNew, & @@ -866,55 +870,88 @@ function mesh_regrid(resNewInput) integer(pInt), dimension(:,:), allocatable :: & sizeStateConst, sizeStateHomog - res = mesh_spectral_getResolution() - geomdim = mesh_spectral_getDimension() - Npoints = res(1)*res(2)*res(3) - - - - allocate(F(res(1),res(2),res(3),3,3)) - call IO_read_jobBinaryFile(777,'convergedSpectralDefgrad',trim(getSolverJobName()),size(F)) - read (777,rec=1) F - close (777) + res = mesh_spectral_getResolution() + geomdim = mesh_spectral_getDimension() + Npoints = res(1)*res(2)*res(3) +!--------------------------------------------------------- + allocate(F(res(1),res(2),res(3),3,3)) + call IO_read_jobBinaryFile(777,'convergedSpectralDefgrad',trim(getSolverJobName()),size(F)) + read (777,rec=1) F + close (777) ! ----Calculate deformed configuration and average-------- - do i= 1_pInt,3_pInt; do j = 1_pInt,3_pInt - Favg(i,j) = sum(F(1:res(1),1:res(2),1:res(3),i,j)) / Npoints - enddo; enddo - allocate(coordinates(res(1),res(2),res(3),3)) - call deformed_fft(res,geomdim,Favg,1.0_pReal,F,coordinates) - deallocate(F) + do i= 1_pInt,3_pInt; do j = 1_pInt,3_pInt + Favg(i,j) = sum(F(1:res(1),1:res(2),1:res(3),i,j)) / Npoints + enddo; enddo + allocate(coordinates(res(1),res(2),res(3),3)) + call deformed_fft(res,geomdim,Favg,1.0_pReal,F,coordinates) + deallocate(F) ! ----Store coordinates into a linear list-------- - allocate(coordinatesLinear(3,Npoints)) - ielem = 0_pInt - do k=1_pInt,res(3); do j=1_pInt, res(2); do i=1_pInt, res(1) - ielem = ielem + 1_pInt - coordinatesLinear(1:3,ielem) = coordinates(i,j,k,1:3) - enddo; enddo; enddo - deallocate(coordinates) + allocate(coordinatesLinear(3,Npoints)) + ielem = 0_pInt + do k=1_pInt,res(3); do j=1_pInt, res(2); do i=1_pInt, res(1) + ielem = ielem + 1_pInt + coordinatesLinear(1:3,ielem) = coordinates(i,j,k,1:3) + enddo; enddo; enddo + deallocate(coordinates) ! ----For 2D /3D case---------------------------------- if (res(3)== 1_pInt) then spatialDim = 2_pInt + if (minRes(3) > 1_pInt) call IO_error(890_pInt) else spatialDim = 3_pInt + if (minRes(3) <= 1_pInt) call IO_error(891_pInt) endif + + geomdimNew = math_mul33x3(Favg,geomdim) +!---- Automatic detection based on current geom ----------------- -! ----determine/calculate new resolution-------------- if (.not. present(resNewInput)) then - resNew = res + ratio = floor(real(res,pReal) * (geomdimNew/geomdim), pInt) + possibleResNew = 1_pInt + + do i = 1_pInt, spatialDim + if (mod(ratio(i),2) == 0_pInt) then + possibleResNew(i,1:2) = [ratio(i),ratio(i) + 2_pInt] + else + possibleResNew(i,1:2) = [ratio(i)-1_pInt, ratio(i) + 1_pInt] + endif + if (.not.present(minRes)) then + possibleResNew = possibleResNew + else + do k = 1_pInt,3_pInt; do j = 1_pInt,3_pInt + possibleResNew(j,k) = max(possibleResNew(j,k), minRes(j)) + enddo; enddo + endif + enddo + + k = huge(1_pInt) + do i = 0_pInt, 2_pInt**spatialDim - 1 + j = possibleResNew(1,iand(i,1_pInt)/1_pInt + 1_pInt) & + * possibleResNew(2,iand(i,2_pInt)/2_pInt + 1_pInt) & + * possibleResNew(3,iand(i,4_pInt)/4_pInt + 1_pInt) & + - Npoints + if (j < k) then + k = j + resNew =[ possibleResNew(1,iand(i,1_pInt)/1_pInt + 1_pInt), & + possibleResNew(2,iand(i,2_pInt)/2_pInt + 1_pInt), & + possibleResNew(3,iand(i,4_pInt)/4_pInt + 1_pInt) ] + endif + enddo else - if(all(resNewInput==0.0_pReal)) then - mesh_regrid =[ 0,0,0] - return !no automatic determination right now + if(all(resNewInput==0_pInt)) then + resNew = res else resNew = resNewInput endif - endif + endif + + mesh_regrid = resNew NpointsNew = resNew(1)*resNew(2)*resNew(3) - -! ----Calculate regular new coordinates------------ + +! ----Calculate regular new coordinates----------------------------- allocate(coordinatesNew(3,NpointsNew)) ielem = 0_pInt do k=1_pInt,resNew(3); do j=1_pInt, resNew(2); do i=1_pInt, resNew(1) @@ -923,14 +960,47 @@ function mesh_regrid(resNewInput) - geomdim/real(2_pInt*resNew,pReal)) enddo; enddo; enddo -!----- Nearest neighbour search ----------------------- +!----- Nearest neighbour search ------------------------------------ allocate(indices(Npoints)) call math_nearestNeighborSearch(spatialDim, Favg, geomdim, NpointsNew, Npoints, & coordinatesNew, coordinatesLinear, indices) - deallocate(coordinatesNew) - indices = indices / 3_pInt**spatialDim +1_pInt + deallocate(coordinatesNew) + +!----- write out indices-------------------------------------------- + write(N_Digits, '(a)') 1_pInt + int(log10(real(maxval(indices),pReal)),pInt) + N_Digits = adjustl(N_Digits) + formatString = '(I'//trim(N_Digits)//'.'//trim(N_Digits)//')' + + call IO_write_jobFile(777,'idx') ! make it a general open-write file + write(777, '(A)') '1 header' + write(777, '(A)') 'Numbered indices as per the large set' + do i = 1_pInt, Npoints + write(777,trim(formatString),advance='no') indices(i) + if(mod(i,res(1)) == 0_pInt) write(777,'(A)') '' + enddo + close(777) -!--------------------------------------------------------------------------- + do i = 1_pInt, Npoints + indices(i) = indices(i) / 3_pInt**spatialDim +1_pInt ! +1 b'coz index count starts from '0' + enddo + +!------Adjusting the point-to-grain association--------------------- + write(N_Digits, '(a)') 1_pInt + int(log10(real(NpointsNew,pReal)),pInt) + N_Digits = adjustl(N_Digits) + formatString = '(I'//trim(N_Digits)//'.'//trim(N_Digits)//')' + + call IO_write_jobFile(777,'geom') + write(777, '(A)') '3 header' + write(777, '(A, I8, A, I8, A, I8)') 'resolution a ', resNew(1), ' b ', resNew(2), ' c ', resNew(3) + write(777, '(A, g17.10, A, g17.10, A, g17.10)') 'dimension x ', geomdim(1), ' y ', geomdim(2), ' z ', geomdim(3) + write(777, '(A)') 'homogenization 1' + do i = 1_pInt, NpointsNew + write(777,trim(formatString),advance='no') mesh_element(4,indices(i)) + if(mod(i,resNew(1)) == 0_pInt) write(777,'(A)') '' + enddo + close(777) + +!------------------------------------------------------------------- allocate(material_phase (1,1, Npoints)) allocate(material_phaseNew (1,1, NpointsNew)) call IO_read_jobBinaryFile(777,'recordedPhase',trim(getSolverJobName()),size(material_phase)) @@ -1097,10 +1167,9 @@ function mesh_regrid(resNewInput) close (777) deallocate(sizeStateHomog) deallocate(stateHomog) - + deallocate(indices) - mesh_regrid = resNew end function mesh_regrid @@ -1247,7 +1316,7 @@ subroutine mesh_spectral_build_elements(myUnit) integer(pInt), parameter :: maxNchunks = 7_pInt integer(pInt), dimension (1_pInt+2_pInt*maxNchunks) :: myPos integer(pInt), dimension(3) :: res - integer(pInt) :: e, i, j, homog = 0_pInt, headerLength = 0_pInt, maxIntCount + integer(pInt) :: e, i, homog = 0_pInt, headerLength = 0_pInt, maxIntCount integer(pInt), dimension(:), allocatable :: microstructures integer(pInt), dimension(1,1) :: dummySet = 0_pInt character(len=65536) :: line,keyword @@ -2509,7 +2578,7 @@ integer(pInt), intent(in) :: myUnit integer(pInt), parameter :: maxNchunks = 5_pInt integer(pInt), dimension (1+2*maxNchunks) :: myPos integer(pInt) chunk, Nchunks -character(len=300) line, keyword, damaskOption, v +character(len=300) line, damaskOption, v mesh_periodicSurface = .false.