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:
Aritra Chakraborty 2016-05-11 14:53:21 -04:00
parent 739fc0f9ad 7ea5ce73ae
commit 4a273f9cb3
207 changed files with 188973 additions and 193107 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
.gitattributes vendored Normal file
View File

@ -0,0 +1,10 @@
# 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
installation/mods_Abaqus/abaqus_v6_windows.env eol=crlf
# Denote all files that are truly binary and should not be modified.
*.png binary
*.jpg binary
*.cae binary

19
DAMASK_env.bat
View File

@ -1,19 +0,0 @@
:: sets up an environment for DAMASK on Windows
:: usage: call DAMASK_env.bat
@echo off
set LOCATION=%~dp0
set DAMASK_ROOT=%LOCATION%\DAMASK
set DAMASK_NUM_THREADS=2
chcp 1252
Title Düsseldorf Advanced Materials Simulation Kit - DAMASK, MPIE Düsseldorf
echo.
echo Düsseldorf Advanced Materials Simulation Kit - DAMASK
echo Max-Planck-Institut für Eisenforschung, Düsseldorf
echo http://damask.mpie.de
echo.
echo Preparing environment ...
echo DAMASK_ROOT=%DAMASK_ROOT%
echo DAMASK_NUM_THREADS=%DAMASK_NUM_THREADS%
set DAMASK_BIN=%DAMASK_ROOT%\bin
set PATH=%PATH%;%DAMASK_BIN%
set PYTHONPATH=%PYTHONPATH%;%DAMASK_ROOT%\lib

