diff --git a/VERSION b/VERSION index 2a89ebb42..37c33ab3c 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.1-130-gf3308db +v2.0.1-150-g5345b42 diff --git a/code/C_routines.c b/code/C_routines.c index 7be6264d7..5bc09745f 100644 --- a/code/C_routines.c +++ b/code/C_routines.c @@ -9,6 +9,14 @@ /* 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 ){ char cwd_tmp[1024]; 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; - if(stat(dir, &statbuf) != 0) - return 0; - return S_ISDIR(statbuf.st_mode); + +void gethostname_c(char hostname[], int *stat ){ + char hostname_tmp[1024]; + if(gethostname(hostname_tmp, sizeof(hostname_tmp)) == 0){ + strcpy(hostname,hostname_tmp); + *stat = 0; + } + else{ + *stat = 1; + } } diff --git a/code/IO.f90 b/code/IO.f90 index f5db5c1f1..060fcf01f 100644 --- a/code/IO.f90 +++ b/code/IO.f90 @@ -1568,8 +1568,6 @@ subroutine IO_error(error_ID,el,ip,g,ext_msg) msg = 'math_check: R*v == q*v failed' case (410_pInt) msg = 'eigenvalues computation error' - case (450_pInt) - msg = 'unknown symmetry type specified' !------------------------------------------------------------------------------------------------- ! homogenization errors diff --git a/code/math.f90 b/code/math.f90 index d6d26a75d..315da2642 100644 --- a/code/math.f90 +++ b/code/math.f90 @@ -1450,7 +1450,7 @@ pure function math_EulerToQ(eulerangles) cos(halfangles(1)-halfangles(3)) * s, & sin(halfangles(1)-halfangles(3)) * s, & 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 @@ -1508,7 +1508,7 @@ pure function math_EulerAxisAngleToR(axis,omega) real(pReal), dimension(3), intent(in) :: axis 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 @@ -1527,7 +1527,7 @@ pure function math_EulerAxisAngleToQ(axis,omega) real(pReal), dimension(3), intent(in) :: axis 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 @@ -1550,7 +1550,7 @@ pure function math_axisAngleToQ(axis,omega) norm = sqrt(math_mul3x3(axis,axis)) 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)] else rotation 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) 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,3) :: math_symmetricEulers - integer(pInt) :: i,j - math_symmetricEulers(1,1) = PI+Euler(1) - math_symmetricEulers(2,1) = Euler(2) - math_symmetricEulers(3,1) = Euler(3) + math_symmetricEulers = transpose(reshape([PI+Euler(1), PI-Euler(1), 2.0_pReal*PI-Euler(1), & + Euler(2), PI-Euler(2), PI -Euler(2), & + 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(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) + math_symmetricEulers = modulo(math_symmetricEulers,2.0_pReal*pi) 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 - case default ! return blank + case default ! triclinic: return blank math_symmetricEulers = 0.0_pReal end select diff --git a/code/spectral_interface.f90 b/code/spectral_interface.f90 index d49a54411..85ba30c51 100644 --- a/code/spectral_interface.f90 +++ b/code/spectral_interface.f90 @@ -20,6 +20,7 @@ module DAMASK_interface geometryFile = '', & !< parameter given for geometry file loadCaseFile = '' !< parameter given for load case file character(len=1024), private :: workingDirectory !< accessed by getSolverWorkingDirectoryName for compatibility reasons + character, private,parameter :: pathSep = '/' public :: & getSolverWorkingDirectoryName, & @@ -31,7 +32,6 @@ module DAMASK_interface getLoadCaseFile, & rectifyPath, & makeRelativePath, & - getPathSep, & IIO_stringValue, & IIO_intValue, & IIO_lc, & @@ -44,6 +44,8 @@ contains !-------------------------------------------------------------------------------------------------- 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 system_routines, only: & + getHostName implicit none character(len=1024) :: & @@ -64,6 +66,7 @@ subroutine DAMASK_interface_init() integer, dimension(8) :: & dateAndTime ! type default integer PetscErrorCode :: ierr + logical :: error external :: & quit,& MPI_Comm_rank,& @@ -116,54 +119,52 @@ subroutine DAMASK_interface_init() tag = IIO_lc(IIO_stringValue(commandLine,chunkPos,i)) ! extract key select case(tag) case ('-h','--help') - mainProcess2: if (worldrank == 0) then - write(6,'(a)') ' #######################################################################' - write(6,'(a)') ' DAMASK_spectral:' - write(6,'(a)') ' The spectral method boundary value problem solver for' - write(6,'(a)') ' the Düsseldorf Advanced Material Simulation Kit' - write(6,'(a,/)')' #######################################################################' - write(6,'(a,/)')' Valid command line switches:' - write(6,'(a)') ' --geom (-g, --geometry)' - write(6,'(a)') ' --load (-l, --loadcase)' - write(6,'(a)') ' --workingdir (-w, --wd, --workingdirectory, -d, --directory)' - write(6,'(a)') ' --restart (-r, --rs)' - write(6,'(a)') ' --help (-h)' - write(6,'(/,a)')' -----------------------------------------------------------------------' - write(6,'(a)') ' Mandatory arguments:' - write(6,'(/,a)')' --geom PathToGeomFile/NameOfGeom.geom' - 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)') ' "PathToGeomFile" will be the working directory if not specified' - write(6,'(a)') ' via --workingdir.' - write(6,'(a)') ' Make sure the file "material.config" exists in the working' - write(6,'(a)') ' directory.' - write(6,'(a)') ' For further configuration place "numerics.config"' - write(6,'(a)')' and "numerics.config" in that directory.' - write(6,'(/,a)')' --load PathToLoadFile/NameOfLoadFile.load' - 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)')' -----------------------------------------------------------------------' - write(6,'(a)') ' Optional arguments:' - write(6,'(/,a)')' --workingdirectory PathToWorkingDirectory' - write(6,'(a)') ' Specifies the working directory and overwrites the default' - write(6,'(a)') ' "PathToGeomFile".' - write(6,'(a)') ' Make sure the file "material.config" exists in the working' - write(6,'(a)') ' directory.' - write(6,'(a)') ' For further configuration place "numerics.config"' - write(6,'(a)')' and "numerics.config" in that directory.' - write(6,'(/,a)')' --restart XX' - write(6,'(a)') ' Reads in total increment No. XX-1 and continues to' - write(6,'(a)') ' calculate total increment No. XX.' - write(6,'(a)') ' Appends to existing results file ' - write(6,'(a)') ' "NameOfGeom_NameOfLoadFile.spectralOut".' - 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)')' -----------------------------------------------------------------------' - write(6,'(a)') ' Help:' - write(6,'(/,a)')' --help' - write(6,'(a,/)')' Prints this message and exits' - call quit(0_pInt) ! normal Termination - endif mainProcess2 + write(6,'(a)') ' #######################################################################' + write(6,'(a)') ' DAMASK_spectral:' + write(6,'(a)') ' The spectral method boundary value problem solver for' + write(6,'(a)') ' the Düsseldorf Advanced Material Simulation Kit' + write(6,'(a,/)')' #######################################################################' + write(6,'(a,/)')' Valid command line switches:' + write(6,'(a)') ' --geom (-g, --geometry)' + write(6,'(a)') ' --load (-l, --loadcase)' + write(6,'(a)') ' --workingdir (-w, --wd, --workingdirectory, -d, --directory)' + write(6,'(a)') ' --restart (-r, --rs)' + write(6,'(a)') ' --help (-h)' + write(6,'(/,a)')' -----------------------------------------------------------------------' + write(6,'(a)') ' Mandatory arguments:' + write(6,'(/,a)')' --geom PathToGeomFile/NameOfGeom.geom' + 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)') ' "PathToGeomFile" will be the working directory if not specified' + write(6,'(a)') ' via --workingdir.' + write(6,'(a)') ' Make sure the file "material.config" exists in the working' + write(6,'(a)') ' directory.' + write(6,'(a)') ' For further configuration place "numerics.config"' + write(6,'(a)')' and "numerics.config" in that directory.' + write(6,'(/,a)')' --load PathToLoadFile/NameOfLoadFile.load' + 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)')' -----------------------------------------------------------------------' + write(6,'(a)') ' Optional arguments:' + write(6,'(/,a)')' --workingdirectory PathToWorkingDirectory' + write(6,'(a)') ' Specifies the working directory and overwrites the default' + write(6,'(a)') ' "PathToGeomFile".' + write(6,'(a)') ' Make sure the file "material.config" exists in the working' + write(6,'(a)') ' directory.' + write(6,'(a)') ' For further configuration place "numerics.config"' + write(6,'(a)')' and "numerics.config" in that directory.' + write(6,'(/,a)')' --restart XX' + write(6,'(a)') ' Reads in total increment No. XX-1 and continues to' + write(6,'(a)') ' calculate total increment No. XX.' + write(6,'(a)') ' Appends to existing results file ' + write(6,'(a)') ' "NameOfGeom_NameOfLoadFile.spectralOut".' + 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)')' -----------------------------------------------------------------------' + write(6,'(a)') ' Help:' + write(6,'(/,a)')' --help' + write(6,'(a,/)')' Prints this message and exits' + call quit(0_pInt) ! normal Termination case ('-l', '--load', '--loadcase') loadcaseArg = IIO_stringValue(commandLine,chunkPos,i+1_pInt) case ('-g', '--geom', '--geometry') @@ -185,25 +186,23 @@ subroutine DAMASK_interface_init() geometryFile = getGeometryFile(geometryArg) loadCaseFile = getLoadCaseFile(loadCaseArg) - call get_environment_variable('HOSTNAME',hostName) call get_environment_variable('USER',userName) - mainProcess3: if (worldrank == 0) then - write(6,'(a,a)') ' Host name: ', trim(hostName) - write(6,'(a,a)') ' User name: ', trim(userName) - write(6,'(a,a)') ' Path separator: ', getPathSep() - write(6,'(a,a)') ' Command line call: ', trim(commandLine) - if (len(trim(workingDirArg))>0) & - write(6,'(a,a)') ' Working dir argument: ', trim(workingDirArg) - write(6,'(a,a)') ' Geometry argument: ', trim(geometryArg) - write(6,'(a,a)') ' Loadcase argument: ', trim(loadcaseArg) - write(6,'(a,a)') ' Working directory: ', trim(getSolverWorkingDirectoryName()) - write(6,'(a,a)') ' Geometry file: ', trim(geometryFile) - write(6,'(a,a)') ' Loadcase file: ', trim(loadCaseFile) - write(6,'(a,a)') ' Solver job name: ', trim(getSolverJobName()) - if (SpectralRestartInc > 1_pInt) & - write(6,'(a,i6.6)') ' Restart at increment: ', spectralRestartInc - write(6,'(a,l1,/)') ' Append to result file: ', appendToOutFile - endif mainProcess3 + error = getHostName(hostName) + write(6,'(a,a)') ' Host name: ', trim(hostName) + write(6,'(a,a)') ' User name: ', trim(userName) + write(6,'(a,a)') ' Path separator: ', pathSep + write(6,'(a,a)') ' Command line call: ', trim(commandLine) + if (len(trim(workingDirArg))>0) & + write(6,'(a,a)') ' Working dir argument: ', trim(workingDirArg) + write(6,'(a,a)') ' Geometry argument: ', trim(geometryArg) + write(6,'(a,a)') ' Loadcase argument: ', trim(loadcaseArg) + write(6,'(a,a)') ' Working directory: ', trim(getSolverWorkingDirectoryName()) + write(6,'(a,a)') ' Geometry file: ', trim(geometryFile) + write(6,'(a,a)') ' Loadcase file: ', trim(loadCaseFile) + write(6,'(a,a)') ' Solver job name: ', trim(getSolverJobName()) + if (SpectralRestartInc > 1_pInt) & + write(6,'(a,i6.6)') ' Restart at increment: ', spectralRestartInc + write(6,'(a,l1,/)') ' Append to result file: ', appendToOutFile 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) :: geometryArg !< geometry argument character(len=1024) :: cwd - character :: pathSep logical :: error external :: quit - pathSep = getPathSep() wdGiven: if (len(workingDirectoryArg)>0) then absolutePath: if (workingDirectoryArg(1:1) == pathSep) then storeWorkingDirectory = workingDirectoryArg @@ -262,6 +259,7 @@ end function storeWorkingDirectory character(len=1024) function getSolverWorkingDirectoryName() implicit none + getSolverWorkingDirectoryName = workingDirectory end function getSolverWorkingDirectoryName @@ -274,10 +272,8 @@ character(len=1024) function getSolverJobName() implicit none integer :: posExt,posSep - character :: pathSep character(len=1024) :: tempString - pathSep = getPathSep() tempString = geometryFile posExt = scan(tempString,'.',back=.true.) @@ -308,11 +304,9 @@ character(len=1024) function getGeometryFile(geometryParameter) cwd integer :: posExt, posSep logical :: error - character :: pathSep external :: quit getGeometryFile = geometryParameter - pathSep = getPathSep() posExt = scan(getGeometryFile,'.',back=.true.) posSep = scan(getGeometryFile,pathSep,back=.true.) @@ -344,11 +338,9 @@ character(len=1024) function getLoadCaseFile(loadCaseParameter) cwd integer :: posExt, posSep logical :: error - character :: pathSep external :: quit getLoadCaseFile = loadcaseParameter - pathSep = getPathSep() posExt = scan(getLoadCaseFile,'.',back=.true.) posSep = scan(getLoadCaseFile,pathSep,back=.true.) @@ -374,11 +366,8 @@ function rectifyPath(path) implicit none character(len=*) :: path character(len=len_trim(path)) :: rectifyPath - character :: pathSep integer :: i,j,k,l ! no pInt - pathSep = getPathSep() - !-------------------------------------------------------------------------------------------------- ! remove /./ from path l = len_trim(path) @@ -415,10 +404,8 @@ character(len=1024) function makeRelativePath(a,b) implicit none character (len=*) :: a,b - character :: pathSep integer :: i,posLastCommonSlash,remainingSlashes !no pInt - pathSep = getPathSep() posLastCommonSlash = 0 remainingSlashes = 0 @@ -434,35 +421,6 @@ character(len=1024) function makeRelativePath(a,b) 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 !-------------------------------------------------------------------------------------------------- diff --git a/code/system_routines.f90 b/code/system_routines.f90 index ab1aae03f..07e12a20b 100644 --- a/code/system_routines.f90 +++ b/code/system_routines.f90 @@ -9,7 +9,8 @@ module system_routines public :: & isDirectory, & - getCWD + getCWD, & + getHostName interface @@ -29,6 +30,14 @@ interface integer(C_INT),intent(out) :: stat 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 @@ -85,5 +94,34 @@ logical function getCWD(str) 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 diff --git a/examples/ConfigFiles/Phase_Nonlocal_Aluminum.config b/examples/ConfigFiles/Phase_Nonlocal_Aluminum.config index b805f09be..6406a47c0 100644 --- a/examples/ConfigFiles/Phase_Nonlocal_Aluminum.config +++ b/examples/ConfigFiles/Phase_Nonlocal_Aluminum.config @@ -8,32 +8,20 @@ plasticity nonlocal (output) rho_edge (output) rho_screw (output) rho_sgl -(output) rho_sgl_edge (output) rho_sgl_edge_pos (output) rho_sgl_edge_neg -(output) rho_sgl_screw (output) rho_sgl_screw_pos (output) rho_sgl_screw_neg -(output) rho_sgl_mobile -(output) rho_sgl_edge_mobile (output) rho_sgl_edge_pos_mobile (output) rho_sgl_edge_neg_mobile -(output) rho_sgl_screw_mobile (output) rho_sgl_screw_pos_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_neg_immobile -(output) rho_sgl_screw_immobile (output) rho_sgl_screw_pos_immobile (output) rho_sgl_screw_neg_immobile -(output) rho_dip (output) rho_dip_edge (output) rho_dip_screw -(output) excess_rho -(output) excess_rho_edge -(output) excess_rho_screw (output) rho_forest (output) delta (output) delta_sgl @@ -47,10 +35,8 @@ plasticity nonlocal (output) rho_dot_sgl (output) rho_dot_sgl_mobile (output) rho_dot_dip -(output) rho_dot_gen (output) rho_dot_gen_edge (output) rho_dot_gen_screw -(output) rho_dot_sgl2dip (output) rho_dot_sgl2dip_edge (output) rho_dot_sgl2dip_screw (output) rho_dot_ann_ath @@ -66,24 +52,6 @@ plasticity nonlocal (output) velocity_edge_neg (output) velocity_screw_pos (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_screw (output) accumulated_shear diff --git a/examples/ConfigFiles/Phase_Nonlocal_Nickel.config b/examples/ConfigFiles/Phase_Nonlocal_Nickel.config old mode 100755 new mode 100644 index 4d4c2d1df..3420b4246 --- a/examples/ConfigFiles/Phase_Nonlocal_Nickel.config +++ b/examples/ConfigFiles/Phase_Nonlocal_Nickel.config @@ -14,8 +14,6 @@ plasticity nonlocal (output) rho_dip_edge (output) rho_dip_screw (output) rho_forest -(output) excess_rho_edge -(output) excess_rho_screw (output) accumulatedshear (output) shearrate (output) resolvedstress @@ -30,24 +28,6 @@ plasticity nonlocal (output) rho_dot_edgejogs (output) rho_dot_flux_edge (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 Nslip 12 # number of slip systems per family diff --git a/examples/ConfigFiles/Phase_Phenopowerlaw_Aluminum.config b/examples/ConfigFiles/Phase_Phenopowerlaw_Aluminum.config index 8fa58557a..63ec8f3c8 100644 --- a/examples/ConfigFiles/Phase_Phenopowerlaw_Aluminum.config +++ b/examples/ConfigFiles/Phase_Phenopowerlaw_Aluminum.config @@ -7,11 +7,6 @@ plasticity phenopowerlaw (output) resolvedstress_slip (output) accumulated_shear_slip (output) totalshear -(output) resistance_twin -(output) shearrate_twin -(output) resolvedstress_twin -(output) accumulated_shear_twin -(output) totalvolfrac_twin lattice_structure fcc Nslip 12 # per family @@ -26,19 +21,6 @@ n_slip 20 tau0_slip 31e6 # per family tausat_slip 63e6 # per family 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_twinslip 0 -h0_twintwin 0 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 diff --git a/lib/damask/orientation.py b/lib/damask/orientation.py index f34f33a53..d5274571b 100644 --- a/lib/damask/orientation.py +++ b/lib/damask/orientation.py @@ -995,11 +995,17 @@ class Orientation: relationModel, direction, 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 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 # NW from H. Kitahara et al./Materials Characterization 54 (2005) 378-386 @@ -1226,14 +1232,14 @@ class Orientation: 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 /= 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 /= np.linalg.norm(otherPlane) otherNormal = [float(i) for i in normals[relationModel][variant,other]] # map(float, planes[...]) does not work in python 3 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 ?? diff --git a/lib/damask/test/test.py b/lib/damask/test/test.py index db58fd4a5..5a06939c6 100644 --- a/lib/damask/test/test.py +++ b/lib/damask/test/test.py @@ -23,6 +23,8 @@ class Test(): 'keep': False, 'accept': False, 'updateRequest': False, + 'show': False, + 'select': None, } for arg in defaults.keys(): setattr(self,arg,kwargs.get(arg) if kwargs.get(arg) else defaults[arg]) @@ -58,10 +60,18 @@ class Test(): action = "store_true", dest = "accept", 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, accept = self.accept, update = self.updateRequest, + show = self.show, + select = self.select, ) @@ -73,21 +83,26 @@ class Test(): self.prepareAll() for variant,name in enumerate(self.variants): - try: - if not self.options.keep: - self.prepare(variant) - self.run(variant) + if self.options.show: + logging.critical('{}: {}'.format(variant,name)) + elif self.options.select is not None and name != self.options.select: + 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.update(variant) != 0: logging.critical('update for "{}" failed.'.format(name)) - elif not (self.options.accept or self.compare(variant)): # no update, do comparison - return variant+1 # return culprit + if self.options.update: + 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 + return variant+1 # return culprit - except Exception as e : - logging.critical('exception during variant execution: {}'.format(e)) - return variant+1 # return culprit + except Exception as e : + logging.critical('exception during variant execution: {}'.format(e)) + return variant+1 # return culprit return 0 def feasible(self): diff --git a/processing/misc/DREAM3D_toTable.py b/processing/misc/DREAM3D_toTable.py new file mode 100755 index 000000000..4cd93cf5b --- /dev/null +++ b/processing/misc/DREAM3D_toTable.py @@ -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() diff --git a/processing/pre/geom_check.sh b/processing/pre/geom_check.sh index 1eb85913c..5a39d4fc7 100755 --- a/processing/pre/geom_check.sh +++ b/processing/pre/geom_check.sh @@ -11,7 +11,7 @@ do < $geom \ | \ vtk_addRectilinearGridData \ - --scalar microstructure \ + --data microstructure \ --inplace \ --vtk ${geom%.*}.vtk rm ${geom%.*}.vtk diff --git a/processing/pre/geom_grainGrowth.py b/processing/pre/geom_grainGrowth.py index ac7af9e85..e67273d38 100755 --- a/processing/pre/geom_grainGrowth.py +++ b/processing/pre/geom_grainGrowth.py @@ -89,9 +89,9 @@ for name in filenames: # 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[:,:,grid[2]/2::] = gauss[:,:,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[0]/2::,:,:] = gauss[round(grid[0]/2.)-1::-1,:,:] + 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[:,int(round(grid[1]/2.))-1::-1,:] + gauss[grid[0]/2::,:,:] = gauss[int(round(grid[0]/2.))-1::-1,:,:] 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 diff --git a/processing/pre/geom_mirror.py b/processing/pre/geom_mirror.py new file mode 100755 index 000000000..cc51749d4 --- /dev/null +++ b/processing/pre/geom_mirror.py @@ -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 = '', + 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 diff --git a/processing/pre/patchFromReconstructedBoundaries.py b/processing/pre/patchFromReconstructedBoundaries.py index ada92b09d..572c929fa 100755 --- a/processing/pre/patchFromReconstructedBoundaries.py +++ b/processing/pre/patchFromReconstructedBoundaries.py @@ -2,6 +2,7 @@ # -*- coding: UTF-8 no BOM -*- import sys,os,math,re +import numpy as np from optparse import OptionParser import damask @@ -63,9 +64,8 @@ def rcbOrientationParser(content,idcolumn): grains = [] myOrientation = [0.0,0.0,0.0] - for line in content: - m = re.match(r'\s*(#|$)',line) - if m: continue # skip comments and blank lines + for j,line in enumerate(content): + if re.match(r'^\s*(#|$)',line): continue # skip comments and blank lines for grain in range(2): myID = int(line.split()[idcolumn+grain]) # get grain id myOrientation = map(float,line.split())[3*grain:3+3*grain] # get orientation @@ -75,8 +75,8 @@ def rcbOrientationParser(content,idcolumn): try: grains[myID-1] = myOrientation # store Euler angles except IndexError: - message = 'You might not have chosen the correct column for the grain IDs! Please check the "--id" option.' - print '\033[1;31m'+message+'\033[0m\n' + damask.util.croak('You might not have chosen the correct column for the grain IDs! '+ + 'Please check the "--id" option.') raise except: raise @@ -91,13 +91,13 @@ def rcbParser(content,M,size,tolerance,idcolumn,segmentcolumn): x = [0.,0.] y = [0.,0.] for line in content: - m = re.match(r'\s*(#|$)',line) + m = re.match(r'^\s*(#|$)',line) if m: continue # skip comments and blank lines try: (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: - message = 'You might not have chosen the correct column for the segment end points! Please check the "--segment" option.' - print '\033[1;31m'+message+'\033[0m\n' + damask.util.croak('You might not have chosen the correct column for the segment end points! '+ + 'Please check the "--segment" option.') raise except: raise @@ -110,6 +110,9 @@ def rcbParser(content,M,size,tolerance,idcolumn,segmentcolumn): dX = boxX[1]-boxX[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 else: scalePatch = 1.0 @@ -122,8 +125,7 @@ def rcbParser(content,M,size,tolerance,idcolumn,segmentcolumn): grainNeighbors = [] for line in content: - m = re.match(r'\s*(#|$)',line) - if m: continue # skip comments and blank lines + if re.match(r'^\s*(#|$)',line): 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]) = (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 @@ -133,8 +135,8 @@ def rcbParser(content,M,size,tolerance,idcolumn,segmentcolumn): y[0] -= boxY[0] y[1] -= boxY[0] 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 - match = False # check whether point is already known (within a small range) + for i in range(2): # store segment to both points + match = False # check whether point is already known (within a small range) for posX in connectivityXY.keys(): if (abs(float(posX)-x[i]) 0: + damask.util.croak(' culling {} duplicate segments...'.format(len(dupSegments))) + for rm in dupSegments: + segments[rm] = None + crappyData = False for pointId,point in enumerate(points): if len(point['segments']) < 2: # point marks a dead end! - print "Dead end at segment %i (%f,%f)"\ - %(1+point['segments'][0],boxX[0]+point['coords'][0]/scalePatch,boxY[0]+point['coords'][1]/scalePatch,) + damask.util.croak('dead end at segment {} for point {} ({},{}).' + .format(point['segments'][0], + pointId, + boxX[0]+point['coords'][0]/scalePatch,boxY[0]+point['coords'][1]/scalePatch,)) crappyData = True grains = {'draw': [], 'legs': []} @@ -249,58 +269,69 @@ def rcbParser(content,M,size,tolerance,idcolumn,segmentcolumn): innerAngleSum = 0.0 myWalk = point['segments'].pop() grainLegs = [myWalk] - if segments[myWalk][0] == myStart: - myEnd = segments[myWalk][1] - else: - myEnd = segments[myWalk][0] - + myEnd = segments[myWalk][1 if segments[myWalk][0] == myStart else 0] while (myEnd != pointId): - myV = [points[myEnd]['coords'][0]-points[myStart]['coords'][0],\ - points[myEnd]['coords'][1]-points[myStart]['coords'][1]] + myV = [points[myEnd]['coords'][0]-points[myStart]['coords'][0], + points[myEnd]['coords'][1]-points[myStart]['coords'][1]] 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} for peek in points[myEnd]['segments']: # trying in turn all segments emanating from current end if peek == myWalk: - continue - peekEnd = segments[peek][1] if segments[peek][0] == myEnd else segments[peek][0] - peekV = [points[myEnd]['coords'][0]-points[peekEnd]['coords'][0], - points[myEnd]['coords'][1]-points[peekEnd]['coords'][1]] + continue # do not go back same path + peekEnd = segments[peek][1 if segments[peek][0] == myEnd else 0] + peekV = [points[peekEnd]['coords'][0]-points[myEnd]['coords'][0], + points[peekEnd]['coords'][1]-points[myEnd]['coords'][1]] peekLen = math.sqrt(peekV[0]**2+peekV[1]**2) - crossproduct = (myV[0]*peekV[1]-myV[1]*peekV[0])/myLen/peekLen - dotproduct = (myV[0]*peekV[0]+myV[1]*peekV[1])/myLen/peekLen - if crossproduct*(dotproduct+1.0) >= best['product']: - best['product'] = crossproduct*(dotproduct+1.0) + if peekLen == 0.0: damask.util.croak('peeklen is zero: peek point {}'.format(peek)) + crossproduct = (myV[0]*peekV[1] - myV[1]*peekV[0])/myLen/peekLen + dotproduct = (myV[0]*peekV[0] + myV[1]*peekV[1])/myLen/peekLen + innerAngle = math.copysign(1.0,crossproduct)*(dotproduct-1.0) + if innerAngle >= best['product']: # takes sharpest left turn + best['product'] = innerAngle best['peek'] = peek best['point'] = peekEnd + innerAngleSum += best['product'] myWalk = best['peek'] myStart = myEnd myEnd = best['point'] + if myWalk in points[myStart]['segments']: points[myStart]['segments'].remove(myWalk) 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']) grainLegs.append(myWalk) + if innerAngleSum > 0.0: grains['draw'].append(grainDraw) grains['legs'].append(grainLegs) else: grains['box'] = grainLegs - # build overall data structure - rcData = {'dimension':[dX,dY], 'point': [],'segment': [], 'grain': [], 'grainMapping': []} - print " dimension %g x %g"%(dX,dY) + rcData = {'dimension':[dX,dY], + 'bounds': [[boxX[0],boxY[0]],[boxX[1],boxY[1]]], + 'scale': scalePatch, + 'point': [], + 'segment': [], + 'neighbors': [], + 'grain': [], + 'grainMapping': [], + } for point in points: 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) - 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 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 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 - if candidateGrains[0][0] not in rcData['grainMapping']: # most frequent one not yet seen? - rcData['grainMapping'].append(candidateGrains[0][0]) # must be me then - else: - rcData['grainMapping'].append(candidateGrains[1][0]) # special case of bi-crystal situation... + # most frequent one not yet seen? + rcData['grainMapping'].append(candidateGrains[0 if candidateGrains[0][0] not in rcData['grainMapping'] else 1][0]) # must be me then + # 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 [] @@ -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)) for id,vertex in enumerate(rcData['segment']): + if vertex: start = rcData['point'][vertex[0]] 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), 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)} frameindex=len(rcData['grain'])+1 # calculate frame index as largest grain index plus one - dx = xsize/(xres+1) # calculate step sizes - dy = ysize/(yres+1) + dx = xsize/(xres) # calculate step sizes + dy = ysize/(yres) grainpoints = [] 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 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 - bestGuess = 0 # assume grain 0 as best guess for i in range(int(xres*yres)): # walk through all points in xy plane 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( ytest < 0 or ytest > maxY): # part of edges 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'] parser.add_option('-o', '--output', action='extend', dest='output', metavar = '', - 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', dest='port', help='Mentat connection port [%default]') parser.add_option('-2', '--twodimensional', action='store_true', @@ -847,16 +877,15 @@ parser.set_defaults(output = [], (options, args) = parser.parse_args() -print '\033[1m'+scriptName+'\033[0m\n' if not len(args): - parser.error('no boundary file specified') + parser.error('no boundary file specified.') try: boundaryFile = open(args[0]) boundarySegments = boundaryFile.readlines() boundaryFile.close() except: - print 'unable to read boundary file "%s"'%args[0] + damask.util.croak('unable to read boundary file "{}".'.format(args[0])) raise 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 myName = os.path.splitext(args[0])[0] -print "%s\n"%myName +damask.util.report(scriptName,myName) orientationData = rcbOrientationParser(boundarySegments,options.idcolumn) 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 ----- if 'image' in options.output and options.imgsize > 0: if ImageCapability: image(myName,options.imgsize,options.xmargin,options.ymargin,rcData) else: - print '...no image drawing possible (PIL missing)...' + damask.util.croak('...no image drawing possible (PIL missing)...') # ----- 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.write('3\theader\n') # write header info - geomFile.write('resolution a %i b %i c %i\n'%(fftdata['resolution'])) # resolution - geomFile.write('dimension x %f y %f z %f\n'%(fftdata['dimension'])) # size + geomFile.write('grid a %i b %i c %i\n'%(fftdata['resolution'])) # grid resolution + geomFile.write('size x %f y %f z %f\n'%(fftdata['dimension'])) # size geomFile.write('homogenization 1\n') # homogenization for z in xrange(fftdata['resolution'][2]): # z repetions 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 geomFile.close() # close geom file - print('assigned %i out of %i (2D) Fourier points.'\ - %(len(fftdata['fftpoints']), int(fftdata['resolution'][0])*int(fftdata['resolution'][1]))) + damask.util.croak('assigned {} out of {} (2D) Fourier points...' + .format(len(fftdata['fftpoints']), + int(fftdata['resolution'][0])*int(fftdata['resolution'][1]))) # ----- write Mentat procedure ----- @@ -936,7 +994,7 @@ if 'mentat' in options.output: if 'procedure' in options.output: output(outputLocals['log'],outputLocals,'Stdout') else: - print '...no interaction with Mentat possible...' + damask.util.croak('...no interaction with Mentat possible...') # ----- write config data to file ----- diff --git a/processing/pre/seeds_check.sh b/processing/pre/seeds_check.sh index f24cd1c64..9bc054406 100755 --- a/processing/pre/seeds_check.sh +++ b/processing/pre/seeds_check.sh @@ -5,7 +5,7 @@ do vtk_pointcloud $seeds vtk_addPointcloudData $seeds \ - --scalar microstructure,weight \ + --data microstructure,weight \ --inplace \ --vtk ${seeds%.*}.vtp \