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

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

This commit is contained in:
Martin Diehl 2016-09-27 18:54:06 +02:00
parent 4a627161d5 3abf0bc2d3
commit aeb633a42d
17 changed files with 494 additions and 295 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
VERSION
View File

@ -1 +1 @@
v2.0.1-130-gf3308db v2.0.1-150-g5345b42

23
code/C_routines.c
View File

@ -9,6 +9,14 @@
/* http://stackoverflow.com/questions/30279228/is-there-an-alternative-to-getcwd-in-fortran-2003-2008 */ /* http://stackoverflow.com/questions/30279228/is-there-an-alternative-to-getcwd-in-fortran-2003-2008 */
int isdirectory_c(const char *dir){
struct stat statbuf;
if(stat(dir, &statbuf) != 0)
return 0;
return S_ISDIR(statbuf.st_mode);
}
void getcurrentworkdir_c(char cwd[], int *stat ){ void getcurrentworkdir_c(char cwd[], int *stat ){
char cwd_tmp[1024]; char cwd_tmp[1024];
if(getcwd(cwd_tmp, sizeof(cwd_tmp)) == cwd_tmp){ if(getcwd(cwd_tmp, sizeof(cwd_tmp)) == cwd_tmp){
@ -20,9 +28,14 @@ void getcurrentworkdir_c(char cwd[], int *stat ){
} }
} }
int isdirectory_c(const char *dir){
struct stat statbuf; void gethostname_c(char hostname[], int *stat ){
if(stat(dir, &statbuf) != 0) char hostname_tmp[1024];
return 0; if(gethostname(hostname_tmp, sizeof(hostname_tmp)) == 0){
return S_ISDIR(statbuf.st_mode); strcpy(hostname,hostname_tmp);
*stat = 0;
}
else{
*stat = 1;
}
} }

2
code/IO.f90
View File

@ -1568,8 +1568,6 @@ subroutine IO_error(error_ID,el,ip,g,ext_msg)
msg = 'math_check: R*v == q*v failed' msg = 'math_check: R*v == q*v failed'
case (410_pInt) case (410_pInt)
msg = 'eigenvalues computation error' msg = 'eigenvalues computation error'
case (450_pInt)
msg = 'unknown symmetry type specified'
!------------------------------------------------------------------------------------------------- !-------------------------------------------------------------------------------------------------
! homogenization errors ! homogenization errors

36
code/math.f90
View File

@ -1450,7 +1450,7 @@ pure function math_EulerToQ(eulerangles)
cos(halfangles(1)-halfangles(3)) * s, & cos(halfangles(1)-halfangles(3)) * s, &
sin(halfangles(1)-halfangles(3)) * s, & sin(halfangles(1)-halfangles(3)) * s, &
sin(halfangles(1)+halfangles(3)) * c ] sin(halfangles(1)+halfangles(3)) * c ]
math_EulerToQ = math_qConj(math_EulerToQ) ! convert to passive rotation math_EulerToQ = math_qConj(math_EulerToQ) ! convert to passive rotation
end function math_EulerToQ end function math_EulerToQ
@ -1508,7 +1508,7 @@ pure function math_EulerAxisAngleToR(axis,omega)
real(pReal), dimension(3), intent(in) :: axis real(pReal), dimension(3), intent(in) :: axis
real(pReal), intent(in) :: omega real(pReal), intent(in) :: omega
math_EulerAxisAngleToR = transpose(math_axisAngleToR(axis,omega)) ! convert to passive rotation math_EulerAxisAngleToR = transpose(math_axisAngleToR(axis,omega)) ! convert to passive rotation
end function math_EulerAxisAngleToR end function math_EulerAxisAngleToR
@ -1527,7 +1527,7 @@ pure function math_EulerAxisAngleToQ(axis,omega)
real(pReal), dimension(3), intent(in) :: axis real(pReal), dimension(3), intent(in) :: axis
real(pReal), intent(in) :: omega real(pReal), intent(in) :: omega
math_EulerAxisAngleToQ = math_qConj(math_axisAngleToQ(axis,omega)) ! convert to passive rotation math_EulerAxisAngleToQ = math_qConj(math_axisAngleToQ(axis,omega)) ! convert to passive rotation
end function math_EulerAxisAngleToQ end function math_EulerAxisAngleToQ
@ -1550,7 +1550,7 @@ pure function math_axisAngleToQ(axis,omega)
norm = sqrt(math_mul3x3(axis,axis)) norm = sqrt(math_mul3x3(axis,axis))
rotation: if (norm > 1.0e-8_pReal) then rotation: if (norm > 1.0e-8_pReal) then
axisNrm = axis/norm ! normalize axis to be sure axisNrm = axis/norm ! normalize axis to be sure
math_axisAngleToQ = [cos(0.5_pReal*omega), sin(0.5_pReal*omega) * axisNrm(1:3)] math_axisAngleToQ = [cos(0.5_pReal*omega), sin(0.5_pReal*omega) * axisNrm(1:3)]
else rotation else rotation
math_axisAngleToQ = [1.0_pReal,0.0_pReal,0.0_pReal,0.0_pReal] math_axisAngleToQ = [1.0_pReal,0.0_pReal,0.0_pReal,0.0_pReal]
@ -1864,37 +1864,29 @@ end function math_sampleGaussVar
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief symmetrically equivalent Euler angles for given sample symmetry 1:triclinic, 2:monoclinic, 4:orthotropic !> @brief symmetrically equivalent Euler angles for given sample symmetry
!> @detail 1 (equivalent to != 2 and !=4):triclinic, 2:monoclinic, 4:orthotropic
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function math_symmetricEulers(sym,Euler) pure function math_symmetricEulers(sym,Euler)
implicit none implicit none
integer(pInt), intent(in) :: sym integer(pInt), intent(in) :: sym !< symmetry Class
real(pReal), dimension(3), intent(in) :: Euler real(pReal), dimension(3), intent(in) :: Euler
real(pReal), dimension(3,3) :: math_symmetricEulers real(pReal), dimension(3,3) :: math_symmetricEulers
integer(pInt) :: i,j
math_symmetricEulers(1,1) = PI+Euler(1) math_symmetricEulers = transpose(reshape([PI+Euler(1), PI-Euler(1), 2.0_pReal*PI-Euler(1), &
math_symmetricEulers(2,1) = Euler(2) Euler(2), PI-Euler(2), PI -Euler(2), &
math_symmetricEulers(3,1) = Euler(3) Euler(3), PI+Euler(3), PI +Euler(3)],[3,3])) ! transpose is needed to have symbolic notation instead of column-major
math_symmetricEulers(1,2) = PI-Euler(1) math_symmetricEulers = modulo(math_symmetricEulers,2.0_pReal*pi)
math_symmetricEulers(2,2) = PI-Euler(2)
math_symmetricEulers(3,2) = PI+Euler(3)
math_symmetricEulers(1,3) = 2.0_pReal*PI-Euler(1)
math_symmetricEulers(2,3) = PI-Euler(2)
math_symmetricEulers(3,3) = PI+Euler(3)
forall (i=1_pInt:3_pInt,j=1_pInt:3_pInt) math_symmetricEulers(j,i) = modulo(math_symmetricEulers(j,i),2.0_pReal*pi)
select case (sym) select case (sym)
case (4_pInt) ! all done case (4_pInt) ! orthotropic: all done
case (2_pInt) ! return only first case (2_pInt) ! monoclinic: return only first
math_symmetricEulers(1:3,2:3) = 0.0_pReal math_symmetricEulers(1:3,2:3) = 0.0_pReal
case default ! return blank case default ! triclinic: return blank
math_symmetricEulers = 0.0_pReal math_symmetricEulers = 0.0_pReal
end select end select

176
code/spectral_interface.f90
View File