13
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
@ -53,8 +53,11 @@ if [ ! -z "$PS1" ]; then
[[ "x$SOLVER" != "x" ]] && echo "Spectral Solver $SOLVER" [[ "x$SOLVER" != "x" ]] && echo "Spectral Solver $SOLVER"
[[ "x$PROCESSING" != "x" ]] && echo "Post Processing $PROCESSING" [[ "x$PROCESSING" != "x" ]] && echo "Post Processing $PROCESSING"
echo "Multithreading DAMASK_NUM_THREADS=$DAMASK_NUM_THREADS" echo "Multithreading DAMASK_NUM_THREADS=$DAMASK_NUM_THREADS"
[[ "x$PETSC_DIR" != "x" ]] && echo "PETSc location $PETSC_DIR" && \ if [ "x$PETSC_DIR" != "x" ]; then
[[ `readlink -f $PETSC_DIR` == $PETSC_DIR ]] || echo " ~~> "`readlink -f $PETSC_DIR` echo "PETSc location $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
[[ "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"
echo echo

6
DAMASK_env.zsh
View File

@ -51,8 +51,10 @@ if [ ! -z "$PS1" ]; then
[[ "x$SOLVER" != "x" ]] && echo "Spectral Solver $SOLVER" [[ "x$SOLVER" != "x" ]] && echo "Spectral Solver $SOLVER"
[[ "x$PROCESSING" != "x" ]] && echo "Post Processing $PROCESSING" [[ "x$PROCESSING" != "x" ]] && echo "Post Processing $PROCESSING"
echo "Multithreading DAMASK_NUM_THREADS=$DAMASK_NUM_THREADS" echo "Multithreading DAMASK_NUM_THREADS=$DAMASK_NUM_THREADS"
[[ "x$PETSC_DIR" != "x" ]] && echo "PETSc location $PETSC_DIR" && \ if [ "x$PETSC_DIR" != "x" ]; then
[[ `readlink -f $PETSC_DIR` == $PETSC_DIR ]] || echo " ~~> "`readlink -f $PETSC_DIR` echo "PETSc location $PETSC_DIR"
[[ `readlink -f $PETSC_DIR` == $PETSC_DIR ]] || echo " ~~> "`readlink -f $PETSC_DIR`
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"
echo echo

2
VERSION
View File

@ -1 +1 @@
v2.0.0-36-g2b7524e v2.0.0-219-g297d5ca

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

2
code/.gitignore vendored
View File

@ -1 +1,3 @@
DAMASK_marc*.f90 DAMASK_marc*.f90
quit__genmod.f90
*.marc

36
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, &
@ -59,8 +60,6 @@ program DAMASK_spectral
materialpoint_sizeResults, & materialpoint_sizeResults, &
materialpoint_results, & materialpoint_results, &
materialpoint_postResults materialpoint_postResults
use material, only: & use material, only: &
thermal_type, & thermal_type, &
damage_type, & damage_type, &
@ -149,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)
@ -341,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
@ -425,28 +426,33 @@ 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
outputIndex=[(i-1)*((maxByteOut/pReal)/materialpoint_sizeResults)+1, & outputIndex=int([(i-1_pInt)*((maxByteOut/pReal)/materialpoint_sizeResults)+1_pInt, &
min(i*((maxByteOut/pReal)/materialpoint_sizeResults),size(materialpoint_results,3))] min(i*((maxByteOut/pReal)/materialpoint_sizeResults),size(materialpoint_results,3))],pLongInt)
call MPI_file_write(resUnit,reshape(materialpoint_results(:,:,outputIndex(1):outputIndex(2)),& call MPI_file_write(resUnit,reshape(materialpoint_results(:,:,outputIndex(1):outputIndex(2)),&
[(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)
fileOffset = fileOffset + sum(outputSize) ! forward to current file position 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
if (worldrank == 0) & if (worldrank == 0) &
write(6,'(1/,a)') ' ... writing initial configuration to file ........................' write(6,'(1/,a)') ' ... writing initial configuration to file ........................'
endif endif
@ -645,15 +651,17 @@ 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=[(i-1)*maxByteOut/pReal/materialpoint_sizeResults+1, & outputIndex=int([(i-1_pInt)*((maxByteOut/pReal)/materialpoint_sizeResults)+1_pInt, &
min(i*maxByteOut/pReal/materialpoint_sizeResults,size(materialpoint_results,3))] min(i*((maxByteOut/pReal)/materialpoint_sizeResults),size(materialpoint_results,3))],pLongInt)
call MPI_file_write(resUnit,reshape(materialpoint_results(:,:,outputIndex(1):outputIndex(2)),& call MPI_file_write(resUnit,reshape(materialpoint_results(:,:,outputIndex(1):outputIndex(2)),&
[(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)
fileOffset = fileOffset + sum(outputSize) ! forward to current file position 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
endif endif
if( loadCases(currentLoadCase)%restartFrequency > 0_pInt .and. & ! at frequency of writing restart information set restart parameter for FEsolving if( loadCases(currentLoadCase)%restartFrequency > 0_pInt .and. & ! at frequency of writing restart information set restart parameter for FEsolving
mod(inc,loadCases(currentLoadCase)%restartFrequency) == 0_pInt) then ! first call to CPFEM_general will write? mod(inc,loadCases(currentLoadCase)%restartFrequency) == 0_pInt) then ! first call to CPFEM_general will write?
@ -700,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 ;)

547
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,7 +1648,9 @@ 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
case (900_pInt) case (900_pInt)
@ -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)

114
code/crystallite.f90
View File

@ -258,7 +258,8 @@ subroutine crystallite_init
allocate(crystallite_orientation(4,cMax,iMax,eMax), source=0.0_pReal) allocate(crystallite_orientation(4,cMax,iMax,eMax), source=0.0_pReal)
allocate(crystallite_orientation0(4,cMax,iMax,eMax), source=0.0_pReal) allocate(crystallite_orientation0(4,cMax,iMax,eMax), source=0.0_pReal)
allocate(crystallite_rotation(4,cMax,iMax,eMax), source=0.0_pReal) allocate(crystallite_rotation(4,cMax,iMax,eMax), source=0.0_pReal)
allocate(crystallite_disorientation(4,nMax,cMax,iMax,eMax), source=0.0_pReal) if (any(plasticState%nonLocal)) &
allocate(crystallite_disorientation(4,nMax,cMax,iMax,eMax),source=0.0_pReal)
allocate(crystallite_localPlasticity(cMax,iMax,eMax), source=.true.) allocate(crystallite_localPlasticity(cMax,iMax,eMax), source=.true.)
allocate(crystallite_requested(cMax,iMax,eMax), source=.false.) allocate(crystallite_requested(cMax,iMax,eMax), source=.false.)
allocate(crystallite_todo(cMax,iMax,eMax), source=.false.) allocate(crystallite_todo(cMax,iMax,eMax), source=.false.)
@ -3569,11 +3570,7 @@ logical function crystallite_integrateStress(&
maxticks maxticks
external :: & external :: &
#if(FLOAT==8)
dgesv dgesv
#elif(FLOAT==4)
sgesv
#endif
!* be pessimistic !* be pessimistic
crystallite_integrateStress = .false. crystallite_integrateStress = .false.
@ -3756,11 +3753,7 @@ logical function crystallite_integrateStress(&
- math_Plain3333to99(math_mul3333xx3333(math_mul3333xx3333(dLp_dT3333,dT_dFe3333),dFe_dLp3333)) - math_Plain3333to99(math_mul3333xx3333(math_mul3333xx3333(dLp_dT3333,dT_dFe3333),dFe_dLp3333))
dRLp_dLp2 = dRLp_dLp ! will be overwritten in first call to LAPACK routine dRLp_dLp2 = dRLp_dLp ! will be overwritten in first call to LAPACK routine
work = math_plain33to9(residuumLp) work = math_plain33to9(residuumLp)
#if(FLOAT==8)
call dgesv(9,1,dRLp_dLp2,9,ipiv,work,9,ierr) ! solve dRLp/dLp * delta Lp = -res for delta Lp call dgesv(9,1,dRLp_dLp2,9,ipiv,work,9,ierr) ! solve dRLp/dLp * delta Lp = -res for delta Lp
#elif(FLOAT==4)
call sgesv(9,1,dRLp_dLp2,9,ipiv,work,9,ierr) ! solve dRLp/dLp * delta Lp = -res for delta Lp
#endif
if (ierr /= 0_pInt) then if (ierr /= 0_pInt) then
#ifndef _OPENMP #ifndef _OPENMP
if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
@ -3849,31 +3842,27 @@ logical function crystallite_integrateStress(&
math_mul3333xx3333(dT_dFi3333, dFi_dLi3333))) & math_mul3333xx3333(dT_dFi3333, dFi_dLi3333))) &
- math_Plain3333to99(math_mul3333xx3333(dLi_dFi3333, dFi_dLi3333)) - math_Plain3333to99(math_mul3333xx3333(dLi_dFi3333, dFi_dLi3333))
work = math_plain33to9(residuumLi) work = math_plain33to9(residuumLi)
#if(FLOAT==8)
call dgesv(9,1,dRLi_dLi,9,ipiv,work,9,ierr) ! solve dRLi/dLp * delta Li = -res for delta Li call dgesv(9,1,dRLi_dLi,9,ipiv,work,9,ierr) ! solve dRLi/dLp * delta Li = -res for delta Li
#elif(FLOAT==4) if (ierr /= 0_pInt) then
call sgesv(9,1,dRLi_dLi,9,ipiv,work,9,ierr) ! solve dRLi/dLp * delta Li = -res for delta Li
#endif
if (ierr /= 0_pInt) then
#ifndef _OPENMP #ifndef _OPENMP
if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3,a,i3)') '<< CRYST >> integrateStress failed on dR/dLi inversion at el ip ipc ', & write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3,a,i3)') '<< CRYST >> integrateStress failed on dR/dLi inversion at el ip ipc ', &
el,mesh_element(1,el),ip,ipc el,mesh_element(1,el),ip,ipc
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt & if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
.and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g)& .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g)&
.or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
write(6,*) write(6,*)
write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dR_dLi',transpose(dRLi_dLi) write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dR_dLi',transpose(dRLi_dLi)
write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dFe_dLi',transpose(math_Plain3333to99(dFe_dLi3333)) write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dFe_dLi',transpose(math_Plain3333to99(dFe_dLi3333))
write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dT_dFi_constitutive',transpose(math_Plain3333to99(dT_dFi3333)) write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dT_dFi_constitutive',transpose(math_Plain3333to99(dT_dFi3333))
write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dLi_dT_constitutive',transpose(math_Plain3333to99(dLi_dT3333)) write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dLi_dT_constitutive',transpose(math_Plain3333to99(dLi_dT3333))
write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> Li_constitutive',math_transpose33(Li_constitutive) write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> Li_constitutive',math_transpose33(Li_constitutive)
write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> Liguess',math_transpose33(Liguess) write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> Liguess',math_transpose33(Liguess)
endif
endif endif
#endif
return
endif endif
#endif
return
endif
deltaLi = - math_plain9to33(work) deltaLi = - math_plain9to33(work)
endif endif
@ -3973,7 +3962,6 @@ subroutine crystallite_orientations
use plastic_nonlocal, only: & use plastic_nonlocal, only: &
plastic_nonlocal_updateCompatibility plastic_nonlocal_updateCompatibility
implicit none implicit none
integer(pInt) & integer(pInt) &
c, & !< counter in integration point component loop c, & !< counter in integration point component loop
@ -3989,50 +3977,51 @@ subroutine crystallite_orientations
! --- CALCULATE ORIENTATION AND LATTICE ROTATION --- ! --- CALCULATE ORIENTATION AND LATTICE ROTATION ---
!$OMP PARALLEL DO PRIVATE(orientation) !$OMP PARALLEL DO PRIVATE(orientation)
do e = FEsolving_execElem(1),FEsolving_execElem(2) do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
do c = 1_pInt,homogenization_Ngrains(mesh_element(3,e)) do c = 1_pInt,homogenization_Ngrains(mesh_element(3,e))
! somehow this subroutine is not threadsafe, so need critical statement here; not clear, what exactly the problem is ! somehow this subroutine is not threadsafe, so need critical statement here; not clear, what exactly the problem is
!$OMP CRITICAL (polarDecomp) !$OMP CRITICAL (polarDecomp)
orientation = math_RtoQ(transpose(math_rotationalPart33(crystallite_Fe(1:3,1:3,c,i,e)))) ! rotational part from polar decomposition as quaternion orientation = math_RtoQ(transpose(math_rotationalPart33(crystallite_Fe(1:3,1:3,c,i,e))))
!$OMP END CRITICAL (polarDecomp) !$OMP END CRITICAL (polarDecomp)
crystallite_rotation(1:4,c,i,e) = lattice_qDisorientation(crystallite_orientation0(1:4,c,i,e), & ! active rotation from ori0 crystallite_rotation(1:4,c,i,e) = lattice_qDisorientation(crystallite_orientation0(1:4,c,i,e), &! active rotation from initial
orientation) ! to current orientation (with no symmetry) orientation) ! to current orientation (with no symmetry)
crystallite_orientation(1:4,c,i,e) = orientation crystallite_orientation(1:4,c,i,e) = orientation
enddo; enddo; enddo enddo; enddo; enddo
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
! --- UPDATE SOME ADDITIONAL VARIABLES THAT ARE NEEDED FOR NONLOCAL MATERIAL --- ! --- UPDATE SOME ADDITIONAL VARIABLES THAT ARE NEEDED FOR NONLOCAL MATERIAL ---
! --- we use crystallite_orientation from above, so need a separate loop ! --- we use crystallite_orientation from above, so need a separate loop
!$OMP PARALLEL DO PRIVATE(myPhase,neighboring_e,neighboring_i,neighboringPhase) nonlocalPresent: if (any(plasticState%nonLocal)) then
!$OMP PARALLEL DO PRIVATE(myPhase,neighboring_e,neighboring_i,neighboringPhase)
do e = FEsolving_execElem(1),FEsolving_execElem(2) do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
myPhase = material_phase(1,i,e) ! get my phase (non-local models make no sense with more than one grain per material point) myPhase = material_phase(1,i,e) ! get my phase (non-local models make no sense with more than one grain per material point)
if (plasticState(myPhase)%nonLocal) then ! if nonlocal model if (plasticState(myPhase)%nonLocal) then ! if nonlocal model
! --- calculate disorientation between me and my neighbor --- ! --- calculate disorientation between me and my neighbor ---
do n = 1_pInt,FE_NipNeighbors(FE_celltype(FE_geomtype(mesh_element(2,e)))) ! loop through my neighbors do n = 1_pInt,FE_NipNeighbors(FE_celltype(FE_geomtype(mesh_element(2,e)))) ! loop through my neighbors
neighboring_e = mesh_ipNeighborhood(1,n,i,e) neighboring_e = mesh_ipNeighborhood(1,n,i,e)
neighboring_i = mesh_ipNeighborhood(2,n,i,e) neighboring_i = mesh_ipNeighborhood(2,n,i,e)
if (neighboring_e > 0 .and. neighboring_i > 0) then ! if neighbor exists if (neighboring_e > 0 .and. neighboring_i > 0) then ! if neighbor exists
neighboringPhase = material_phase(1,neighboring_i,neighboring_e) ! get my neighbor's phase neighboringPhase = material_phase(1,neighboring_i,neighboring_e) ! get my neighbor's phase
if (plasticState(neighboringPhase)%nonLocal) then ! neighbor got also nonlocal plasticity if (plasticState(neighboringPhase)%nonLocal) then ! neighbor got also nonlocal plasticity
if (lattice_structure(myPhase) == lattice_structure(neighboringPhase)) then ! if my neighbor has same crystal structure like me if (lattice_structure(myPhase) == lattice_structure(neighboringPhase)) then ! if my neighbor has same crystal structure like me
crystallite_disorientation(:,n,1,i,e) = & crystallite_disorientation(:,n,1,i,e) = &
lattice_qDisorientation( crystallite_orientation(1:4,1,i,e), & lattice_qDisorientation( crystallite_orientation(1:4,1,i,e), &
crystallite_orientation(1:4,1,neighboring_i,neighboring_e), & crystallite_orientation(1:4,1,neighboring_i,neighboring_e), &
lattice_structure(myPhase)) ! calculate disorientation for given symmetry lattice_structure(myPhase)) ! calculate disorientation for given symmetry
else ! for neighbor with different phase else ! for neighbor with different phase
crystallite_disorientation(:,n,1,i,e) = [0.0_pReal, 1.0_pReal, 0.0_pReal, 0.0_pReal] ! 180 degree rotation about 100 axis crystallite_disorientation(:,n,1,i,e) = [0.0_pReal, 1.0_pReal, 0.0_pReal, 0.0_pReal]! 180 degree rotation about 100 axis
endif endif
else ! for neighbor with local plasticity else ! for neighbor with local plasticity
crystallite_disorientation(:,n,1,i,e) = [-1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal] ! homomorphic identity crystallite_disorientation(:,n,1,i,e) = [-1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal]! homomorphic identity
endif endif
else ! no existing neighbor else ! no existing neighbor
crystallite_disorientation(:,n,1,i,e) = [-1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal] ! homomorphic identity crystallite_disorientation(:,n,1,i,e) = [-1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal] ! homomorphic identity
endif endif
enddo enddo
@ -4043,7 +4032,8 @@ subroutine crystallite_orientations
endif endif
enddo; enddo enddo; enddo
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
endif nonlocalPresent
end subroutine crystallite_orientations end subroutine crystallite_orientations

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

165
code/lattice.f90
View File

@ -17,13 +17,7 @@ module lattice
LATTICE_maxNslipFamily = 13_pInt, & !< max # of slip system families over lattice structures LATTICE_maxNslipFamily = 13_pInt, & !< max # of slip system families over lattice structures
LATTICE_maxNtwinFamily = 4_pInt, & !< max # of twin system families over lattice structures LATTICE_maxNtwinFamily = 4_pInt, & !< max # of twin system families over lattice structures
LATTICE_maxNtransFamily = 2_pInt, & !< max # of transformation system families over lattice structures LATTICE_maxNtransFamily = 2_pInt, & !< max # of transformation system families over lattice structures
LATTICE_maxNcleavageFamily = 3_pInt, & !< max # of transformation system families over lattice structures LATTICE_maxNcleavageFamily = 3_pInt !< max # of transformation system families over lattice structures
LATTICE_maxNslip = 52_pInt, & !< max # of slip systems over lattice structures
LATTICE_maxNtwin = 24_pInt, & !< max # of twin systems over lattice structures
LATTICE_maxNinteraction = 182_pInt, & !< max # of interaction types (in hardening matrix part)
LATTICE_maxNnonSchmid = 6_pInt, & !< max # of non schmid contributions over lattice structures
LATTICE_maxNtrans = 12_pInt, & !< max # of transformations over lattice structures
LATTICE_maxNcleavage = 9_pInt !< max # of cleavage over lattice structures
integer(pInt), allocatable, dimension(:,:), protected, public :: & integer(pInt), allocatable, dimension(:,:), protected, public :: &
lattice_NslipSystem, & !< total # of slip systems in each family lattice_NslipSystem, & !< total # of slip systems in each family
@ -80,25 +74,25 @@ module lattice
lattice_NnonSchmid !< total # of non-Schmid contributions for each structure lattice_NnonSchmid !< total # of non-Schmid contributions for each structure
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! fcc ! face centered cubic
integer(pInt), dimension(LATTICE_maxNslipFamily), parameter, public :: & integer(pInt), dimension(LATTICE_maxNslipFamily), parameter, public :: &
LATTICE_fcc_NslipSystem = int([12, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],pInt) !< total # of slip systems per family for fcc LATTICE_fcc_NslipSystem = int([12, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],pInt) !< # of slip systems per family for fcc
integer(pInt), dimension(LATTICE_maxNtwinFamily), parameter, public :: & integer(pInt), dimension(LATTICE_maxNtwinFamily), parameter, public :: &
LATTICE_fcc_NtwinSystem = int([12, 0, 0, 0],pInt) !< total # of twin systems per family for fcc LATTICE_fcc_NtwinSystem = int([12, 0, 0, 0],pInt) !< # of twin systems per family for fcc
integer(pInt), dimension(LATTICE_maxNtransFamily), parameter, public :: & integer(pInt), dimension(LATTICE_maxNtransFamily), parameter, public :: &
LATTICE_fcc_NtransSystem = int([12, 0],pInt) !< total # of transformation systems per family for fcc LATTICE_fcc_NtransSystem = int([12, 0],pInt) !< # of transformation systems per family for fcc
integer(pInt), dimension(LATTICE_maxNcleavageFamily), parameter, public :: & integer(pInt), dimension(LATTICE_maxNcleavageFamily), parameter, public :: &
LATTICE_fcc_NcleavageSystem = int([3, 4, 0],pInt) !< total # of cleavage systems per family for fcc LATTICE_fcc_NcleavageSystem = int([3, 4, 0],pInt) !< # of cleavage systems per family for fcc
integer(pInt), parameter, private :: & integer(pInt), parameter, private :: &
LATTICE_fcc_Nslip = 12_pInt, & ! sum(lattice_fcc_NslipSystem), & !< total # of slip systems for fcc LATTICE_fcc_Nslip = 12_pInt, & !sum(lattice_fcc_NslipSystem), & !< total # of slip systems for fcc
LATTICE_fcc_Ntwin = 12_pInt, & ! sum(lattice_fcc_NtwinSystem) !< total # of twin systems for fcc LATTICE_fcc_Ntwin = 12_pInt, & !sum(lattice_fcc_NtwinSystem), & !< total # of twin systems for fcc
LATTICE_fcc_NnonSchmid = 0_pInt, & !< total # of non-Schmid contributions for fcc LATTICE_fcc_NnonSchmid = 0_pInt, & !< total # of non-Schmid contributions for fcc
LATTICE_fcc_Ntrans = 12_pInt, & !< total # of transformations for fcc LATTICE_fcc_Ntrans = 12_pInt, & !sum(lattice_fcc_NtransSystem), & !< total # of transformation systems for fcc
LATTICE_fcc_Ncleavage = 7_pInt !< total # of cleavage systems for fcc LATTICE_fcc_Ncleavage = 7_pInt !sum(lattice_fcc_NcleavageSystem) !< total # of cleavage systems for fcc
real(pReal), dimension(3+3,LATTICE_fcc_Nslip), parameter, private :: & real(pReal), dimension(3+3,LATTICE_fcc_Nslip), parameter, private :: &
LATTICE_fcc_systemSlip = reshape(real([& LATTICE_fcc_systemSlip = reshape(real([&
@ -312,8 +306,8 @@ module lattice
0.0, 0.0, 1.0, 45.0 & 0.0, 0.0, 1.0, 45.0 &
],[ 4_pInt,LATTICE_fcc_Ntrans]) ],[ 4_pInt,LATTICE_fcc_Ntrans])
real(pReal), dimension(LATTICE_fcc_Ntrans,LATTICE_fcc_Ntrans), parameter, private :: & ! Matrix for projection of shear from slip system to fault-band (twin) systems real(pReal), dimension(LATTICE_fcc_Ntrans,LATTICE_fcc_Ntrans), parameter, private :: & ! Matrix for projection of shear from slip system to fault-band (twin) systems
LATTICE_fccTobcc_projectionTrans = reshape(real([& ! For ns = nt = nr LATTICE_fccTobcc_projectionTrans = reshape(real([& ! For ns = nt = nr
0, 1,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, & 0, 1,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
-1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, & -1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
1,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & 1,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
@ -363,27 +357,26 @@ module lattice
],pReal),[ 3_pInt + 3_pInt,LATTICE_fcc_Ncleavage]) ],pReal),[ 3_pInt + 3_pInt,LATTICE_fcc_Ncleavage])
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! bcc ! body centered cubic
integer(pInt), dimension(LATTICE_maxNslipFamily), parameter, public :: & integer(pInt), dimension(LATTICE_maxNslipFamily), parameter, public :: &
LATTICE_bcc_NslipSystem = int([ 12, 12, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], pInt) !< total # of slip systems per family for bcc LATTICE_bcc_NslipSystem = int([ 12, 12, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], pInt) !< # of slip systems per family for bcc
integer(pInt), dimension(LATTICE_maxNtwinFamily), parameter, public :: & integer(pInt), dimension(LATTICE_maxNtwinFamily), parameter, public :: &
LATTICE_bcc_NtwinSystem = int([ 12, 0, 0, 0], pInt) !< total # of twin systems per family for bcc LATTICE_bcc_NtwinSystem = int([ 12, 0, 0, 0], pInt) !< # of twin systems per family for bcc
integer(pInt), dimension(LATTICE_maxNtransFamily), parameter, public :: & integer(pInt), dimension(LATTICE_maxNtransFamily), parameter, public :: &
LATTICE_bcc_NtransSystem = int([0,0],pInt) !< total # of transformation systems per family for bcc LATTICE_bcc_NtransSystem = int([0,0],pInt) !< # of transformation systems per family for bcc
integer(pInt), dimension(LATTICE_maxNcleavageFamily), parameter, public :: & integer(pInt), dimension(LATTICE_maxNcleavageFamily), parameter, public :: &
LATTICE_bcc_NcleavageSystem = int([3,6,0],pInt) !< total # of cleavage systems per family for bcc LATTICE_bcc_NcleavageSystem = int([3,6,0],pInt) !< # of cleavage systems per family for bcc
integer(pInt), parameter, private :: & integer(pInt), parameter, private :: &
LATTICE_bcc_Nslip = 24_pInt, & ! sum(lattice_bcc_NslipSystem), & !< total # of slip systems for bcc LATTICE_bcc_Nslip = 24_pInt, & !sum(lattice_bcc_NslipSystem), & !< total # of slip systems for bcc
LATTICE_bcc_Ntwin = 12_pInt, & ! sum(lattice_bcc_NtwinSystem) !< total # of twin systems for bcc LATTICE_bcc_Ntwin = 12_pInt, & !sum(lattice_bcc_NtwinSystem), & !< total # of twin systems for bcc
LATTICE_bcc_NnonSchmid = 6_pInt, & !< # of non-Schmid contributions for bcc. 6 known non schmid contributions for BCC (A. Koester, A. Ma, A. Hartmaier 2012) LATTICE_bcc_NnonSchmid = 6_pInt, & !< total # of non-Schmid contributions for bcc (A. Koester, A. Ma, A. Hartmaier 2012)
LATTICE_bcc_Ntrans = 0_pInt, & !< total # of transformations for bcc LATTICE_bcc_Ntrans = 0_pInt, & !sum(lattice_bcc_NtransSystem), & !< total # of transformation systems for bcc
LATTICE_bcc_Ncleavage = 9_pInt !< total # of cleavage systems for bcc LATTICE_bcc_Ncleavage = 9_pInt !sum(lattice_bcc_NcleavageSystem) !< total # of cleavage systems for bcc
real(pReal), dimension(3+3,LATTICE_bcc_Nslip), parameter, private :: & real(pReal), dimension(3+3,LATTICE_bcc_Nslip), parameter, private :: &
LATTICE_bcc_systemSlip = reshape(real([& LATTICE_bcc_systemSlip = reshape(real([&
! Slip direction Plane normal ! Slip direction Plane normal
@ -561,25 +554,25 @@ module lattice
],pReal),[ 3_pInt + 3_pInt,LATTICE_bcc_Ncleavage]) ],pReal),[ 3_pInt + 3_pInt,LATTICE_bcc_Ncleavage])
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! hex ! hexagonal
integer(pInt), dimension(LATTICE_maxNslipFamily), parameter, public :: & integer(pInt), dimension(LATTICE_maxNslipFamily), parameter, public :: &
lattice_hex_NslipSystem = int([ 3, 3, 3, 6, 12, 6, 0, 0, 0, 0, 0, 0, 0],pInt) !< # of slip systems per family for hex lattice_hex_NslipSystem = int([ 3, 3, 3, 6, 12, 6, 0, 0, 0, 0, 0, 0, 0],pInt) !< # of slip systems per family for hex
integer(pInt), dimension(LATTICE_maxNtwinFamily), parameter, public :: & integer(pInt), dimension(LATTICE_maxNtwinFamily), parameter, public :: &
lattice_hex_NtwinSystem = int([ 6, 6, 6, 6],pInt) !< # of slip systems per family for hex lattice_hex_NtwinSystem = int([ 6, 6, 6, 6],pInt) !< # of slip systems per family for hex
integer(pInt), dimension(LATTICE_maxNtransFamily), parameter, public :: & integer(pInt), dimension(LATTICE_maxNtransFamily), parameter, public :: &
LATTICE_hex_NtransSystem = int([0,0],pInt) !< total # of transformation systems per family for hex LATTICE_hex_NtransSystem = int([0,0],pInt) !< # of transformation systems per family for hex
integer(pInt), dimension(LATTICE_maxNcleavageFamily), parameter, public :: & integer(pInt), dimension(LATTICE_maxNcleavageFamily), parameter, public :: &
LATTICE_hex_NcleavageSystem = int([3,0,0],pInt) !< total # of cleavage systems per family for hex LATTICE_hex_NcleavageSystem = int([3,0,0],pInt) !< # of cleavage systems per family for hex
integer(pInt), parameter , private :: & integer(pInt), parameter, private :: &
LATTICE_hex_Nslip = 33_pInt, & ! sum(lattice_hex_NslipSystem), !< total # of slip systems for hex LATTICE_hex_Nslip = 33_pInt, & !sum(lattice_hex_NslipSystem), & !< total # of slip systems for hex
LATTICE_hex_Ntwin = 24_pInt, & ! sum(lattice_hex_NtwinSystem) !< total # of twin systems for hex LATTICE_hex_Ntwin = 24_pInt, & !sum(lattice_hex_NtwinSystem), & !< total # of twin systems for hex
LATTICE_hex_NnonSchmid = 0_pInt, & !< # of non-Schmid contributions for hex LATTICE_hex_NnonSchmid = 0_pInt, & !< total # of non-Schmid contributions for hex
LATTICE_hex_Ntrans = 0_pInt, & !< total # of transformations for hex LATTICE_hex_Ntrans = 0_pInt, & !sum(lattice_hex_NtransSystem), & !< total # of transformation systems for hex
LATTICE_hex_Ncleavage = 3_pInt !< total # of transformations for hex LATTICE_hex_Ncleavage = 3_pInt !sum(lattice_hex_NcleavageSystem) !< total # of cleavage systems for hex
real(pReal), dimension(4+4,LATTICE_hex_Nslip), parameter, private :: & real(pReal), dimension(4+4,LATTICE_hex_Nslip), parameter, private :: &
LATTICE_hex_systemSlip = reshape(real([& LATTICE_hex_systemSlip = reshape(real([&
@ -842,28 +835,26 @@ module lattice
],pReal),[ 4_pInt + 4_pInt,LATTICE_hex_Ncleavage]) ],pReal),[ 4_pInt + 4_pInt,LATTICE_hex_Ncleavage])
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! bct ! body centered tetragonal
integer(pInt), dimension(LATTICE_maxNslipFamily), parameter, public :: & integer(pInt), dimension(LATTICE_maxNslipFamily), parameter, public :: &
LATTICE_bct_NslipSystem = int([2, 2, 2, 4, 2, 4, 2, 2, 4, 8, 4, 8, 8 ],pInt) !< # of slip systems per family for bct (Sn) Bieler J. Electr Mater 2009 LATTICE_bct_NslipSystem = int([2, 2, 2, 4, 2, 4, 2, 2, 4, 8, 4, 8, 8 ],pInt) !< # of slip systems per family for bct (Sn) Bieler J. Electr Mater 2009
integer(pInt), dimension(LATTICE_maxNtwinFamily), parameter, public :: & integer(pInt), dimension(LATTICE_maxNtwinFamily), parameter, public :: &
LATTICE_bct_NtwinSystem = int([0, 0, 0, 0], pInt) !< total # of twin systems per family for bct-example LATTICE_bct_NtwinSystem = int([0, 0, 0, 0], pInt) !< # of twin systems per family for bct
integer(pInt), dimension(LATTICE_maxNtransFamily), parameter, public :: & integer(pInt), dimension(LATTICE_maxNtransFamily), parameter, public :: &
LATTICE_bct_NtransSystem = int([0,0],pInt) !< total # of transformation systems per family for bct LATTICE_bct_NtransSystem = int([0,0],pInt) !< # of transformation systems per family for bct
integer(pInt), dimension(LATTICE_maxNcleavageFamily), parameter, public :: & integer(pInt), dimension(LATTICE_maxNcleavageFamily), parameter, public :: &
LATTICE_bct_NcleavageSystem = int([0,0,0],pInt) !< total # of cleavage systems per family for bct LATTICE_bct_NcleavageSystem = int([0,0,0],pInt) !< # of cleavage systems per family for bct
integer(pInt), parameter , private :: & integer(pInt), parameter, private :: &
LATTICE_bct_Nslip = 52_pInt, & ! sum(lattice_bct_NslipSystem), !< total # of slip systems for bct LATTICE_bct_Nslip = 52_pInt, & !sum(lattice_bct_NslipSystem), & !< total # of slip systems for bct
LATTICE_bct_Ntwin = 0_pInt, & ! sum(lattice_bcc_NtwinSystem) !< total # of twin systems for bct LATTICE_bct_Ntwin = 0_pInt, & !sum(lattice_bct_NtwinSystem), & !< total # of twin systems for bct
LATTICE_bct_NnonSchmid = 0_pInt, & !< # of non-Schmid contributions for bct LATTICE_bct_NnonSchmid = 0_pInt, & !< total # of non-Schmid contributions for bct
LATTICE_bct_Ntrans = 0_pInt, & !< total # of transformations for bct LATTICE_bct_Ntrans = 0_pInt, & !sum(lattice_bct_NtransSystem), & !< total # of transformation systems for bct
LATTICE_bct_Ncleavage = 0_pInt !< total # of transformations for bct LATTICE_bct_Ncleavage = 0_pInt !sum(lattice_bct_NcleavageSystem) !< total # of cleavage systems for bct
real(pReal), dimension(3+3,LATTICE_bct_Nslip), parameter, private :: & real(pReal), dimension(3+3,LATTICE_bct_Nslip), parameter, private :: &
LATTICE_bct_systemSlip = reshape(real([& LATTICE_bct_systemSlip = reshape(real([&
@ -1006,12 +997,25 @@ module lattice
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! isotropic ! isotropic
integer(pInt), dimension(LATTICE_maxNslipFamily), parameter, public :: &
LATTICE_iso_NslipSystem = int([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ],pInt) !< # of slip systems per family for iso
integer(pInt), dimension(LATTICE_maxNtwinFamily), parameter, public :: &
LATTICE_iso_NtwinSystem = int([0, 0, 0, 0], pInt) !< # of twin systems per family for iso
integer(pInt), dimension(LATTICE_maxNtransFamily), parameter, public :: &
LATTICE_iso_NtransSystem = int([0, 0],pInt) !< # of transformation systems per family for iso
integer(pInt), dimension(LATTICE_maxNcleavageFamily), parameter, public :: & integer(pInt), dimension(LATTICE_maxNcleavageFamily), parameter, public :: &
LATTICE_iso_NcleavageSystem = int([3,0,0],pInt) !< total # of cleavage systems per family for isotropic LATTICE_iso_NcleavageSystem = int([3,0,0],pInt) !< # of cleavage systems per family for iso
integer(pInt), parameter, private :: & integer(pInt), parameter, private :: &
LATTICE_iso_Ncleavage = 3_pInt !< total # of cleavage systems for bcc LATTICE_iso_Nslip = 0_pInt, & !sum(lattice_iso_NslipSystem), & !< total # of slip systems for iso
LATTICE_iso_Ntwin = 0_pInt, & !sum(lattice_iso_NtwinSystem), & !< total # of twin systems for iso
LATTICE_iso_NnonSchmid = 0_pInt, & !< total # of non-Schmid contributions for iso
LATTICE_iso_Ntrans = 0_pInt, & !sum(lattice_iso_NtransSystem), & !< total # of transformation systems for iso
LATTICE_iso_Ncleavage = 3_pInt !sum(lattice_iso_NcleavageSystem) !< total # of cleavage systems for iso
real(pReal), dimension(3+3,LATTICE_iso_Ncleavage), parameter, private :: & real(pReal), dimension(3+3,LATTICE_iso_Ncleavage), parameter, private :: &
LATTICE_iso_systemCleavage = reshape(real([& LATTICE_iso_systemCleavage = reshape(real([&
! Cleavage direction Plane normal ! Cleavage direction Plane normal
@ -1022,12 +1026,25 @@ module lattice
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! orthorhombic ! orthorhombic
integer(pInt), dimension(LATTICE_maxNslipFamily), parameter, public :: &
LATTICE_ortho_NslipSystem = int([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ],pInt) !< # of slip systems per family for ortho
integer(pInt), dimension(LATTICE_maxNtwinFamily), parameter, public :: &
LATTICE_ortho_NtwinSystem = int([0, 0, 0, 0], pInt) !< # of twin systems per family for ortho
integer(pInt), dimension(LATTICE_maxNtransFamily), parameter, public :: &
LATTICE_ortho_NtransSystem = int([0, 0],pInt) !< # of transformation systems per family for ortho
integer(pInt), dimension(LATTICE_maxNcleavageFamily), parameter, public :: & integer(pInt), dimension(LATTICE_maxNcleavageFamily), parameter, public :: &
LATTICE_ortho_NcleavageSystem = int([1,1,1],pInt) !< total # of cleavage systems per family for orthotropic LATTICE_ortho_NcleavageSystem = int([1,1,1],pInt) !< # of cleavage systems per family for ortho
integer(pInt), parameter, private :: & integer(pInt), parameter, private :: &
LATTICE_ortho_Ncleavage = 3_pInt !< total # of cleavage systems for bcc LATTICE_ortho_Nslip = 0_pInt, & !sum(lattice_ortho_NslipSystem), & !< total # of slip systems for ortho
LATTICE_ortho_Ntwin = 0_pInt, & !sum(lattice_ortho_NtwinSystem), & !< total # of twin systems for ortho
LATTICE_ortho_NnonSchmid = 0_pInt, & !< total # of non-Schmid contributions for ortho
LATTICE_ortho_Ntrans = 0_pInt, & !sum(lattice_ortho_NtransSystem), & !< total # of transformation systems for ortho
LATTICE_ortho_Ncleavage = 3_pInt !sum(lattice_ortho_NcleavageSystem) !< total # of cleavage systems for ortho
real(pReal), dimension(3+3,LATTICE_ortho_Ncleavage), parameter, private :: & real(pReal), dimension(3+3,LATTICE_ortho_Ncleavage), parameter, private :: &
LATTICE_ortho_systemCleavage = reshape(real([& LATTICE_ortho_systemCleavage = reshape(real([&
! Cleavage direction Plane normal ! Cleavage direction Plane normal
@ -1036,16 +1053,36 @@ module lattice
1, 0, 0, 0, 0, 1 & 1, 0, 0, 0, 0, 1 &
],pReal),[ 3_pInt + 3_pInt,LATTICE_ortho_Ncleavage]) ],pReal),[ 3_pInt + 3_pInt,LATTICE_ortho_Ncleavage])
real(pReal), dimension(:,:,:), allocatable, public, protected :: & integer(pInt), parameter, public :: &
LATTICE_maxNslip = 52_pInt, &
!LATTICE_maxNslip = maxval([LATTICE_fcc_Nslip,LATTICE_bcc_Nslip,LATTICE_hex_Nslip,\
! LATTICE_bct_Nslip,LATTICE_iso_Nslip,LATTICE_ortho_Nslip]), & !< max # of slip systems over lattice structures
LATTICE_maxNtwin = 24_pInt, &
!LATTICE_maxNtwin = maxval([LATTICE_fcc_Ntwin,LATTICE_bcc_Ntwin,LATTICE_hex_Ntwin,\
! LATTICE_bct_Ntwin,LATTICE_iso_Ntwin,LATTICE_ortho_Ntwin]), & !< max # of twin systems over lattice structures
LATTICE_maxNnonSchmid = 6_pInt, &
!LATTICE_maxNtwin = maxval([LATTICE_fcc_NnonSchmid,LATTICE_bcc_NnonSchmid,\
! LATTICE_hex_NnonSchmid,LATTICE_bct_NnonSchmid,\
! LATTICE_iso_NnonSchmid,LATTICE_ortho_NnonSchmid]), & !< max # of non-Schmid contributions over lattice structures
LATTICE_maxNtrans = 12_pInt, &
!LATTICE_maxNtrans = maxval([LATTICE_fcc_Ntrans,LATTICE_bcc_Ntrans,LATTICE_hex_Ntrans,\
! LATTICE_bct_Ntrans,LATTICE_iso_Ntrans,LATTICE_ortho_Ntrans]),&!< max # of transformation systems over lattice structures
LATTICE_maxNcleavage = 9_pInt, &
!LATTICE_maxNcleavage = maxval([LATTICE_fcc_Ncleavage,LATTICE_bcc_Ncleavage,\
! LATTICE_hex_Ncleavage,LATTICE_bct_Ncleavage,\
! LATTICE_iso_Ncleavage,LATTICE_ortho_Ncleavage]) !< max # of cleavage systems over lattice structures
LATTICE_maxNinteraction = 182_pInt !< max # of interaction types (in hardening matrix part)
real(pReal), dimension(:,:,:), allocatable, public, protected :: &
lattice_C66, lattice_trans_C66 lattice_C66, lattice_trans_C66
real(pReal), dimension(:,:,:,:,:), allocatable, public, protected :: & real(pReal), dimension(:,:,:,:,:), allocatable, public, protected :: &
lattice_C3333, lattice_trans_C3333 lattice_C3333, lattice_trans_C3333
real(pReal), dimension(:), allocatable, public, protected :: & real(pReal), dimension(:), allocatable, public, protected :: &
lattice_mu, & lattice_mu, &
lattice_nu, & lattice_nu, &
lattice_trans_mu, & lattice_trans_mu, &
lattice_trans_nu lattice_trans_nu
real(pReal), dimension(:,:,:), allocatable, public, protected :: & real(pReal), dimension(:,:,:), allocatable, public, protected :: &
lattice_thermalConductivity33, & lattice_thermalConductivity33, &
lattice_thermalExpansion33, & lattice_thermalExpansion33, &
lattice_damageDiffusion33, & lattice_damageDiffusion33, &
@ -1054,7 +1091,7 @@ module lattice
lattice_porosityDiffusion33, & lattice_porosityDiffusion33, &
lattice_hydrogenfluxDiffusion33, & lattice_hydrogenfluxDiffusion33, &
lattice_hydrogenfluxMobility33 lattice_hydrogenfluxMobility33
real(pReal), dimension(:), allocatable, public, protected :: & real(pReal), dimension(:), allocatable, public, protected :: &
lattice_damageMobility, & lattice_damageMobility, &
lattice_porosityMobility, & lattice_porosityMobility, &
lattice_massDensity, & lattice_massDensity, &
@ -1253,7 +1290,7 @@ subroutine lattice_init
endif mainProcess endif mainProcess
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! consistency checks ! consistency checks (required since ifort 15.0 does not support sum/maxval in parameter definition)
if (LATTICE_maxNslip /= maxval([LATTICE_fcc_Nslip,LATTICE_bcc_Nslip,LATTICE_hex_Nslip,LATTICE_bct_Nslip])) & if (LATTICE_maxNslip /= maxval([LATTICE_fcc_Nslip,LATTICE_bcc_Nslip,LATTICE_hex_Nslip,LATTICE_bct_Nslip])) &
call IO_error(0_pInt,ext_msg = 'LATTICE_maxNslip') call IO_error(0_pInt,ext_msg = 'LATTICE_maxNslip')
@ -2182,7 +2219,7 @@ pure function lattice_qDisorientation(Q1, Q2, struct)
real(pReal), dimension(4) :: lattice_qDisorientation real(pReal), dimension(4) :: lattice_qDisorientation
real(pReal), dimension(4), intent(in) :: & real(pReal), dimension(4), intent(in) :: &
Q1, & ! 1st orientation Q1, & ! 1st orientation
Q2 ! 2nd orientation Q2 ! 2nd orientation
integer(kind(LATTICE_undefined_ID)), optional, intent(in) :: & ! if given, symmetries between the two orientation will be considered integer(kind(LATTICE_undefined_ID)), optional, intent(in) :: & ! if given, symmetries between the two orientation will be considered
struct struct

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

42
code/math.f90
View File

@ -186,10 +186,6 @@ module math
halton_seed_set, & halton_seed_set, &
i_to_halton, & i_to_halton, &
prime prime
external :: &
dsyev, &
dgetrf, &
dgetri
contains contains
@ -811,15 +807,13 @@ function math_invSym3333(A)
integer(pInt), dimension(6) :: ipiv6 integer(pInt), dimension(6) :: ipiv6
real(pReal), dimension(6,6) :: temp66_Real real(pReal), dimension(6,6) :: temp66_Real
real(pReal), dimension(6) :: work6 real(pReal), dimension(6) :: work6
external :: &
dgetrf, &
dgetri
temp66_real = math_Mandel3333to66(A) temp66_real = math_Mandel3333to66(A)
#if(FLOAT==8)
call dgetrf(6,6,temp66_real,6,ipiv6,ierr) call dgetrf(6,6,temp66_real,6,ipiv6,ierr)
call dgetri(6,temp66_real,6,ipiv6,work6,6,ierr) call dgetri(6,temp66_real,6,ipiv6,work6,6,ierr)
#elif(FLOAT==4)
call sgetrf(6,6,temp66_real,6,ipiv6,ierr)
call sgetri(6,temp66_real,6,ipiv6,work6,6,ierr)
#endif
if (ierr == 0_pInt) then if (ierr == 0_pInt) then
math_invSym3333 = math_Mandel66to3333(temp66_real) math_invSym3333 = math_Mandel66to3333(temp66_real)
else else
@ -845,15 +839,13 @@ subroutine math_invert(myDim,A, InvA, error)
real(pReal), dimension(myDim,myDim), intent(out) :: invA real(pReal), dimension(myDim,myDim), intent(out) :: invA
logical, intent(out) :: error logical, intent(out) :: error
external :: &
dgetrf, &
dgetri
invA = A invA = A
#if(FLOAT==8)
call dgetrf(myDim,myDim,invA,myDim,ipiv,ierr) call dgetrf(myDim,myDim,invA,myDim,ipiv,ierr)
call dgetri(myDim,InvA,myDim,ipiv,work,myDim,ierr) call dgetri(myDim,InvA,myDim,ipiv,work,myDim,ierr)
#elif(FLOAT==4)
call sgetrf(myDim,myDim,invA,myDim,ipiv,ierr)
call sgetri(myDim,InvA,myDim,ipiv,work,myDim,ierr)
#endif
error = merge(.true.,.false., ierr /= 0_pInt) ! http://fortraninacworld.blogspot.de/2012/12/ternary-operator.html error = merge(.true.,.false., ierr /= 0_pInt) ! http://fortraninacworld.blogspot.de/2012/12/ternary-operator.html
end subroutine math_invert end subroutine math_invert
@ -1656,14 +1648,14 @@ pure function math_qToAxisAngle(Q)
real(pReal) :: halfAngle, sinHalfAngle real(pReal) :: halfAngle, sinHalfAngle
real(pReal), dimension(4) :: math_qToAxisAngle real(pReal), dimension(4) :: math_qToAxisAngle
halfAngle = acos(max(-1.0_pReal, min(1.0_pReal, Q(1)))) ! limit to [-1,1] --> 0 to 180 deg halfAngle = acos(math_limit(Q(1),-1.0_pReal,1.0_pReal))
sinHalfAngle = sin(halfAngle) sinHalfAngle = sin(halfAngle)
if (sinHalfAngle <= 1.0e-4_pReal) then ! very small rotation angle? smallRotation: if (sinHalfAngle <= 1.0e-4_pReal) then
math_qToAxisAngle = 0.0_pReal math_qToAxisAngle = 0.0_pReal
else else smallRotation
math_qToAxisAngle= [ Q(2:4)/sinHalfAngle, halfAngle*2.0_pReal] math_qToAxisAngle= [ Q(2:4)/sinHalfAngle, halfAngle*2.0_pReal]
endif endif smallRotation
end function math_qToAxisAngle end function math_qToAxisAngle
@ -1937,16 +1929,13 @@ subroutine math_eigenValuesVectorsSym(m,values,vectors,error)
real(pReal), dimension(size(m,1)), intent(out) :: values real(pReal), dimension(size(m,1)), intent(out) :: values
real(pReal), dimension(size(m,1),size(m,1)), intent(out) :: vectors real(pReal), dimension(size(m,1),size(m,1)), intent(out) :: vectors
logical, intent(out) :: error logical, intent(out) :: error
integer(pInt) :: info integer(pInt) :: info
real(pReal), dimension((64+2)*size(m,1)) :: work ! block size of 64 taken from http://www.netlib.org/lapack/double/dsyev.f real(pReal), dimension((64+2)*size(m,1)) :: work ! block size of 64 taken from http://www.netlib.org/lapack/double/dsyev.f
external :: &
dsyev
vectors = m ! copy matrix to input (doubles as output) array vectors = m ! copy matrix to input (doubles as output) array
#if(FLOAT==8)
call dsyev('V','U',size(m,1),vectors,size(m,1),values,work,(64+2)*size(m,1),info) call dsyev('V','U',size(m,1),vectors,size(m,1),values,work,(64+2)*size(m,1),info)
#elif(FLOAT==4)
call ssyev('V','U',size(m,1),vectors,size(m,1),values,work,(64+2)*size(m,1),info)
#endif
error = (info == 0_pInt) error = (info == 0_pInt)
end subroutine math_eigenValuesVectorsSym end subroutine math_eigenValuesVectorsSym
@ -2135,16 +2124,13 @@ function math_eigenvaluesSym(m)
real(pReal), dimension(:,:), intent(in) :: m real(pReal), dimension(:,:), intent(in) :: m
real(pReal), dimension(size(m,1)) :: math_eigenvaluesSym real(pReal), dimension(size(m,1)) :: math_eigenvaluesSym
real(pReal), dimension(size(m,1),size(m,1)) :: vectors real(pReal), dimension(size(m,1),size(m,1)) :: vectors
integer(pInt) :: info integer(pInt) :: info
real(pReal), dimension((64+2)*size(m,1)) :: work ! block size of 64 taken from http://www.netlib.org/lapack/double/dsyev.f real(pReal), dimension((64+2)*size(m,1)) :: work ! block size of 64 taken from http://www.netlib.org/lapack/double/dsyev.f
external :: &
dsyev
vectors = m ! copy matrix to input (doubles as output) array vectors = m ! copy matrix to input (doubles as output) array
#if(FLOAT==8)
call dsyev('N','U',size(m,1),vectors,size(m,1),math_eigenvaluesSym,work,(64+2)*size(m,1),info) call dsyev('N','U',size(m,1),vectors,size(m,1),math_eigenvaluesSym,work,(64+2)*size(m,1),info)
#elif(FLOAT==4)
call ssyev('N','U',size(m,1),vectors,size(m,1),math_eigenvaluesSym,work,(64+2)*size(m,1),info)
#endif
if (info /= 0_pInt) math_eigenvaluesSym = DAMASK_NaN if (info /= 0_pInt) math_eigenvaluesSym = DAMASK_NaN
end function math_eigenvaluesSym end function math_eigenvaluesSym

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', &

823
code/plastic_disloUCLA.f90
View File

File diff suppressed because it is too large Load Diff

104
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,22 +36,22 @@ 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, & 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, & tausat_SinhFitA= 0.0_pReal, &
tausat_SinhFitB, & tausat_SinhFitB= 0.0_pReal, &
tausat_SinhFitC, & tausat_SinhFitC= 0.0_pReal, &
tausat_SinhFitD tausat_SinhFitD= 0.0_pReal
logical :: & logical :: &
dilatation dilatation = .false.
end type end type
type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance) type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance)
@ -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)
@ -466,15 +462,15 @@ subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dTstar_3333,Tstar_v,ipc,ip,e
math_spherical33, & math_spherical33, &
math_mul33xx33 math_mul33xx33
use material, only: & use material, only: &
phaseAt, 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,9 +480,7 @@ 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) :: & real(pReal) :: &
dLi_dTstar_3333 !< derivative of Li with respect to Tstar as 4th order tensor
real(pReal) :: &
gamma_dot, & !< strainrate gamma_dot, & !< strainrate
norm_Tstar_sph, & !< euclidean norm of Tstar_sph norm_Tstar_sph, & !< euclidean norm of Tstar_sph
squarenorm_Tstar_sph !< square of the euclidean norm of Tstar_sph squarenorm_Tstar_sph !< square of the euclidean norm of Tstar_sph
@ -495,34 +489,32 @@ subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dTstar_3333,Tstar_v,ipc,ip,e
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_sph_33 = math_spherical33(math_Mandel6to33(Tstar_v)) ! spherical part of 2nd Piola-Kirchhoff stress Tstar_sph_33 = math_spherical33(math_Mandel6to33(Tstar_v)) ! spherical part of 2nd Piola-Kirchhoff stress
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 gamma_dot = param(instance)%gdot0 &
Li = 0.0_pReal * (sqrt(1.5_pReal) * norm_Tstar_sph / param(instance)%fTaylor / state(instance)%flowstress(of) ) &
dLi_dTstar_3333 = 0.0_pReal **param(instance)%n
else
gamma_dot = param(instance)%gdot0 &
* (sqrt(1.5_pReal) * norm_Tstar_sph / param(instance)%fTaylor / state(instance)%flowstress(of) ) &
**param(instance)%n
Li = Tstar_sph_33/norm_Tstar_sph * gamma_dot/param(instance)%fTaylor Li = Tstar_sph_33/norm_Tstar_sph * gamma_dot/param(instance)%fTaylor
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! Calculation of the tangent of Li ! Calculation of the tangent of Li
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
dLi_dTstar_3333(k,l,m,n) = (param(instance)%n-1.0_pReal) * & dLi_dTstar_3333(k,l,m,n) = (param(instance)%n-1.0_pReal) * &
Tstar_sph_33(k,l)*Tstar_sph_33(m,n) / squarenorm_Tstar_sph Tstar_sph_33(k,l)*Tstar_sph_33(m,n) / squarenorm_Tstar_sph
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt) & forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt) &
dLi_dTstar_3333(k,l,k,l) = dLi_dTstar_3333(k,l,k,l) + 1.0_pReal dLi_dTstar_3333(k,l,k,l) = dLi_dTstar_3333(k,l,k,l) + 1.0_pReal
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
@ -535,8 +527,7 @@ subroutine plastic_isotropic_dotState(Tstar_v,ipc,ip,el)
use math, only: & use math, only: &
math_mul6x6 math_mul6x6
use material, only: & use material, only: &
phaseAt, phasememberAt, & phasememberAt, &
plasticState, &
material_phase, & material_phase, &
phase_plasticityInstance phase_plasticityInstance
@ -559,7 +550,7 @@ subroutine plastic_isotropic_dotState(Tstar_v,ipc,ip,el)
of !< shortcut notation for offset position in state array of !< shortcut notation for offset position in state array
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))
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! norm of (deviatoric) 2nd Piola-Kirchhoff stress ! norm of (deviatoric) 2nd Piola-Kirchhoff stress
@ -615,8 +606,7 @@ 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, &
phaseAt, phasememberAt, &
phase_plasticityInstance phase_plasticityInstance
implicit none implicit none
@ -640,7 +630,7 @@ function plastic_isotropic_postResults(Tstar_v,ipc,ip,el)
o o
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))
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! norm of (deviatoric) 2nd Piola-Kirchhoff stress ! norm of (deviatoric) 2nd Piola-Kirchhoff stress

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

51
code/prec.f90
View File

@ -4,9 +4,9 @@
!> @author Christoph Kords, Max-Planck-Institut für Eisenforschung GmbH !> @author Christoph Kords, Max-Planck-Institut für Eisenforschung GmbH
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @author Luv Sharma, Max-Planck-Institut für Eisenforschung GmbH !> @author Luv Sharma, Max-Planck-Institut für Eisenforschung GmbH
!> @brief setting precision for real and int type depending on makros "FLOAT" and "INT" !> @brief setting precision for real and int type
!> @details setting precision for real and int type and for DAMASK_NaN. Definition is made !> @details setting precision for real and int type and for DAMASK_NaN. Definition is made
!! depending on makros "FLOAT" and "INT" defined during compilation !! depending on makro "INT" defined during compilation
!! for details on NaN see https://software.intel.com/en-us/forums/topic/294680 !! for details on NaN see https://software.intel.com/en-us/forums/topic/294680
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module prec module prec
@ -18,18 +18,7 @@ module prec
implicit none implicit none
private private
#if (FLOAT==4) #if (FLOAT==8)
#if defined(Spectral) || defined(FEM)
SPECTRAL SOLVER AND OWN FEM DO NOT SUPPORT SINGLE PRECISION, STOPPING COMPILATION
#endif
integer, parameter, public :: pReal = 4 !< floating point single precition (was selected_real_kind(6,37), number with 6 significant digits, up to 1e+-37)
#ifdef __INTEL_COMPILER
real(pReal), parameter, public :: DAMASK_NaN = Z'7F800001' !< quiet NaN for single precision (from http://www.hpc.unimelb.edu.au/doc/f90lrm/dfum_035.html, copy can be found in documentation/Code/Fortran)
#endif
#ifdef __GFORTRAN__
real(pReal), parameter, public :: DAMASK_NaN = real(Z'7F800001', pReal) !< quiet NaN for single precision (from http://www.hpc.unimelb.edu.au/doc/f90lrm/dfum_035.html, copy can be found in documentation/Code/Fortran)
#endif
#elif (FLOAT==8)
integer, parameter, public :: pReal = 8 !< floating point double precision (was selected_real_kind(15,300), number with 15 significant digits, up to 1e+-300) integer, parameter, public :: pReal = 8 !< floating point double precision (was selected_real_kind(15,300), number with 15 significant digits, up to 1e+-300)
#ifdef __INTEL_COMPILER #ifdef __INTEL_COMPILER
real(pReal), parameter, public :: DAMASK_NaN = Z'7FF8000000000000' !< quiet NaN for double precision (from http://www.hpc.unimelb.edu.au/doc/f90lrm/dfum_035.html, copy can be found in documentation/Code/Fortran) real(pReal), parameter, public :: DAMASK_NaN = Z'7FF8000000000000' !< quiet NaN for double precision (from http://www.hpc.unimelb.edu.au/doc/f90lrm/dfum_035.html, copy can be found in documentation/Code/Fortran)
@ -124,7 +113,9 @@ module prec
public :: & public :: &
prec_init, & prec_init, &
prec_isNaN prec_isNaN, &
dEq, &
dNeq
contains contains
@ -172,9 +163,9 @@ end subroutine prec_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief figures out if a floating point number is NaN !> @brief figures out if a floating point number is NaN
! basically just a small wrapper, because gfortran < 4.9 does not have the IEEE module ! basically just a small wrapper, because gfortran < 5.0 does not have the IEEE module
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
logical elemental function prec_isNaN(a) logical elemental pure function prec_isNaN(a)
implicit none implicit none
real(pReal), intent(in) :: a real(pReal), intent(in) :: a
@ -187,4 +178,30 @@ logical elemental function prec_isNaN(a)
#endif #endif
end function prec_isNaN end function prec_isNaN
!--------------------------------------------------------------------------------------------------
!> @brief equality comparison for double precision
! replaces "==" but for certain (relative) tolerance. Counterpart to dNeq
! http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm
!--------------------------------------------------------------------------------------------------
logical elemental pure function dEq(a,b,tol)
real(pReal), intent(in) :: a,b
real(pReal), intent(in), optional :: tol
real(pReal), parameter :: eps = 2.2204460492503131E-16 ! DBL_EPSILON in C
dEq = merge(.True., .False.,abs(a-b) <= merge(tol,eps,present(tol))*maxval(abs([a,b])))
end function dEq
!--------------------------------------------------------------------------------------------------
!> @brief inequality comparison for double precision
! replaces "!=" but for certain (relative) tolerance. Counterpart to dEq
! http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm
!--------------------------------------------------------------------------------------------------
logical elemental pure function dNeq(a,b,tol)
real(pReal), intent(in) :: a,b
real(pReal), intent(in), optional :: tol
real(pReal), parameter :: eps = 2.2204460492503131E-16 ! DBL_EPSILON in C
dNeq = merge(.False., .True.,abs(a-b) <= merge(tol,eps,present(tol))*maxval(abs([a,b])))
end function dNeq
end module prec end module prec

1
config/Crystallite_All.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
[all] [all]
(output) phase (output) phase
(output) texture (output) texture

1
config/Crystallite_None.config
View File

@ -1,2 +1 @@
### $Id$ ###
[none] [none]

1
config/Crystallite_aLittleSomething.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
[aLittleSomething] [aLittleSomething]
(output) f (output) f
(output) p (output) p

1
config/Homogenization_Damage_NonLocal.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
damage nonlocal damage nonlocal
initialDamage 1.0 initialDamage 1.0
(output) damage (output) damage

1
config/Homogenization_HydrogenFlux_CahnHilliard.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
hydrogenflux cahnhilliard hydrogenflux cahnhilliard
initialHydrogenConc 0.0 initialHydrogenConc 0.0
(output) hydrogenconc (output) hydrogenconc

1
config/Homogenization_Isostrain_Parallel3.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
[Parallel3] [Parallel3]
type isostrain type isostrain
Ngrains 3 Ngrains 3

1
config/Homogenization_Isostrain_SX.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
[SX] [SX]
type isostrain type isostrain
Ngrains 1 Ngrains 1

1
config/Homogenization_Isostrain_Taylor2.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
[Taylor2] [Taylor2]
type isostrain type isostrain
Ngrains 2 Ngrains 2

1
config/Homogenization_None_Dummy.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
[directSX] [directSX]
type none type none

1
config/Homogenization_Porosity_PhaseField.config
View File

@ -1,3 +1,2 @@
### $Id$ ###
porosity phasefield porosity phasefield
(output) porosity (output) porosity

1
config/Homogenization_RGC_8Grains.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
[8Grains] [8Grains]
type RGC type RGC
Ngrains 8 Ngrains 8

1
config/Homogenization_Thermal_Conduction.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
thermal conduction thermal conduction
initialT 300.0 initialT 300.0
(output) temperature (output) temperature

1
config/Homogenization_VacancyFlux_CahnHilliard.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
vacancyflux cahnhilliard vacancyflux cahnhilliard
initialVacancyConc 1e-6 initialVacancyConc 1e-6
(output) vacancyconc (output) vacancyconc

1
config/Homogenization_multiField.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
[SX] [SX]
type isostrain type isostrain
Ngrains 1 Ngrains 1

1
config/Kinematics_Hydrogen_Strain.config
View File

@ -1,3 +1,2 @@
### $Id$ ###
(kinematics) vacancy_strain (kinematics) vacancy_strain
vacancy_strain_coeff 0.006 vacancy_strain_coeff 0.006

1
config/Kinematics_Thermal_Expansion.config
View File

@ -1,3 +1,2 @@
### $Id$ ###
(kinematics) thermal_expansion (kinematics) thermal_expansion
thermal_expansion11 0.00231 thermal_expansion11 0.00231

1
config/Kinematics_Vacancy_Strain.config
View File

@ -1,3 +1,2 @@
### $Id$ ###
(kinematics) hydrogen_strain (kinematics) hydrogen_strain
hydrogen_strain_coeff 0.06 hydrogen_strain_coeff 0.06

1
config/Microstructure_DP_Steel.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
[DP_Steel] [DP_Steel]
/elementhomogeneous/ /elementhomogeneous/
crystallite 1 crystallite 1

1
config/Microstructure_ElementHomogeneous.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
[ElementHomogeneous] [ElementHomogeneous]
/elementhomogeneous/ # put this flag to set ips identical in one element (something like reduced integration) /elementhomogeneous/ # put this flag to set ips identical in one element (something like reduced integration)
crystallite 1 crystallite 1

1
config/Phase_Damage.config
View File

@ -1,3 +1,2 @@
### $Id$ ###
damage_diffusion11 1.0 damage_diffusion11 1.0
damage_mobility 0.001 damage_mobility 0.001

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

1
config/Phase_Dislotwin_TWIP-Steel-FeMnC.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
[TWIP_Steel_FeMnC] [TWIP_Steel_FeMnC]
elasticity hooke elasticity hooke

1
config/Phase_Dislotwin_Tungsten.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
[Tungsten] [Tungsten]
elasticity hooke elasticity hooke

1
config/Phase_Hydrogen.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
hydrogenflux_diffusion11 1.0 hydrogenflux_diffusion11 1.0
hydrogenflux_mobility11 1.0 hydrogenflux_mobility11 1.0
hydrogenVolume 1e-28 hydrogenVolume 1e-28

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

@ -1,9 +1,8 @@
### $Id$ ###
[Aluminum_Isotropic] [Aluminum_Isotropic]
# 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

1
config/Phase_None_IsotropicVolumePreservation.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
[IsotropicVolumePreservation] [IsotropicVolumePreservation]
elasticity hooke elasticity hooke

1
config/Phase_None_Orthorombic.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
[Orthorombic] [Orthorombic]
elasticity hooke elasticity hooke

1
config/Phase_Nonlocal_Aluminum.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
[Aluminum] [Aluminum]
elasticity hooke elasticity hooke

1
config/Phase_Phenopowerlaw_Aluminum.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
[Aluminum] [Aluminum]
elasticity hooke elasticity hooke
plasticity phenopowerlaw plasticity phenopowerlaw

1
config/Phase_Phenopowerlaw_BCC-Ferrite.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
# Tasan et.al. 2015 Acta Materalia # Tasan et.al. 2015 Acta Materalia
# Tasan et.al. 2015 International Journal of Plasticity # Tasan et.al. 2015 International Journal of Plasticity
# Diehl et.al. 2015 Meccanica # Diehl et.al. 2015 Meccanica

1
config/Phase_Phenopowerlaw_BCC-Martensite.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
# Tasan et.al. 2015 Acta Materalia # Tasan et.al. 2015 Acta Materalia
# Tasan et.al. 2015 International Journal of Plasticity # Tasan et.al. 2015 International Journal of Plasticity
# Diehl et.al. 2015 Meccanica # Diehl et.al. 2015 Meccanica

1
config/Phase_Phenopowerlaw_Gold.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
# parameters fitted by D. Ma to: # parameters fitted by D. Ma to:

1
config/Phase_Phenopowerlaw_cpTi-alpha.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
[cpTi-alpha] [cpTi-alpha]
plasticity phenopowerlaw plasticity phenopowerlaw
elasticity hooke elasticity hooke

1
config/Phase_Phenopowerlaw_multiField.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
[Aluminum] [Aluminum]
elasticity hooke elasticity hooke
plasticity phenopowerlaw plasticity phenopowerlaw

1
config/Phase_Porosity.config
View File

@ -1,3 +1,2 @@
### $Id$ ###
porosity_diffusion11 1.0 porosity_diffusion11 1.0
porosity_mobility 0.001 porosity_mobility 0.001

1
config/Phase_Thermal.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
thermal_conductivity11 237.0 thermal_conductivity11 237.0
specific_heat 910.0 specific_heat 910.0
mass_density 2700.0 mass_density 2700.0

1
config/Phase_Vacancy.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
vacancyflux_diffusion11 1.0 vacancyflux_diffusion11 1.0
vacancyflux_mobility11 1.0 vacancyflux_mobility11 1.0
vacancyFormationEnergy 1e-19 vacancyFormationEnergy 1e-19

1
config/Source_Damage_IsoBrittle.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
(source) damage_isoBrittle (source) damage_isoBrittle
isobrittle_criticalStrainEnergy 1400000.0 isobrittle_criticalStrainEnergy 1400000.0
isobrittle_atol 0.01 isobrittle_atol 0.01

1
config/Source_Thermal_Dissipation.config
View File

@ -1,3 +1,2 @@
### $Id$ ###
(source) thermal_dissipation (source) thermal_dissipation
dissipation_ColdWorkCoeff 0.95 dissipation_ColdWorkCoeff 0.95

1
config/Source_Vacancy_Irradiation.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
(source) vacancy_irradiation (source) vacancy_irradiation
irradiation_cascadeprobability 0.00001 irradiation_cascadeprobability 0.00001
irradiation_cascadevolume 1000.0 irradiation_cascadevolume 1000.0

1
config/Source_Vacancy_PhenoPlasticity.config
View File

@ -1,3 +1,2 @@
### $Id$ ###
(source) vacancy_phenoplasticity (source) vacancy_phenoplasticity
phenoplasticity_ratecoeff 0.01 phenoplasticity_ratecoeff 0.01

1
config/Texture_FiberExample.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
[FiberExample] [FiberExample]
axes x y -z # model coordinate x-, y-, z-axes correspond to which axes during texture measurement? this was a left handed coordinate system! axes x y -z # model coordinate x-, y-, z-axes correspond to which axes during texture measurement? this was a left handed coordinate system!
# fiber axis in spherical coordinates: alpha crystal system, beta sample system # fiber axis in spherical coordinates: alpha crystal system, beta sample system

1
config/Texture_Gauss_001.config
View File

@ -1,3 +1,2 @@
### $Id$ ###
[001] [001]
(gauss) phi1 0.000 Phi 0.000 phi2 0.000 scatter 0.000 fraction 1.000 (gauss) phi1 0.000 Phi 0.000 phi2 0.000 scatter 0.000 fraction 1.000

1
config/Texture_Gauss_101.config
View File

@ -1,3 +1,2 @@
### $Id$ ###
[101] [101]
(gauss) phi1 0.000 Phi 45.000 phi2 90.000 scatter 0.000 fraction 1.000 (gauss) phi1 0.000 Phi 45.000 phi2 90.000 scatter 0.000 fraction 1.000

1
config/Texture_Gauss_111.config
View File

@ -1,3 +1,2 @@
### $Id$ ###
[111] [111]
(gauss) phi1 0.000 Phi 54.7356 phi2 45.000 scatter 0.000 fraction 1.000 (gauss) phi1 0.000 Phi 54.7356 phi2 45.000 scatter 0.000 fraction 1.000

1
config/Texture_Gauss_123.config
View File

@ -1,3 +1,2 @@
### $Id$ ###
[123] [123]
(gauss) phi1 209.805 Phi 29.206 phi2 63.435 scatter 0.000 fraction 1.000 (gauss) phi1 209.805 Phi 29.206 phi2 63.435 scatter 0.000 fraction 1.000

1
config/Texture_RandomSingleCrystals.config
View File

@ -1,3 +1,2 @@
### $Id$ ###
[RandomSingleCrystals] [RandomSingleCrystals]
(random) scatter 0.000 fraction 1.000 (random) scatter 0.000 fraction 1.000

1
config/Texture_Rolling.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
[Rolling] [Rolling]
hybridIA rollingTexture.linearODF hybridIA rollingTexture.linearODF
symmetry orthotropic # or monoclinic symmetry orthotropic # or monoclinic

1
config/debug.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
### debugging parameters ### ### debugging parameters ###
# example: # example:

1
config/material.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
### numerical parameters ### ### numerical parameters ###
# The material.config file needs to specify five parts: # The material.config file needs to specify five parts:

1
config/numerics.config
View File

@ -1,4 +1,3 @@
### $Id$ ###
### numerical parameters ### ### numerical parameters ###
relevantStrain 1.0e-7 # strain increment considered significant (used by crystallite to determine whether strain inc is considered significant) relevantStrain 1.0e-7 # strain increment considered significant (used by crystallite to determine whether strain inc is considered significant)

373258
config/rollingTexture.linearODF
View File

File diff suppressed because it is too large Load Diff

1
examples/.gitignore vendored Normal file
View File

@ -0,0 +1 @@
postProc

4
examples/AbaqusStandard/material.config
View File

@ -1,7 +1,3 @@
#####################
# $Id$
#####################
#-------------------# #-------------------#
<homogenization> <homogenization>
#-------------------# #-------------------#

4
examples/MSC.Marc/material.config
View File

@ -1,7 +1,3 @@
#####################
# $Id$
#####################
#-------------------# #-------------------#
<homogenization> <homogenization>
#-------------------# #-------------------#

33
examples/SpectralMethod/example_postprocess.py
View File

@ -1,33 +0,0 @@
#!/usr/bin/env python
# -*- coding: UTF-8 no BOM -*-
import os
import glob
from subprocess import call
geom_name = '20grains16x16x16_tensionX'
postResults = 'postResults --cr f,p --split --separation x,y,z '+geom_name+'.spectralOut'
sts = call(postResults, shell=True)
os.chdir('./postProc/')
ascii_files = glob.glob(geom_name+'_inc*.txt')
print ascii_files
showTable = "showTable -a "
addCauchy = 'addCauchy '
addMises = 'addMises -s Cauchy '
addStrainTensors = "addStrainTensors -0 -v "
visualize3D = "3Dvisualize -s 'Mises(Cauchy)',1_p Cauchy "
postProc = [addCauchy, addMises, addStrainTensors, visualize3D]
for f in ascii_files:
print f
for p in postProc:
p = p+f
print p
sts = call(p,shell=True)

17
examples/SpectralMethod/example_run.py
View File

@ -1,17 +0,0 @@
#!/usr/bin/env python
# -*- coding: UTF-8 no BOM -*-
import sys
resolutions = [16,32,64]
resolution = resolutions[0]
try:
resolution = int(sys.argv[1])
except:
pass
if resolution not in resolutions:
resolution = resolutions[0]
from subprocess import call
call('make run%s'%('x'.join([str(resolution)]*3)), shell=True)

12
installation/compile_CoreModule.py
View File

@ -1,7 +1,7 @@
#!/usr/bin/env python #!/usr/bin/env python2
# -*- coding: UTF-8 no BOM -*- # -*- coding: UTF-8 no BOM -*-
import os,sys,glob,string,subprocess,shlex import os,sys,glob,subprocess,shlex
from damask import Environment from damask import Environment
from damask import version as DAMASKVERSION from damask import version as DAMASKVERSION
@ -47,9 +47,13 @@ compileOptions = ' -DSpectral -DFLOAT=8 -DINT=4 -I%s/lib -DDAMASKVERSION=\\\\\"\
#--- this saves the path of libraries to core.so, hence it is known during runtime ---------------- #--- this saves the path of libraries to core.so, hence it is known during runtime ----------------
if options['F90'] == 'gfortran': if options['F90'] == 'gfortran':
LDFLAGS = '-shared -Wl,-undefined,dynamic_lookup' # solved error: Undefined symbols for architecture x86_64: "_PyArg_ParseTupleAndKeywords" as found on https://lists.macosforge.org/pipermail/macports-dev/2013-May/022735.html # solved error: Undefined symbols for architecture x86_64: "_PyArg_ParseTupleAndKeywords"
# as found on https://lists.macosforge.org/pipermail/macports-dev/2013-May/022735.html
LDFLAGS = '-shared -Wl,-undefined,dynamic_lookup'
else: else:
LDFLAGS = ' -openmp -Wl' # some f2py versions/configurations compile with openMP, so linking against openMP is needed to prevent errors during loading of core module # some f2py versions/configurations compile with openMP, so linking against openMP is needed
# to prevent errors during loading of core module
LDFLAGS = ' -openmp -Wl'
#--- run path of for fftw during runtime ---------------------------------------------------------- #--- run path of for fftw during runtime ----------------------------------------------------------
LDFLAGS += ',-rpath,%s/lib,-rpath,%s/lib64'%(options['FFTW_ROOT'],options['FFTW_ROOT']) LDFLAGS += ',-rpath,%s/lib,-rpath,%s/lib64'%(options['FFTW_ROOT'],options['FFTW_ROOT'])

1
installation/mods_Abaqus/abaqus_v6.env
View File

@ -1,7 +1,6 @@
# #
# DAMASK Abaqus Environment File # DAMASK Abaqus Environment File
# #
# $Id$
# ------------------------------------ # ------------------------------------
# originally taken from Abaqus ver. 6.11.1 # originally taken from Abaqus ver. 6.11.1
# #

1
installation/mods_Abaqus/abaqus_v6_serial.env
View File

@ -1,7 +1,6 @@
# #
# DAMASK Abaqus Environment File # DAMASK Abaqus Environment File
# #
# $Id$
# ------------------------------------ # ------------------------------------
# originally taken from Abaqus ver. 6.11.1 # originally taken from Abaqus ver. 6.11.1
# #

354
installation/mods_Abaqus/abaqus_v6_windows.env
View File

@ -1,177 +1,177 @@
# #
# System-Wide Abaqus Environment File # System-Wide Abaqus Environment File
# ------------------------------------- # -------------------------------------
standard_parallel = ALL standard_parallel = ALL
mp_mode = MPI mp_mode = MPI
mp_file_system = (DETECT,DETECT) mp_file_system = (DETECT,DETECT)
mp_num_parallel_ftps = (4, 4) mp_num_parallel_ftps = (4, 4)
mp_environment_export = ('MPI_PROPAGATE_TSTP', mp_environment_export = ('MPI_PROPAGATE_TSTP',
'ABA_CM_BUFFERING', 'ABA_CM_BUFFERING',
'ABA_CM_BUFFERING_LIMIT', 'ABA_CM_BUFFERING_LIMIT',
'ABA_ITERATIVE_SOLVER_VERBOSE', 'ABA_ITERATIVE_SOLVER_VERBOSE',
'ABA_DMPSOLVER_BWDPARALLELOFF', 'ABA_DMPSOLVER_BWDPARALLELOFF',
'ABA_ELP_SURFACE_SPLIT', 'ABA_ELP_SURFACE_SPLIT',
'ABA_ELP_SUSPEND', 'ABA_ELP_SUSPEND',
'ABA_HOME', 'ABA_HOME',
'ABA_MEMORY_MODE', 'ABA_MEMORY_MODE',
'ABA_MPI_MESSAGE_TRACKING', 'ABA_MPI_MESSAGE_TRACKING',
'ABA_MPI_VERBOSE_LEVEL', 'ABA_MPI_VERBOSE_LEVEL',
'ABA_PATH', 'ABA_PATH',
'ABAQUS_CSE_RELTIMETOLERANCE', 'ABAQUS_CSE_RELTIMETOLERANCE',
'ABA_RESOURCE_MONITOR', 'ABA_RESOURCE_MONITOR',
'ABA_RESOURCE_USEMALLINFO', 'ABA_RESOURCE_USEMALLINFO',
'ABAQUS_LANG', 'ABAQUS_LANG',
'ABAQUS_CSE_CURRCONFIGMAPPING', 'ABAQUS_CSE_CURRCONFIGMAPPING',
'ABAQUS_MPF_DIAGNOSTIC_LEVEL', 'ABAQUS_MPF_DIAGNOSTIC_LEVEL',
'ABAQUSLM_LICENSE_FILE', 'ABAQUSLM_LICENSE_FILE',
'ABQ_CRTMALLOC', 'ABQ_CRTMALLOC',
'ABQ_DATACHECK', 'ABQ_DATACHECK',
'ABQ_RECOVER', 'ABQ_RECOVER',
'ABQ_RESTART', 'ABQ_RESTART',
'ABQ_SPLITFILE', 'ABQ_SPLITFILE',
'ABQ_XPL_WINDOWDUMP', 'ABQ_XPL_WINDOWDUMP',
'ABQ_XPL_PARTITIONSIZE', 'ABQ_XPL_PARTITIONSIZE',
'ABQLMHANGLIMIT', 'ABQLMHANGLIMIT',
'ABQLMQUEUE', 'ABQLMQUEUE',
'ABQLMUSER', 'ABQLMUSER',
'CCI_RENDEZVOUS', 'CCI_RENDEZVOUS',
'DOMAIN', 'DOMAIN',
'DOMAIN_CPUS', 'DOMAIN_CPUS',
'DOUBLE_PRECISION', 'DOUBLE_PRECISION',
'FLEXLM_DIAGNOSTICS', 'FLEXLM_DIAGNOSTICS',
'FOR0006', 'FOR0006',
'FOR0064', 'FOR0064',
'FOR_IGNORE_EXCEPTIONS', 'FOR_IGNORE_EXCEPTIONS',
'FOR_DISABLE_DIAGNOSTIC_DISPLAY', 'FOR_DISABLE_DIAGNOSTIC_DISPLAY',
'LD_PRELOAD', 'LD_PRELOAD',
'MP_NUMBER_OF_THREADS', 'MP_NUMBER_OF_THREADS',
'MPC_GANG', 'MPC_GANG',
'MPI_FLAGS', 'MPI_FLAGS',
'MPI_FLUSH_FCACHE', 'MPI_FLUSH_FCACHE',
'MPI_RDMA_NENVELOPE', 'MPI_RDMA_NENVELOPE',
'MPI_SOCKBUFSIZE', 'MPI_SOCKBUFSIZE',
'MPI_USE_MALLOPT_MMAP_MAX', 'MPI_USE_MALLOPT_MMAP_MAX',
'MPI_USE_MALLOPT_MMAP_THRESHOLD', 'MPI_USE_MALLOPT_MMAP_THRESHOLD',
'MPI_USE_MALLOPT_SBRK_PROTECTION', 'MPI_USE_MALLOPT_SBRK_PROTECTION',
'MPI_WORKDIR', 'MPI_WORKDIR',
'MPCCI_DEBUG', 'MPCCI_DEBUG',
'MPCCI_CODEID', 'MPCCI_CODEID',
'MPCCI_JOBID', 'MPCCI_JOBID',
'MPCCI_NETDEVICE', 'MPCCI_NETDEVICE',
'MPCCI_TINFO', 'MPCCI_TINFO',
'MPCCI_SERVER', 'MPCCI_SERVER',
'ABAQUS_CCI_DEBUG', 'ABAQUS_CCI_DEBUG',
'NCPUS', 'NCPUS',
'OMP_DYNAMIC', 'OMP_DYNAMIC',
'OMP_NUM_THREADS', 'OMP_NUM_THREADS',
'OUTDIR', 'OUTDIR',
'PAIDUP', 'PAIDUP',
'PARALLEL_METHOD', 'PARALLEL_METHOD',
'RAIDEV_NDREG_LAZYMEM', 'RAIDEV_NDREG_LAZYMEM',
'ABA_SYMBOLIC_GENERALCOLLAPSE', 'ABA_SYMBOLIC_GENERALCOLLAPSE',
'ABA_SYMBOLIC_GENERAL_MAXCLIQUERANK', 'ABA_SYMBOLIC_GENERAL_MAXCLIQUERANK',
'ABA_ADM_MINIMUMINCREASE', 'ABA_ADM_MINIMUMINCREASE',
'ABA_ADM_MINIMUMDECREASE', 'ABA_ADM_MINIMUMDECREASE',
'IPATH_NO_CPUAFFINITY', 'IPATH_NO_CPUAFFINITY',
'MALLOC_MMAP_THRESHOLD_', 'MALLOC_MMAP_THRESHOLD_',
'ABA_EXT_SIMOUTPUT', 'ABA_EXT_SIMOUTPUT',
'SMA_WS', 'SMA_WS',
'SMA_PARENT', 'SMA_PARENT',
'SMA_PLATFORM', 'SMA_PLATFORM',
'ABA_PRE_DECOMPOSITION', 'ABA_PRE_DECOMPOSITION',
'ACML_FAST_MALLOC', 'ACML_FAST_MALLOC',
'ACML_FAST_MALLOC_CHUNK_SIZE', 'ACML_FAST_MALLOC_CHUNK_SIZE',
'ACML_FAST_MALLOC_MAX_CHUNKS', 'ACML_FAST_MALLOC_MAX_CHUNKS',
'ACML_FAST_MALLOC_DEBUG') 'ACML_FAST_MALLOC_DEBUG')
import driverUtils, os import driverUtils, os
#-*- mode: python -*- #-*- mode: python -*-
# # # #
# Compile and Link command settings for the Windows 64 Platform # # Compile and Link command settings for the Windows 64 Platform #
# ( AMD Opteron / Intel EM64T ) # # ( AMD Opteron / Intel EM64T ) #
# # # #
compile_fortran=['ifort', compile_fortran=['ifort',
'/c','/DABQ_WIN86_64', '/u', '/c','/DABQ_WIN86_64', '/u',
'/iface:cref', '/recursive', '/Qauto-scalar', '/iface:cref', '/recursive', '/Qauto-scalar',
'/QxSSE3', '/QaxAVX', '/QxSSE3', '/QaxAVX',
'/heap-arrays:1', '/heap-arrays:1',
# '/Od', '/Ob0' # <-- Optimization # '/Od', '/Ob0' # <-- Optimization
# '/Zi', # <-- Debugging # '/Zi', # <-- Debugging
'/include:%I', '/free', '/O1', '/fpp', '/openmp', '/Qmkl'] '/include:%I', '/free', '/O1', '/fpp', '/openmp', '/Qmkl']
link_sl=['LINK', link_sl=['LINK',
'/nologo', '/NOENTRY', '/INCREMENTAL:NO', '/subsystem:console', '/machine:AMD64', '/nologo', '/NOENTRY', '/INCREMENTAL:NO', '/subsystem:console', '/machine:AMD64',
'/NODEFAULTLIB:LIBC.LIB', '/NODEFAULTLIB:LIBCMT.LIB', '/NODEFAULTLIB:LIBC.LIB', '/NODEFAULTLIB:LIBCMT.LIB',
'/DEFAULTLIB:OLDNAMES.LIB', '/DEFAULTLIB:LIBIFCOREMD.LIB', '/DEFAULTLIB:LIBIFPORTMD', '/DEFAULTLIB:LIBMMD.LIB', '/DEFAULTLIB:OLDNAMES.LIB', '/DEFAULTLIB:LIBIFCOREMD.LIB', '/DEFAULTLIB:LIBIFPORTMD', '/DEFAULTLIB:LIBMMD.LIB',
'/DEFAULTLIB:kernel32.lib', '/DEFAULTLIB:user32.lib', '/DEFAULTLIB:advapi32.lib', '/DEFAULTLIB:kernel32.lib', '/DEFAULTLIB:user32.lib', '/DEFAULTLIB:advapi32.lib',
'/FIXED:NO', '/dll', '/FIXED:NO', '/dll',
'/def:%E', '/out:%U', '%F', '%A', '%L', '%B', '/def:%E', '/out:%U', '%F', '%A', '%L', '%B',
'oldnames.lib', 'user32.lib', 'ws2_32.lib', 'netapi32.lib', 'advapi32.lib'] 'oldnames.lib', 'user32.lib', 'ws2_32.lib', 'netapi32.lib', 'advapi32.lib']
link_exe=['LINK', link_exe=['LINK',
'/nologo', '/INCREMENTAL:NO', '/subsystem:console', '/machine:AMD64', '/STACK:20000000', '/nologo', '/INCREMENTAL:NO', '/subsystem:console', '/machine:AMD64', '/STACK:20000000',
'/NODEFAULTLIB:LIBC.LIB', '/NODEFAULTLIB:LIBCMT.LIB', '/DEFAULTLIB:OLDNAMES.LIB', '/DEFAULTLIB:LIBIFCOREMD.LIB', '/NODEFAULTLIB:LIBC.LIB', '/NODEFAULTLIB:LIBCMT.LIB', '/DEFAULTLIB:OLDNAMES.LIB', '/DEFAULTLIB:LIBIFCOREMD.LIB',
'/DEFAULTLIB:LIBIFPORTMD', '/DEFAULTLIB:LIBMMD.LIB', '/DEFAULTLIB:kernel32.lib', '/DEFAULTLIB:LIBIFPORTMD', '/DEFAULTLIB:LIBMMD.LIB', '/DEFAULTLIB:kernel32.lib',
'/DEFAULTLIB:user32.lib', '/DEFAULTLIB:advapi32.lib', '/DEFAULTLIB:user32.lib', '/DEFAULTLIB:advapi32.lib',
'/FIXED:NO', '/LARGEADDRESSAWARE', '/FIXED:NO', '/LARGEADDRESSAWARE',
'/out:%J', '%F', '%M', '%L', '%B', '%O', '/out:%J', '%F', '%M', '%L', '%B', '%O',
'oldnames.lib', 'user32.lib', 'ws2_32.lib', 'netapi32.lib', 'advapi32.lib'] 'oldnames.lib', 'user32.lib', 'ws2_32.lib', 'netapi32.lib', 'advapi32.lib']
# Link command to be used for MAKE w/o fortran compiler. # Link command to be used for MAKE w/o fortran compiler.
# remove the pound signs in order to remove the comments and have the file take effect. # remove the pound signs in order to remove the comments and have the file take effect.
# #
#link_exe=['LINK', '/nologo', 'INCREMENTAL:NO', '/subsystem:console', '/machine:AMD64', '/NODEFAULTLIB:LIBC.LIB', '/NODEFAULTLIB:LIBCMT.LIB', #link_exe=['LINK', '/nologo', 'INCREMENTAL:NO', '/subsystem:console', '/machine:AMD64', '/NODEFAULTLIB:LIBC.LIB', '/NODEFAULTLIB:LIBCMT.LIB',
# '/DEFAULTLIB:OLDNAMES.LIB', '/DEFAULTLIB:MSVCRT.LIB', '/DEFAULTLIB:kernel32.lib', 'DEFAULTLIB:user32.lib', '/DEFAULTLIB:advapi32.lib', # '/DEFAULTLIB:OLDNAMES.LIB', '/DEFAULTLIB:MSVCRT.LIB', '/DEFAULTLIB:kernel32.lib', 'DEFAULTLIB:user32.lib', '/DEFAULTLIB:advapi32.lib',
# '/FIXED:NO', '/LARGEADDRESSAWARE', '/DEBUG', '/out:%J', '%F', '%M', '%L', '%B', '%O', 'oldnames.lib', 'user32.lib', 'ws2_32.lib', # '/FIXED:NO', '/LARGEADDRESSAWARE', '/DEBUG', '/out:%J', '%F', '%M', '%L', '%B', '%O', 'oldnames.lib', 'user32.lib', 'ws2_32.lib',
# 'netapi32.lib', 'advapi32.lib] # 'netapi32.lib', 'advapi32.lib]
# MPI Configuration # MPI Configuration
mp_mode = THREADS mp_mode = THREADS
mp_mpi_implementation = NATIVE mp_mpi_implementation = NATIVE
mp_rsh_command = 'dummy %H -l %U -n %C' mp_rsh_command = 'dummy %H -l %U -n %C'
mp_mpirun_path = {} mp_mpirun_path = {}
mpirun = '' mpirun = ''
progDir = os.environ.get('ProgramFiles','C:\\Program Files') progDir = os.environ.get('ProgramFiles','C:\\Program Files')
for mpiDir in ('Microsoft HPC Pack', 'Microsoft HPC Pack 2008 R2', 'Microsoft HPC Pack 2008', 'Microsoft HPC Pack 2008 SDK'): for mpiDir in ('Microsoft HPC Pack', 'Microsoft HPC Pack 2008 R2', 'Microsoft HPC Pack 2008', 'Microsoft HPC Pack 2008 SDK'):
mpirun = progDir + os.sep + mpiDir + os.sep + 'bin' + os.sep + 'mpiexec.exe' mpirun = progDir + os.sep + mpiDir + os.sep + 'bin' + os.sep + 'mpiexec.exe'
if os.path.exists(mpirun): if os.path.exists(mpirun):
mp_mpirun_path[NATIVE] = mpirun mp_mpirun_path[NATIVE] = mpirun
mp_mpirun_path[MSSDK] = os.path.join(progDir, mpiDir) mp_mpirun_path[MSSDK] = os.path.join(progDir, mpiDir)
break break
if os.environ.has_key('CCP_HOME'): if os.environ.has_key('CCP_HOME'):
from queueCCS import QueueCCS from queueCCS import QueueCCS
queues['default'] = QueueCCS(queueName='share') queues['default'] = QueueCCS(queueName='share')
queues['share'] = QueueCCS(queueName='share') queues['share'] = QueueCCS(queueName='share')
queues['local'] = QueueCCS(queueName='local') queues['local'] = QueueCCS(queueName='local')
queues['genxmlshare'] = QueueCCS(queueName='genxmlshare') queues['genxmlshare'] = QueueCCS(queueName='genxmlshare')
queues['genxmllocal'] = QueueCCS(queueName='genxmllocal') queues['genxmllocal'] = QueueCCS(queueName='genxmllocal')
del QueueCCS del QueueCCS
mpirun = os.path.join(os.environ['CCP_HOME'], 'bin', 'mpiexec.exe') mpirun = os.path.join(os.environ['CCP_HOME'], 'bin', 'mpiexec.exe')
if os.path.exists(mpirun): if os.path.exists(mpirun):
mp_mpirun_path[NATIVE] = mpirun mp_mpirun_path[NATIVE] = mpirun
run_mode=BATCH run_mode=BATCH
if mp_mpirun_path: if mp_mpirun_path:
mp_mode=MPI mp_mode=MPI
del progDir, mpiDir, mpirun del progDir, mpiDir, mpirun
graphicsEnv = driverUtils.locateFile(os.environ['ABA_PATH'],'site','graphicsConfig','env') graphicsEnv = driverUtils.locateFile(os.environ['ABA_PATH'],'site','graphicsConfig','env')
if graphicsEnv: if graphicsEnv:
execfile(graphicsEnv) execfile(graphicsEnv)
else: else:
raise 'Cannot find the graphics configuration environment file (graphicsConfig.env)' raise 'Cannot find the graphics configuration environment file (graphicsConfig.env)'
del driverUtils, os, graphicsEnv del driverUtils, os, graphicsEnv
license_server_type=FLEXNET license_server_type=FLEXNET
abaquslm_license_file="" abaquslm_license_file=""
doc_root=" doc_root="
doc_root_type="html" doc_root_type="html"
academic=RESEARCH academic=RESEARCH

2
installation/symlink_Code.py
View File

@ -1,4 +1,4 @@
#!/usr/bin/env python #!/usr/bin/env python2
# -*- coding: UTF-8 no BOM -*- # -*- coding: UTF-8 no BOM -*-
import os,sys import os,sys

6
installation/symlink_Processing.py
View File

@ -1,4 +1,4 @@
#!/usr/bin/env python #!/usr/bin/env python2
# -*- coding: UTF-8 no BOM -*- # -*- coding: UTF-8 no BOM -*-
# Makes postprocessing routines acessible from everywhere. # Makes postprocessing routines acessible from everywhere.
@ -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]))

2
lib/damask/.gitignore vendored Normal file
View File

@ -0,0 +1,2 @@
core.so
corientation.so

87
lib/damask/asciitable.py
View File

@ -1,6 +1,5 @@
# -*- coding: UTF-8 no BOM -*- # -*- coding: UTF-8 no BOM -*-
# $Id$
import os,sys import os,sys
import numpy as np import numpy as np
@ -25,7 +24,7 @@ class ASCIItable():
readonly = False, # no reading from file readonly = False, # no reading from file
): ):
self.__IO__ = {'output': [], self.__IO__ = {'output': [],
'buffered': buffered, 'buffered': buffered,
'labeled': labeled, # header contains labels 'labeled': labeled, # header contains labels
'labels': [], # labels according to file info 'labels': [], # labels according to file info
'readBuffer': [], # buffer to hold non-advancing reads 'readBuffer': [], # buffer to hold non-advancing reads
@ -35,18 +34,18 @@ class ASCIItable():
self.__IO__['inPlace'] = not outname and name and not readonly self.__IO__['inPlace'] = not outname and name and not readonly
if self.__IO__['inPlace']: outname = name + self.tmpext # transparently create tmp file if self.__IO__['inPlace']: outname = name + self.tmpext # transparently create tmp file
try: try:
self.__IO__['in'] = (open( name,'r') if os.access( name, os.R_OK) else None) if name else sys.stdin self.__IO__['in'] = (open( name,'r') if os.access( name, os.R_OK) else None) if name else sys.stdin
except TypeError: except TypeError:
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 +65,24 @@ 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 _quote(self,
what):
"""quote empty or white space-containing output"""
import re
return '{quote}{content}{quote}'.format(
quote = ('"' if str(what)=='' or re.search(r"\s",str(what)) else ''),
content = what)
# ------------------------------------------------------------------ # ------------------------------------------------------------------
def close(self, def close(self,
dismiss = False): dismiss = False):
@ -128,7 +145,7 @@ class ASCIItable():
the first row or, if keyword "head[*]" is present, the first row or, if keyword "head[*]" is present,
the last line of the header the last line of the header
""" """
import re import re,shlex
try: try:
self.__IO__['in'].seek(0) self.__IO__['in'].seek(0)
@ -143,7 +160,7 @@ class ASCIItable():
if self.__IO__['labeled']: # table features labels if self.__IO__['labeled']: # table features labels
self.info = [self.__IO__['in'].readline().strip() for i in xrange(1,int(m.group(1)))] self.info = [self.__IO__['in'].readline().strip() for i in xrange(1,int(m.group(1)))]
self.labels = self.__IO__['in'].readline().split() # store labels found in last line self.labels = shlex.split(self.__IO__['in'].readline()) # store labels found in last line
else: else:
@ -179,7 +196,7 @@ class ASCIItable():
"""write current header information (info + labels)""" """write current header information (info + labels)"""
head = ['{}\theader'.format(len(self.info)+self.__IO__['labeled'])] if header else [] head = ['{}\theader'.format(len(self.info)+self.__IO__['labeled'])] if header else []
head.append(self.info) head.append(self.info)
if self.__IO__['labeled']: head.append('\t'.join(self.labels)) if self.__IO__['labeled']: head.append('\t'.join(map(self._quote,self.labels)))
return self.output_write(head) return self.output_write(head)
@ -243,9 +260,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 +289,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 +300,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 +310,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 +329,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 +348,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 +362,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,9 +378,10 @@ 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))) \
else xrange(start,start+dim) if isinstance(labels, Iterable) and not isinstance(labels, str) \
else xrange(start,start+dim)
# ------------------------------------------------------------------ # ------------------------------------------------------------------
def info_append(self, def info_append(self,
what): what):
@ -372,9 +390,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):
@ -402,6 +420,8 @@ class ASCIItable():
advance = True, advance = True,
respectLabels = True): respectLabels = True):
"""read next line (possibly buffered) and parse it into data array""" """read next line (possibly buffered) and parse it into data array"""
import shlex
self.line = self.__IO__['readBuffer'].pop(0) if len(self.__IO__['readBuffer']) > 0 \ self.line = self.__IO__['readBuffer'].pop(0) if len(self.__IO__['readBuffer']) > 0 \
else self.__IO__['in'].readline().strip() # take buffered content or get next data row from file else self.__IO__['in'].readline().strip() # take buffered content or get next data row from file
@ -411,10 +431,10 @@ class ASCIItable():
self.line = self.line.rstrip('\n') self.line = self.line.rstrip('\n')
if self.__IO__['labeled'] and respectLabels: # if table has labels if self.__IO__['labeled'] and respectLabels: # if table has labels
items = self.line.split()[:len(self.__IO__['labels'])] # use up to label count (from original file info) items = shlex.split(self.line)[:len(self.__IO__['labels'])] # use up to label count (from original file info)
self.data = items if len(items) == len(self.__IO__['labels']) else [] # take entries if label count matches self.data = items if len(items) == len(self.__IO__['labels']) else [] # take entries if label count matches
else: else:
self.data = self.line.split() # otherwise take all self.data = shlex.split(self.line) # otherwise take all
return self.data != [] return self.data != []
@ -462,9 +482,9 @@ class ASCIItable():
if len(self.data) == 0: return True if len(self.data) == 0: return True
if isinstance(self.data[0],list): if isinstance(self.data[0],list):
return self.output_write([delimiter.join(map(str,items)) for items in self.data]) return self.output_write([delimiter.join(map(self._quote,items)) for items in self.data])
else: else:
return self.output_write(delimiter.join(map(str,self.data))) return self.output_write( delimiter.join(map(self._quote,self.data)))
# ------------------------------------------------------------------ # ------------------------------------------------------------------
def data_writeArray(self, def data_writeArray(self,
@ -518,7 +538,8 @@ class ASCIItable():
# ------------------------------------------------------------------ # ------------------------------------------------------------------
def microstructure_read(self, def microstructure_read(self,
grid, grid,
type = 'i'): type = 'i',
strict = False):
"""read microstructure data (from .geom format)""" """read microstructure data (from .geom format)"""
def datatype(item): def datatype(item):
return int(item) if type.lower() == 'i' else float(item) return int(item) if type.lower() == 'i' else float(item)
@ -537,6 +558,6 @@ class ASCIItable():
s = min(len(items), N-i) # prevent overflow of microstructure array s = min(len(items), N-i) # prevent overflow of microstructure array
microstructure[i:i+s] = items[:s] microstructure[i:i+s] = items[:s]
i += s i += len(items)
return microstructure return (microstructure, i == N and not self.data_read()) if strict else microstructure # check for proper point count and end of file

40
lib/damask/colormaps.py
View File

@ -303,36 +303,45 @@ class Colormap():
'interpolate', 'interpolate',
] ]
__predefined__ = { __predefined__ = {
'gray': {'left': Color('HSL',[0,1,1]), 'gray': {'left': Color('HSL',[0,1,1]),
'right': Color('HSL',[0,0,0.15]), 'right': Color('HSL',[0,0,0.15]),
'interpolate': 'perceptualuniform'}, 'interpolate': 'perceptualuniform'},
'grey': {'left': Color('HSL',[0,1,1]), 'grey': {'left': Color('HSL',[0,1,1]),
'right': Color('HSL',[0,0,0.15]), 'right': Color('HSL',[0,0,0.15]),
'interpolate': 'perceptualuniform'}, 'interpolate': 'perceptualuniform'},
'red': {'left': Color('HSL',[0,1,0.14]), 'red': {'left': Color('HSL',[0,1,0.14]),
'right': Color('HSL',[0,0.35,0.91]), 'right': Color('HSL',[0,0.35,0.91]),
'interpolate': 'perceptualuniform'}, 'interpolate': 'perceptualuniform'},
'green': {'left': Color('HSL',[0.33333,1,0.14]), 'green': {'left': Color('HSL',[0.33333,1,0.14]),
'right': Color('HSL',[0.33333,0.35,0.91]), 'right': Color('HSL',[0.33333,0.35,0.91]),
'interpolate': 'perceptualuniform'}, 'interpolate': 'perceptualuniform'},
'blue': {'left': Color('HSL',[0.66,1,0.14]), 'blue': {'left': Color('HSL',[0.66,1,0.14]),
'right': Color('HSL',[0.66,0.35,0.91]), 'right': Color('HSL',[0.66,0.35,0.91]),
'interpolate': 'perceptualuniform'}, 'interpolate': 'perceptualuniform'},
'seaweed': {'left': Color('HSL',[0.78,1.0,0.1]), 'seaweed': {'left': Color('HSL',[0.78,1.0,0.1]),
'right': Color('HSL',[0.40000,0.1,0.9]), 'right': Color('HSL',[0.40000,0.1,0.9]),
'interpolate': 'perceptualuniform'}, 'interpolate': 'perceptualuniform'},
'bluebrown': {'left': Color('HSL',[0.65,0.53,0.49]), 'bluebrown': {'left': Color('HSL',[0.65,0.53,0.49]),
'right': Color('HSL',[0.11,0.75,0.38]), 'right': Color('HSL',[0.11,0.75,0.38]),
'interpolate': 'perceptualuniform'}, 'interpolate': 'perceptualuniform'},
'redgreen': {'left': Color('HSL',[0.97,0.96,0.36]), 'redgreen': {'left': Color('HSL',[0.97,0.96,0.36]),
'right': Color('HSL',[0.33333,1.0,0.14]), 'right': Color('HSL',[0.33333,1.0,0.14]),
'interpolate': 'perceptualuniform'}, 'interpolate': 'perceptualuniform'},
'bluered': {'left': Color('HSL',[0.65,0.53,0.49]), 'bluered': {'left': Color('HSL',[0.65,0.53,0.49]),
'right': Color('HSL',[0.97,0.96,0.36]), 'right': Color('HSL',[0.97,0.96,0.36]),
'interpolate': 'perceptualuniform'}, 'interpolate': 'perceptualuniform'},
'blueredrainbow':{'left': Color('HSL',[2.0/3.0,1,0.5]), 'blueredrainbow':{'left': Color('HSL',[2.0/3.0,1,0.5]),
'right': Color('HSL',[0,1,0.5]), 'right': Color('HSL',[0,1,0.5]),
'interpolate': 'linear' }, 'interpolate': 'linear' },
'orientation': {'left': Color('RGB',[0.933334,0.878432,0.878431]),
'right': Color('RGB',[0.250980,0.007843,0.000000]),
'interpolate': 'perceptualuniform'},
'strain': {'left': Color('RGB',[0.941177,0.941177,0.870588]),
'right': Color('RGB',[0.266667,0.266667,0.000000]),
'interpolate': 'perceptualuniform'},
'stress': {'left': Color('RGB',[0.878432,0.874511,0.949019]),
'right': Color('RGB',[0.000002,0.000000,0.286275]),
'interpolate': 'perceptualuniform'},
} }
@ -344,7 +353,7 @@ class Colormap():
predefined = None predefined = None
): ):
if str(predefined).lower() in self.__predefined__: if predefined is not None:
left = self.__predefined__[predefined.lower()]['left'] left = self.__predefined__[predefined.lower()]['left']
right= self.__predefined__[predefined.lower()]['right'] right= self.__predefined__[predefined.lower()]['right']
interpolate = self.__predefined__[predefined.lower()]['interpolate'] interpolate = self.__predefined__[predefined.lower()]['interpolate']
@ -442,11 +451,12 @@ class Colormap():
format = format.lower() # consistent comparison basis format = format.lower() # consistent comparison basis
frac = 0.5*(np.array(crop) + 1.0) # rescale crop range to fractions frac = 0.5*(np.array(crop) + 1.0) # rescale crop range to fractions
colors = [self.color(float(i)/(steps-1)*(frac[1]-frac[0])+frac[0]).expressAs(model).color for i in xrange(steps)] colors = [self.color(float(i)/(steps-1)*(frac[1]-frac[0])+frac[0]).expressAs(model).color for i in xrange(steps)]
if format == 'paraview': if format == 'paraview':
colormap = ['<ColorMap name="'+str(name)+'" space="Diverging">'] \ colormap = ['[\n {{\n "ColorSpace" : "RGB", "Name" : "{}",\n "RGBPoints" : ['.format(name)] \
+ ['<Point x="%i"'%i + ' o="1" r="%g" g="%g" b="%g"/>'%(color[0],color[1],color[2],) for i,color in colors] \ + [' {:4d},{:8.6f},{:8.6f},{:8.6f},'.format(i,color[0],color[1],color[2],)
+ ['</ColorMap>'] for i,color in enumerate(colors[:-1])]\
+ [' {:4d},{:8.6f},{:8.6f},{:8.6f} '.format(i+1,colors[-1][0],colors[-1][1],colors[-1][2],)]\
+ [' ]\n }\n]']
elif format == 'gmsh': elif format == 'gmsh':
colormap = ['View.ColorTable = {'] \ colormap = ['View.ColorTable = {'] \

