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

Merge remote-tracking branch 'origin/development' into restructure-numerics

This commit is contained in:
Sharan Roongta 2023-07-16 15:41:37 +02:00
parent 915d03f2bf 6770bedb7d
commit 054f548f4c
23 changed files with 300 additions and 303 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

27
.gitlab-ci.yml
View File

@ -196,7 +196,7 @@ Marc:
################################################################################################### ###################################################################################################
grid_runtime: grid_performance:
stage: statistics stage: statistics
before_script: before_script:
- ${LOCAL_HOME}/bin/queue ${CI_JOB_ID} --blocking - ${LOCAL_HOME}/bin/queue ${CI_JOB_ID} --blocking
@ -210,22 +210,27 @@ grid_runtime:
- export PATH=${PWD}/bin:${PATH} - export PATH=${PWD}/bin:${PATH}
- cd $(mktemp -d) - cd $(mktemp -d)
- git clone -q git@git.damask.mpie.de:damask/performance.git . - git clone -q git@git.damask.mpie.de:damask/performance.git .
- ./measure_performance.py --input_dir ${CI_PROJECT_DIR}/examples/grid --tag ${CI_COMMIT_SHA}
- > - >
${CI_PROJECT_DIR}/PRIVATE/testing/runtime.py if [ ${CI_COMMIT_BRANCH} == development ]; then
--input_dir ${CI_PROJECT_DIR}/examples/grid git add performance.txt
--output_dir ./ git commit -m ${CI_PIPELINE_ID}_${CI_COMMIT_SHA}
--tag ${CI_COMMIT_SHA} git push
- if [ ${CI_COMMIT_BRANCH} == development ]; then git commit -am ${CI_PIPELINE_ID}_${CI_COMMIT_SHA}; git push; fi fi
commit_history: update_plots:
stage: statistics stage: statistics
script: script:
- cd $(mktemp -d) - cd $(mktemp -d)
- ${CI_PROJECT_DIR}/PRIVATE/testing/plot_commithistory.py --color green -n 5 -N 100 - git clone -q git@git.damask.mpie.de:damask/performance.git .
- ${CI_PROJECT_DIR}/PRIVATE/testing/plot_commithistory.py --color green -n 5 -N 1000 - ./plot_commithistory.py --color green -n 5 -N 100
- ${CI_PROJECT_DIR}/PRIVATE/testing/plot_commithistory.py --color green -n 5 -N 10000 - ./plot_commithistory.py --color green -n 5 -N 1000
- ./plot_commithistory.py --color green -n 5 -N 10000
- scp -r ./commits_*.html damask3.mpie.de:~/ - scp -r ./commits_*.html damask3.mpie.de:~/
- ssh damask3.mpie.de "./update_statistics.sh" - ssh damask3.mpie.de "./update_statistics_commits.sh"
- ./plot_performance.py --template=xgridoff
- scp -r ./runtime.html ./memory.html damask3.mpie.de:~/
- ssh damask3.mpie.de "./update_statistics_performance.sh"
only: only:
- development - development

2
VERSION
View File

@ -1 +1 @@
3.0.0-alpha7-614-gad6220c26 3.0.0-alpha7-630-ga63fe8d02

57
src/CLI.f90
View File