@ -20,6 +20,7 @@ module DAMASK_interface
geometryFile = '', & !< parameter given for geometry file geometryFile = '', & !< parameter given for geometry file
loadCaseFile = '' !< parameter given for load case file loadCaseFile = '' !< parameter given for load case file
character(len=1024), private :: workingDirectory !< accessed by getSolverWorkingDirectoryName for compatibility reasons character(len=1024), private :: workingDirectory !< accessed by getSolverWorkingDirectoryName for compatibility reasons
character, private,parameter :: pathSep = '/'
public :: & public :: &
getSolverWorkingDirectoryName, & getSolverWorkingDirectoryName, &
@ -31,7 +32,6 @@ module DAMASK_interface
getLoadCaseFile, & getLoadCaseFile, &
rectifyPath, & rectifyPath, &
makeRelativePath, & makeRelativePath, &
getPathSep, &
IIO_stringValue, & IIO_stringValue, &
IIO_intValue, & IIO_intValue, &
IIO_lc, & IIO_lc, &
@ -44,6 +44,8 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine DAMASK_interface_init() subroutine DAMASK_interface_init()
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 system_routines, only: &
getHostName
implicit none implicit none
character(len=1024) :: & character(len=1024) :: &
@ -64,6 +66,7 @@ subroutine DAMASK_interface_init()
integer, dimension(8) :: & integer, dimension(8) :: &
dateAndTime ! type default integer dateAndTime ! type default integer
PetscErrorCode :: ierr PetscErrorCode :: ierr
logical :: error
external :: & external :: &
quit,& quit,&
MPI_Comm_rank,& MPI_Comm_rank,&
@ -116,54 +119,52 @@ subroutine DAMASK_interface_init()
tag = IIO_lc(IIO_stringValue(commandLine,chunkPos,i)) ! extract key tag = IIO_lc(IIO_stringValue(commandLine,chunkPos,i)) ! extract key
select case(tag) select case(tag)
case ('-h','--help') case ('-h','--help')
mainProcess2: if (worldrank == 0) then write(6,'(a)') ' #######################################################################'
write(6,'(a)') ' #######################################################################' write(6,'(a)') ' DAMASK_spectral:'
write(6,'(a)') ' DAMASK_spectral:' write(6,'(a)') ' The spectral method boundary value problem solver for'
write(6,'(a)') ' The spectral method boundary value problem solver for' write(6,'(a)') ' the Düsseldorf Advanced Material Simulation Kit'
write(6,'(a)') ' the Düsseldorf Advanced Material Simulation Kit' write(6,'(a,/)')' #######################################################################'
write(6,'(a,/)')' #######################################################################' write(6,'(a,/)')' Valid command line switches:'
write(6,'(a,/)')' Valid command line switches:' write(6,'(a)') ' --geom (-g, --geometry)'
write(6,'(a)') ' --geom (-g, --geometry)' write(6,'(a)') ' --load (-l, --loadcase)'
write(6,'(a)') ' --load (-l, --loadcase)' write(6,'(a)') ' --workingdir (-w, --wd, --workingdirectory, -d, --directory)'
write(6,'(a)') ' --workingdir (-w, --wd, --workingdirectory, -d, --directory)' write(6,'(a)') ' --restart (-r, --rs)'
write(6,'(a)') ' --restart (-r, --rs)' write(6,'(a)') ' --help (-h)'
write(6,'(a)') ' --help (-h)' write(6,'(/,a)')' -----------------------------------------------------------------------'
write(6,'(/,a)')' -----------------------------------------------------------------------' write(6,'(a)') ' Mandatory arguments:'
write(6,'(a)') ' Mandatory arguments:' write(6,'(/,a)')' --geom PathToGeomFile/NameOfGeom.geom'
write(6,'(/,a)')' --geom PathToGeomFile/NameOfGeom.geom' write(6,'(a)') ' Specifies the location of the geometry definition file,'
write(6,'(a)') ' Specifies the location of the geometry definition file,' write(6,'(a)') ' if no extension is given, .geom will be appended.'
write(6,'(a)') ' if no extension is given, .geom will be appended.' write(6,'(a)') ' "PathToGeomFile" will be the working directory if not specified'
write(6,'(a)') ' "PathToGeomFile" will be the working directory if not specified' write(6,'(a)') ' via --workingdir.'
write(6,'(a)') ' via --workingdir.' write(6,'(a)') ' Make sure the file "material.config" exists in the working'
write(6,'(a)') ' Make sure the file "material.config" exists in the working' write(6,'(a)') ' directory.'
write(6,'(a)') ' directory.' write(6,'(a)') ' For further configuration place "numerics.config"'
write(6,'(a)') ' For further configuration place "numerics.config"' write(6,'(a)')' and "numerics.config" in that directory.'
write(6,'(a)')' and "numerics.config" in that directory.' write(6,'(/,a)')' --load PathToLoadFile/NameOfLoadFile.load'
write(6,'(/,a)')' --load PathToLoadFile/NameOfLoadFile.load' write(6,'(a)') ' Specifies the location of the load case definition file,'
write(6,'(a)') ' Specifies the location of the load case definition file,' write(6,'(a)') ' if no extension is given, .load will be appended.'
write(6,'(a)') ' if no extension is given, .load will be appended.' write(6,'(/,a)')' -----------------------------------------------------------------------'
write(6,'(/,a)')' -----------------------------------------------------------------------' write(6,'(a)') ' Optional arguments:'
write(6,'(a)') ' Optional arguments:' write(6,'(/,a)')' --workingdirectory PathToWorkingDirectory'
write(6,'(/,a)')' --workingdirectory PathToWorkingDirectory' write(6,'(a)') ' Specifies the working directory and overwrites the default'
write(6,'(a)') ' Specifies the working directory and overwrites the default' write(6,'(a)') ' "PathToGeomFile".'
write(6,'(a)') ' "PathToGeomFile".' write(6,'(a)') ' Make sure the file "material.config" exists in the working'
write(6,'(a)') ' Make sure the file "material.config" exists in the working' write(6,'(a)') ' directory.'
write(6,'(a)') ' directory.' write(6,'(a)') ' For further configuration place "numerics.config"'
write(6,'(a)') ' For further configuration place "numerics.config"' write(6,'(a)')' and "numerics.config" in that directory.'
write(6,'(a)')' and "numerics.config" in that directory.' write(6,'(/,a)')' --restart XX'
write(6,'(/,a)')' --restart XX' write(6,'(a)') ' Reads in total increment No. XX-1 and continues to'
write(6,'(a)') ' Reads in total increment No. XX-1 and continues to' write(6,'(a)') ' calculate total increment No. XX.'
write(6,'(a)') ' calculate total increment No. XX.' write(6,'(a)') ' Appends to existing results file '
write(6,'(a)') ' Appends to existing results file ' write(6,'(a)') ' "NameOfGeom_NameOfLoadFile.spectralOut".'
write(6,'(a)') ' "NameOfGeom_NameOfLoadFile.spectralOut".' write(6,'(a)') ' Works only if the restart information for total increment'
write(6,'(a)') ' Works only if the restart information for total increment' write(6,'(a)') ' No. XX-1 is available in the working directory.'
write(6,'(a)') ' No. XX-1 is available in the working directory.' write(6,'(/,a)')' -----------------------------------------------------------------------'
write(6,'(/,a)')' -----------------------------------------------------------------------' write(6,'(a)') ' Help:'
write(6,'(a)') ' Help:' write(6,'(/,a)')' --help'
write(6,'(/,a)')' --help' write(6,'(a,/)')' Prints this message and exits'
write(6,'(a,/)')' Prints this message and exits' call quit(0_pInt) ! normal Termination
call quit(0_pInt) ! normal Termination
endif mainProcess2
case ('-l', '--load', '--loadcase') case ('-l', '--load', '--loadcase')
loadcaseArg = IIO_stringValue(commandLine,chunkPos,i+1_pInt) loadcaseArg = IIO_stringValue(commandLine,chunkPos,i+1_pInt)
case ('-g', '--geom', '--geometry') case ('-g', '--geom', '--geometry')
@ -185,25 +186,23 @@ subroutine DAMASK_interface_init()
geometryFile = getGeometryFile(geometryArg) geometryFile = getGeometryFile(geometryArg)
loadCaseFile = getLoadCaseFile(loadCaseArg) loadCaseFile = getLoadCaseFile(loadCaseArg)
call get_environment_variable('HOSTNAME',hostName)
call get_environment_variable('USER',userName) call get_environment_variable('USER',userName)
mainProcess3: if (worldrank == 0) then error = getHostName(hostName)
write(6,'(a,a)') ' Host name: ', trim(hostName) write(6,'(a,a)') ' Host name: ', trim(hostName)
write(6,'(a,a)') ' User name: ', trim(userName) write(6,'(a,a)') ' User name: ', trim(userName)
write(6,'(a,a)') ' Path separator: ', getPathSep() write(6,'(a,a)') ' Path separator: ', pathSep
write(6,'(a,a)') ' Command line call: ', trim(commandLine) write(6,'(a,a)') ' Command line call: ', trim(commandLine)
if (len(trim(workingDirArg))>0) & if (len(trim(workingDirArg))>0) &
write(6,'(a,a)') ' Working dir argument: ', trim(workingDirArg) write(6,'(a,a)') ' Working dir argument: ', trim(workingDirArg)
write(6,'(a,a)') ' Geometry argument: ', trim(geometryArg) write(6,'(a,a)') ' Geometry argument: ', trim(geometryArg)
write(6,'(a,a)') ' Loadcase argument: ', trim(loadcaseArg) write(6,'(a,a)') ' Loadcase argument: ', trim(loadcaseArg)
write(6,'(a,a)') ' Working directory: ', trim(getSolverWorkingDirectoryName()) write(6,'(a,a)') ' Working directory: ', trim(getSolverWorkingDirectoryName())
write(6,'(a,a)') ' Geometry file: ', trim(geometryFile) write(6,'(a,a)') ' Geometry file: ', trim(geometryFile)
write(6,'(a,a)') ' Loadcase file: ', trim(loadCaseFile) write(6,'(a,a)') ' Loadcase file: ', trim(loadCaseFile)
write(6,'(a,a)') ' Solver job name: ', trim(getSolverJobName()) write(6,'(a,a)') ' Solver job name: ', trim(getSolverJobName())
if (SpectralRestartInc > 1_pInt) & if (SpectralRestartInc > 1_pInt) &
write(6,'(a,i6.6)') ' Restart at increment: ', spectralRestartInc write(6,'(a,i6.6)') ' Restart at increment: ', spectralRestartInc
write(6,'(a,l1,/)') ' Append to result file: ', appendToOutFile write(6,'(a,l1,/)') ' Append to result file: ', appendToOutFile
endif mainProcess3
end subroutine DAMASK_interface_init end subroutine DAMASK_interface_init
@ -222,11 +221,9 @@ character(len=1024) function storeWorkingDirectory(workingDirectoryArg,geometryA
character(len=*), intent(in) :: workingDirectoryArg !< working directory argument character(len=*), intent(in) :: workingDirectoryArg !< working directory argument
character(len=*), intent(in) :: geometryArg !< geometry argument character(len=*), intent(in) :: geometryArg !< geometry argument
character(len=1024) :: cwd character(len=1024) :: cwd
character :: pathSep
logical :: error logical :: error
external :: quit external :: quit
pathSep = getPathSep()
wdGiven: if (len(workingDirectoryArg)>0) then wdGiven: if (len(workingDirectoryArg)>0) then
absolutePath: if (workingDirectoryArg(1:1) == pathSep) then absolutePath: if (workingDirectoryArg(1:1) == pathSep) then
storeWorkingDirectory = workingDirectoryArg storeWorkingDirectory = workingDirectoryArg
@ -262,6 +259,7 @@ end function storeWorkingDirectory
character(len=1024) function getSolverWorkingDirectoryName() character(len=1024) function getSolverWorkingDirectoryName()
implicit none implicit none
getSolverWorkingDirectoryName = workingDirectory getSolverWorkingDirectoryName = workingDirectory
end function getSolverWorkingDirectoryName end function getSolverWorkingDirectoryName
@ -274,10 +272,8 @@ character(len=1024) function getSolverJobName()
implicit none implicit none
integer :: posExt,posSep integer :: posExt,posSep
character :: pathSep
character(len=1024) :: tempString character(len=1024) :: tempString
pathSep = getPathSep()
tempString = geometryFile tempString = geometryFile
posExt = scan(tempString,'.',back=.true.) posExt = scan(tempString,'.',back=.true.)
@ -308,11 +304,9 @@ character(len=1024) function getGeometryFile(geometryParameter)
cwd cwd
integer :: posExt, posSep integer :: posExt, posSep
logical :: error logical :: error
character :: pathSep
external :: quit external :: quit
getGeometryFile = geometryParameter getGeometryFile = geometryParameter
pathSep = getPathSep()
posExt = scan(getGeometryFile,'.',back=.true.) posExt = scan(getGeometryFile,'.',back=.true.)
posSep = scan(getGeometryFile,pathSep,back=.true.) posSep = scan(getGeometryFile,pathSep,back=.true.)
@ -344,11 +338,9 @@ character(len=1024) function getLoadCaseFile(loadCaseParameter)
cwd cwd
integer :: posExt, posSep integer :: posExt, posSep
logical :: error logical :: error
character :: pathSep
external :: quit external :: quit
getLoadCaseFile = loadcaseParameter getLoadCaseFile = loadcaseParameter
pathSep = getPathSep()
posExt = scan(getLoadCaseFile,'.',back=.true.) posExt = scan(getLoadCaseFile,'.',back=.true.)
posSep = scan(getLoadCaseFile,pathSep,back=.true.) posSep = scan(getLoadCaseFile,pathSep,back=.true.)
@ -374,11 +366,8 @@ function rectifyPath(path)
implicit none implicit none
character(len=*) :: path character(len=*) :: path
character(len=len_trim(path)) :: rectifyPath character(len=len_trim(path)) :: rectifyPath
character :: pathSep
integer :: i,j,k,l ! no pInt integer :: i,j,k,l ! no pInt
pathSep = getPathSep()
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! remove /./ from path ! remove /./ from path
l = len_trim(path) l = len_trim(path)
@ -415,10 +404,8 @@ character(len=1024) function makeRelativePath(a,b)
implicit none implicit none
character (len=*) :: a,b character (len=*) :: a,b
character :: pathSep
integer :: i,posLastCommonSlash,remainingSlashes !no pInt integer :: i,posLastCommonSlash,remainingSlashes !no pInt
pathSep = getPathSep()
posLastCommonSlash = 0 posLastCommonSlash = 0
remainingSlashes = 0 remainingSlashes = 0
@ -434,35 +421,6 @@ character(len=1024) function makeRelativePath(a,b)
end function makeRelativePath end function makeRelativePath
!--------------------------------------------------------------------------------------------------
!> @brief counting / and \ in $PATH System variable the character occuring more often is assumed
! to be the path separator
!--------------------------------------------------------------------------------------------------
character function getPathSep()
implicit none
character(len=2048) :: &
path
integer(pInt) :: &
backslash = 0_pInt, &
slash = 0_pInt
integer :: i
call get_environment_variable('PATH',path)
do i=1, len(trim(path))
if (path(i:i)=='/') slash = slash + 1_pInt
if (path(i:i)=='\') backslash = backslash + 1_pInt
enddo
if (backslash>slash) then
getPathSep = '\'
else
getPathSep = '/'
endif
end function getPathSep
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief taken from IO, check IO_stringValue for documentation !> @brief taken from IO, check IO_stringValue for documentation
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------

40
code/system_routines.f90
View File

@ -9,7 +9,8 @@ module system_routines
public :: & public :: &
isDirectory, & isDirectory, &
getCWD getCWD, &
getHostName
interface interface
@ -29,6 +30,14 @@ interface
integer(C_INT),intent(out) :: stat integer(C_INT),intent(out) :: stat
end subroutine getCurrentWorkDir_C end subroutine getCurrentWorkDir_C
subroutine getHostName_C(str, stat) bind(C)
use, intrinsic :: ISO_C_Binding, only: &
C_INT, &
C_CHAR
character(kind=C_CHAR), dimension(1024), intent(out) :: str ! C string is an array
integer(C_INT),intent(out) :: stat
end subroutine getHostName_C
end interface end interface
@ -85,5 +94,34 @@ logical function getCWD(str)
end function getCWD end function getCWD
!--------------------------------------------------------------------------------------------------
!> @brief gets the current host name
!--------------------------------------------------------------------------------------------------
logical function getHostName(str)
use, intrinsic :: ISO_C_Binding, only: &
C_INT, &
C_CHAR, &
C_NULL_CHAR
implicit none
character(len=*), intent(out) :: str
character(kind=C_CHAR), dimension(1024) :: strFixedLength ! C string is an array
integer(C_INT) :: stat
integer :: i
str = repeat('',len(str))
call getHostName_C(strFixedLength,stat)
do i=1,1024 ! copy array components until Null string is found
if (strFixedLength(i) /= C_NULL_CHAR) then
str(i:i)=strFixedLength(i)
else
exit
endif
enddo
getHostName=merge(.True.,.False.,stat /= 0_C_INT)
end function getHostName
end module system_routines end module system_routines

32
examples/ConfigFiles/Phase_Nonlocal_Aluminum.config
View File

@ -8,32 +8,20 @@ plasticity nonlocal
(output) rho_edge (output) rho_edge
(output) rho_screw (output) rho_screw
(output) rho_sgl (output) rho_sgl
(output) rho_sgl_edge
(output) rho_sgl_edge_pos (output) rho_sgl_edge_pos
(output) rho_sgl_edge_neg (output) rho_sgl_edge_neg
(output) rho_sgl_screw
(output) rho_sgl_screw_pos (output) rho_sgl_screw_pos
(output) rho_sgl_screw_neg (output) rho_sgl_screw_neg
(output) rho_sgl_mobile
(output) rho_sgl_edge_mobile
(output) rho_sgl_edge_pos_mobile (output) rho_sgl_edge_pos_mobile
(output) rho_sgl_edge_neg_mobile (output) rho_sgl_edge_neg_mobile
(output) rho_sgl_screw_mobile
(output) rho_sgl_screw_pos_mobile (output) rho_sgl_screw_pos_mobile
(output) rho_sgl_screw_neg_mobile (output) rho_sgl_screw_neg_mobile
(output) rho_sgl_immobile
(output) rho_sgl_edge_immobile
(output) rho_sgl_edge_pos_immobile (output) rho_sgl_edge_pos_immobile
(output) rho_sgl_edge_neg_immobile (output) rho_sgl_edge_neg_immobile
(output) rho_sgl_screw_immobile
(output) rho_sgl_screw_pos_immobile (output) rho_sgl_screw_pos_immobile
(output) rho_sgl_screw_neg_immobile (output) rho_sgl_screw_neg_immobile
(output) rho_dip
(output) rho_dip_edge (output) rho_dip_edge
(output) rho_dip_screw (output) rho_dip_screw
(output) excess_rho
(output) excess_rho_edge
(output) excess_rho_screw
(output) rho_forest (output) rho_forest
(output) delta (output) delta
(output) delta_sgl (output) delta_sgl
@ -47,10 +35,8 @@ plasticity nonlocal
(output) rho_dot_sgl (output) rho_dot_sgl
(output) rho_dot_sgl_mobile (output) rho_dot_sgl_mobile
(output) rho_dot_dip (output) rho_dot_dip
(output) rho_dot_gen
(output) rho_dot_gen_edge (output) rho_dot_gen_edge
(output) rho_dot_gen_screw (output) rho_dot_gen_screw
(output) rho_dot_sgl2dip
(output) rho_dot_sgl2dip_edge (output) rho_dot_sgl2dip_edge
(output) rho_dot_sgl2dip_screw (output) rho_dot_sgl2dip_screw
(output) rho_dot_ann_ath (output) rho_dot_ann_ath
@ -66,24 +52,6 @@ plasticity nonlocal
(output) velocity_edge_neg (output) velocity_edge_neg
(output) velocity_screw_pos (output) velocity_screw_pos
(output) velocity_screw_neg (output) velocity_screw_neg
(output) slipdirection.x
(output) slipdirection.y
(output) slipdirection.z
(output) slipnormal.x
(output) slipnormal.y
(output) slipnormal.z
(output) fluxDensity_edge_pos.x
(output) fluxDensity_edge_pos.y
(output) fluxDensity_edge_pos.z
(output) fluxDensity_edge_neg.x
(output) fluxDensity_edge_neg.y
(output) fluxDensity_edge_neg.z
(output) fluxDensity_screw_pos.x
(output) fluxDensity_screw_pos.y
(output) fluxDensity_screw_pos.z
(output) fluxDensity_screw_neg.x
(output) fluxDensity_screw_neg.y
(output) fluxDensity_screw_neg.z
(output) maximumDipoleHeight_edge (output) maximumDipoleHeight_edge
(output) maximumDipoleHeight_screw (output) maximumDipoleHeight_screw
(output) accumulated_shear (output) accumulated_shear

20
examples/ConfigFiles/Phase_Nonlocal_Nickel.config Executable file → Normal file
View File

@ -14,8 +14,6 @@ plasticity nonlocal
(output) rho_dip_edge (output) rho_dip_edge
(output) rho_dip_screw (output) rho_dip_screw
(output) rho_forest (output) rho_forest
(output) excess_rho_edge
(output) excess_rho_screw
(output) accumulatedshear (output) accumulatedshear
(output) shearrate (output) shearrate
(output) resolvedstress (output) resolvedstress
@ -30,24 +28,6 @@ plasticity nonlocal
(output) rho_dot_edgejogs (output) rho_dot_edgejogs
(output) rho_dot_flux_edge (output) rho_dot_flux_edge
(output) rho_dot_flux_screw (output) rho_dot_flux_screw
(output) slipdirection.x
(output) slipdirection.y
(output) slipdirection.z
(output) slipnormal.x
(output) slipnormal.y
(output) slipnormal.z
(output) fluxdensity_edge_pos.x
(output) fluxdensity_edge_pos.y
(output) fluxdensity_edge_pos.z
(output) fluxdensity_edge_neg.x
(output) fluxdensity_edge_neg.y
(output) fluxdensity_edge_neg.z
(output) fluxdensity_screw_pos.x
(output) fluxdensity_screw_pos.y
(output) fluxdensity_screw_pos.z
(output) fluxdensity_screw_neg.x
(output) fluxdensity_screw_neg.y
(output) fluxdensity_screw_neg.z
lattice_structure fcc lattice_structure fcc
Nslip 12 # number of slip systems per family Nslip 12 # number of slip systems per family

18
examples/ConfigFiles/Phase_Phenopowerlaw_Aluminum.config
View File

@ -7,11 +7,6 @@ plasticity phenopowerlaw
(output) resolvedstress_slip (output) resolvedstress_slip
(output) accumulated_shear_slip (output) accumulated_shear_slip
(output) totalshear (output) totalshear
(output) resistance_twin
(output) shearrate_twin
(output) resolvedstress_twin
(output) accumulated_shear_twin
(output) totalvolfrac_twin
lattice_structure fcc lattice_structure fcc
Nslip 12 # per family Nslip 12 # per family
@ -26,19 +21,6 @@ n_slip 20
tau0_slip 31e6 # per family tau0_slip 31e6 # per family
tausat_slip 63e6 # per family tausat_slip 63e6 # per family
a_slip 2.25 a_slip 2.25
gdot0_twin 0.001
n_twin 20
tau0_twin 31e6 # per family
s_pr 0 # push-up factor for slip saturation due to twinning
twin_b 0
twin_c 0
twin_d 0
twin_e 0
h0_slipslip 75e6 h0_slipslip 75e6
h0_twinslip 0
h0_twintwin 0
interaction_slipslip 1 1 1.4 1.4 1.4 1.4 interaction_slipslip 1 1 1.4 1.4 1.4 1.4
interaction_sliptwin 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
interaction_twinslip 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
interaction_twintwin 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
atol_resistance 1 atol_resistance 1

14
lib/damask/orientation.py
View File

@ -995,11 +995,17 @@ class Orientation:
relationModel, relationModel,
direction, direction,
targetSymmetry = None): targetSymmetry = None):
"""
orientation relationship
positive number: fcc --> bcc
negative number: bcc --> fcc
"""
if relationModel not in ['KS','GT','GTdash','NW','Pitsch','Bain']: return None if relationModel not in ['KS','GT','GTdash','NW','Pitsch','Bain']: return None
if int(direction) == 0: return None if int(direction) == 0: return None
# KS from S. Morito et al./Journal of Alloys and Compounds 5775 (2013) S587-S592 DOES THIS PAPER EXISTS? # KS from S. Morito et al./Journal of Alloys and Compounds 5775 (2013) S587-S592
# for KS rotation matrices also check K. Kitahara et al./Acta Materialia 54 (2006) 1279-1288
# GT from Y. He et al./Journal of Applied Crystallography (2006). 39, 72-81 # GT from Y. He et al./Journal of Applied Crystallography (2006). 39, 72-81
# GT' from Y. He et al./Journal of Applied Crystallography (2006). 39, 72-81 # GT' from Y. He et al./Journal of Applied Crystallography (2006). 39, 72-81
# NW from H. Kitahara et al./Materials Characterization 54 (2005) 378-386 # NW from H. Kitahara et al./Materials Characterization 54 (2005) 378-386
@ -1226,14 +1232,14 @@ class Orientation:
myPlane /= np.linalg.norm(myPlane) myPlane /= np.linalg.norm(myPlane)
myNormal = [float(i) for i in normals[relationModel][variant,me]] # map(float, planes[...]) does not work in python 3 myNormal = [float(i) for i in normals[relationModel][variant,me]] # map(float, planes[...]) does not work in python 3
myNormal /= np.linalg.norm(myNormal) myNormal /= np.linalg.norm(myNormal)
myMatrix = np.array([myPlane,myNormal,np.cross(myPlane,myNormal)]) myMatrix = np.array([myNormal,np.cross(myPlane,myNormal),myPlane]).T
otherPlane = [float(i) for i in planes[relationModel][variant,other]] # map(float, planes[...]) does not work in python 3 otherPlane = [float(i) for i in planes[relationModel][variant,other]] # map(float, planes[...]) does not work in python 3
otherPlane /= np.linalg.norm(otherPlane) otherPlane /= np.linalg.norm(otherPlane)
otherNormal = [float(i) for i in normals[relationModel][variant,other]] # map(float, planes[...]) does not work in python 3 otherNormal = [float(i) for i in normals[relationModel][variant,other]] # map(float, planes[...]) does not work in python 3
otherNormal /= np.linalg.norm(otherNormal) otherNormal /= np.linalg.norm(otherNormal)
otherMatrix = np.array([otherPlane,otherNormal,np.cross(otherPlane,otherNormal)]) otherMatrix = np.array([otherNormal,np.cross(otherPlane,otherNormal),otherPlane]).T
rot=np.dot(otherMatrix.T,myMatrix) rot=np.dot(otherMatrix,myMatrix.T)
return Orientation(matrix=np.dot(rot,self.asMatrix())) # no symmetry information ?? return Orientation(matrix=np.dot(rot,self.asMatrix())) # no symmetry information ??

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