1
lib/damask/config/material.py
View File

@ -1,6 +1,5 @@
# -*- coding: UTF-8 no BOM -*- # -*- coding: UTF-8 no BOM -*-
# $Id$
import re import re
class Section(): class Section():

1
lib/damask/environment.py
View File

@ -1,6 +1,5 @@
# -*- coding: UTF-8 no BOM -*- # -*- coding: UTF-8 no BOM -*-
# $Id$
import os,subprocess,shlex import os,subprocess,shlex

1
lib/damask/geometry/geometry.py
View File

@ -1,6 +1,5 @@
# -*- coding: UTF-8 no BOM -*- # -*- coding: UTF-8 no BOM -*-
# $Id$
import damask.geometry import damask.geometry

1
lib/damask/geometry/marc.py
View File

@ -1,6 +1,5 @@
# -*- coding: UTF-8 no BOM -*- # -*- coding: UTF-8 no BOM -*-
# $Id$
from .geometry import Geometry from .geometry import Geometry

1
lib/damask/geometry/spectral.py
View File

@ -1,6 +1,5 @@
# -*- coding: UTF-8 no BOM -*- # -*- coding: UTF-8 no BOM -*-
# $Id$
from .geometry import Geometry from .geometry import Geometry

1
lib/damask/result.py
View File