@ -74,7 +74,7 @@ subroutine CLI_init()
print'(a)', achar(27)//'[31m' print'(a)', achar(27)//'[31m'
print'(1x,a,/)', 'debug version - debug version - debug version - debug version - debug version' print'(1x,a,/)', 'debug version - debug version - debug version - debug version - debug version'
#else #else
print'(a)', achar(27)//'[94m' print'(a)', achar(27)//'[1;94m'
#endif #endif
print'(1x,a)', ' _/_/_/ _/_/ _/ _/ _/_/ _/_/_/ _/ _/ _/_/_/' print'(1x,a)', ' _/_/_/ _/_/ _/ _/ _/_/ _/_/_/ _/ _/ _/_/_/'
print'(1x,a)', ' _/ _/ _/ _/ _/_/ _/_/ _/ _/ _/ _/ _/ _/' print'(1x,a)', ' _/ _/ _/ _/ _/_/ _/_/ _/ _/ _/ _/ _/ _/'
@ -158,48 +158,34 @@ subroutine CLI_init()
print'(1x,a,/)',' Prints this message and exits' print'(1x,a,/)',' Prints this message and exits'
call quit(0) ! normal Termination call quit(0) ! normal Termination
case ('-g', '--geom', '--geometry') case ('-g', '--geom', '--geometry')
if (.not. hasArg) print'(/,1x,a)', 'ERROR: missing argument for --geom' if (.not. hasArg) call IO_error(610,ext_msg='--geom')
geomArg = getArg(i+1) geomArg = getArg(i+1)
case ('-l', '--load', '--loadcase') case ('-l', '--load', '--loadcase')
if (.not. hasArg) print'(/,1x,a)', 'ERROR: missing argument for --load' if (.not. hasArg) call IO_error(610,ext_msg='--load')
loadArg = getArg(i+1) loadArg = getArg(i+1)
case ('-m', '--material', '--materialconfig') case ('-m', '--material', '--materialconfig')
if (.not. hasArg) print'(/,1x,a)', 'ERROR: missing argument for --material' if (.not. hasArg) call IO_error(610,ext_msg='--material')
materialArg = getArg(i+1) materialArg = getArg(i+1)
case ('-n', '--numerics', '--numericsconfig') case ('-n', '--numerics', '--numericsconfig')
if (.not. hasArg) print'(/,1x,a)', 'ERROR: missing argument for --numerics' if (.not. hasArg) call IO_error(610,ext_msg='--numerics')
numericsArg = getArg(i+1) numericsArg = getArg(i+1)
case ('-j', '--job', '--jobname') case ('-j', '--job', '--jobname')
if (.not. hasArg) print'(/,1x,a)', 'ERROR: missing argument for --jobname' if (.not. hasArg) call IO_error(610,ext_msg='--jobname')
solverJobname = getArg(i+1) solverJobname = getArg(i+1)
case ('-w', '--wd', '--workingdir', '--workingdirectory') case ('-w', '--wd', '--workingdir', '--workingdirectory')
if (.not. hasArg) print'(/,1x,a)', 'ERROR: missing argument for --workingdirectory' if (.not. hasArg) call IO_error(610,ext_msg='--workingdirectory')
workingDirArg = getArg(i+1) workingDirArg = getArg(i+1)
case ('-r', '--rs', '--restart') case ('-r', '--rs', '--restart')
if (.not. hasArg) print'(/,1x,a)', 'ERROR: missing argument for --restart' if (.not. hasArg) call IO_error(610,ext_msg='--jobname')
arg = getArg(i+1) arg = getArg(i+1)
read(arg,*,iostat=stat) CLI_restartInc read(arg,*,iostat=stat) CLI_restartInc
if (CLI_restartInc < 0 .or. stat /= 0) then if (CLI_restartInc < 0 .or. stat /= 0) call IO_error(611,ext_msg=arg)
print'(/,1x,a)', 'ERROR: could not parse restart increment: '//trim(arg)
call quit(1)
end if
end select end select
end do end do
if (.not. allocated(loadArg)) then if (.not. allocated(geomArg)) call IO_error(612,ext_msg='--geom')
print'(/,1x,a)', 'Error: no load case specified (-h for help)' if (.not. allocated(loadArg)) call IO_error(612,ext_msg='--load')
call quit(1) if (.not. allocated(materialArg)) call IO_error(612,ext_msg='--material')
end if
if (.not. allocated(geomArg)) then
print'(/,1x,a)', 'Error: no geometry specified (-h for help)'
call quit(1)
end if
if (.not. allocated(materialArg)) then
print'(/,1x,a)', 'Error: no material configuration specified (-h for help)'
call quit(1)
end if
call setWorkingDirectory(trim(workingDirArg)) call setWorkingDirectory(trim(workingDirArg))
CLI_geomFile = getPathRelCWD(geomArg,'geometry') CLI_geomFile = getPathRelCWD(geomArg,'geometry')
@ -211,8 +197,7 @@ subroutine CLI_init()
if (.not. allocated(solverJobname)) then if (.not. allocated(solverJobname)) then
solverJobname = jobname(CLI_geomFile,CLI_loadFile,CLI_materialFile,CLI_numericsFile) solverJobname = jobname(CLI_geomFile,CLI_loadFile,CLI_materialFile,CLI_numericsFile)
elseif (scan(solverJobname,'/') > 0) then elseif (scan(solverJobname,'/') > 0) then
print'(/,1x,a)', 'ERROR: JOBNAME must not contain any slashes' call IO_error(630)
call quit(1)
endif endif
commandLine = getArg(-1) commandLine = getArg(-1)
@ -272,9 +257,6 @@ subroutine setWorkingDirectory(workingDirectoryArg)
character(len=*), intent(in) :: workingDirectoryArg !< working directory argument character(len=*), intent(in) :: workingDirectoryArg !< working directory argument
character(len=:), allocatable :: workingDirectory character(len=:), allocatable :: workingDirectory
logical :: error
external :: quit
absolutePath: if (workingDirectoryArg(1:1) == '/') then absolutePath: if (workingDirectoryArg(1:1) == '/') then
workingDirectory = workingDirectoryArg workingDirectory = workingDirectoryArg
@ -284,11 +266,7 @@ subroutine setWorkingDirectory(workingDirectoryArg)
end if absolutePath end if absolutePath
workingDirectory = trim(normpath(workingDirectory)) workingDirectory = trim(normpath(workingDirectory))
error = setCWD(trim(workingDirectory)) if (setCWD(trim(workingDirectory))) call IO_error(640,ext_msg=workingDirectory)
if (error) then
print'(1x,a)', 'ERROR: invalid working directory: '//trim(workingDirectory)
call quit(1)
end if
end subroutine setWorkingDirectory end subroutine setWorkingDirectory
@ -344,7 +322,6 @@ function getPathRelCWD(path,fileType)
character(len=*), intent(in) :: fileType character(len=*), intent(in) :: fileType
logical :: file_exists logical :: file_exists
external :: quit
getPathRelCWD = trim(path) getPathRelCWD = trim(path)
@ -352,10 +329,7 @@ function getPathRelCWD(path,fileType)
getPathRelCWD = trim(relpath(getPathRelCWD,getCWD())) getPathRelCWD = trim(relpath(getPathRelCWD,getCWD()))
inquire(file=getPathRelCWD, exist=file_exists) inquire(file=getPathRelCWD, exist=file_exists)
if (.not. file_exists) then if (.not. file_exists) call IO_error(100,ext_msg=fileType//' "'//trim(getPathRelCWD)//'"')
print'(/,1x,a)', 'ERROR: '//fileType//' file does not exist: '//trim(getPathRelCWD)
call quit(1)
end if
end function getPathRelCWD end function getPathRelCWD
@ -438,4 +412,5 @@ function relpath(path,start)
end function relpath end function relpath
end module CLI end module CLI

2
src/CMakeLists.txt
View File

@ -1,7 +1,7 @@
# special flags for some files # special flags for some files
if(CMAKE_Fortran_COMPILER_ID STREQUAL "GNU") if(CMAKE_Fortran_COMPILER_ID STREQUAL "GNU")
# long lines for interaction matrix # long lines for interaction matrix
set_source_files_properties("lattice.f90" PROPERTIES COMPILE_FLAGS "-ffree-line-length-240") set_source_files_properties("crystal.f90" PROPERTIES COMPILE_FLAGS "-ffree-line-length-240")
set_source_files_properties("parallelization.f90" PROPERTIES COMPILE_FLAGS "-ffree-line-length-none") set_source_files_properties("parallelization.f90" PROPERTIES COMPILE_FLAGS "-ffree-line-length-none")
endif() endif()

47
src/IO.f90
View File

@ -500,7 +500,7 @@ subroutine IO_error(error_ID,ext_msg,label1,ID1,label2,ID2)
case (131) case (131)
msg = 'hex lattice structure with invalid c/a ratio' msg = 'hex lattice structure with invalid c/a ratio'
case (132) case (132)
msg = 'trans_lattice_structure not possible' msg = 'invalid parameters for transformation'
case (134) case (134)
msg = 'negative lattice parameter' msg = 'negative lattice parameter'
case (135) case (135)
@ -579,6 +579,18 @@ subroutine IO_error(error_ID,ext_msg,label1,ID1,label2,ID2)
! user errors ! user errors
case (603) case (603)
msg = 'invalid data for table' msg = 'invalid data for table'
case (610)
msg = 'missing argument for option'
case (611)
msg = 'could not parse restart increment'
case (612)
msg = 'missing option'
case (630)
msg = 'JOBNAME must not contain any slashes'
case (640)
msg = 'invalid working directory'
!------------------------------------------------------------------------------------------------ !------------------------------------------------------------------------------------------------
! errors related to YAML data ! errors related to YAML data
@ -728,38 +740,43 @@ subroutine panel(paneltype,ID,msg,ext_msg,label1,ID1,label2,ID2)
character(len=pSTRLEN) :: formatString character(len=pSTRLEN) :: formatString
integer, parameter :: panelwidth = 69 integer, parameter :: panelwidth = 69
character(len=:), allocatable :: msg_,ID_,msg1,msg2
character(len=*), parameter :: DIVIDER = repeat('─',panelwidth) character(len=*), parameter :: DIVIDER = repeat('─',panelwidth)
if (.not. present(label1) .and. present(ID1)) error stop 'missing label for value 1' if (.not. present(label1) .and. present(ID1)) error stop 'missing label for value 1'
if (.not. present(label2) .and. present(ID2)) error stop 'missing label for value 2' if (.not. present(label2) .and. present(ID2)) error stop 'missing label for value 2'
if ( present(label1) .and. .not. present(ID1)) error stop 'missing value for label 1'
if ( present(label2) .and. .not. present(ID2)) error stop 'missing value for label 2'
ID_ = IO_intAsStr(ID)
if (present(label1)) msg1 = label1
if (present(label2)) msg2 = label2
if (present(ID1)) msg1 = msg1//' '//IO_intAsStr(ID1)
if (present(ID2)) msg2 = msg2//' '//IO_intAsStr(ID2)
if (paneltype == 'error') msg_ = achar(27)//'[31m'//trim(msg)//achar(27)//'[0m'
if (paneltype == 'warning') msg_ = achar(27)//'[33m'//trim(msg)//achar(27)//'[0m'
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(IO_STDERR,'(/,a)') ' ┌'//DIVIDER//'┐' write(IO_STDERR,'(/,a)') ' ┌'//DIVIDER//'┐'
write(formatString,'(a,i2,a)') '(a,24x,a,',max(1,panelwidth-24-len_trim(paneltype)),'x,a)' write(formatString,'(a,i2,a)') '(a,24x,a,1x,i0,',max(1,panelwidth-24-len_trim(paneltype)-1-len_trim(ID_)),'x,a)'
write(IO_STDERR,formatString) ' │',trim(paneltype), '│' write(IO_STDERR,formatString) ' │',trim(paneltype),ID, '│'
write(formatString,'(a,i2,a)') '(a,24x,i3,',max(1,panelwidth-24-3),'x,a)'
write(IO_STDERR,formatString) ' │',ID, '│'
write(IO_STDERR,'(a)') ' ├'//DIVIDER//'┤' write(IO_STDERR,'(a)') ' ├'//DIVIDER//'┤'
write(formatString,'(a,i3.3,a,i3.3,a)') '(1x,a4,a',max(1,len_trim(msg)),',',& write(formatString,'(a,i3.3,a,i3.3,a)') '(1x,a4,a',max(1,len_trim(msg_)),',',&
max(1,panelwidth+3-len_trim(msg)-4),'x,a)' max(1,panelwidth+3-len_trim(msg)-4),'x,a)'
write(IO_STDERR,formatString) '│ ',trim(msg), '│' write(IO_STDERR,formatString) '│ ',trim(msg_), '│'
if (present(ext_msg)) then if (present(ext_msg)) then
write(formatString,'(a,i3.3,a,i3.3,a)') '(1x,a4,a',max(1,len_trim(ext_msg)),',',& write(formatString,'(a,i3.3,a,i3.3,a)') '(1x,a4,a',max(1,len_trim(ext_msg)),',',&
max(1,panelwidth+3-len_trim(ext_msg)-4),'x,a)' max(1,panelwidth+3-len_trim(ext_msg)-4),'x,a)'
write(IO_STDERR,formatString) '│ ',trim(ext_msg), '│' write(IO_STDERR,formatString) '│ ',trim(ext_msg), '│'
end if end if
if (present(label1)) then if (present(label1)) then
write(formatString,'(a,i3.3,a,i3.3,a)') '(1x,a7,a',max(1,len_trim(label1)),',i9,',& write(formatString,'(a,i3.3,a,i3.3,a)') '(1x,a7,a',max(1,len_trim(msg1)),',',&
max(1,panelwidth+3-len_trim(label1)-9-7),'x,a)' max(1,panelwidth+3-len_trim(msg1)-7),'x,a)'
write(IO_STDERR,formatString) '│ at ',trim(label1),ID1, '│' write(IO_STDERR,formatString) '│ at ',trim(msg1), '│'
end if end if
if (present(label2)) then if (present(label2)) then
write(formatString,'(a,i3.3,a,i3.3,a)') '(1x,a7,a',max(1,len_trim(label2)),',i9,',& write(formatString,'(a,i3.3,a,i3.3,a)') '(1x,a7,a',max(1,len_trim(msg2)),',',&
max(1,panelwidth+3-len_trim(label2)-9-7),'x,a)' max(1,panelwidth+3-len_trim(msg2)-7),'x,a)'
write(IO_STDERR,formatString) '│ at ',trim(label2),ID2, '│' write(IO_STDERR,formatString) '│ at ',trim(msg2), '│'
end if end if
write(formatString,'(a,i2.2,a)') '(a,',max(1,panelwidth),'x,a)' write(formatString,'(a,i2.2,a)') '(a,',max(1,panelwidth),'x,a)'
write(IO_STDERR,formatString) ' │', '│' write(IO_STDERR,formatString) ' │', '│'

2
src/Marc/DAMASK_Marc.f90
View File

@ -155,7 +155,7 @@ end module DAMASK_interface
#include "../rotations.f90" #include "../rotations.f90"
#include "../polynomials.f90" #include "../polynomials.f90"
#include "../tables.f90" #include "../tables.f90"
#include "../lattice.f90" #include "../crystal.f90"
#include "element.f90" #include "element.f90"
#include "../geometry_plastic_nonlocal.f90" #include "../geometry_plastic_nonlocal.f90"
#include "../discretization.f90" #include "../discretization.f90"

4
src/Marc/materialpoint_Marc.f90
View File

@ -16,7 +16,7 @@ module materialpoint_Marc
use rotations use rotations
use polynomials use polynomials
use tables use tables
use lattice use crystal
use material use material
use phase use phase
use homogenization use homogenization
@ -75,7 +75,7 @@ subroutine materialpoint_initAll()
call rotations_init() call rotations_init()
call polynomials_init() call polynomials_init()
call tables_init() call tables_init()
call lattice_init() call crystal_init()
call discretization_Marc_init() call discretization_Marc_init()
call material_init(.false.) call material_init(.false.)
call phase_init() call phase_init()

316
src/lattice.f90 → src/crystal.f90
View File

@ -3,10 +3,10 @@
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH !> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief contains lattice definitions including Schmid matrices for slip, twin, trans, !> @brief Contains crystal definitions including Schmid matrices for slip, twin, trans,
! and cleavage as well as interaction among the various systems ! and cleavage as well as interaction among the various systems.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module lattice module crystal
use prec use prec
use misc use misc
use IO use IO
@ -80,7 +80,7 @@ module lattice
],pREAL),shape(CF_SYSTEMTWIN)) !< cF twin systems ],pREAL),shape(CF_SYSTEMTWIN)) !< cF twin systems
integer, dimension(2,CF_NTWIN), parameter, public :: & integer, dimension(2,CF_NTWIN), parameter, public :: &
lattice_CF_TWINNUCLEATIONSLIPPAIR = reshape( [& crystal_CF_TWINNUCLEATIONSLIPPAIR = reshape( [&
2,3, & 2,3, &
1,3, & 1,3, &
1,2, & 1,2, &
@ -93,7 +93,7 @@ module lattice
11,12, & 11,12, &
10,12, & 10,12, &
10,11 & 10,11 &
],shape(lattice_CF_TWINNUCLEATIONSLIPPAIR)) ],shape(crystal_CF_TWINNUCLEATIONSLIPPAIR))
real(pREAL), dimension(3+3,CF_NCLEAVAGE), parameter :: & real(pREAL), dimension(3+3,CF_NCLEAVAGE), parameter :: &
CF_SYSTEMCLEAVAGE = reshape(real([& CF_SYSTEMCLEAVAGE = reshape(real([&
@ -367,60 +367,60 @@ module lattice
],pREAL),shape(TI_SYSTEMSLIP)) !< tI slip systems for c/a = 0.5456 (Sn), sorted by Bieler 2009 (https://doi.org/10.1007/s11664-009-0909-x) ],pREAL),shape(TI_SYSTEMSLIP)) !< tI slip systems for c/a = 0.5456 (Sn), sorted by Bieler 2009 (https://doi.org/10.1007/s11664-009-0909-x)
interface lattice_forestProjection_edge interface crystal_forestProjection_edge
module procedure slipProjection_transverse module procedure slipProjection_transverse
end interface lattice_forestProjection_edge end interface crystal_forestProjection_edge
interface lattice_forestProjection_screw interface crystal_forestProjection_screw
module procedure slipProjection_direction module procedure slipProjection_direction
end interface lattice_forestProjection_screw end interface crystal_forestProjection_screw
public :: & public :: &
lattice_init, & crystal_init, &
lattice_isotropic_nu, & crystal_isotropic_nu, &
lattice_isotropic_mu, & crystal_isotropic_mu, &
lattice_symmetrize_33, & crystal_symmetrize_33, &
lattice_symmetrize_C66, & crystal_symmetrize_C66, &
lattice_SchmidMatrix_slip, & crystal_SchmidMatrix_slip, &
lattice_SchmidMatrix_twin, & crystal_SchmidMatrix_twin, &
lattice_SchmidMatrix_trans, & crystal_SchmidMatrix_trans, &
lattice_SchmidMatrix_cleavage, & crystal_SchmidMatrix_cleavage, &
lattice_nonSchmidMatrix, & crystal_nonSchmidMatrix, &
lattice_interaction_SlipBySlip, & crystal_interaction_SlipBySlip, &
lattice_interaction_TwinByTwin, & crystal_interaction_TwinByTwin, &
lattice_interaction_TransByTrans, & crystal_interaction_TransByTrans, &
lattice_interaction_SlipByTwin, & crystal_interaction_SlipByTwin, &
lattice_interaction_SlipByTrans, & crystal_interaction_SlipByTrans, &
lattice_interaction_TwinBySlip, & crystal_interaction_TwinBySlip, &
lattice_characteristicShear_Twin, & crystal_characteristicShear_Twin, &
lattice_C66_twin, & crystal_C66_twin, &
lattice_C66_trans, & crystal_C66_trans, &
lattice_forestProjection_edge, & crystal_forestProjection_edge, &
lattice_forestProjection_screw, & crystal_forestProjection_screw, &
lattice_slip_normal, & crystal_slip_normal, &
lattice_slip_direction, & crystal_slip_direction, &
lattice_slip_transverse, & crystal_slip_transverse, &
lattice_labels_slip, & crystal_labels_slip, &
lattice_labels_twin crystal_labels_twin
contains contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Run self test. !> @brief Run self test.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine lattice_init() subroutine crystal_init()
print'(/,1x,a)', '<<<+- lattice init -+>>>'; flush(IO_STDOUT) print'(/,1x,a)', '<<<+- crystal init -+>>>'; flush(IO_STDOUT)
call selfTest() call selfTest()
end subroutine lattice_init end subroutine crystal_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Characteristic shear for twinning !> @brief Characteristic shear for twinning
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_characteristicShear_Twin(Ntwin,lattice,CoverA) result(characteristicShear) function crystal_characteristicShear_Twin(Ntwin,lattice,CoverA) result(characteristicShear)
integer, dimension(:), intent(in) :: Ntwin !< number of active twin systems per family integer, dimension(:), intent(in) :: Ntwin !< number of active twin systems per family
character(len=*), intent(in) :: lattice !< Bravais lattice (Pearson symbol) character(len=*), intent(in) :: lattice !< Bravais lattice (Pearson symbol)
@ -470,7 +470,7 @@ function lattice_characteristicShear_Twin(Ntwin,lattice,CoverA) result(character
characteristicShear(a) = 0.5_pREAL*sqrt(2.0_pREAL) characteristicShear(a) = 0.5_pREAL*sqrt(2.0_pREAL)
case('hP') case('hP')
if (cOverA < 1.0_pREAL .or. cOverA > 2.0_pREAL) & if (cOverA < 1.0_pREAL .or. cOverA > 2.0_pREAL) &
call IO_error(131,ext_msg='lattice_characteristicShear_Twin') call IO_error(131,ext_msg='crystal_characteristicShear_Twin')
p = sum(HP_NTWINSYSTEM(1:f-1))+s p = sum(HP_NTWINSYSTEM(1:f-1))+s
select case(HP_SHEARTWIN(p)) ! from Christian & Mahajan 1995 p.29 select case(HP_SHEARTWIN(p)) ! from Christian & Mahajan 1995 p.29
case (1) ! <-10.1>{10.2} case (1) ! <-10.1>{10.2}
@ -483,24 +483,24 @@ function lattice_characteristicShear_Twin(Ntwin,lattice,CoverA) result(character
characteristicShear(a) = 2.0_pREAL*(cOverA**2-2.0_pREAL)/3.0_pREAL/cOverA characteristicShear(a) = 2.0_pREAL*(cOverA**2-2.0_pREAL)/3.0_pREAL/cOverA
end select end select
case default case default
call IO_error(137,ext_msg='lattice_characteristicShear_Twin: '//trim(lattice)) call IO_error(137,ext_msg='crystal_characteristicShear_Twin: '//trim(lattice))
end select end select
end do mySystems end do mySystems
end do myFamilies end do myFamilies
end function lattice_characteristicShear_Twin end function crystal_characteristicShear_Twin
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Rotated elasticity matrices for twinning in 6x6-matrix notation !> @brief Rotated elasticity matrices for twinning in 6x6-matrix notation
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_C66_twin(Ntwin,C66,lattice,CoverA) function crystal_C66_twin(Ntwin,C66,lattice,CoverA)
integer, dimension(:), intent(in) :: Ntwin !< number of active twin systems per family integer, dimension(:), intent(in) :: Ntwin !< number of active twin systems per family
character(len=*), intent(in) :: lattice !< Bravais lattice (Pearson symbol) character(len=*), intent(in) :: lattice !< Bravais lattice (Pearson symbol)
real(pREAL), dimension(6,6), intent(in) :: C66 !< unrotated parent stiffness matrix real(pREAL), dimension(6,6), intent(in) :: C66 !< unrotated parent stiffness matrix
real(pREAL), intent(in) :: cOverA !< c/a ratio real(pREAL), intent(in) :: cOverA !< c/a ratio
real(pREAL), dimension(6,6,sum(Ntwin)) :: lattice_C66_twin real(pREAL), dimension(6,6,sum(Ntwin)) :: crystal_C66_twin
real(pREAL), dimension(3,3,sum(Ntwin)):: coordinateSystem real(pREAL), dimension(3,3,sum(Ntwin)):: coordinateSystem
type(tRotation) :: R type(tRotation) :: R
@ -518,28 +518,28 @@ function lattice_C66_twin(Ntwin,C66,lattice,CoverA)
coordinateSystem = buildCoordinateSystem(Ntwin,HP_NSLIPSYSTEM,HP_SYSTEMTWIN,& coordinateSystem = buildCoordinateSystem(Ntwin,HP_NSLIPSYSTEM,HP_SYSTEMTWIN,&
lattice,cOverA) lattice,cOverA)
case default case default
call IO_error(137,ext_msg='lattice_C66_twin: '//trim(lattice)) call IO_error(137,ext_msg='crystal_C66_twin: '//trim(lattice))
end select end select
do i = 1, sum(Ntwin) do i = 1, sum(Ntwin)
call R%fromAxisAngle([coordinateSystem(1:3,2,i),PI],P=1) ! ToDo: Why always 180 deg? call R%fromAxisAngle([coordinateSystem(1:3,2,i),PI],P=1) ! ToDo: Why always 180 deg?
lattice_C66_twin(1:6,1:6,i) = R%rotStiffness(C66) crystal_C66_twin(1:6,1:6,i) = R%rotStiffness(C66)
end do end do
end function lattice_C66_twin end function crystal_C66_twin
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Rotated elasticity matrices for transformation in 6x6-matrix notation !> @brief Rotated elasticity matrices for transformation in 6x6-matrix notation
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_C66_trans(Ntrans,C_parent66,lattice_target, & function crystal_C66_trans(Ntrans,C_parent66,crystal_target, &
cOverA_trans,a_cF,a_cI) cOverA_trans,a_cF,a_cI)
integer, dimension(:), intent(in) :: Ntrans !< number of active twin systems per family integer, dimension(:), intent(in) :: Ntrans !< number of active twin systems per family
character(len=*), intent(in) :: lattice_target !< Bravais lattice (Pearson symbol) character(len=*), intent(in) :: crystal_target !< Bravais lattice (Pearson symbol)
real(pREAL), dimension(6,6), intent(in) :: C_parent66 real(pREAL), dimension(6,6), intent(in) :: C_parent66
real(pREAL), optional, intent(in) :: cOverA_trans, a_cF, a_cI real(pREAL), optional, intent(in) :: cOverA_trans, a_cF, a_cI
real(pREAL), dimension(6,6,sum(Ntrans)) :: lattice_C66_trans real(pREAL), dimension(6,6,sum(Ntrans)) :: crystal_C66_trans
real(pREAL), dimension(6,6) :: C_bar66, C_target_unrotated66 real(pREAL), dimension(6,6) :: C_bar66, C_target_unrotated66
real(pREAL), dimension(3,3,sum(Ntrans)) :: Q,S real(pREAL), dimension(3,3,sum(Ntrans)) :: Q,S
@ -548,11 +548,11 @@ function lattice_C66_trans(Ntrans,C_parent66,lattice_target, &
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! elasticity matrix of the target phase in cube orientation ! elasticity matrix of the target phase in cube orientation
if (lattice_target == 'hP' .and. present(cOverA_trans)) then if (crystal_target == 'hP' .and. present(cOverA_trans)) then
! https://doi.org/10.1063/1.1663858 eq. (16), eq. (18), eq. (19) ! https://doi.org/10.1063/1.1663858 eq. (16), eq. (18), eq. (19)
! https://doi.org/10.1016/j.actamat.2016.07.032 eq. (47), eq. (48) ! https://doi.org/10.1016/j.actamat.2016.07.032 eq. (47), eq. (48)
if (cOverA_trans < 1.0_pREAL .or. cOverA_trans > 2.0_pREAL) & if (cOverA_trans < 1.0_pREAL .or. cOverA_trans > 2.0_pREAL) &
call IO_error(131,ext_msg='lattice_C66_trans: '//trim(lattice_target)) call IO_error(131,ext_msg='crystal_C66_trans: '//trim(crystal_target))
C_bar66(1,1) = (C_parent66(1,1) + C_parent66(1,2) + 2.0_pREAL*C_parent66(4,4))/2.0_pREAL C_bar66(1,1) = (C_parent66(1,1) + C_parent66(1,2) + 2.0_pREAL*C_parent66(4,4))/2.0_pREAL
C_bar66(1,2) = (C_parent66(1,1) + 5.0_pREAL*C_parent66(1,2) - 2.0_pREAL*C_parent66(4,4))/6.0_pREAL C_bar66(1,2) = (C_parent66(1,1) + 5.0_pREAL*C_parent66(1,2) - 2.0_pREAL*C_parent66(4,4))/6.0_pREAL
C_bar66(3,3) = (C_parent66(1,1) + 2.0_pREAL*C_parent66(1,2) + 4.0_pREAL*C_parent66(4,4))/3.0_pREAL C_bar66(3,3) = (C_parent66(1,1) + 2.0_pREAL*C_parent66(1,2) + 4.0_pREAL*C_parent66(4,4))/3.0_pREAL
@ -566,13 +566,13 @@ function lattice_C66_trans(Ntrans,C_parent66,lattice_target, &
C_target_unrotated66(1,3) = C_bar66(1,3) C_target_unrotated66(1,3) = C_bar66(1,3)
C_target_unrotated66(3,3) = C_bar66(3,3) C_target_unrotated66(3,3) = C_bar66(3,3)
C_target_unrotated66(4,4) = C_bar66(4,4) - C_bar66(1,4)**2/(0.5_pREAL*(C_bar66(1,1) - C_bar66(1,2))) C_target_unrotated66(4,4) = C_bar66(4,4) - C_bar66(1,4)**2/(0.5_pREAL*(C_bar66(1,1) - C_bar66(1,2)))
C_target_unrotated66 = lattice_symmetrize_C66(C_target_unrotated66,'hP') C_target_unrotated66 = crystal_symmetrize_C66(C_target_unrotated66,'hP')
elseif (lattice_target == 'cI' .and. present(a_cF) .and. present(a_cI)) then elseif (crystal_target == 'cI' .and. present(a_cF) .and. present(a_cI)) then
if (a_cI <= 0.0_pREAL .or. a_cF <= 0.0_pREAL) & if (a_cI <= 0.0_pREAL .or. a_cF <= 0.0_pREAL) &
call IO_error(134,ext_msg='lattice_C66_trans: '//trim(lattice_target)) call IO_error(134,ext_msg='crystal_C66_trans: '//trim(crystal_target))
C_target_unrotated66 = C_parent66 C_target_unrotated66 = C_parent66
else else
call IO_error(137,ext_msg='lattice_C66_trans : '//trim(lattice_target)) call IO_error(137,ext_msg='crystal_C66_trans : '//trim(crystal_target))
end if end if
do i = 1,6 do i = 1,6
@ -584,10 +584,10 @@ function lattice_C66_trans(Ntrans,C_parent66,lattice_target, &
do i = 1,sum(Ntrans) do i = 1,sum(Ntrans)
call R%fromMatrix(Q(1:3,1:3,i)) call R%fromMatrix(Q(1:3,1:3,i))
lattice_C66_trans(1:6,1:6,i) = R%rotStiffness(C_target_unrotated66) crystal_C66_trans(1:6,1:6,i) = R%rotStiffness(C_target_unrotated66)
end do end do
end function lattice_C66_trans end function crystal_C66_trans
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -595,7 +595,7 @@ function lattice_C66_trans(Ntrans,C_parent66,lattice_target, &
! https://doi.org/10.1016/j.actamat.2012.03.053, eq. (17) ! https://doi.org/10.1016/j.actamat.2012.03.053, eq. (17)
! https://doi.org/10.1016/j.actamat.2008.07.037, table 1 ! https://doi.org/10.1016/j.actamat.2008.07.037, table 1
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_nonSchmidMatrix(Nslip,nonSchmidCoefficients,sense) result(nonSchmidMatrix) function crystal_nonSchmidMatrix(Nslip,nonSchmidCoefficients,sense) result(nonSchmidMatrix)
integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family
real(pREAL), dimension(:), intent(in) :: nonSchmidCoefficients !< non-Schmid coefficients for projections real(pREAL), dimension(:), intent(in) :: nonSchmidCoefficients !< non-Schmid coefficients for projections
@ -608,11 +608,11 @@ function lattice_nonSchmidMatrix(Nslip,nonSchmidCoefficients,sense) result(nonSc
integer :: i integer :: i
if (abs(sense) /= 1) error stop 'Sense in lattice_nonSchmidMatrix' if (abs(sense) /= 1) error stop 'Sense in crystal_nonSchmidMatrix'
coordinateSystem = buildCoordinateSystem(Nslip,CI_NSLIPSYSTEM,CI_SYSTEMSLIP,'cI',0.0_pREAL) coordinateSystem = buildCoordinateSystem(Nslip,CI_NSLIPSYSTEM,CI_SYSTEMSLIP,'cI',0.0_pREAL)
coordinateSystem(1:3,1,1:sum(Nslip)) = coordinateSystem(1:3,1,1:sum(Nslip))*real(sense,pREAL) ! convert unidirectional coordinate system coordinateSystem(1:3,1,1:sum(Nslip)) = coordinateSystem(1:3,1,1:sum(Nslip))*real(sense,pREAL) ! convert unidirectional coordinate system
nonSchmidMatrix = lattice_SchmidMatrix_slip(Nslip,'cI',0.0_pREAL) ! Schmid contribution nonSchmidMatrix = crystal_SchmidMatrix_slip(Nslip,'cI',0.0_pREAL) ! Schmid contribution
do i = 1,sum(Nslip) do i = 1,sum(Nslip)
direction = coordinateSystem(1:3,1,i) direction = coordinateSystem(1:3,1,i)
@ -635,7 +635,7 @@ function lattice_nonSchmidMatrix(Nslip,nonSchmidCoefficients,sense) result(nonSc
+ nonSchmidCoefficients(6) * math_outer(direction, direction) + nonSchmidCoefficients(6) * math_outer(direction, direction)
end do end do
end function lattice_nonSchmidMatrix end function crystal_nonSchmidMatrix
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -644,7 +644,7 @@ end function lattice_nonSchmidMatrix
!> @details https://doi.org/10.1016/j.actamat.2016.12.040 (cF: Tab S4-1, cI: Tab S5-1) !> @details https://doi.org/10.1016/j.actamat.2016.12.040 (cF: Tab S4-1, cI: Tab S5-1)
!> @details https://doi.org/10.1016/j.ijplas.2014.06.010 (hP: Tab 3b) !> @details https://doi.org/10.1016/j.ijplas.2014.06.010 (hP: Tab 3b)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_interaction_SlipBySlip(Nslip,interactionValues,lattice) result(interactionMatrix) function crystal_interaction_SlipBySlip(Nslip,interactionValues,lattice) result(interactionMatrix)
integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family
real(pREAL), dimension(:), intent(in) :: interactionValues !< values for slip-slip interaction real(pREAL), dimension(:), intent(in) :: interactionValues !< values for slip-slip interaction
@ -950,19 +950,19 @@ function lattice_interaction_SlipBySlip(Nslip,interactionValues,lattice) result(
interactionTypes = TI_INTERACTIONSLIPSLIP interactionTypes = TI_INTERACTIONSLIPSLIP
NslipMax = TI_NSLIPSYSTEM NslipMax = TI_NSLIPSYSTEM
case default case default
call IO_error(137,ext_msg='lattice_interaction_SlipBySlip: '//trim(lattice)) call IO_error(137,ext_msg='crystal_interaction_SlipBySlip: '//trim(lattice))
end select end select
interactionMatrix = buildInteraction(Nslip,Nslip,NslipMax,NslipMax,interactionValues,interactionTypes) interactionMatrix = buildInteraction(Nslip,Nslip,NslipMax,NslipMax,interactionValues,interactionTypes)
end function lattice_interaction_SlipBySlip end function crystal_interaction_SlipBySlip
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Twin-twin interaction matrix !> @brief Twin-twin interaction matrix
!> details only active twin systems are considered !> details only active twin systems are considered
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_interaction_TwinByTwin(Ntwin,interactionValues,lattice) result(interactionMatrix) function crystal_interaction_TwinByTwin(Ntwin,interactionValues,lattice) result(interactionMatrix)
integer, dimension(:), intent(in) :: Ntwin !< number of active twin systems per family integer, dimension(:), intent(in) :: Ntwin !< number of active twin systems per family
real(pREAL), dimension(:), intent(in) :: interactionValues !< values for twin-twin interaction real(pREAL), dimension(:), intent(in) :: interactionValues !< values for twin-twin interaction
@ -1049,19 +1049,19 @@ function lattice_interaction_TwinByTwin(Ntwin,interactionValues,lattice) result(
interactionTypes = HP_INTERACTIONTWINTWIN interactionTypes = HP_INTERACTIONTWINTWIN
NtwinMax = HP_NTWINSYSTEM NtwinMax = HP_NTWINSYSTEM
case default case default
call IO_error(137,ext_msg='lattice_interaction_TwinByTwin: '//trim(lattice)) call IO_error(137,ext_msg='crystal_interaction_TwinByTwin: '//trim(lattice))
end select end select
interactionMatrix = buildInteraction(Ntwin,Ntwin,NtwinMax,NtwinMax,interactionValues,interactionTypes) interactionMatrix = buildInteraction(Ntwin,Ntwin,NtwinMax,NtwinMax,interactionValues,interactionTypes)
end function lattice_interaction_TwinByTwin end function crystal_interaction_TwinByTwin
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Trans-trans interaction matrix !> @brief Trans-trans interaction matrix
!> details only active trans systems are considered !> details only active trans systems are considered
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_interaction_TransByTrans(Ntrans,interactionValues,lattice) result(interactionMatrix) function crystal_interaction_TransByTrans(Ntrans,interactionValues,lattice) result(interactionMatrix)
integer, dimension(:), intent(in) :: Ntrans !< number of active trans systems per family integer, dimension(:), intent(in) :: Ntrans !< number of active trans systems per family
real(pREAL), dimension(:), intent(in) :: interactionValues !< values for trans-trans interaction real(pREAL), dimension(:), intent(in) :: interactionValues !< values for trans-trans interaction
@ -1091,19 +1091,19 @@ function lattice_interaction_TransByTrans(Ntrans,interactionValues,lattice) resu
interactionTypes = CF_INTERACTIONTRANSTRANS interactionTypes = CF_INTERACTIONTRANSTRANS
NtransMax = CF_NTRANSSYSTEM NtransMax = CF_NTRANSSYSTEM
else else
call IO_error(137,ext_msg='lattice_interaction_TransByTrans: '//trim(lattice)) call IO_error(137,ext_msg='crystal_interaction_TransByTrans: '//trim(lattice))
end if end if
interactionMatrix = buildInteraction(Ntrans,Ntrans,NtransMax,NtransMax,interactionValues,interactionTypes) interactionMatrix = buildInteraction(Ntrans,Ntrans,NtransMax,NtransMax,interactionValues,interactionTypes)
end function lattice_interaction_TransByTrans end function crystal_interaction_TransByTrans
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Slip-twin interaction matrix !> @brief Slip-twin interaction matrix
!> details only active slip and twin systems are considered !> details only active slip and twin systems are considered
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_interaction_SlipByTwin(Nslip,Ntwin,interactionValues,lattice) result(interactionMatrix) function crystal_interaction_SlipByTwin(Nslip,Ntwin,interactionValues,lattice) result(interactionMatrix)
integer, dimension(:), intent(in) :: Nslip, & !< number of active slip systems per family integer, dimension(:), intent(in) :: Nslip, & !< number of active slip systems per family
Ntwin !< number of active twin systems per family Ntwin !< number of active twin systems per family
@ -1251,19 +1251,19 @@ function lattice_interaction_SlipByTwin(Nslip,Ntwin,interactionValues,lattice) r
NslipMax = HP_NSLIPSYSTEM NslipMax = HP_NSLIPSYSTEM
NtwinMax = HP_NTWINSYSTEM NtwinMax = HP_NTWINSYSTEM
case default case default
call IO_error(137,ext_msg='lattice_interaction_SlipByTwin: '//trim(lattice)) call IO_error(137,ext_msg='crystal_interaction_SlipByTwin: '//trim(lattice))
end select end select
interactionMatrix = buildInteraction(Nslip,Ntwin,NslipMax,NtwinMax,interactionValues,interactionTypes) interactionMatrix = buildInteraction(Nslip,Ntwin,NslipMax,NtwinMax,interactionValues,interactionTypes)
end function lattice_interaction_SlipByTwin end function crystal_interaction_SlipByTwin
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Slip-trans interaction matrix !> @brief Slip-trans interaction matrix
!> details only active slip and trans systems are considered !> details only active slip and trans systems are considered
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_interaction_SlipByTrans(Nslip,Ntrans,interactionValues,lattice) result(interactionMatrix) function crystal_interaction_SlipByTrans(Nslip,Ntrans,interactionValues,lattice) result(interactionMatrix)
integer, dimension(:), intent(in) :: Nslip, & !< number of active slip systems per family integer, dimension(:), intent(in) :: Nslip, & !< number of active slip systems per family
Ntrans !< number of active trans systems per family Ntrans !< number of active trans systems per family
@ -1304,19 +1304,19 @@ function lattice_interaction_SlipByTrans(Nslip,Ntrans,interactionValues,lattice)
NslipMax = CF_NSLIPSYSTEM NslipMax = CF_NSLIPSYSTEM
NtransMax = CF_NTRANSSYSTEM NtransMax = CF_NTRANSSYSTEM
case default case default
call IO_error(137,ext_msg='lattice_interaction_SlipByTrans: '//trim(lattice)) call IO_error(137,ext_msg='crystal_interaction_SlipByTrans: '//trim(lattice))
end select end select
interactionMatrix = buildInteraction(Nslip,Ntrans,NslipMax,NtransMax,interactionValues,interactionTypes) interactionMatrix = buildInteraction(Nslip,Ntrans,NslipMax,NtransMax,interactionValues,interactionTypes)
end function lattice_interaction_SlipByTrans end function crystal_interaction_SlipByTrans
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Twin-slip interaction matrix !> @brief Twin-slip interaction matrix
!> details only active twin and slip systems are considered !> details only active twin and slip systems are considered
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_interaction_TwinBySlip(Ntwin,Nslip,interactionValues,lattice) result(interactionMatrix) function crystal_interaction_TwinBySlip(Ntwin,Nslip,interactionValues,lattice) result(interactionMatrix)
integer, dimension(:), intent(in) :: Ntwin, & !< number of active twin systems per family integer, dimension(:), intent(in) :: Ntwin, & !< number of active twin systems per family
Nslip !< number of active slip systems per family Nslip !< number of active slip systems per family
@ -1380,19 +1380,19 @@ function lattice_interaction_TwinBySlip(Ntwin,Nslip,interactionValues,lattice) r
NtwinMax = HP_NTWINSYSTEM NtwinMax = HP_NTWINSYSTEM
NslipMax = HP_NSLIPSYSTEM NslipMax = HP_NSLIPSYSTEM
case default case default
call IO_error(137,ext_msg='lattice_interaction_TwinBySlip: '//trim(lattice)) call IO_error(137,ext_msg='crystal_interaction_TwinBySlip: '//trim(lattice))
end select end select
interactionMatrix = buildInteraction(Ntwin,Nslip,NtwinMax,NslipMax,interactionValues,interactionTypes) interactionMatrix = buildInteraction(Ntwin,Nslip,NtwinMax,NslipMax,interactionValues,interactionTypes)
end function lattice_interaction_TwinBySlip end function crystal_interaction_TwinBySlip
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Schmid matrix for slip !> @brief Schmid matrix for slip
!> details only active slip systems are considered !> details only active slip systems are considered
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_SchmidMatrix_slip(Nslip,lattice,cOverA) result(SchmidMatrix) function crystal_SchmidMatrix_slip(Nslip,lattice,cOverA) result(SchmidMatrix)
integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family
character(len=*), intent(in) :: lattice !< Bravais lattice (Pearson symbol) character(len=*), intent(in) :: lattice !< Bravais lattice (Pearson symbol)
@ -1419,7 +1419,7 @@ function lattice_SchmidMatrix_slip(Nslip,lattice,cOverA) result(SchmidMatrix)
slipSystems = TI_SYSTEMSLIP slipSystems = TI_SYSTEMSLIP
case default case default
allocate(NslipMax(0)) allocate(NslipMax(0))
call IO_error(137,ext_msg='lattice_SchmidMatrix_slip: '//trim(lattice)) call IO_error(137,ext_msg='crystal_SchmidMatrix_slip: '//trim(lattice))
end select end select
if (any(NslipMax(1:size(Nslip)) - Nslip < 0)) & if (any(NslipMax(1:size(Nslip)) - Nslip < 0)) &
@ -1435,14 +1435,14 @@ function lattice_SchmidMatrix_slip(Nslip,lattice,cOverA) result(SchmidMatrix)
error stop 'dilatational Schmid matrix for slip' error stop 'dilatational Schmid matrix for slip'
end do end do
end function lattice_SchmidMatrix_slip end function crystal_SchmidMatrix_slip
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Schmid matrix for twinning !> @brief Schmid matrix for twinning
!> details only active twin systems are considered !> details only active twin systems are considered
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_SchmidMatrix_twin(Ntwin,lattice,cOverA) result(SchmidMatrix) function crystal_SchmidMatrix_twin(Ntwin,lattice,cOverA) result(SchmidMatrix)
integer, dimension(:), intent(in) :: Ntwin !< number of active twin systems per family integer, dimension(:), intent(in) :: Ntwin !< number of active twin systems per family
character(len=*), intent(in) :: lattice !< Bravais lattice (Pearson symbol) character(len=*), intent(in) :: lattice !< Bravais lattice (Pearson symbol)
@ -1466,7 +1466,7 @@ function lattice_SchmidMatrix_twin(Ntwin,lattice,cOverA) result(SchmidMatrix)
twinSystems = HP_SYSTEMTWIN twinSystems = HP_SYSTEMTWIN
case default case default
allocate(NtwinMax(0)) allocate(NtwinMax(0))
call IO_error(137,ext_msg='lattice_SchmidMatrix_twin: '//trim(lattice)) call IO_error(137,ext_msg='crystal_SchmidMatrix_twin: '//trim(lattice))
end select end select
if (any(NtwinMax(1:size(Ntwin)) - Ntwin < 0)) & if (any(NtwinMax(1:size(Ntwin)) - Ntwin < 0)) &
@ -1482,43 +1482,43 @@ function lattice_SchmidMatrix_twin(Ntwin,lattice,cOverA) result(SchmidMatrix)
error stop 'dilatational Schmid matrix for twin' error stop 'dilatational Schmid matrix for twin'
end do end do
end function lattice_SchmidMatrix_twin end function crystal_SchmidMatrix_twin
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Schmid matrix for transformation !> @brief Schmid matrix for transformation
!> details only active twin systems are considered !> details only active twin systems are considered
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_SchmidMatrix_trans(Ntrans,lattice_target,cOverA,a_cF,a_cI) result(SchmidMatrix) function crystal_SchmidMatrix_trans(Ntrans,crystal_target,cOverA,a_cF,a_cI) result(SchmidMatrix)
integer, dimension(:), intent(in) :: Ntrans !< number of active twin systems per family integer, dimension(:), intent(in) :: Ntrans !< number of active twin systems per family
character(len=*), intent(in) :: lattice_target !< Bravais lattice (Pearson symbol) character(len=*), intent(in) :: crystal_target !< Bravais lattice (Pearson symbol)
real(pREAL), optional, intent(in) :: cOverA, a_cI, a_cF real(pREAL), optional, intent(in) :: cOverA, a_cI, a_cF
real(pREAL), dimension(3,3,sum(Ntrans)) :: SchmidMatrix real(pREAL), dimension(3,3,sum(Ntrans)) :: SchmidMatrix
real(pREAL), dimension(3,3,sum(Ntrans)) :: devNull real(pREAL), dimension(3,3,sum(Ntrans)) :: devNull
if (lattice_target == 'hP' .and. present(cOverA)) then if (crystal_target == 'hP' .and. present(cOverA)) then
if (cOverA < 1.0_pREAL .or. cOverA > 2.0_pREAL) & if (cOverA < 1.0_pREAL .or. cOverA > 2.0_pREAL) &
call IO_error(131,ext_msg='lattice_SchmidMatrix_trans: '//trim(lattice_target)) call IO_error(131,ext_msg='crystal_SchmidMatrix_trans: '//trim(crystal_target))
call buildTransformationSystem(devNull,SchmidMatrix,Ntrans,cOverA=cOverA) call buildTransformationSystem(devNull,SchmidMatrix,Ntrans,cOverA=cOverA)
else if (lattice_target == 'cI' .and. present(a_cF) .and. present(a_cI)) then else if (crystal_target == 'cI' .and. present(a_cF) .and. present(a_cI)) then
if (a_cI <= 0.0_pREAL .or. a_cF <= 0.0_pREAL) & if (a_cI <= 0.0_pREAL .or. a_cF <= 0.0_pREAL) &
call IO_error(134,ext_msg='lattice_SchmidMatrix_trans: '//trim(lattice_target)) call IO_error(134,ext_msg='crystal_SchmidMatrix_trans: '//trim(crystal_target))
call buildTransformationSystem(devNull,SchmidMatrix,Ntrans,a_cF=a_cF,a_cI=a_cI) call buildTransformationSystem(devNull,SchmidMatrix,Ntrans,a_cF=a_cF,a_cI=a_cI)
else else
call IO_error(131,ext_msg='lattice_SchmidMatrix_trans: '//trim(lattice_target)) call IO_error(131,ext_msg='crystal_SchmidMatrix_trans: '//trim(crystal_target))
end if end if
end function lattice_SchmidMatrix_trans end function crystal_SchmidMatrix_trans
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Schmid matrix for cleavage !> @brief Schmid matrix for cleavage
!> details only active cleavage systems are considered !> details only active cleavage systems are considered
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_SchmidMatrix_cleavage(Ncleavage,lattice,cOverA) result(SchmidMatrix) function crystal_SchmidMatrix_cleavage(Ncleavage,lattice,cOverA) result(SchmidMatrix)
integer, dimension(:), intent(in) :: Ncleavage !< number of active cleavage systems per family integer, dimension(:), intent(in) :: Ncleavage !< number of active cleavage systems per family
character(len=*), intent(in) :: lattice !< Bravais lattice (Pearson symbol) character(len=*), intent(in) :: lattice !< Bravais lattice (Pearson symbol)
@ -1539,7 +1539,7 @@ function lattice_SchmidMatrix_cleavage(Ncleavage,lattice,cOverA) result(SchmidMa
cleavageSystems = CI_SYSTEMCLEAVAGE cleavageSystems = CI_SYSTEMCLEAVAGE
case default case default
allocate(NcleavageMax(0)) allocate(NcleavageMax(0))
call IO_error(137,ext_msg='lattice_SchmidMatrix_cleavage: '//trim(lattice)) call IO_error(137,ext_msg='crystal_SchmidMatrix_cleavage: '//trim(lattice))
end select end select
if (any(NcleavageMax(1:size(Ncleavage)) - Ncleavage < 0)) & if (any(NcleavageMax(1:size(Ncleavage)) - Ncleavage < 0)) &
@ -1555,13 +1555,13 @@ function lattice_SchmidMatrix_cleavage(Ncleavage,lattice,cOverA) result(SchmidMa
SchmidMatrix(1:3,1:3,3,i) = math_outer(coordinateSystem(1:3,2,i),coordinateSystem(1:3,2,i)) SchmidMatrix(1:3,1:3,3,i) = math_outer(coordinateSystem(1:3,2,i),coordinateSystem(1:3,2,i))
end do end do
end function lattice_SchmidMatrix_cleavage end function crystal_SchmidMatrix_cleavage
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Slip direction of slip systems (|| b) !> @brief Slip direction of slip systems (|| b)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_slip_direction(Nslip,lattice,cOverA) result(d) function crystal_slip_direction(Nslip,lattice,cOverA) result(d)
integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family
character(len=*), intent(in) :: lattice !< Bravais lattice (Pearson symbol) character(len=*), intent(in) :: lattice !< Bravais lattice (Pearson symbol)
@ -1573,13 +1573,13 @@ function lattice_slip_direction(Nslip,lattice,cOverA) result(d)
coordinateSystem = coordinateSystem_slip(Nslip,lattice,cOverA) coordinateSystem = coordinateSystem_slip(Nslip,lattice,cOverA)
d = coordinateSystem(1:3,1,1:sum(Nslip)) d = coordinateSystem(1:3,1,1:sum(Nslip))
end function lattice_slip_direction end function crystal_slip_direction
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Normal direction of slip systems (|| n) !> @brief Normal direction of slip systems (|| n)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_slip_normal(Nslip,lattice,cOverA) result(n) function crystal_slip_normal(Nslip,lattice,cOverA) result(n)
integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family
character(len=*), intent(in) :: lattice !< Bravais lattice (Pearson symbol) character(len=*), intent(in) :: lattice !< Bravais lattice (Pearson symbol)
@ -1591,13 +1591,13 @@ function lattice_slip_normal(Nslip,lattice,cOverA) result(n)
coordinateSystem = coordinateSystem_slip(Nslip,lattice,cOverA) coordinateSystem = coordinateSystem_slip(Nslip,lattice,cOverA)
n = coordinateSystem(1:3,2,1:sum(Nslip)) n = coordinateSystem(1:3,2,1:sum(Nslip))
end function lattice_slip_normal end function crystal_slip_normal
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Transverse direction of slip systems (|| t = b x n) !> @brief Transverse direction of slip systems (|| t = b x n)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_slip_transverse(Nslip,lattice,cOverA) result(t) function crystal_slip_transverse(Nslip,lattice,cOverA) result(t)
integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family
character(len=*), intent(in) :: lattice !< Bravais lattice (Pearson symbol) character(len=*), intent(in) :: lattice !< Bravais lattice (Pearson symbol)
@ -1609,14 +1609,14 @@ function lattice_slip_transverse(Nslip,lattice,cOverA) result(t)
coordinateSystem = coordinateSystem_slip(Nslip,lattice,cOverA) coordinateSystem = coordinateSystem_slip(Nslip,lattice,cOverA)
t = coordinateSystem(1:3,3,1:sum(Nslip)) t = coordinateSystem(1:3,3,1:sum(Nslip))
end function lattice_slip_transverse end function crystal_slip_transverse
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Labels of slip systems !> @brief Labels of slip systems
!> details only active slip systems are considered !> details only active slip systems are considered
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_labels_slip(Nslip,lattice) result(labels) function crystal_labels_slip(Nslip,lattice) result(labels)
integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family
character(len=*), intent(in) :: lattice !< Bravais lattice (Pearson symbol) character(len=*), intent(in) :: lattice !< Bravais lattice (Pearson symbol)
@ -1640,7 +1640,7 @@ function lattice_labels_slip(Nslip,lattice) result(labels)
NslipMax = TI_NSLIPSYSTEM NslipMax = TI_NSLIPSYSTEM
slipSystems = TI_SYSTEMSLIP slipSystems = TI_SYSTEMSLIP
case default case default
call IO_error(137,ext_msg='lattice_labels_slip: '//trim(lattice)) call IO_error(137,ext_msg='crystal_labels_slip: '//trim(lattice))
end select end select
if (any(NslipMax(1:size(Nslip)) - Nslip < 0)) & if (any(NslipMax(1:size(Nslip)) - Nslip < 0)) &
@ -1650,13 +1650,13 @@ function lattice_labels_slip(Nslip,lattice) result(labels)
labels = getLabels(Nslip,NslipMax,slipSystems) labels = getLabels(Nslip,NslipMax,slipSystems)
end function lattice_labels_slip end function crystal_labels_slip
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Return 3x3 tensor with symmetry according to given Bravais lattice !> @brief Return 3x3 tensor with symmetry according to given Bravais lattice
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function lattice_symmetrize_33(T,lattice) result(T_sym) pure function crystal_symmetrize_33(T,lattice) result(T_sym)
real(pREAL), dimension(3,3) :: T_sym real(pREAL), dimension(3,3) :: T_sym
@ -1677,14 +1677,14 @@ pure function lattice_symmetrize_33(T,lattice) result(T_sym)
T_sym(3,3) = T(3,3) T_sym(3,3) = T(3,3)
end select end select
end function lattice_symmetrize_33 end function crystal_symmetrize_33
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Return stiffness matrix in 6x6 notation with symmetry according to given Bravais lattice !> @brief Return stiffness matrix in 6x6 notation with symmetry according to given Bravais lattice
!> @details J. A. Rayne and B. S. Chandrasekhar Phys. Rev. 120, 1658 Erratum Phys. Rev. 122, 1962 !> @details J. A. Rayne and B. S. Chandrasekhar Phys. Rev. 120, 1658 Erratum Phys. Rev. 122, 1962
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function lattice_symmetrize_C66(C66,lattice) result(C66_sym) pure function crystal_symmetrize_C66(C66,lattice) result(C66_sym)
real(pREAL), dimension(6,6) :: C66_sym real(pREAL), dimension(6,6) :: C66_sym
@ -1723,14 +1723,14 @@ pure function lattice_symmetrize_C66(C66,lattice) result(C66_sym)
end do end do
end do end do
end function lattice_symmetrize_C66 end function crystal_symmetrize_C66
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Labels for twin systems !> @brief Labels for twin systems
!> details only active twin systems are considered !> details only active twin systems are considered
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_labels_twin(Ntwin,lattice) result(labels) function crystal_labels_twin(Ntwin,lattice) result(labels)
integer, dimension(:), intent(in) :: Ntwin !< number of active slip systems per family integer, dimension(:), intent(in) :: Ntwin !< number of active slip systems per family
character(len=*), intent(in) :: lattice !< Bravais lattice (Pearson symbol) character(len=*), intent(in) :: lattice !< Bravais lattice (Pearson symbol)
@ -1751,7 +1751,7 @@ function lattice_labels_twin(Ntwin,lattice) result(labels)
NtwinMax = HP_NTWINSYSTEM NtwinMax = HP_NTWINSYSTEM
twinSystems = HP_SYSTEMTWIN twinSystems = HP_SYSTEMTWIN
case default case default
call IO_error(137,ext_msg='lattice_labels_twin: '//trim(lattice)) call IO_error(137,ext_msg='crystal_labels_twin: '//trim(lattice))
end select end select
if (any(NtwinMax(1:size(Ntwin)) - Ntwin < 0)) & if (any(NtwinMax(1:size(Ntwin)) - Ntwin < 0)) &
@ -1761,7 +1761,7 @@ function lattice_labels_twin(Ntwin,lattice) result(labels)
labels = getLabels(Ntwin,NtwinMax,twinSystems) labels = getLabels(Ntwin,NtwinMax,twinSystems)
end function lattice_labels_twin end function crystal_labels_twin
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -1778,8 +1778,8 @@ function slipProjection_transverse(Nslip,lattice,cOverA) result(projection)
real(pREAL), dimension(3,sum(Nslip)) :: n, t real(pREAL), dimension(3,sum(Nslip)) :: n, t
integer :: i, j integer :: i, j
n = lattice_slip_normal (Nslip,lattice,cOverA) n = crystal_slip_normal (Nslip,lattice,cOverA)
t = lattice_slip_transverse(Nslip,lattice,cOverA) t = crystal_slip_transverse(Nslip,lattice,cOverA)
do i=1, sum(Nslip); do j=1, sum(Nslip) do i=1, sum(Nslip); do j=1, sum(Nslip)
projection(i,j) = abs(math_inner(n(:,i),t(:,j))) projection(i,j) = abs(math_inner(n(:,i),t(:,j)))
@ -1802,8 +1802,8 @@ function slipProjection_direction(Nslip,lattice,cOverA) result(projection)
real(pREAL), dimension(3,sum(Nslip)) :: n, d real(pREAL), dimension(3,sum(Nslip)) :: n, d
integer :: i, j integer :: i, j
n = lattice_slip_normal (Nslip,lattice,cOverA) n = crystal_slip_normal (Nslip,lattice,cOverA)
d = lattice_slip_direction(Nslip,lattice,cOverA) d = crystal_slip_direction(Nslip,lattice,cOverA)
do i=1, sum(Nslip); do j=1, sum(Nslip) do i=1, sum(Nslip); do j=1, sum(Nslip)
projection(i,j) = abs(math_inner(n(:,i),d(:,j))) projection(i,j) = abs(math_inner(n(:,i),d(:,j)))
@ -2150,7 +2150,7 @@ end function getlabels
!> @brief Equivalent Poisson's ratio (ν) !> @brief Equivalent Poisson's ratio (ν)
!> @details https://doi.org/10.1143/JPSJ.20.635 !> @details https://doi.org/10.1143/JPSJ.20.635
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function lattice_isotropic_nu(C,assumption,lattice) result(nu) pure function crystal_isotropic_nu(C,assumption,lattice) result(nu)
real(pREAL), dimension(6,6), intent(in) :: C !< Stiffness tensor (Voigt notation) real(pREAL), dimension(6,6), intent(in) :: C !< Stiffness tensor (Voigt notation)
character(len=*), intent(in) :: assumption !< Assumption (isostrain = 'Voigt', isostress = 'Reuss') character(len=*), intent(in) :: assumption !< Assumption (isostrain = 'Voigt', isostress = 'Reuss')
@ -2172,10 +2172,10 @@ pure function lattice_isotropic_nu(C,assumption,lattice) result(nu)
error stop 'invalid assumption' error stop 'invalid assumption'
end if end if
mu = lattice_isotropic_mu(C,assumption,lattice) mu = crystal_isotropic_mu(C,assumption,lattice)
nu = (1.5_pREAL*K-mu)/(3.0_pREAL*K+mu) nu = (1.5_pREAL*K-mu)/(3.0_pREAL*K+mu)
end function lattice_isotropic_nu end function crystal_isotropic_nu
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -2183,7 +2183,7 @@ end function lattice_isotropic_nu
!> @details https://doi.org/10.1143/JPSJ.20.635 !> @details https://doi.org/10.1143/JPSJ.20.635
!> @details Nonlinear Mechanics of Crystals 10.1007/978-94-007-0350-6, pp 563 !> @details Nonlinear Mechanics of Crystals 10.1007/978-94-007-0350-6, pp 563
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function lattice_isotropic_mu(C,assumption,lattice) result(mu) pure function crystal_isotropic_mu(C,assumption,lattice) result(mu)
real(pREAL), dimension(6,6), intent(in) :: C !< Stiffness tensor (Voigt notation) real(pREAL), dimension(6,6), intent(in) :: C !< Stiffness tensor (Voigt notation)
character(len=*), intent(in) :: assumption !< Assumption (isostrain = 'Voigt', isostress = 'Reuss') character(len=*), intent(in) :: assumption !< Assumption (isostrain = 'Voigt', isostress = 'Reuss')
@ -2220,11 +2220,11 @@ pure function lattice_isotropic_mu(C,assumption,lattice) result(mu)
error stop 'invalid assumption' error stop 'invalid assumption'
end if end if
end function lattice_isotropic_mu end function crystal_isotropic_mu
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Check correctness of some lattice functions. !> @brief Check correctness of some crystal functions.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine selfTest subroutine selfTest
@ -2246,10 +2246,10 @@ subroutine selfTest
do i = 1, 10 do i = 1, 10
call random_number(C) call random_number(C)
C_cF = lattice_symmetrize_C66(C,'cI') C_cF = crystal_symmetrize_C66(C,'cI')
C_cI = lattice_symmetrize_C66(C,'cF') C_cI = crystal_symmetrize_C66(C,'cF')
C_hP = lattice_symmetrize_C66(C,'hP') C_hP = crystal_symmetrize_C66(C,'hP')
C_tI = lattice_symmetrize_C66(C,'tI') C_tI = crystal_symmetrize_C66(C,'tI')
if (any(dNeq(C_cI,transpose(C_cF)))) error stop 'SymmetryC66/cI-cF' if (any(dNeq(C_cI,transpose(C_cF)))) error stop 'SymmetryC66/cI-cF'
if (any(dNeq(C_cF,transpose(C_cI)))) error stop 'SymmetryC66/cF-cI' if (any(dNeq(C_cF,transpose(C_cI)))) error stop 'SymmetryC66/cF-cI'
@ -2269,10 +2269,10 @@ subroutine selfTest
if (any(dNeq(C(4,4),[C_tI(4,4),C_tI(5,5)]))) error stop 'SymmetryC_44-55/tI' if (any(dNeq(C(4,4),[C_tI(4,4),C_tI(5,5)]))) error stop 'SymmetryC_44-55/tI'
call random_number(T) call random_number(T)
T_cF = lattice_symmetrize_33(T,'cI') T_cF = crystal_symmetrize_33(T,'cI')
T_cI = lattice_symmetrize_33(T,'cF') T_cI = crystal_symmetrize_33(T,'cF')
T_hP = lattice_symmetrize_33(T,'hP') T_hP = crystal_symmetrize_33(T,'hP')
T_tI = lattice_symmetrize_33(T,'tI') T_tI = crystal_symmetrize_33(T,'tI')
if (any(dNeq0(T_cF) .and. math_I3<1.0_pREAL)) error stop 'Symmetry33/c' if (any(dNeq0(T_cF) .and. math_I3<1.0_pREAL)) error stop 'Symmetry33/c'
if (any(dNeq0(T_hP) .and. math_I3<1.0_pREAL)) error stop 'Symmetry33/hP' if (any(dNeq0(T_hP) .and. math_I3<1.0_pREAL)) error stop 'Symmetry33/hP'
@ -2291,48 +2291,48 @@ subroutine selfTest
C(4,4) = 0.5_pREAL * (C(1,1) - C(1,2)) C(4,4) = 0.5_pREAL * (C(1,1) - C(1,2))
C(6,6) = C(4,4) C(6,6) = C(4,4)
C_cI = lattice_symmetrize_C66(C,'cI') C_cI = crystal_symmetrize_C66(C,'cI')
if (dNeq(C_cI(4,4),lattice_isotropic_mu(C_cI,'isostrain','cI'),1.0e-12_pREAL)) error stop 'isotropic_mu/isostrain/cI' if (dNeq(C_cI(4,4),crystal_isotropic_mu(C_cI,'isostrain','cI'),1.0e-12_pREAL)) error stop 'isotropic_mu/isostrain/cI'
if (dNeq(C_cI(4,4),lattice_isotropic_mu(C_cI,'isostress','cI'),1.0e-12_pREAL)) error stop 'isotropic_mu/isostress/cI' if (dNeq(C_cI(4,4),crystal_isotropic_mu(C_cI,'isostress','cI'),1.0e-12_pREAL)) error stop 'isotropic_mu/isostress/cI'
lambda = C_cI(1,2) lambda = C_cI(1,2)
if (dNeq(lambda*0.5_pREAL/(lambda+lattice_isotropic_mu(C_cI,'isostrain','cI')), & if (dNeq(lambda*0.5_pREAL/(lambda+crystal_isotropic_mu(C_cI,'isostrain','cI')), &
lattice_isotropic_nu(C_cI,'isostrain','cI'),1.0e-12_pREAL)) error stop 'isotropic_nu/isostrain/cI' crystal_isotropic_nu(C_cI,'isostrain','cI'),1.0e-12_pREAL)) error stop 'isotropic_nu/isostrain/cI'
if (dNeq(lambda*0.5_pREAL/(lambda+lattice_isotropic_mu(C_cI,'isostress','cI')), & if (dNeq(lambda*0.5_pREAL/(lambda+crystal_isotropic_mu(C_cI,'isostress','cI')), &
lattice_isotropic_nu(C_cI,'isostress','cI'),1.0e-12_pREAL)) error stop 'isotropic_nu/isostress/cI' crystal_isotropic_nu(C_cI,'isostress','cI'),1.0e-12_pREAL)) error stop 'isotropic_nu/isostress/cI'
C_hP = lattice_symmetrize_C66(C,'hP') C_hP = crystal_symmetrize_C66(C,'hP')
if (dNeq(C(4,4),lattice_isotropic_mu(C_hP,'isostrain','hP'),1.0e-12_pREAL)) error stop 'isotropic_mu/isostrain/hP' if (dNeq(C(4,4),crystal_isotropic_mu(C_hP,'isostrain','hP'),1.0e-12_pREAL)) error stop 'isotropic_mu/isostrain/hP'
if (dNeq(C(4,4),lattice_isotropic_mu(C_hP,'isostress','hP'),1.0e-12_pREAL)) error stop 'isotropic_mu/isostress/hP' if (dNeq(C(4,4),crystal_isotropic_mu(C_hP,'isostress','hP'),1.0e-12_pREAL)) error stop 'isotropic_mu/isostress/hP'
lambda = C_hP(1,2) lambda = C_hP(1,2)
if (dNeq(lambda*0.5_pREAL/(lambda+lattice_isotropic_mu(C_hP,'isostrain','hP')), & if (dNeq(lambda*0.5_pREAL/(lambda+crystal_isotropic_mu(C_hP,'isostrain','hP')), &
lattice_isotropic_nu(C_hP,'isostrain','hP'),1.0e-12_pREAL)) error stop 'isotropic_nu/isostrain/hP' crystal_isotropic_nu(C_hP,'isostrain','hP'),1.0e-12_pREAL)) error stop 'isotropic_nu/isostrain/hP'
if (dNeq(lambda*0.5_pREAL/(lambda+lattice_isotropic_mu(C_hP,'isostress','hP')), & if (dNeq(lambda*0.5_pREAL/(lambda+crystal_isotropic_mu(C_hP,'isostress','hP')), &
lattice_isotropic_nu(C_hP,'isostress','hP'),1.0e-12_pREAL)) error stop 'isotropic_nu/isostress/hP' crystal_isotropic_nu(C_hP,'isostress','hP'),1.0e-12_pREAL)) error stop 'isotropic_nu/isostress/hP'
C_tI = lattice_symmetrize_C66(C,'tI') C_tI = crystal_symmetrize_C66(C,'tI')
if (dNeq(C(6,6),lattice_isotropic_mu(C_tI,'isostrain','tI'),1.0e-12_pREAL)) error stop 'isotropic_mu/isostrain/tI' if (dNeq(C(6,6),crystal_isotropic_mu(C_tI,'isostrain','tI'),1.0e-12_pREAL)) error stop 'isotropic_mu/isostrain/tI'
if (dNeq(C(6,6),lattice_isotropic_mu(C_tI,'isostress','tI'),1.0e-12_pREAL)) error stop 'isotropic_mu/isostress/tI' if (dNeq(C(6,6),crystal_isotropic_mu(C_tI,'isostress','tI'),1.0e-12_pREAL)) error stop 'isotropic_mu/isostress/tI'
lambda = C_tI(1,2) lambda = C_tI(1,2)
if (dNeq(lambda*0.5_pREAL/(lambda+lattice_isotropic_mu(C_tI,'isostrain','tI')), & if (dNeq(lambda*0.5_pREAL/(lambda+crystal_isotropic_mu(C_tI,'isostrain','tI')), &
lattice_isotropic_nu(C_tI,'isostrain','tI'),1.0e-12_pREAL)) error stop 'isotropic_nu/isostrain/tI' crystal_isotropic_nu(C_tI,'isostrain','tI'),1.0e-12_pREAL)) error stop 'isotropic_nu/isostrain/tI'
if (dNeq(lambda*0.5_pREAL/(lambda+lattice_isotropic_mu(C_tI,'isostress','tI')), & if (dNeq(lambda*0.5_pREAL/(lambda+crystal_isotropic_mu(C_tI,'isostress','tI')), &
lattice_isotropic_nu(C_tI,'isostress','tI'),1.0e-12_pREAL)) error stop 'isotropic_nu/isostress/tI' crystal_isotropic_nu(C_tI,'isostress','tI'),1.0e-12_pREAL)) error stop 'isotropic_nu/isostress/tI'
call random_number(C) call random_number(C)
C = lattice_symmetrize_C66(C+math_eye(6),'cI') C = crystal_symmetrize_C66(C+math_eye(6),'cI')
if (dNeq(lattice_isotropic_mu(C,'isostrain','cI'), lattice_isotropic_mu(C,'isostrain','hP'), 1.0e-12_pREAL)) & if (dNeq(crystal_isotropic_mu(C,'isostrain','cI'), crystal_isotropic_mu(C,'isostrain','hP'), 1.0e-12_pREAL)) &
error stop 'isotropic_mu/isostrain/cI-hP' error stop 'isotropic_mu/isostrain/cI-hP'
if (dNeq(lattice_isotropic_nu(C,'isostrain','cF'), lattice_isotropic_nu(C,'isostrain','cI'), 1.0e-12_pREAL)) & if (dNeq(crystal_isotropic_nu(C,'isostrain','cF'), crystal_isotropic_nu(C,'isostrain','cI'), 1.0e-12_pREAL)) &
error stop 'isotropic_nu/isostrain/cF-tI' error stop 'isotropic_nu/isostrain/cF-tI'
if (dNeq(lattice_isotropic_mu(C,'isostress','cI'), lattice_isotropic_mu(C,'isostress'), 1.0e-12_pREAL)) & if (dNeq(crystal_isotropic_mu(C,'isostress','cI'), crystal_isotropic_mu(C,'isostress'), 1.0e-12_pREAL)) &
error stop 'isotropic_mu/isostress/cI-hP' error stop 'isotropic_mu/isostress/cI-hP'
if (dNeq(lattice_isotropic_nu(C,'isostress','cF'), lattice_isotropic_nu(C,'isostress'), 1.0e-12_pREAL)) & if (dNeq(crystal_isotropic_nu(C,'isostress','cF'), crystal_isotropic_nu(C,'isostress'), 1.0e-12_pREAL)) &
error stop 'isotropic_nu/isostress/cF-tI' error stop 'isotropic_nu/isostress/cF-tI'
end subroutine selfTest end subroutine selfTest
end module lattice end module crystal

6
src/grid/grid_mech_spectral_basic.f90
View File

@ -381,7 +381,7 @@ end subroutine grid_mechanical_spectral_basic_forward
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Update coordinates !> @brief Update coordinates.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine grid_mechanical_spectral_basic_updateCoords() subroutine grid_mechanical_spectral_basic_updateCoords()
@ -390,7 +390,7 @@ subroutine grid_mechanical_spectral_basic_updateCoords()
call DMDAVecGetArrayReadF90(da,solution_vec,F,err_PETSc) call DMDAVecGetArrayReadF90(da,solution_vec,F,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call utilities_updateCoords(F) call utilities_updateCoords(reshape(F,[3,3,size(F,2),size(F,3),size(F,4)]))
call DMDAVecRestoreArrayReadF90(da,solution_vec,F,err_PETSc) call DMDAVecRestoreArrayReadF90(da,solution_vec,F,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
@ -398,7 +398,7 @@ end subroutine grid_mechanical_spectral_basic_updateCoords
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Write current solver and constitutive data for restart to file !> @brief Write current solver and constitutive data for restart to file.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine grid_mechanical_spectral_basic_restartWrite() subroutine grid_mechanical_spectral_basic_restartWrite()

6
src/grid/grid_mech_spectral_polarisation.f90
View File

@ -438,7 +438,7 @@ end subroutine grid_mechanical_spectral_polarisation_forward
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Update coordinates !> @brief Update coordinates.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine grid_mechanical_spectral_polarisation_updateCoords() subroutine grid_mechanical_spectral_polarisation_updateCoords()
@ -447,7 +447,7 @@ subroutine grid_mechanical_spectral_polarisation_updateCoords()
call DMDAVecGetArrayReadF90(da,solution_vec,FandF_tau,err_PETSc) call DMDAVecGetArrayReadF90(da,solution_vec,FandF_tau,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call utilities_updateCoords(FandF_tau(0:8,:,:,:)) call utilities_updateCoords(reshape(FandF_tau(0:8,:,:,:),[3,3,size(FandF_tau,2),size(FandF_tau,3),size(FandF_tau,4)]))
call DMDAVecRestoreArrayReadF90(da,solution_vec,FandF_tau,err_PETSc) call DMDAVecRestoreArrayReadF90(da,solution_vec,FandF_tau,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
@ -455,7 +455,7 @@ end subroutine grid_mechanical_spectral_polarisation_updateCoords
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Write current solver and constitutive data for restart to file !> @brief Write current solver and constitutive data for restart to file.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine grid_mechanical_spectral_polarisation_restartWrite() subroutine grid_mechanical_spectral_polarisation_restartWrite()

2
src/homogenization.f90
View File

@ -16,7 +16,7 @@ module homogenization
use HDF5 use HDF5
use HDF5_utilities use HDF5_utilities
use result use result
use lattice use crystal
implicit none(type,external) implicit none(type,external)
private private

4
src/homogenization_mechanical_RGC.f90
View File

@ -8,7 +8,7 @@
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
submodule(homogenization:mechanical) RGC submodule(homogenization:mechanical) RGC
use rotations use rotations
use lattice use crystal
type :: tParameters type :: tParameters
integer, dimension(:), allocatable :: & integer, dimension(:), allocatable :: &
@ -654,7 +654,7 @@ module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy)
C = phase_homogenizedC66(material_ID_phase(co,ce),material_entry_phase(co,ce)) ! damage not included! C = phase_homogenizedC66(material_ID_phase(co,ce),material_entry_phase(co,ce)) ! damage not included!
equivalentMu = lattice_isotropic_mu(C,'isostrain') equivalentMu = crystal_isotropic_mu(C,'isostrain')
end function equivalentMu end function equivalentMu

4
src/materialpoint.f90
View File

@ -20,7 +20,7 @@ module materialpoint
use rotations use rotations
use polynomials use polynomials
use tables use tables
use lattice use crystal
use material use material
use phase use phase
use homogenization use homogenization
@ -64,7 +64,7 @@ subroutine materialpoint_initAll()
call rotations_init() call rotations_init()
call polynomials_init() call polynomials_init()
call tables_init() call tables_init()
call lattice_init() call crystal_init()
#if defined(MESH) #if defined(MESH)
call discretization_mesh_init(restart=CLI_restartInc>0) call discretization_mesh_init(restart=CLI_restartInc>0)
#elif defined(GRID) #elif defined(GRID)

4
src/phase.f90
View File

@ -14,7 +14,7 @@ module phase
use config use config
use material use material
use result use result
use lattice use crystal
use discretization use discretization
use parallelization use parallelization
use HDF5 use HDF5
@ -336,7 +336,7 @@ module phase
config, & config, &
material, & material, &
result, & result, &
lattice, & crystal, &
discretization, & discretization, &
HDF5_utilities HDF5_utilities
#endif #endif

2
src/phase_damage_anisobrittle.f90
View File

@ -77,7 +77,7 @@ module function anisobrittle_init() result(mySources)
prm%s_crit = src%get_as1dReal('s_crit',requiredSize=size(N_cl)) prm%s_crit = src%get_as1dReal('s_crit',requiredSize=size(N_cl))
prm%g_crit = src%get_as1dReal('g_crit',requiredSize=size(N_cl)) prm%g_crit = src%get_as1dReal('g_crit',requiredSize=size(N_cl))
prm%cleavage_systems = lattice_SchmidMatrix_cleavage(N_cl,phase_lattice(ph),phase_cOverA(ph)) prm%cleavage_systems = crystal_SchmidMatrix_cleavage(N_cl,phase_lattice(ph),phase_cOverA(ph))
! expand: family => system ! expand: family => system
prm%s_crit = math_expand(prm%s_crit,N_cl) prm%s_crit = math_expand(prm%s_crit,N_cl)

2
src/phase_mechanical_eigen_thermalexpansion.f90
View File

@ -92,7 +92,7 @@ module subroutine thermalexpansion_LiAndItsTangent(Li, dLi_dTstar, ph,me)
Alpha = 0.0_pREAL Alpha = 0.0_pREAL
Alpha(1,1) = prm%Alpha_11%at(T) Alpha(1,1) = prm%Alpha_11%at(T)
if (any(phase_lattice(ph) == ['hP','tI'])) Alpha(3,3) = prm%Alpha_33%at(T) if (any(phase_lattice(ph) == ['hP','tI'])) Alpha(3,3) = prm%Alpha_33%at(T)
Alpha = lattice_symmetrize_33(Alpha,phase_lattice(ph)) Alpha = crystal_symmetrize_33(Alpha,phase_lattice(ph))
Li = dot_T * Alpha Li = dot_T * Alpha
end associate end associate

6
src/phase_mechanical_elastic.f90
View File

@ -97,7 +97,7 @@ pure module function elastic_C66(ph,en) result(C66)
if (phase_lattice(ph) == 'tI') C66(6,6) = prm%C_66%at(T) if (phase_lattice(ph) == 'tI') C66(6,6) = prm%C_66%at(T)
C66 = lattice_symmetrize_C66(C66,phase_lattice(ph)) C66 = crystal_symmetrize_C66(C66,phase_lattice(ph))
end associate end associate
@ -119,7 +119,7 @@ pure module function elastic_mu(ph,en,isotropic_bound) result(mu)
associate(prm => param(ph)) associate(prm => param(ph))
mu = lattice_isotropic_mu(elastic_C66(ph,en),isotropic_bound,phase_lattice(ph)) mu = crystal_isotropic_mu(elastic_C66(ph,en),isotropic_bound,phase_lattice(ph))
end associate end associate
@ -141,7 +141,7 @@ pure module function elastic_nu(ph,en,isotropic_bound) result(nu)
associate(prm => param(ph)) associate(prm => param(ph))
nu = lattice_isotropic_nu(elastic_C66(ph,en),isotropic_bound,phase_lattice(ph)) nu = crystal_isotropic_nu(elastic_C66(ph,en),isotropic_bound,phase_lattice(ph))
end associate end associate

14
src/phase_mechanical_plastic_dislotungsten.f90
View File

@ -149,13 +149,13 @@ module function plastic_dislotungsten_init() result(myPlasticity)
N_sl = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray) N_sl = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray)
prm%sum_N_sl = sum(abs(N_sl)) prm%sum_N_sl = sum(abs(N_sl))
slipActive: if (prm%sum_N_sl > 0) then slipActive: if (prm%sum_N_sl > 0) then
prm%systems_sl = lattice_labels_slip(N_sl,phase_lattice(ph)) prm%systems_sl = crystal_labels_slip(N_sl,phase_lattice(ph))
prm%P_sl = lattice_SchmidMatrix_slip(N_sl,phase_lattice(ph),phase_cOverA(ph)) prm%P_sl = crystal_SchmidMatrix_slip(N_sl,phase_lattice(ph),phase_cOverA(ph))
if (phase_lattice(ph) == 'cI') then if (phase_lattice(ph) == 'cI') then
a = pl%get_as1dReal('a_nonSchmid',defaultVal = emptyRealArray) a = pl%get_as1dReal('a_nonSchmid',defaultVal = emptyRealArray)
prm%P_nS_pos = lattice_nonSchmidMatrix(N_sl,a,+1) prm%P_nS_pos = crystal_nonSchmidMatrix(N_sl,a,+1)
prm%P_nS_neg = lattice_nonSchmidMatrix(N_sl,a,-1) prm%P_nS_neg = crystal_nonSchmidMatrix(N_sl,a,-1)
else else
prm%P_nS_pos = prm%P_sl prm%P_nS_pos = prm%P_sl
prm%P_nS_neg = prm%P_sl prm%P_nS_neg = prm%P_sl
@ -184,13 +184,13 @@ module function plastic_dislotungsten_init() result(myPlasticity)
prm%d_caron = prm%b_sl * pl%get_asReal('D_a') prm%d_caron = prm%b_sl * pl%get_asReal('D_a')
prm%f_at = prm%b_sl**3*pl%get_asReal('f_at') prm%f_at = prm%b_sl**3*pl%get_asReal('f_at')
prm%h_sl_sl = lattice_interaction_SlipBySlip(N_sl,pl%get_as1dReal('h_sl-sl'), & prm%h_sl_sl = crystal_interaction_SlipBySlip(N_sl,pl%get_as1dReal('h_sl-sl'), &
phase_lattice(ph)) phase_lattice(ph))
prm%forestProjection = spread( f_edge,1,prm%sum_N_sl) & prm%forestProjection = spread( f_edge,1,prm%sum_N_sl) &
* lattice_forestProjection_edge (N_sl,phase_lattice(ph),phase_cOverA(ph)) & * crystal_forestProjection_edge (N_sl,phase_lattice(ph),phase_cOverA(ph)) &
+ spread(1.0_pREAL-f_edge,1,prm%sum_N_sl) & + spread(1.0_pREAL-f_edge,1,prm%sum_N_sl) &
* lattice_forestProjection_screw(N_sl,phase_lattice(ph),phase_cOverA(ph)) * crystal_forestProjection_screw(N_sl,phase_lattice(ph),phase_cOverA(ph))
! sanity checks ! sanity checks
if ( prm%D_0 < 0.0_pREAL) extmsg = trim(extmsg)//' D_0' if ( prm%D_0 < 0.0_pREAL) extmsg = trim(extmsg)//' D_0'

36
src/phase_mechanical_plastic_dislotwin.f90
View File

@ -73,7 +73,7 @@ submodule(phase:plastic) dislotwin
integer, allocatable, dimension(:,:) :: & integer, allocatable, dimension(:,:) :: &
fcc_twinNucleationSlipPair ! ToDo: Better name? Is also used for trans fcc_twinNucleationSlipPair ! ToDo: Better name? Is also used for trans
character(len=:), allocatable :: & character(len=:), allocatable :: &
lattice_tr, & crystal_tr, &
isotropic_bound isotropic_bound
character(len=pSTRLEN), allocatable, dimension(:) :: & character(len=pSTRLEN), allocatable, dimension(:) :: &
output output
@ -202,9 +202,9 @@ module function plastic_dislotwin_init() result(myPlasticity)
N_sl = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray) N_sl = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray)
prm%sum_N_sl = sum(abs(N_sl)) prm%sum_N_sl = sum(abs(N_sl))
slipActive: if (prm%sum_N_sl > 0) then slipActive: if (prm%sum_N_sl > 0) then
prm%systems_sl = lattice_labels_slip(N_sl,phase_lattice(ph)) prm%systems_sl = crystal_labels_slip(N_sl,phase_lattice(ph))
prm%P_sl = lattice_SchmidMatrix_slip(N_sl,phase_lattice(ph),phase_cOverA(ph)) prm%P_sl = crystal_SchmidMatrix_slip(N_sl,phase_lattice(ph),phase_cOverA(ph))
prm%n0_sl = lattice_slip_normal(N_sl,phase_lattice(ph),phase_cOverA(ph)) prm%n0_sl = crystal_slip_normal(N_sl,phase_lattice(ph),phase_cOverA(ph))
prm%extendedDislocations = pl%get_asBool('extend_dislocations',defaultVal=.false.) prm%extendedDislocations = pl%get_asBool('extend_dislocations',defaultVal=.false.)
prm%omitDipoles = pl%get_asBool('omit_dipoles', defaultVal=.false.) prm%omitDipoles = pl%get_asBool('omit_dipoles', defaultVal=.false.)
@ -226,15 +226,15 @@ module function plastic_dislotwin_init() result(myPlasticity)
defaultVal=[(0.0_pREAL,i=1,size(N_sl))]),N_sl) defaultVal=[(0.0_pREAL,i=1,size(N_sl))]),N_sl)
prm%d_caron = prm%b_sl * pl%get_asReal('D_a') prm%d_caron = prm%b_sl * pl%get_asReal('D_a')
prm%h_sl_sl = lattice_interaction_SlipBySlip(N_sl,pl%get_as1dReal('h_sl-sl'),phase_lattice(ph)) prm%h_sl_sl = crystal_interaction_SlipBySlip(N_sl,pl%get_as1dReal('h_sl-sl'),phase_lattice(ph))
prm%forestProjection = spread( f_edge,1,prm%sum_N_sl) & prm%forestProjection = spread( f_edge,1,prm%sum_N_sl) &
* lattice_forestProjection_edge (N_sl,phase_lattice(ph),phase_cOverA(ph)) & * crystal_forestProjection_edge (N_sl,phase_lattice(ph),phase_cOverA(ph)) &
+ spread(1.0_pREAL-f_edge,1,prm%sum_N_sl) & + spread(1.0_pREAL-f_edge,1,prm%sum_N_sl) &
* lattice_forestProjection_screw(N_sl,phase_lattice(ph),phase_cOverA(ph)) * crystal_forestProjection_screw(N_sl,phase_lattice(ph),phase_cOverA(ph))
prm%fccTwinTransNucleation = phase_lattice(ph) == 'cF' .and. N_sl(1) == 12 prm%fccTwinTransNucleation = phase_lattice(ph) == 'cF' .and. N_sl(1) == 12
if (prm%fccTwinTransNucleation) prm%fcc_twinNucleationSlipPair = lattice_CF_TWINNUCLEATIONSLIPPAIR if (prm%fccTwinTransNucleation) prm%fcc_twinNucleationSlipPair = crystal_CF_TWINNUCLEATIONSLIPPAIR
! multiplication factor according to crystal structure (nearest neighbors bcc vs fcc/hex) ! multiplication factor according to crystal structure (nearest neighbors bcc vs fcc/hex)
! details: Argon & Moffat, Acta Metallurgica, Vol. 29, pg 293 to 299, 1981 ! details: Argon & Moffat, Acta Metallurgica, Vol. 29, pg 293 to 299, 1981
@ -274,9 +274,9 @@ module function plastic_dislotwin_init() result(myPlasticity)
prm%N_tw = pl%get_as1dInt('N_tw', defaultVal=emptyIntArray) prm%N_tw = pl%get_as1dInt('N_tw', defaultVal=emptyIntArray)
prm%sum_N_tw = sum(abs(prm%N_tw)) prm%sum_N_tw = sum(abs(prm%N_tw))
twinActive: if (prm%sum_N_tw > 0) then twinActive: if (prm%sum_N_tw > 0) then
prm%systems_tw = lattice_labels_twin(prm%N_tw,phase_lattice(ph)) prm%systems_tw = crystal_labels_twin(prm%N_tw,phase_lattice(ph))
prm%P_tw = lattice_SchmidMatrix_twin(prm%N_tw,phase_lattice(ph),phase_cOverA(ph)) prm%P_tw = crystal_SchmidMatrix_twin(prm%N_tw,phase_lattice(ph),phase_cOverA(ph))
prm%gamma_char_tw = lattice_characteristicShear_Twin(prm%N_tw,phase_lattice(ph),phase_cOverA(ph)) prm%gamma_char_tw = crystal_characteristicShear_Twin(prm%N_tw,phase_lattice(ph),phase_cOverA(ph))
prm%L_tw = pl%get_asReal('L_tw') prm%L_tw = pl%get_asReal('L_tw')
prm%i_tw = pl%get_asReal('i_tw') prm%i_tw = pl%get_asReal('i_tw')
@ -285,7 +285,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
prm%t_tw = math_expand(pl%get_as1dReal('t_tw', requiredSize=size(prm%N_tw)),prm%N_tw) prm%t_tw = math_expand(pl%get_as1dReal('t_tw', requiredSize=size(prm%N_tw)),prm%N_tw)
prm%r = math_expand(pl%get_as1dReal('p_tw', requiredSize=size(prm%N_tw)),prm%N_tw) prm%r = math_expand(pl%get_as1dReal('p_tw', requiredSize=size(prm%N_tw)),prm%N_tw)
prm%h_tw_tw = lattice_interaction_TwinByTwin(prm%N_tw,pl%get_as1dReal('h_tw-tw'), & prm%h_tw_tw = crystal_interaction_TwinByTwin(prm%N_tw,pl%get_as1dReal('h_tw-tw'), &
phase_lattice(ph)) phase_lattice(ph))
! sanity checks ! sanity checks
@ -309,7 +309,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
prm%N_tr = pl%get_as1dInt('N_tr', defaultVal=emptyIntArray) prm%N_tr = pl%get_as1dInt('N_tr', defaultVal=emptyIntArray)
prm%sum_N_tr = sum(abs(prm%N_tr)) prm%sum_N_tr = sum(abs(prm%N_tr))
transActive: if (prm%sum_N_tr > 0) then transActive: if (prm%sum_N_tr > 0) then
prm%P_tr = lattice_SchmidMatrix_trans(prm%N_tr,'hP',prm%cOverA_hP) prm%P_tr = crystal_SchmidMatrix_trans(prm%N_tr,'hP',prm%cOverA_hP)
prm%Delta_G = polynomial(pl,'Delta_G','T') prm%Delta_G = polynomial(pl,'Delta_G','T')
prm%i_tr = pl%get_asReal('i_tr') prm%i_tr = pl%get_asReal('i_tr')
@ -324,7 +324,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
a_cF = prm%b_tr(1)*sqrt(6.0_pREAL) ! b_tr is Shockley partial a_cF = prm%b_tr(1)*sqrt(6.0_pREAL) ! b_tr is Shockley partial
prm%h = 5.0_pREAL * a_cF/sqrt(3.0_pREAL) prm%h = 5.0_pREAL * a_cF/sqrt(3.0_pREAL)
prm%rho = 4.0_pREAL/(sqrt(3.0_pREAL)*a_cF**2)/N_A prm%rho = 4.0_pREAL/(sqrt(3.0_pREAL)*a_cF**2)/N_A
prm%h_tr_tr = lattice_interaction_TransByTrans(prm%N_tr,pl%get_as1dReal('h_tr-tr'),& prm%h_tr_tr = crystal_interaction_TransByTrans(prm%N_tr,pl%get_as1dReal('h_tr-tr'),&
phase_lattice(ph)) phase_lattice(ph))
@ -372,13 +372,13 @@ module function plastic_dislotwin_init() result(myPlasticity)
prm%Gamma_sf = polynomial(pl,'Gamma_sf','T') prm%Gamma_sf = polynomial(pl,'Gamma_sf','T')
slipAndTwinActive: if (prm%sum_N_sl * prm%sum_N_tw > 0) then slipAndTwinActive: if (prm%sum_N_sl * prm%sum_N_tw > 0) then
prm%h_sl_tw = lattice_interaction_SlipByTwin(N_sl,prm%N_tw,pl%get_as1dReal('h_sl-tw'), & prm%h_sl_tw = crystal_interaction_SlipByTwin(N_sl,prm%N_tw,pl%get_as1dReal('h_sl-tw'), &
phase_lattice(ph)) phase_lattice(ph))
if (prm%fccTwinTransNucleation .and. size(prm%N_tw) /= 1) extmsg = trim(extmsg)//' N_tw: nucleation' if (prm%fccTwinTransNucleation .and. size(prm%N_tw) /= 1) extmsg = trim(extmsg)//' N_tw: nucleation'
end if slipAndTwinActive end if slipAndTwinActive
slipAndTransActive: if (prm%sum_N_sl * prm%sum_N_tr > 0) then slipAndTransActive: if (prm%sum_N_sl * prm%sum_N_tr > 0) then
prm%h_sl_tr = lattice_interaction_SlipByTrans(N_sl,prm%N_tr,pl%get_as1dReal('h_sl-tr'), & prm%h_sl_tr = crystal_interaction_SlipByTrans(N_sl,prm%N_tr,pl%get_as1dReal('h_sl-tr'), &
phase_lattice(ph)) phase_lattice(ph))
if (prm%fccTwinTransNucleation .and. size(prm%N_tr) /= 1) extmsg = trim(extmsg)//' N_tr: nucleation' if (prm%fccTwinTransNucleation .and. size(prm%N_tr) /= 1) extmsg = trim(extmsg)//' N_tr: nucleation'
end if slipAndTransActive end if slipAndTransActive
@ -480,7 +480,7 @@ module function plastic_dislotwin_homogenizedC(ph,en) result(homogenizedC)
homogenizedC = f_matrix * C homogenizedC = f_matrix * C
twinActive: if (prm%sum_N_tw > 0) then twinActive: if (prm%sum_N_tw > 0) then
C66_tw = lattice_C66_twin(prm%N_tw,C,phase_lattice(ph),phase_cOverA(ph)) C66_tw = crystal_C66_twin(prm%N_tw,C,phase_lattice(ph),phase_cOverA(ph))
do i = 1, prm%sum_N_tw do i = 1, prm%sum_N_tw
homogenizedC = homogenizedC & homogenizedC = homogenizedC &
+ stt%f_tw(i,en)*C66_tw(1:6,1:6,i) + stt%f_tw(i,en)*C66_tw(1:6,1:6,i)
@ -488,7 +488,7 @@ module function plastic_dislotwin_homogenizedC(ph,en) result(homogenizedC)
end if twinActive end if twinActive
transActive: if (prm%sum_N_tr > 0) then transActive: if (prm%sum_N_tr > 0) then
C66_tr = lattice_C66_trans(prm%N_tr,C,'hP',prm%cOverA_hP) C66_tr = crystal_C66_trans(prm%N_tr,C,'hP',prm%cOverA_hP)
do i = 1, prm%sum_N_tr do i = 1, prm%sum_N_tr
homogenizedC = homogenizedC & homogenizedC = homogenizedC &
+ stt%f_tr(i,en)*C66_tr(1:6,1:6,i) + stt%f_tr(i,en)*C66_tr(1:6,1:6,i)

10
src/phase_mechanical_plastic_kinehardening.f90
View File

@ -139,14 +139,14 @@ module function plastic_kinehardening_init() result(myPlasticity)
N_sl = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray) N_sl = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray)
prm%sum_N_sl = sum(abs(N_sl)) prm%sum_N_sl = sum(abs(N_sl))
slipActive: if (prm%sum_N_sl > 0) then slipActive: if (prm%sum_N_sl > 0) then
prm%systems_sl = lattice_labels_slip(N_sl,phase_lattice(ph)) prm%systems_sl = crystal_labels_slip(N_sl,phase_lattice(ph))
prm%P = lattice_SchmidMatrix_slip(N_sl,phase_lattice(ph),phase_cOverA(ph)) prm%P = crystal_SchmidMatrix_slip(N_sl,phase_lattice(ph),phase_cOverA(ph))
if (phase_lattice(ph) == 'cI') then if (phase_lattice(ph) == 'cI') then
a = pl%get_as1dReal('a_nonSchmid',defaultVal=emptyRealArray) a = pl%get_as1dReal('a_nonSchmid',defaultVal=emptyRealArray)
prm%nonSchmidActive = size(a) > 0 prm%nonSchmidActive = size(a) > 0
prm%P_nS_pos = lattice_nonSchmidMatrix(N_sl,a,+1) prm%P_nS_pos = crystal_nonSchmidMatrix(N_sl,a,+1)
prm%P_nS_neg = lattice_nonSchmidMatrix(N_sl,a,-1) prm%P_nS_neg = crystal_nonSchmidMatrix(N_sl,a,-1)
else else
prm%P_nS_pos = prm%P prm%P_nS_pos = prm%P
prm%P_nS_neg = prm%P prm%P_nS_neg = prm%P
@ -155,7 +155,7 @@ module function plastic_kinehardening_init() result(myPlasticity)
prm%dot_gamma_0 = pl%get_asReal('dot_gamma_0') prm%dot_gamma_0 = pl%get_asReal('dot_gamma_0')
prm%n = pl%get_asReal('n') prm%n = pl%get_asReal('n')
prm%h_sl_sl = lattice_interaction_SlipBySlip(N_sl,pl%get_as1dReal('h_sl-sl'), & prm%h_sl_sl = crystal_interaction_SlipBySlip(N_sl,pl%get_as1dReal('h_sl-sl'), &
phase_lattice(ph)) phase_lattice(ph))
xi_0 = math_expand(pl%get_as1dReal('xi_0', requiredSize=size(N_sl)),N_sl) xi_0 = math_expand(pl%get_as1dReal('xi_0', requiredSize=size(N_sl)),N_sl)

20
src/phase_mechanical_plastic_nonlocal.f90
View File

@ -249,30 +249,30 @@ module function plastic_nonlocal_init() result(myPlasticity)
ini%N_sl = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray) ini%N_sl = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray)
prm%sum_N_sl = sum(abs(ini%N_sl)) prm%sum_N_sl = sum(abs(ini%N_sl))
slipActive: if (prm%sum_N_sl > 0) then slipActive: if (prm%sum_N_sl > 0) then
prm%systems_sl = lattice_labels_slip(ini%N_sl,phase_lattice(ph)) prm%systems_sl = crystal_labels_slip(ini%N_sl,phase_lattice(ph))
prm%P_sl = lattice_SchmidMatrix_slip(ini%N_sl,phase_lattice(ph), phase_cOverA(ph)) prm%P_sl = crystal_SchmidMatrix_slip(ini%N_sl,phase_lattice(ph), phase_cOverA(ph))
if (phase_lattice(ph) == 'cI') then if (phase_lattice(ph) == 'cI') then
a = pl%get_as1dReal('a_nonSchmid',defaultVal = emptyRealArray) a = pl%get_as1dReal('a_nonSchmid',defaultVal = emptyRealArray)
if (size(a) > 0) prm%nonSchmidActive = .true. if (size(a) > 0) prm%nonSchmidActive = .true.
prm%P_nS_pos = lattice_nonSchmidMatrix(ini%N_sl,a,+1) prm%P_nS_pos = crystal_nonSchmidMatrix(ini%N_sl,a,+1)
prm%P_nS_neg = lattice_nonSchmidMatrix(ini%N_sl,a,-1) prm%P_nS_neg = crystal_nonSchmidMatrix(ini%N_sl,a,-1)
else else
prm%P_nS_pos = prm%P_sl prm%P_nS_pos = prm%P_sl
prm%P_nS_neg = prm%P_sl prm%P_nS_neg = prm%P_sl
end if end if
prm%h_sl_sl = lattice_interaction_SlipBySlip(ini%N_sl,pl%get_as1dReal('h_sl-sl'), & prm%h_sl_sl = crystal_interaction_SlipBySlip(ini%N_sl,pl%get_as1dReal('h_sl-sl'), &
phase_lattice(ph)) phase_lattice(ph))
prm%forestProjection_edge = lattice_forestProjection_edge (ini%N_sl,phase_lattice(ph),& prm%forestProjection_edge = crystal_forestProjection_edge (ini%N_sl,phase_lattice(ph),&
phase_cOverA(ph)) phase_cOverA(ph))
prm%forestProjection_screw = lattice_forestProjection_screw(ini%N_sl,phase_lattice(ph),& prm%forestProjection_screw = crystal_forestProjection_screw(ini%N_sl,phase_lattice(ph),&
phase_cOverA(ph)) phase_cOverA(ph))
prm%slip_direction = lattice_slip_direction (ini%N_sl,phase_lattice(ph),phase_cOverA(ph)) prm%slip_direction = crystal_slip_direction (ini%N_sl,phase_lattice(ph),phase_cOverA(ph))
prm%slip_transverse = lattice_slip_transverse(ini%N_sl,phase_lattice(ph),phase_cOverA(ph)) prm%slip_transverse = crystal_slip_transverse(ini%N_sl,phase_lattice(ph),phase_cOverA(ph))
prm%slip_normal = lattice_slip_normal (ini%N_sl,phase_lattice(ph),phase_cOverA(ph)) prm%slip_normal = crystal_slip_normal (ini%N_sl,phase_lattice(ph),phase_cOverA(ph))
! collinear systems (only for octahedral slip systems in fcc) ! collinear systems (only for octahedral slip systems in fcc)
allocate(prm%colinearSystem(prm%sum_N_sl), source = -1) allocate(prm%colinearSystem(prm%sum_N_sl), source = -1)

22
src/phase_mechanical_plastic_phenopowerlaw.f90
View File

@ -149,21 +149,21 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
prm%h_int = math_expand(pl%get_as1dReal('h_int', requiredSize=size(N_sl), & prm%h_int = math_expand(pl%get_as1dReal('h_int', requiredSize=size(N_sl), &
defaultVal=[(0.0_pREAL,i=1,size(N_sl))]),N_sl) defaultVal=[(0.0_pREAL,i=1,size(N_sl))]),N_sl)
prm%h_sl_sl = lattice_interaction_SlipBySlip(N_sl,pl%get_as1dReal('h_sl-sl'),phase_lattice(ph)) prm%h_sl_sl = crystal_interaction_SlipBySlip(N_sl,pl%get_as1dReal('h_sl-sl'),phase_lattice(ph))
prm%P_sl = lattice_SchmidMatrix_slip(N_sl,phase_lattice(ph),phase_cOverA(ph)) prm%P_sl = crystal_SchmidMatrix_slip(N_sl,phase_lattice(ph),phase_cOverA(ph))
if (phase_lattice(ph) == 'cI') then if (phase_lattice(ph) == 'cI') then
a = pl%get_as1dReal('a_nonSchmid',defaultVal=emptyRealArray) a = pl%get_as1dReal('a_nonSchmid',defaultVal=emptyRealArray)
if (size(a) > 0) prm%nonSchmidActive = .true. if (size(a) > 0) prm%nonSchmidActive = .true.
prm%P_nS_pos = lattice_nonSchmidMatrix(N_sl,a,+1) prm%P_nS_pos = crystal_nonSchmidMatrix(N_sl,a,+1)
prm%P_nS_neg = lattice_nonSchmidMatrix(N_sl,a,-1) prm%P_nS_neg = crystal_nonSchmidMatrix(N_sl,a,-1)
else else
prm%P_nS_pos = prm%P_sl prm%P_nS_pos = prm%P_sl
prm%P_nS_neg = prm%P_sl prm%P_nS_neg = prm%P_sl
end if end if
prm%systems_sl = lattice_labels_slip(N_sl,phase_lattice(ph)) prm%systems_sl = crystal_labels_slip(N_sl,phase_lattice(ph))
! sanity checks ! sanity checks
if ( prm%dot_gamma_0_sl <= 0.0_pREAL) extmsg = trim(extmsg)//' dot_gamma_0_sl' if ( prm%dot_gamma_0_sl <= 0.0_pREAL) extmsg = trim(extmsg)//' dot_gamma_0_sl'
@ -196,11 +196,11 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
xi_0_tw = math_expand(pl%get_as1dReal('xi_0_tw',requiredSize=size(N_tw)),N_tw) xi_0_tw = math_expand(pl%get_as1dReal('xi_0_tw',requiredSize=size(N_tw)),N_tw)
prm%gamma_char = lattice_characteristicShear_twin(N_tw,phase_lattice(ph),phase_cOverA(ph)) prm%gamma_char = crystal_characteristicShear_twin(N_tw,phase_lattice(ph),phase_cOverA(ph))
prm%h_tw_tw = lattice_interaction_TwinByTwin(N_tw,pl%get_as1dReal('h_tw-tw'),phase_lattice(ph)) prm%h_tw_tw = crystal_interaction_TwinByTwin(N_tw,pl%get_as1dReal('h_tw-tw'),phase_lattice(ph))
prm%P_tw = lattice_SchmidMatrix_twin(N_tw,phase_lattice(ph),phase_cOverA(ph)) prm%P_tw = crystal_SchmidMatrix_twin(N_tw,phase_lattice(ph),phase_cOverA(ph))
prm%systems_tw = lattice_labels_twin(N_tw,phase_lattice(ph)) prm%systems_tw = crystal_labels_twin(N_tw,phase_lattice(ph))
! sanity checks ! sanity checks
if (prm%dot_gamma_0_tw <= 0.0_pREAL) extmsg = trim(extmsg)//' dot_gamma_0_tw' if (prm%dot_gamma_0_tw <= 0.0_pREAL) extmsg = trim(extmsg)//' dot_gamma_0_tw'
@ -216,9 +216,9 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
! slip-twin related parameters ! slip-twin related parameters
slipAndTwinActive: if (prm%sum_N_sl > 0 .and. prm%sum_N_tw > 0) then slipAndTwinActive: if (prm%sum_N_sl > 0 .and. prm%sum_N_tw > 0) then
prm%h_0_tw_sl = pl%get_asReal('h_0_tw-sl') prm%h_0_tw_sl = pl%get_asReal('h_0_tw-sl')
prm%h_sl_tw = lattice_interaction_SlipByTwin(N_sl,N_tw,pl%get_as1dReal('h_sl-tw'), & prm%h_sl_tw = crystal_interaction_SlipByTwin(N_sl,N_tw,pl%get_as1dReal('h_sl-tw'), &
phase_lattice(ph)) phase_lattice(ph))
prm%h_tw_sl = lattice_interaction_TwinBySlip(N_tw,N_sl,pl%get_as1dReal('h_tw-sl'), & prm%h_tw_sl = crystal_interaction_TwinBySlip(N_tw,N_sl,pl%get_as1dReal('h_tw-sl'), &
phase_lattice(ph)) phase_lattice(ph))
else slipAndTwinActive else slipAndTwinActive
allocate(prm%h_sl_tw(prm%sum_N_sl,prm%sum_N_tw)) ! at least one dimension is 0 allocate(prm%h_sl_tw(prm%sum_N_sl,prm%sum_N_tw)) ! at least one dimension is 0

2
src/phase_thermal.f90
View File

@ -112,7 +112,7 @@ module subroutine thermal_init(phases)
param(ph)%C_p = thermal%get_asReal('C_p') param(ph)%C_p = thermal%get_asReal('C_p')
param(ph)%K(1,1) = thermal%get_asReal('K_11') param(ph)%K(1,1) = thermal%get_asReal('K_11')
if (any(phase_lattice(ph) == ['hP','tI'])) param(ph)%K(3,3) = thermal%get_asReal('K_33') if (any(phase_lattice(ph) == ['hP','tI'])) param(ph)%K(3,3) = thermal%get_asReal('K_33')
param(ph)%K = lattice_symmetrize_33(param(ph)%K,phase_lattice(ph)) param(ph)%K = crystal_symmetrize_33(param(ph)%K,phase_lattice(ph))
#if defined(__GFORTRAN__) #if defined(__GFORTRAN__)
param(ph)%output = output_as1dStr(thermal) param(ph)%output = output_as1dStr(thermal)