@ -23,6 +23,8 @@ class Test():
'keep': False, 'keep': False,
'accept': False, 'accept': False,
'updateRequest': False, 'updateRequest': False,
'show': False,
'select': None,
} }
for arg in defaults.keys(): for arg in defaults.keys():
setattr(self,arg,kwargs.get(arg) if kwargs.get(arg) else defaults[arg]) setattr(self,arg,kwargs.get(arg) if kwargs.get(arg) else defaults[arg])
@ -58,10 +60,18 @@ class Test():
action = "store_true", action = "store_true",
dest = "accept", dest = "accept",
help = "calculate results but always consider test as successfull") help = "calculate results but always consider test as successfull")
self.parser.add_option("-l", "--list",
action = "store_true",
dest = "show",
help = "show all test variants and do no calculation")
self.parser.add_option("-s", "--select",
dest = "select",
help = "run test of given name only")
self.parser.set_defaults(keep = self.keep, self.parser.set_defaults(keep = self.keep,
accept = self.accept, accept = self.accept,
update = self.updateRequest, update = self.updateRequest,
show = self.show,
select = self.select,
) )
@ -73,21 +83,26 @@ class Test():
self.prepareAll() self.prepareAll()
for variant,name in enumerate(self.variants): for variant,name in enumerate(self.variants):
try: if self.options.show:
if not self.options.keep: logging.critical('{}: {}'.format(variant,name))
self.prepare(variant) elif self.options.select is not None and name != self.options.select:
self.run(variant) pass
else:
try:
if not self.options.keep:
self.prepare(variant)
self.run(variant)
self.postprocess(variant) self.postprocess(variant)
if self.options.update: if self.options.update:
if self.update(variant) != 0: logging.critical('update for "{}" failed.'.format(name)) if self.update(variant) != 0: logging.critical('update for "{}" failed.'.format(name))
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('exception during variant execution: {}'.format(e)) logging.critical('exception during variant execution: {}'.format(e))
return variant+1 # return culprit return variant+1 # return culprit
return 0 return 0
def feasible(self): def feasible(self):