@ -1,6 +1,5 @@
# -*- coding: UTF-8 no BOM -*- # -*- coding: UTF-8 no BOM -*-
# $Id$
import numpy as np import numpy as np
import sys import sys

2
lib/damask/setup_corientation.py
View File

@ -1,4 +1,4 @@
#!/usr/bin/env python #!/usr/bin/env python2
# -*- coding: UTF-8 no BOM -*- # -*- coding: UTF-8 no BOM -*-
from distutils.core import setup from distutils.core import setup

1
lib/damask/solver/abaqus.py
View File

@ -1,6 +1,5 @@
# -*- coding: UTF-8 no BOM -*- # -*- coding: UTF-8 no BOM -*-
# $Id$
from .solver import Solver from .solver import Solver

1
lib/damask/solver/marc.py
View File

@ -1,6 +1,5 @@
# -*- coding: UTF-8 no BOM -*- # -*- coding: UTF-8 no BOM -*-
# $Id$
from .solver import Solver from .solver import Solver

1
lib/damask/solver/solver.py
View File

@ -1,6 +1,5 @@
# -*- coding: UTF-8 no BOM -*- # -*- coding: UTF-8 no BOM -*-
# $Id$
import damask.solver import damask.solver

1
lib/damask/solver/spectral.py
View File

