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

Merge branch 'development' of magit1.mpie.de:damask/DAMASK into development

This commit is contained in:
Martin Diehl 2016-04-25 12:35:31 +02:00
parent ce99a071c2 43e74f00cf
commit dadf840c95
62 changed files with 188001 additions and 190901 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

9
.gitattributes vendored Normal file
View File

@ -0,0 +1,9 @@
# from https://help.github.com/articles/dealing-with-line-endings/
#
# always use LF, even if the files are edited on windows, they need to be compiled/used on unix
* text eol=lf
# Denote all files that are truly binary and should not be modified.
*.png binary
*.jpg binary
*.cae binary

9
DAMASK_env.sh
View File

@ -2,7 +2,7 @@
# usage: source DAMASK_env.sh # usage: source DAMASK_env.sh
if [ "$OSTYPE" == "linux-gnu" ] || [ "$OSTYPE" == 'linux' ]; then if [ "$OSTYPE" == "linux-gnu" ] || [ "$OSTYPE" == 'linux' ]; then
DAMASK_ROOT=$(readlink -f "`dirname $BASH_SOURCE`") DAMASK_ROOT=$(python -c "import os,sys; print(os.path.realpath(os.path.expanduser(sys.argv[1])))" "`dirname $BASH_SOURCE`")
else else
[[ "${BASH_SOURCE::1}" == "/" ]] && BASE="" || BASE="`pwd`/" [[ "${BASH_SOURCE::1}" == "/" ]] && BASE="" || BASE="`pwd`/"
STAT=$(stat "`dirname $BASE$BASH_SOURCE`") STAT=$(stat "`dirname $BASE$BASH_SOURCE`")
@ -18,11 +18,11 @@ fi
SOLVER=`which DAMASK_spectral 2>/dev/null` SOLVER=`which DAMASK_spectral 2>/dev/null`
if [ "x$SOLVER" == "x" ]; then if [ "x$SOLVER" == "x" ]; then
export SOLVER='Not found!' SOLVER='Not found!'
fi fi
PROCESSING=`which postResults 2>/dev/null` PROCESSING=`which postResults 2>/dev/null`
if [ "x$PROCESSING" == "x" ]; then if [ "x$PROCESSING" == "x" ]; then
export PROCESSING='Not found!' PROCESSING='Not found!'
fi fi
# according to http://software.intel.com/en-us/forums/topic/501500 # according to http://software.intel.com/en-us/forums/topic/501500
@ -55,7 +55,8 @@ if [ ! -z "$PS1" ]; then
echo "Multithreading DAMASK_NUM_THREADS=$DAMASK_NUM_THREADS" echo "Multithreading DAMASK_NUM_THREADS=$DAMASK_NUM_THREADS"
if [ "x$PETSC_DIR" != "x" ]; then if [ "x$PETSC_DIR" != "x" ]; then
echo "PETSc location $PETSC_DIR" echo "PETSc location $PETSC_DIR"
[[ `readlink -f $PETSC_DIR` == $PETSC_DIR ]] || echo " ~~> "`readlink -f $PETSC_DIR` [[ `python -c "import os,sys; print(os.path.realpath(os.path.expanduser(sys.argv[1])))" "$PETSC_DIR"` == $PETSC_DIR ]] \
|| echo " ~~> "`python -c "import os,sys; print(os.path.realpath(os.path.expanduser(sys.argv[1])))" "$PETSC_DIR"`
fi fi
[[ "x$PETSC_ARCH" != "x" ]] && echo "PETSc architecture $PETSC_ARCH" [[ "x$PETSC_ARCH" != "x" ]] && echo "PETSc architecture $PETSC_ARCH"
echo "MSC.Marc/Mentat $MSC_ROOT" echo "MSC.Marc/Mentat $MSC_ROOT"

2
VERSION
View File

@ -1 +1 @@
v2.0.0-43-ge39441f v2.0.0-160-g26f437b

8
code/.gitattributes vendored Normal file
View File

@ -0,0 +1,8 @@
# from https://help.github.com/articles/dealing-with-line-endings/
#
# always use LF, even if the files are edited on windows, they need to be compiled/used on unix
* text eol=lf
# Denote all files that are truly binary and should not be modified.
*.png binary
*.jpg binary

22
code/DAMASK_spectral.f90
View File

@ -13,7 +13,8 @@ program DAMASK_spectral
pInt, & pInt, &
pLongInt, & pLongInt, &
pReal, & pReal, &
tol_math_check tol_math_check, &
dNeq
use DAMASK_interface, only: & use DAMASK_interface, only: &
DAMASK_interface_init, & DAMASK_interface_init, &
loadCaseFile, & loadCaseFile, &
@ -147,7 +148,9 @@ program DAMASK_spectral
MPI_file_seek, & MPI_file_seek, &
MPI_file_get_position, & MPI_file_get_position, &
MPI_file_write, & MPI_file_write, &
MPI_allreduce MPI_abort, &
MPI_allreduce, &
PETScFinalize
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! init DAMASK (all modules) ! init DAMASK (all modules)
@ -339,7 +342,7 @@ program DAMASK_spectral
reshape(spread(tol_math_check,1,9),[ 3,3]))& reshape(spread(tol_math_check,1,9),[ 3,3]))&
.or. abs(math_det33(loadCases(currentLoadCase)%rotation)) > & .or. abs(math_det33(loadCases(currentLoadCase)%rotation)) > &
1.0_pReal + tol_math_check) errorID = 846_pInt ! given rotation matrix contains strain 1.0_pReal + tol_math_check) errorID = 846_pInt ! given rotation matrix contains strain
if (any(loadCases(currentLoadCase)%rotation /= math_I3)) & if (any(dNeq(loadCases(currentLoadCase)%rotation, math_I3))) &
write(6,'(2x,a,/,3(3(3x,f12.7,1x)/))',advance='no') 'rotation of loadframe:',& write(6,'(2x,a,/,3(3(3x,f12.7,1x)/))',advance='no') 'rotation of loadframe:',&
math_transpose33(loadCases(currentLoadCase)%rotation) math_transpose33(loadCases(currentLoadCase)%rotation)
if (loadCases(currentLoadCase)%time < 0.0_pReal) errorID = 834_pInt ! negative time increment if (loadCases(currentLoadCase)%time < 0.0_pReal) errorID = 834_pInt ! negative time increment
@ -423,17 +426,21 @@ program DAMASK_spectral
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! prepare MPI parallel out (including opening of file) ! prepare MPI parallel out (including opening of file)
allocate(outputSize(worldsize), source = 0_MPI_OFFSET_KIND) allocate(outputSize(worldsize), source = 0_MPI_OFFSET_KIND)
outputSize(worldrank+1) = int(size(materialpoint_results)*pReal,MPI_OFFSET_KIND) outputSize(worldrank+1) = size(materialpoint_results,kind=MPI_OFFSET_KIND)*int(pReal,MPI_OFFSET_KIND)
call MPI_allreduce(MPI_IN_PLACE,outputSize,worldsize,MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr) ! get total output size over each process call MPI_allreduce(MPI_IN_PLACE,outputSize,worldsize,MPI_LONG,MPI_SUM,PETSC_COMM_WORLD,ierr) ! get total output size over each process
if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='MPI_allreduce')
call MPI_file_open(PETSC_COMM_WORLD, & call MPI_file_open(PETSC_COMM_WORLD, &
trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//'.spectralOut', & trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//'.spectralOut', &
MPI_MODE_WRONLY + MPI_MODE_APPEND, & MPI_MODE_WRONLY + MPI_MODE_APPEND, &
MPI_INFO_NULL, & MPI_INFO_NULL, &
resUnit, & resUnit, &
ierr) ierr)
if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='MPI_file_open')
call MPI_file_get_position(resUnit,fileOffset,ierr) ! get offset from header call MPI_file_get_position(resUnit,fileOffset,ierr) ! get offset from header
if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='MPI_file_get_position')
fileOffset = fileOffset + sum(outputSize(1:worldrank)) ! offset of my process in file (header + processes before me) fileOffset = fileOffset + sum(outputSize(1:worldrank)) ! offset of my process in file (header + processes before me)
call MPI_file_seek (resUnit,fileOffset,MPI_SEEK_SET,ierr) call MPI_file_seek (resUnit,fileOffset,MPI_SEEK_SET,ierr)
if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='MPI_file_seek')
if (.not. appendToOutFile) then ! if not restarting, write 0th increment if (.not. appendToOutFile) then ! if not restarting, write 0th increment
do i=1, size(materialpoint_results,3)/(maxByteOut/(materialpoint_sizeResults*pReal))+1 ! slice the output of my process in chunks not exceeding the limit for one output do i=1, size(materialpoint_results,3)/(maxByteOut/(materialpoint_sizeResults*pReal))+1 ! slice the output of my process in chunks not exceeding the limit for one output
@ -443,6 +450,7 @@ program DAMASK_spectral
[(outputIndex(2)-outputIndex(1)+1)*materialpoint_sizeResults]), & [(outputIndex(2)-outputIndex(1)+1)*materialpoint_sizeResults]), &
(outputIndex(2)-outputIndex(1)+1)*materialpoint_sizeResults,& (outputIndex(2)-outputIndex(1)+1)*materialpoint_sizeResults,&
MPI_DOUBLE, MPI_STATUS_IGNORE, ierr) MPI_DOUBLE, MPI_STATUS_IGNORE, ierr)
if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='MPI_file_write')
enddo enddo
fileOffset = fileOffset + sum(outputSize) ! forward to current file position fileOffset = fileOffset + sum(outputSize) ! forward to current file position
if (worldrank == 0) & if (worldrank == 0) &
@ -643,6 +651,7 @@ program DAMASK_spectral
write(6,'(1/,a)') ' ... writing results to file ......................................' write(6,'(1/,a)') ' ... writing results to file ......................................'
call materialpoint_postResults() call materialpoint_postResults()
call MPI_file_seek (resUnit,fileOffset,MPI_SEEK_SET,ierr) call MPI_file_seek (resUnit,fileOffset,MPI_SEEK_SET,ierr)
if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='MPI_file_seek')
do i=1, size(materialpoint_results,3)/(maxByteOut/(materialpoint_sizeResults*pReal))+1 ! slice the output of my process in chunks not exceeding the limit for one output do i=1, size(materialpoint_results,3)/(maxByteOut/(materialpoint_sizeResults*pReal))+1 ! slice the output of my process in chunks not exceeding the limit for one output
outputIndex=int([(i-1_pInt)*((maxByteOut/pReal)/materialpoint_sizeResults)+1_pInt, & outputIndex=int([(i-1_pInt)*((maxByteOut/pReal)/materialpoint_sizeResults)+1_pInt, &
min(i*((maxByteOut/pReal)/materialpoint_sizeResults),size(materialpoint_results,3))],pLongInt) min(i*((maxByteOut/pReal)/materialpoint_sizeResults),size(materialpoint_results,3))],pLongInt)
@ -650,6 +659,7 @@ program DAMASK_spectral
[(outputIndex(2)-outputIndex(1)+1)*materialpoint_sizeResults]), & [(outputIndex(2)-outputIndex(1)+1)*materialpoint_sizeResults]), &
(outputIndex(2)-outputIndex(1)+1)*materialpoint_sizeResults,& (outputIndex(2)-outputIndex(1)+1)*materialpoint_sizeResults,&
MPI_DOUBLE, MPI_STATUS_IGNORE, ierr) MPI_DOUBLE, MPI_STATUS_IGNORE, ierr)
if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='MPI_file_write')
enddo enddo
fileOffset = fileOffset + sum(outputSize) ! forward to current file position fileOffset = fileOffset + sum(outputSize) ! forward to current file position
endif endif
@ -698,7 +708,7 @@ program DAMASK_spectral
enddo enddo
call utilities_destroy() call utilities_destroy()
call PetscFinalize(ierr); CHKERRQ(ierr) call PETScFinalize(ierr); CHKERRQ(ierr)
if (notConvergedCounter > 0_pInt) call quit(3_pInt) ! error if some are not converged if (notConvergedCounter > 0_pInt) call quit(3_pInt) ! error if some are not converged
call quit(0_pInt) ! no complains ;) call quit(0_pInt) ! no complains ;)

545
code/IO.f90
View File