74
processing/misc/DREAM3D_toTable.py Executable file
View File

@ -0,0 +1,74 @@
#!/usr/bin/env python2.7
# -*- coding: UTF-8 no BOM -*-
import os,h5py
import numpy as np
from optparse import OptionParser
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
#--------------------------------------------------------------------------------------------------
# MAIN
#--------------------------------------------------------------------------------------------------
parser = OptionParser(option_class=damask.extendableOption, usage='%prog [geomfile[s]]', description = """
Convert DREAM3D file to ASCIItable
""", version = scriptID)
(options, filenames) = parser.parse_args()
rootDir ='DataContainers/ImageDataContainer'
# --- loop over input files -------------------------------------------------------------------------
if filenames == []: parser.error('no input file specified.')
for name in filenames:
try:
table = damask.ASCIItable(outname = os.path.splitext(name)[0]+'.txt',
buffered = False
)
except: continue
damask.util.report(scriptName,name)
inFile = h5py.File(name, 'r')
grid = inFile[rootDir+'/_SIMPL_GEOMETRY/DIMENSIONS'][...]
# --- read comments --------------------------------------------------------------------------------
coords = (np.mgrid[0:grid[2], 0:grid[1], 0: grid[0]]).reshape(3, -1).T
coords = (np.fliplr(coords)*inFile[rootDir+'/_SIMPL_GEOMETRY/SPACING'][...] \
+ inFile[rootDir+'/_SIMPL_GEOMETRY/ORIGIN'][...] \
+ inFile[rootDir+'/_SIMPL_GEOMETRY/SPACING'][...]*0.5)
table.data = np.hstack( (coords,
inFile[rootDir+'/CellData/EulerAngles'][...].reshape(grid.prod(),3),
inFile[rootDir+'/CellData/Phases'][...].reshape(grid.prod(),1),
inFile[rootDir+'/CellData/Confidence Index'][...].reshape(grid.prod(),1),
inFile[rootDir+'/CellData/Fit'][...].reshape(grid.prod(),1),
inFile[rootDir+'/CellData/Image Quality'][...].reshape(grid.prod(),1)))
labels = ['1_pos','2_pos','3_pos',
'1_Euler','2_Euler','3_Euler',
'PhaseID','CI','Fit','IQ']
try:
table.data = np.hstack((table.data, inFile[rootDir+'/CellData/FeatureIds'][...].reshape(grid.prod(),1)))
labels.append(['FeatureID'])
except Exception:
pass
# ------------------------------------------ assemble header ---------------------------------------
table.labels_clear()
table.labels_append(labels,reset = True)
table.head_write()
# ------------------------------------------ finalize output ---------------------------------------
table.data_writeArray() #(fmt='%e2.2')
table.close()