@ -1,6 +1,5 @@
# -*- coding: UTF-8 no BOM -*- # -*- coding: UTF-8 no BOM -*-
# $Id$
from .solver import Solver from .solver import Solver

148
lib/damask/test/test.py
View File

@ -1,6 +1,5 @@
# -*- coding: UTF-8 no BOM -*- # -*- coding: UTF-8 no BOM -*-
# $Id$
import os,sys,shutil import os,sys,shutil
import logging,logging.config import logging,logging.config
@ -22,11 +21,11 @@ class Test():
logger = logging.getLogger() logger = logging.getLogger()
logger.setLevel(0) logger.setLevel(0)
fh = logging.FileHandler('test.log') # create file handler which logs even debug messages fh = logging.FileHandler('test.log') # create file handler which logs even debug messages
fh.setLevel(logging.DEBUG) fh.setLevel(logging.DEBUG)
full = logging.Formatter('%(asctime)s - %(levelname)s: \n%(message)s') full = logging.Formatter('%(asctime)s - %(levelname)s: \n%(message)s')
fh.setFormatter(full) fh.setFormatter(full)
ch = logging.StreamHandler(stream=sys.stdout) # create console handler with a higher log level ch = logging.StreamHandler(stream=sys.stdout) # create console handler with a higher log level
ch.setLevel(logging.INFO) ch.setLevel(logging.INFO)
# create formatter and add it to the handlers # create formatter and add it to the handlers
plain = logging.Formatter('%(message)s') plain = logging.Formatter('%(message)s')
@ -41,7 +40,7 @@ class Test():
+'----------------------------------------------------------------') +'----------------------------------------------------------------')
self.dirBase = os.path.dirname(os.path.realpath(sys.modules[self.__class__.__module__].__file__)) self.dirBase = os.path.dirname(os.path.realpath(sys.modules[self.__class__.__module__].__file__))
self.parser = OptionParser( self.parser = OptionParser(
description = test_description+' (using class: $Id$)', description = test_description+' (using class: {})'.format(damask.version),
usage='./test.py [options]') usage='./test.py [options]')
self.updateRequested = False self.updateRequested = False
self.parser.add_option("-d", "--debug", action="store_true",\ self.parser.add_option("-d", "--debug", action="store_true",\
@ -60,10 +59,10 @@ class Test():
try: try:
self.postprocess(variant) self.postprocess(variant)
if not self.compare(variant): if not self.compare(variant):
return variant+1 # return culprit return variant+1 # return culprit
except Exception as e : except Exception as e :
logging.critical('\nWARNING:\n %s\n'%e) logging.critical('\nWARNING:\n {}\n'.format(e))
return variant+1 # return culprit return variant+1 # return culprit
return 0 return 0
else: else:
if not self.testPossible(): return -1 if not self.testPossible(): return -1
@ -74,13 +73,13 @@ class Test():
self.prepare(variant) self.prepare(variant)
self.run(variant) self.run(variant)
self.postprocess(variant) self.postprocess(variant)
if self.updateRequested: # update requested if self.updateRequested: # update requested
self.update(variant) self.update(variant)
elif not (self.options.accept or self.compare(variant)): # no update, do comparison elif not (self.options.accept or self.compare(variant)): # no update, do comparison
return variant+1 # return culprit return variant+1 # return culprit
except Exception as e : except Exception as e :
logging.critical('\nWARNING:\n %s\n'%e) logging.critical('\nWARNING:\n {}\n'.format(e))
return variant+1 # return culprit return variant+1 # return culprit
return 0 return 0
def testPossible(self): def testPossible(self):
@ -94,13 +93,13 @@ class Test():
try: try:
shutil.rmtree(self.dirCurrent()) shutil.rmtree(self.dirCurrent())
except: except:
logging.warning('removal of directory "%s" not possible...'%(self.dirCurrent())) logging.warning('removal of directory "{}" not possible...'.format(self.dirCurrent()))
status = status and False status = status and False
try: try:
os.mkdir(self.dirCurrent()) os.mkdir(self.dirCurrent())
except: except:
logging.critical('creation of directory "%s" failed...'%(self.dirCurrent())) logging.critical('creation of directory "{}" failed...'.format(self.dirCurrent()))
status = status and False status = status and False
return status return status
@ -193,19 +192,19 @@ class Test():
try: try:
shutil.copy2(self.fileInReference(file),self.fileInCurrent(targetfiles[i])) shutil.copy2(self.fileInReference(file),self.fileInCurrent(targetfiles[i]))
except: except:
logging.critical('Reference2Current: Unable to copy file %s'%file) logging.critical('Reference2Current: Unable to copy file "{}"'.format(file))
def copy_Base2Current(self,sourceDir,sourcefiles=[],targetfiles=[]): def copy_Base2Current(self,sourceDir,sourcefiles=[],targetfiles=[]):
source=os.path.normpath(os.path.join(self.dirBase,'../../../'+sourceDir)) source=os.path.normpath(os.path.join(self.dirBase,'../../..',sourceDir))
if len(targetfiles) == 0: targetfiles = sourcefiles if len(targetfiles) == 0: targetfiles = sourcefiles
for i,file in enumerate(sourcefiles): for i,file in enumerate(sourcefiles):
try: try:
shutil.copy2(os.path.join(source,file),self.fileInCurrent(targetfiles[i])) shutil.copy2(os.path.join(source,file),self.fileInCurrent(targetfiles[i]))
except: except:
logging.error(os.path.join(source,file)) logging.error(os.path.join(source,file))
logging.critical('Base2Current: Unable to copy file %s'%file) logging.critical('Base2Current: Unable to copy file "{}"'.format(file))
def copy_Current2Reference(self,sourcefiles=[],targetfiles=[]): def copy_Current2Reference(self,sourcefiles=[],targetfiles=[]):
@ -215,7 +214,7 @@ class Test():
try: try:
shutil.copy2(self.fileInCurrent(file),self.fileInReference(targetfiles[i])) shutil.copy2(self.fileInCurrent(file),self.fileInReference(targetfiles[i]))
except: except:
logging.critical('Current2Reference: Unable to copy file %s'%file) logging.critical('Current2Reference: Unable to copy file "{}"'.format(file))
def copy_Proof2Current(self,sourcefiles=[],targetfiles=[]): def copy_Proof2Current(self,sourcefiles=[],targetfiles=[]):
@ -225,7 +224,7 @@ class Test():
try: try:
shutil.copy2(self.fileInProof(file),self.fileInCurrent(targetfiles[i])) shutil.copy2(self.fileInProof(file),self.fileInCurrent(targetfiles[i]))
except: except:
logging.critical('Proof2Current: Unable to copy file %s'%file) logging.critical('Proof2Current: Unable to copy file "{}"'.format(file))
def copy_Current2Current(self,sourcefiles=[],targetfiles=[]): def copy_Current2Current(self,sourcefiles=[],targetfiles=[]):
@ -234,7 +233,7 @@ class Test():
try: try:
shutil.copy2(self.fileInReference(file),self.fileInCurrent(targetfiles[i])) shutil.copy2(self.fileInReference(file),self.fileInCurrent(targetfiles[i]))
except: except:
logging.critical('Current2Current: Unable to copy file %s'%file) logging.critical('Current2Current: Unable to copy file "{}"'.format(file))
def execute_inCurrentDir(self,cmd,streamIn=None): def execute_inCurrentDir(self,cmd,streamIn=None):
@ -252,7 +251,7 @@ class Test():
def compare_Array(self,File1,File2): def compare_Array(self,File1,File2):
import numpy as np import numpy as np
logging.info('comparing\n '+File1+'\n '+File2) logging.info('\n '.join(['comparing',File1,File2]))
table1 = damask.ASCIItable(name=File1,readonly=True) table1 = damask.ASCIItable(name=File1,readonly=True)
table1.head_read() table1.head_read()
len1=len(table1.info)+2 len1=len(table1.info)+2
@ -270,8 +269,9 @@ class Test():
max_loc=np.argmax(abs(refArrayNonZero[curArray.nonzero()]/curArray[curArray.nonzero()]-1.)) max_loc=np.argmax(abs(refArrayNonZero[curArray.nonzero()]/curArray[curArray.nonzero()]-1.))
refArrayNonZero = refArrayNonZero[curArray.nonzero()] refArrayNonZero = refArrayNonZero[curArray.nonzero()]
curArray = curArray[curArray.nonzero()] curArray = curArray[curArray.nonzero()]
print(' ********\n * maximum relative error %e for %e and %e\n ********' print(' ********\n * maximum relative error {} between {} and {}\n ********'.format(max_err,
%(max_err, refArrayNonZero[max_loc],curArray[max_loc])) refArrayNonZero[max_loc],
curArray[max_loc]))
return max_err return max_err
else: else:
raise Exception('mismatch in array size to compare') raise Exception('mismatch in array size to compare')
@ -295,7 +295,7 @@ class Test():
absoluteTolerance=False,perLine=False,skipLines=[]): absoluteTolerance=False,perLine=False,skipLines=[]):
import numpy as np import numpy as np
logging.info('comparing ASCII Tables\n %s \n %s'%(file0,file1)) logging.info('\n '.join(['comparing ASCII Tables',file0,file1]))
if normHeadings == '': normHeadings = headings0 if normHeadings == '': normHeadings = headings0
# check if comparison is possible and determine lenght of columns # check if comparison is possible and determine lenght of columns
@ -315,7 +315,7 @@ class Test():
for i in xrange(dataLength): for i in xrange(dataLength):
if headings0[i]['shape'] != headings1[i]['shape']: if headings0[i]['shape'] != headings1[i]['shape']:
raise Exception('shape mismatch when comparing %s with %s '%(headings0[i]['label'],headings1[i]['label'])) raise Exception('shape mismatch between {} and {} '.format(headings0[i]['label'],headings1[i]['label']))
shape[i] = headings0[i]['shape'] shape[i] = headings0[i]['shape']
for j in xrange(np.shape(shape[i])[0]): for j in xrange(np.shape(shape[i])[0]):
length[i] *= shape[i][j] length[i] *= shape[i][j]
@ -323,7 +323,9 @@ class Test():
for j in xrange(np.shape(normShape[i])[0]): for j in xrange(np.shape(normShape[i])[0]):
normLength[i] *= normShape[i][j] normLength[i] *= normShape[i][j]
else: else:
raise Exception('trying to compare %i with %i normed by %i data sets'%(len(headings0),len(headings1),len(normHeadings))) raise Exception('trying to compare {} with {} normed by {} data sets'.format(len(headings0),
len(headings1),
len(normHeadings)))
table0 = damask.ASCIItable(name=file0,readonly=True) table0 = damask.ASCIItable(name=file0,readonly=True)
table0.head_read() table0.head_read()
@ -331,37 +333,34 @@ class Test():
table1.head_read() table1.head_read()
for i in xrange(dataLength): for i in xrange(dataLength):
key0 = {True :'1_%s', key0 = ('1_' if length[i]>1 else '') + headings0[i]['label']
False:'%s' }[length[i]>1]%headings0[i]['label'] key1 = ('1_' if length[i]>1 else '') + headings1[i]['label']
key1 = {True :'1_%s', normKey = ('1_' if normLength[i]>1 else '') + normHeadings[i]['label']
False:'%s' }[length[i]>1]%headings1[i]['label']
normKey = {True :'1_%s',
False:'%s' }[normLength[i]>1]%normHeadings[i]['label']
if key0 not in table0.labels: if key0 not in table0.labels:
raise Exception('column %s not found in 1. table...\n'%key0) raise Exception('column {} not found in 1. table...\n'.format(key0))
elif key1 not in table1.labels: elif key1 not in table1.labels:
raise Exception('column %s not found in 2. table...\n'%key1) raise Exception('column {} not found in 2. table...\n'.format(key1))
elif normKey not in table0.labels: elif normKey not in table0.labels:
raise Exception('column %s not found in 1. table...\n'%normKey) raise Exception('column {} not found in 1. table...\n'.format(normKey))
else: else:
column[0][i] = table0.labels.index(key0) # remember columns of requested data column[0][i] = table0.labels.index(key0)
column[1][i] = table1.labels.index(key1) # remember columns of requested data in second column column[1][i] = table1.labels.index(key1)
normColumn[i] = table0.labels.index(normKey) # remember columns of requested data in second column normColumn[i] = table0.labels.index(normKey)
line0 = 0 line0 = 0
while table0.data_read(): # read next data line of ASCII table while table0.data_read(): # read next data line of ASCII table
if line0 not in skipLines: if line0 not in skipLines:
for i in xrange(dataLength): for i in xrange(dataLength):
myData = np.array(map(float,table0.data[column[0][i]:\ myData = np.array(map(float,table0.data[column[0][i]:\
column[0][i]+length[i]]),'d') column[0][i]+length[i]]),'d')
normData = np.array(map(float,table0.data[normColumn[i]:\ normData = np.array(map(float,table0.data[normColumn[i]:\
normColumn[i]+normLength[i]]),'d') normColumn[i]+normLength[i]]),'d')
data[i] = np.append(data[i],np.reshape(myData,shape[i])) data[i] = np.append(data[i],np.reshape(myData,shape[i]))
if normType == 'pInf': if normType == 'pInf':
norm[i] = np.append(norm[i],np.max(np.abs(normData))) norm[i] = np.append(norm[i],np.max(np.abs(normData)))
else: else:
norm[i] = np.append(norm[i],np.linalg.norm(np.reshape(normData,normShape[i]),normType)) norm[i] = np.append(norm[i],np.linalg.norm(np.reshape(normData,normShape[i]),normType))
line0 +=1 line0 += 1
for i in xrange(dataLength): for i in xrange(dataLength):
if not perLine: norm[i] = [np.max(norm[i]) for j in xrange(line0-len(skipLines))] if not perLine: norm[i] = [np.max(norm[i]) for j in xrange(line0-len(skipLines))]
@ -370,12 +369,12 @@ class Test():
norm[i] = [1.0 for j in xrange(line0-len(skipLines))] norm[i] = [1.0 for j in xrange(line0-len(skipLines))]
absTol[i] = True absTol[i] = True
if perLine: if perLine:
logging.warning('At least one norm of %s in 1. table is 0.0, using absolute tolerance'%headings0[i]['label']) logging.warning('At least one norm of {} in 1. table is 0.0, using absolute tolerance'.format(headings0[i]['label']))
else: else:
logging.warning('Maximum norm of %s in 1. table is 0.0, using absolute tolerance'%headings0[i]['label']) logging.warning('Maximum norm of {} in 1. table is 0.0, using absolute tolerance'.format(headings0[i]['label']))
line1 = 0 line1 = 0
while table1.data_read(): # read next data line of ASCII table while table1.data_read(): # read next data line of ASCII table
if line1 not in skipLines: if line1 not in skipLines:
for i in xrange(dataLength): for i in xrange(dataLength):
myData = np.array(map(float,table1.data[column[1][i]:\ myData = np.array(map(float,table1.data[column[1][i]:\
@ -384,21 +383,25 @@ class Test():
norm[i][line1-len(skipLines)]) norm[i][line1-len(skipLines)])
line1 +=1 line1 +=1
if (line0 != line1): raise Exception('found %s lines in 1. table and %s in 2. table'%(line0,line1)) if (line0 != line1): raise Exception('found {} lines in 1. table but {} in 2. table'.format(line0,line1))
logging.info(' ********') logging.info(' ********')
for i in xrange(dataLength): for i in xrange(dataLength):
if absTol[i]: if absTol[i]:
logging.info(' * maximum absolute error %e for %s and %s'%(maxError[i],headings0[i]['label'],headings1[i]['label'])) logging.info(' * maximum absolute error {} between {} and {}'.format(maxError[i],
headings0[i]['label'],
headings1[i]['label']))
else: else:
logging.info(' * maximum relative error %e for %s and %s'%(maxError[i],headings0[i]['label'],headings1[i]['label'])) logging.info(' * maximum relative error {} between {} and {}'.format(maxError[i],
headings0[i]['label'],
headings1[i]['label']))
logging.info(' ********') logging.info(' ********')
return maxError return maxError
def compare_TablesStatistically(self, def compare_TablesStatistically(self,
files = [None,None], # list of file names files = [None,None], # list of file names
columns = [None], # list of list of column labels (per file) columns = [None], # list of list of column labels (per file)
meanTol = 1.0e-4, meanTol = 1.0e-4,
stdTol = 1.0e-6, stdTol = 1.0e-6,
preFilter = 1.0e-9): preFilter = 1.0e-9):
@ -407,18 +410,18 @@ class Test():
threshold can be used to ignore small values (a negative number disables this feature) threshold can be used to ignore small values (a negative number disables this feature)
""" """
if not (isinstance(files, Iterable) and not isinstance(files, str)): # check whether list of files is requested if not (isinstance(files, Iterable) and not isinstance(files, str)): # check whether list of files is requested
files = [str(files)] files = [str(files)]
tables = [damask.ASCIItable(name = filename,readonly = True) for filename in files] tables = [damask.ASCIItable(name = filename,readonly = True) for filename in files]
for table in tables: for table in tables:
table.head_read() table.head_read()
columns += [columns[0]]*(len(files)-len(columns)) # extend to same length as files columns += [columns[0]]*(len(files)-len(columns)) # extend to same length as files
columns = columns[:len(files)] # truncate to same length as files columns = columns[:len(files)] # truncate to same length as files
for i,column in enumerate(columns): for i,column in enumerate(columns):
if column is None: columns[i] = tables[i].labels # if no column is given, read all if column is None: columns[i] = tables[i].labels # if no column is given, read all
logging.info('comparing ASCIItables statistically') logging.info('comparing ASCIItables statistically')
for i in xrange(len(columns)): for i in xrange(len(columns)):
@ -428,7 +431,7 @@ class Test():
) )
logging.info(files[i]+':'+','.join(columns[i])) logging.info(files[i]+':'+','.join(columns[i]))
if len(files) < 2: return True # single table is always close to itself... if len(files) < 2: return True # single table is always close to itself...
data = [] data = []
for table,labels in zip(tables,columns): for table,labels in zip(tables,columns):
@ -443,16 +446,16 @@ class Test():
normedDelta = np.where(normBy>preFilter,delta/normBy,0.0) normedDelta = np.where(normBy>preFilter,delta/normBy,0.0)
mean = np.amax(np.abs(np.mean(normedDelta,0))) mean = np.amax(np.abs(np.mean(normedDelta,0)))
std = np.amax(np.std(normedDelta,0)) std = np.amax(np.std(normedDelta,0))
logging.info('mean: %f'%mean) logging.info('mean: {:f}'.format(mean))
logging.info('std: %f'%std) logging.info('std: {:f}'.format(std))
return (mean<meanTol) & (std < stdTol) return (mean<meanTol) & (std < stdTol)
def compare_Tables(self, def compare_Tables(self,
files = [None,None], # list of file names files = [None,None], # list of file names
columns = [None], # list of list of column labels (per file) columns = [None], # list of list of column labels (per file)
rtol = 1e-5, rtol = 1e-5,
atol = 1e-8, atol = 1e-8,
preFilter = -1.0, preFilter = -1.0,
@ -463,18 +466,18 @@ class Test():
threshold can be used to ignore small values (a negative number disables this feature) threshold can be used to ignore small values (a negative number disables this feature)
""" """
if not (isinstance(files, Iterable) and not isinstance(files, str)): # check whether list of files is requested if not (isinstance(files, Iterable) and not isinstance(files, str)): # check whether list of files is requested
files = [str(files)] files = [str(files)]
tables = [damask.ASCIItable(name = filename,readonly = True) for filename in files] tables = [damask.ASCIItable(name = filename,readonly = True) for filename in files]
for table in tables: for table in tables:
table.head_read() table.head_read()
columns += [columns[0]]*(len(files)-len(columns)) # extend to same length as files columns += [columns[0]]*(len(files)-len(columns)) # extend to same length as files
columns = columns[:len(files)] # truncate to same length as files columns = columns[:len(files)] # truncate to same length as files
for i,column in enumerate(columns): for i,column in enumerate(columns):
if column is None: columns[i] = tables[i].labels # if no column is given, read all if column is None: columns[i] = tables[i].labels # if no column is given, read all
logging.info('comparing ASCIItables') logging.info('comparing ASCIItables')
for i in xrange(len(columns)): for i in xrange(len(columns)):
@ -484,7 +487,7 @@ class Test():
) )
logging.info(files[i]+':'+','.join(columns[i])) logging.info(files[i]+':'+','.join(columns[i]))
if len(files) < 2: return True # single table is always close to itself... if len(files) < 2: return True # single table is always close to itself...
maximum = np.zeros(len(columns[0]),dtype='f') maximum = np.zeros(len(columns[0]),dtype='f')
data = [] data = []
@ -495,26 +498,26 @@ class Test():
table.close() table.close()
maximum /= len(tables) maximum /= len(tables)
maximum = np.where(maximum >0.0, maximum, 1) # do not devide by zero for empty columns maximum = np.where(maximum >0.0, maximum, 1) # avoid div by zero for empty columns
for i in xrange(len(data)): for i in xrange(len(data)):
data[i] /= maximum data[i] /= maximum
mask = np.zeros_like(table.data,dtype='bool') mask = np.zeros_like(table.data,dtype='bool')
for table in data: for table in data:
mask |= np.where(np.abs(table)<postFilter,True,False) # mask out (all) tiny values mask |= np.where(np.abs(table)<postFilter,True,False) # mask out (all) tiny values
allclose = True # start optimistic allclose = True # start optimistic
for i in xrange(1,len(data)): for i in xrange(1,len(data)):
if debug: if debug:
t0 = np.where(mask,0.0,data[i-1]) t0 = np.where(mask,0.0,data[i-1])
t1 = np.where(mask,0.0,data[i ]) t1 = np.where(mask,0.0,data[i ])
j = np.argmin(np.abs(t1)*rtol+atol-np.abs(t0-t1)) j = np.argmin(np.abs(t1)*rtol+atol-np.abs(t0-t1))
logging.info('%f'%np.amax(np.abs(t0-t1)/(np.abs(t1)*rtol+atol))) logging.info('{:f}'.format(np.amax(np.abs(t0-t1)/(np.abs(t1)*rtol+atol))))
logging.info('%f %f'%((t0*maximum).flatten()[j],(t1*maximum).flatten()[j])) logging.info('{:f} {:f}'.format((t0*maximum).flatten()[j],(t1*maximum).flatten()[j]))
allclose &= np.allclose(np.where(mask,0.0,data[i-1]), allclose &= np.allclose(np.where(mask,0.0,data[i-1]),
np.where(mask,0.0,data[i ]),rtol,atol) # accumulate "pessimism" np.where(mask,0.0,data[i ]),rtol,atol) # accumulate "pessimism"
return allclose return allclose
@ -543,14 +546,13 @@ class Test():
def report_Success(self,culprit): def report_Success(self,culprit):
if culprit == 0: if culprit == 0:
logging.critical('%s passed.'%({False: 'The test', logging.critical(('The test' if len(self.variants) == 1 else 'All {} tests'.format(len(self.variants))) + ' passed')
True: 'All %i tests'%(len(self.variants))}[len(self.variants) > 1]))
logging.critical('\n!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n') logging.critical('\n!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n')
return 0 return 0
if culprit == -1: if culprit == -1:
logging.warning('Warning: Could not start test') logging.warning('Warning: Could not start test')
return 0 return 0
else: else:
logging.critical(' ********\n * Test %i failed...\n ********'%(culprit)) logging.critical(' ********\n * Test {} failed...\n ********'.format(culprit))
logging.critical('\n!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n') logging.critical('\n!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n')
return culprit return culprit

Some files were not shown because too many files have changed in this diff Show More