@ -6,10 +6,6 @@
!> @brief input/output functions, partly depending on chosen solver !> @brief input/output functions, partly depending on chosen solver
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module IO module IO
#ifdef HDF
use hdf5, only: &
HID_T
#endif
use prec, only: & use prec, only: &
pInt, & pInt, &
pReal pReal
@ -18,22 +14,8 @@ module IO
private private
character(len=5), parameter, public :: & character(len=5), parameter, public :: &
IO_EOF = '#EOF#' !< end of file string IO_EOF = '#EOF#' !< end of file string
#ifdef HDF
integer(HID_T), public, protected :: tempCoordinates, tempResults
integer(HID_T), private :: resultsFile, tempFile
integer(pInt), private :: currentInc
#endif
public :: & public :: &
#ifdef HDF
HDF5_mappingConstitutive, &
HDF5_mappingHomogenization, &
HDF5_mappingCells, &
HDF5_addGroup ,&
HDF5_forwardResults, &
HDF5_addScalarDataset, &
IO_formatIntToString ,&
#endif
IO_init, & IO_init, &
IO_read, & IO_read, &
IO_checkAndRewind, & IO_checkAndRewind, &
@ -117,9 +99,6 @@ subroutine IO_init
#include "compilation_info.f90" #include "compilation_info.f90"
endif mainProcess endif mainProcess
#ifdef HDF
call HDF5_createJobFile
#endif
end subroutine IO_init end subroutine IO_init
@ -1669,6 +1648,8 @@ subroutine IO_error(error_ID,el,ip,g,ext_msg)
msg = 'unknown filter type selected' msg = 'unknown filter type selected'
case (893_pInt) case (893_pInt)
msg = 'PETSc: SNES_DIVERGED_FNORM_NAN' msg = 'PETSc: SNES_DIVERGED_FNORM_NAN'
case (894_pInt)
msg = 'MPI error'
!------------------------------------------------------------------------------------------------- !-------------------------------------------------------------------------------------------------
! error messages related to parsing of Abaqus input file ! error messages related to parsing of Abaqus input file
@ -1942,526 +1923,4 @@ recursive function abaqus_assembleInputFile(unit1,unit2) result(createSuccess)
end function abaqus_assembleInputFile end function abaqus_assembleInputFile
#endif #endif
#ifdef HDF
!--------------------------------------------------------------------------------------------------
!> @brief creates and initializes HDF5 output files
!--------------------------------------------------------------------------------------------------
subroutine HDF5_createJobFile
use hdf5
use DAMASK_interface, only: &
getSolverWorkingDirectoryName, &
getSolverJobName
implicit none
integer :: hdferr
integer(SIZE_T) :: typeSize
character(len=1024) :: path
integer(HID_T) :: prp_id
integer(SIZE_T), parameter :: increment = 104857600 ! increase temp file in memory in 100MB steps
!--------------------------------------------------------------------------------------------------
! initialize HDF5 library and check if integer and float type size match
call h5open_f(hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_createJobFile: h5open_f')
call h5tget_size_f(H5T_NATIVE_INTEGER,typeSize, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_createJobFile: h5tget_size_f (int)')
if (int(pInt,SIZE_T)/=typeSize) call IO_error(0_pInt,ext_msg='pInt does not match H5T_NATIVE_INTEGER')
call h5tget_size_f(H5T_NATIVE_DOUBLE,typeSize, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_createJobFile: h5tget_size_f (double)')
if (int(pReal,SIZE_T)/=typeSize) call IO_error(0_pInt,ext_msg='pReal does not match H5T_NATIVE_DOUBLE')
!--------------------------------------------------------------------------------------------------
! open file
path = trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//'.'//'DAMASKout'
call h5fcreate_f(path,H5F_ACC_TRUNC_F,resultsFile,hdferr)
if (hdferr < 0) call IO_error(100_pInt,ext_msg=path)
call HDF5_addStringAttribute(resultsFile,'createdBy','$Id$')
!--------------------------------------------------------------------------------------------------
! open temp file
path = trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//'.'//'DAMASKoutTemp'
call h5pcreate_f(H5P_FILE_ACCESS_F, prp_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_createJobFile: h5pcreate_f')
call h5pset_fapl_core_f(prp_id, increment, .false., hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_createJobFile: h5pset_fapl_core_f')
call h5fcreate_f(path,H5F_ACC_TRUNC_F,tempFile,hdferr)
if (hdferr < 0) call IO_error(100_pInt,ext_msg=path)
!--------------------------------------------------------------------------------------------------
! create mapping groups in out file
call HDF5_closeGroup(HDF5_addGroup("mapping"))
call HDF5_closeGroup(HDF5_addGroup("results"))
call HDF5_closeGroup(HDF5_addGroup("coordinates"))
!--------------------------------------------------------------------------------------------------
! create results group in temp file
tempResults = HDF5_addGroup("results",tempFile)
tempCoordinates = HDF5_addGroup("coordinates",tempFile)
end subroutine HDF5_createJobFile
!--------------------------------------------------------------------------------------------------
!> @brief creates and initializes HDF5 output file
!--------------------------------------------------------------------------------------------------
subroutine HDF5_closeJobFile()
use hdf5
implicit none
integer :: hdferr
call h5fclose_f(resultsFile,hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_closeJobFile: h5fclose_f')
end subroutine HDF5_closeJobFile
!--------------------------------------------------------------------------------------------------
!> @brief adds a new group to the results file, or if loc is present at the given location
!--------------------------------------------------------------------------------------------------
integer(HID_T) function HDF5_addGroup(path,loc)
use hdf5
implicit none
character(len=*), intent(in) :: path
integer(HID_T), intent(in),optional :: loc
integer :: hdferr
if (present(loc)) then
call h5gcreate_f(loc, trim(path), HDF5_addGroup, hdferr)
else
call h5gcreate_f(resultsFile, trim(path), HDF5_addGroup, hdferr)
endif
if (hdferr < 0) call IO_error(1_pInt,ext_msg = 'HDF5_addGroup: h5gcreate_f ('//trim(path)//' )')
end function HDF5_addGroup
!--------------------------------------------------------------------------------------------------
!> @brief adds a new group to the results file
!--------------------------------------------------------------------------------------------------
integer(HID_T) function HDF5_openGroup(path)
use hdf5
implicit none
character(len=*), intent(in) :: path
integer :: hdferr
call h5gopen_f(resultsFile, trim(path), HDF5_openGroup, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg = 'HDF5_openGroup: h5gopen_f ('//trim(path)//' )')
end function HDF5_openGroup
!--------------------------------------------------------------------------------------------------
!> @brief closes a group
!--------------------------------------------------------------------------------------------------
subroutine HDF5_closeGroup(ID)
use hdf5
implicit none
integer(HID_T), intent(in) :: ID
integer :: hdferr
call h5gclose_f(ID, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg = 'HDF5_closeGroup: h5gclose_f')
end subroutine HDF5_closeGroup
!--------------------------------------------------------------------------------------------------
!> @brief adds a new group to the results file
!--------------------------------------------------------------------------------------------------
subroutine HDF5_addStringAttribute(entity,attrLabel,attrValue)
use hdf5
implicit none
integer(HID_T), intent(in) :: entity
character(len=*), intent(in) :: attrLabel, attrValue
integer :: hdferr
integer(HID_T) :: attr_id, space_id, type_id
call h5screate_f(H5S_SCALAR_F,space_id,hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_addStringAttribute: h5screate_f')
call h5tcopy_f(H5T_NATIVE_CHARACTER, type_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_addStringAttribute: h5tcopy_f')
call h5tset_size_f(type_id, int(len(trim(attrValue)),HSIZE_T), hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_addStringAttribute: h5tset_size_f')
call h5acreate_f(entity, trim(attrLabel),type_id,space_id,attr_id,hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_addStringAttribute: h5acreate_f')
call h5awrite_f(attr_id, type_id, trim(attrValue), int([1],HSIZE_T), hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_addStringAttribute: h5awrite_f')
call h5aclose_f(attr_id,hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_addStringAttribute: h5aclose_f')
call h5sclose_f(space_id,hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_addStringAttribute: h5sclose_f')
end subroutine HDF5_addStringAttribute
!--------------------------------------------------------------------------------------------------
!> @brief adds the unique mapping from spatial position and constituent ID to results
!--------------------------------------------------------------------------------------------------
subroutine HDF5_mappingConstitutive(mapping)
use hdf5
implicit none
integer(pInt), intent(in), dimension(:,:,:) :: mapping
integer :: hdferr, NmatPoints,Nconstituents
integer(HID_T) :: mapping_id, dtype_id, dset_id, space_id,instance_id,position_id
Nconstituents=size(mapping,1)
NmatPoints=size(mapping,2)
mapping_ID = HDF5_openGroup("mapping")
!--------------------------------------------------------------------------------------------------
! create dataspace
call h5screate_simple_f(2, int([Nconstituents,NmatPoints],HSIZE_T), space_id, hdferr, &
int([Nconstituents,NmatPoints],HSIZE_T))
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingConstitutive')
!--------------------------------------------------------------------------------------------------
! compound type
call h5tcreate_f(H5T_COMPOUND_F, 6_SIZE_T, dtype_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingConstitutive: h5tcreate_f dtype_id')
call h5tinsert_f(dtype_id, "Constitutive Instance", 0_SIZE_T, H5T_STD_U16LE, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingConstitutive: h5tinsert_f 0')
call h5tinsert_f(dtype_id, "Position in Instance Results", 2_SIZE_T, H5T_STD_U32LE, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingConstitutive: h5tinsert_f 2')
!--------------------------------------------------------------------------------------------------
! create Dataset
call h5dcreate_f(mapping_id, "Constitutive", dtype_id, space_id, dset_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingConstitutive')
!--------------------------------------------------------------------------------------------------
! Create memory types (one compound datatype for each member)
call h5tcreate_f(H5T_COMPOUND_F, int(pInt,SIZE_T), instance_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingConstitutive: h5tcreate_f instance_id')
call h5tinsert_f(instance_id, "Constitutive Instance", 0_SIZE_T, H5T_NATIVE_INTEGER, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingConstitutive: h5tinsert_f instance_id')
call h5tcreate_f(H5T_COMPOUND_F, int(pInt,SIZE_T), position_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingConstitutive: h5tcreate_f position_id')
call h5tinsert_f(position_id, "Position in Instance Results", 0_SIZE_T, H5T_NATIVE_INTEGER, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingConstitutive: h5tinsert_f position_id')
!--------------------------------------------------------------------------------------------------
! write data by fields in the datatype. Fields order is not important.
call h5dwrite_f(dset_id, position_id, mapping(1:Nconstituents,1:NmatPoints,1), &
int([Nconstituents, NmatPoints],HSIZE_T), hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingConstitutive: h5dwrite_f position_id')
call h5dwrite_f(dset_id, instance_id, mapping(1:Nconstituents,1:NmatPoints,2), &
int([Nconstituents, NmatPoints],HSIZE_T), hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingConstitutive: h5dwrite_f instance_id')
!--------------------------------------------------------------------------------------------------
!close types, dataspaces
call h5tclose_f(dtype_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingConstitutive: h5tclose_f dtype_id')
call h5tclose_f(position_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingConstitutive: h5tclose_f position_id')
call h5tclose_f(instance_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingConstitutive: h5tclose_f instance_id')
call h5dclose_f(dset_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingConstitutive: h5dclose_f')
call h5sclose_f(space_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingConstitutive: h5sclose_f')
call HDF5_closeGroup(mapping_ID)
end subroutine HDF5_mappingConstitutive
!--------------------------------------------------------------------------------------------------
!> @brief adds the unique mapping from spatial position and constituent ID to results
!--------------------------------------------------------------------------------------------------
subroutine HDF5_mappingCrystallite(mapping)
use hdf5
implicit none
integer(pInt), intent(in), dimension(:,:,:) :: mapping
integer :: hdferr, NmatPoints,Nconstituents
integer(HID_T) :: mapping_id, dtype_id, dset_id, space_id,instance_id,position_id
Nconstituents=size(mapping,1)
NmatPoints=size(mapping,2)
mapping_ID = HDF5_openGroup("mapping")
!--------------------------------------------------------------------------------------------------
! create dataspace
call h5screate_simple_f(2, int([Nconstituents,NmatPoints],HSIZE_T), space_id, hdferr, &
int([Nconstituents,NmatPoints],HSIZE_T))
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite')
!--------------------------------------------------------------------------------------------------
! compound type
call h5tcreate_f(H5T_COMPOUND_F, 6_SIZE_T, dtype_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5tcreate_f dtype_id')
call h5tinsert_f(dtype_id, "Crystallite Instance", 0_SIZE_T, H5T_STD_U16LE, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5tinsert_f 0')
call h5tinsert_f(dtype_id, "Position in Instance Results", 2_SIZE_T, H5T_STD_U32LE, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5tinsert_f 2')
!--------------------------------------------------------------------------------------------------
! create Dataset
call h5dcreate_f(mapping_id, "Crystallite", dtype_id, space_id, dset_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite')
!--------------------------------------------------------------------------------------------------
! Create memory types (one compound datatype for each member)
call h5tcreate_f(H5T_COMPOUND_F, int(pInt,SIZE_T), instance_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5tcreate_f instance_id')
call h5tinsert_f(instance_id, "Crystallite Instance", 0_SIZE_T, H5T_NATIVE_INTEGER, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5tinsert_f instance_id')
call h5tcreate_f(H5T_COMPOUND_F, int(pInt,SIZE_T), position_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5tcreate_f position_id')
call h5tinsert_f(position_id, "Position in Instance Results", 0_SIZE_T, H5T_NATIVE_INTEGER, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5tinsert_f position_id')
!--------------------------------------------------------------------------------------------------
! write data by fields in the datatype. Fields order is not important.
call h5dwrite_f(dset_id, position_id, mapping(1:Nconstituents,1:NmatPoints,1), &
int([Nconstituents, NmatPoints],HSIZE_T), hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5dwrite_f position_id')
call h5dwrite_f(dset_id, instance_id, mapping(1:Nconstituents,1:NmatPoints,2), &
int([Nconstituents, NmatPoints],HSIZE_T), hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5dwrite_f instance_id')
!--------------------------------------------------------------------------------------------------
!close types, dataspaces
call h5tclose_f(dtype_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5tclose_f dtype_id')
call h5tclose_f(position_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5tclose_f position_id')
call h5tclose_f(instance_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5tclose_f instance_id')
call h5dclose_f(dset_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5dclose_f')
call h5sclose_f(space_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5sclose_f')
call HDF5_closeGroup(mapping_ID)
end subroutine HDF5_mappingCrystallite
!--------------------------------------------------------------------------------------------------
!> @brief adds the unique mapping from spatial position to results
!--------------------------------------------------------------------------------------------------
subroutine HDF5_mappingHomogenization(mapping)
use hdf5
implicit none
integer(pInt), intent(in), dimension(:,:) :: mapping
integer :: hdferr, NmatPoints
integer(HID_T) :: mapping_id, dtype_id, dset_id, space_id,instance_id,position_id,elem_id,ip_id
NmatPoints=size(mapping,1)
mapping_ID = HDF5_openGroup("mapping")
!--------------------------------------------------------------------------------------------------
! create dataspace
call h5screate_simple_f(1, int([NmatPoints],HSIZE_T), space_id, hdferr, &
int([NmatPoints],HSIZE_T))
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization')
!--------------------------------------------------------------------------------------------------
! compound type
call h5tcreate_f(H5T_COMPOUND_F, 11_SIZE_T, dtype_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5tcreate_f dtype_id')
call h5tinsert_f(dtype_id, "Homogenization Instance", 0_SIZE_T, H5T_STD_U16LE, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5tinsert_f 0')
call h5tinsert_f(dtype_id, "Position in Instance Results", 2_SIZE_T, H5T_STD_U32LE, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5tinsert_f 2')
call h5tinsert_f(dtype_id, "Element Number", 6_SIZE_T, H5T_STD_U32LE, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5tinsert_f 6')
call h5tinsert_f(dtype_id, "Material Point Number", 10_SIZE_T, H5T_STD_U8LE, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5tinsert_f 10')
!--------------------------------------------------------------------------------------------------
! create Dataset
call h5dcreate_f(mapping_id, "Homogenization", dtype_id, space_id, dset_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization')
!--------------------------------------------------------------------------------------------------
! Create memory types (one compound datatype for each member)
call h5tcreate_f(H5T_COMPOUND_F, int(pInt,SIZE_T), instance_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5tcreate_f instance_id')
call h5tinsert_f(instance_id, "Homogenization Instance", 0_SIZE_T, H5T_NATIVE_INTEGER, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5tinsert_f instance_id')
call h5tcreate_f(H5T_COMPOUND_F, int(pInt,SIZE_T), position_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5tcreate_f position_id')
call h5tinsert_f(position_id, "Position in Instance Results", 0_SIZE_T, H5T_NATIVE_INTEGER, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5tinsert_f position_id')
call h5tcreate_f(H5T_COMPOUND_F, int(pInt,SIZE_T), elem_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5tcreate_f elem_id')
call h5tinsert_f(elem_id, "Element Number", 0_SIZE_T, H5T_NATIVE_INTEGER, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5tinsert_f elem_id')
call h5tcreate_f(H5T_COMPOUND_F, int(pInt,SIZE_T), ip_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5tcreate_f ip_id')
call h5tinsert_f(ip_id, "Material Point Number", 0_SIZE_T, H5T_NATIVE_INTEGER, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5tinsert_f ip_id')
!--------------------------------------------------------------------------------------------------
! write data by fields in the datatype. Fields order is not important.
call h5dwrite_f(dset_id, position_id, mapping(1:NmatPoints,1), &
int([NmatPoints],HSIZE_T), hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5dwrite_f position_id')
call h5dwrite_f(dset_id, instance_id, mapping(1:NmatPoints,2), &
int([NmatPoints],HSIZE_T), hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5dwrite_f position_id')
call h5dwrite_f(dset_id, elem_id, mapping(1:NmatPoints,3), &
int([NmatPoints],HSIZE_T), hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5dwrite_f elem_id')
call h5dwrite_f(dset_id, ip_id, mapping(1:NmatPoints,4), &
int([NmatPoints],HSIZE_T), hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5dwrite_f ip_id')
!--------------------------------------------------------------------------------------------------
!close types, dataspaces
call h5tclose_f(dtype_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5tclose_f dtype_id')
call h5tclose_f(position_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5tclose_f position_id')
call h5tclose_f(instance_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5tclose_f instance_id')
call h5tclose_f(ip_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5tclose_f ip_id')
call h5tclose_f(elem_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5tclose_f elem_id')
call h5dclose_f(dset_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5dclose_f')
call h5sclose_f(space_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5sclose_f')
call HDF5_closeGroup(mapping_ID)
end subroutine HDF5_mappingHomogenization
!--------------------------------------------------------------------------------------------------
!> @brief adds the unique cell to node mapping
!--------------------------------------------------------------------------------------------------
subroutine HDF5_mappingCells(mapping)
use hdf5
implicit none
integer(pInt), intent(in), dimension(:) :: mapping
integer :: hdferr, Nnodes
integer(HID_T) :: mapping_id, dset_id, space_id
Nnodes=size(mapping)
mapping_ID = HDF5_openGroup("mapping")
!--------------------------------------------------------------------------------------------------
! create dataspace
call h5screate_simple_f(1, int([Nnodes],HSIZE_T), space_id, hdferr, &
int([Nnodes],HSIZE_T))
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCells: h5screate_simple_f')
!--------------------------------------------------------------------------------------------------
! create Dataset
call h5dcreate_f(mapping_id, "Cell",H5T_NATIVE_INTEGER, space_id, dset_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCells')
!--------------------------------------------------------------------------------------------------
! write data by fields in the datatype. Fields order is not important.
call h5dwrite_f(dset_id, H5T_NATIVE_INTEGER, mapping, int([Nnodes],HSIZE_T), hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCells: h5dwrite_f instance_id')
!--------------------------------------------------------------------------------------------------
!close types, dataspaces
call h5dclose_f(dset_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingConstitutive: h5dclose_f')
call h5sclose_f(space_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingConstitutive: h5sclose_f')
call HDF5_closeGroup(mapping_ID)
end subroutine HDF5_mappingCells
!--------------------------------------------------------------------------------------------------
!> @brief creates a new scalar dataset in the given group location
!--------------------------------------------------------------------------------------------------
subroutine HDF5_addScalarDataset(group,nnodes,label,SIunit)
use hdf5
implicit none
integer(HID_T), intent(in) :: group
integer(pInt), intent(in) :: nnodes
character(len=*), intent(in) :: SIunit,label
integer :: hdferr
integer(HID_T) :: dset_id, space_id
!--------------------------------------------------------------------------------------------------
! create dataspace
call h5screate_simple_f(1, int([Nnodes],HSIZE_T), space_id, hdferr, &
int([Nnodes],HSIZE_T))
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addScalarDataset: h5screate_simple_f')
!--------------------------------------------------------------------------------------------------
! create Dataset
call h5dcreate_f(group, trim(label),H5T_NATIVE_DOUBLE, space_id, dset_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addScalarDataset: h5dcreate_f')
call HDF5_addStringAttribute(dset_id,'unit',trim(SIunit))
!--------------------------------------------------------------------------------------------------
!close types, dataspaces
call h5dclose_f(dset_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addScalarDataset: h5dclose_f')
call h5sclose_f(space_id, hdferr)
end subroutine HDF5_addScalarDataset
!--------------------------------------------------------------------------------------------------
!> @brief returns nicely formatted string of integer value
!--------------------------------------------------------------------------------------------------
function IO_formatIntToString(myInt)
implicit none
integer(pInt), intent(in) :: myInt
character(len=1_pInt + int(log10(real(myInt)),pInt)) :: IO_formatIntToString
write(IO_formatIntToString,'('//IO_intOut(myInt)//')') myInt
end function
!--------------------------------------------------------------------------------------------------
!> @brief copies the current temp results to the actual results file
!--------------------------------------------------------------------------------------------------
subroutine HDF5_forwardResults
use hdf5
implicit none
integer :: hdferr
integer(HID_T) :: new_loc_id
new_loc_id = HDF5_openGroup("results")
currentInc = currentInc + 1_pInt
call h5ocopy_f(tempFile, 'results', new_loc_id,dst_name=IO_formatIntToString(currentInc), hdferr=hdferr)
if (hdferr < 0_pInt) call IO_error(1_pInt,ext_msg='HDF5_forwardResults: h5ocopy_f')
call HDF5_closeGroup(new_loc_id)
end subroutine HDF5_forwardResults
#endif
end module IO end module IO

5
code/Makefile
View File

@ -309,7 +309,7 @@ KINEMATICS_FILES = \
kinematics_vacancy_strain.o kinematics_hydrogen_strain.o kinematics_vacancy_strain.o kinematics_hydrogen_strain.o
PLASTIC_FILES = \ PLASTIC_FILES = \
plastic_dislotwin.o plastic_disloUCLA.o plastic_isotropic.o plastic_j2.o \ plastic_dislotwin.o plastic_disloUCLA.o plastic_isotropic.o \
plastic_phenopowerlaw.o plastic_titanmod.o plastic_nonlocal.o plastic_none.o \ plastic_phenopowerlaw.o plastic_titanmod.o plastic_nonlocal.o plastic_none.o \
plastic_phenoplus.o plastic_phenoplus.o
@ -579,9 +579,6 @@ plastic_phenoplus.o: plastic_phenoplus.f90 \
plastic_isotropic.o: plastic_isotropic.f90 \ plastic_isotropic.o: plastic_isotropic.f90 \
lattice.o lattice.o
plastic_j2.o: plastic_j2.f90 \
lattice.o
plastic_none.o: plastic_none.f90 \ plastic_none.o: plastic_none.f90 \
lattice.o lattice.o
ifeq "$(F90)" "gfortran" ifeq "$(F90)" "gfortran"

1
code/commercialFEM_fileList.f90
View File

@ -28,7 +28,6 @@
#include "kinematics_hydrogen_strain.f90" #include "kinematics_hydrogen_strain.f90"
#include "plastic_none.f90" #include "plastic_none.f90"
#include "plastic_isotropic.f90" #include "plastic_isotropic.f90"
#include "plastic_j2.f90"
#include "plastic_phenopowerlaw.f90" #include "plastic_phenopowerlaw.f90"
#include "plastic_phenoplus.f90" #include "plastic_phenoplus.f90"
#include "plastic_titanmod.f90" #include "plastic_titanmod.f90"

28
code/constitutive.f90
View File

@ -69,7 +69,6 @@ subroutine constitutive_init()
ELASTICITY_hooke_ID, & ELASTICITY_hooke_ID, &
PLASTICITY_none_ID, & PLASTICITY_none_ID, &
PLASTICITY_isotropic_ID, & PLASTICITY_isotropic_ID, &
PLASTICITY_j2_ID, &
PLASTICITY_phenopowerlaw_ID, & PLASTICITY_phenopowerlaw_ID, &
PLASTICITY_phenoplus_ID, & PLASTICITY_phenoplus_ID, &
PLASTICITY_dislotwin_ID, & PLASTICITY_dislotwin_ID, &
@ -93,7 +92,6 @@ subroutine constitutive_init()
ELASTICITY_HOOKE_label, & ELASTICITY_HOOKE_label, &
PLASTICITY_NONE_label, & PLASTICITY_NONE_label, &
PLASTICITY_ISOTROPIC_label, & PLASTICITY_ISOTROPIC_label, &
PLASTICITY_J2_label, &
PLASTICITY_PHENOPOWERLAW_label, & PLASTICITY_PHENOPOWERLAW_label, &
PLASTICITY_PHENOPLUS_label, & PLASTICITY_PHENOPLUS_label, &
PLASTICITY_DISLOTWIN_label, & PLASTICITY_DISLOTWIN_label, &
@ -114,7 +112,6 @@ subroutine constitutive_init()
use plastic_none use plastic_none
use plastic_isotropic use plastic_isotropic
use plastic_j2
use plastic_phenopowerlaw use plastic_phenopowerlaw
use plastic_phenoplus use plastic_phenoplus
use plastic_dislotwin use plastic_dislotwin
@ -160,7 +157,6 @@ subroutine constitutive_init()
! parse plasticities from config file ! parse plasticities from config file
if (any(phase_plasticity == PLASTICITY_NONE_ID)) call plastic_none_init if (any(phase_plasticity == PLASTICITY_NONE_ID)) call plastic_none_init
if (any(phase_plasticity == PLASTICITY_ISOTROPIC_ID)) call plastic_isotropic_init(FILEUNIT) if (any(phase_plasticity == PLASTICITY_ISOTROPIC_ID)) call plastic_isotropic_init(FILEUNIT)
if (any(phase_plasticity == PLASTICITY_J2_ID)) call plastic_j2_init(FILEUNIT)
if (any(phase_plasticity == PLASTICITY_PHENOPOWERLAW_ID)) call plastic_phenopowerlaw_init(FILEUNIT) if (any(phase_plasticity == PLASTICITY_PHENOPOWERLAW_ID)) call plastic_phenopowerlaw_init(FILEUNIT)
if (any(phase_plasticity == PLASTICITY_PHENOPLUS_ID)) call plastic_phenoplus_init(FILEUNIT) if (any(phase_plasticity == PLASTICITY_PHENOPLUS_ID)) call plastic_phenoplus_init(FILEUNIT)
if (any(phase_plasticity == PLASTICITY_DISLOTWIN_ID)) call plastic_dislotwin_init(FILEUNIT) if (any(phase_plasticity == PLASTICITY_DISLOTWIN_ID)) call plastic_dislotwin_init(FILEUNIT)
@ -217,11 +213,6 @@ subroutine constitutive_init()
thisNoutput => plastic_isotropic_Noutput thisNoutput => plastic_isotropic_Noutput
thisOutput => plastic_isotropic_output thisOutput => plastic_isotropic_output
thisSize => plastic_isotropic_sizePostResult thisSize => plastic_isotropic_sizePostResult
case (PLASTICITY_J2_ID) plasticityType
outputName = PLASTICITY_J2_label
thisNoutput => plastic_j2_Noutput
thisOutput => plastic_j2_output
thisSize => plastic_j2_sizePostResult
case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType
outputName = PLASTICITY_PHENOPOWERLAW_label outputName = PLASTICITY_PHENOPOWERLAW_label
thisNoutput => plastic_phenopowerlaw_Noutput thisNoutput => plastic_phenopowerlaw_Noutput
@ -408,8 +399,6 @@ function constitutive_homogenizedC(ipc,ip,el)
plastic_titanmod_homogenizedC plastic_titanmod_homogenizedC
use plastic_dislotwin, only: & use plastic_dislotwin, only: &
plastic_dislotwin_homogenizedC plastic_dislotwin_homogenizedC
use plastic_disloucla, only: &
plastic_disloucla_homogenizedC
use lattice, only: & use lattice, only: &
lattice_C66 lattice_C66
@ -423,8 +412,6 @@ function constitutive_homogenizedC(ipc,ip,el)
plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el))) plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
case (PLASTICITY_DISLOTWIN_ID) plasticityType case (PLASTICITY_DISLOTWIN_ID) plasticityType
constitutive_homogenizedC = plastic_dislotwin_homogenizedC(ipc,ip,el) constitutive_homogenizedC = plastic_dislotwin_homogenizedC(ipc,ip,el)
case (PLASTICITY_DISLOUCLA_ID) plasticityType
constitutive_homogenizedC = plastic_disloucla_homogenizedC(ipc,ip,el)
case (PLASTICITY_TITANMOD_ID) plasticityType case (PLASTICITY_TITANMOD_ID) plasticityType
constitutive_homogenizedC = plastic_titanmod_homogenizedC (ipc,ip,el) constitutive_homogenizedC = plastic_titanmod_homogenizedC (ipc,ip,el)
case default plasticityType case default plasticityType
@ -513,7 +500,6 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar3333, dLp_dFi3333, Tstar_v
thermalMapping, & thermalMapping, &
PLASTICITY_NONE_ID, & PLASTICITY_NONE_ID, &
PLASTICITY_ISOTROPIC_ID, & PLASTICITY_ISOTROPIC_ID, &
PLASTICITY_J2_ID, &
PLASTICITY_PHENOPOWERLAW_ID, & PLASTICITY_PHENOPOWERLAW_ID, &
PLASTICITY_PHENOPLUS_ID, & PLASTICITY_PHENOPLUS_ID, &
PLASTICITY_DISLOTWIN_ID, & PLASTICITY_DISLOTWIN_ID, &
@ -522,8 +508,6 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar3333, dLp_dFi3333, Tstar_v
PLASTICITY_NONLOCAL_ID PLASTICITY_NONLOCAL_ID
use plastic_isotropic, only: & use plastic_isotropic, only: &
plastic_isotropic_LpAndItsTangent plastic_isotropic_LpAndItsTangent
use plastic_j2, only: &
plastic_j2_LpAndItsTangent
use plastic_phenopowerlaw, only: & use plastic_phenopowerlaw, only: &
plastic_phenopowerlaw_LpAndItsTangent plastic_phenopowerlaw_LpAndItsTangent
use plastic_phenoplus, only: & use plastic_phenoplus, only: &
@ -574,8 +558,6 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar3333, dLp_dFi3333, Tstar_v
dLp_dMstar = 0.0_pReal dLp_dMstar = 0.0_pReal
case (PLASTICITY_ISOTROPIC_ID) plasticityType case (PLASTICITY_ISOTROPIC_ID) plasticityType
call plastic_isotropic_LpAndItsTangent(Lp,dLp_dMstar,Mstar_v,ipc,ip,el) call plastic_isotropic_LpAndItsTangent(Lp,dLp_dMstar,Mstar_v,ipc,ip,el)
case (PLASTICITY_J2_ID) plasticityType
call plastic_j2_LpAndItsTangent(Lp,dLp_dMstar,Mstar_v,ipc,ip,el)
case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType
call plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMstar,Mstar_v,ipc,ip,el) call plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMstar,Mstar_v,ipc,ip,el)
case (PLASTICITY_PHENOPLUS_ID) plasticityType case (PLASTICITY_PHENOPLUS_ID) plasticityType
@ -903,7 +885,6 @@ subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, subdt, subfra
homogenization_maxNgrains, & homogenization_maxNgrains, &
PLASTICITY_none_ID, & PLASTICITY_none_ID, &
PLASTICITY_isotropic_ID, & PLASTICITY_isotropic_ID, &
PLASTICITY_j2_ID, &
PLASTICITY_phenopowerlaw_ID, & PLASTICITY_phenopowerlaw_ID, &
PLASTICITY_phenoplus_ID, & PLASTICITY_phenoplus_ID, &
PLASTICITY_dislotwin_ID, & PLASTICITY_dislotwin_ID, &
@ -916,8 +897,6 @@ subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, subdt, subfra
SOURCE_thermal_externalheat_ID SOURCE_thermal_externalheat_ID
use plastic_isotropic, only: & use plastic_isotropic, only: &
plastic_isotropic_dotState plastic_isotropic_dotState
use plastic_j2, only: &
plastic_j2_dotState
use plastic_phenopowerlaw, only: & use plastic_phenopowerlaw, only: &
plastic_phenopowerlaw_dotState plastic_phenopowerlaw_dotState
use plastic_phenoplus, only: & use plastic_phenoplus, only: &
@ -971,8 +950,6 @@ subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, subdt, subfra
plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el))) plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
case (PLASTICITY_ISOTROPIC_ID) plasticityType case (PLASTICITY_ISOTROPIC_ID) plasticityType
call plastic_isotropic_dotState (Tstar_v,ipc,ip,el) call plastic_isotropic_dotState (Tstar_v,ipc,ip,el)
case (PLASTICITY_J2_ID) plasticityType
call plastic_j2_dotState (Tstar_v,ipc,ip,el)
case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType
call plastic_phenopowerlaw_dotState(Tstar_v,ipc,ip,el) call plastic_phenopowerlaw_dotState(Tstar_v,ipc,ip,el)
case (PLASTICITY_PHENOPLUS_ID) plasticityType case (PLASTICITY_PHENOPLUS_ID) plasticityType
@ -1117,7 +1094,6 @@ function constitutive_postResults(Tstar_v, FeArray, ipc, ip, el)
homogenization_maxNgrains, & homogenization_maxNgrains, &
PLASTICITY_NONE_ID, & PLASTICITY_NONE_ID, &
PLASTICITY_ISOTROPIC_ID, & PLASTICITY_ISOTROPIC_ID, &
PLASTICITY_J2_ID, &
PLASTICITY_PHENOPOWERLAW_ID, & PLASTICITY_PHENOPOWERLAW_ID, &
PLASTICITY_PHENOPLUS_ID, & PLASTICITY_PHENOPLUS_ID, &
PLASTICITY_DISLOTWIN_ID, & PLASTICITY_DISLOTWIN_ID, &
@ -1130,8 +1106,6 @@ function constitutive_postResults(Tstar_v, FeArray, ipc, ip, el)
SOURCE_damage_anisoDuctile_ID SOURCE_damage_anisoDuctile_ID
use plastic_isotropic, only: & use plastic_isotropic, only: &
plastic_isotropic_postResults plastic_isotropic_postResults
use plastic_j2, only: &
plastic_j2_postResults
use plastic_phenopowerlaw, only: & use plastic_phenopowerlaw, only: &
plastic_phenopowerlaw_postResults plastic_phenopowerlaw_postResults
use plastic_phenoplus, only: & use plastic_phenoplus, only: &
@ -1185,8 +1159,6 @@ function constitutive_postResults(Tstar_v, FeArray, ipc, ip, el)
constitutive_postResults(startPos:endPos) = plastic_titanmod_postResults(ipc,ip,el) constitutive_postResults(startPos:endPos) = plastic_titanmod_postResults(ipc,ip,el)
case (PLASTICITY_ISOTROPIC_ID) plasticityType case (PLASTICITY_ISOTROPIC_ID) plasticityType
constitutive_postResults(startPos:endPos) = plastic_isotropic_postResults(Tstar_v,ipc,ip,el) constitutive_postResults(startPos:endPos) = plastic_isotropic_postResults(Tstar_v,ipc,ip,el)
case (PLASTICITY_J2_ID) plasticityType
constitutive_postResults(startPos:endPos) = plastic_j2_postResults(Tstar_v,ipc,ip,el)
case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType
constitutive_postResults(startPos:endPos) = & constitutive_postResults(startPos:endPos) = &
plastic_phenopowerlaw_postResults(Tstar_v,ipc,ip,el) plastic_phenopowerlaw_postResults(Tstar_v,ipc,ip,el)

23
code/homogenization.f90
View File

@ -71,12 +71,6 @@ contains
!> @brief module initialization !> @brief module initialization
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine homogenization_init subroutine homogenization_init
#ifdef HDF
use hdf5, only: &
HID_T
use IO, only : &
HDF5_mappingHomogenization
#endif
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use math, only: & use math, only: &
math_I3 math_I3
@ -131,12 +125,6 @@ subroutine homogenization_init
character(len=64), dimension(:,:), pointer :: thisOutput character(len=64), dimension(:,:), pointer :: thisOutput
character(len=32) :: outputName !< name of output, intermediate fix until HDF5 output is ready character(len=32) :: outputName !< name of output, intermediate fix until HDF5 output is ready
logical :: knownHomogenization, knownThermal, knownDamage, knownVacancyflux, knownPorosity, knownHydrogenflux logical :: knownHomogenization, knownThermal, knownDamage, knownVacancyflux, knownPorosity, knownHydrogenflux
#ifdef HDF
integer(pInt), dimension(:,:), allocatable :: mapping
integer(pInt), dimension(:), allocatable :: InstancePosition
allocate(mapping(mesh_ncpelems,4),source=0_pInt)
allocate(InstancePosition(material_Nhomogenization),source=0_pInt)
#endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -396,17 +384,6 @@ subroutine homogenization_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate and initialize global state and postresutls variables ! allocate and initialize global state and postresutls variables
#ifdef HDF
elementLooping: do e = 1,mesh_NcpElems
myInstance = homogenization_typeInstance(mesh_element(3,e))
IpLooping: do i = 1,FE_Nips(FE_geomtype(mesh_element(2,e)))
InstancePosition(myInstance) = InstancePosition(myInstance)+1_pInt
mapping(e,1:4) = [instancePosition(myinstance),myinstance,e,i]
enddo IpLooping
enddo elementLooping
call HDF5_mappingHomogenization(mapping)
#endif
homogenization_maxSizePostResults = 0_pInt homogenization_maxSizePostResults = 0_pInt
thermal_maxSizePostResults = 0_pInt thermal_maxSizePostResults = 0_pInt
damage_maxSizePostResults = 0_pInt damage_maxSizePostResults = 0_pInt

27
code/material.f90
View File

@ -24,7 +24,6 @@ module material
ELASTICITY_hooke_label = 'hooke', & ELASTICITY_hooke_label = 'hooke', &
PLASTICITY_none_label = 'none', & PLASTICITY_none_label = 'none', &
PLASTICITY_isotropic_label = 'isotropic', & PLASTICITY_isotropic_label = 'isotropic', &
PLASTICITY_j2_label = 'j2', &
PLASTICITY_phenopowerlaw_label = 'phenopowerlaw', & PLASTICITY_phenopowerlaw_label = 'phenopowerlaw', &
PLASTICITY_phenoplus_label = 'phenoplus', & PLASTICITY_phenoplus_label = 'phenoplus', &
PLASTICITY_dislotwin_label = 'dislotwin', & PLASTICITY_dislotwin_label = 'dislotwin', &
@ -74,7 +73,6 @@ module material
enumerator :: PLASTICITY_undefined_ID, & enumerator :: PLASTICITY_undefined_ID, &
PLASTICITY_none_ID, & PLASTICITY_none_ID, &
PLASTICITY_isotropic_ID, & PLASTICITY_isotropic_ID, &
PLASTICITY_j2_ID, &
PLASTICITY_phenopowerlaw_ID, & PLASTICITY_phenopowerlaw_ID, &
PLASTICITY_phenoplus_ID, & PLASTICITY_phenoplus_ID, &
PLASTICITY_dislotwin_ID, & PLASTICITY_dislotwin_ID, &
@ -313,7 +311,6 @@ module material
ELASTICITY_hooke_ID ,& ELASTICITY_hooke_ID ,&
PLASTICITY_none_ID, & PLASTICITY_none_ID, &
PLASTICITY_isotropic_ID, & PLASTICITY_isotropic_ID, &
PLASTICITY_J2_ID, &
PLASTICITY_phenopowerlaw_ID, & PLASTICITY_phenopowerlaw_ID, &
PLASTICITY_phenoplus_ID, & PLASTICITY_phenoplus_ID, &
PLASTICITY_dislotwin_ID, & PLASTICITY_dislotwin_ID, &
@ -351,9 +348,6 @@ module material
HYDROGENFLUX_cahnhilliard_ID, & HYDROGENFLUX_cahnhilliard_ID, &
HOMOGENIZATION_none_ID, & HOMOGENIZATION_none_ID, &
HOMOGENIZATION_isostrain_ID, & HOMOGENIZATION_isostrain_ID, &
#ifdef HDF
material_NconstituentsPhase, &
#endif
HOMOGENIZATION_RGC_ID HOMOGENIZATION_RGC_ID
private :: & private :: &
@ -982,8 +976,6 @@ subroutine material_parsePhase(fileUnit,myPart)
phase_plasticity(section) = PLASTICITY_NONE_ID phase_plasticity(section) = PLASTICITY_NONE_ID
case (PLASTICITY_ISOTROPIC_label) case (PLASTICITY_ISOTROPIC_label)
phase_plasticity(section) = PLASTICITY_ISOTROPIC_ID phase_plasticity(section) = PLASTICITY_ISOTROPIC_ID
case (PLASTICITY_J2_label)
phase_plasticity(section) = PLASTICITY_J2_ID
case (PLASTICITY_PHENOPOWERLAW_label) case (PLASTICITY_PHENOPOWERLAW_label)
phase_plasticity(section) = PLASTICITY_PHENOPOWERLAW_ID phase_plasticity(section) = PLASTICITY_PHENOPOWERLAW_ID
case (PLASTICITY_PHENOPLUS_label) case (PLASTICITY_PHENOPLUS_label)
@ -1280,7 +1272,7 @@ subroutine material_populateGrains
integer(pInt) :: t,e,i,g,j,m,c,r,homog,micro,sgn,hme, myDebug, & integer(pInt) :: t,e,i,g,j,m,c,r,homog,micro,sgn,hme, myDebug, &
phaseID,textureID,dGrains,myNgrains,myNorientations,myNconstituents, & phaseID,textureID,dGrains,myNgrains,myNorientations,myNconstituents, &
grain,constituentGrain,ipGrain,symExtension, ip grain,constituentGrain,ipGrain,symExtension, ip
real(pReal) :: extreme,rnd real(pReal) :: deviation,extreme,rnd
integer(pInt), dimension (:,:), allocatable :: Nelems ! counts number of elements in homog, micro array integer(pInt), dimension (:,:), allocatable :: Nelems ! counts number of elements in homog, micro array
type(p_intvec), dimension (:,:), allocatable :: elemsOfHomogMicro ! lists element number in homog, micro array type(p_intvec), dimension (:,:), allocatable :: elemsOfHomogMicro ! lists element number in homog, micro array
@ -1407,8 +1399,11 @@ subroutine material_populateGrains
extreme = 0.0_pReal extreme = 0.0_pReal
t = 0_pInt t = 0_pInt
do i = 1_pInt,myNconstituents ! find largest deviator do i = 1_pInt,myNconstituents ! find largest deviator
if (real(sgn,pReal)*log(NgrainsOfConstituent(i)/myNgrains/microstructure_fraction(i,micro)) > extreme) then deviation = real(sgn,pReal)*log( microstructure_fraction(i,micro) / &
extreme = real(sgn,pReal)*log(NgrainsOfConstituent(i)/myNgrains/microstructure_fraction(i,micro)) !-------------------------------- &
(real(NgrainsOfConstituent(i),pReal)/real(myNgrains,pReal) ) )
if (deviation > extreme) then
extreme = deviation
t = i t = i
endif endif
enddo enddo
@ -1600,14 +1595,4 @@ subroutine material_populateGrains
end subroutine material_populateGrains end subroutine material_populateGrains
#ifdef HDF
integer(pInt) pure function material_NconstituentsPhase(matID)
implicit none
integer(pInt), intent(in) :: matID
material_NconstituentsPhase = count(microstructure_phase == matID)
end function
#endif
end module material end module material

16
code/mesh.f90
View File

@ -4525,17 +4525,9 @@ subroutine mesh_write_cellGeom
VTK_geo, & VTK_geo, &
VTK_con, & VTK_con, &
VTK_end VTK_end
#ifdef HDF
use IO, only: &
HDF5_mappingCells
#endif
implicit none implicit none
integer(I4P), dimension(1:mesh_Ncells) :: celltype integer(I4P), dimension(1:mesh_Ncells) :: celltype
integer(I4P), dimension(mesh_Ncells*(1_pInt+FE_maxNcellnodesPerCell)) :: cellconnection integer(I4P), dimension(mesh_Ncells*(1_pInt+FE_maxNcellnodesPerCell)) :: cellconnection
#ifdef HDF
integer(pInt), dimension(mesh_Ncells*FE_maxNcellnodesPerCell) :: cellconnectionHDF5
integer(pInt) :: j2=0_pInt
#endif
integer(I4P):: error integer(I4P):: error
integer(I4P):: g, c, e, CellID, i, j integer(I4P):: g, c, e, CellID, i, j
@ -4550,16 +4542,8 @@ subroutine mesh_write_cellGeom
cellconnection(j+1_pInt:j+FE_NcellnodesPerCell(c)+1_pInt) & cellconnection(j+1_pInt:j+FE_NcellnodesPerCell(c)+1_pInt) &
= [FE_NcellnodesPerCell(c),mesh_cell(1:FE_NcellnodesPerCell(c),i,e)-1_pInt] ! number of cellnodes per cell & list of global cellnode IDs belnging to this cell (cellnode counting starts at 0) = [FE_NcellnodesPerCell(c),mesh_cell(1:FE_NcellnodesPerCell(c),i,e)-1_pInt] ! number of cellnodes per cell & list of global cellnode IDs belnging to this cell (cellnode counting starts at 0)
j = j + FE_NcellnodesPerCell(c) + 1_pInt j = j + FE_NcellnodesPerCell(c) + 1_pInt
#ifdef HDF
cellconnectionHDF5(j2+1_pInt:j2+FE_NcellnodesPerCell(c)) &
= mesh_cell(1:FE_NcellnodesPerCell(c),i,e)-1_pInt
j2=j2 + FE_ncellnodesPerCell(c)
#endif
enddo enddo
enddo enddo
#ifdef HDF
call HDF5_mappingCells(cellconnectionHDF5(1:j2))
#endif
error=VTK_ini(output_format = 'ASCII', & error=VTK_ini(output_format = 'ASCII', &
title=trim(getSolverJobName())//' cell mesh', & title=trim(getSolverJobName())//' cell mesh', &

737
code/plastic_disloUCLA.f90
View File

File diff suppressed because it is too large Load Diff

52
code/plastic_isotropic.f90
View File

@ -7,14 +7,10 @@
!! untextured polycrystal !! untextured polycrystal
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module plastic_isotropic module plastic_isotropic
#ifdef HDF
use hdf5, only: &
HID_T
#endif
use prec, only: & use prec, only: &
pReal,& pReal,&
pInt pInt, &
DAMASK_NaN
implicit none implicit none
private private
@ -40,16 +36,16 @@ module plastic_isotropic
integer(kind(undefined_ID)), allocatable, dimension(:) :: & integer(kind(undefined_ID)), allocatable, dimension(:) :: &
outputID outputID
real(pReal) :: & real(pReal) :: &
fTaylor, & fTaylor = DAMASK_NaN, &
tau0, & tau0 = DAMASK_NaN, &
gdot0, & gdot0 = DAMASK_NaN, &
n, & n = DAMASK_NaN, &
h0, & h0 = DAMASK_NaN, &
h0_slopeLnRate = 0.0_pReal, & h0_slopeLnRate = 0.0_pReal, &
tausat, & tausat = DAMASK_NaN, &
a, & a = DAMASK_NaN, &
aTolFlowstress, & aTolFlowstress = 1.0_pReal, &
aTolShear , & aTolShear = 1.0e-6_pReal, &
tausat_SinhFitA= 0.0_pReal, & tausat_SinhFitA= 0.0_pReal, &
tausat_SinhFitB= 0.0_pReal, & tausat_SinhFitB= 0.0_pReal, &
tausat_SinhFitC= 0.0_pReal, & tausat_SinhFitC= 0.0_pReal, &
@ -143,9 +139,10 @@ subroutine plastic_isotropic_init(fileUnit)
sizeDeltaState sizeDeltaState
character(len=65536) :: & character(len=65536) :: &
tag = '', & tag = '', &
outputtag = '', &
line = '', & line = '', &
extmsg = '' extmsg = ''
character(len=64) :: &
outputtag = ''
integer(pInt) :: NipcMyPhase integer(pInt) :: NipcMyPhase
mainProcess: if (worldrank == 0) then mainProcess: if (worldrank == 0) then
@ -385,8 +382,7 @@ subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el)
math_mul33xx33, & math_mul33xx33, &
math_transpose33 math_transpose33
use material, only: & use material, only: &
phaseAt, phasememberAt, & phasememberAt, &
plasticState, &
material_phase, & material_phase, &
phase_plasticityInstance phase_plasticityInstance
@ -416,7 +412,7 @@ subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el)
k, l, m, n k, l, m, n
of = phasememberAt(ipc,ip,el) ! phasememberAt should be tackled by material and be renamed to material_phasemember of = phasememberAt(ipc,ip,el) ! phasememberAt should be tackled by material and be renamed to material_phasemember
instance = phase_plasticityInstance(phaseAt(ipc,ip,el)) ! "phaseAt" equivalent to "material_phase" !! instance = phase_plasticityInstance(material_phase(ipc,ip,el))
Tstar_dev_33 = math_deviatoric33(math_Mandel6to33(Tstar_v)) ! deviatoric part of 2nd Piola-Kirchhoff stress Tstar_dev_33 = math_deviatoric33(math_Mandel6to33(Tstar_v)) ! deviatoric part of 2nd Piola-Kirchhoff stress
squarenorm_Tstar_dev = math_mul33xx33(Tstar_dev_33,Tstar_dev_33) squarenorm_Tstar_dev = math_mul33xx33(Tstar_dev_33,Tstar_dev_33)
@ -467,14 +463,14 @@ subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dTstar_3333,Tstar_v,ipc,ip,e
math_mul33xx33 math_mul33xx33
use material, only: & use material, only: &
phasememberAt, & phasememberAt, &
plasticState, &
material_phase, & material_phase, &
phase_plasticityInstance phase_plasticityInstance
implicit none implicit none
real(pReal), dimension(3,3), intent(out) :: & real(pReal), dimension(3,3), intent(out) :: &
Li !< plastic velocity gradient Li !< plastic velocity gradient
real(pReal), dimension(3,3,3,3), intent(out) :: &
dLi_dTstar_3333 !< derivative of Li with respect to Tstar as 4th order tensor
real(pReal), dimension(6), intent(in) :: & real(pReal), dimension(6), intent(in) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
integer(pInt), intent(in) :: & integer(pInt), intent(in) :: &
@ -484,8 +480,6 @@ subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dTstar_3333,Tstar_v,ipc,ip,e
real(pReal), dimension(3,3) :: & real(pReal), dimension(3,3) :: &
Tstar_sph_33 !< sphiatoric part of the 2nd Piola Kirchhoff stress tensor as 2nd order tensor Tstar_sph_33 !< sphiatoric part of the 2nd Piola Kirchhoff stress tensor as 2nd order tensor
real(pReal), dimension(3,3,3,3), intent(out) :: &
dLi_dTstar_3333 !< derivative of Li with respect to Tstar as 4th order tensor
real(pReal) :: & real(pReal) :: &
gamma_dot, & !< strainrate gamma_dot, & !< strainrate
norm_Tstar_sph, & !< euclidean norm of Tstar_sph norm_Tstar_sph, & !< euclidean norm of Tstar_sph
@ -501,11 +495,7 @@ subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dTstar_3333,Tstar_v,ipc,ip,e
squarenorm_Tstar_sph = math_mul33xx33(Tstar_sph_33,Tstar_sph_33) squarenorm_Tstar_sph = math_mul33xx33(Tstar_sph_33,Tstar_sph_33)
norm_Tstar_sph = sqrt(squarenorm_Tstar_sph) norm_Tstar_sph = sqrt(squarenorm_Tstar_sph)
if (param(instance)%dilatation) then if (param(instance)%dilatation .and. norm_Tstar_sph > 0.0_pReal) then ! Tstar == 0 or J2 plascitiy --> both Li and dLi_dTstar are zero
if (norm_Tstar_sph <= 0.0_pReal) then ! Tstar == 0 --> both Li and dLi_dTstar are zero
Li = 0.0_pReal
dLi_dTstar_3333 = 0.0_pReal
else
gamma_dot = param(instance)%gdot0 & gamma_dot = param(instance)%gdot0 &
* (sqrt(1.5_pReal) * norm_Tstar_sph / param(instance)%fTaylor / state(instance)%flowstress(of) ) & * (sqrt(1.5_pReal) * norm_Tstar_sph / param(instance)%fTaylor / state(instance)%flowstress(of) ) &
**param(instance)%n **param(instance)%n
@ -522,7 +512,9 @@ subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dTstar_3333,Tstar_v,ipc,ip,e
dLi_dTstar_3333 = gamma_dot / param(instance)%fTaylor * & dLi_dTstar_3333 = gamma_dot / param(instance)%fTaylor * &
dLi_dTstar_3333 / norm_Tstar_sph dLi_dTstar_3333 / norm_Tstar_sph
endif else
Li = 0.0_pReal
dLi_dTstar_3333 = 0.0_pReal
endif endif
end subroutine plastic_isotropic_LiAndItsTangent end subroutine plastic_isotropic_LiAndItsTangent
@ -536,7 +528,6 @@ subroutine plastic_isotropic_dotState(Tstar_v,ipc,ip,el)
math_mul6x6 math_mul6x6
use material, only: & use material, only: &
phasememberAt, & phasememberAt, &
plasticState, &
material_phase, & material_phase, &
phase_plasticityInstance phase_plasticityInstance
@ -615,7 +606,6 @@ function plastic_isotropic_postResults(Tstar_v,ipc,ip,el)
math_mul6x6 math_mul6x6
use material, only: & use material, only: &
material_phase, & material_phase, &
plasticState, &
phasememberAt, & phasememberAt, &
phase_plasticityInstance phase_plasticityInstance

577
code/plastic_j2.f90
View File

@ -1,577 +0,0 @@
!--------------------------------------------------------------------------------------------------
!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @brief material subroutine for isotropic (J2) plasticity
!> @details Isotropic (J2) Plasticity which resembles the phenopowerlaw plasticity without
!! resolving the stress on the slip systems. Will give the response of phenopowerlaw for an
!! untextured polycrystal
!--------------------------------------------------------------------------------------------------
module plastic_j2
#ifdef HDF
use hdf5, only: &
HID_T
#endif
use prec, only: &
pReal,&
pInt
implicit none
private
integer(pInt), dimension(:), allocatable, public, protected :: &
plastic_j2_sizePostResults !< cumulative size of post results
integer(pInt), dimension(:,:), allocatable, target, public :: &
plastic_j2_sizePostResult !< size of each post result output
character(len=64), dimension(:,:), allocatable, target, public :: &
plastic_j2_output !< name of each post result output
integer(pInt), dimension(:), allocatable, target, public :: &
plastic_j2_Noutput !< number of outputs per instance
real(pReal), dimension(:), allocatable, private :: &
plastic_j2_fTaylor, & !< Taylor factor
plastic_j2_tau0, & !< initial plastic stress
plastic_j2_gdot0, & !< reference velocity
plastic_j2_n, & !< Visco-plastic parameter
!--------------------------------------------------------------------------------------------------
! h0 as function of h0 = A + B log (gammadot)
plastic_j2_h0, &
plastic_j2_h0_slopeLnRate, &
plastic_j2_tausat, & !< final plastic stress
plastic_j2_a, &
plastic_j2_aTolResistance, &
plastic_j2_aTolShear, &
!--------------------------------------------------------------------------------------------------
! tausat += (asinh((gammadot / SinhFitA)**(1 / SinhFitD)))**(1 / SinhFitC) / (SinhFitB * (gammadot / gammadot0)**(1/n))
plastic_j2_tausat_SinhFitA, & !< fitting parameter for normalized strain rate vs. stress function
plastic_j2_tausat_SinhFitB, & !< fitting parameter for normalized strain rate vs. stress function
plastic_j2_tausat_SinhFitC, & !< fitting parameter for normalized strain rate vs. stress function
plastic_j2_tausat_SinhFitD !< fitting parameter for normalized strain rate vs. stress function
enum, bind(c)
enumerator :: undefined_ID, &
flowstress_ID, &
strainrate_ID
end enum
integer(kind(undefined_ID)), dimension(:,:), allocatable, private :: &
plastic_j2_outputID !< ID of each post result output
#ifdef HDF
type plastic_j2_tOutput
real(pReal), dimension(:), allocatable, private :: &
flowstress, &
strainrate
logical :: flowstressActive = .false., strainrateActive = .false. ! if we can write the output block wise, this is not needed anymore because we can do an if(allocated(xxx))
end type plastic_j2_tOutput
type(plastic_j2_tOutput), allocatable, dimension(:) :: plastic_j2_Output2
integer(HID_T), allocatable, dimension(:) :: outID
#endif
public :: &
plastic_j2_init, &
plastic_j2_LpAndItsTangent, &
plastic_j2_dotState, &
plastic_j2_postResults
contains
!--------------------------------------------------------------------------------------------------
!> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
subroutine plastic_j2_init(fileUnit)
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
#ifdef HDF
use hdf5
#endif
use debug, only: &
debug_level, &
debug_constitutive, &
debug_levelBasic
use numerics, only: &
analyticJaco, &
worldrank, &
numerics_integrator
use math, only: &
math_Mandel3333to66, &
math_Voigt66to3333
use IO, only: &
IO_read, &
IO_lc, &
IO_getTag, &
IO_isBlank, &
IO_stringPos, &
IO_stringValue, &
IO_floatValue, &
IO_error, &
IO_timeStamp, &
#ifdef HDF
tempResults, &
HDF5_addGroup, &
HDF5_addScalarDataset,&
#endif
IO_EOF
use material, only: &
phase_plasticity, &
phase_plasticityInstance, &
phase_Noutput, &
PLASTICITY_J2_label, &
PLASTICITY_J2_ID, &
material_phase, &
plasticState, &
MATERIAL_partPhase
use lattice
implicit none
integer(pInt), intent(in) :: fileUnit
integer(pInt), allocatable, dimension(:) :: chunkPos
integer(pInt) :: &
o, &
phase, &
maxNinstance, &
instance, &
mySize, &
sizeDotState, &
sizeState, &
sizeDeltaState
character(len=65536) :: &
tag = '', &
line = ''
integer(pInt) :: NofMyPhase
#ifdef HDF
character(len=5) :: &
str1
integer(HID_T) :: ID,ID2,ID4
#endif
mainProcess: if (worldrank == 0) then
write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_J2_label//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
endif mainProcess
maxNinstance = int(count(phase_plasticity == PLASTICITY_J2_ID),pInt)
if (maxNinstance == 0_pInt) return
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
#ifdef HDF
allocate(plastic_j2_Output2(maxNinstance))
allocate(outID(maxNinstance))
#endif
allocate(plastic_j2_sizePostResults(maxNinstance), source=0_pInt)
allocate(plastic_j2_sizePostResult(maxval(phase_Noutput), maxNinstance),source=0_pInt)
allocate(plastic_j2_output(maxval(phase_Noutput), maxNinstance))
plastic_j2_output = ''
allocate(plastic_j2_outputID(maxval(phase_Noutput),maxNinstance), source=undefined_ID)
allocate(plastic_j2_Noutput(maxNinstance), source=0_pInt)
allocate(plastic_j2_fTaylor(maxNinstance), source=0.0_pReal)
allocate(plastic_j2_tau0(maxNinstance), source=0.0_pReal)
allocate(plastic_j2_gdot0(maxNinstance), source=0.0_pReal)
allocate(plastic_j2_n(maxNinstance), source=0.0_pReal)
allocate(plastic_j2_h0(maxNinstance), source=0.0_pReal)
allocate(plastic_j2_h0_slopeLnRate(maxNinstance), source=0.0_pReal)
allocate(plastic_j2_tausat(maxNinstance), source=0.0_pReal)
allocate(plastic_j2_a(maxNinstance), source=0.0_pReal)
allocate(plastic_j2_aTolResistance(maxNinstance), source=0.0_pReal)
allocate(plastic_j2_aTolShear (maxNinstance), source=0.0_pReal)
allocate(plastic_j2_tausat_SinhFitA(maxNinstance), source=0.0_pReal)
allocate(plastic_j2_tausat_SinhFitB(maxNinstance), source=0.0_pReal)
allocate(plastic_j2_tausat_SinhFitC(maxNinstance), source=0.0_pReal)
allocate(plastic_j2_tausat_SinhFitD(maxNinstance), source=0.0_pReal)
rewind(fileUnit)
phase = 0_pInt
do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= material_partPhase) ! wind forward to <phase>
line = IO_read(fileUnit)
enddo
parsingFile: do while (trim(line) /= IO_EOF) ! read through sections of phase part
line = IO_read(fileUnit)
if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') then ! stop at next part
line = IO_read(fileUnit, .true.) ! reset IO_read
exit
endif
if (IO_getTag(line,'[',']') /= '') then ! next section
phase = phase + 1_pInt ! advance section counter
if (phase_plasticity(phase) == PLASTICITY_J2_ID) then
instance = phase_plasticityInstance(phase)
#ifdef HDF
outID(instance)=HDF5_addGroup(str1,tempResults)
#endif
endif
cycle ! skip to next line
endif
if (phase > 0_pInt ) then; if (phase_plasticity(phase) == PLASTICITY_J2_ID) then ! one of my phases. Do not short-circuit here (.and. between if-statements), it's not safe in Fortran
instance = phase_plasticityInstance(phase) ! which instance of my plasticity is present phase
chunkPos = IO_stringPos(line)
tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key
select case(tag)
case ('(output)')
select case(IO_lc(IO_stringValue(line,chunkPos,2_pInt)))
case ('flowstress')
plastic_j2_Noutput(instance) = plastic_j2_Noutput(instance) + 1_pInt
plastic_j2_outputID(plastic_j2_Noutput(instance),instance) = flowstress_ID
plastic_j2_output(plastic_j2_Noutput(instance),instance) = &
IO_lc(IO_stringValue(line,chunkPos,2_pInt))
#ifdef HDF
call HDF5_addScalarDataset(outID(instance),myConstituents,'flowstress','MPa')
allocate(plastic_j2_Output2(instance)%flowstress(myConstituents))
plastic_j2_Output2(instance)%flowstressActive = .true.
#endif
case ('strainrate')
plastic_j2_Noutput(instance) = plastic_j2_Noutput(instance) + 1_pInt
plastic_j2_outputID(plastic_j2_Noutput(instance),instance) = strainrate_ID
plastic_j2_output(plastic_j2_Noutput(instance),instance) = &
IO_lc(IO_stringValue(line,chunkPos,2_pInt))
#ifdef HDF
call HDF5_addScalarDataset(outID(instance),myConstituents,'strainrate','1/s')
allocate(plastic_j2_Output2(instance)%strainrate(myConstituents))
plastic_j2_Output2(instance)%strainrateActive = .true.
#endif
case default
end select
case ('tau0')
plastic_j2_tau0(instance) = IO_floatValue(line,chunkPos,2_pInt)
if (plastic_j2_tau0(instance) < 0.0_pReal) &
call IO_error(211_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_J2_label//')')
case ('gdot0')
plastic_j2_gdot0(instance) = IO_floatValue(line,chunkPos,2_pInt)
if (plastic_j2_gdot0(instance) <= 0.0_pReal) &
call IO_error(211_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_J2_label//')')
case ('n')
plastic_j2_n(instance) = IO_floatValue(line,chunkPos,2_pInt)
if (plastic_j2_n(instance) <= 0.0_pReal) &
call IO_error(211_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_J2_label//')')
case ('h0')
plastic_j2_h0(instance) = IO_floatValue(line,chunkPos,2_pInt)
case ('h0_slope','slopelnrate')
plastic_j2_h0_slopeLnRate(instance) = IO_floatValue(line,chunkPos,2_pInt)
case ('tausat')
plastic_j2_tausat(instance) = IO_floatValue(line,chunkPos,2_pInt)
if (plastic_j2_tausat(instance) <= 0.0_pReal) &
call IO_error(211_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_J2_label//')')
case ('tausat_sinhfita')
plastic_j2_tausat_SinhFitA(instance) = IO_floatValue(line,chunkPos,2_pInt)
case ('tausat_sinhfitb')
plastic_j2_tausat_SinhFitB(instance) = IO_floatValue(line,chunkPos,2_pInt)
case ('tausat_sinhfitc')
plastic_j2_tausat_SinhFitC(instance) = IO_floatValue(line,chunkPos,2_pInt)
case ('tausat_sinhfitd')
plastic_j2_tausat_SinhFitD(instance) = IO_floatValue(line,chunkPos,2_pInt)
case ('a', 'w0')
plastic_j2_a(instance) = IO_floatValue(line,chunkPos,2_pInt)
if (plastic_j2_a(instance) <= 0.0_pReal) &
call IO_error(211_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_J2_label//')')
case ('taylorfactor')
plastic_j2_fTaylor(instance) = IO_floatValue(line,chunkPos,2_pInt)
if (plastic_j2_fTaylor(instance) <= 0.0_pReal) &
call IO_error(211_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_J2_label//')')
case ('atol_resistance')
plastic_j2_aTolResistance(instance) = IO_floatValue(line,chunkPos,2_pInt)
if (plastic_j2_aTolResistance(instance) <= 0.0_pReal) &
call IO_error(211_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_J2_label//')')
case ('atol_shear')
plastic_j2_aTolShear(instance) = IO_floatValue(line,chunkPos,2_pInt)
case default
end select
endif; endif
enddo parsingFile
initializeInstances: do phase = 1_pInt, size(phase_plasticity)
myPhase: if (phase_plasticity(phase) == PLASTICITY_j2_ID) then
NofMyPhase=count(material_phase==phase)
instance = phase_plasticityInstance(phase)
!--------------------------------------------------------------------------------------------------
! sanity checks
if (plastic_j2_aTolShear(instance) <= 0.0_pReal) &
plastic_j2_aTolShear(instance) = 1.0e-6_pReal ! default absolute tolerance 1e-6
!--------------------------------------------------------------------------------------------------
! Determine size of postResults array
outputsLoop: do o = 1_pInt,plastic_j2_Noutput(instance)
select case(plastic_j2_outputID(o,instance))
case(flowstress_ID,strainrate_ID)
mySize = 1_pInt
case default
end select
outputFound: if (mySize > 0_pInt) then
plastic_j2_sizePostResult(o,instance) = mySize
plastic_j2_sizePostResults(instance) = &
plastic_j2_sizePostResults(instance) + mySize
endif outputFound
enddo outputsLoop
!--------------------------------------------------------------------------------------------------
! allocate state arrays
sizeState = 2_pInt
sizeDotState = sizeState
sizeDeltaState = 0_pInt
plasticState(phase)%sizeState = sizeState
plasticState(phase)%sizeDotState = sizeDotState
plasticState(phase)%sizeDeltaState = sizeDeltaState
plasticState(phase)%sizePostResults = plastic_j2_sizePostResults(instance)
plasticState(phase)%nSlip = 1
plasticState(phase)%nTwin = 0
plasticState(phase)%nTrans= 0
allocate(plasticState(phase)%aTolState ( sizeState))
plasticState(phase)%aTolState(1) = plastic_j2_aTolResistance(instance)
plasticState(phase)%aTolState(2) = plastic_j2_aTolShear(instance)
allocate(plasticState(phase)%state0 ( sizeState,NofMyPhase))
plasticState(phase)%state0(1,1:NofMyPhase) = plastic_j2_tau0(instance)
plasticState(phase)%state0(2,1:NofMyPhase) = 0.0_pReal
allocate(plasticState(phase)%partionedState0 ( sizeState,NofMyPhase),source=0.0_pReal)
allocate(plasticState(phase)%subState0 ( sizeState,NofMyPhase),source=0.0_pReal)
allocate(plasticState(phase)%state ( sizeState,NofMyPhase),source=0.0_pReal)
allocate(plasticState(phase)%dotState (sizeDotState,NofMyPhase),source=0.0_pReal)
allocate(plasticState(phase)%deltaState (sizeDeltaState,NofMyPhase),source=0.0_pReal)
if (.not. analyticJaco) then
allocate(plasticState(phase)%state_backup ( sizeState,NofMyPhase),source=0.0_pReal)
allocate(plasticState(phase)%dotState_backup (sizeDotState,NofMyPhase),source=0.0_pReal)
endif
if (any(numerics_integrator == 1_pInt)) then
allocate(plasticState(phase)%previousDotState (sizeDotState,NofMyPhase),source=0.0_pReal)
allocate(plasticState(phase)%previousDotState2(sizeDotState,NofMyPhase),source=0.0_pReal)
endif
if (any(numerics_integrator == 4_pInt)) &
allocate(plasticState(phase)%RK4dotState (sizeDotState,NofMyPhase),source=0.0_pReal)
if (any(numerics_integrator == 5_pInt)) &
allocate(plasticState(phase)%RKCK45dotState (6,sizeDotState,NofMyPhase),source=0.0_pReal)
plasticState(phase)%slipRate => plasticState(phase)%dotState(2:2,1:NofMyPhase)
plasticState(phase)%accumulatedSlip => plasticState(phase)%state (2:2,1:NofMyPhase)
endif myPhase
enddo initializeInstances
end subroutine plastic_j2_init
!--------------------------------------------------------------------------------------------------
!> @brief calculates plastic velocity gradient and its tangent
!--------------------------------------------------------------------------------------------------
subroutine plastic_j2_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el)
use math, only: &
math_mul6x6, &
math_Mandel6to33, &
math_Plain3333to99, &
math_deviatoric33, &
math_mul33xx33
use material, only: &
phaseAt, phasememberAt, &
plasticState, &
material_phase, &
phase_plasticityInstance
implicit none
real(pReal), dimension(3,3), intent(out) :: &
Lp !< plastic velocity gradient
real(pReal), dimension(9,9), intent(out) :: &
dLp_dTstar99 !< derivative of Lp with respect to 2nd Piola Kirchhoff stress
real(pReal), dimension(6), intent(in) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
real(pReal), dimension(3,3) :: &
Tstar_dev_33 !< deviatoric part of the 2nd Piola Kirchhoff stress tensor as 2nd order tensor
real(pReal), dimension(3,3,3,3) :: &
dLp_dTstar_3333 !< derivative of Lp with respect to Tstar as 4th order tensor
real(pReal) :: &
gamma_dot, & !< strainrate
norm_Tstar_dev, & !< euclidean norm of Tstar_dev
squarenorm_Tstar_dev !< square of the euclidean norm of Tstar_dev
integer(pInt) :: &
instance, &
k, l, m, n
instance = phase_plasticityInstance(material_phase(ipc,ip,el))
Tstar_dev_33 = math_deviatoric33(math_Mandel6to33(Tstar_v)) ! deviatoric part of 2nd Piola-Kirchhoff stress
squarenorm_Tstar_dev = math_mul33xx33(Tstar_dev_33,Tstar_dev_33)
norm_Tstar_dev = sqrt(squarenorm_Tstar_dev)
if (norm_Tstar_dev <= 0.0_pReal) then ! Tstar == 0 --> both Lp and dLp_dTstar are zero
Lp = 0.0_pReal
dLp_dTstar99 = 0.0_pReal
else
gamma_dot = plastic_j2_gdot0(instance) &
* (sqrt(1.5_pReal) * norm_Tstar_dev / (plastic_j2_fTaylor(instance) * &
plasticState(phaseAt(ipc,ip,el))%state(1,phasememberAt(ipc,ip,el)))) &
**plastic_j2_n(instance)
Lp = Tstar_dev_33/norm_Tstar_dev * gamma_dot/plastic_j2_fTaylor(instance)
!--------------------------------------------------------------------------------------------------
! Calculation of the tangent of Lp
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
dLp_dTstar_3333(k,l,m,n) = (plastic_j2_n(instance)-1.0_pReal) * &
Tstar_dev_33(k,l)*Tstar_dev_33(m,n) / squarenorm_Tstar_dev
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt) &
dLp_dTstar_3333(k,l,k,l) = dLp_dTstar_3333(k,l,k,l) + 1.0_pReal
forall (k=1_pInt:3_pInt,m=1_pInt:3_pInt) &
dLp_dTstar_3333(k,k,m,m) = dLp_dTstar_3333(k,k,m,m) - 1.0_pReal/3.0_pReal
dLp_dTstar99 = math_Plain3333to99(gamma_dot / plastic_j2_fTaylor(instance) * &
dLp_dTstar_3333 / norm_Tstar_dev)
end if
end subroutine plastic_j2_LpAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief calculates the rate of change of microstructure
!--------------------------------------------------------------------------------------------------
subroutine plastic_j2_dotState(Tstar_v,ipc,ip,el)
use math, only: &
math_mul6x6
use material, only: &
phaseAt, phasememberAt, &
plasticState, &
material_phase, &
phase_plasticityInstance
implicit none
real(pReal), dimension(6), intent(in):: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
real(pReal), dimension(6) :: &
Tstar_dev_v !< deviatoric part of the 2nd Piola Kirchhoff stress tensor in Mandel notation
real(pReal) :: &
gamma_dot, & !< strainrate
hardening, & !< hardening coefficient
saturation, & !< saturation resistance
norm_Tstar_dev !< euclidean norm of Tstar_dev
integer(pInt) :: &
instance, & !< instance of my instance (unique number of my constitutive model)
of, & !< shortcut notation for offset position in state array
ph !< shortcut notation for phase ID (unique number of all phases, regardless of constitutive model)
of = phasememberAt(ipc,ip,el)
ph = phaseAt(ipc,ip,el)
instance = phase_plasticityInstance(material_phase(ipc,ip,el))
!--------------------------------------------------------------------------------------------------
! norm of deviatoric part of 2nd Piola-Kirchhoff stress
Tstar_dev_v(1:3) = Tstar_v(1:3) - sum(Tstar_v(1:3))/3.0_pReal
Tstar_dev_v(4:6) = Tstar_v(4:6)
norm_Tstar_dev = sqrt(math_mul6x6(Tstar_dev_v,Tstar_dev_v))
!--------------------------------------------------------------------------------------------------
! strain rate
gamma_dot = plastic_j2_gdot0(instance) * ( sqrt(1.5_pReal) * norm_Tstar_dev &
/ &!-----------------------------------------------------------------------------------
(plastic_j2_fTaylor(instance)*plasticState(ph)%state(1,of)) )**plastic_j2_n(instance)
!--------------------------------------------------------------------------------------------------
! hardening coefficient
if (abs(gamma_dot) > 1e-12_pReal) then
if (abs(plastic_j2_tausat_SinhFitA(instance)) <= tiny(0.0_pReal)) then
saturation = plastic_j2_tausat(instance)
else
saturation = ( plastic_j2_tausat(instance) &
+ ( log( ( gamma_dot / plastic_j2_tausat_SinhFitA(instance)&
)**(1.0_pReal / plastic_j2_tausat_SinhFitD(instance))&
+ sqrt( ( gamma_dot / plastic_j2_tausat_SinhFitA(instance) &
)**(2.0_pReal / plastic_j2_tausat_SinhFitD(instance)) &
+ 1.0_pReal ) &
) & ! asinh(K) = ln(K + sqrt(K^2 +1))
)**(1.0_pReal / plastic_j2_tausat_SinhFitC(instance)) &
/ ( plastic_j2_tausat_SinhFitB(instance) &
* (gamma_dot / plastic_j2_gdot0(instance))**(1.0_pReal / plastic_j2_n(instance)) &
) &
)
endif
hardening = ( plastic_j2_h0(instance) + plastic_j2_h0_slopeLnRate(instance) * log(gamma_dot) ) &
* abs( 1.0_pReal - plasticState(ph)%state(1,of)/saturation )**plastic_j2_a(instance) &
* sign(1.0_pReal, 1.0_pReal - plasticState(ph)%state(1,of)/saturation)
else
hardening = 0.0_pReal
endif
plasticState(ph)%dotState(1,of) = hardening * gamma_dot
plasticState(ph)%dotState(2,of) = gamma_dot
end subroutine plastic_j2_dotState
!--------------------------------------------------------------------------------------------------
!> @brief return array of constitutive results
!--------------------------------------------------------------------------------------------------
function plastic_j2_postResults(Tstar_v,ipc,ip,el)
use math, only: &
math_mul6x6
use material, only: &
material_phase, &
plasticState, &
phaseAt, phasememberAt, &
phase_plasticityInstance
implicit none
real(pReal), dimension(6), intent(in) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
real(pReal), dimension(plastic_j2_sizePostResults(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
plastic_j2_postResults
real(pReal), dimension(6) :: &
Tstar_dev_v ! deviatoric part of the 2nd Piola Kirchhoff stress tensor in Mandel notation
real(pReal) :: &
norm_Tstar_dev ! euclidean norm of Tstar_dev
integer(pInt) :: &
instance, & !< instance of my instance (unique number of my constitutive model)
of, & !< shortcut notation for offset position in state array
ph, & !< shortcut notation for phase ID (unique number of all phases, regardless of constitutive model)
c, &
o
of = phasememberAt(ipc,ip,el)
ph = phaseAt(ipc,ip,el)
instance = phase_plasticityInstance(material_phase(ipc,ip,el))
!--------------------------------------------------------------------------------------------------
! calculate deviatoric part of 2nd Piola-Kirchhoff stress and its norm
Tstar_dev_v(1:3) = Tstar_v(1:3) - sum(Tstar_v(1:3))/3.0_pReal
Tstar_dev_v(4:6) = Tstar_v(4:6)
norm_Tstar_dev = sqrt(math_mul6x6(Tstar_dev_v,Tstar_dev_v))
c = 0_pInt
plastic_j2_postResults = 0.0_pReal
outputsLoop: do o = 1_pInt,plastic_j2_Noutput(instance)
select case(plastic_j2_outputID(o,instance))
case (flowstress_ID)
plastic_j2_postResults(c+1_pInt) = plasticState(ph)%state(1,of)
c = c + 1_pInt
case (strainrate_ID)
plastic_j2_postResults(c+1_pInt) = &
plastic_j2_gdot0(instance) * ( sqrt(1.5_pReal) * norm_Tstar_dev &
/ &!----------------------------------------------------------------------------------
(plastic_j2_fTaylor(instance) * plasticState(ph)%state(1,of)) ) ** plastic_j2_n(instance)
c = c + 1_pInt
end select
enddo outputsLoop
end function plastic_j2_postResults
end module plastic_j2

4
code/prec.f90
View File

@ -113,7 +113,9 @@ module prec
public :: & public :: &
prec_init, & prec_init, &
prec_isNaN prec_isNaN, &
dEq, &
dNeq
contains contains

37
config/Phase_DisloKMC_Tungsten.config
View File

@ -1,37 +0,0 @@
### $Id$ ###
[Tungsten]
elasticity hooke
plasticity dislokmc
### Material parameters ###
lattice_structure bcc
C11 523.0e9 # From Marinica et al. Journal of Physics: Condensed Matter(2013)
C12 202.0e9
C44 161.0e9
grainsize 2.0e-5 # Average grain size [m] 2.0e-5
SolidSolutionStrength 0.0 # Strength due to elements in solid solution
### Dislocation glide parameters ###
#per family
Nslip 12 0
slipburgers 2.72e-10 # Burgers vector of slip system [m]
rhoedge0 1.0e12 # Initial edge dislocation density [m/m**3]
rhoedgedip0 1.0 # Initial edged dipole dislocation density [m/m**3]
Qedge 2.725e-19 # Activation energy for dislocation glide [J]
v0 3560.3 # Initial glide velocity [m/s](kmC)
p_slip 0.16 # p-exponent in glide velocity
q_slip 1.00 # q-exponent in glide velocity
u_slip 2.47 # u-exponent of stress pre-factor (kmC)
s_slip 0.97 # self hardening in glide velocity (kmC)
tau_peierls 2.03e9 # peierls stress [Pa]
#hardening
dipoleformationfactor 0 # to have hardening due to dipole formation off
CLambdaSlip 10.0 # Adj. parameter controlling dislocation mean free path
D0 4.0e-5 # Vacancy diffusion prefactor [m**2/s]
Qsd 4.5e-19 # Activation energy for climb [J]
Catomicvolume 1.0 # Adj. parameter controlling the atomic volume [in b]
Cedgedipmindistance 1.0 # Adj. parameter controlling the minimum dipole distance [in b]
interaction_slipslip 0.2 0.11 0.19 0.15 0.11 0.17

2
config/Phase_J2_AluminumIsotropic.config → config/Phase_Isotropic_AluminumIsotropic.config
View File

@ -3,7 +3,7 @@
# Kuo, J. C., Mikrostrukturmechanik von Bikristallen mit Kippkorngrenzen. Shaker-Verlag 2004. http://edoc.mpg.de/204079 # Kuo, J. C., Mikrostrukturmechanik von Bikristallen mit Kippkorngrenzen. Shaker-Verlag 2004. http://edoc.mpg.de/204079
elasticity hooke elasticity hooke
plasticity j2 plasticity isotropic
(output) flowstress (output) flowstress
(output) strainrate (output) strainrate

4
installation/symlink_Processing.py
View File

@ -21,13 +21,13 @@ if not os.path.isdir(binDir):
#define ToDo list #define ToDo list
processing_subDirs = ['pre','post','misc',] processing_subDirs = ['pre','post','misc',]
processing_extensions = ['.py',] processing_extensions = ['.py','.sh',]
for subDir in processing_subDirs: for subDir in processing_subDirs:
theDir = os.path.abspath(os.path.join(baseDir,subDir)) theDir = os.path.abspath(os.path.join(baseDir,subDir))
for theFile in os.listdir(theDir): for theFile in os.listdir(theDir):
if os.path.splitext(theFile)[1] in processing_extensions: # omit anything not fitting our script extensions (skip .py.bak, .py~, and the like) if os.path.splitext(theFile)[1] in processing_extensions: # only consider files with proper extensions
src = os.path.abspath(os.path.join(theDir,theFile)) src = os.path.abspath(os.path.join(theDir,theFile))
sym_link = os.path.abspath(os.path.join(binDir,os.path.splitext(theFile)[0])) sym_link = os.path.abspath(os.path.join(binDir,os.path.splitext(theFile)[0]))

43
lib/damask/asciitable.py
View File

@ -40,12 +40,12 @@ class ASCIItable():
self.__IO__['in'] = name self.__IO__['in'] = name
try: try:
self.__IO__['out'] = (open(outname,'w') if (not os.path.isfile(outname) \ self.__IO__['out'] = (open(outname,'w') if (not os.path.isfile(outname) or
or os.access( outname, os.W_OK) \ os.access( outname, os.W_OK)
) \ ) and
and (not self.__IO__['inPlace'] \ (not self.__IO__['inPlace'] or
or not os.path.isfile(name) \ not os.path.isfile(name) or
or os.access( name, os.W_OK) \ os.access( name, os.W_OK)
) else None) if outname else sys.stdout ) else None) if outname else sys.stdout
except TypeError: except TypeError:
self.__IO__['out'] = outname self.__IO__['out'] = outname
@ -66,6 +66,14 @@ class ASCIItable():
except: except:
return 0.0 return 0.0
# ------------------------------------------------------------------
def _removeCRLF(self,
string):
try:
return string.replace('\n','').replace('\r','')
except:
return string
# ------------------------------------------------------------------ # ------------------------------------------------------------------
def close(self, def close(self,
dismiss = False): dismiss = False):
@ -243,9 +251,9 @@ class ASCIItable():
try: try:
for item in what: self.labels_append(item) for item in what: self.labels_append(item)
except: except:
self.labels += [str(what)] self.labels += [self._removeCRLF(str(what))]
else: else:
self.labels += [what] self.labels += [self._removeCRLF(what)]
self.__IO__['labeled'] = True # switch on processing (in particular writing) of labels self.__IO__['labeled'] = True # switch on processing (in particular writing) of labels
if reset: self.__IO__['labels'] = list(self.labels) # subsequent data_read uses current labels as data size if reset: self.__IO__['labels'] = list(self.labels) # subsequent data_read uses current labels as data size
@ -272,7 +280,7 @@ class ASCIItable():
for label in labels: for label in labels:
if label is not None: if label is not None:
try: try:
idx.append(int(label)) # column given as integer number? idx.append(int(label)-1) # column given as integer number?
except ValueError: except ValueError:
try: try:
idx.append(self.labels.index(label)) # locate string in label list idx.append(self.labels.index(label)) # locate string in label list
@ -283,7 +291,7 @@ class ASCIItable():
idx.append(-1) # not found... idx.append(-1) # not found...
else: else:
try: try:
idx = int(labels) idx = int(labels)-1 # offset for python array indexing
except ValueError: except ValueError:
try: try:
idx = self.labels.index(labels) idx = self.labels.index(labels)
@ -293,7 +301,7 @@ class ASCIItable():
except ValueError: except ValueError:
idx = None if labels is None else -1 idx = None if labels is None else -1
return np.array(idx) if isinstance(idx,list) else idx return np.array(idx) if isinstance(idx,Iterable) else idx
# ------------------------------------------------------------------ # ------------------------------------------------------------------
def label_dimension(self, def label_dimension(self,
@ -312,7 +320,7 @@ class ASCIItable():
if label is not None: if label is not None:
myDim = -1 myDim = -1
try: # column given as number? try: # column given as number?
idx = int(label) idx = int(label)-1
myDim = 1 # if found has at least dimension 1 myDim = 1 # if found has at least dimension 1
if self.labels[idx].startswith('1_'): # column has multidim indicator? if self.labels[idx].startswith('1_'): # column has multidim indicator?
while idx+myDim < len(self.labels) and self.labels[idx+myDim].startswith("%i_"%(myDim+1)): while idx+myDim < len(self.labels) and self.labels[idx+myDim].startswith("%i_"%(myDim+1)):
@ -331,7 +339,7 @@ class ASCIItable():
dim = -1 # assume invalid label dim = -1 # assume invalid label
idx = -1 idx = -1
try: # column given as number? try: # column given as number?
idx = int(labels) idx = int(labels)-1
dim = 1 # if found has at least dimension 1 dim = 1 # if found has at least dimension 1
if self.labels[idx].startswith('1_'): # column has multidim indicator? if self.labels[idx].startswith('1_'): # column has multidim indicator?
while idx+dim < len(self.labels) and self.labels[idx+dim].startswith("%i_"%(dim+1)): while idx+dim < len(self.labels) and self.labels[idx+dim].startswith("%i_"%(dim+1)):
@ -345,7 +353,7 @@ class ASCIItable():
while idx+dim < len(self.labels) and self.labels[idx+dim].startswith("%i_"%(dim+1)): while idx+dim < len(self.labels) and self.labels[idx+dim].startswith("%i_"%(dim+1)):
dim += 1 # keep adding while going through object dim += 1 # keep adding while going through object
return np.array(dim) if isinstance(dim,list) else dim return np.array(dim) if isinstance(dim,Iterable) else dim
# ------------------------------------------------------------------ # ------------------------------------------------------------------
def label_indexrange(self, def label_indexrange(self,
@ -361,7 +369,8 @@ class ASCIItable():
start = self.label_index(labels) start = self.label_index(labels)
dim = self.label_dimension(labels) dim = self.label_dimension(labels)
return map(lambda a,b: xrange(a,a+b), zip(start,dim)) if isinstance(labels, Iterable) and not isinstance(labels, str) \ return np.hstack(map(lambda c: xrange(c[0],c[0]+c[1]), zip(start,dim))) \
if isinstance(labels, Iterable) and not isinstance(labels, str) \
else xrange(start,start+dim) else xrange(start,start+dim)
# ------------------------------------------------------------------ # ------------------------------------------------------------------
@ -372,9 +381,9 @@ class ASCIItable():
try: try:
for item in what: self.info_append(item) for item in what: self.info_append(item)
except: except:
self.info += [str(what)] self.info += [self._removeCRLF(str(what))]
else: else:
self.info += [what] self.info += [self._removeCRLF(what)]
# ------------------------------------------------------------------ # ------------------------------------------------------------------
def info_clear(self): def info_clear(self):

14
lib/damask/util.py
View File

@ -53,6 +53,20 @@ def report(who,what):
"""reports script and file name""" """reports script and file name"""
croak( (emph(who) if who else '') + (': '+what if what else '') ) croak( (emph(who) if who else '') + (': '+what if what else '') )
# -----------------------------
def report_geom(info,
what = ['grid','size','origin','homogenization','microstructures']):
"""reports (selected) geometry information"""
output = {
'grid' : 'grid a b c: {}'.format(' x '.join(map(str,info['grid' ]))),
'size' : 'size x y z: {}'.format(' x '.join(map(str,info['size' ]))),
'origin' : 'origin x y z: {}'.format(' : '.join(map(str,info['origin']))),
'homogenization' : 'homogenization: {}'.format(info['homogenization']),
'microstructures' : 'microstructures: {}'.format(info['microstructures']),
}
for item in what: croak(output[item.lower()])
# ----------------------------- # -----------------------------
def emph(what): def emph(what):
"""emphasizes string on screen""" """emphasizes string on screen"""

BIN
misc/DAMASK QR-Code.png
View File

Binary file not shown.
Side by Side Swipe Overlay

Before

Width:  |  Height:  |  Size: 80 KiB

After

Width:  |  Height:  |  Size: 34 KiB

61
processing/post/addCauchyX.py
View File

@ -1,61 +0,0 @@
#!/usr/bin/env python
# -*- coding: UTF-8 no BOM -*-
import os,string,h5py
import numpy as np
from optparse import OptionParser
import damask
# --------------------------------------------------------------------
# MAIN
# --------------------------------------------------------------------
parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """
Add column(s) containing Cauchy stress based on given column(s) of
deformation gradient and first Piola--Kirchhoff stress.
""" + string.replace('$Id$','\n','\\n')
)
parser.add_option('-f','--defgrad', dest='defgrad', \
help='heading of columns containing deformation gradient [%default]')
parser.add_option('-p','--stress', dest='stress', \
help='heading of columns containing first Piola--Kirchhoff stress [%default]')
parser.add_option('-o','--output', dest='output', \
help='group containing requested data [%default]')
parser.set_defaults(defgrad = 'f')
parser.set_defaults(stress = 'p')
parser.set_defaults(output = 'crystallite')
(options,filenames) = parser.parse_args()
if options.defgrad is None or options.stress is None or options.output is None:
parser.error('missing data column...')
# ------------------------------------------ setup file handles ---------------------------------------
files = []
for name in filenames:
if os.path.exists(name):
files.append({'name':name, 'file':h5py.File(name,"a")})
# ------------------------------------------ loop over input files ------------------------------------
for myFile in files:
print(myFile['name'])
# ------------------------------------------ loop over increments -------------------------------------
for inc in myFile['file']['increments'].keys():
print("Current Increment: "+inc)
for instance in myFile['file']['increments/'+inc+'/'+options.output].keys():
dsets = myFile['file']['increments/'+inc+'/'+options.output+'/'+instance].keys()
if (options.defgrad in dsets and options.stress in dsets):
defgrad = myFile['file']['increments/'+inc+'/'+options.output+'/'+instance+'/'+options.defgrad]
stress = myFile['file']['increments/'+inc+'/'+options.output+'/'+instance+'/'+options.stress]
cauchy=np.zeros(np.shape(stress),'f')
for p in range(stress.shape[0]):
cauchy[p,...] = 1.0/np.linalg.det(defgrad[p,...])*np.dot(stress[p,...],defgrad[p,...].T) # [Cauchy] = (1/det(F)) * [P].[F_transpose]
cauchyFile = myFile['file']['increments/'+inc+'/'+options.output+'/'+instance].create_dataset('cauchy', data=cauchy)
cauchyFile.attrs['units'] = 'Pa'

2
processing/post/addCompatibilityMismatch.py
View File

@ -36,7 +36,7 @@ parser.add_option('--no-volume','-v',
dest = 'volume', dest = 'volume',
action = 'store_false', action = 'store_false',
help = 'omit volume mismatch') help = 'omit volume mismatch')
parser.set_defaults(coords = 'ipinitialcoord', parser.set_defaults(coords = 'pos',
defgrad = 'f', defgrad = 'f',
shape = True, shape = True,
volume = True, volume = True,

8
processing/post/addCurl.py
View File

@ -72,17 +72,17 @@ Deals with both vector- and tensor-valued fields.
parser.add_option('-c','--coordinates', parser.add_option('-c','--coordinates',
dest = 'coords', dest = 'coords',
type = 'string', metavar = 'string', type = 'string', metavar = 'string',
help = 'column heading for coordinates [%default]') help = 'column label of coordinates [%default]')
parser.add_option('-v','--vector', parser.add_option('-v','--vector',
dest = 'vector', dest = 'vector',
action = 'extend', metavar = '<string LIST>', action = 'extend', metavar = '<string LIST>',
help = 'heading of columns containing vector field values') help = 'column label(s) of vector field values')
parser.add_option('-t','--tensor', parser.add_option('-t','--tensor',
dest = 'tensor', dest = 'tensor',
action = 'extend', metavar = '<string LIST>', action = 'extend', metavar = '<string LIST>',
help = 'heading of columns containing tensor field values') help = 'column label(s) of tensor field values')
parser.set_defaults(coords = 'ipinitialcoord', parser.set_defaults(coords = 'pos',
) )
(options,filenames) = parser.parse_args() (options,filenames) = parser.parse_args()

164
processing/post/addDeformedConfiguration.py
View File

@ -1,164 +0,0 @@
#!/usr/bin/env python
# -*- coding: UTF-8 no BOM -*-
import os,sys,math
import numpy as np
from optparse import OptionParser
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
#--------------------------------------------------------------------------------------------------
def deformedCoordsFFT(F,undeformed=False):
wgt = 1.0/grid.prod()
integrator = np.array([0.+1.j,0.+1.j,0.+1.j],'c16') * size/ 2.0 / math.pi
step = size/grid
F_fourier = np.fft.rfftn(F,axes=(0,1,2))
coords_fourier = np.zeros(F_fourier.shape[0:4],'c16')
if undeformed:
Favg=np.eye(3)
else:
Favg=np.real(F_fourier[0,0,0,:,:])*wgt
#--------------------------------------------------------------------------------------------------
# integration in Fourier space
k_s = np.zeros([3],'i')
for i in xrange(grid[2]):
k_s[2] = i
if(i > grid[2]//2 ): k_s[2] = k_s[2] - grid[2]
for j in xrange(grid[1]):
k_s[1] = j
if(j > grid[1]//2 ): k_s[1] = k_s[1] - grid[1]
for k in xrange(grid[0]//2+1):
k_s[0] = k
for m in xrange(3):
coords_fourier[i,j,k,m] = sum(F_fourier[i,j,k,m,0:3]*k_s*integrator)
if (any(k_s != 0)):
coords_fourier[i,j,k,0:3] /= -sum(k_s*k_s)
#--------------------------------------------------------------------------------------------------
# add average to scaled fluctuation and put (0,0,0) on (0,0,0)
coords = np.fft.irfftn(coords_fourier,F.shape[0:3],axes=(0,1,2))
offset_coords = np.dot(F[0,0,0,:,:],step/2.0) - scaling*coords[0,0,0,0:3]
for z in xrange(grid[2]):
for y in xrange(grid[1]):
for x in xrange(grid[0]):
coords[z,y,x,0:3] = scaling*coords[z,y,x,0:3] \
+ offset_coords \
+ np.dot(Favg,step*np.array([x,y,z]))
return coords
# --------------------------------------------------------------------
# MAIN
# --------------------------------------------------------------------
parser = OptionParser(option_class=damask.extendableOption, usage='%prog options file[s]', description = """
Add deformed configuration of given initial coordinates.
Operates on periodic three-dimensional x,y,z-ordered data sets.
""", version = scriptID)
parser.add_option('-f', '--defgrad',dest='defgrad', metavar = 'string',
help='heading of deformation gradient columns [%default]')
parser.add_option('--reference', dest='undeformed', action='store_true',
help='map results to reference (undeformed) average configuration [%default]')
parser.add_option('--scaling', dest='scaling', action='extend', metavar = '<float LIST>',
help='scaling of fluctuation')
parser.add_option('-u', '--unitlength', dest='unitlength', type='float', metavar = 'float',
help='set unit length for 2D model [%default]')
parser.add_option('--coordinates', dest='coords', metavar='string',
help='column heading for coordinates [%default]')
parser.set_defaults(defgrad = 'f')
parser.set_defaults(coords = 'ipinitialcoord')
parser.set_defaults(scaling = [])
parser.set_defaults(undeformed = False)
parser.set_defaults(unitlength = 0.0)
(options,filenames) = parser.parse_args()
options.scaling += [1.0 for i in xrange(max(0,3-len(options.scaling)))]
scaling = map(float, options.scaling)
# --- loop over input files -------------------------------------------------------------------------
if filenames == []: filenames = [None]
for name in filenames:
try:
table = damask.ASCIItable(name = name,
buffered = False)
except: continue
damask.util.report(scriptName,name)
# ------------------------------------------ read header ------------------------------------------
table.head_read()
# ------------------------------------------ sanity checks ----------------------------------------
errors = []
remarks = []
if table.label_dimension(options.coords) != 3: errors.append('coordinates {} are not a vector.'.format(options.coords))
else: colCoord = table.label_index(options.coords)
if table.label_dimension(options.defgrad) != 9: errors.append('deformation gradient {} is not a tensor.'.format(options.defgrad))
else: colF = table.label_index(options.defgrad)
if remarks != []: damask.util.croak(remarks)
if errors != []:
damask.util.croak(errors)
table.close(dismiss = True)
continue
# --------------- figure out size and grid ---------------------------------------------------------
table.data_readArray()
coords = [np.unique(table.data[:,colCoord+i]) for i in xrange(3)]
mincorner = np.array(map(min,coords))
maxcorner = np.array(map(max,coords))
grid = np.array(map(len,coords),'i')
size = grid/np.maximum(np.ones(3,'d'), grid-1.0) * (maxcorner-mincorner) # size from edge to edge = dim * n/(n-1)
size = np.where(grid > 1, size, min(size[grid > 1]/grid[grid > 1])) # spacing for grid==1 set to smallest among other spacings
N = grid.prod()
if N != len(table.data): errors.append('data count {} does not match grid {}x{}x{}.'.format(N,*grid))
if errors != []:
damask.util.croak(errors)
table.close(dismiss = True)
continue
# ------------------------------------------ assemble header ---------------------------------------
table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:]))
for coord in xrange(3):
label = '{}_{}_{}'.format(coord+1,options.defgrad,options.coords)
if np.any(scaling) != 1.0: label+='_{}_{}_{}'.format(scaling)
if options.undeformed: label+='_undeformed'
table.labels_append([label]) # extend ASCII header with new labels
table.head_write()
# ------------------------------------------ read deformation gradient field -----------------------
centroids = deformedCoordsFFT(table.data[:,colF:colF+9].reshape(grid[2],grid[1],grid[0],3,3),
options.undeformed)
# ------------------------------------------ process data ------------------------------------------
table.data_rewind()
for z in xrange(grid[2]):
for y in xrange(grid[1]):
for x in xrange(grid[0]):
table.data_read()
table.data_append(list(centroids[z,y,x,:]))
table.data_write()
# ------------------------------------------ output finalization -----------------------------------
table.close() # close ASCII tables

221
processing/post/addDisplacement.py Executable file
View File

@ -0,0 +1,221 @@
#!/usr/bin/env python
# -*- coding: UTF-8 no BOM -*-
import os,sys,math
import numpy as np
import scipy.ndimage
from optparse import OptionParser
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
#--------------------------------------------------------------------------------------------------
def cell2node(cellData,grid):
nodeData = 0.0
datalen = np.array(cellData.shape[3:]).prod()
for i in xrange(datalen):
node = scipy.ndimage.convolve(cellData.reshape(tuple(grid)+(datalen,))[...,i],
np.ones((2,2,2))/8., # 2x2x2 neighborhood of cells
mode = 'wrap',
origin = -1, # offset to have cell origin as center
) # now averaged at cell origins
node = np.append(node,node[np.newaxis,0,:,:,...],axis=0) # wrap along z
node = np.append(node,node[:,0,np.newaxis,:,...],axis=1) # wrap along y
node = np.append(node,node[:,:,0,np.newaxis,...],axis=2) # wrap along x
nodeData = node[...,np.newaxis] if i==0 else np.concatenate((nodeData,node[...,np.newaxis]),axis=-1)
return nodeData
#--------------------------------------------------------------------------------------------------
def displacementAvgFFT(F,grid,size,nodal=False,transformed=False):
"""calculate average cell center (or nodal) displacement for deformation gradient field specified in each grid cell"""
if nodal:
x, y, z = np.meshgrid(np.linspace(0,size[0],1+grid[0]),
np.linspace(0,size[1],1+grid[1]),
np.linspace(0,size[2],1+grid[2]),
indexing = 'ij')
else:
x, y, z = np.meshgrid(np.linspace(0,size[0],grid[0],endpoint=False),
np.linspace(0,size[1],grid[1],endpoint=False),
np.linspace(0,size[2],grid[2],endpoint=False),
indexing = 'ij')
origCoords = np.concatenate((z[:,:,:,None],y[:,:,:,None],x[:,:,:,None]),axis = 3)
F_fourier = F if transformed else np.fft.rfftn(F,axes=(0,1,2)) # transform or use provided data
Favg = np.real(F_fourier[0,0,0,:,:])/grid.prod() # take zero freq for average
avgDisplacement = np.einsum('ml,ijkl->ijkm',Favg-np.eye(3),origCoords) # dX = Favg.X
return avgDisplacement
#--------------------------------------------------------------------------------------------------
def displacementFluctFFT(F,grid,size,nodal=False,transformed=False):
"""calculate cell center (or nodal) displacement for deformation gradient field specified in each grid cell"""
integrator = 0.5j * size / math.pi
kk, kj, ki = np.meshgrid(np.where(np.arange(grid[2])>grid[2]//2,np.arange(grid[2])-grid[2],np.arange(grid[2])),
np.where(np.arange(grid[1])>grid[1]//2,np.arange(grid[1])-grid[1],np.arange(grid[1])),
np.arange(grid[0]//2+1),
indexing = 'ij')
k_s = np.concatenate((ki[:,:,:,None],kj[:,:,:,None],kk[:,:,:,None]),axis = 3)
k_sSquared = np.einsum('...l,...l',k_s,k_s)
k_sSquared[0,0,0] = 1.0 # ignore global average frequency
#--------------------------------------------------------------------------------------------------
# integration in Fourier space
displacement_fourier = -np.einsum('ijkml,ijkl,l->ijkm',
F if transformed else np.fft.rfftn(F,axes=(0,1,2)),
k_s,
integrator,
) / k_sSquared[...,np.newaxis]
#--------------------------------------------------------------------------------------------------
# backtransformation to real space
displacement = np.fft.irfftn(displacement_fourier,grid,axes=(0,1,2))
return cell2node(displacement,grid) if nodal else displacement
# --------------------------------------------------------------------
# MAIN
# --------------------------------------------------------------------
parser = OptionParser(option_class=damask.extendableOption, usage='%prog options file[s]', description = """
Add displacments resulting from deformation gradient field.
Operates on periodic three-dimensional x,y,z-ordered data sets.
Outputs at cell centers or cell nodes (into separate file).
""", version = scriptID)
parser.add_option('-f',
'--defgrad',
dest = 'defgrad',
metavar = 'string',
help = 'column label of deformation gradient [%default]')
parser.add_option('-p',
'--pos', '--position',
dest = 'coords',
metavar = 'string',
help = 'label of coordinates [%default]')
parser.add_option('--nodal',
dest = 'nodal',
action = 'store_true',
help = 'output nodal (instad of cell-centered) displacements')
parser.set_defaults(defgrad = 'f',
coords = 'pos',
nodal = False,
)
(options,filenames) = parser.parse_args()
# --- loop over input files -------------------------------------------------------------------------
if filenames == []: filenames = [None]
for name in filenames:
try: table = damask.ASCIItable(name = name,
outname = (os.path.splitext(name)[0] +
'_nodal' +
os.path.splitext(name)[1]) if (options.nodal and name) else None,
buffered = False)
except: continue
damask.util.report(scriptName,name)
# ------------------------------------------ read header ------------------------------------------
table.head_read()
# ------------------------------------------ sanity checks ----------------------------------------
errors = []
remarks = []
if table.label_dimension(options.defgrad) != 9:
errors.append('deformation gradient "{}" is not a 3x3 tensor.'.format(options.defgrad))
coordDim = table.label_dimension(options.coords)
if not 3 >= coordDim >= 1:
errors.append('coordinates "{}" need to have one, two, or three dimensions.'.format(options.coords))
elif coordDim < 3:
remarks.append('appending {} dimension{} to coordinates "{}"...'.format(3-coordDim,
's' if coordDim < 2 else '',
options.coords))
if remarks != []: damask.util.croak(remarks)
if errors != []:
damask.util.croak(errors)
table.close(dismiss=True)
continue
# --------------- figure out size and grid ---------------------------------------------------------
table.data_readArray([options.defgrad,options.coords])
table.data_rewind()
if len(table.data.shape) < 2: table.data.shape += (1,) # expand to 2D shape
if table.data[:,9:].shape[1] < 3:
table.data = np.hstack((table.data,
np.zeros((table.data.shape[0],
3-table.data[:,9:].shape[1]),dtype='f'))) # fill coords up to 3D with zeros
coords = [np.unique(table.data[:,9+i]) for i in xrange(3)]
mincorner = np.array(map(min,coords))
maxcorner = np.array(map(max,coords))
grid = np.array(map(len,coords),'i')
size = grid/np.maximum(np.ones(3,'d'), grid-1.0) * (maxcorner-mincorner) # size from edge to edge = dim * n/(n-1)
size = np.where(grid > 1, size, min(size[grid > 1]/grid[grid > 1])) # spacing for grid==1 set to smallest among other spacings
N = grid.prod()
if N != len(table.data): errors.append('data count {} does not match grid {}x{}x{}.'.format(N,*grid))
if errors != []:
damask.util.croak(errors)
table.close(dismiss = True)
continue
# ------------------------------------------ process data ------------------------------------------
F_fourier = np.fft.rfftn(table.data[:,:9].reshape(grid[2],grid[1],grid[0],3,3),axes=(0,1,2)) # perform transform only once...
displacement = displacementFluctFFT(F_fourier,grid,size,options.nodal,transformed=True)
avgDisplacement = displacementAvgFFT (F_fourier,grid,size,options.nodal,transformed=True)
# ------------------------------------------ assemble header ---------------------------------------
if options.nodal:
table.info_clear()
table.labels_clear()
table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:]))
table.labels_append((['{}_pos' .format(i+1) for i in xrange(3)] if options.nodal else []) +
['{}_avg({}).{}' .format(i+1,options.defgrad,options.coords) for i in xrange(3)] +
['{}_fluct({}).{}'.format(i+1,options.defgrad,options.coords) for i in xrange(3)] )
table.head_write()
# ------------------------------------------ output data -------------------------------------------
zrange = np.linspace(0,size[2],1+grid[2]) if options.nodal else xrange(grid[2])
yrange = np.linspace(0,size[1],1+grid[1]) if options.nodal else xrange(grid[1])
xrange = np.linspace(0,size[0],1+grid[0]) if options.nodal else xrange(grid[0])
for i,z in enumerate(zrange):
for j,y in enumerate(yrange):
for k,x in enumerate(xrange):
if options.nodal: table.data_clear()
else: table.data_read()
table.data_append([x,y,z] if options.nodal else [])
table.data_append(list(avgDisplacement[i,j,k,:]))
table.data_append(list( displacement[i,j,k,:]))
table.data_write()
# ------------------------------------------ output finalization -----------------------------------
table.close() # close ASCII tables

8
processing/post/addDivergence.py
View File

@ -58,17 +58,17 @@ Deals with both vector- and tensor-valued fields.
parser.add_option('-c','--coordinates', parser.add_option('-c','--coordinates',
dest = 'coords', dest = 'coords',
type = 'string', metavar = 'string', type = 'string', metavar = 'string',
help = 'column heading for coordinates [%default]') help = 'column label of coordinates [%default]')
parser.add_option('-v','--vector', parser.add_option('-v','--vector',
dest = 'vector', dest = 'vector',
action = 'extend', metavar = '<string LIST>', action = 'extend', metavar = '<string LIST>',
help = 'heading of columns containing vector field values') help = 'column label(s) of vector field values')
parser.add_option('-t','--tensor', parser.add_option('-t','--tensor',
dest = 'tensor', dest = 'tensor',
action = 'extend', metavar = '<string LIST>', action = 'extend', metavar = '<string LIST>',
help = 'heading of columns containing tensor field values') help = 'column label(s) of tensor field values')
parser.set_defaults(coords = 'ipinitialcoord', parser.set_defaults(coords = 'pos',
) )
(options,filenames) = parser.parse_args() (options,filenames) = parser.parse_args()

72
processing/post/addEuclideanDistance.py
View File

@ -88,20 +88,32 @@ Add column(s) containing Euclidean distance to grain structural features: bounda
""", version = scriptID) """, version = scriptID)
parser.add_option('-c','--coordinates', dest='coords', metavar='string', parser.add_option('-p',
help='column heading for coordinates [%default]') '--pos', '--position',
parser.add_option('-i','--identifier', dest='id', metavar = 'string', dest = 'coords', metavar = 'string',
help='heading of column containing grain identifier [%default]') help = 'label of coordinates [%default]')
parser.add_option('-t','--type', dest = 'type', action = 'extend', metavar = '<string LIST>', parser.add_option('-i',
help = 'feature type {%s} '%(', '.join(map(lambda x:'/'.join(x['names']),features))) ) '--id', '--identifier',
parser.add_option('-n','--neighborhood',dest='neighborhood', choices = neighborhoods.keys(), metavar = 'string', dest = 'id', metavar = 'string',
help = 'type of neighborhood [neumann] {%s}'%(', '.join(neighborhoods.keys()))) help='label of grain identifier [%default]')
parser.add_option('-s', '--scale', dest = 'scale', type = 'float', metavar='float', parser.add_option('-t',
'--type',
dest = 'type', action = 'extend', metavar = '<string LIST>',
help = 'feature type {{{}}} '.format(', '.join(map(lambda x:'/'.join(x['names']),features))) )
parser.add_option('-n',
'--neighborhood',
dest = 'neighborhood', choices = neighborhoods.keys(), metavar = 'string',
help = 'neighborhood type [neumann] {{{}}}'.format(', '.join(neighborhoods.keys())))
parser.add_option('-s',
'--scale',
dest = 'scale', type = 'float', metavar = 'float',
help = 'voxel size [%default]') help = 'voxel size [%default]')
parser.set_defaults(coords = 'ipinitialcoord')
parser.set_defaults(id = 'texture') parser.set_defaults(coords = 'pos',
parser.set_defaults(neighborhood = 'neumann') id = 'texture',
parser.set_defaults(scale = 1.0) neighborhood = 'neumann',
scale = 1.0,
)
(options,filenames) = parser.parse_args() (options,filenames) = parser.parse_args()
@ -110,7 +122,7 @@ if options.type is None:
if not set(options.type).issubset(set(list(itertools.chain(*map(lambda x: x['names'],features))))): if not set(options.type).issubset(set(list(itertools.chain(*map(lambda x: x['names'],features))))):
parser.error('type must be chosen from (%s).'%(', '.join(map(lambda x:'|'.join(x['names']),features))) ) parser.error('type must be chosen from (%s).'%(', '.join(map(lambda x:'|'.join(x['names']),features))) )
if 'biplane' in options.type and 'boundary' in options.type: if 'biplane' in options.type and 'boundary' in options.type:
parser.error("only one from aliases 'biplane' and 'boundary' possible.") parser.error('only one from aliases "biplane" and "boundary" possible.')
feature_list = [] feature_list = []
for i,feature in enumerate(features): for i,feature in enumerate(features):
@ -125,10 +137,8 @@ for i,feature in enumerate(features):
if filenames == []: filenames = [None] if filenames == []: filenames = [None]
for name in filenames: for name in filenames:
try: try: table = damask.ASCIItable(name = name, buffered = False)
table = damask.ASCIItable(name = name, buffered = False) except: continue
except:
continue
damask.util.report(scriptName,name) damask.util.report(scriptName,name)
# ------------------------------------------ read header ------------------------------------------ # ------------------------------------------ read header ------------------------------------------
@ -141,7 +151,9 @@ for name in filenames:
remarks = [] remarks = []
column = {} column = {}
if table.label_dimension(options.coords) != 3: errors.append('coordinates {} are not a vector.'.format(options.coords)) coordDim = table.label_dimension(options.coords)
if not 3 >= coordDim >= 1:
errors.append('coordinates "{}" need to have one, two, or three dimensions.'.format(options.coords))
else: coordCol = table.label_index(options.coords) else: coordCol = table.label_index(options.coords)
if table.label_dimension(options.id) != 1: errors.append('grain identifier {} not found.'.format(options.id)) if table.label_dimension(options.id) != 1: errors.append('grain identifier {} not found.'.format(options.id))
@ -164,18 +176,18 @@ for name in filenames:
table.data_readArray() table.data_readArray()
coords = [{},{},{}] coords = [np.unique(table.data[:,coordCol+i]) for i in xrange(coordDim)]
for i in xrange(len(table.data)): mincorner = np.array(map(min,coords))
for j in xrange(3): maxcorner = np.array(map(max,coords))
coords[j][str(table.data[i,coordCol+j])] = True grid = np.array(map(len,coords)+[1]*(3-len(coords)),'i')
grid = np.array(map(len,coords),'i')
size = grid/np.maximum(np.ones(3,'d'),grid-1.0)* \
np.array([max(map(float,coords[0].keys()))-min(map(float,coords[0].keys())),\
max(map(float,coords[1].keys()))-min(map(float,coords[1].keys())),\
max(map(float,coords[2].keys()))-min(map(float,coords[2].keys())),\
],'d') # size from bounding box, corrected for cell-centeredness
size = np.where(grid > 1, size, min(size[grid > 1]/grid[grid > 1])) # spacing for grid==1 set to smallest among other spacings N = grid.prod()
if N != len(table.data): errors.append('data count {} does not match grid {}.'.format(N,'x'.join(map(str,grid))))
if errors != []:
damask.util.croak(errors)
table.close(dismiss = True)
continue
# ------------------------------------------ process value field ----------------------------------- # ------------------------------------------ process value field -----------------------------------

10
processing/post/addGradient.py
View File

@ -9,7 +9,9 @@ import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0] scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version]) scriptID = ' '.join([scriptName,damask.version])
#--------------------------------------------------------------------------------------------------
def gradFFT(geomdim,field): def gradFFT(geomdim,field):
grid = np.array(np.shape(field)[2::-1]) grid = np.array(np.shape(field)[2::-1])
N = grid.prod() # field size N = grid.prod() # field size
n = np.array(np.shape(field)[3:]).prod() # data size n = np.array(np.shape(field)[3:]).prod() # data size
@ -61,17 +63,17 @@ Deals with both vector- and scalar fields.
parser.add_option('-c','--coordinates', parser.add_option('-c','--coordinates',
dest = 'coords', dest = 'coords',
type = 'string', metavar='string', type = 'string', metavar='string',
help = 'column heading for coordinates [%default]') help = 'column label of coordinates [%default]')
parser.add_option('-v','--vector', parser.add_option('-v','--vector',
dest = 'vector', dest = 'vector',
action = 'extend', metavar = '<string LIST>', action = 'extend', metavar = '<string LIST>',
help = 'heading of columns containing vector field values') help = 'column label(s) of vector field values')
parser.add_option('-s','--scalar', parser.add_option('-s','--scalar',
dest = 'scalar', dest = 'scalar',
action = 'extend', metavar = '<string LIST>', action = 'extend', metavar = '<string LIST>',
help = 'heading of columns containing scalar field values') help = 'column label(s) of scalar field values')
parser.set_defaults(coords = 'ipinitialcoord', parser.set_defaults(coords = 'pos',
) )
(options,filenames) = parser.parse_args() (options,filenames) = parser.parse_args()

139
processing/post/addGrainID.py
View File

@ -5,7 +5,6 @@ import numpy as np
import damask import damask
from optparse import OptionParser from optparse import OptionParser
from scipy import spatial from scipy import spatial
from collections import defaultdict
scriptName = os.path.splitext(os.path.basename(__file__))[0] scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version]) scriptID = ' '.join([scriptName,damask.version])
@ -16,53 +15,61 @@ Add grain index based on similiarity of crystal lattice orientation.
""", version = scriptID) """, version = scriptID)
parser.add_option('-r', '--radius', parser.add_option('-r',
'--radius',
dest = 'radius', dest = 'radius',
type = 'float', metavar = 'float', type = 'float', metavar = 'float',
help = 'search radius') help = 'search radius')
parser.add_option('-d', '--disorientation', parser.add_option('-d',
'--disorientation',
dest = 'disorientation', dest = 'disorientation',
type = 'float', metavar = 'float', type = 'float', metavar = 'float',
help = 'disorientation threshold per grain [%default] (degrees)') help = 'disorientation threshold in degrees [%default]')
parser.add_option('-s', '--symmetry', parser.add_option('-s',
'--symmetry',
dest = 'symmetry', dest = 'symmetry',
type = 'string', metavar = 'string', type = 'string', metavar = 'string',
help = 'crystal symmetry [%default]') help = 'crystal symmetry [%default]')
parser.add_option('-e', '--eulers', parser.add_option('-e',
'--eulers',
dest = 'eulers', dest = 'eulers',
type = 'string', metavar = 'string', type = 'string', metavar = 'string',
help = 'Euler angles') help = 'label of Euler angles')
parser.add_option('--degrees', parser.add_option('--degrees',
dest = 'degrees', dest = 'degrees',
action = 'store_true', action = 'store_true',
help = 'Euler angles are given in degrees [%default]') help = 'Euler angles are given in degrees [%default]')
parser.add_option('-m', '--matrix', parser.add_option('-m',
'--matrix',
dest = 'matrix', dest = 'matrix',
type = 'string', metavar = 'string', type = 'string', metavar = 'string',
help = 'orientation matrix') help = 'label of orientation matrix')
parser.add_option('-a', parser.add_option('-a',
dest = 'a', dest = 'a',
type = 'string', metavar = 'string', type = 'string', metavar = 'string',
help = 'crystal frame a vector') help = 'label of crystal frame a vector')
parser.add_option('-b', parser.add_option('-b',
dest = 'b', dest = 'b',
type = 'string', metavar = 'string', type = 'string', metavar = 'string',
help = 'crystal frame b vector') help = 'label of crystal frame b vector')
parser.add_option('-c', parser.add_option('-c',
dest = 'c', dest = 'c',
type = 'string', metavar = 'string', type = 'string', metavar = 'string',
help = 'crystal frame c vector') help = 'label of crystal frame c vector')
parser.add_option('-q', '--quaternion', parser.add_option('-q',
'--quaternion',
dest = 'quaternion', dest = 'quaternion',
type = 'string', metavar = 'string', type = 'string', metavar = 'string',
help = 'quaternion') help = 'label of quaternion')
parser.add_option('-p', '--position', parser.add_option('-p',
dest = 'coords', '--pos', '--position',
dest = 'pos',
type = 'string', metavar = 'string', type = 'string', metavar = 'string',
help = 'spatial position of voxel [%default]') help = 'label of coordinates [%default]')
parser.set_defaults(symmetry = 'cubic', parser.set_defaults(disorientation = 5,
coords = 'pos', symmetry = 'cubic',
pos = 'pos',
degrees = False, degrees = False,
) )
@ -87,15 +94,14 @@ if np.sum(input) != 1: parser.error('needs exactly one input format.')
(options.quaternion,4,'quaternion'), (options.quaternion,4,'quaternion'),
][np.where(input)[0][0]] # select input label that was requested ][np.where(input)[0][0]] # select input label that was requested
toRadians = np.pi/180.0 if options.degrees else 1.0 # rescale degrees to radians toRadians = np.pi/180.0 if options.degrees else 1.0 # rescale degrees to radians
cos_disorientation = np.cos(options.disorientation/2.*toRadians) cos_disorientation = np.cos(np.radians(options.disorientation/2.)) # cos of half the disorientation angle
# --- loop over input files ------------------------------------------------------------------------- # --- loop over input files -------------------------------------------------------------------------
if filenames == []: filenames = [None] if filenames == []: filenames = [None]
for name in filenames: for name in filenames:
try: try: table = damask.ASCIItable(name = name,
table = damask.ASCIItable(name = name,
buffered = False) buffered = False)
except: continue except: continue
damask.util.report(scriptName,name) damask.util.report(scriptName,name)
@ -109,8 +115,10 @@ for name in filenames:
errors = [] errors = []
remarks = [] remarks = []
if table.label_dimension(options.coords) != 3: errors.append('coordinates {} are not a vector.'.format(options.coords)) if not 3 >= table.label_dimension(options.pos) >= 1:
if not np.all(table.label_dimension(label) == dim): errors.append('input {} does not have dimension {}.'.format(label,dim)) errors.append('coordinates "{}" need to have one, two, or three dimensions.'.format(options.pos))
if not np.all(table.label_dimension(label) == dim):
errors.append('input "{}" does not have dimension {}.'.format(label,dim))
else: column = table.label_index(label) else: column = table.label_index(label)
if remarks != []: damask.util.croak(remarks) if remarks != []: damask.util.croak(remarks)
@ -122,8 +130,10 @@ for name in filenames:
# ------------------------------------------ assemble header --------------------------------------- # ------------------------------------------ assemble header ---------------------------------------
table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:])) table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:]))
table.labels_append('grainID_{}@{}'.format(label, table.labels_append('grainID_{}@{:g}'.format('+'.join(label)
options.disorientation if options.degrees else np.degrees(options.disorientation))) # report orientation source and disorientation in degrees if isinstance(label, (list,tuple))
else label,
options.disorientation)) # report orientation source and disorientation
table.head_write() table.head_write()
# ------------------------------------------ process data ------------------------------------------ # ------------------------------------------ process data ------------------------------------------
@ -137,7 +147,7 @@ for name in filenames:
bg.set_message('reading positions...') bg.set_message('reading positions...')
table.data_readArray(options.coords) # read position vectors table.data_readArray(options.pos) # read position vectors
grainID = -np.ones(len(table.data),dtype=int) grainID = -np.ones(len(table.data),dtype=int)
start = tick = time.clock() start = tick = time.clock()
@ -162,7 +172,7 @@ for name in filenames:
time_delta = (time.clock()-tick) * (len(grainID) - p) / p time_delta = (time.clock()-tick) * (len(grainID) - p) / p
bg.set_message('(%02i:%02i:%02i) processing point %i of %i (grain count %i)...'\ bg.set_message('(%02i:%02i:%02i) processing point %i of %i (grain count %i)...'\
%(time_delta//3600,time_delta%3600//60,time_delta%60,p,len(grainID),len(orientations))) %(time_delta//3600,time_delta%3600//60,time_delta%60,p,len(grainID),np.count_nonzero(memberCounts)))
if inputtype == 'eulers': if inputtype == 'eulers':
o = damask.Orientation(Eulers = np.array(map(float,table.data[column:column+3]))*toRadians, o = damask.Orientation(Eulers = np.array(map(float,table.data[column:column+3]))*toRadians,
@ -180,83 +190,50 @@ for name in filenames:
symmetry = options.symmetry).reduced() symmetry = options.symmetry).reduced()
matched = False matched = False
# check against last matched needs to be really picky. best would be to exclude jumps across the poke (checking distance between last and me?)
# when walking through neighborhood first check whether grainID of that point has already been tested, if yes, skip!
if matchedID != -1: # has matched before?
matched = (o.quaternion.conjugated() * orientations[matchedID].quaternion).w > cos_disorientation
if not matched:
alreadyChecked = {} alreadyChecked = {}
candidates = []
bestDisorientation = damask.Quaternion([0,0,0,1]) # initialize to 180 deg rotation as worst case bestDisorientation = damask.Quaternion([0,0,0,1]) # initialize to 180 deg rotation as worst case
for i in kdtree.query_ball_point(kdtree.data[p],options.radius): # check all neighboring points for i in kdtree.query_ball_point(kdtree.data[p],options.radius): # check all neighboring points
gID = grainID[i] gID = grainID[i]
if gID != -1 and gID not in alreadyChecked: # indexed point belonging to a grain not yet tested? if gID != -1 and gID not in alreadyChecked: # indexed point belonging to a grain not yet tested?
alreadyChecked[gID] = True # remember not to check again alreadyChecked[gID] = True # remember not to check again
disorientation = o.disorientation(orientations[gID],SST = False)[0] # compare against other orientation disorientation = o.disorientation(orientations[gID],SST = False)[0] # compare against other orientation
if disorientation.quaternion.w > cos_disorientation and \ if disorientation.quaternion.w > cos_disorientation: # within threshold ...
disorientation.quaternion.w >= bestDisorientation.w: # within threshold and betterthan current best? candidates.append(gID) # remember as potential candidate
if disorientation.quaternion.w >= bestDisorientation.w: # ... and better than current best?
matched = True matched = True
matchedID = gID # remember that grain matchedID = gID # remember that grain
bestDisorientation = disorientation.quaternion bestDisorientation = disorientation.quaternion
if not matched: # no match -> new grain found if matched: # did match existing grain
memberCounts += [1] # start new membership counter memberCounts[matchedID] += 1
if len(candidates) > 1: # ambiguity in grain identification?
largestGrain = sorted(candidates,key=lambda x:memberCounts[x])[-1] # find largest among potential candidate grains
matchedID = largestGrain
for c in [c for c in candidates if c != largestGrain]: # loop over smaller candidates
memberCounts[largestGrain] += memberCounts[c] # reassign member count of smaller to largest
memberCounts[c] = 0
grainID = np.where(np.in1d(grainID,candidates), largestGrain, grainID) # relabel grid points of smaller candidates as largest one
else: # no match -> new grain found
orientations += [o] # initialize with current orientation orientations += [o] # initialize with current orientation
memberCounts += [1] # start new membership counter
matchedID = g matchedID = g
g += 1 # increment grain counter g += 1 # increment grain counter
else: # did match existing grain
memberCounts[matchedID] += 1
grainID[p] = matchedID # remember grain index assigned to point grainID[p] = matchedID # remember grain index assigned to point
p += 1 # increment point p += 1 # increment point
bg.set_message('identifying similar orientations among {} grains...'.format(len(orientations))) grainIDs = np.where(np.array(memberCounts) > 0)[0] # identify "live" grain identifiers
packingMap = dict(zip(list(grainIDs),range(len(grainIDs)))) # map to condense into consecutive IDs
memberCounts = np.array(memberCounts)
similarOrientations = [[] for i in xrange(len(orientations))]
for i,orientation in enumerate(orientations[:-1]): # compare each identified orientation...
for j in xrange(i+1,len(orientations)): # ...against all others that were defined afterwards
if orientation.disorientation(orientations[j],SST = False)[0].quaternion.w > cos_disorientation: # similar orientations in both grainIDs?
similarOrientations[i].append(j) # remember in upper triangle...
similarOrientations[j].append(i) # ...and lower triangle of matrix
if similarOrientations[i] != []:
bg.set_message('grainID {} is as: {}'.format(i,' '.join(map(str,similarOrientations[i]))))
stillShifting = True
while stillShifting:
stillShifting = False
tick = time.clock()
for p,gID in enumerate(grainID): # walk through all points
if p > 0 and p % 1000 == 0:
time_delta = (time.clock()-tick) * (len(grainID) - p) / p
bg.set_message('(%02i:%02i:%02i) shifting ID of point %i out of %i (grain count %i)...'
%(time_delta//3600,time_delta%3600//60,time_delta%60,p,len(grainID),len(orientations)))
if similarOrientations[gID] != []: # orientation of my grainID is similar to someone else?
similarNeighbors = defaultdict(int) # frequency of neighboring grainIDs sharing my orientation
for i in kdtree.query_ball_point(kdtree.data[p],options.radius): # check all neighboring point
if grainID[i] in similarOrientations[gID]: # neighboring point shares my orientation?
similarNeighbors[grainID[i]] += 1 # remember its grainID
if similarNeighbors != {}: # found similar orientation(s) in neighborhood
candidates = np.array([gID]+similarNeighbors.keys()) # possible replacement grainIDs for me
grainID[p] = candidates[np.argsort(memberCounts[candidates])[-1]] # adopt ID that is most frequent in overall dataset
memberCounts[gID] -= 1 # my former ID loses one fellow
memberCounts[grainID[p]] += 1 # my new ID gains one fellow
bg.set_message('{}:{} --> {}'.format(p,gID,grainID[p])) # report switch of grainID
stillShifting = True
table.data_rewind() table.data_rewind()
outputAlive = True outputAlive = True
p = 0 p = 0
while outputAlive and table.data_read(): # read next data line of ASCII table while outputAlive and table.data_read(): # read next data line of ASCII table
table.data_append(1+grainID[p]) # add grain ID table.data_append(1+packingMap[grainID[p]]) # add (condensed) grain ID
outputAlive = table.data_write() # output processed line outputAlive = table.data_write() # output processed line
p += 1 p += 1

56
processing/post/addInfo.py Executable file
View File

@ -0,0 +1,56 @@
#!/usr/bin/env python
# -*- coding: UTF-8 no BOM -*-
import os
from optparse import OptionParser
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
# --------------------------------------------------------------------
# MAIN
# --------------------------------------------------------------------
parser = OptionParser(option_class=damask.extendableOption, usage='%prog options file[s]', description = """
Add info lines to ASCIItable header.
""", version = scriptID)
parser.add_option('-i',
'--info',
dest = 'info', action = 'extend', metavar = '<string LIST>',
help = 'items to add')
parser.set_defaults(info = [],
)
(options,filenames) = parser.parse_args()
# --- loop over input files ------------------------------------------------------------------------
if filenames == []: filenames = [None]
for name in filenames:
try: table = damask.ASCIItable(name = name,
buffered = False)
except: continue
damask.util.report(scriptName,name)
# ------------------------------------------ assemble header ---------------------------------------
table.head_read()
table.info_append(options.info)
table.head_write()
# ------------------------------------------ pass through data -------------------------------------
outputAlive = True
while outputAlive and table.data_read(): # read next data line of ASCII table
outputAlive = table.data_write() # output processed line
# ------------------------------------------ output finalization -----------------------------------
table.close() # close ASCII tables

4
processing/post/addSchmidfactors.py
View File

@ -239,7 +239,9 @@ for name in filenames:
# ------------------------------------------ assemble header --------------------------------------- # ------------------------------------------ assemble header ---------------------------------------
table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:])) table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:]))
table.labels_append(['{id}_S[{direction[0]:.1g}_{direction[1]:.1g}_{direction[2]:.1g}]({normal[0]:.1g}_{normal[1]:.1g}_{normal[2]:.1g})'\ table.labels_append(['{id}_'
'S[{direction[0]:.1g}_{direction[1]:.1g}_{direction[2]:.1g}]'
'({normal[0]:.1g}_{normal[1]:.1g}_{normal[2]:.1g})'\
.format( id = i+1, .format( id = i+1,
normal = theNormal, normal = theNormal,
direction = theDirection, direction = theDirection,

9
processing/post/averageDown.py
View File

@ -22,7 +22,7 @@ Average each data block of size 'packing' into single values thus reducing the f
parser.add_option('-c','--coordinates', parser.add_option('-c','--coordinates',
dest = 'coords', dest = 'coords',
type = 'string', metavar = 'string', type = 'string', metavar = 'string',
help = 'column heading for coordinates [%default]') help = 'column label of coordinates [%default]')
parser.add_option('-p','--packing', parser.add_option('-p','--packing',
dest = 'packing', dest = 'packing',
type = 'int', nargs = 3, metavar = 'int int int', type = 'int', nargs = 3, metavar = 'int int int',
@ -39,7 +39,7 @@ parser.add_option('-s', '--size',
dest = 'size', dest = 'size',
type = 'float', nargs = 3, metavar = 'float float float', type = 'float', nargs = 3, metavar = 'float float float',
help = 'size in x,y,z [autodetect]') help = 'size in x,y,z [autodetect]')
parser.set_defaults(coords = 'ipinitialcoord', parser.set_defaults(coords = 'pos',
packing = (2,2,2), packing = (2,2,2),
shift = (0,0,0), shift = (0,0,0),
grid = (0,0,0), grid = (0,0,0),
@ -59,8 +59,7 @@ if any(shift != 0): prefix += 'shift{:+}{:+}{:+}_'.format(*shift)
if filenames == []: filenames = [None] if filenames == []: filenames = [None]
for name in filenames: for name in filenames:
try: try: table = damask.ASCIItable(name = name,
table = damask.ASCIItable(name = name,
outname = os.path.join(os.path.dirname(name), outname = os.path.join(os.path.dirname(name),
prefix+os.path.basename(name)) if name else name, prefix+os.path.basename(name)) if name else name,
buffered = False) buffered = False)
@ -75,7 +74,6 @@ for name in filenames:
errors = [] errors = []
remarks = [] remarks = []
colCoord = None
if table.label_dimension(options.coords) != 3: errors.append('coordinates {} are not a vector.'.format(options.coords)) if table.label_dimension(options.coords) != 3: errors.append('coordinates {} are not a vector.'.format(options.coords))
else: colCoord = table.label_index(options.coords) else: colCoord = table.label_index(options.coords)
@ -86,7 +84,6 @@ for name in filenames:
table.close(dismiss = True) table.close(dismiss = True)
continue continue
# ------------------------------------------ assemble header --------------------------------------- # ------------------------------------------ assemble header ---------------------------------------
table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:])) table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:]))

10
processing/post/averageTable.py
View File

@ -37,9 +37,13 @@ if options.label is None:
if filenames == []: filenames = [None] if filenames == []: filenames = [None]
for name in filenames: for name in filenames:
try: damask.util.croak(name)
table = damask.ASCIItable(name = name,
outname = options.label+'_averaged_'+name if name else name, try: table = damask.ASCIItable(name = name,
outname = os.path.join(
os.path.split(name)[0],
options.label+'_averaged_'+os.path.split(name)[1]
) if name else name,
buffered = False) buffered = False)
except: continue except: continue
damask.util.report(scriptName,name) damask.util.report(scriptName,name)

81
processing/post/blowUp.py
View File

@ -19,35 +19,37 @@ to resolution*packing.
""", version = scriptID) """, version = scriptID)
parser.add_option('-c','--coordinates', dest='coords', metavar='string', parser.add_option('-c','--coordinates',
help='column heading for coordinates [%default]') dest = 'coords', metavar = 'string',
parser.add_option('-p','--packing', dest='packing', type='int', nargs=3, metavar='int int int', help = 'column label of coordinates [%default]')
parser.add_option('-p','--packing',
dest = 'packing', type = 'int', nargs = 3, metavar = 'int int int',
help = 'dimension of packed group [%default]') help = 'dimension of packed group [%default]')
parser.add_option('-g','--grid', dest='resolution', type='int', nargs=3, metavar='int int int', parser.add_option('-g','--grid',
dest = 'resolution', type = 'int', nargs = 3, metavar = 'int int int',
help = 'resolution in x,y,z [autodetect]') help = 'resolution in x,y,z [autodetect]')
parser.add_option('-s','--size', dest='dimension', type='float', nargs=3, metavar='int int int', parser.add_option('-s','--size',
dest = 'dimension', type = 'float', nargs = 3, metavar = 'int int int',
help = 'dimension in x,y,z [autodetect]') help = 'dimension in x,y,z [autodetect]')
parser.set_defaults(coords = 'ipinitialcoord') parser.set_defaults(coords = 'pos',
parser.set_defaults(packing = (2,2,2)) packing = (2,2,2),
parser.set_defaults(grid = (0,0,0)) grid = (0,0,0),
parser.set_defaults(size = (0.0,0.0,0.0)) size = (0.0,0.0,0.0),
)
(options,filenames) = parser.parse_args() (options,filenames) = parser.parse_args()
options.packing = np.array(options.packing) options.packing = np.array(options.packing)
prefix = 'blowUp%ix%ix%i_'%(options.packing[0],options.packing[1],options.packing[2]) prefix = 'blowUp{}x{}x{}_'.format(*options.packing)
# --- loop over input files ------------------------------------------------------------------------- # --- loop over input files -------------------------------------------------------------------------
if filenames == []: filenames = [None] if filenames == []: filenames = [None]
for name in filenames: for name in filenames:
try: try: table = damask.ASCIItable(name = name,
table = damask.ASCIItable(name = name,
outname = os.path.join(os.path.dirname(name), outname = os.path.join(os.path.dirname(name),
prefix+ \ prefix+os.path.basename(name)) if name else name,
os.path.basename(name)) if name else name,
buffered = False) buffered = False)
except: continue except: continue
damask.util.report(scriptName,name) damask.util.report(scriptName,name)
@ -55,48 +57,45 @@ for name in filenames:
# ------------------------------------------ read header ------------------------------------------ # ------------------------------------------ read header ------------------------------------------
table.head_read() table.head_read()
errors = []
# ------------------------------------------ sanity checks ---------------------------------------- # ------------------------------------------ sanity checks ----------------------------------------
if table.label_dimension(options.coords) != 3: errors = []
damask.util.croak('coordinates {} are not a vector.'.format(options.coords)) remarks = []
if table.label_dimension(options.coords) != 3: errors.append('coordinates {} are not a vector.'.format(options.coords))
else: colCoord = table.label_index(options.coords)
colElem = table.label_index('elem')
if remarks != []: damask.util.croak(remarks)
if errors != []:
damask.util.croak(errors)
table.close(dismiss = True) table.close(dismiss = True)
continue continue
else:
coordCol = table.label_index(options.coords)
# ------------------------------------------ assemble header --------------------------------------
table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:]))
# --------------- figure out size and grid --------------------------------------------------------- # --------------- figure out size and grid ---------------------------------------------------------
table.data_readArray() table.data_readArray(options.coords)
table.data_rewind()
coords = [{},{},{}] coords = [np.unique(table.data[:,i]) for i in xrange(3)]
for i in xrange(len(table.data)): mincorner = np.array(map(min,coords))
for j in xrange(3): maxcorner = np.array(map(max,coords))
coords[j][str(table.data[i,coordCol+j])] = True
grid = np.array(map(len,coords),'i') grid = np.array(map(len,coords),'i')
size = grid/np.maximum(np.ones(3,'d'),grid-1.0)* \ size = grid/np.maximum(np.ones(3,'d'), grid-1.0) * (maxcorner-mincorner) # size from edge to edge = dim * n/(n-1)
np.array([max(map(float,coords[0].keys()))-min(map(float,coords[0].keys())),\
max(map(float,coords[1].keys()))-min(map(float,coords[1].keys())),\
max(map(float,coords[2].keys()))-min(map(float,coords[2].keys())),\
],'d') # size from bounding box, corrected for cell-centeredness
size = np.where(grid > 1, size, min(size[grid > 1]/grid[grid > 1])) # spacing for grid==1 set to smallest among other spacings size = np.where(grid > 1, size, min(size[grid > 1]/grid[grid > 1])) # spacing for grid==1 set to smallest among other spacings
packing = np.array(options.packing,'i') packing = np.array(options.packing,'i')
outSize = grid*packing outSize = grid*packing
# ------------------------------------------ assemble header --------------------------------------- # ------------------------------------------ assemble header --------------------------------------
table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:]))
table.head_write() table.head_write()
# ------------------------------------------ process data ------------------------------------------- # ------------------------------------------ process data -------------------------------------------
table.data_rewind()
data = np.zeros(outSize.tolist()+[len(table.labels)]) data = np.zeros(outSize.tolist()+[len(table.labels)])
p = np.zeros(3,'i') p = np.zeros(3,'i')
@ -114,8 +113,8 @@ for name in filenames:
for c in xrange(outSize[2]): for c in xrange(outSize[2]):
for b in xrange(outSize[1]): for b in xrange(outSize[1]):
for a in xrange(outSize[0]): for a in xrange(outSize[0]):
data[a,b,c,coordCol:coordCol+3] = [a+0.5,b+0.5,c+0.5]*elementSize data[a,b,c,colCoord:colCoord+3] = [a+0.5,b+0.5,c+0.5]*elementSize
data[a,b,c,table.label_index('elem')] = elem if colElem != -1: data[a,b,c,colElem] = elem
table.data = data[a,b,c,:].tolist() table.data = data[a,b,c,:].tolist()
outputAlive = table.data_write() # output processed line outputAlive = table.data_write() # output processed line
elem += 1 elem += 1

142
processing/post/histogram.py Executable file
View File

@ -0,0 +1,142 @@
#!/usr/bin/env python
# -*- coding: UTF-8 no BOM -*-
import os,sys
import numpy as np
from optparse import OptionParser
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
# --------------------------------------------------------------------
# MAIN
# --------------------------------------------------------------------
parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """
Generate histogram of N bins in given data range.
""", version = scriptID)
parser.add_option('-d','--data',
dest = 'data',
type = 'string', metavar = 'string',
help = 'column heading for data')
parser.add_option('-w','--weights',
dest = 'weights',
type = 'string', metavar = 'string',
help = 'column heading for weights')
parser.add_option('--range',
dest = 'range',
type = 'float', nargs = 2, metavar = 'float float',
help = 'data range of histogram [min - max]')
parser.add_option('-N',
dest = 'N',
type = 'int', metavar = 'int',
help = 'number of bins')
parser.add_option('--density',
dest = 'density',
action = 'store_true',
help = 'report probability density')
parser.add_option('--logarithmic',
dest = 'log',
action = 'store_true',
help = 'logarithmically spaced bins')
parser.set_defaults(data = None,
weights = None,
range = None,
N = None,
density = False,
log = False,
)
(options,filenames) = parser.parse_args()
if not options.data: parser.error('no data specified.')
if not options.N: parser.error('no bin number specified.')
if options.log:
def forward(x):
return np.log(x)
def reverse(x):
return np.exp(x)
else:
def forward(x):
return x
def reverse(x):
return x
# --- loop over input files ------------------------------------------------------------------------
if filenames == []: filenames = [None]
for name in filenames:
try: table = damask.ASCIItable(name = name,
buffered = False,
readonly = True)
except: continue
damask.util.report(scriptName,name)
# ------------------------------------------ read header ------------------------------------------
table.head_read()
# ------------------------------------------ sanity checks ----------------------------------------
errors = []
remarks = []
if table.label_dimension(options.data) != 1: errors.append('data {} are not scalar.'.format(options.data))
if options.weights and \
table.label_dimension(options.data) != 1: errors.append('weights {} are not scalar.'.format(options.weights))
if remarks != []: damask.util.croak(remarks)
if errors != []:
damask.util.croak(errors)
table.close(dismiss = True)
continue
# --------------- read data ----------------------------------------------------------------
table.data_readArray([options.data,options.weights])
# --------------- auto range ---------------------------------------------------------------
if options.range is None:
rangeMin,rangeMax = min(table.data[:,0]),max(table.data[:,0])
else:
rangeMin,rangeMax = min(options.range),max(options.range)
# --------------- bin data ----------------------------------------------------------------
count,edges = np.histogram(table.data[:,0],
bins = reverse(forward(rangeMin) + np.arange(options.N+1) *
(forward(rangeMax)-forward(rangeMin))/options.N),
range = (rangeMin,rangeMax),
weights = None if options.weights is None else table.data[:,1],
density = options.density,
)
bincenter = reverse(forward(rangeMin) + (0.5+np.arange(options.N)) *
(forward(rangeMax)-forward(rangeMin))/options.N) # determine center of bins
# ------------------------------------------ assemble header ---------------------------------------
table.info_clear()
table.info_append([scriptID + '\t' + ' '.join(sys.argv[1:]),
scriptID + ':\t' +
'data range {} -- {}'.format(rangeMin,rangeMax) +
(' (log)' if options.log else ''),
])
table.labels_clear()
table.labels_append(['bincenter','count'])
table.head_write()
# ------------------------------------------ output result -----------------------------------------
table.data = np.squeeze(np.dstack((bincenter,count)))
table.data_writeArray()
# ------------------------------------------ output finalization -----------------------------------
table.close() # close ASCII tables

144
processing/post/marc_deformedGeometry.py
View File

@ -1,144 +0,0 @@
#!/usr/bin/env python
# -*- coding: UTF-8 no BOM -*-
import os,sys,shutil
import numpy as np
import damask
from optparse import OptionParser
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
# -----------------------------
# MAIN FUNCTION STARTS HERE
# -----------------------------
# --- input parsing
parser = OptionParser(usage='%prog [options] resultfile', description = """
Create vtk files for the (deformed) geometry that belongs to a .t16 (MSC.Marc) results file.
""", version = scriptID)
parser.add_option('-d','--dir', dest='dir', \
help='name of subdirectory to hold output [%default]')
parser.add_option('-r','--range', dest='range', type='int', nargs=3, \
help='range of positions (or increments) to output (start, end, step) [all]')
parser.add_option('--increments', action='store_true', dest='getIncrements', \
help='switch to increment range [%default]')
parser.add_option('-t','--type', dest='type', type='choice', choices=['ipbased','nodebased'], \
help='processed geometry type [ipbased and nodebased]')
parser.set_defaults(dir = 'vtk')
parser.set_defaults(getIncrements= False)
(options, files) = parser.parse_args()
# --- basic sanity checks
if files == []:
parser.print_help()
parser.error('no file specified...')
filename = os.path.splitext(files[0])[0]
if not os.path.exists(filename+'.t16'):
parser.print_help()
parser.error('invalid file "%s" specified...'%filename+'.t16')
if not options.type :
options.type = ['nodebased', 'ipbased']
else:
options.type = [options.type]
# --- more sanity checks
sys.path.append(damask.solver.Marc().libraryPath('../../'))
try:
import py_post
except:
print('error: no valid Mentat release found')
sys.exit(-1)
# --------------------------- open results file and initialize mesh ----------
p = py_post.post_open(filename+'.t16')
p.moveto(0)
Nnodes = p.nodes()
Nincrements = p.increments() - 1 # t16 contains one "virtual" increment (at 0)
if damask.core.mesh.mesh_init_postprocessing(filename+'.mesh') > 0:
print('error: init not successful')
sys.exit(-1)
Ncellnodes = damask.core.mesh.mesh_get_Ncellnodes()
unitlength = damask.core.mesh.mesh_get_unitlength()
# --------------------------- create output dir --------------------------------
dirname = os.path.abspath(os.path.join(os.path.dirname(filename),options.dir))
if not os.path.isdir(dirname):
os.mkdir(dirname,0755)
# --------------------------- get positions --------------------------------
incAtPosition = {}
positionOfInc = {}
for position in range(Nincrements):
p.moveto(position+1)
incAtPosition[position] = p.increment # remember "real" increment at this position
positionOfInc[p.increment] = position # remember position of "real" increment
if not options.range:
options.getIncrements = False
locations = range(Nincrements) # process all positions
else:
options.range = list(options.range) # convert to list
if options.getIncrements:
locations = [positionOfInc[x] for x in range(options.range[0],options.range[1]+1,options.range[2])
if x in positionOfInc]
else:
locations = range( max(0,options.range[0]),
min(Nincrements,options.range[1]+1),
options.range[2] )
increments = [incAtPosition[x] for x in locations] # build list of increments to process
# --------------------------- loop over positions --------------------------------
for incCount,position in enumerate(locations): # walk through locations
p.moveto(position+1) # wind to correct position
# --- get displacements
node_displacement = [[0,0,0] for i in range(Nnodes)]
for n in range(Nnodes):
if p.node_displacements():
node_displacement[n] = map(lambda x:x*unitlength,list(p.node_displacement(n)))
c = damask.core.mesh.mesh_build_cellnodes(np.array(node_displacement).T,Ncellnodes)
cellnode_displacement = [[c[i][n] for i in range(3)] for n in range(Ncellnodes)]
# --- append displacements to corresponding files
for geomtype in options.type:
outFilename = eval('"'+eval("'%%s_%%s_inc%%0%ii.vtk'%(math.log10(max(increments+[1]))+1)")\
+'"%(dirname + os.sep + os.path.split(filename)[1],geomtype,increments[incCount])')
print outFilename
shutil.copyfile('%s_%s.vtk'%(filename,geomtype),outFilename)
with open(outFilename,'a') as myfile:
myfile.write("POINT_DATA %i\n"%{'nodebased':Nnodes,'ipbased':Ncellnodes}[geomtype])
myfile.write("VECTORS displacement double\n")
coordinates = {'nodebased':node_displacement,'ipbased':cellnode_displacement}[geomtype]
for n in range({'nodebased':Nnodes,'ipbased':Ncellnodes}[geomtype]):
myfile.write("%.8e %.8e %.8e\n"%(coordinates[n][0],coordinates[n][1],coordinates[n][2]))
# --------------------------- DONE --------------------------------

421
processing/post/marc_extractData.py
View File

@ -1,421 +0,0 @@
#!/usr/bin/env python
# -*- coding: UTF-8 no BOM -*-
import os,sys,string,re,time
from optparse import OptionParser, OptionGroup
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
# -----------------------------
def ParseOutputFormat(filename,homogID,crystID,phaseID):
"""parse .output* files in order to get a list of outputs"""
myID = {'Homogenization': homogID,
'Crystallite': crystID,
'Constitutive': phaseID,
}
format = {}
for what in ['Homogenization','Crystallite','Constitutive']:
content = []
format[what] = {'outputs':{},'specials':{'brothers':[]}}
for prefix in ['']+map(str,range(1,17)):
if os.path.exists(prefix+filename+'.output'+what):
try:
file = open(prefix+filename+'.output'+what)
content = file.readlines()
file.close()
break
except:
pass
if content == []: continue # nothing found...
tag = ''
tagID = 0
for line in content:
if re.match("\s*$",line) or re.match("#",line): # skip blank lines and comments
continue
m = re.match("\[(.+)\]",line) # look for block indicator
if m: # next section
tag = m.group(1)
tagID += 1
format[what]['specials']['brothers'].append(tag)
if tag == myID[what] or (myID[what].isdigit() and tagID == int(myID[what])):
format[what]['specials']['_id'] = tagID
format[what]['outputs'] = []
tag = myID[what]
else: # data from section
if tag == myID[what]:
(output,length) = line.split()
output.lower()
if length.isdigit():
length = int(length)
if re.match("\((.+)\)",output): # special data, e.g. (Ngrains)
format[what]['specials'][output] = length
elif length > 0:
format[what]['outputs'].append([output,length])
if '_id' not in format[what]['specials']:
print "\nsection '%s' not found in <%s>"%(myID[what], what)
print '\n'.join(map(lambda x:' [%s]'%x, format[what]['specials']['brothers']))
return format
# -----------------------------
def ParsePostfile(p,filename, outputFormat, legacyFormat):
"""
parse postfile in order to get position and labels of outputs
needs "outputFormat" for mapping of output names to postfile output indices
"""
startVar = {True: 'GrainCount',
False:'HomogenizationCount'}
# --- build statistics
stat = { \
'IndexOfLabel': {}, \
'Title': p.title(), \
'Extrapolation': p.extrapolate, \
'NumberOfIncrements': p.increments() - 1, \
'NumberOfNodes': p.nodes(), \
'NumberOfNodalScalars': p.node_scalars(), \
'LabelOfNodalScalar': [None]*p.node_scalars() , \
'NumberOfElements': p.elements(), \
'NumberOfElementalScalars': p.element_scalars(), \
'LabelOfElementalScalar': [None]*p.element_scalars() , \
'NumberOfElementalTensors': p.element_tensors(), \
'LabelOfElementalTensor': [None]*p.element_tensors(), \
}
# --- find labels
for labelIndex in range(stat['NumberOfNodalScalars']):
label = p.node_scalar_label(labelIndex)
stat['IndexOfLabel'][label] = labelIndex
stat['LabelOfNodalScalar'][labelIndex] = label
for labelIndex in range(stat['NumberOfElementalScalars']):
label = p.element_scalar_label(labelIndex)
stat['IndexOfLabel'][label] = labelIndex
stat['LabelOfElementalScalar'][labelIndex] = label
for labelIndex in range(stat['NumberOfElementalTensors']):
label = p.element_tensor_label(labelIndex)
stat['IndexOfLabel'][label] = labelIndex
stat['LabelOfElementalTensor'][labelIndex] = label
if 'User Defined Variable 1' in stat['IndexOfLabel']: # output format without dedicated names?
stat['IndexOfLabel'][startVar[legacyFormat]] = stat['IndexOfLabel']['User Defined Variable 1'] # adjust first named entry
if startVar[legacyFormat] in stat['IndexOfLabel']: # does the result file contain relevant user defined output at all?
startIndex = stat['IndexOfLabel'][startVar[legacyFormat]]
stat['LabelOfElementalScalar'][startIndex] = startVar[legacyFormat]
# We now have to find a mapping for each output label as defined in the .output* files to the output position in the post file
# Since we know where the user defined outputs start ("startIndex"), we can simply assign increasing indices to the labels
# given in the .output* file
offset = 1
if legacyFormat:
stat['LabelOfElementalScalar'][startIndex + offset] = startVar[not legacyFormat] # add HomogenizationCount as second
offset += 1
for (name,N) in outputFormat['Homogenization']['outputs']:
for i in range(N):
label = {False: '%s'%( name),
True:'%i_%s'%(i+1,name)}[N > 1]
stat['IndexOfLabel'][label] = startIndex + offset
stat['LabelOfElementalScalar'][startIndex + offset] = label
offset += 1
if not legacyFormat:
stat['IndexOfLabel'][startVar[not legacyFormat]] = startIndex + offset
stat['LabelOfElementalScalar'][startIndex + offset] = startVar[not legacyFormat] # add GrainCount
offset += 1
if '(ngrains)' in outputFormat['Homogenization']['specials']:
for grain in range(outputFormat['Homogenization']['specials']['(ngrains)']):
stat['IndexOfLabel']['%i_CrystalliteCount'%(grain+1)] = startIndex + offset # report crystallite count
stat['LabelOfElementalScalar'][startIndex + offset] = '%i_CrystalliteCount'%(grain+1) # add GrainCount
offset += 1
for (name,N) in outputFormat['Crystallite']['outputs']: # add crystallite outputs
for i in range(N):
label = {False: '%i_%s'%(grain+1, name),
True:'%i_%i_%s'%(grain+1,i+1,name)}[N > 1]
stat['IndexOfLabel'][label] = startIndex + offset
stat['LabelOfElementalScalar'][startIndex + offset] = label
offset += 1
stat['IndexOfLabel']['%i_ConstitutiveCount'%(grain+1)] = startIndex + offset # report constitutive count
stat['LabelOfElementalScalar'][startIndex + offset] = '%i_ConstitutiveCount'%(grain+1) # add GrainCount
offset += 1
for (name,N) in outputFormat['Constitutive']['outputs']: # add constitutive outputs
for i in range(N):
label = {False: '%i_%s'%(grain+1, name),
True:'%i_%i_%s'%(grain+1,i+1,name)}[N > 1]
stat['IndexOfLabel'][label] = startIndex + offset
try:
stat['LabelOfElementalScalar'][startIndex + offset] = label
except IndexError:
print 'trying to assign %s at position %i+%i'%(label,startIndex,offset)
sys.exit(1)
offset += 1
return stat
# -----------------------------
def GetIncrementLocations(p,Nincrements,options):
"""get mapping between positions in postfile and increment number"""
incAtPosition = {}
positionOfInc = {}
for position in range(Nincrements):
p.moveto(position+1)
incAtPosition[position] = p.increment # remember "real" increment at this position
positionOfInc[p.increment] = position # remember position of "real" increment
if not options.range:
options.getIncrements = False
locations = range(Nincrements) # process all positions
else:
options.range = list(options.range) # convert to list
if options.getIncrements:
locations = [positionOfInc[x] for x in range(options.range[0],options.range[1]+1,options.range[2])
if x in positionOfInc]
else:
locations = range( max(0,options.range[0]),
min(Nincrements,options.range[1]+1),
options.range[2] )
increments = [incAtPosition[x] for x in locations] # build list of increments to process
return [increments,locations]
# -----------------------------
def SummarizePostfile(stat,where=sys.stdout):
where.write('\n\n')
where.write('title:\t%s'%stat['Title'] + '\n\n')
where.write('extraplation:\t%s'%stat['Extrapolation'] + '\n\n')
where.write('increments:\t%i'%(stat['NumberOfIncrements']) + '\n\n')
where.write('nodes:\t%i'%stat['NumberOfNodes'] + '\n\n')
where.write('elements:\t%i'%stat['NumberOfElements'] + '\n\n')
where.write('nodal scalars:\t%i'%stat['NumberOfNodalScalars'] + '\n\n '\
+'\n '.join(stat['LabelOfNodalScalar']) + '\n\n')
where.write('elemental scalars:\t%i'%stat['NumberOfElementalScalars'] + '\n\n '\
+ '\n '.join(stat['LabelOfElementalScalar']) + '\n\n')
where.write('elemental tensors:\t%i'%stat['NumberOfElementalTensors'] + '\n\n '\
+ '\n '.join(stat['LabelOfElementalTensor']) + '\n\n')
return True
# -----------------------------
def SummarizeOutputfile(format,where=sys.stdout):
where.write('\nUser Defined Outputs')
for what in format.keys():
where.write('\n\n %s:'%what)
for output in format[what]['outputs']:
where.write('\n %s'%output)
return True
# -----------------------------
def writeHeader(myfile,stat,geomtype):
myfile.write('2\theader\n')
myfile.write(string.replace('$Id$','\n','\\n')+
'\t' + ' '.join(sys.argv[1:]) + '\n')
if geomtype == 'nodebased':
myfile.write('node')
for i in range(stat['NumberOfNodalScalars']):
myfile.write('\t%s'%''.join(stat['LabelOfNodalScalar'][i].split()))
elif geomtype == 'ipbased':
myfile.write('elem\tip')
for i in range(stat['NumberOfElementalScalars']):
myfile.write('\t%s'%''.join(stat['LabelOfElementalScalar'][i].split()))
myfile.write('\n')
return True
# --------------------------------------------------------------------
# MAIN
# --------------------------------------------------------------------
parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """
Extract data from a .t16 (MSC.Marc) results file.
""", version = scriptID)
parser.add_option('-i','--info', action='store_true', dest='info', \
help='list contents of resultfile [%default]')
parser.add_option('-l','--legacy', action='store_true', dest='legacy', \
help='legacy user result block (starts with GrainCount) [%default]')
parser.add_option('-d','--dir', dest='dir', \
help='name of subdirectory to hold output [%default]')
parser.add_option('-r','--range', dest='range', type='int', nargs=3, \
help='range of positions (or increments) to output (start, end, step) [all]')
parser.add_option('--increments', action='store_true', dest='getIncrements', \
help='switch to increment range [%default]')
parser.add_option('-t','--type', dest='type', type='choice', choices=['ipbased','nodebased'], \
help='processed geometry type [ipbased and nodebased]')
group_material = OptionGroup(parser,'Material identifier')
group_material.add_option('--homogenization', dest='homog', \
help='homogenization identifier (as string or integer [%default])', metavar='<ID>')
group_material.add_option('--crystallite', dest='cryst', \
help='crystallite identifier (as string or integer [%default])', metavar='<ID>')
group_material.add_option('--phase', dest='phase', \
help='phase identifier (as string or integer [%default])', metavar='<ID>')
parser.add_option_group(group_material)
parser.set_defaults(info = False)
parser.set_defaults(legacy = False)
parser.set_defaults(dir = 'vtk')
parser.set_defaults(getIncrements= False)
parser.set_defaults(homog = '1')
parser.set_defaults(cryst = '1')
parser.set_defaults(phase = '1')
(options, files) = parser.parse_args()
# --- sanity checks
if files == []:
parser.print_help()
parser.error('no file specified...')
filename = os.path.splitext(files[0])[0]
if not os.path.exists(filename+'.t16'):
parser.print_help()
parser.error('invalid file "%s" specified...'%filename+'.t16')
sys.path.append(damask.solver.Marc().libraryPath('../../'))
try:
import py_post
except:
print('error: no valid Mentat release found')
sys.exit(-1)
if not options.type :
options.type = ['nodebased', 'ipbased']
else:
options.type = [options.type]
# --- initialize mesh data
if damask.core.mesh.mesh_init_postprocessing(filename+'.mesh'):
print('error: init not successful')
sys.exit(-1)
# --- check if ip data available for all elements; if not, then .t19 file is required
p = py_post.post_open(filename+'.t16')
asciiFile = False
p.moveto(1)
for e in range(p.elements()):
if not damask.core.mesh.mesh_get_nodeAtIP(str(p.element(e).type),1):
if os.path.exists(filename+'.t19'):
p.close()
p = py_post.post_open(filename+'.t19')
asciiFile = True
break
# --- parse *.output and *.t16 file
outputFormat = ParseOutputFormat(filename,options.homog,options.cryst,options.phase)
p.moveto(1)
p.extrapolation('translate')
stat = ParsePostfile(p,filename,outputFormat,options.legacy)
# --- output info
if options.info:
print '\n\nMentat release %s'%damask.solver.Marc().version('../../')
SummarizePostfile(stat)
SummarizeOutputfile(outputFormat)
sys.exit(0)
# --- create output dir
dirname = os.path.abspath(os.path.join(os.path.dirname(filename),options.dir))
if not os.path.isdir(dirname):
os.mkdir(dirname,0755)
# --- get positions
[increments,locations] = GetIncrementLocations(p,stat['NumberOfIncrements'],options)
# --- loop over positions
time_start = time.time()
for incCount,position in enumerate(locations): # walk through locations
p.moveto(position+1) # wind to correct position
time_delta = (float(len(locations)) / float(incCount+1) - 1.0) * (time.time() - time_start)
sys.stdout.write("\r(%02i:%02i:%02i) processing increment %i of %i..."\
%(time_delta//3600,time_delta%3600//60,time_delta%60,incCount+1,len(locations)))
sys.stdout.flush()
# --- write header
outFilename = {}
for geomtype in options.type:
outFilename[geomtype] = eval('"'+eval("'%%s_%%s_inc%%0%ii.txt'%(math.log10(max(increments+[1]))+1)")\
+'"%(dirname + os.sep + os.path.split(filename)[1],geomtype,increments[incCount])')
with open(outFilename[geomtype],'w') as myfile:
writeHeader(myfile,stat,geomtype)
# --- write node based data
if geomtype == 'nodebased':
for n in range(stat['NumberOfNodes']):
myfile.write(str(n))
for l in range(stat['NumberOfNodalScalars']):
myfile.write('\t'+str(p.node_scalar(n,l)))
myfile.write('\n')
# --- write ip based data
elif geomtype == 'ipbased':
for e in range(stat['NumberOfElements']):
if asciiFile:
print 'ascii postfile not yet supported'
sys.exit(-1)
else:
ipData = [[]]
for l in range(stat['NumberOfElementalScalars']):
data = p.element_scalar(e,l)
for i in range(len(data)): # at least as many nodes as ips
node = damask.core.mesh.mesh_get_nodeAtIP(str(p.element(e).type),i+1) # fortran indexing starts at 1
if not node: break # no more ips
while i >= len(ipData): ipData.append([])
ipData[i].extend([data[node-1].value]) # python indexing starts at 0
for i in range(len(ipData)):
myfile.write('\t'.join(map(str,[e,i]+ipData[i]))+'\n')
p.close()
sys.stdout.write("\n")

7
processing/post/postResults.py
View File

@ -1003,11 +1003,8 @@ fileOpen = False
assembleHeader = True assembleHeader = True
header = [] header = []
standard = ['inc'] + \ standard = ['inc'] + \
{True: ['time'], (['time'] if options.time else []) + \
False:[]}[options.time] + \ ['elem','node','ip','grain','1_pos','2_pos','3_pos']
['elem','node','ip','grain'] + \
{True: ['1_nodeinitialcoord','2_nodeinitialcoord','3_nodeinitialcoord'],
False:['1_ipinitialcoord','2_ipinitialcoord','3_ipinitialcoord']}[options.nodalScalar != []]
# --------------------------- loop over positions -------------------------------- # --------------------------- loop over positions --------------------------------

12
processing/post/reLabel.py
View File

@ -1,7 +1,7 @@
#!/usr/bin/env python #!/usr/bin/env python
# -*- coding: UTF-8 no BOM -*- # -*- coding: UTF-8 no BOM -*-
import os,sys import os,sys,re
import damask import damask
from optparse import OptionParser from optparse import OptionParser
@ -32,13 +32,16 @@ parser.set_defaults(label = [],
(options,filenames) = parser.parse_args() (options,filenames) = parser.parse_args()
pattern = [re.compile('^()(.+)$'), # label pattern for scalar
re.compile('^(\d+_)?(.+)$'), # label pattern for multidimension
]
# --- loop over input files ------------------------------------------------------------------------- # --- loop over input files -------------------------------------------------------------------------
if filenames == []: filenames = [None] if filenames == []: filenames = [None]
for name in filenames: for name in filenames:
try: try: table = damask.ASCIItable(name = name,
table = damask.ASCIItable(name = name,
buffered = False) buffered = False)
except: continue except: continue
damask.util.report(scriptName,name) damask.util.report(scriptName,name)
@ -63,8 +66,9 @@ for name in filenames:
for i,index in enumerate(indices): for i,index in enumerate(indices):
if index == -1: remarks.append('label {} not present...'.format(options.label[i])) if index == -1: remarks.append('label {} not present...'.format(options.label[i]))
else: else:
m = pattern[dimensions[i]>1].match(table.labels[index]) # isolate label name
for j in xrange(dimensions[i]): for j in xrange(dimensions[i]):
table.labels[index+j] = table.labels[index+j].replace(options.label[i],options.substitute[i]) table.labels[index+j] = table.labels[index+j].replace(m.group(2),options.substitute[i]) # replace name with substitute
if remarks != []: damask.util.croak(remarks) if remarks != []: damask.util.croak(remarks)
if errors != []: if errors != []:

82
processing/post/vtk_addData.py
View File

@ -1,82 +0,0 @@
#!/usr/bin/env python
# -*- coding: UTF-8 no BOM -*-
import os,glob,re
import damask
from optparse import OptionParser
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
# -----------------------------
def findTag(filename,tag):
with open(filename,'r') as myfile:
mypattern = re.compile(str(tag))
for line in myfile:
if mypattern.search(line): return True
return False
# -----------------------------
# MAIN FUNCTION STARTS HERE
# -----------------------------
# --- input parsing
parser = OptionParser(usage='%prog [options] directory', description = """
Add data from an ASCII table to a VTK geometry file.
""", version = scriptID)
parser.add_option('-s','--sub', action='store_true', dest='subdir', \
help='include files in subdirectories [%default]')
parser.set_defaults(subdir = False)
(options, dirname) = parser.parse_args()
# --- sanity checks
if dirname == []:
parser.print_help()
parser.error('no directory specified...')
else:
dirname = os.path.abspath(dirname[0]) # only use first argument
if not os.path.isdir(dirname):
parser.print_help()
parser.error('invalid directory "%s" specified...'%dirname)
# --- loop over "nodebased" and "ipbased" data files and
# copy data to corresponding geometry files
dataSetTag = {'nodebased':'POINT_DATA', 'ipbased':'CELL_DATA'}
for geomtype in ['nodebased','ipbased']:
for vtkfilename in glob.iglob(dirname+os.sep+'*'+geomtype+'*.vtk'):
if not os.path.dirname(vtkfilename) == dirname and not options.subdir: continue # include files in subdir?
datafilename = os.path.splitext(vtkfilename)[0] + '.txt'
if not os.path.exists(datafilename): continue # no corresponding datafile found
# --- read data from datafile
with open(datafilename,'r') as datafile: # open datafile in read mode
table = damask.ASCIItable(fileIn=datafile) # use ASCIItable class to read data file
table.head_read() # read ASCII header info
myData = []
while table.data_read(): # read line in datafile
myData.append(table.data)
myData = zip(*myData) # reorder data: first index now label, not node
# --- append data to vtkfile
with open(vtkfilename,'a') as vtkfile: # open vtkfile in append mode
print vtkfilename
if not findTag(vtkfilename,dataSetTag[geomtype]): # check if data set is already present...
vtkfile.write(dataSetTag[geomtype] + ' %i'%len(myData[0])) # ... if not, write keyword
for idx,label in enumerate(table.labels): # write data
vtkfile.write('\nSCALARS '+label+' double 1\nLOOKUP_TABLE default\n') # all scalar data
vtkfile.write('\n'.join(map(str,myData[idx])))

196
processing/post/vtk_addPointcloudData.py
View File

@ -1,8 +1,9 @@
#!/usr/bin/env python #!/usr/bin/env python
# -*- coding: UTF-8 no BOM -*- # -*- coding: UTF-8 no BOM -*-
import os,sys,vtk import os,vtk
import damask import damask
from collections import defaultdict
from optparse import OptionParser from optparse import OptionParser
scriptName = os.path.splitext(os.path.basename(__file__))[0] scriptName = os.path.splitext(os.path.basename(__file__))[0]
@ -17,100 +18,115 @@ Add scalar and RGB tuples from ASCIItable to existing VTK point cloud (.vtp).
""", version = scriptID) """, version = scriptID)
parser.add_option('-v', '--vtk', dest='vtk', \ parser.add_option( '--vtk',
dest = 'vtk',
type = 'string', metavar = 'string',
help = 'VTK file name') help = 'VTK file name')
parser.add_option( '--inplace',
dest = 'inplace',
action = 'store_true',
help = 'modify VTK file in-place')
parser.add_option('-r', '--render',
dest = 'render',
action = 'store_true',
help = 'open output in VTK render window')
parser.add_option('-s', '--scalar', dest='scalar', action='extend', \ parser.add_option('-s', '--scalar', dest='scalar', action='extend', \
help = 'scalar values') help = 'scalar values')
parser.add_option('-v', '--vector',
dest = 'vector',
action = 'extend', metavar = '<string LIST>',
help = 'vector value label(s)')
parser.add_option('-c', '--color', dest='color', action='extend', \ parser.add_option('-c', '--color', dest='color', action='extend', \
help = 'RGB color tuples') help = 'RGB color tuples')
parser.set_defaults(scalar = []) parser.set_defaults(scalar = [],
parser.set_defaults(color = []) vector = [],
color = [],
inplace = False,
render = False,
)
(options, filenames) = parser.parse_args() (options, filenames) = parser.parse_args()
datainfo = { # list of requested labels per datatype if not options.vtk: parser.error('No VTK file specified.')
'scalar': {'len':1, if not os.path.exists(options.vtk): parser.error('VTK file does not exist.')
'label':[]},
'color': {'len':3,
'label':[]},
}
if not os.path.exists(options.vtk):
parser.error('VTK file does not exist'); sys.exit()
if os.path.splitext(options.vtk)[1] == '.vtp':
reader = vtk.vtkXMLPolyDataReader() reader = vtk.vtkXMLPolyDataReader()
reader.SetFileName(options.vtk) reader.SetFileName(options.vtk)
reader.Update() reader.Update()
Npoints = reader.GetNumberOfPoints()
Ncells = reader.GetNumberOfCells()
Nvertices = reader.GetNumberOfVerts()
Polydata = reader.GetOutput() Polydata = reader.GetOutput()
elif os.path.splitext(options.vtk)[1] == '.vtk':
reader = vtk.vtkGenericDataObjectReader()
reader.SetFileName(options.vtk)
reader.Update()
Polydata = reader.GetPolyDataOutput()
else:
parser.error('Unsupported VTK file type extension.')
Npoints = Polydata.GetNumberOfPoints()
Ncells = Polydata.GetNumberOfCells()
Nvertices = Polydata.GetNumberOfVerts()
if Npoints != Ncells or Npoints != Nvertices: if Npoints != Ncells or Npoints != Nvertices:
parser.error('Number of points, cells, and vertices in VTK differ from each other'); sys.exit() parser.error('Number of points, cells, and vertices in VTK differ from each other.')
if options.scalar is not None: datainfo['scalar']['label'] += options.scalar
if options.color is not None: datainfo['color']['label'] += options.color
# ------------------------------------------ setup file handles --------------------------------------- damask.util.croak('{}: {} points, {} vertices, and {} cells...'.format(options.vtk,Npoints,Nvertices,Ncells))
# --- loop over input files -------------------------------------------------------------------------
if filenames == []: filenames = [None]
files = []
if filenames == []:
files.append({'name':'STDIN', 'input':sys.stdin, 'output':sys.stdout, 'croak':sys.stderr})
else:
for name in filenames: for name in filenames:
if os.path.exists(name): try: table = damask.ASCIItable(name = name,
files.append({'name':name, 'input':open(name), 'output':sys.stderr, 'croak':sys.stderr}) buffered = False,
readonly = True)
except: continue
damask.util.report(scriptName, name)
#--- loop over input files ------------------------------------------------------------------------ # --- interpret header ----------------------------------------------------------------------------
for file in files:
if file['name'] != 'STDIN': file['croak'].write('\033[1m'+scriptName+'\033[0m: '+file['name']+'\n')
else: file['croak'].write('\033[1m'+scriptName+'\033[0m\n')
table = damask.ASCIItable(file['input'],file['output'],False) # make unbuffered ASCII_table table.head_read()
table.head_read() # read ASCII header info
# --------------- figure out columns to process remarks = []
active = {} errors = []
column = {} VTKarray = {}
active = defaultdict(list)
array = {} for datatype,dimension,label in [['scalar',1,options.scalar],
['vector',3,options.vector],
['color',3,options.color],
]:
for i,dim in enumerate(table.label_dimension(label)):
me = label[i]
if dim == -1: remarks.append('{} "{}" not found...'.format(datatype,me))
elif dim > dimension: remarks.append('"{}" not of dimension {}...'.format(me,dimension))
else:
remarks.append('adding {} "{}"...'.format(datatype,me))
active[datatype].append(me)
for datatype,info in datainfo.items(): if datatype in ['scalar','vector']: VTKarray[me] = vtk.vtkDoubleArray()
for label in info['label']: elif datatype == 'color': VTKarray[me] = vtk.vtkUnsignedCharArray()
foundIt = False
for key in ['1_'+label,label]: VTKarray[me].SetNumberOfComponents(dimension)
if key in table.labels: VTKarray[me].SetName(label[i])
foundIt = True
if datatype not in active: active[datatype] = [] if remarks != []: damask.util.croak(remarks)
if datatype not in column: column[datatype] = {} if errors != []:
if datatype not in array: array[datatype] = {} damask.util.croak(errors)
active[datatype].append(label) table.close(dismiss = True)
column[datatype][label] = table.labels.index(key) # remember columns of requested data continue
if datatype == 'scalar':
array[datatype][label] = vtk.vtkDoubleArray()
array[datatype][label].SetNumberOfComponents(1)
array[datatype][label].SetName(label)
elif datatype == 'color':
array[datatype][label] = vtk.vtkUnsignedCharArray()
array[datatype][label].SetNumberOfComponents(3)
array[datatype][label].SetName(label)
if not foundIt:
file['croak'].write('column %s not found...\n'%label)
# ------------------------------------------ process data --------------------------------------- # ------------------------------------------ process data ---------------------------------------
while table.data_read(): # read next data line of ASCII table while table.data_read(): # read next data line of ASCII table
for datatype,labels in active.items(): # loop over scalar,color for datatype,labels in active.items(): # loop over scalar,color
for label in labels: # loop over all requested items for me in labels: # loop over all requested items
theData = table.data[column[datatype][label]:\ theData = [table.data[i] for i in table.label_indexrange(me)] # read strings
column[datatype][label]+datainfo[datatype]['len']] # read strings if datatype == 'color': VTKarray[me].InsertNextTuple3(*map(lambda x: int(255.*float(x)),theData))
if datatype == 'color': elif datatype == 'vector': VTKarray[me].InsertNextTuple3(*map(float,theData))
theData = map(lambda x: int(255.*float(x)),theData) elif datatype == 'scalar': VTKarray[me].InsertNextValue(float(theData[0]))
array[datatype][label].InsertNextTuple3(theData[0],theData[1],theData[2],)
elif datatype == 'scalar':
array[datatype][label].InsertNextValue(float(theData[0]))
table.input_close() # close input ASCII table table.input_close() # close input ASCII table
@ -118,24 +134,50 @@ for file in files:
for datatype,labels in active.items(): # loop over scalar,color for datatype,labels in active.items(): # loop over scalar,color
if datatype == 'color': if datatype == 'color':
Polydata.GetPointData().SetScalars(array[datatype][labels[0]]) Polydata.GetPointData().SetScalars(VTKarray[active['color'][0]])
Polydata.GetCellData().SetScalars(array[datatype][labels[0]]) Polydata.GetCellData().SetScalars(VTKarray[active['color'][0]])
for label in labels: # loop over all requested items for me in labels: # loop over all requested items
Polydata.GetPointData().AddArray(array[datatype][label]) Polydata.GetPointData().AddArray(VTKarray[me])
Polydata.GetCellData().AddArray(array[datatype][label]) Polydata.GetCellData().AddArray(VTKarray[me])
Polydata.Modified() Polydata.Modified()
if vtk.VTK_MAJOR_VERSION <= 5: if vtk.VTK_MAJOR_VERSION <= 5: Polydata.Update()
Polydata.Update()
# ------------------------------------------ output result --------------------------------------- # ------------------------------------------ output result ---------------------------------------
writer = vtk.vtkXMLPolyDataWriter() writer = vtk.vtkXMLPolyDataWriter()
writer.SetDataModeToBinary() writer.SetDataModeToBinary()
writer.SetCompressorTypeToZLib() writer.SetCompressorTypeToZLib()
writer.SetFileName(os.path.splitext(options.vtk)[0]+'_added.vtp') writer.SetFileName(os.path.splitext(options.vtk)[0]+('.vtp' if options.inplace else '_added.vtp'))
if vtk.VTK_MAJOR_VERSION <= 5: if vtk.VTK_MAJOR_VERSION <= 5: writer.SetInput(Polydata)
writer.SetInput(Polydata) else: writer.SetInputData(Polydata)
else:
writer.SetInputData(Polydata)
writer.Write() writer.Write()
# ------------------------------------------ render result ---------------------------------------
if options.render:
mapper = vtk.vtkDataSetMapper()
mapper.SetInputData(Polydata)
actor = vtk.vtkActor()
actor.SetMapper(mapper)
# Create the graphics structure. The renderer renders into the
# render window. The render window interactor captures mouse events
# and will perform appropriate camera or actor manipulation
# depending on the nature of the events.
ren = vtk.vtkRenderer()
renWin = vtk.vtkRenderWindow()
renWin.AddRenderer(ren)
ren.AddActor(actor)
ren.SetBackground(1, 1, 1)
renWin.SetSize(200, 200)
iren = vtk.vtkRenderWindowInteractor()
iren.SetRenderWindow(renWin)
iren.Initialize()
renWin.Render()
iren.Start()

32
processing/post/vtk_addRectilinearGridData.py
View File

@ -30,10 +30,6 @@ parser.add_option('-r', '--render',
dest = 'render', dest = 'render',
action = 'store_true', action = 'store_true',
help = 'open output in VTK render window') help = 'open output in VTK render window')
parser.add_option('-m', '--mode',
dest = 'mode',
type = 'choice', metavar = 'string', choices = ['cell', 'point'],
help = 'cell-centered or point-centered data')
parser.add_option('-s', '--scalar', parser.add_option('-s', '--scalar',
dest = 'scalar', dest = 'scalar',
action = 'extend', metavar = '<string LIST>', action = 'extend', metavar = '<string LIST>',
@ -56,7 +52,6 @@ parser.set_defaults(scalar = [],
(options, filenames) = parser.parse_args() (options, filenames) = parser.parse_args()
if not options.mode: parser.error('No data mode specified.')
if not options.vtk: parser.error('No VTK file specified.') if not options.vtk: parser.error('No VTK file specified.')
if not os.path.exists(options.vtk): parser.error('VTK file does not exist.') if not os.path.exists(options.vtk): parser.error('VTK file does not exist.')
@ -83,9 +78,9 @@ damask.util.croak('{}: {} points and {} cells...'.format(options.vtk,Npoints,Nce
if filenames == []: filenames = [None] if filenames == []: filenames = [None]
for name in filenames: for name in filenames:
try: try: table = damask.ASCIItable(name = name,
table = damask.ASCIItable(name = name, buffered = False,
buffered = False, readonly = True) readonly = True)
except: continue except: continue
damask.util.report(scriptName, name) damask.util.report(scriptName, name)
@ -124,8 +119,11 @@ for name in filenames:
# ------------------------------------------ process data --------------------------------------- # ------------------------------------------ process data ---------------------------------------
datacount = 0
while table.data_read(): # read next data line of ASCII table while table.data_read(): # read next data line of ASCII table
datacount += 1 # count data lines
for datatype,labels in active.items(): # loop over scalar,color for datatype,labels in active.items(): # loop over scalar,color
for me in labels: # loop over all requested items for me in labels: # loop over all requested items
theData = [table.data[i] for i in table.label_indexrange(me)] # read strings theData = [table.data[i] for i in table.label_indexrange(me)] # read strings
@ -133,15 +131,25 @@ for name in filenames:
elif datatype == 'vector': VTKarray[me].InsertNextTuple3(*map(float,theData)) elif datatype == 'vector': VTKarray[me].InsertNextTuple3(*map(float,theData))
elif datatype == 'scalar': VTKarray[me].InsertNextValue(float(theData[0])) elif datatype == 'scalar': VTKarray[me].InsertNextValue(float(theData[0]))
table.close() # close input ASCII table
# ------------------------------------------ add data --------------------------------------- # ------------------------------------------ add data ---------------------------------------
if datacount == Npoints: mode = 'point'
elif datacount == Ncells: mode = 'cell'
else:
damask.util.croak('Data count is incompatible with grid...')
continue
damask.util.croak('{} mode...'.format(mode))
for datatype,labels in active.items(): # loop over scalar,color for datatype,labels in active.items(): # loop over scalar,color
if datatype == 'color': if datatype == 'color':
if options.mode == 'cell': rGrid.GetCellData().SetScalars(VTKarray[active['color'][0]]) if mode == 'cell': rGrid.GetCellData().SetScalars(VTKarray[active['color'][0]])
elif options.mode == 'point': rGrid.GetPointData().SetScalars(VTKarray[active['color'][0]]) elif mode == 'point': rGrid.GetPointData().SetScalars(VTKarray[active['color'][0]])
for me in labels: # loop over all requested items for me in labels: # loop over all requested items
if options.mode == 'cell': rGrid.GetCellData().AddArray(VTKarray[me]) if mode == 'cell': rGrid.GetCellData().AddArray(VTKarray[me])
elif options.mode == 'point': rGrid.GetPointData().AddArray(VTKarray[me]) elif mode == 'point': rGrid.GetPointData().AddArray(VTKarray[me])
rGrid.Modified() rGrid.Modified()
if vtk.VTK_MAJOR_VERSION <= 5: rGrid.Update() if vtk.VTK_MAJOR_VERSION <= 5: rGrid.Update()

135
processing/post/vtk_addVoxelcloudData.py
View File

@ -1,135 +0,0 @@
#!/usr/bin/env python
# -*- coding: UTF-8 no BOM -*-
import os,sys,vtk
import damask
from collections import defaultdict
from optparse import OptionParser
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
# --------------------------------------------------------------------
# MAIN
# --------------------------------------------------------------------
parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """
Add scalar and RGB tuples from ASCIItable to existing VTK voxel cloud (.vtu/.vtk).
""", version = scriptID)
parser.add_option('-v', '--vtk', dest='vtk', \
help = 'VTK file name')
parser.add_option('-s', '--scalar', dest='scalar', action='extend', \
help = 'scalar values')
parser.add_option('-c', '--color', dest='color', action='extend', \
help = 'RGB color tuples')
parser.set_defaults(scalar = [],
color = [],
render = False,
)
(options, filenames) = parser.parse_args()
if options.vtk is None or not os.path.exists(options.vtk):
parser.error('VTK file does not exist')
if os.path.splitext(options.vtk)[1] == '.vtu':
reader = vtk.vtkXMLUnstructuredGridReader()
reader.SetFileName(options.vtk)
reader.Update()
uGrid = reader.GetOutput()
elif os.path.splitext(options.vtk)[1] == '.vtk':
reader = vtk.vtkGenericDataObjectReader()
reader.SetFileName(options.vtk)
reader.Update()
uGrid = reader.GetUnstructuredGridOutput()
else:
parser.error('unsupported VTK file type extension')
Npoints = uGrid.GetNumberOfPoints()
Ncells = uGrid.GetNumberOfCells()
sys.stderr.write('{}: {} points and {} cells...\n'.format(damask.util.emph(options.vtk),Npoints,Ncells))
# --- loop over input files -------------------------------------------------------------------------
if filenames == []: filenames = [None]
for name in filenames:
try:
table = damask.ASCIItable(name = name,
buffered = False, readonly = True)
except: continue
damask.util.croak(damask.util.emph(scriptName)+(': '+name if name else ''))
# --- interpret header ----------------------------------------------------------------------------
table.head_read()
remarks = []
errors = []
VTKarray = {}
active = defaultdict(list)
for datatype,dimension,label in [['scalar',1,options.scalar],
['color',3,options.color],
]:
for i,dim in enumerate(table.label_dimension(label)):
me = label[i]
if dim == -1: remarks.append('{} "{}" not found...'.format(datatype,me))
elif dim > dimension: remarks.append('"{}" not of dimension{}...'.format(me,dimension))
else:
damask.util.croak('adding {} {}'.format(datatype,me))
active[datatype].append(me)
if datatype in ['scalar']:
VTKarray[me] = vtk.vtkDoubleArray()
elif datatype == 'color':
VTKarray[me] = vtk.vtkUnsignedCharArray()
VTKarray[me].SetNumberOfComponents(dimension)
VTKarray[me].SetName(label[i])
if remarks != []: damask.util.croak(remarks)
if errors != []:
damask.util.croak(errors)
table.close(dismiss=True)
continue
# ------------------------------------------ process data ---------------------------------------
while table.data_read(): # read next data line of ASCII table
for datatype,labels in active.items(): # loop over scalar,color
for me in labels: # loop over all requested items
theData = [table.data[i] for i in table.label_indexrange(me)] # read strings
if datatype == 'color':
VTKarray[me].InsertNextTuple3(*map(lambda x: int(255.*float(x)),theData))
elif datatype == 'scalar':
VTKarray[me].InsertNextValue(float(theData[0]))
# ------------------------------------------ add data ---------------------------------------
for datatype,labels in active.items(): # loop over scalar,color
if datatype == 'color':
uGrid.GetCellData().SetScalars(VTKarray[active['color'][0]])
for label in labels: # loop over all requested items
uGrid.GetCellData().AddArray(VTKarray[me])
uGrid.Modified()
if vtk.VTK_MAJOR_VERSION <= 5:
uGrid.Update()
# ------------------------------------------ output result ---------------------------------------
writer = vtk.vtkXMLUnstructuredGridWriter()
writer.SetDataModeToBinary()
writer.SetCompressorTypeToZLib()
writer.SetFileName(os.path.splitext(options.vtk)[0]+'_added.vtu')
if vtk.VTK_MAJOR_VERSION <= 5:
writer.SetInput(uGrid)
else:
writer.SetInputData(uGrid)
writer.Write()

31
processing/post/vtk_pointcloud.py
View File

@ -18,12 +18,13 @@ Produce a VTK point cloud dataset based on coordinates given in an ASCIItable.
""", version = scriptID) """, version = scriptID)
parser.add_option('-d', '--deformed', parser.add_option('-p',
dest = 'deformed', '--pos', '--position',
dest = 'pos',
type = 'string', metavar = 'string', type = 'string', metavar = 'string',
help = 'deformed coordinate label [%default]') help = 'label of coordinates [%default]')
parser.set_defaults(deformed = 'ipdeformedcoord' parser.set_defaults(pos = 'pos',
) )
(options, filenames) = parser.parse_args() (options, filenames) = parser.parse_args()
@ -46,9 +47,9 @@ for name in filenames:
errors = [] errors = []
remarks = [] remarks = []
coordDim = table.label_dimension(options.deformed) coordDim = table.label_dimension(options.pos)
if not 3 >= coordDim >= 1: errors.append('coordinates "{}" need to have one, two, or three dimensions.'.format(options.deformed)) if not 3 >= coordDim >= 1: errors.append('coordinates "{}" need to have one, two, or three dimensions.'.format(options.pos))
elif coordDim < 3: remarks.append('appending {} dimensions to coordinates "{}"...'.format(3-coordDim,options.deformed)) elif coordDim < 3: remarks.append('appending {} dimensions to coordinates "{}"...'.format(3-coordDim,options.pos))
if remarks != []: damask.util.croak(remarks) if remarks != []: damask.util.croak(remarks)
if errors != []: if errors != []:
@ -58,7 +59,7 @@ for name in filenames:
# ------------------------------------------ process data --------------------------------------- # ------------------------------------------ process data ---------------------------------------
table.data_readArray(options.deformed) table.data_readArray(options.pos)
if len(table.data.shape) < 2: table.data.shape += (1,) # expand to 2D shape if len(table.data.shape) < 2: table.data.shape += (1,) # expand to 2D shape
if table.data.shape[1] < 3: if table.data.shape[1] < 3:
table.data = np.hstack((table.data, table.data = np.hstack((table.data,
@ -83,19 +84,21 @@ for name in filenames:
if name: if name:
writer = vtk.vtkXMLPolyDataWriter() writer = vtk.vtkXMLPolyDataWriter()
(directory,filename) = os.path.split(name)
writer.SetDataModeToBinary()
writer.SetCompressorTypeToZLib() writer.SetCompressorTypeToZLib()
writer.SetFileName(os.path.join(directory,os.path.splitext(filename)[0]\ writer.SetDataModeToBinary()
+'.'+writer.GetDefaultFileExtension())) writer.SetFileName(os.path.join(os.path.split(name)[0],
os.path.splitext(os.path.split(name)[1])[0] +
'.' + writer.GetDefaultFileExtension()))
else: else:
writer = vtk.vtkDataSetWriter() writer = vtk.vtkDataSetWriter()
writer.WriteToOutputStringOn()
writer.SetHeader('# powered by '+scriptID) writer.SetHeader('# powered by '+scriptID)
writer.WriteToOutputStringOn()
if vtk.VTK_MAJOR_VERSION <= 5: writer.SetInput(Polydata) if vtk.VTK_MAJOR_VERSION <= 5: writer.SetInput(Polydata)
else: writer.SetInputData(Polydata) else: writer.SetInputData(Polydata)
writer.Write() writer.Write()
if name is None: sys.stdout.write(writer.GetOutputString()[0:writer.GetOutputStringLength()])
if name is None: sys.stdout.write(writer.GetOutputString()[:writer.GetOutputStringLength()]) # limiting of outputString is fix for vtk <7.0
table.close() table.close()

80
processing/post/vtk_rectilinearGrid.py
View File

@ -6,7 +6,6 @@ import numpy as np
import damask import damask
from optparse import OptionParser from optparse import OptionParser
scriptName = os.path.splitext(os.path.basename(__file__))[0] scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version]) scriptID = ' '.join([scriptName,damask.version])
@ -15,20 +14,29 @@ scriptID = ' '.join([scriptName,damask.version])
# -------------------------------------------------------------------- # --------------------------------------------------------------------
parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """ parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """
Create regular voxel grid from points in an ASCIItable. Create regular voxel grid from points in an ASCIItable (or geom file).
""", version = scriptID) """, version = scriptID)
parser.add_option('-m', '--mode', parser.add_option('-m',
'--mode',
dest = 'mode', dest = 'mode',
type = 'choice', choices = ['cell','point'], type = 'choice', choices = ['cell','point'],
help = 'cell-centered or point-centered coordinates') help = 'cell-centered or point-centered coordinates')
parser.add_option('-c', '--coordinates', parser.add_option('-p',
dest = 'position', '--pos', '--position',
dest = 'pos',
type = 'string', metavar = 'string', type = 'string', metavar = 'string',
help = 'coordinate label [%default]') help = 'label of coordinates [%default]')
parser.set_defaults(position ='ipinitialcoord', parser.add_option('-g',
mode ='cell' '--geom',
dest = 'geom',
action = 'store_true',
help = 'geom input format')
parser.set_defaults(mode = 'cell',
pos = 'pos',
geom = False,
) )
(options, filenames) = parser.parse_args() (options, filenames) = parser.parse_args()
@ -38,9 +46,12 @@ parser.set_defaults(position ='ipinitialcoord',
if filenames == []: filenames = [None] if filenames == []: filenames = [None]
for name in filenames: for name in filenames:
try: isGeom = options.geom or (name is not None and name.endswith('.geom'))
table = damask.ASCIItable(name = name, try: table = damask.ASCIItable(name = name,
buffered = False, readonly = True) buffered = False,
labeled = not isGeom,
readonly = True,
)
except: continue except: continue
damask.util.report(scriptName,name) damask.util.report(scriptName,name)
@ -48,10 +59,13 @@ for name in filenames:
table.head_read() table.head_read()
remarks = []
errors = [] errors = []
if table.label_dimension(options.position) != 3: coordDim = 3 if isGeom else table.label_dimension(options.pos)
errors.append('coordinates {} are not a vector.'.format(options.position)) if not 3 >= coordDim >= 1: errors.append('coordinates "{}" need to have one, two, or three dimensions.'.format(options.pos))
elif coordDim < 3: remarks.append('appending {} dimensions to coordinates "{}"...'.format(3-coordDim,options.pos))
if remarks != []: damask.util.croak(remarks)
if errors != []: if errors != []:
damask.util.croak(errors) damask.util.croak(errors)
table.close(dismiss=True) table.close(dismiss=True)
@ -59,17 +73,33 @@ for name in filenames:
# --------------- figure out size and grid --------------------------------------------------------- # --------------- figure out size and grid ---------------------------------------------------------
table.data_readArray(options.position) if isGeom:
info,extra_header = table.head_getGeom()
coords = [np.linspace(info['origin'][i],
info['origin'][i]+info['size'][i],
num = info['grid'][i]+1,
endpoint = True,
) for i in xrange(3)]
else:
table.data_readArray(options.pos)
if len(table.data.shape) < 2: table.data.shape += (1,) # expand to 2D shape
if table.data.shape[1] < 3:
table.data = np.hstack((table.data,
np.zeros((table.data.shape[0],
3-table.data.shape[1]),dtype='f'))) # fill coords up to 3D with zeros
coords = [np.unique(table.data[:,i]) for i in xrange(3)] coords = [np.unique(table.data[:,i]) for i in xrange(3)]
if options.mode == 'cell': if options.mode == 'cell':
coords = [0.5 * np.array([3.0 * coords[i][0] - coords[i][0 + len(coords[i]) > 1]] + \ coords = [0.5 * np.array([3.0 * coords[i][0] - coords[i][0 + len(coords[i]) > 1]] + \
[coords[i][j-1] + coords[i][j] for j in xrange(1,len(coords[i]))] + \ [coords[i][j-1] + coords[i][j] for j in xrange(1,len(coords[i]))] + \
[3.0 * coords[i][-1] - coords[i][-1 - (len(coords[i]) > 1)]]) for i in xrange(3)] [3.0 * coords[i][-1] - coords[i][-1 - (len(coords[i]) > 1)]]) for i in xrange(3)]
grid = np.array(map(len,coords),'i')
N = grid.prod() if options.mode == 'point' else (grid-1).prod()
if N != len(table.data): errors.append('data count {} does not match grid {}x{}x{}.'.format(N,*(grid - options.mode == 'cell') )) grid = np.array(map(len,coords),'i')
N = grid.prod() if options.mode == 'point' or isGeom else (grid-1).prod()
if not isGeom and N != len(table.data):
errors.append('data count {} does not match grid {}x{}x{}.'.format(N,*(grid - (options.mode == 'cell')) ))
if errors != []: if errors != []:
damask.util.croak(errors) damask.util.croak(errors)
table.close(dismiss = True) table.close(dismiss = True)
@ -97,20 +127,22 @@ for name in filenames:
if name: if name:
writer = vtk.vtkXMLRectilinearGridWriter() writer = vtk.vtkXMLRectilinearGridWriter()
(directory,filename) = os.path.split(name)
writer.SetDataModeToBinary()
writer.SetCompressorTypeToZLib() writer.SetCompressorTypeToZLib()
writer.SetFileName(os.path.join(directory,os.path.splitext(filename)[0] \ writer.SetDataModeToBinary()
+'_{}({})'.format(options.position, options.mode) \ writer.SetFileName(os.path.join(os.path.split(name)[0],
+'.'+writer.GetDefaultFileExtension())) os.path.splitext(os.path.split(name)[1])[0] +
('' if isGeom else '_{}({})'.format(options.pos, options.mode)) +
'.' + writer.GetDefaultFileExtension()))
else: else:
writer = vtk.vtkDataSetWriter() writer = vtk.vtkDataSetWriter()
writer.WriteToOutputStringOn()
writer.SetHeader('# powered by '+scriptID) writer.SetHeader('# powered by '+scriptID)
writer.WriteToOutputStringOn()
if vtk.VTK_MAJOR_VERSION <= 5: writer.SetInput(rGrid) if vtk.VTK_MAJOR_VERSION <= 5: writer.SetInput(rGrid)
else: writer.SetInputData(rGrid) else: writer.SetInputData(rGrid)
writer.Write() writer.Write()
if name is None: sys.stdout.write(writer.GetOutputString()[0:writer.GetOutputStringLength()])
if name is None: sys.stdout.write(writer.GetOutputString()[:writer.GetOutputStringLength()]) # limiting of outputString is fix for vtk <7.0
table.close() table.close()

122
processing/post/vtk_scalars2vectors.py
View File

@ -1,122 +0,0 @@
#!/usr/bin/env python
# -*- coding: UTF-8 no BOM -*-
import os,sys,shutil
import damask
from optparse import OptionParser
import vtk
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
# --------------------------------------------------------------------
# MAIN
# --------------------------------------------------------------------
parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """
""", version = scriptID)
parser.add_option('-v','--vector', nargs=3, dest='vector', \
help='suffices indicating vector components [%default]')
parser.add_option('-s','--separator', dest='separator', \
help='separator between label and suffix [%default]')
parser.set_defaults(vector = ['x','y','z'])
parser.set_defaults(separator = '.')
(options, filenames) = parser.parse_args()
# --- sanity checks
if filenames == []:
parser.print_help()
parser.error('no file specified...')
for filename in filenames:
if not os.path.isfile(filename):
parser.print_help()
parser.error('invalid file "%s" specified...'%filename)
# --- ITERATE OVER FILES AND PROCESS THEM
for filename in filenames:
sys.stdout.write('read file "%s" ...'%filename)
sys.stdout.flush()
suffix = os.path.splitext(filename)[1]
if suffix == '.vtk':
reader = vtk.vtkUnstructuredGridReader()
reader.ReadAllScalarsOn()
reader.ReadAllVectorsOn()
reader.ReadAllTensorsOn()
elif suffix == '.vtu':
reader = vtk.vtkXMLUnstructuredGridReader()
else:
parser.error('filetype "%s" not supported'%suffix)
reader.SetFileName(filename)
reader.Update()
uGrid = reader.GetOutput()
sys.stdout.write(' done\n')
sys.stdout.flush()
# Read the scalar data
scalarData = {}
scalarsToBeRemoved = []
Nscalars = uGrid.GetCellData().GetNumberOfArrays()
for i in range(Nscalars):
sys.stdout.write("\rread scalar data %d%%" %(100*i/Nscalars))
sys.stdout.flush()
scalarName = uGrid.GetCellData().GetArrayName(i)
if scalarName.split(options.separator)[-1] in options.vector:
label,suffix = scalarName.split(options.separator)
if label not in scalarData:
scalarData[label] = [[],[],[]]
uGrid.GetCellData().SetActiveScalars(scalarName)
scalarData[label][options.vector.index(suffix)] = uGrid.GetCellData().GetScalars(scalarName)
scalarsToBeRemoved.append(scalarName)
for scalarName in scalarsToBeRemoved:
uGrid.GetCellData().RemoveArray(scalarName)
sys.stdout.write('\rread scalar data done\n')
sys.stdout.flush()
# Convert the scalar data to vector data
NscalarData = len(scalarData)
for n,label in enumerate(scalarData):
sys.stdout.write("\rconvert to vector data %d%%" %(100*n/NscalarData))
sys.stdout.flush()
Nvalues = scalarData[label][0].GetNumberOfTuples()
vectorData = vtk.vtkDoubleArray()
vectorData.SetName(label)
vectorData.SetNumberOfComponents(3) # set this before NumberOfTuples !!!
vectorData.SetNumberOfTuples(Nvalues)
for i in range(Nvalues):
for j in range(3):
vectorData.SetComponent(i,j,scalarData[label][j].GetValue(i))
uGrid.GetCellData().AddArray(vectorData)
sys.stdout.write('\rconvert to vector data done\n')
# Write to new vtk file
outfilename = os.path.splitext(filename)[0]+'.vtu'
sys.stdout.write('write to file "%s" ...'%outfilename)
sys.stdout.flush()
writer = vtk.vtkXMLUnstructuredGridWriter()
writer.SetFileName(outfilename+'_tmp')
writer.SetDataModeToAscii()
writer.SetInput(uGrid)
writer.Write()
sys.stdout.write(' done\n')
sys.stdout.flush()
shutil.move(outfilename+'_tmp',outfilename)
# --------------------------- DONE --------------------------------

118
processing/post/vtk_voxelcloud.py
View File

@ -1,118 +0,0 @@
#!/usr/bin/env python
# -*- coding: UTF-8 no BOM -*-
import os,sys,vtk
import numpy as np
from optparse import OptionParser
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
# --------------------------------------------------------------------
# MAIN
# --------------------------------------------------------------------
parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """
Create hexahedral voxels around points in an ASCIItable.
""", version = scriptID)
parser.add_option('-d', '--deformed',
dest = 'deformed',
type = 'string', metavar = 'string',
help = 'deformed coordinate label [%default]')
parser.add_option('-c','--coordinates',
dest = 'coords',
type = 'string', metavar='string',
help = 'undeformed coordinates label [%default]')
parser.set_defaults(deformed = 'ipdeformedcoord',
coords = 'ipinitialcoord',
)
(options, filenames) = parser.parse_args()
# --- loop over input files -------------------------------------------------------------------------
if filenames == []: filenames = [None]
for name in filenames:
try:
table = damask.ASCIItable(name = name,
buffered = False, readonly = True)
except: continue
damask.util.report(scriptName,name)
# --------------- interprete header -----------------------------------------------------------------
table.head_read()
errors=[]
if table.label_dimension(options.deformed) != 3:
errors.append('columns "{}" have dimension {}'.format(options.deformed,table.label_dimension(options.deformed)))
if table.label_dimension(options.coords) != 3:
errors.append('coordinates {} are not a vector.'.format(options.coords))
table.data_readArray([options.coords,options.deformed])
# --------------- figure out size and grid ---------------------------------------------------------
coords = [{},{},{}]
for i in xrange(len(table.data)):
for j in xrange(3):
coords[j][str(table.data[i,table.label_index(options.coords)+j])] = True
grid = np.array(map(len,coords),'i')
size = grid/np.maximum(np.ones(3,'d'),grid-1.0)* \
np.array([max(map(float,coords[0].keys()))-min(map(float,coords[0].keys())),\
max(map(float,coords[1].keys()))-min(map(float,coords[1].keys())),\
max(map(float,coords[2].keys()))-min(map(float,coords[2].keys())),\
],'d') # size from bounding box, corrected for cell-centeredness
size = np.where(grid > 1, size, min(size[grid > 1]/grid[grid > 1])) # spacing for grid==1 set to smallest among other spacings
# ------------------------------------------ process data ---------------------------------------
hexPoints = np.array([[-1,-1,-1],
[ 1,-1,-1],
[ 1, 1,-1],
[-1, 1,-1],
[-1,-1, 1],
[ 1,-1, 1],
[ 1, 1, 1],
[-1, 1, 1],
])
halfDelta = 0.5*size/grid
Points = vtk.vtkPoints()
Hex = vtk.vtkHexahedron()
uGrid = vtk.vtkUnstructuredGrid()
for p in table.data:
for i,h in enumerate(hexPoints):
pointID = Points.InsertNextPoint(p[table.label_index(options.deformed):table.label_index(options.deformed)+3]+h*halfDelta)
Hex.GetPointIds().SetId(i,pointID)
uGrid.InsertNextCell(Hex.GetCellType(), Hex.GetPointIds())
uGrid.SetPoints(Points)
# ------------------------------------------ output result ---------------------------------------
if name:
writer = vtk.vtkXMLUnstructuredGridWriter()
(directory,filename) = os.path.split(name)
writer.SetDataModeToBinary()
writer.SetCompressorTypeToZLib()
writer.SetFileName(os.path.join(directory,os.path.splitext(filename)[0]\
+'.'+writer.GetDefaultFileExtension()))
else:
writer = vtk.vtkDataSetWriter()
writer.WriteToOutputStringOn()
writer.SetHeader('# powered by '+scriptID)
if vtk.VTK_MAJOR_VERSION <= 5: writer.SetInput(uGrid)
else: writer.SetInputData(uGrid)
writer.Write()
if name is None: sys.stdout.write(writer.GetOutputString()[0:writer.GetOutputStringLength()])
table.close() # close input ASCII table

42
processing/pre/geom_canvas.py
View File

@ -15,25 +15,29 @@ scriptID = ' '.join([scriptName,damask.version])
parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """ parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """
Changes the (three-dimensional) canvas of a spectral geometry description. Changes the (three-dimensional) canvas of a spectral geometry description.
Grid can be given as absolute or relative values, e.g. 16 16 16 or 2x 0.5x 32.
""", version = scriptID) """, version = scriptID)
parser.add_option('-g', '--grid', parser.add_option('-g',
'--grid',
dest = 'grid', dest = 'grid',
type = 'string', nargs = 3, metavar = ' '.join(['string']*3), type = 'string', nargs = 3, metavar = ' '.join(['string']*3),
help = 'a,b,c grid of hexahedral box [auto]') help = 'a,b,c grid of hexahedral box. [auto]')
parser.add_option('-o', '--offset', parser.add_option('-o',
'--offset',
dest = 'offset', dest = 'offset',
type = 'int', nargs = 3, metavar = ' '.join(['int']*3), type = 'int', nargs = 3, metavar = ' '.join(['int']*3),
help = 'a,b,c offset from old to new origin of grid [%default]') help = 'a,b,c offset from old to new origin of grid [%default]')
parser.add_option('-f', '--fill', parser.add_option('-f',
'--fill',
dest = 'fill', dest = 'fill',
type = 'float', metavar = 'float', type = 'float', metavar = 'float',
help = '(background) canvas grain index. "0" selects maximum microstructure index + 1 [%default]') help = '(background) canvas grain index. "0" selects maximum microstructure index + 1 [%default]')
parser.add_option('--float', parser.add_option('--float',
dest = 'real', dest = 'real',
action = 'store_true', action = 'store_true',
help = 'input data is float [%default]') help = 'use float input')
parser.set_defaults(grid = ['0','0','0'], parser.set_defaults(grid = ['0','0','0'],
offset = (0,0,0), offset = (0,0,0),
@ -61,13 +65,7 @@ for name in filenames:
table.head_read() table.head_read()
info,extra_header = table.head_getGeom() info,extra_header = table.head_getGeom()
damask.util.report_geom(info)
damask.util.croak(['grid a b c: %s'%(' x '.join(map(str,info['grid']))),
'size x y z: %s'%(' x '.join(map(str,info['size']))),
'origin x y z: %s'%(' : '.join(map(str,info['origin']))),
'homogenization: %i'%info['homogenization'],
'microstructures: %i'%info['microstructures'],
])
errors = [] errors = []
if np.any(info['grid'] < 1): errors.append('invalid grid a b c.') if np.any(info['grid'] < 1): errors.append('invalid grid a b c.')
@ -105,7 +103,7 @@ for name in filenames:
translate_y = [i - options.offset[1] for i in yindex] translate_y = [i - options.offset[1] for i in yindex]
translate_z = [i - options.offset[2] for i in zindex] translate_z = [i - options.offset[2] for i in zindex]
if 0 in map(len,[xindex,yindex,zindex,translate_x,translate_y,translate_z]): if 0 in map(len,[xindex,yindex,zindex,translate_x,translate_y,translate_z]):
damask.util.croak('Invaldid grid-offset comination') damask.util.croak('invaldid grid-offset combination.')
table.close(dismiss = True) table.close(dismiss = True)
continue continue
microstructure_cropped[min(translate_x):(max(translate_x)+1),\ microstructure_cropped[min(translate_x):(max(translate_x)+1),\
@ -125,13 +123,13 @@ for name in filenames:
errors = [] errors = []
if (any(newInfo['grid'] != info['grid'])): if (any(newInfo['grid'] != info['grid'])):
remarks.append('--> grid a b c: %s'%(' x '.join(map(str,newInfo['grid'])))) remarks.append('--> grid a b c: {}'.format(' x '.join(map(str,newInfo['grid']))))
if (any(newInfo['size'] != info['size'])): if (any(newInfo['size'] != info['size'])):
remarks.append('--> size x y z: %s'%(' x '.join(map(str,newInfo['size'])))) remarks.append('--> size x y z: {}'.format(' x '.join(map(str,newInfo['size']))))
if (any(newInfo['origin'] != info['origin'])): if (any(newInfo['origin'] != info['origin'])):
remarks.append('--> origin x y z: %s'%(' : '.join(map(str,newInfo['origin'])))) remarks.append('--> origin x y z: {}'.format(' : '.join(map(str,newInfo['origin']))))
if ( newInfo['microstructures'] != info['microstructures']): if ( newInfo['microstructures'] != info['microstructures']):
remarks.append('--> microstructures: %i'%newInfo['microstructures']) remarks.append('--> microstructures: {}'.format(newInfo['microstructures']))
if np.any(newInfo['grid'] < 1): errors.append('invalid new grid a b c.') if np.any(newInfo['grid'] < 1): errors.append('invalid new grid a b c.')
if np.any(newInfo['size'] <= 0.0): errors.append('invalid new size x y z.') if np.any(newInfo['size'] <= 0.0): errors.append('invalid new size x y z.')
@ -147,11 +145,11 @@ for name in filenames:
table.info_clear() table.info_clear()
table.info_append([ table.info_append([
scriptID + ' ' + ' '.join(sys.argv[1:]), scriptID + ' ' + ' '.join(sys.argv[1:]),
"grid\ta {grid[0]}\tb {grid[1]}\tc {grid[2]}".format(grid=newInfo['grid']), "grid\ta {}\tb {}\tc {}".format(*newInfo['grid']),
"size\tx {size[0]}\ty {size[1]}\tz {size[2]}".format(size=newInfo['size']), "size\tx {}\ty {}\tz {}".format(*newInfo['size']),
"origin\tx {origin[0]}\ty {origin[1]}\tz {origin[2]}".format(origin=newInfo['origin']), "origin\tx {}\ty {}\tz {}".format(*newInfo['origin']),
"homogenization\t{homog}".format(homog=info['homogenization']), "homogenization\t{}".format(info['homogenization']),
"microstructures\t{microstructures}".format(microstructures=newInfo['microstructures']), "microstructures\t{}".format(newInfo['microstructures']),
extra_header extra_header
]) ])
table.labels_clear() table.labels_clear()

108
processing/pre/geom_check.py
View File

@ -1,108 +0,0 @@
#!/usr/bin/env python
# -*- coding: UTF-8 no BOM -*-
import os,sys,vtk
import numpy as np
from optparse import OptionParser
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
#--------------------------------------------------------------------------------------------------
# MAIN
#--------------------------------------------------------------------------------------------------
parser = OptionParser(option_class=damask.extendableOption, usage='%prog [file[s]]', description = """
Produce VTK rectilinear mesh of structure data from geom description
""", version = scriptID)
parser.add_option('-m','--nodata',
dest = 'data',
action = 'store_false',
help = 'generate mesh without microstructure index data')
parser.set_defaults(data = True,
)
(options, filenames) = parser.parse_args()
# --- loop over input files -------------------------------------------------------------------------
if filenames == []: filenames = [None]
for name in filenames:
try:
table = damask.ASCIItable(name = name,
buffered = False, labeled = False, readonly = True)
except: continue
damask.util.report(scriptName,name)
# --- interpret header ----------------------------------------------------------------------------
table.head_read()
info,extra_header = table.head_getGeom()
damask.util.croak(['grid a b c: %s'%(' x '.join(map(str,info['grid']))),
'size x y z: %s'%(' x '.join(map(str,info['size']))),
'origin x y z: %s'%(' : '.join(map(str,info['origin']))),
'homogenization: %i'%info['homogenization'],
'microstructures: %i'%info['microstructures'],
])
errors = []
if np.any(info['grid'] < 1): errors.append('invalid grid a b c.')
if np.any(info['size'] <= 0.0): errors.append('invalid size x y z.')
#--- read microstructure information --------------------------------------------------------------
if options.data:
microstructure,ok = table.microstructure_read(info['grid'],strict = True) # read microstructure
if ok:
structure = vtk.vtkIntArray()
structure.SetName('Microstructures')
for idx in microstructure: structure.InsertNextValue(idx)
else: errors.append('mismatch between data and grid dimension.')
if errors != []:
damask.util.croak(errors)
table.close(dismiss = True)
continue
# --- generate VTK rectilinear grid ---------------------------------------------------------------
grid = vtk.vtkRectilinearGrid()
grid.SetDimensions([x+1 for x in info['grid']])
for i in xrange(3):
temp = vtk.vtkDoubleArray()
temp.SetNumberOfTuples(info['grid'][i]+1)
for j in xrange(info['grid'][i]+1):
temp.InsertTuple1(j,j*info['size'][i]/info['grid'][i]+info['origin'][i])
if i == 0: grid.SetXCoordinates(temp)
elif i == 1: grid.SetYCoordinates(temp)
elif i == 2: grid.SetZCoordinates(temp)
if options.data: grid.GetCellData().AddArray(structure)
# --- write data -----------------------------------------------------------------------------------
if name:
writer = vtk.vtkXMLRectilinearGridWriter()
(directory,filename) = os.path.split(name)
writer.SetDataModeToBinary()
writer.SetCompressorTypeToZLib()
writer.SetFileName(os.path.join(directory,os.path.splitext(filename)[0]+'.'+writer.GetDefaultFileExtension()))
else:
writer = vtk.vtkDataSetWriter()
writer.WriteToOutputStringOn()
writer.SetHeader('# powered by '+scriptID)
if vtk.VTK_MAJOR_VERSION <= 5: writer.SetInput(grid)
else: writer.SetInputData(grid)
writer.Write()
if name is None: sys.stdout.write(writer.GetOutputString()[0:writer.GetOutputStringLength()])
table.close()

15
processing/pre/geom_check.sh Executable file
View File

@ -0,0 +1,15 @@
#!/bin/bash
for geom in "$@"
do
vtk_rectilinearGrid \
--geom $geom
geom_toTable \
< $geom \
| \
vtk_addRectilinearGridData \
--scalar microstructure \
--inplace \
--vtk ${geom%.*}.vtr
done

8
processing/pre/geom_clean.py
View File

@ -6,15 +6,12 @@ import numpy as np
import damask import damask
from scipy import ndimage from scipy import ndimage
from optparse import OptionParser from optparse import OptionParser
from collections import defaultdict
scriptName = os.path.splitext(os.path.basename(__file__))[0] scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version]) scriptID = ' '.join([scriptName,damask.version])
def mostFrequent(arr): def mostFrequent(arr):
d = defaultdict(int) return np.argmax(np.bincount(arr))
for i in arr: d[i] += 1
return sorted(d.iteritems(), key=lambda x: x[1], reverse=True)[0][0] # return value of most frequent microstructure
#-------------------------------------------------------------------------------------------------- #--------------------------------------------------------------------------------------------------
@ -43,8 +40,7 @@ parser.set_defaults(stencil = 3,
if filenames == []: filenames = [None] if filenames == []: filenames = [None]
for name in filenames: for name in filenames:
try: try: table = damask.ASCIItable(name = name,
table = damask.ASCIItable(name = name,
buffered = False, buffered = False,
labeled = False) labeled = False)
except: continue except: continue

110
processing/pre/geom_fromImage.py
View File

@ -1,110 +0,0 @@
#!/usr/bin/env python
# -*- coding: UTF-8 no BOM -*-
import os,sys,math
import numpy as np
from optparse import OptionParser
from PIL import Image
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
#--------------------------------------------------------------------------------------------------
# MAIN
#--------------------------------------------------------------------------------------------------
parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """
Generate geometry description from (multilayer) images.
Microstructure index is based on gray scale value (1..256).
""", version = scriptID)
parser.add_option('--homogenization',
dest = 'homogenization',
type = 'int', metavar = 'int',
help = 'homogenization index [%default]')
parser.set_defaults(homogenization = 1,
)
(options,filenames) = parser.parse_args()
# --- loop over input files -------------------------------------------------------------------------
if filenames == []: filenames = [None]
for name in filenames:
try:
table = damask.ASCIItable(name = name,
outname = os.path.splitext(name)[0]+'.geom' if name else name,
buffered = False, labeled = False)
except: continue
table.croak('\033[1m'+scriptName+'\033[0m'+(': '+name if name else ''))
# --- read image ------------------------------------------------------------------------------------
img = Image.open(name).convert(mode = 'L') # open and convert to grayscale 8bit
slice = 0
while True:
try:
img.seek(slice) # advance to slice
layer = np.expand_dims(1+np.array(img,dtype = 'uint16'),axis = 0) # read image layer
microstructure = layer if slice == 0 else np.vstack((microstructure,layer)) # noqa
slice += 1 # advance to next slice
except EOFError:
break
# http://docs.scipy.org/doc/scipy/reference/ndimage.html
# http://scipy-lectures.github.io/advanced/image_processing/
info = {
'grid': np.array(microstructure.shape,'i')[::-1],
'size': np.array(microstructure.shape,'d')[::-1],
'origin': np.zeros(3,'d'),
'microstructures': len(np.unique(microstructure)),
'homogenization': options.homogenization,
}
# --- report ---------------------------------------------------------------------------------------
table.croak(['grid a b c: %s'%(' x '.join(map(str,info['grid']))),
'size x y z: %s'%(' x '.join(map(str,info['size']))),
'origin x y z: %s'%(' : '.join(map(str,info['origin']))),
'homogenization: %i'%info['homogenization'],
'microstructures: %i'%info['microstructures'],
])
errors = []
if np.any(info['grid'] < 1): errors.append('invalid grid a b c.')
if np.any(info['size'] <= 0.0): errors.append('invalid size x y z.')
if errors != []:
table.croak(errors)
table.close(dismiss = True)
continue
# --- write header ---------------------------------------------------------------------------------
table.info_clear()
table.info_append([
scriptID + ' ' + ' '.join(sys.argv[1:]),
"grid\ta {grid[0]}\tb {grid[1]}\tc {grid[2]}".format(grid=info['grid']),
"size\tx {size[0]}\ty {size[1]}\tz {size[2]}".format(size=info['size']),
"origin\tx {origin[0]}\ty {origin[1]}\tz {origin[2]}".format(origin=info['origin']),
"homogenization\t{homog}".format(homog=info['homogenization']),
"microstructures\t{microstructures}".format(microstructures=info['microstructures']),
])
table.labels_clear()
table.head_write()
table.output_flush()
# --- write microstructure information ------------------------------------------------------------
formatwidth = int(math.floor(math.log10(microstructure.max())+1))
table.data = microstructure.reshape((info['grid'][1]*info['grid'][2],info['grid'][0]),order='C')
table.data_writeArray('%%%ii'%(formatwidth),delimiter = ' ')
# --- output finalization --------------------------------------------------------------------------
table.close() # close ASCII table

211
processing/pre/geom_fromVoronoiTessellation.py
View File

@ -4,7 +4,7 @@
import os,sys,math import os,sys,math
import numpy as np import numpy as np
import multiprocessing import multiprocessing
from optparse import OptionParser from optparse import OptionParser,OptionGroup
from scipy import spatial from scipy import spatial
import damask import damask
@ -109,97 +109,111 @@ Generate geometry description and material configuration by standard Voronoi tes
""", version = scriptID) """, version = scriptID)
parser.add_option('-g', '--grid',
dest = 'grid', group = OptionGroup(parser, "Tessellation","")
type = 'int', nargs = 3, metavar = ' '.join(['int']*3),
help = 'a,b,c grid of hexahedral box [auto]') group.add_option('-l',
parser.add_option('-s', '--size', '--laguerre',
dest = 'size',
type = 'float', nargs = 3, metavar=' '.join(['float']*3),
help = 'x,y,z size of hexahedral box [auto]')
parser.add_option('-o', '--origin',
dest = 'origin',
type = 'float', nargs = 3, metavar=' '.join(['float']*3),
help = 'offset from old to new origin of grid')
parser.add_option('-p', '--position',
dest = 'position',
type = 'string', metavar = 'string',
help = 'column label for seed positions [%default]')
parser.add_option('-w', '--weight',
dest = 'weight',
type = 'string', metavar = 'string',
help = 'column label for seed weights [%default]')
parser.add_option('-m', '--microstructure',
dest = 'microstructure',
type = 'string', metavar = 'string',
help = 'column label for seed microstructures [%default]')
parser.add_option('-e', '--eulers',
dest = 'eulers',
type = 'string', metavar = 'string',
help = 'column label for seed Euler angles [%default]')
parser.add_option('--axes',
dest = 'axes',
type = 'string', nargs = 3, metavar = ' '.join(['string']*3),
help = 'orientation coordinate frame in terms of position coordinate frame')
parser.add_option('--homogenization',
dest = 'homogenization',
type = 'int', metavar = 'int',
help = 'homogenization index to be used [%default]')
parser.add_option('--crystallite',
dest = 'crystallite',
type = 'int', metavar = 'int',
help = 'crystallite index to be used [%default]')
parser.add_option('--phase',
dest = 'phase',
type = 'int', metavar = 'int',
help = 'phase index to be used [%default]')
parser.add_option('-r', '--rnd',
dest = 'randomSeed',
type = 'int', metavar='int',
help = 'seed of random number generator for second phase distribution [%default]')
parser.add_option('--secondphase',
dest = 'secondphase',
type = 'float', metavar= 'float',
help = 'volume fraction of randomly distribute second phase [%default]')
parser.add_option('-l', '--laguerre',
dest = 'laguerre', dest = 'laguerre',
action = 'store_true', action = 'store_true',
help = 'use Laguerre (weighted Voronoi) tessellation') help = 'use Laguerre (weighted Voronoi) tessellation')
parser.add_option('--cpus', group.add_option('--cpus',
dest = 'cpus', dest = 'cpus',
type = 'int', metavar = 'int', type = 'int', metavar = 'int',
help = 'number of parallel processes to use for Laguerre tessellation [%default]') help = 'number of parallel processes to use for Laguerre tessellation [%default]')
parser.add_option('--nonperiodic', group.add_option('--nonperiodic',
dest = 'nonperiodic', dest = 'nonperiodic',
action = 'store_true', action = 'store_true',
help = 'use nonperiodic tessellation') help = 'nonperiodic tessellation')
parser.set_defaults(position = 'pos', parser.add_option_group(group)
group = OptionGroup(parser, "Geometry","")
group.add_option('-g',
'--grid',
dest = 'grid',
type = 'int', nargs = 3, metavar = ' '.join(['int']*3),
help = 'a,b,c grid of hexahedral box [auto]')
group.add_option('-s',
'--size',
dest = 'size',
type = 'float', nargs = 3, metavar=' '.join(['float']*3),
help = 'x,y,z size of hexahedral box [auto]')
group.add_option('-o',
'--origin',
dest = 'origin',
type = 'float', nargs = 3, metavar=' '.join(['float']*3),
help = 'origin of grid')
parser.add_option_group(group)
group = OptionGroup(parser, "Seeds","")
group.add_option('-p',
'--pos', '--seedposition',
dest = 'pos',
type = 'string', metavar = 'string',
help = 'label of coordinates [%default]')
group.add_option('-w',
'--weight',
dest = 'weight',
type = 'string', metavar = 'string',
help = 'label of weights [%default]')
group.add_option('-m',
'--microstructure',
dest = 'microstructure',
type = 'string', metavar = 'string',
help = 'label of microstructures [%default]')
group.add_option('-e',
'--eulers',
dest = 'eulers',
type = 'string', metavar = 'string',
help = 'label of Euler angles [%default]')
group.add_option('--axes',
dest = 'axes',
type = 'string', nargs = 3, metavar = ' '.join(['string']*3),
help = 'orientation coordinate frame in terms of position coordinate frame')
parser.add_option_group(group)
group = OptionGroup(parser, "Configuration","")
group.add_option('--homogenization',
dest = 'homogenization',
type = 'int', metavar = 'int',
help = 'homogenization index to be used [%default]')
group.add_option('--crystallite',
dest = 'crystallite',
type = 'int', metavar = 'int',
help = 'crystallite index to be used [%default]')
group.add_option('--phase',
dest = 'phase',
type = 'int', metavar = 'int',
help = 'phase index to be used [%default]')
parser.add_option_group(group)
parser.set_defaults(pos = 'pos',
weight = 'weight', weight = 'weight',
microstructure = 'microstructure', microstructure = 'microstructure',
eulers = 'eulerangles', eulers = 'euler',
homogenization = 1, homogenization = 1,
crystallite = 1, crystallite = 1,
phase = 1, phase = 1,
secondphase = 0.0,
cpus = 2, cpus = 2,
laguerre = False, laguerre = False,
nonperiodic = False, nonperiodic = False,
randomSeed = None,
) )
(options,filenames) = parser.parse_args() (options,filenames) = parser.parse_args()
if options.secondphase > 1.0 or options.secondphase < 0.0:
parser.error('volume fraction of second phase ({}) out of bounds.'.format(options.secondphase))
# --- loop over input files ------------------------------------------------------------------------- # --- loop over input files -------------------------------------------------------------------------
if filenames == []: filenames = [None] if filenames == []: filenames = [None]
for name in filenames: for name in filenames:
try: try: table = damask.ASCIItable(name = name,
table = damask.ASCIItable(name = name, outname = os.path.splitext(name)[0]+'.geom' if name else name,
outname = os.path.splitext(name)[-2]+'.geom' if name else name,
buffered = False) buffered = False)
except: continue except: continue
damask.util.report(scriptName,name) damask.util.report(scriptName,name)
@ -232,11 +246,11 @@ for name in filenames:
info['size'][i] = float(info['grid'][i])/max(info['grid']) # normalize to grid info['size'][i] = float(info['grid'][i])/max(info['grid']) # normalize to grid
remarks.append('rescaling size {} to {}...'.format({0:'x',1:'y',2:'z'}[i],info['size'][i])) remarks.append('rescaling size {} to {}...'.format({0:'x',1:'y',2:'z'}[i],info['size'][i]))
if table.label_dimension(options.position) != 3: if table.label_dimension(options.pos) != 3:
errors.append('position columns "{}" have dimension {}.'.format(options.position, errors.append('seed positions "{}" have dimension {}.'.format(options.pos,
table.label_dimension(options.position))) table.label_dimension(options.pos)))
else: else:
labels += [options.position] labels += [options.pos]
if not hasEulers: remarks.append('missing seed orientations...') if not hasEulers: remarks.append('missing seed orientations...')
else: labels += [options.eulers] else: labels += [options.eulers]
@ -254,14 +268,13 @@ for name in filenames:
# ------------------------------------------ read seeds --------------------------------------- # ------------------------------------------ read seeds ---------------------------------------
table.data_readArray(labels) table.data_readArray(labels)
coords = table.data[:,table.label_index(options.position):table.label_index(options.position)+3]\ coords = table.data[:,table.label_indexrange(options.pos)] * info['size']
* info['size'] eulers = table.data[:,table.label_indexrange(options.eulers)] if hasEulers \
eulers = table.data[:,table.label_index(options.eulers ):table.label_index(options.eulers )+3]\ else np.zeros(3*len(coords))
if hasEulers else np.zeros(3*len(coords)) grains = table.data[:,table.label_indexrange(options.microstructure)].astype('i') if hasGrains \
grains = table.data[:,table.label_index(options.microstructure)].astype('i')\ else 1+np.arange(len(coords))
if hasGrains else 1+np.arange(len(coords)) weights = table.data[:,table.label_indexrange(options.weight)] if hasWeights \
weights = table.data[:,table.label_index(options.weight)]\ else np.zeros(len(coords))
if hasWeights else np.zeros(len(coords))
grainIDs = np.unique(grains).astype('i') grainIDs = np.unique(grains).astype('i')
NgrainIDs = len(grainIDs) NgrainIDs = len(grainIDs)
@ -283,50 +296,38 @@ for name in filenames:
if info['homogenization'] == 0: info['homogenization'] = options.homogenization if info['homogenization'] == 0: info['homogenization'] = options.homogenization
damask.util.croak(['grid a b c: %s'%(' x '.join(map(str,info['grid']))), damask.util.report_geom(info,['grid','size','origin','homogenization',])
'size x y z: %s'%(' x '.join(map(str,info['size']))), damask.util.croak(['microstructures: {}{}'.format(info['microstructures'],
'origin x y z: %s'%(' : '.join(map(str,info['origin']))), (' out of {}'.format(NgrainIDs) if NgrainIDs != info['microstructures'] else '')),
'homogenization: %i'%info['homogenization'],
'microstructures: %i%s'%(info['microstructures'],
(' out of %i'%NgrainIDs if NgrainIDs != info['microstructures'] else '')),
]) ])
config_header = [] config_header = []
formatwidth = 1+int(math.log10(NgrainIDs)) formatwidth = 1+int(math.log10(NgrainIDs))
phase = options.phase * np.ones(NgrainIDs,'i')
if int(options.secondphase*info['microstructures']) > 0:
phase[0:int(options.secondphase*info['microstructures'])] += 1 # alter fraction 'options.secondphase' of used grainIDs
randomSeed = options.randomSeed if options.randomSeed \
else int(os.urandom(4).encode('hex'), 16) # random seed for second phase
np.random.seed(randomSeed)
np.random.shuffle(phase)
config_header += ['# random seed (phase shuffling): {}'.format(randomSeed)]
config_header += ['<microstructure>'] config_header += ['<microstructure>']
for i,ID in enumerate(grainIDs): for i,ID in enumerate(grainIDs):
config_header += ['[Grain%s]'%(str(ID).zfill(formatwidth)), config_header += ['[Grain{}]'.format(str(ID).zfill(formatwidth)),
'crystallite %i'%options.crystallite, 'crystallite {}'.format(options.crystallite),
'(constituent)\tphase %i\ttexture %s\tfraction 1.0'%(phase[i],str(ID).rjust(formatwidth)), '(constituent)\tphase {}\ttexture {}\tfraction 1.0'.format(options.phase,str(ID).rjust(formatwidth)),
] ]
if hasEulers: if hasEulers:
config_header += ['<texture>'] config_header += ['<texture>']
for ID in grainIDs: for ID in grainIDs:
eulerID = np.nonzero(grains == ID)[0][0] # find first occurrence of this grain id eulerID = np.nonzero(grains == ID)[0][0] # find first occurrence of this grain id
config_header += ['[Grain%s]'%(str(ID).zfill(formatwidth)), config_header += ['[Grain{}]'.format(str(ID).zfill(formatwidth)),
'(gauss)\tphi1 %g\tPhi %g\tphi2 %g\tscatter 0.0\tfraction 1.0'%tuple(eulers[eulerID]) '(gauss)\tphi1 {:g}\tPhi {:g}\tphi2 {:g}\tscatter 0.0\tfraction 1.0'.format(*eulers[eulerID])
] ]
if options.axes is not None: config_header.append('axes\t%s %s %s'%tuple(options.axes)) if options.axes is not None: config_header.append('axes\t{} {} {}'.format(*options.axes))
table.labels_clear() table.labels_clear()
table.info_clear() table.info_clear()
table.info_append([ table.info_append([
scriptID + ' ' + ' '.join(sys.argv[1:]), scriptID + ' ' + ' '.join(sys.argv[1:]),
"grid\ta {grid[0]}\tb {grid[1]}\tc {grid[2]}".format(grid=info['grid']), "grid\ta {}\tb {}\tc {}".format(*info['grid']),
"size\tx {size[0]}\ty {size[1]}\tz {size[2]}".format(size=info['size']), "size\tx {}\ty {}\tz {}".format(*info['size']),
"origin\tx {origin[0]}\ty {origin[1]}\tz {origin[2]}".format(origin=info['origin']), "origin\tx {}\ty {}\tz {}".format(*info['origin']),
"homogenization\t{homog}".format(homog=info['homogenization']), "homogenization\t{}".format(info['homogenization']),
"microstructures\t{microstructures}".format(microstructures=info['microstructures']), "microstructures\t{}".format(info['microstructures']),
config_header, config_header,
]) ])
table.head_write() table.head_write()

23
processing/pre/geom_toTable.py
View File

@ -19,16 +19,18 @@ Translate geom description into ASCIItable containing 1/2/3_pos and microstructu
""", version = scriptID) """, version = scriptID)
parser.add_option('-p','--position', parser.add_option('-p',
dest = 'position', '--pos', '--position',
dest = 'pos',
type = 'string', metavar = 'string', type = 'string', metavar = 'string',
help = 'column label for position [%default]') help = 'label of coordinates [%default]')
parser.add_option('-m','--microstructure', parser.add_option('-m',
'--microstructure',
dest = 'microstructure', dest = 'microstructure',
type = 'string', metavar = 'string', type = 'string', metavar = 'string',
help = 'column label for microstructure index [%default]') help = 'label of microstructure index [%default]')
parser.set_defaults(position = 'pos', parser.set_defaults(pos = 'pos',
microstructure = 'microstructure', microstructure = 'microstructure',
) )
@ -39,10 +41,11 @@ parser.set_defaults(position = 'pos',
if filenames == []: filenames = [None] if filenames == []: filenames = [None]
for name in filenames: for name in filenames:
try: try: table = damask.ASCIItable(name = name,
table = damask.ASCIItable(name = name,
outname = os.path.splitext(name)[0]+'.txt' if name else name, outname = os.path.splitext(name)[0]+'.txt' if name else name,
buffered = False, labeled = False) buffered = False,
labeled = False,
)
except: continue except: continue
damask.util.report(scriptName,name) damask.util.report(scriptName,name)
@ -75,7 +78,7 @@ for name in filenames:
table.info_clear() table.info_clear()
table.info_append(extra_header + [scriptID + '\t' + ' '.join(sys.argv[1:])]) table.info_append(extra_header + [scriptID + '\t' + ' '.join(sys.argv[1:])])
table.labels_clear() table.labels_clear()
table.labels_append(['{}_{}'.format(1+i,options.position) for i in xrange(3)]+[options.microstructure]) table.labels_append(['{}_{}'.format(1+i,options.pos) for i in xrange(3)]+[options.microstructure])
table.head_write() table.head_write()
table.output_flush() table.output_flush()

125
processing/pre/seeds_check.py
View File

@ -1,125 +0,0 @@
#!/usr/bin/env python
# -*- coding: UTF-8 no BOM -*-
import os,sys,vtk
import numpy as np
import damask
from optparse import OptionParser
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
#--------------------------------------------------------------------------------------------------
# MAIN
#--------------------------------------------------------------------------------------------------
parser = OptionParser(option_class=damask.extendableOption, usage='%prog [seedsfile[s]]', description = """
Produce VTK point mesh from seeds file
""", version = scriptID)
parser.add_option('-s', '--size',
dest = 'size',
type = 'float', nargs = 3, metavar = 'float float float',
help = 'x,y,z size of hexahedral box [1.0 along largest grid point number]')
parser.add_option('-p','--position',
dest = 'position',
type = 'string', metavar = 'string',
help = 'column label for coordinates [%default]')
parser.set_defaults(size = [0.0,0.0,0.0],
position = 'pos',
)
(options, filenames) = parser.parse_args()
# --- loop over input files -------------------------------------------------------------------------
if filenames == []: filenames = [None]
for name in filenames:
try:
table = damask.ASCIItable(name = name,
buffered = False, readonly = True)
except: continue
damask.util.report(scriptName,name)
# --- interpret header ----------------------------------------------------------------------------
table.head_read()
info,extra_header = table.head_getGeom()
damask.util.croak(['grid a b c: %s'%(' x '.join(map(str,info['grid']))),
'size x y z: %s'%(' x '.join(map(str,info['size']))),
'origin x y z: %s'%(' : '.join(map(str,info['origin']))),
'homogenization: %i'%info['homogenization'],
'microstructures: %i'%info['microstructures'],
])
remarks = []
errors = []
if np.any(info['grid'] < 1): remarks.append('invalid grid a b c.')
if np.any(info['size'] <= 0.0) \
and np.all(info['grid'] < 1): errors.append('invalid size x y z.')
else:
for i in xrange(3):
if info['size'][i] <= 0.0: # any invalid size?
info['size'][i] = float(info['grid'][i])/max(info['grid']) # normalize to grid
remarks.append('rescaling size {} to {}...'.format({0:'x',1:'y',2:'z'}[i],info['size'][i]))
if table.label_dimension(options.position) != 3:
errors.append('columns "{}" have dimension {}'.format(options.position,table.label_dimension(options.position)))
if remarks != []: damask.util.croak(remarks)
if errors != []:
damask.util.croak(errors)
table.close(dismiss=True)
continue
labels = ['{dim}_{label}'.format(dim = 1+i,label = options.position) for i in xrange(3)]
hasGrains = table.label_index('microstructure') != -1
labels += ['microstructure'] if hasGrains else []
table.data_readArray(labels) # read ASCIItable columns
coords = table.data[:,:3]*info['size'] # assign coordinates (rescaled to box size)
grain = table.data[:,3].astype('i') if hasGrains else 1+np.arange(len(coords),dtype='i') # assign grains
# --- generate grid --------------------------------------------------------------------------------
grid = vtk.vtkUnstructuredGrid()
pts = vtk.vtkPoints()
# --- process microstructure information -----------------------------------------------------------
IDs = vtk.vtkIntArray()
IDs.SetNumberOfComponents(1)
IDs.SetName("GrainID")
for i,item in enumerate(coords):
IDs.InsertNextValue(grain[i])
pid = pts.InsertNextPoint(item[0:3])
pointIds = vtk.vtkIdList()
pointIds.InsertId(0, pid)
grid.InsertNextCell(1, pointIds)
grid.SetPoints(pts)
grid.GetCellData().AddArray(IDs)
#--- write data -----------------------------------------------------------------------------------
if name:
writer = vtk.vtkXMLRectilinearGridWriter()
(directory,filename) = os.path.split(name)
writer.SetDataModeToBinary()
writer.SetCompressorTypeToZLib()
writer.SetFileName(os.path.join(directory,os.path.splitext(filename)[0]+'.'+writer.GetDefaultFileExtension()))
else:
writer = vtk.vtkDataSetWriter()
writer.WriteToOutputStringOn()
writer.SetHeader('# powered by '+scriptID)
if vtk.VTK_MAJOR_VERSION <= 5: writer.SetInput(grid)
else: writer.SetInputData(grid)
writer.Write()
if name is None: sys.stdout.write(writer.GetOutputString()[0:writer.GetOutputStringLength()])
table.close()

12
processing/pre/seeds_check.sh Executable file
View File

@ -0,0 +1,12 @@
#!/bin/bash
for seeds in "$@"
do
vtk_pointcloud $seeds
vtk_addPointCloudData $seeds \
--scalar microstructure,weight \
--inplace \
--vtk ${seeds%.*}.vtp \
done

45
processing/pre/seeds_fromGeom.py
View File

@ -14,26 +14,30 @@ scriptID = ' '.join([scriptName,damask.version])
#-------------------------------------------------------------------------------------------------- #--------------------------------------------------------------------------------------------------
parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """ parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """
Create seed file taking microstructure indices from given geom file but excluding black-listed grains. Create seed file taking microstructure indices from given geom file.
Indices can be black-listed or white-listed.
""", version = scriptID) """, version = scriptID)
parser.add_option('-w','--white', parser.add_option('-w',
'--white',
action = 'extend', metavar = '<int LIST>', action = 'extend', metavar = '<int LIST>',
dest = 'whitelist', dest = 'whitelist',
help = 'whitelist of grain IDs') help = 'whitelist of grain IDs')
parser.add_option('-b','--black', parser.add_option('-b',
'--black',
action = 'extend', metavar = '<int LIST>', action = 'extend', metavar = '<int LIST>',
dest = 'blacklist', dest = 'blacklist',
help = 'blacklist of grain IDs') help = 'blacklist of grain IDs')
parser.add_option('-p','--position', parser.add_option('-p',
'--pos', '--seedposition',
dest = 'position', dest = 'position',
type = 'string', metavar = 'string', type = 'string', metavar = 'string',
help = 'column label for coordinates [%default]') help = 'label of coordinates [%default]')
parser.set_defaults(whitelist = [], parser.set_defaults(whitelist = [],
blacklist = [], blacklist = [],
position = 'pos', pos = 'pos',
) )
(options,filenames) = parser.parse_args() (options,filenames) = parser.parse_args()
@ -46,25 +50,18 @@ options.blacklist = map(int,options.blacklist)
if filenames == []: filenames = [None] if filenames == []: filenames = [None]
for name in filenames: for name in filenames:
try: try: table = damask.ASCIItable(name = name,
table = damask.ASCIItable(name = name,
outname = os.path.splitext(name)[0]+'.seeds' if name else name, outname = os.path.splitext(name)[0]+'.seeds' if name else name,
buffered = False, labeled = False) buffered = False,
except: labeled = False)
continue except: continue
damask.util.report(scriptName,name) damask.util.report(scriptName,name)
# --- interpret header ---------------------------------------------------------------------------- # --- interpret header ----------------------------------------------------------------------------
table.head_read() table.head_read()
info,extra_header = table.head_getGeom() info,extra_header = table.head_getGeom()
damask.util.report_geom(info)
damask.util.croak(['grid a b c: %s'%(' x '.join(map(str,info['grid']))),
'size x y z: %s'%(' x '.join(map(str,info['size']))),
'origin x y z: %s'%(' : '.join(map(str,info['origin']))),
'homogenization: %i'%info['homogenization'],
'microstructures: %i'%info['microstructures'],
])
errors = [] errors = []
if np.any(info['grid'] < 1): errors.append('invalid grid a b c.') if np.any(info['grid'] < 1): errors.append('invalid grid a b c.')
@ -98,14 +95,14 @@ for name in filenames:
table.info_clear() table.info_clear()
table.info_append(extra_header+[ table.info_append(extra_header+[
scriptID + ' ' + ' '.join(sys.argv[1:]), scriptID + ' ' + ' '.join(sys.argv[1:]),
"grid\ta {grid[0]}\tb {grid[1]}\tc {grid[2]}".format(grid=info['grid']), "grid\ta {}\tb {}\tc {}".format(*info['grid']),
"size\tx {size[0]}\ty {size[1]}\tz {size[2]}".format(size=info['size']), "size\tx {}\ty {}\tz {}".format(*info['size']),
"origin\tx {origin[0]}\ty {origin[1]}\tz {origin[2]}".format(origin=info['origin']), "origin\tx {}\ty {}\tz {}".format(*info['origin']),
"homogenization\t{homog}".format(homog=info['homogenization']), "homogenization\t{}".format(info['homogenization']),
"microstructures\t{microstructures}".format(microstructures=info['microstructures']), "microstructures\t{}".format(info['microstructures']),
]) ])
table.labels_clear() table.labels_clear()
table.labels_append(['{dim}_{label}'.format(dim = 1+i,label = options.position) for i in range(3)]+['microstructure']) table.labels_append(['{dim}_{label}'.format(dim = 1+i,label = options.pos) for i in range(3)]+['microstructure'])
table.head_write() table.head_write()
table.output_flush() table.output_flush()

57
processing/pre/seeds_fromRandom.py
View File

@ -34,28 +34,36 @@ Reports positions with random crystal orientations in seeds file format to STDOU
""", version = scriptID) """, version = scriptID)
parser.add_option('-N', dest='N', parser.add_option('-N',
dest = 'N',
type = 'int', metavar = 'int', type = 'int', metavar = 'int',
help = 'number of seed points to distribute [%default]') help = 'number of seed points [%default]')
parser.add_option('-g','--grid', parser.add_option('-g',
'--grid',
dest = 'grid', dest = 'grid',
type = 'int', nargs = 3, metavar = 'int int int', type = 'int', nargs = 3, metavar = 'int int int',
help='min a,b,c grid of hexahedral box %default') help='min a,b,c grid of hexahedral box %default')
parser.add_option('-m', '--microstructure', parser.add_option('-m',
'--microstructure',
dest = 'microstructure', dest = 'microstructure',
type = 'int', metavar = 'int', type = 'int', metavar = 'int',
help = 'first microstructure index [%default]') help = 'first microstructure index [%default]')
parser.add_option('-r', '--rnd', parser.add_option('-r',
'--rnd',
dest = 'randomSeed', type = 'int', metavar = 'int', dest = 'randomSeed', type = 'int', metavar = 'int',
help = 'seed of random number generator [%default]') help = 'seed of random number generator [%default]')
parser.add_option('--format',
dest = 'format', type = 'string', metavar = 'string',
help = 'output number format [auto]')
group = OptionGroup(parser, "Laguerre Tessellation Options", group = OptionGroup(parser, "Laguerre Tessellation",
"Parameters determining shape of weight distribution of seed points" "Parameters determining shape of weight distribution of seed points"
) )
group.add_option('-w', '--weights', group.add_option( '-w',
'--weights',
action = 'store_true', action = 'store_true',
dest = 'weights', dest = 'weights',
help = 'assign random weigts to seed points for Laguerre tessellation [%default]') help = 'assign random weights to seed points for Laguerre tessellation [%default]')
group.add_option( '--max', group.add_option( '--max',
dest = 'max', dest = 'max',
type = 'float', metavar = 'float', type = 'float', metavar = 'float',
@ -70,21 +78,23 @@ group.add_option('--sigma',
help='standard deviation of normal distribution for weights [%default]') help='standard deviation of normal distribution for weights [%default]')
parser.add_option_group(group) parser.add_option_group(group)
group = OptionGroup(parser, "Selective Seeding Options", group = OptionGroup(parser, "Selective Seeding",
"More uniform distribution of seed points using Mitchell\'s Best Candidate Algorithm" "More uniform distribution of seed points using Mitchell's Best Candidate Algorithm"
) )
group.add_option('-s','--selective', group.add_option( '-s',
'--selective',
action = 'store_true', action = 'store_true',
dest = 'selective', dest = 'selective',
help = 'selective picking of seed points from random seed points [%default]') help = 'selective picking of seed points from random seed points [%default]')
group.add_option('-f','--force', group.add_option( '-f',
'--force',
action = 'store_true', action = 'store_true',
dest = 'force', dest = 'force',
help = 'try selective picking despite large seed point number [%default]') help = 'try selective picking despite large seed point number [%default]')
group.add_option( '--distance', group.add_option( '--distance',
dest = 'distance', dest = 'distance',
type = 'float', metavar = 'float', type = 'float', metavar = 'float',
help = 'minimum distance to the next neighbor [%default]') help = 'minimum distance to next neighbor [%default]')
group.add_option( '--numCandidates', group.add_option( '--numCandidates',
dest = 'numCandidates', dest = 'numCandidates',
type = 'int', metavar = 'int', type = 'int', metavar = 'int',
@ -103,6 +113,7 @@ parser.set_defaults(randomSeed = None,
force = False, force = False,
distance = 0.2, distance = 0.2,
numCandidates = 10, numCandidates = 10,
format = None,
) )
(options,filenames) = parser.parse_args() (options,filenames) = parser.parse_args()
@ -120,11 +131,9 @@ random.seed(options.randomSeed)
if filenames == []: filenames = [None] if filenames == []: filenames = [None]
for name in filenames: for name in filenames:
try: try: table = damask.ASCIItable(outname = name,
table = damask.ASCIItable(outname = name,
buffered = False) buffered = False)
except: except: continue
continue
damask.util.report(scriptName,name) damask.util.report(scriptName,name)
# --- sanity checks ------------------------------------------------------------------------- # --- sanity checks -------------------------------------------------------------------------
@ -132,7 +141,7 @@ for name in filenames:
remarks = [] remarks = []
errors = [] errors = []
if gridSize == 0: if gridSize == 0:
errors.append('zero grid dimension for %s.'%(', '.join([['a','b','c'][x] for x in np.where(options.grid == 0)[0]]))) errors.append('zero grid dimension for {}.'.format(', '.join([['a','b','c'][x] for x in np.where(options.grid == 0)[0]])))
if options.N > gridSize/10.: errors.append('seed count exceeds 0.1 of grid points.') if options.N > gridSize/10.: errors.append('seed count exceeds 0.1 of grid points.')
if options.selective and 4./3.*math.pi*(options.distance/2.)**3*options.N > 0.5: if options.selective and 4./3.*math.pi*(options.distance/2.)**3*options.N > 0.5:
(remarks if options.force else errors).append('maximum recommended seed point count for given distance is {}.{}'. (remarks if options.force else errors).append('maximum recommended seed point count for given distance is {}.{}'.
@ -182,10 +191,8 @@ for name in filenames:
seeds = seeds.T # prepare shape for stacking seeds = seeds.T # prepare shape for stacking
if options.weights: if options.weights:
if options.max > 0.0: weights = [np.random.uniform(low = 0, high = options.max, size = options.N)] if options.max > 0.0 \
weights = [np.random.uniform(low = 0, high = options.max, size = options.N)] else [np.random.normal(loc = options.mean, scale = options.sigma, size = options.N)]
else:
weights = [np.random.normal(loc = options.mean, scale = options.sigma, size = options.N)]
else: else:
weights = [] weights = []
seeds = np.transpose(np.vstack(tuple([seeds, seeds = np.transpose(np.vstack(tuple([seeds,
@ -200,13 +207,13 @@ for name in filenames:
table.info_clear() table.info_clear()
table.info_append([ table.info_append([
scriptID + ' ' + ' '.join(sys.argv[1:]), scriptID + ' ' + ' '.join(sys.argv[1:]),
"grid\ta {grid[0]}\tb {grid[1]}\tc {grid[2]}".format(grid=options.grid), "grid\ta {}\tb {}\tc {}".format(*options.grid),
"microstructures\t{}".format(options.N), "microstructures\t{}".format(options.N),
"randomSeed\t{}".format(options.randomSeed), "randomSeed\t{}".format(options.randomSeed),
]) ])
table.labels_clear() table.labels_clear()
table.labels_append( ['{dim}_{label}'.format(dim = 1+k,label = 'pos') for k in xrange(3)] + table.labels_append( ['{dim}_{label}'.format(dim = 1+k,label = 'pos') for k in xrange(3)] +
['{dim}_{label}'.format(dim = 1+k,label = 'eulerangles') for k in xrange(3)] + ['{dim}_{label}'.format(dim = 1+k,label = 'euler') for k in xrange(3)] +
['microstructure'] + ['microstructure'] +
(['weight'] if options.weights else [])) (['weight'] if options.weights else []))
table.head_write() table.head_write()
@ -215,7 +222,7 @@ for name in filenames:
# --- write seeds information ------------------------------------------------------------ # --- write seeds information ------------------------------------------------------------
table.data = seeds table.data = seeds
table.data_writeArray() table.data_writeArray(fmt = options.format)
# --- output finalization -------------------------------------------------------------------------- # --- output finalization --------------------------------------------------------------------------

42
processing/pre/seeds_fromTable.py
View File

@ -23,29 +23,38 @@ Examples:
""", version = scriptID) """, version = scriptID)
parser.add_option('-p', '--positions', parser.add_option('-p',
'--pos', '--seedposition',
dest = 'pos', dest = 'pos',
type = 'string', metavar = 'string', type = 'string', metavar = 'string',
help = 'coordinate label [%default]') help = 'label of coordinates [%default]')
parser.add_option('--boundingbox', parser.add_option('--boundingbox',
dest = 'box', dest = 'box',
type = 'float', nargs = 6, metavar = ' '.join(['float']*6), type = 'float', nargs = 6, metavar = ' '.join(['float']*6),
help = 'min (x,y,z) and max (x,y,z) coordinates of bounding box [tight]') help = 'min (x,y,z) and max (x,y,z) coordinates of bounding box [tight]')
parser.add_option('-i', '--index', parser.add_option('-m',
dest = 'index', '--microstructure',
dest = 'microstructure',
type = 'string', metavar = 'string', type = 'string', metavar = 'string',
help = 'microstructure index label [%default]') help = 'label of microstructures [%default]')
parser.add_option('-w','--white', parser.add_option('--weight',
dest = 'weight',
type = 'string', metavar = 'string',
help = 'label of weights [%default]')
parser.add_option('-w',
'--white',
dest = 'whitelist', dest = 'whitelist',
action = 'extend', metavar = '<int LIST>', action = 'extend', metavar = '<int LIST>',
help = 'whitelist of microstructure indices') help = 'whitelist of microstructure indices')
parser.add_option('-b','--black', parser.add_option('-b',
'--black',
dest = 'blacklist', dest = 'blacklist',
action = 'extend', metavar = '<int LIST>', action = 'extend', metavar = '<int LIST>',
help = 'blacklist of microstructure indices') help = 'blacklist of microstructure indices')
parser.set_defaults(pos = 'pos', parser.set_defaults(pos = 'pos',
index ='microstructure', microstructure = 'microstructure',
weight = None,
) )
(options,filenames) = parser.parse_args() (options,filenames) = parser.parse_args()
@ -58,8 +67,7 @@ if options.blacklist is not None: options.blacklist = map(int,options.blacklist)
if filenames == []: filenames = [None] if filenames == []: filenames = [None]
for name in filenames: for name in filenames:
try: try: table = damask.ASCIItable(name = name,
table = damask.ASCIItable(name = name,
outname = os.path.splitext(name)[0]+'.seeds' if name else name, outname = os.path.splitext(name)[0]+'.seeds' if name else name,
buffered = False) buffered = False)
except: continue except: continue
@ -69,14 +77,17 @@ for name in filenames:
# ------------------------------------------ sanity checks --------------------------------------- # ------------------------------------------ sanity checks ---------------------------------------
missing_labels = table.data_readArray([options.pos,options.index]) missing_labels = table.data_readArray([options.pos,options.microstructure] +
([options.weight] if options.weight else []))
errors = [] errors = []
if len(missing_labels) > 0: if len(missing_labels) > 0:
errors.append('column{} {} not found'.format('s' if len(missing_labels) > 1 else '', errors.append('column{} {} not found'.format('s' if len(missing_labels) > 1 else '',
', '.join(missing_labels))) ', '.join(missing_labels)))
for label, dim in {options.pos: 3, input = {options.pos: 3,
options.index: 1}.iteritems(): options.microstructure: 1,}
if options.weight: input.update({options.weight: 1})
for label, dim in input.iteritems():
if table.label_dimension(label) != dim: if table.label_dimension(label) != dim:
errors.append('column {} has wrong dimension'.format(label)) errors.append('column {} has wrong dimension'.format(label))
@ -113,11 +124,12 @@ for name in filenames:
table.info = [ table.info = [
scriptID, scriptID,
'size %s'%(' '.join(list(itertools.chain.from_iterable(zip(['x','y','z'], 'size {}'.format(' '.join(list(itertools.chain.from_iterable(zip(['x','y','z'],
map(str,boundingBox[1,:]-boundingBox[0,:])))))), map(str,boundingBox[1,:]-boundingBox[0,:])))))),
] ]
table.labels_clear() table.labels_clear()
table.labels_append(['1_pos','2_pos','3_pos','microstructure']) # implicitly switching label processing/writing on table.labels_append(['1_pos','2_pos','3_pos','microstructure'] +
['weight'] if options.weight else []) # implicitly switching label processing/writing on
table.head_write() table.head_write()
# ------------------------------------------ output result --------------------------------------- # ------------------------------------------ output result ---------------------------------------