2
processing/pre/geom_check.sh
View File

@ -11,7 +11,7 @@ do
< $geom \ < $geom \
| \ | \
vtk_addRectilinearGridData \ vtk_addRectilinearGridData \
--scalar microstructure \ --data microstructure \
--inplace \ --inplace \
--vtk ${geom%.*}.vtk --vtk ${geom%.*}.vtk
rm ${geom%.*}.vtk rm ${geom%.*}.vtk

6
processing/pre/geom_grainGrowth.py
View File

@ -89,9 +89,9 @@ for name in filenames:
# Calculates gaussian weights for simulating 3d diffusion # Calculates gaussian weights for simulating 3d diffusion
gauss = np.exp(-(X*X + Y*Y + Z*Z)/(2.0*options.d*options.d))/math.pow(2.0*np.pi*options.d*options.d,1.5) gauss = np.exp(-(X*X + Y*Y + Z*Z)/(2.0*options.d*options.d))/math.pow(2.0*np.pi*options.d*options.d,1.5)
gauss[:,:,grid[2]/2::] = gauss[:,:,round(grid[2]/2.)-1::-1] # trying to cope with uneven (odd) grid size gauss[:,:,grid[2]/2::] = gauss[:,:,int(round(grid[2]/2.))-1::-1] # trying to cope with uneven (odd) grid size
gauss[:,grid[1]/2::,:] = gauss[:,round(grid[1]/2.)-1::-1,:] gauss[:,grid[1]/2::,:] = gauss[:,int(round(grid[1]/2.))-1::-1,:]
gauss[grid[0]/2::,:,:] = gauss[round(grid[0]/2.)-1::-1,:,:] gauss[grid[0]/2::,:,:] = gauss[int(round(grid[0]/2.))-1::-1,:,:]
gauss = np.fft.rfftn(gauss) gauss = np.fft.rfftn(gauss)
interfacialEnergy = lambda A,B: (A*B != 0)*(A != B)*1.0 #1.0 if A & B are distinct & nonzero, 0.0 otherwise interfacialEnergy = lambda A,B: (A*B != 0)*(A != B)*1.0 #1.0 if A & B are distinct & nonzero, 0.0 otherwise

117
processing/pre/geom_mirror.py Executable file
View File

@ -0,0 +1,117 @@
#!/usr/bin/env python2.7
# -*- coding: UTF-8 no BOM -*-
import os,sys,math
import numpy as np
import damask
from optparse import OptionParser
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
#--------------------------------------------------------------------------------------------------
# MAIN
#--------------------------------------------------------------------------------------------------
validDirections = ['x','y','z']
parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [geomfile(s)]', description = """
Mirrors spectral geometry description along given directions.
""", version=scriptID)
parser.add_option('-d','--direction',
dest = 'directions',
action = 'extend', metavar = '<string LIST>',
help = "directions in which to mirror {'x','y','z'}")
(options, filenames) = parser.parse_args()
if options.directions is None:
parser.error('no direction given.')
if not set(options.directions).issubset(validDirections):
invalidDirections = [str(e) for e in set(options.directions).difference(validDirections)]
parser.error('invalid directions {}. '.format(*invalidDirections))
# --- loop over input files -------------------------------------------------------------------------
if filenames == []: filenames = [None]
for name in filenames:
try:
table = damask.ASCIItable(name = name,
buffered = False, labeled = False)
except: continue
damask.util.report(scriptName,name)
# --- interpret header ----------------------------------------------------------------------------
table.head_read()
info,extra_header = table.head_getGeom()
damask.util.croak(['grid a b c: %s'%(' x '.join(map(str,info['grid']))),
'size x y z: %s'%(' x '.join(map(str,info['size']))),
'origin x y z: %s'%(' : '.join(map(str,info['origin']))),
'homogenization: %i'%info['homogenization'],
'microstructures: %i'%info['microstructures'],
])
errors = []
if np.any(info['grid'] < 1): errors.append('invalid grid a b c.')
if np.any(info['size'] <= 0.0): errors.append('invalid size x y z.')
if errors != []:
damask.util.croak(errors)
table.close(dismiss = True)
continue
# --- read data ------------------------------------------------------------------------------------
microstructure = table.microstructure_read(info['grid']).reshape(info['grid'],order='F') # read microstructure
if 'z' in options.directions:
microstructure = np.concatenate([microstructure,microstructure[:,:,::-1]],2)
if 'y' in options.directions:
microstructure = np.concatenate([microstructure,microstructure[:,::-1,:]],1)
if 'x' in options.directions:
microstructure = np.concatenate([microstructure,microstructure[::-1,:,:]],0)
# --- do work ------------------------------------------------------------------------------------
newInfo = {
'size': microstructure.shape*info['size']/info['grid'],
'grid': microstructure.shape,
}
# --- report ---------------------------------------------------------------------------------------
remarks = []
if (any(newInfo['grid'] != info['grid'])):
remarks.append('--> grid a b c: %s'%(' x '.join(map(str,newInfo['grid']))))
if (any(newInfo['size'] != info['size'])):
remarks.append('--> size x y z: %s'%(' x '.join(map(str,newInfo['size']))))
if remarks != []: damask.util.croak(remarks)
# --- write header ---------------------------------------------------------------------------------
table.labels_clear()
table.info_clear()
table.info_append([
scriptID + ' ' + ' '.join(sys.argv[1:]),
"grid\ta {grid[0]}\tb {grid[1]}\tc {grid[2]}".format(grid=newInfo['grid']),
"size\tx {size[0]}\ty {size[1]}\tz {size[2]}".format(size=newInfo['size']),
"origin\tx {origin[0]}\ty {origin[1]}\tz {origin[2]}".format(origin=info['origin']),
"homogenization\t{homog}".format(homog=info['homogenization']),
"microstructures\t{microstructures}".format(microstructures=info['microstructures']),
extra_header
])
table.head_write()
# --- write microstructure information ------------------------------------------------------------
formatwidth = int(math.floor(math.log10(microstructure.max())+1))
table.data = microstructure.reshape((newInfo['grid'][0],np.prod(newInfo['grid'][1:])),order='F').transpose()
table.data_writeArray('%%%ii'%(formatwidth),delimiter = ' ')
# --- output finalization --------------------------------------------------------------------------
table.close() # close ASCII table

184
processing/pre/patchFromReconstructedBoundaries.py
View File

@ -2,6 +2,7 @@
# -*- coding: UTF-8 no BOM -*- # -*- coding: UTF-8 no BOM -*-
import sys,os,math,re import sys,os,math,re
import numpy as np
from optparse import OptionParser from optparse import OptionParser
import damask import damask
@ -63,9 +64,8 @@ def rcbOrientationParser(content,idcolumn):
grains = [] grains = []
myOrientation = [0.0,0.0,0.0] myOrientation = [0.0,0.0,0.0]
for line in content: for j,line in enumerate(content):
m = re.match(r'\s*(#|$)',line) if re.match(r'^\s*(#|$)',line): continue # skip comments and blank lines
if m: continue # skip comments and blank lines
for grain in range(2): for grain in range(2):
myID = int(line.split()[idcolumn+grain]) # get grain id myID = int(line.split()[idcolumn+grain]) # get grain id
myOrientation = map(float,line.split())[3*grain:3+3*grain] # get orientation myOrientation = map(float,line.split())[3*grain:3+3*grain] # get orientation
@ -75,8 +75,8 @@ def rcbOrientationParser(content,idcolumn):
try: try:
grains[myID-1] = myOrientation # store Euler angles grains[myID-1] = myOrientation # store Euler angles
except IndexError: except IndexError:
message = 'You might not have chosen the correct column for the grain IDs! Please check the "--id" option.' damask.util.croak('You might not have chosen the correct column for the grain IDs! '+
print '\033[1;31m'+message+'\033[0m\n' 'Please check the "--id" option.')
raise raise
except: except:
raise raise
@ -91,13 +91,13 @@ def rcbParser(content,M,size,tolerance,idcolumn,segmentcolumn):
x = [0.,0.] x = [0.,0.]
y = [0.,0.] y = [0.,0.]
for line in content: for line in content:
m = re.match(r'\s*(#|$)',line) m = re.match(r'^\s*(#|$)',line)
if m: continue # skip comments and blank lines if m: continue # skip comments and blank lines
try: try:
(x[0],y[0],x[1],y[1]) = map(float,line.split())[segmentcolumn:segmentcolumn+4] # get start and end coordinates of each segment. (x[0],y[0],x[1],y[1]) = map(float,line.split())[segmentcolumn:segmentcolumn+4] # get start and end coordinates of each segment.
except IndexError: except IndexError:
message = 'You might not have chosen the correct column for the segment end points! Please check the "--segment" option.' damask.util.croak('You might not have chosen the correct column for the segment end points! '+
print '\033[1;31m'+message+'\033[0m\n' 'Please check the "--segment" option.')
raise raise
except: except:
raise raise
@ -110,6 +110,9 @@ def rcbParser(content,M,size,tolerance,idcolumn,segmentcolumn):
dX = boxX[1]-boxX[0] dX = boxX[1]-boxX[0]
dY = boxY[1]-boxY[0] dY = boxY[1]-boxY[0]
damask.util.croak(' bounding box {},{} -- {},{}'.format(boxX[0],boxY[0],boxX[1],boxY[1]))
damask.util.croak(' dimension {} x {}'.format(dX,dY))
if size > 0.0: scalePatch = size/dX if size > 0.0: scalePatch = size/dX
else: scalePatch = 1.0 else: scalePatch = 1.0
@ -122,8 +125,7 @@ def rcbParser(content,M,size,tolerance,idcolumn,segmentcolumn):
grainNeighbors = [] grainNeighbors = []
for line in content: for line in content:
m = re.match(r'\s*(#|$)',line) if re.match(r'^\s*(#|$)',line): continue # skip comments and blank lines
if m: continue # skip comments and blank lines
(x[0],y[0],x[1],y[1]) = map(float,line.split())[segmentcolumn:segmentcolumn+4] # get start and end coordinates of each segment. (x[0],y[0],x[1],y[1]) = map(float,line.split())[segmentcolumn:segmentcolumn+4] # get start and end coordinates of each segment.
(x[0],y[0]) = (M[0]*x[0]+M[1]*y[0],M[2]*x[0]+M[3]*y[0]) # apply transformation to coordinates (x[0],y[0]) = (M[0]*x[0]+M[1]*y[0],M[2]*x[0]+M[3]*y[0]) # apply transformation to coordinates
(x[1],y[1]) = (M[0]*x[1]+M[1]*y[1],M[2]*x[1]+M[3]*y[1]) # to get rcb --> Euler system (x[1],y[1]) = (M[0]*x[1]+M[1]*y[1],M[2]*x[1]+M[3]*y[1]) # to get rcb --> Euler system
@ -133,8 +135,8 @@ def rcbParser(content,M,size,tolerance,idcolumn,segmentcolumn):
y[0] -= boxY[0] y[0] -= boxY[0]
y[1] -= boxY[0] y[1] -= boxY[0]
grainNeighbors.append(map(int,line.split()[idcolumn:idcolumn+2])) # remember right and left grain per segment grainNeighbors.append(map(int,line.split()[idcolumn:idcolumn+2])) # remember right and left grain per segment
for i in range(2): # store segment to both points for i in range(2): # store segment to both points
match = False # check whether point is already known (within a small range) match = False # check whether point is already known (within a small range)
for posX in connectivityXY.keys(): for posX in connectivityXY.keys():
if (abs(float(posX)-x[i])<dX*tolerance): if (abs(float(posX)-x[i])<dX*tolerance):
for posY in connectivityXY[posX].keys(): for posY in connectivityXY[posX].keys():
@ -144,7 +146,7 @@ def rcbParser(content,M,size,tolerance,idcolumn,segmentcolumn):
match = True match = True
break break
break break
# force to boundary if inside tolerance to it # force onto boundary if inside tolerance to it
if (not match): if (not match):
if (abs(x[i])<dX*tolerance): if (abs(x[i])<dX*tolerance):
x[i] = 0 x[i] = 0
@ -168,7 +170,6 @@ def rcbParser(content,M,size,tolerance,idcolumn,segmentcolumn):
connectivityYX[keyY][keyX].append(segment) connectivityYX[keyY][keyX].append(segment)
segment += 1 segment += 1
# top border # top border
keyId = "0" keyId = "0"
boundary = connectivityYX[keyId].keys() boundary = connectivityYX[keyId].keys()
@ -225,17 +226,36 @@ def rcbParser(content,M,size,tolerance,idcolumn,segmentcolumn):
for keyY in allkeysY: for keyY in allkeysY:
points.append({'coords': [float(keyX)*scalePatch,float(keyY)*scalePatch], 'segments': connectivityXY[keyX][keyY]}) points.append({'coords': [float(keyX)*scalePatch,float(keyY)*scalePatch], 'segments': connectivityXY[keyX][keyY]})
for segment in connectivityXY[keyX][keyY]: for segment in connectivityXY[keyX][keyY]:
if (segments[segment] is None): segments[segment].append(pointId)
segments[segment] = pointId
else:
segments[segment].append(pointId)
pointId += 1 pointId += 1
dupSegments = []
for pointId,point in enumerate(points):
ends = []
goners = []
for segment in point['segments']:
end = segments[segment][1 if segments[segment][0] == pointId else 0]
if end in ends:
goners.append(segment)
dupSegments.append(segment)
else:
ends.append(end)
for item in goners:
point['segments'].remove(item)
if len(dupSegments) > 0:
damask.util.croak(' culling {} duplicate segments...'.format(len(dupSegments)))
for rm in dupSegments:
segments[rm] = None
crappyData = False crappyData = False
for pointId,point in enumerate(points): for pointId,point in enumerate(points):
if len(point['segments']) < 2: # point marks a dead end! if len(point['segments']) < 2: # point marks a dead end!
print "Dead end at segment %i (%f,%f)"\ damask.util.croak('dead end at segment {} for point {} ({},{}).'
%(1+point['segments'][0],boxX[0]+point['coords'][0]/scalePatch,boxY[0]+point['coords'][1]/scalePatch,) .format(point['segments'][0],
pointId,
boxX[0]+point['coords'][0]/scalePatch,boxY[0]+point['coords'][1]/scalePatch,))
crappyData = True crappyData = True
grains = {'draw': [], 'legs': []} grains = {'draw': [], 'legs': []}
@ -249,58 +269,69 @@ def rcbParser(content,M,size,tolerance,idcolumn,segmentcolumn):
innerAngleSum = 0.0 innerAngleSum = 0.0
myWalk = point['segments'].pop() myWalk = point['segments'].pop()
grainLegs = [myWalk] grainLegs = [myWalk]
if segments[myWalk][0] == myStart: myEnd = segments[myWalk][1 if segments[myWalk][0] == myStart else 0]
myEnd = segments[myWalk][1]
else:
myEnd = segments[myWalk][0]
while (myEnd != pointId): while (myEnd != pointId):
myV = [points[myEnd]['coords'][0]-points[myStart]['coords'][0],\ myV = [points[myEnd]['coords'][0]-points[myStart]['coords'][0],
points[myEnd]['coords'][1]-points[myStart]['coords'][1]] points[myEnd]['coords'][1]-points[myStart]['coords'][1]]
myLen = math.sqrt(myV[0]**2+myV[1]**2) myLen = math.sqrt(myV[0]**2+myV[1]**2)
if myLen == 0.0: damask.util.croak('mylen is zero: point {} --> {}'.format(myStart,myEnd))
best = {'product': -2.0, 'peek': -1, 'len': -1, 'point': -1} best = {'product': -2.0, 'peek': -1, 'len': -1, 'point': -1}
for peek in points[myEnd]['segments']: # trying in turn all segments emanating from current end for peek in points[myEnd]['segments']: # trying in turn all segments emanating from current end
if peek == myWalk: if peek == myWalk:
continue continue # do not go back same path
peekEnd = segments[peek][1] if segments[peek][0] == myEnd else segments[peek][0] peekEnd = segments[peek][1 if segments[peek][0] == myEnd else 0]
peekV = [points[myEnd]['coords'][0]-points[peekEnd]['coords'][0], peekV = [points[peekEnd]['coords'][0]-points[myEnd]['coords'][0],
points[myEnd]['coords'][1]-points[peekEnd]['coords'][1]] points[peekEnd]['coords'][1]-points[myEnd]['coords'][1]]
peekLen = math.sqrt(peekV[0]**2+peekV[1]**2) peekLen = math.sqrt(peekV[0]**2+peekV[1]**2)
crossproduct = (myV[0]*peekV[1]-myV[1]*peekV[0])/myLen/peekLen if peekLen == 0.0: damask.util.croak('peeklen is zero: peek point {}'.format(peek))
dotproduct = (myV[0]*peekV[0]+myV[1]*peekV[1])/myLen/peekLen crossproduct = (myV[0]*peekV[1] - myV[1]*peekV[0])/myLen/peekLen
if crossproduct*(dotproduct+1.0) >= best['product']: dotproduct = (myV[0]*peekV[0] + myV[1]*peekV[1])/myLen/peekLen
best['product'] = crossproduct*(dotproduct+1.0) innerAngle = math.copysign(1.0,crossproduct)*(dotproduct-1.0)
if innerAngle >= best['product']: # takes sharpest left turn
best['product'] = innerAngle
best['peek'] = peek best['peek'] = peek
best['point'] = peekEnd best['point'] = peekEnd
innerAngleSum += best['product'] innerAngleSum += best['product']
myWalk = best['peek'] myWalk = best['peek']
myStart = myEnd myStart = myEnd
myEnd = best['point'] myEnd = best['point']
if myWalk in points[myStart]['segments']: if myWalk in points[myStart]['segments']:
points[myStart]['segments'].remove(myWalk) points[myStart]['segments'].remove(myWalk)
else: else:
sys.stderr.write(str(myWalk)+' not in segments of point '+str(myStart)+'\n') damask.util.croak('{} not in segments of point {}'.format(myWalk,myStart))
grainDraw.append(points[myStart]['coords']) grainDraw.append(points[myStart]['coords'])
grainLegs.append(myWalk) grainLegs.append(myWalk)
if innerAngleSum > 0.0: if innerAngleSum > 0.0:
grains['draw'].append(grainDraw) grains['draw'].append(grainDraw)
grains['legs'].append(grainLegs) grains['legs'].append(grainLegs)
else: else:
grains['box'] = grainLegs grains['box'] = grainLegs
# build overall data structure # build overall data structure
rcData = {'dimension':[dX,dY], 'point': [],'segment': [], 'grain': [], 'grainMapping': []} rcData = {'dimension':[dX,dY],
print " dimension %g x %g"%(dX,dY) 'bounds': [[boxX[0],boxY[0]],[boxX[1],boxY[1]]],
'scale': scalePatch,
'point': [],
'segment': [],
'neighbors': [],
'grain': [],
'grainMapping': [],
}
for point in points: for point in points:
rcData['point'].append(point['coords']) rcData['point'].append(point['coords'])
print " found %i points"%(len(rcData['point'])) damask.util.croak(' found {} points'.format(len(rcData['point'])))
for segment in segments: for segment in segments:
rcData['segment'].append(segment) rcData['segment'].append(segment)
print " built %i segments"%(len(rcData['segment'])) damask.util.croak(' built {} segments'.format(len(rcData['segment'])))
for neighbors in grainNeighbors:
rcData['neighbors'].append(neighbors)
for legs in grains['legs']: # loop over grains for legs in grains['legs']: # loop over grains
rcData['grain'].append(legs) # store list of boundary segments rcData['grain'].append(legs) # store list of boundary segments
@ -314,12 +345,11 @@ def rcbParser(content,M,size,tolerance,idcolumn,segmentcolumn):
myNeighbors[grainNeighbors[leg][side]] = 1 myNeighbors[grainNeighbors[leg][side]] = 1
if myNeighbors: # do I have any neighbors (i.e., non-bounding box segment) if myNeighbors: # do I have any neighbors (i.e., non-bounding box segment)
candidateGrains = sorted(myNeighbors.iteritems(), key=lambda (k,v): (v,k), reverse=True) # sort grain counting candidateGrains = sorted(myNeighbors.iteritems(), key=lambda (k,v): (v,k), reverse=True) # sort grain counting
if candidateGrains[0][0] not in rcData['grainMapping']: # most frequent one not yet seen? # most frequent one not yet seen?
rcData['grainMapping'].append(candidateGrains[0][0]) # must be me then rcData['grainMapping'].append(candidateGrains[0 if candidateGrains[0][0] not in rcData['grainMapping'] else 1][0]) # must be me then
else: # special case of bi-crystal situation...
rcData['grainMapping'].append(candidateGrains[1][0]) # special case of bi-crystal situation...
print " found %i grains\n"%(len(rcData['grain'])) damask.util.croak(' found {} grains'.format(len(rcData['grain'])))
rcData['box'] = grains['box'] if 'box' in grains else [] rcData['box'] = grains['box'] if 'box' in grains else []
@ -670,9 +700,10 @@ def image(name,imgsize,marginX,marginY,rcData):
draw.text([offsetX+point[0]*scaleImg,sizeY-(offsetY+point[1]*scaleImg)],"%i"%id,fill=(0,0,0)) draw.text([offsetX+point[0]*scaleImg,sizeY-(offsetY+point[1]*scaleImg)],"%i"%id,fill=(0,0,0))
for id,vertex in enumerate(rcData['segment']): for id,vertex in enumerate(rcData['segment']):
if vertex:
start = rcData['point'][vertex[0]] start = rcData['point'][vertex[0]]
end = rcData['point'][vertex[1]] end = rcData['point'][vertex[1]]
draw.text([offsetX+(start[0]+end[0])/2.0*scaleImg,sizeY-(offsetY+(start[1]+end[1])/2.0*scaleImg)],"%i"%id,fill=(0,0,128)) draw.text([offsetX+(start[0]+end[0])/2.0*scaleImg,sizeY-(offsetY+(start[1]+end[1])/2.0*scaleImg)],"%i"%id,fill=(255,0,128))
draw.line([offsetX+start[0]*scaleImg,sizeY-(offsetY+start[1]*scaleImg), draw.line([offsetX+start[0]*scaleImg,sizeY-(offsetY+start[1]*scaleImg),
offsetX+ end[0]*scaleImg,sizeY-(offsetY+ end[1]*scaleImg)],width=1,fill=(128,128,128)) offsetX+ end[0]*scaleImg,sizeY-(offsetY+ end[1]*scaleImg)],width=1,fill=(128,128,128))
@ -739,8 +770,8 @@ def fftbuild(rcData,height,xframe,yframe,resolution,extrusion):
'dimension':(xsize,ysize,zsize)} 'dimension':(xsize,ysize,zsize)}
frameindex=len(rcData['grain'])+1 # calculate frame index as largest grain index plus one frameindex=len(rcData['grain'])+1 # calculate frame index as largest grain index plus one
dx = xsize/(xres+1) # calculate step sizes dx = xsize/(xres) # calculate step sizes
dy = ysize/(yres+1) dy = ysize/(yres)
grainpoints = [] grainpoints = []
for segments in rcData['grain']: # get segments of each grain for segments in rcData['grain']: # get segments of each grain
@ -755,11 +786,10 @@ def fftbuild(rcData,height,xframe,yframe,resolution,extrusion):
grainpoints.append([]) # start out blank for current grain grainpoints.append([]) # start out blank for current grain
for p in sorted(points, key=points.get): # loop thru set of sorted points for p in sorted(points, key=points.get): # loop thru set of sorted points
grainpoints[-1].append([rcData['point'][p][0],rcData['point'][p][1]]) # append x,y of point grainpoints[-1].append([rcData['point'][p][0],rcData['point'][p][1]]) # append x,y of point
bestGuess = 0 # assume grain 0 as best guess bestGuess = 0 # assume grain 0 as best guess
for i in range(int(xres*yres)): # walk through all points in xy plane for i in range(int(xres*yres)): # walk through all points in xy plane
xtest = -xframe+((i%xres)+0.5)*dx # calculate coordinates xtest = -xframe+((i%xres)+0.5)*dx # calculate coordinates
ytest = -yframe+(int(i/xres)+0.5)*dy ytest = -yframe+((i//xres)+0.5)*dy
if(xtest < 0 or xtest > maxX): # check wether part of frame if(xtest < 0 or xtest > maxX): # check wether part of frame
if( ytest < 0 or ytest > maxY): # part of edges if( ytest < 0 or ytest > maxY): # part of edges
fftdata['fftpoints'].append(frameindex+2) # append frameindex to result array fftdata['fftpoints'].append(frameindex+2) # append frameindex to result array
@ -789,7 +819,7 @@ reconstructed boundary file
meshes=['dt_planar_trimesh','af_planar_trimesh','af_planar_quadmesh'] meshes=['dt_planar_trimesh','af_planar_trimesh','af_planar_quadmesh']
parser.add_option('-o', '--output', action='extend', dest='output', metavar = '<string LIST>', parser.add_option('-o', '--output', action='extend', dest='output', metavar = '<string LIST>',
help='types of output {image, mentat, procedure, spectral}') help='types of output {rcb, image, mentat, procedure, spectral}')
parser.add_option('-p', '--port', type='int', metavar = 'int', parser.add_option('-p', '--port', type='int', metavar = 'int',
dest='port', help='Mentat connection port [%default]') dest='port', help='Mentat connection port [%default]')
parser.add_option('-2', '--twodimensional', action='store_true', parser.add_option('-2', '--twodimensional', action='store_true',
@ -847,16 +877,15 @@ parser.set_defaults(output = [],
(options, args) = parser.parse_args() (options, args) = parser.parse_args()
print '\033[1m'+scriptName+'\033[0m\n'
if not len(args): if not len(args):
parser.error('no boundary file specified') parser.error('no boundary file specified.')
try: try:
boundaryFile = open(args[0]) boundaryFile = open(args[0])
boundarySegments = boundaryFile.readlines() boundarySegments = boundaryFile.readlines()
boundaryFile.close() boundaryFile.close()
except: except:
print 'unable to read boundary file "%s"'%args[0] damask.util.croak('unable to read boundary file "{}".'.format(args[0]))
raise raise
options.output = [s.lower() for s in options.output] # lower case options.output = [s.lower() for s in options.output] # lower case
@ -864,18 +893,46 @@ options.idcolumn -= 1 # py
options.segmentcolumn -= 1 # python indexing starts with 0 options.segmentcolumn -= 1 # python indexing starts with 0
myName = os.path.splitext(args[0])[0] myName = os.path.splitext(args[0])[0]
print "%s\n"%myName damask.util.report(scriptName,myName)
orientationData = rcbOrientationParser(boundarySegments,options.idcolumn) orientationData = rcbOrientationParser(boundarySegments,options.idcolumn)
rcData = rcbParser(boundarySegments,options.M,options.size,options.tolerance,options.idcolumn,options.segmentcolumn) rcData = rcbParser(boundarySegments,options.M,options.size,options.tolerance,options.idcolumn,options.segmentcolumn)
# ----- write corrected RCB -----
Minv = np.linalg.inv(np.array(options.M).reshape(2,2))
if 'rcb' in options.output:
print """# Header:
#
# Column 1-3: right hand average orientation (phi1, PHI, phi2 in radians)
# Column 4-6: left hand average orientation (phi1, PHI, phi2 in radians)
# Column 7: length (in microns)
# Column 8: trace angle (in degrees)
# Column 9-12: x,y coordinates of endpoints (in microns)
# Column 13-14: IDs of right hand and left hand grains"""
for i,(left,right) in enumerate(rcData['neighbors']):
if rcData['segment'][i]:
first = np.dot(Minv,np.array([rcData['bounds'][0][0]+rcData['point'][rcData['segment'][i][0]][0]/rcData['scale'],
rcData['bounds'][0][1]+rcData['point'][rcData['segment'][i][0]][1]/rcData['scale'],
]))
second = np.dot(Minv,np.array([rcData['bounds'][0][0]+rcData['point'][rcData['segment'][i][1]][0]/rcData['scale'],
rcData['bounds'][0][1]+rcData['point'][rcData['segment'][i][1]][1]/rcData['scale'],
]))
print ' '.join(map(str,orientationData[left-1]+orientationData[right-1])),
print np.linalg.norm(first-second),
print '0',
print ' '.join(map(str,first)),
print ' '.join(map(str,second)),
print ' '.join(map(str,[left,right]))
# ----- write image ----- # ----- write image -----
if 'image' in options.output and options.imgsize > 0: if 'image' in options.output and options.imgsize > 0:
if ImageCapability: if ImageCapability:
image(myName,options.imgsize,options.xmargin,options.ymargin,rcData) image(myName,options.imgsize,options.xmargin,options.ymargin,rcData)
else: else:
print '...no image drawing possible (PIL missing)...' damask.util.croak('...no image drawing possible (PIL missing)...')
# ----- write spectral geom ----- # ----- write spectral geom -----
@ -884,8 +941,8 @@ if 'spectral' in options.output:
geomFile = open(myName+'_'+str(int(fftdata['resolution'][0]))+'.geom','w') # open geom file for writing geomFile = open(myName+'_'+str(int(fftdata['resolution'][0]))+'.geom','w') # open geom file for writing
geomFile.write('3\theader\n') # write header info geomFile.write('3\theader\n') # write header info
geomFile.write('resolution a %i b %i c %i\n'%(fftdata['resolution'])) # resolution geomFile.write('grid a %i b %i c %i\n'%(fftdata['resolution'])) # grid resolution
geomFile.write('dimension x %f y %f z %f\n'%(fftdata['dimension'])) # size geomFile.write('size x %f y %f z %f\n'%(fftdata['dimension'])) # size
geomFile.write('homogenization 1\n') # homogenization geomFile.write('homogenization 1\n') # homogenization
for z in xrange(fftdata['resolution'][2]): # z repetions for z in xrange(fftdata['resolution'][2]): # z repetions
for y in xrange(fftdata['resolution'][1]): # each x-row separately for y in xrange(fftdata['resolution'][1]): # each x-row separately
@ -893,8 +950,9 @@ if 'spectral' in options.output:
(y+1)*fftdata['resolution'][0]]))+'\n') # grain indexes, x-row per line (y+1)*fftdata['resolution'][0]]))+'\n') # grain indexes, x-row per line
geomFile.close() # close geom file geomFile.close() # close geom file
print('assigned %i out of %i (2D) Fourier points.'\ damask.util.croak('assigned {} out of {} (2D) Fourier points...'
%(len(fftdata['fftpoints']), int(fftdata['resolution'][0])*int(fftdata['resolution'][1]))) .format(len(fftdata['fftpoints']),
int(fftdata['resolution'][0])*int(fftdata['resolution'][1])))
# ----- write Mentat procedure ----- # ----- write Mentat procedure -----
@ -936,7 +994,7 @@ if 'mentat' in options.output:
if 'procedure' in options.output: if 'procedure' in options.output:
output(outputLocals['log'],outputLocals,'Stdout') output(outputLocals['log'],outputLocals,'Stdout')
else: else:
print '...no interaction with Mentat possible...' damask.util.croak('...no interaction with Mentat possible...')
# ----- write config data to file ----- # ----- write config data to file -----

2
processing/pre/seeds_check.sh
View File

@ -5,7 +5,7 @@ do
vtk_pointcloud $seeds vtk_pointcloud $seeds
vtk_addPointcloudData $seeds \ vtk_addPointcloudData $seeds \
--scalar microstructure,weight \ --data microstructure,weight \
--inplace \ --inplace \
--vtk ${seeds%.*}.vtp \ --vtk ${seeds%.*}.vtp \