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

Merge branch 'cleaning-pInt-use' into 'development' 

Remove use-only statements

See merge request damask/DAMASK!83
This commit is contained in:
Martin Diehl 2019-05-30 08:48:23 +02:00
parent 144361295d 358272eb2e
commit b25c64d10d
30 changed files with 1049 additions and 1959 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
cmake/Compiler-Intel.cmake
View File

@ -32,6 +32,8 @@
# disables warnings ...
set (COMPILE_FLAGS "${COMPILE_FLAGS} 5268")
# ... the text exceeds right hand column allowed on the line (we have only comments there)
set (COMPILE_FLAGS "${COMPILE_FLAGS},7624")
# ... about deprecated forall (has nice syntax and most likely a performance advantage)
set (COMPILE_FLAGS "${COMPILE_FLAGS} -warn")
# enables warnings ...

7
src/DAMASK_interface.f90
View File

@ -14,7 +14,11 @@
#define PETSC_MAJOR 3
#define PETSC_MINOR_MIN 10
#define PETSC_MINOR_MAX 11
module DAMASK_interface
use, intrinsic :: iso_fortran_env
use PETScSys
use prec
use system_routines
@ -50,9 +54,6 @@ contains
!! information on computation to screen
!--------------------------------------------------------------------------------------------------
subroutine DAMASK_interface_init
use, intrinsic :: iso_fortran_env
use PETScSys
#include <petsc/finclude/petscsys.h>
#if defined(__GFORTRAN__) && __GNUC__<GCC_MIN
===================================================================================================

10
src/HDF5_utilities.f90
View File

@ -8,6 +8,8 @@ module HDF5_utilities
use prec
use IO
use HDF5
use rotations
use numerics
#ifdef PETSc
use PETSC
#endif
@ -1676,8 +1678,6 @@ end subroutine HDF5_write_int7
! ToDo: We could optionally write out other representations (axis angle, euler, ...)
!--------------------------------------------------------------------------------------------------
subroutine HDF5_write_rotation(loc_id,dataset,datasetName,parallel)
use rotations, only: &
rotation
type(rotation), intent(in), dimension(:) :: dataset
integer(HID_T), intent(in) :: loc_id !< file or group handle
@ -1754,9 +1754,6 @@ end subroutine HDF5_write_rotation
subroutine initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, &
myStart, globalShape, &
loc_id,localShape,datasetName,parallel)
use numerics, only: &
worldrank, &
worldsize
integer(HID_T), intent(in) :: loc_id !< file or group handle
character(len=*), intent(in) :: datasetName !< name of the dataset in the file
@ -1850,9 +1847,6 @@ end subroutine finalize_read
subroutine initialize_write(dset_id, filespace_id, memspace_id, plist_id, &
myStart, totalShape, &
loc_id,myShape,datasetName,datatype,parallel)
use numerics, only: &
worldrank, &
worldsize
integer(HID_T), intent(in) :: loc_id !< file or group handle
character(len=*), intent(in) :: datasetName !< name of the dataset in the file

18
src/Lambert.f90
View File

@ -38,11 +38,13 @@
!> Modeling and Simulations in Materials Science and Engineering 22, 075013 (2014).
!--------------------------------------------------------------------------
module Lambert
use prec
use math
implicit none
private
real(pReal), parameter, private :: &
real(pReal), parameter :: &
SPI = sqrt(PI), &
PREF = sqrt(6.0_pReal/PI), &
A = PI**(5.0_pReal/6.0_pReal)/6.0_pReal**(1.0_pReal/6.0_pReal), &
@ -55,10 +57,8 @@ module Lambert
PREK = R1 * 2.0_pReal**(1.0_pReal/4.0_pReal)/BETA
public :: &
LambertCubeToBall, &
LambertBallToCube
private :: &
GetPyramidOrder
Lambert_CubeToBall, &
Lambert_BallToCube
contains
@ -68,7 +68,7 @@ contains
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief map from 3D cubic grid to 3D ball
!--------------------------------------------------------------------------
function LambertCubeToBall(cube) result(ball)
function Lambert_CubeToBall(cube) result(ball)
real(pReal), intent(in), dimension(3) :: cube
real(pReal), dimension(3) :: ball, LamXYZ, XYZ
@ -116,7 +116,7 @@ function LambertCubeToBall(cube) result(ball)
endif center
end function LambertCubeToBall
end function Lambert_CubeToBall
!--------------------------------------------------------------------------
@ -124,7 +124,7 @@ end function LambertCubeToBall
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief map from 3D ball to 3D cubic grid
!--------------------------------------------------------------------------
pure function LambertBallToCube(xyz) result(cube)
pure function Lambert_BallToCube(xyz) result(cube)
real(pReal), intent(in), dimension(3) :: xyz
real(pReal), dimension(3) :: cube, xyz1, xyz3
@ -170,7 +170,7 @@ pure function LambertBallToCube(xyz) result(cube)
endif center
end function LambertBallToCube
end function Lambert_BallToCube
!--------------------------------------------------------------------------

40
src/config.f90
View File

@ -7,10 +7,14 @@
!--------------------------------------------------------------------------------------------------
module config
use prec
use DAMASK_interface
use IO
use debug
use list
implicit none
private
type(tPartitionedStringList), public, protected, allocatable, dimension(:) :: &
config_phase, &
config_microstructure, &
@ -18,10 +22,11 @@ module config
config_texture, &
config_crystallite
type(tPartitionedStringList), public, protected :: &
type(tPartitionedStringList), public, protected :: &
config_numerics, &
config_debug
!ToDo: bad names (how should one know that those variables are defined in config?)
character(len=64), dimension(:), allocatable, public, protected :: &
phase_name, & !< name of each phase
homogenization_name, & !< name of each homogenization
@ -45,19 +50,9 @@ contains
!> @brief reads material.config and stores its content per part
!--------------------------------------------------------------------------------------------------
subroutine config_init
use DAMASK_interface, only: &
getSolverJobName
use IO, only: &
IO_read_ASCII, &
IO_error, &
IO_lc, &
IO_getTag
use debug, only: &
debug_level, &
debug_material, &
debug_levelBasic
integer :: myDebug,i
integer :: i
logical :: verbose
character(len=pStringLen) :: &
line, &
@ -67,7 +62,7 @@ subroutine config_init
write(6,'(/,a)') ' <<<+- config init -+>>>'
myDebug = debug_level(debug_material)
verbose = iand(debug_level(debug_material),debug_levelBasic) /= 0
inquire(file=trim(getSolverJobName())//'.materialConfig',exist=fileExists)
if(fileExists) then
@ -87,23 +82,23 @@ subroutine config_init
case (trim('phase'))
call parse_materialConfig(phase_name,config_phase,line,fileContent(i+1:))
if (iand(myDebug,debug_levelBasic) /= 0) write(6,'(a)') ' Phase parsed'; flush(6)
if (verbose) write(6,'(a)') ' Phase parsed'; flush(6)
case (trim('microstructure'))
call parse_materialConfig(microstructure_name,config_microstructure,line,fileContent(i+1:))
if (iand(myDebug,debug_levelBasic) /= 0) write(6,'(a)') ' Microstructure parsed'; flush(6)
if (verbose) write(6,'(a)') ' Microstructure parsed'; flush(6)
case (trim('crystallite'))
call parse_materialConfig(crystallite_name,config_crystallite,line,fileContent(i+1:))
if (iand(myDebug,debug_levelBasic) /= 0) write(6,'(a)') ' Crystallite parsed'; flush(6)
if (verbose) write(6,'(a)') ' Crystallite parsed'; flush(6)
case (trim('homogenization'))
call parse_materialConfig(homogenization_name,config_homogenization,line,fileContent(i+1:))
if (iand(myDebug,debug_levelBasic) /= 0) write(6,'(a)') ' Homogenization parsed'; flush(6)
if (verbose) write(6,'(a)') ' Homogenization parsed'; flush(6)
case (trim('texture'))
call parse_materialConfig(texture_name,config_texture,line,fileContent(i+1:))
if (iand(myDebug,debug_levelBasic) /= 0) write(6,'(a)') ' Texture parsed'; flush(6)
if (verbose) write(6,'(a)') ' Texture parsed'; flush(6)
end select
@ -141,8 +136,6 @@ contains
!! Recursion is triggered by "{path/to/inputfile}" in a line
!--------------------------------------------------------------------------------------------------
recursive function read_materialConfig(fileName,cnt) result(fileContent)
use IO, only: &
IO_warning
character(len=*), intent(in) :: fileName
integer, intent(in), optional :: cnt !< recursion counter
@ -226,9 +219,6 @@ end function read_materialConfig
subroutine parse_materialConfig(sectionNames,part,line, &
fileContent)
use IO, only: &
IO_intOut
character(len=64), allocatable, dimension(:), intent(out) :: sectionNames
type(tPartitionedStringList), allocatable, dimension(:), intent(inout) :: part
character(len=pStringLen), intent(inout) :: line
@ -298,8 +288,6 @@ end subroutine config_init
!> @brief deallocates the linked lists that store the content of the configuration files
!--------------------------------------------------------------------------------------------------
subroutine config_deallocate(what)
use IO, only: &
IO_error
character(len=*), intent(in) :: what

328
src/crystallite.f90
View File

@ -9,36 +9,43 @@
!--------------------------------------------------------------------------------------------------
module crystallite
use prec, only: &
pReal, &
pStringLen
use rotations, only: &
rotation
use FEsolving, only: &
FEsolving_execElem, &
FEsolving_execIP
use material, only: &
homogenization_Ngrains
use prec
use IO
use config
use debug
use numerics
use rotations
use math
use mesh
use FEsolving
use material
use constitutive
use lattice
use future
use plastic_nonlocal
#if defined(PETSc) || defined(DAMASK_HDF5)
use HDF5_utilities
use results
#endif
implicit none
private
character(len=64), dimension(:,:), allocatable, private :: &
character(len=64), dimension(:,:), allocatable :: &
crystallite_output !< name of each post result output
integer, public, protected :: &
crystallite_maxSizePostResults !< description not available
integer, dimension(:), allocatable, public, protected :: &
crystallite_sizePostResults !< description not available
integer, dimension(:,:), allocatable, private :: &
integer, dimension(:,:), allocatable :: &
crystallite_sizePostResult !< description not available
real(pReal), dimension(:,:,:), allocatable, public :: &
crystallite_dt !< requested time increment of each grain
real(pReal), dimension(:,:,:), allocatable, private :: &
real(pReal), dimension(:,:,:), allocatable :: &
crystallite_subdt, & !< substepped time increment of each grain
crystallite_subFrac, & !< already calculated fraction of increment
crystallite_subStep !< size of next integration step
type(rotation), dimension(:,:,:), allocatable, private :: &
type(rotation), dimension(:,:,:), allocatable :: &
crystallite_orientation, & !< orientation
crystallite_orientation0 !< initial orientation
real(pReal), dimension(:,:,:,:,:), allocatable, public, protected :: &
@ -63,7 +70,7 @@ module crystallite
crystallite_Li, & !< current intermediate velocitiy grad (end of converged time step)
crystallite_Li0, & !< intermediate velocitiy grad at start of FE inc
crystallite_partionedLi0 !< intermediate velocity grad at start of homog inc
real(pReal), dimension(:,:,:,:,:), allocatable, private :: &
real(pReal), dimension(:,:,:,:,:), allocatable :: &
crystallite_subS0, & !< 2nd Piola-Kirchhoff stress vector at start of crystallite inc
crystallite_invFp, & !< inverse of current plastic def grad (end of converged time step)
crystallite_subFp0,& !< plastic def grad at start of crystallite inc
@ -77,7 +84,7 @@ module crystallite
crystallite_dPdF !< current individual dPdF per grain (end of converged time step)
logical, dimension(:,:,:), allocatable, public :: &
crystallite_requested !< used by upper level (homogenization) to request crystallite calculation
logical, dimension(:,:,:), allocatable, private :: &
logical, dimension(:,:,:), allocatable :: &
crystallite_converged, & !< convergence flag
crystallite_todo, & !< flag to indicate need for further computation
crystallite_localPlasticity !< indicates this grain to have purely local constitutive law
@ -101,16 +108,16 @@ module crystallite
neighboringip_ID, &
neighboringelement_ID
end enum
integer(kind(undefined_ID)),dimension(:,:), allocatable, private :: &
integer(kind(undefined_ID)),dimension(:,:), allocatable :: &
crystallite_outputID !< ID of each post result output
type, private :: tOutput !< new requested output (per phase)
type :: tOutput !< new requested output (per phase)
character(len=65536), allocatable, dimension(:) :: &
label
end type tOutput
type(tOutput), allocatable, dimension(:), private :: output_constituent
type(tOutput), allocatable, dimension(:) :: output_constituent
type, private :: tNumerics
type :: tNumerics
integer :: &
iJacoLpresiduum, & !< frequency of Jacobian update of residuum in Lp
nState, & !< state loop limit
@ -138,15 +145,6 @@ module crystallite
crystallite_push33ToRef, &
crystallite_postResults, &
crystallite_results
private :: &
integrateStress, &
integrateState, &
integrateStateFPI, &
integrateStateEuler, &
integrateStateAdaptiveEuler, &
integrateStateRK4, &
integrateStateRKCK45, &
stateJump
contains
@ -155,39 +153,6 @@ contains
!> @brief allocates and initialize per grain variables
!--------------------------------------------------------------------------------------------------
subroutine crystallite_init
#ifdef DEBUG
use debug, only: &
debug_info, &
debug_reset, &
debug_level, &
debug_crystallite, &
debug_levelBasic
#endif
use numerics, only: &
numerics_integrator, &
worldrank, &
usePingPong
use math, only: &
math_I3, &
math_EulerToR, &
math_inv33
use mesh, only: &
theMesh, &
mesh_element
use IO, only: &
IO_stringValue, &
IO_write_jobFile, &
IO_error
use material
use config, only: &
config_deallocate, &
config_crystallite, &
config_numerics, &
config_phase, &
crystallite_name
use constitutive, only: &
constitutive_initialFi, &
constitutive_microstructure ! derived (shortcut) quantities of given state
integer, parameter :: FILEUNIT=434
logical, dimension(:,:), allocatable :: devNull
@ -478,34 +443,6 @@ end subroutine crystallite_init
!> @brief calculate stress (P)
!--------------------------------------------------------------------------------------------------
function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC)
use prec, only: &
tol_math_check, &
dNeq0
#ifdef DEBUG
use debug, only: &
debug_level, &
debug_crystallite, &
debug_levelBasic, &
debug_levelExtensive, &
debug_levelSelective, &
debug_e, &
debug_i, &
debug_g
#endif
use IO, only: &
IO_warning, &
IO_error
use math, only: &
math_inv33
use mesh, only: &
theMesh, &
mesh_element
use material, only: &
homogenization_Ngrains, &
plasticState, &
sourceState, &
phase_Nsources, &
phaseAt, phasememberAt
logical, dimension(theMesh%elem%nIPs,theMesh%Nelems) :: crystallite_stress
real(pReal), intent(in), optional :: &
@ -746,30 +683,6 @@ end function crystallite_stress
!> @brief calculate tangent (dPdF)
!--------------------------------------------------------------------------------------------------
subroutine crystallite_stressTangent
use prec, only: &
tol_math_check, &
dNeq0
use IO, only: &
IO_warning, &
IO_error
use math, only: &
math_inv33, &
math_identity2nd, &
math_3333to99, &
math_99to3333, &
math_I3, &
math_mul3333xx3333, &
math_mul33xx33, &
math_invert2, &
math_det33
use mesh, only: &
mesh_element
use material, only: &
homogenization_Ngrains
use constitutive, only: &
constitutive_SandItsTangents, &
constitutive_LpAndItsTangents, &
constitutive_LiAndItsTangents
integer :: &
c, & !< counter in integration point component loop
@ -910,19 +823,6 @@ end subroutine crystallite_stressTangent
!> @brief calculates orientations
!--------------------------------------------------------------------------------------------------
subroutine crystallite_orientations
use math, only: &
math_rotationalPart33, &
math_RtoQ
use material, only: &
plasticState, &
material_phase, &
homogenization_Ngrains
use mesh, only: &
mesh_element
use lattice, only: &
lattice_qDisorientation
use plastic_nonlocal, only: &
plastic_nonlocal_updateCompatibility
integer &
c, & !< counter in integration point component loop
@ -979,28 +879,6 @@ end function crystallite_push33ToRef
!> @brief return results of particular grain
!--------------------------------------------------------------------------------------------------
function crystallite_postResults(ipc, ip, el)
use math, only: &
math_det33, &
math_I3, &
inDeg
use mesh, only: &
theMesh, &
mesh_element, &
mesh_ipVolume, &
mesh_ipNeighborhood
use material, only: &
plasticState, &
sourceState, &
microstructure_crystallite, &
crystallite_Noutput, &
material_phase, &
material_texture, &
homogenization_Ngrains
use constitutive, only: &
constitutive_homogenizedC, &
constitutive_postResults
use rotations, only: &
rotation
integer, intent(in):: &
el, & !< element index
@ -1118,16 +996,9 @@ end function crystallite_postResults
!--------------------------------------------------------------------------------------------------
subroutine crystallite_results
#if defined(PETSc) || defined(DAMASK_HDF5)
use lattice
use results
use HDF5_utilities
use rotations
use config, only: &
config_name_phase => phase_name ! anticipate logical name
use material, only: &
material_phase_plasticity_type => phase_plasticity
integer :: p,o
real(pReal), allocatable, dimension(:,:,:) :: selected_tensors
type(rotation), allocatable, dimension(:) :: selected_rotations
@ -1267,33 +1138,6 @@ end subroutine crystallite_results
!> intermediate acceleration of the Newton-Raphson correction
!--------------------------------------------------------------------------------------------------
logical function integrateStress(ipc,ip,el,timeFraction)
use, intrinsic :: &
IEEE_arithmetic
use prec, only: tol_math_check, &
dEq0
#ifdef DEBUG
use debug, only: debug_level, &
debug_e, &
debug_i, &
debug_g, &
debug_crystallite, &
debug_levelBasic, &
debug_levelExtensive, &
debug_levelSelective
#endif
use constitutive, only: constitutive_LpAndItsTangents, &
constitutive_LiAndItsTangents, &
constitutive_SandItsTangents
use math, only: math_mul33xx33, &
math_mul3333xx3333, &
math_inv33, &
math_det33, &
math_I3, &
math_identity2nd, &
math_3333to99, &
math_33to9, &
math_9to33
integer, intent(in):: el, & ! element index
ip, & ! integration point index
@ -1693,27 +1537,6 @@ end function integrateStress
!> using Fixed Point Iteration to adapt the stepsize
!--------------------------------------------------------------------------------------------------
subroutine integrateStateFPI
#ifdef DEBUG
use debug, only: debug_level, &
debug_e, &
debug_i, &
debug_g, &
debug_crystallite, &
debug_levelBasic, &
debug_levelExtensive, &
debug_levelSelective
#endif
use mesh, only: &
mesh_element
use material, only: &
plasticState, &
sourceState, &
phaseAt, phasememberAt, &
phase_Nsources, &
homogenization_Ngrains
use constitutive, only: &
constitutive_plasticity_maxSizeDotState, &
constitutive_source_maxSizeDotState
integer :: &
NiterationState, & !< number of iterations in state loop
@ -1901,8 +1724,6 @@ end subroutine integrateStateFPI
!> @brief integrate state with 1st order explicit Euler method
!--------------------------------------------------------------------------------------------------
subroutine integrateStateEuler
use material, only: &
plasticState
call update_dotState(1.0_pReal)
call update_state(1.0_pReal)
@ -1919,19 +1740,6 @@ end subroutine integrateStateEuler
!> @brief integrate stress, state with 1st order Euler method with adaptive step size
!--------------------------------------------------------------------------------------------------
subroutine integrateStateAdaptiveEuler
use mesh, only: &
theMesh, &
mesh_element
use material, only: &
homogenization_Ngrains, &
plasticState, &
sourceState, &
phaseAt, phasememberAt, &
phase_Nsources, &
homogenization_maxNgrains
use constitutive, only: &
constitutive_plasticity_maxSizeDotState, &
constitutive_source_maxSizeDotState
integer :: &
e, & ! element index in element loop
@ -2025,14 +1833,6 @@ end subroutine integrateStateAdaptiveEuler
! ToDo: This is totally BROKEN: RK4dotState is never used!!!
!--------------------------------------------------------------------------------------------------
subroutine integrateStateRK4
use mesh, only: &
mesh_element
use material, only: &
homogenization_Ngrains, &
plasticState, &
sourceState, &
phase_Nsources, &
phaseAt, phasememberAt
real(pReal), dimension(4), parameter :: &
TIMESTEPFRACTION = [0.5_pReal, 0.5_pReal, 1.0_pReal, 1.0_pReal] ! factor giving the fraction of the original timestep used for Runge Kutta Integration
@ -2092,19 +1892,6 @@ end subroutine integrateStateRK4
!> adaptive step size (use 5th order solution to advance = "local extrapolation")
!--------------------------------------------------------------------------------------------------
subroutine integrateStateRKCK45
use mesh, only: &
mesh_element, &
theMesh
use material, only: &
homogenization_Ngrains, &
plasticState, &
sourceState, &
phase_Nsources, &
phaseAt, phasememberAt, &
homogenization_maxNgrains
use constitutive, only: &
constitutive_plasticity_maxSizeDotState, &
constitutive_source_maxSizeDotState
real(pReal), dimension(5,5), parameter :: &
A = reshape([&
@ -2287,8 +2074,6 @@ end subroutine nonlocalConvergenceCheck
!> @details: For explicitEuler, RK4 and RKCK45, adaptive Euler and FPI have their on criteria
!--------------------------------------------------------------------------------------------------
subroutine setConvergenceFlag
use mesh, only: &
mesh_element
integer :: &
e, & !< element index in element loop
@ -2327,8 +2112,6 @@ end subroutine setConvergenceFlag
!> @brief Standard forwarding of state as state = state0 + dotState * (delta t)
!--------------------------------------------------------------------------------------------------
subroutine update_stress(timeFraction)
use mesh, only: &
mesh_element
real(pReal), intent(in) :: &
timeFraction
@ -2360,8 +2143,6 @@ end subroutine update_stress
!> @brief tbd
!--------------------------------------------------------------------------------------------------
subroutine update_dependentState
use mesh, only: &
mesh_element
use constitutive, only: &
constitutive_dependentState => constitutive_microstructure
@ -2387,13 +2168,6 @@ end subroutine update_dependentState
!> @brief Standard forwarding of state as state = state0 + dotState * (delta t)
!--------------------------------------------------------------------------------------------------
subroutine update_state(timeFraction)
use material, only: &
plasticState, &
sourceState, &
phase_Nsources, &
phaseAt, phasememberAt
use mesh, only: &
mesh_element
real(pReal), intent(in) :: &
timeFraction
@ -2435,17 +2209,6 @@ end subroutine update_state
!> if NaN occurs, crystallite_todo is set to FALSE. Any NaN in a nonlocal propagates to all others
!--------------------------------------------------------------------------------------------------
subroutine update_dotState(timeFraction)
use, intrinsic :: &
IEEE_arithmetic
use material, only: &
plasticState, &
sourceState, &
phaseAt, phasememberAt, &
phase_Nsources
use mesh, only: &
mesh_element
use constitutive, only: &
constitutive_collectDotState
real(pReal), intent(in) :: &
timeFraction
@ -2492,19 +2255,7 @@ end subroutine update_DotState
subroutine update_deltaState
use, intrinsic :: &
IEEE_arithmetic
use prec, only: &
dNeq0
use mesh, only: &
mesh_element
use material, only: &
plasticState, &
sourceState, &
phase_Nsources, &
phaseAt, phasememberAt
use constitutive, only: &
constitutive_collectDeltaState
integer :: &
e, & !< element index in element loop
i, & !< integration point index in ip loop
@ -2569,29 +2320,6 @@ end subroutine update_deltaState
!> returns true, if state jump was successfull or not needed. false indicates NaN in delta state
!--------------------------------------------------------------------------------------------------
logical function stateJump(ipc,ip,el)
use, intrinsic :: &
IEEE_arithmetic
use prec, only: &
dNeq0
#ifdef DEBUG
use debug, only: &
debug_e, &
debug_i, &
debug_g, &
debug_level, &
debug_crystallite, &
debug_levelExtensive, &
debug_levelSelective
#endif
use material, only: &
plasticState, &
sourceState, &
phase_Nsources, &
phaseAt, phasememberAt
use mesh, only: &
mesh_element
use constitutive, only: &
constitutive_collectDeltaState
integer, intent(in):: &
el, & ! element index

57
src/damage_local.f90
View File

@ -4,9 +4,13 @@
!--------------------------------------------------------------------------------------------------
module damage_local
use prec
use material
use numerics
use config
implicit none
private
integer, dimension(:,:), allocatable, target, public :: &
damage_local_sizePostResult !< size of each post result output
@ -20,23 +24,22 @@ module damage_local
enumerator :: undefined_ID, &
damage_ID
end enum
integer(kind(undefined_ID)), dimension(:,:), allocatable, private :: &
integer(kind(undefined_ID)), dimension(:,:), allocatable :: &
damage_local_outputID !< ID of each post result output
type, private :: tParameters
type :: tParameters
integer(kind(undefined_ID)), dimension(:), allocatable :: &
outputID
end type tParameters
type(tparameters), dimension(:), allocatable, private :: &
type(tparameters), dimension(:), allocatable :: &
param
public :: &
damage_local_init, &
damage_local_updateState, &
damage_local_postResults
private :: &
damage_local_getSourceAndItsTangent
contains
@ -45,23 +48,8 @@ contains
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
subroutine damage_local_init
use material, only: &
damage_type, &
damage_typeInstance, &
homogenization_Noutput, &
DAMAGE_local_label, &
DAMAGE_local_ID, &
material_homogenizationAt, &
mappingHomogenization, &
damageState, &
damageMapping, &
damage, &
damage_initialPhi
use config, only: &
config_homogenization
integer :: maxNinstance,homog,instance,o,i
integer :: maxNinstance,homog,instance,i
integer :: sizeState
integer :: NofMyHomog, h
integer(kind(undefined_ID)) :: &
@ -72,7 +60,7 @@ subroutine damage_local_init
write(6,'(/,a)') ' <<<+- damage_'//DAMAGE_local_label//' init -+>>>'
maxNinstance = int(count(damage_type == DAMAGE_local_ID),pInt)
maxNinstance = count(damage_type == DAMAGE_local_ID)
if (maxNinstance == 0) return
allocate(damage_local_sizePostResult (maxval(homogenization_Noutput),maxNinstance),source=0)
@ -135,14 +123,6 @@ end subroutine damage_local_init
!> @brief calculates local change in damage field
!--------------------------------------------------------------------------------------------------
function damage_local_updateState(subdt, ip, el)
use numerics, only: &
residualStiffness, &
err_damage_tolAbs, &
err_damage_tolRel
use material, only: &
material_homogenizationAt, &
mappingHomogenization, &
damageState
integer, intent(in) :: &
ip, & !< integration point number
@ -177,17 +157,6 @@ end function damage_local_updateState
!> @brief calculates homogenized local damage driving forces
!--------------------------------------------------------------------------------------------------
subroutine damage_local_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el)
use material, only: &
homogenization_Ngrains, &
material_homogenizationAt, &
phaseAt, &
phasememberAt, &
phase_source, &
phase_Nsources, &
SOURCE_damage_isoBrittle_ID, &
SOURCE_damage_isoDuctile_ID, &
SOURCE_damage_anisoBrittle_ID, &
SOURCE_damage_anisoDuctile_ID
use source_damage_isoBrittle, only: &
source_damage_isobrittle_getRateAndItsTangent
use source_damage_isoDuctile, only: &
@ -244,15 +213,11 @@ subroutine damage_local_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el
end subroutine damage_local_getSourceAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief return array of damage results
!--------------------------------------------------------------------------------------------------
function damage_local_postResults(ip,el)
use material, only: &
material_homogenizationAt, &
damage_typeInstance, &
damageMapping, &
damage
integer, intent(in) :: &
ip, & !< integration point

16
src/damage_none.f90
View File

@ -3,6 +3,8 @@
!> @brief material subroutine for constant damage field
!--------------------------------------------------------------------------------------------------
module damage_none
use config
use material
implicit none
private
@ -15,18 +17,8 @@ contains
!--------------------------------------------------------------------------------------------------
!> @brief allocates all neccessary fields, reads information from material configuration file
!--------------------------------------------------------------------------------------------------
subroutine damage_none_init()
use config, only: &
config_homogenization
use material, only: &
damage_initialPhi, &
damage, &
damage_type, &
material_homogenizationAt, &
damageState, &
DAMAGE_NONE_LABEL, &
DAMAGE_NONE_ID
subroutine damage_none_init
integer :: &
homog, &
NofMyHomog

450
src/damage_nonlocal.f90
View File

@ -4,39 +4,50 @@
!> @details to be done
!--------------------------------------------------------------------------------------------------
module damage_nonlocal
use prec
use prec
use material
use numerics
use config
use crystallite
use lattice
use mesh
use source_damage_isoBrittle
use source_damage_isoDuctile
use source_damage_anisoBrittle
use source_damage_anisoDuctile
implicit none
private
integer, dimension(:,:), allocatable, target, public :: &
damage_nonlocal_sizePostResult !< size of each post result output
implicit none
private
integer, dimension(:,:), allocatable, target, public :: &
damage_nonlocal_sizePostResult !< size of each post result output
character(len=64), dimension(:,:), allocatable, target, public :: &
damage_nonlocal_output !< name of each post result output
integer, dimension(:), allocatable, target, public :: &
damage_nonlocal_Noutput !< number of outputs per instance of this damage
character(len=64), dimension(:,:), allocatable, target, public :: &
damage_nonlocal_output !< name of each post result output
integer, dimension(:), allocatable, target, public :: &
damage_nonlocal_Noutput !< number of outputs per instance of this damage
enum, bind(c)
enumerator :: undefined_ID, &
damage_ID
end enum
enum, bind(c)
enumerator :: undefined_ID, &
damage_ID
end enum
type, private :: tParameters
integer(kind(undefined_ID)), dimension(:), allocatable :: &
outputID
end type tParameters
type(tparameters), dimension(:), allocatable, private :: &
param
type :: tParameters
integer(kind(undefined_ID)), dimension(:), allocatable :: &
outputID
end type tParameters
type(tparameters), dimension(:), allocatable :: &
param
public :: &
damage_nonlocal_init, &
damage_nonlocal_getSourceAndItsTangent, &
damage_nonlocal_getDiffusion33, &
damage_nonlocal_getMobility, &
damage_nonlocal_putNonLocalDamage, &
damage_nonlocal_postResults
public :: &
damage_nonlocal_init, &
damage_nonlocal_getSourceAndItsTangent, &
damage_nonlocal_getDiffusion33, &
damage_nonlocal_getMobility, &
damage_nonlocal_putNonLocalDamage, &
damage_nonlocal_postResults
contains
@ -45,283 +56,228 @@ contains
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
subroutine damage_nonlocal_init
use material, only: &
damage_type, &
damage_typeInstance, &
homogenization_Noutput, &
DAMAGE_nonlocal_label, &
DAMAGE_nonlocal_ID, &
material_homogenizationAt, &
mappingHomogenization, &
damageState, &
damageMapping, &
damage, &
damage_initialPhi
use config, only: &
config_homogenization
integer :: maxNinstance,homog,instance,o,i
integer :: sizeState
integer :: NofMyHomog, h
integer(kind(undefined_ID)) :: &
outputID
character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::]
character(len=65536), dimension(:), allocatable :: &
outputs
integer :: maxNinstance,homog,instance,o,i
integer :: sizeState
integer :: NofMyHomog, h
integer(kind(undefined_ID)) :: &
outputID
character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::]
character(len=65536), dimension(:), allocatable :: &
outputs
write(6,'(/,a)') ' <<<+- damage_'//DAMAGE_nonlocal_label//' init -+>>>'
maxNinstance = int(count(damage_type == DAMAGE_nonlocal_ID))
if (maxNinstance == 0) return
allocate(damage_nonlocal_sizePostResult (maxval(homogenization_Noutput),maxNinstance),source=0)
allocate(damage_nonlocal_output (maxval(homogenization_Noutput),maxNinstance))
damage_nonlocal_output = ''
allocate(damage_nonlocal_Noutput (maxNinstance), source=0)
allocate(param(maxNinstance))
write(6,'(/,a)') ' <<<+- damage_'//DAMAGE_nonlocal_label//' init -+>>>'
do h = 1, size(damage_type)
if (damage_type(h) /= DAMAGE_NONLOCAL_ID) cycle
associate(prm => param(damage_typeInstance(h)), &
config => config_homogenization(h))
instance = damage_typeInstance(h)
outputs = config%getStrings('(output)',defaultVal=emptyStringArray)
allocate(prm%outputID(0))
maxNinstance = count(damage_type == DAMAGE_nonlocal_ID)
if (maxNinstance == 0) return
allocate(damage_nonlocal_sizePostResult (maxval(homogenization_Noutput),maxNinstance),source=0)
allocate(damage_nonlocal_output (maxval(homogenization_Noutput),maxNinstance))
damage_nonlocal_output = ''
allocate(damage_nonlocal_Noutput (maxNinstance), source=0)
allocate(param(maxNinstance))
do i=1, size(outputs)
outputID = undefined_ID
select case(outputs(i))
case ('damage')
damage_nonlocal_output(i,damage_typeInstance(h)) = outputs(i)
damage_nonlocal_Noutput(instance) = damage_nonlocal_Noutput(instance) + 1
damage_nonlocal_sizePostResult(i,damage_typeInstance(h)) = 1
prm%outputID = [prm%outputID , damage_ID]
end select
enddo
do h = 1, size(damage_type)
if (damage_type(h) /= DAMAGE_NONLOCAL_ID) cycle
associate(prm => param(damage_typeInstance(h)), &
config => config_homogenization(h))
instance = damage_typeInstance(h)
outputs = config%getStrings('(output)',defaultVal=emptyStringArray)
allocate(prm%outputID(0))
do i=1, size(outputs)
outputID = undefined_ID
select case(outputs(i))
case ('damage')
damage_nonlocal_output(i,damage_typeInstance(h)) = outputs(i)
damage_nonlocal_Noutput(instance) = damage_nonlocal_Noutput(instance) + 1
damage_nonlocal_sizePostResult(i,damage_typeInstance(h)) = 1
prm%outputID = [prm%outputID , damage_ID]
end select
enddo
homog = h
homog = h
NofMyHomog = count(material_homogenizationAt == homog)
instance = damage_typeInstance(homog)
NofMyHomog = count(material_homogenizationAt == homog)
instance = damage_typeInstance(homog)
! allocate state arrays
sizeState = 1
damageState(homog)%sizeState = sizeState
damageState(homog)%sizePostResults = sum(damage_nonlocal_sizePostResult(:,instance))
allocate(damageState(homog)%state0 (sizeState,NofMyHomog), source=damage_initialPhi(homog))
allocate(damageState(homog)%subState0(sizeState,NofMyHomog), source=damage_initialPhi(homog))
allocate(damageState(homog)%state (sizeState,NofMyHomog), source=damage_initialPhi(homog))
! allocate state arrays
sizeState = 1
damageState(homog)%sizeState = sizeState
damageState(homog)%sizePostResults = sum(damage_nonlocal_sizePostResult(:,instance))
allocate(damageState(homog)%state0 (sizeState,NofMyHomog), source=damage_initialPhi(homog))
allocate(damageState(homog)%subState0(sizeState,NofMyHomog), source=damage_initialPhi(homog))
allocate(damageState(homog)%state (sizeState,NofMyHomog), source=damage_initialPhi(homog))
nullify(damageMapping(homog)%p)
damageMapping(homog)%p => mappingHomogenization(1,:,:)
deallocate(damage(homog)%p)
damage(homog)%p => damageState(homog)%state(1,:)
end associate
enddo
nullify(damageMapping(homog)%p)
damageMapping(homog)%p => mappingHomogenization(1,:,:)
deallocate(damage(homog)%p)
damage(homog)%p => damageState(homog)%state(1,:)
end associate
enddo
end subroutine damage_nonlocal_init
!--------------------------------------------------------------------------------------------------
!> @brief calculates homogenized damage driving forces
!--------------------------------------------------------------------------------------------------
subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el)
use material, only: &
homogenization_Ngrains, &
material_homogenizationAt, &
phaseAt, &
phasememberAt, &
phase_source, &
phase_Nsources, &
SOURCE_damage_isoBrittle_ID, &
SOURCE_damage_isoDuctile_ID, &
SOURCE_damage_anisoBrittle_ID, &
SOURCE_damage_anisoDuctile_ID
use source_damage_isoBrittle, only: &
source_damage_isobrittle_getRateAndItsTangent
use source_damage_isoDuctile, only: &
source_damage_isoductile_getRateAndItsTangent
use source_damage_anisoBrittle, only: &
source_damage_anisobrittle_getRateAndItsTangent
use source_damage_anisoDuctile, only: &
source_damage_anisoductile_getRateAndItsTangent
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
real(pReal), intent(in) :: &
phi
integer :: &
phase, &
grain, &
source, &
constituent
real(pReal) :: &
phiDot, dPhiDot_dPhi, localphiDot, dLocalphiDot_dPhi
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
real(pReal), intent(in) :: &
phi
integer :: &
phase, &
grain, &
source, &
constituent
real(pReal) :: &
phiDot, dPhiDot_dPhi, localphiDot, dLocalphiDot_dPhi
phiDot = 0.0_pReal
dPhiDot_dPhi = 0.0_pReal
do grain = 1, homogenization_Ngrains(material_homogenizationAt(el))
phase = phaseAt(grain,ip,el)
constituent = phasememberAt(grain,ip,el)
do source = 1, phase_Nsources(phase)
select case(phase_source(source,phase))
case (SOURCE_damage_isoBrittle_ID)
call source_damage_isobrittle_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
phiDot = 0.0_pReal
dPhiDot_dPhi = 0.0_pReal
do grain = 1, homogenization_Ngrains(material_homogenizationAt(el))
phase = phaseAt(grain,ip,el)
constituent = phasememberAt(grain,ip,el)
do source = 1, phase_Nsources(phase)
select case(phase_source(source,phase))
case (SOURCE_damage_isoBrittle_ID)
call source_damage_isobrittle_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
case (SOURCE_damage_isoDuctile_ID)
call source_damage_isoductile_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
case (SOURCE_damage_isoDuctile_ID)
call source_damage_isoductile_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
case (SOURCE_damage_anisoBrittle_ID)
call source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
case (SOURCE_damage_anisoBrittle_ID)
call source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
case (SOURCE_damage_anisoDuctile_ID)
call source_damage_anisoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
case (SOURCE_damage_anisoDuctile_ID)
call source_damage_anisoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
case default
localphiDot = 0.0_pReal
dLocalphiDot_dPhi = 0.0_pReal
case default
localphiDot = 0.0_pReal
dLocalphiDot_dPhi = 0.0_pReal
end select
phiDot = phiDot + localphiDot
dPhiDot_dPhi = dPhiDot_dPhi + dLocalphiDot_dPhi
enddo
enddo
phiDot = phiDot/real(homogenization_Ngrains(material_homogenizationAt(el)),pReal)
dPhiDot_dPhi = dPhiDot_dPhi/real(homogenization_Ngrains(material_homogenizationAt(el)),pReal)
end select
phiDot = phiDot + localphiDot
dPhiDot_dPhi = dPhiDot_dPhi + dLocalphiDot_dPhi
enddo
enddo
phiDot = phiDot/real(homogenization_Ngrains(material_homogenizationAt(el)),pReal)
dPhiDot_dPhi = dPhiDot_dPhi/real(homogenization_Ngrains(material_homogenizationAt(el)),pReal)
end subroutine damage_nonlocal_getSourceAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief returns homogenized non local damage diffusion tensor in reference configuration
!--------------------------------------------------------------------------------------------------
function damage_nonlocal_getDiffusion33(ip,el)
use numerics, only: &
charLength
use lattice, only: &
lattice_DamageDiffusion33
use material, only: &
homogenization_Ngrains, &
material_phase, &
material_homogenizationAt
use crystallite, only: &
crystallite_push33ToRef
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
real(pReal), dimension(3,3) :: &
damage_nonlocal_getDiffusion33
integer :: &
homog, &
grain
homog = material_homogenizationAt(el)
damage_nonlocal_getDiffusion33 = 0.0_pReal
do grain = 1, homogenization_Ngrains(homog)
damage_nonlocal_getDiffusion33 = damage_nonlocal_getDiffusion33 + &
crystallite_push33ToRef(grain,ip,el,lattice_DamageDiffusion33(1:3,1:3,material_phase(grain,ip,el)))
enddo
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
real(pReal), dimension(3,3) :: &
damage_nonlocal_getDiffusion33
integer :: &
homog, &
grain
homog = material_homogenizationAt(el)
damage_nonlocal_getDiffusion33 = 0.0_pReal
do grain = 1, homogenization_Ngrains(homog)
damage_nonlocal_getDiffusion33 = damage_nonlocal_getDiffusion33 + &
crystallite_push33ToRef(grain,ip,el,lattice_DamageDiffusion33(1:3,1:3,material_phase(grain,ip,el)))
enddo
damage_nonlocal_getDiffusion33 = &
charLength**2*damage_nonlocal_getDiffusion33/real(homogenization_Ngrains(homog),pReal)
damage_nonlocal_getDiffusion33 = &
charLength**2*damage_nonlocal_getDiffusion33/real(homogenization_Ngrains(homog),pReal)
end function damage_nonlocal_getDiffusion33
!--------------------------------------------------------------------------------------------------
!> @brief Returns homogenized nonlocal damage mobility
!--------------------------------------------------------------------------------------------------
real(pReal) function damage_nonlocal_getMobility(ip,el)
use mesh, only: &
mesh_element
use lattice, only: &
lattice_damageMobility
use material, only: &
material_phase, &
homogenization_Ngrains
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
integer :: &
ipc
damage_nonlocal_getMobility = 0.0_pReal
do ipc = 1, homogenization_Ngrains(mesh_element(3,el))
damage_nonlocal_getMobility = damage_nonlocal_getMobility + lattice_DamageMobility(material_phase(ipc,ip,el))
enddo
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
integer :: &
ipc
damage_nonlocal_getMobility = 0.0_pReal
do ipc = 1, homogenization_Ngrains(mesh_element(3,el))
damage_nonlocal_getMobility = damage_nonlocal_getMobility + lattice_DamageMobility(material_phase(ipc,ip,el))
enddo
damage_nonlocal_getMobility = damage_nonlocal_getMobility/&
real(homogenization_Ngrains(mesh_element(3,el)),pReal)
damage_nonlocal_getMobility = damage_nonlocal_getMobility/&
real(homogenization_Ngrains(mesh_element(3,el)),pReal)
end function damage_nonlocal_getMobility
!--------------------------------------------------------------------------------------------------
!> @brief updated nonlocal damage field with solution from damage phase field PDE
!--------------------------------------------------------------------------------------------------
subroutine damage_nonlocal_putNonLocalDamage(phi,ip,el)
use material, only: &
material_homogenizationAt, &
damageMapping, &
damage
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
real(pReal), intent(in) :: &
phi
integer :: &
homog, &
offset
homog = material_homogenizationAt(el)
offset = damageMapping(homog)%p(ip,el)
damage(homog)%p(offset) = phi
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
real(pReal), intent(in) :: &
phi
integer :: &
homog, &
offset
homog = material_homogenizationAt(el)
offset = damageMapping(homog)%p(ip,el)
damage(homog)%p(offset) = phi
end subroutine damage_nonlocal_putNonLocalDamage
!--------------------------------------------------------------------------------------------------
!> @brief return array of damage results
!--------------------------------------------------------------------------------------------------
function damage_nonlocal_postResults(ip,el)
use material, only: &
material_homogenizationAt, &
damage_typeInstance, &
damageMapping, &
damage
integer, intent(in) :: &
ip, & !< integration point
el !< element
real(pReal), dimension(sum(damage_nonlocal_sizePostResult(:,damage_typeInstance(material_homogenizationAt(el))))) :: &
damage_nonlocal_postResults
integer, intent(in) :: &
ip, & !< integration point
el !< element
real(pReal), dimension(sum(damage_nonlocal_sizePostResult(:,damage_typeInstance(material_homogenizationAt(el))))) :: &
damage_nonlocal_postResults
integer :: &
instance, homog, offset, o, c
homog = material_homogenizationAt(el)
offset = damageMapping(homog)%p(ip,el)
instance = damage_typeInstance(homog)
associate(prm => param(instance))
c = 0
integer :: &
instance, homog, offset, o, c
homog = material_homogenizationAt(el)
offset = damageMapping(homog)%p(ip,el)
instance = damage_typeInstance(homog)
associate(prm => param(instance))
c = 0
outputsLoop: do o = 1,size(prm%outputID)
select case(prm%outputID(o))
case (damage_ID)
damage_nonlocal_postResults(c+1) = damage(homog)%p(offset)
c = c + 1
end select
enddo outputsLoop
outputsLoop: do o = 1,size(prm%outputID)
select case(prm%outputID(o))
case (damage_ID)
damage_nonlocal_postResults(c+1) = damage(homog)%p(offset)
c = c + 1
end select
enddo outputsLoop
end associate
end associate
end function damage_nonlocal_postResults
end module damage_nonlocal

22
src/debug.f90
View File

@ -6,12 +6,12 @@
!> @brief Reading in and interpretating the debugging settings for the various modules
!--------------------------------------------------------------------------------------------------
module debug
use prec, only: &
pInt, &
pReal
use prec
use IO
implicit none
private
integer(pInt), parameter, public :: &
debug_LEVELSELECTIVE = 2_pInt**0_pInt, &
debug_LEVELBASIC = 2_pInt**1_pInt, &
@ -78,19 +78,7 @@ contains
!> @brief reads in parameters from debug.config and allocates arrays
!--------------------------------------------------------------------------------------------------
subroutine debug_init
use prec, only: &
pStringLen
use IO, only: &
IO_read_ASCII, &
IO_error, &
IO_isBlank, &
IO_stringPos, &
IO_stringValue, &
IO_lc, &
IO_floatValue, &
IO_intValue
implicit none
character(len=pStringLen), dimension(:), allocatable :: fileContent
integer :: i, what, j
@ -253,8 +241,6 @@ end subroutine debug_init
!--------------------------------------------------------------------------------------------------
subroutine debug_reset
implicit none
debug_stressMaxLocation = 0_pInt
debug_stressMinLocation = 0_pInt
debug_jacobianMaxLocation = 0_pInt
@ -272,8 +258,6 @@ end subroutine debug_reset
!--------------------------------------------------------------------------------------------------
subroutine debug_info
implicit none
!$OMP CRITICAL (write2out)
debugOutputCPFEM: if (iand(debug_level(debug_CPFEM),debug_LEVELBASIC) /= 0 &
.and. any(debug_stressMinLocation /= 0_pInt) &

24
src/homogenization.f90
View File

@ -16,6 +16,12 @@ module homogenization
use crystallite
use mesh
use FEsolving
use thermal_isothermal
use thermal_adiabatic
use thermal_conduction
use damage_none
use damage_local
use damage_nonlocal
#if defined(PETSc) || defined(DAMASK_HDF5)
use results
use HDF5_utilities
@ -131,12 +137,6 @@ contains
!> @brief module initialization
!--------------------------------------------------------------------------------------------------
subroutine homogenization_init
use thermal_isothermal
use thermal_adiabatic
use thermal_conduction
use damage_none
use damage_local
use damage_nonlocal
integer, parameter :: FILEUNIT = 200
integer :: e,i,p
@ -668,10 +668,6 @@ end subroutine partitionDeformation
!> "happy" with result
!--------------------------------------------------------------------------------------------------
function updateState(ip,el)
use thermal_adiabatic, only: &
thermal_adiabatic_updateState
use damage_local, only: &
damage_local_updateState
integer, intent(in) :: &
ip, & !< integration point
@ -753,14 +749,6 @@ end subroutine averageStressAndItsTangent
!> if homogenization_sizePostResults(i,e) > 0 !!
!--------------------------------------------------------------------------------------------------
function postResults(ip,el)
use thermal_adiabatic, only: &
thermal_adiabatic_postResults
use thermal_conduction, only: &
thermal_conduction_postResults
use damage_local, only: &
damage_local_postResults
use damage_nonlocal, only: &
damage_nonlocal_postResults
integer, intent(in) :: &
ip, & !< integration point

335
src/kinematics_cleavage_opening.f90
View File

@ -5,44 +5,51 @@
!> @details to be done
!--------------------------------------------------------------------------------------------------
module kinematics_cleavage_opening
use prec
use prec
use IO
use config
use debug
use math
use lattice
use material
implicit none
private
integer, dimension(:), allocatable, private :: kinematics_cleavage_opening_instance
implicit none
private
type, private :: tParameters !< container type for internal constitutive parameters
integer :: &
totalNcleavage
integer, dimension(:), allocatable :: &
Ncleavage !< active number of cleavage systems per family
real(pReal) :: &
sdot0, &
n
real(pReal), dimension(:), allocatable :: &
critDisp, &
critLoad
end type
integer, dimension(:), allocatable :: kinematics_cleavage_opening_instance
type :: tParameters !< container type for internal constitutive parameters
integer :: &
totalNcleavage
integer, dimension(:), allocatable :: &
Ncleavage !< active number of cleavage systems per family
real(pReal) :: &
sdot0, &
n
real(pReal), dimension(:), allocatable :: &
critDisp, &
critLoad
end type
! Begin Deprecated
integer, dimension(:), allocatable, private :: &
kinematics_cleavage_opening_totalNcleavage !< total number of cleavage systems
integer, dimension(:,:), allocatable, private :: &
kinematics_cleavage_opening_Ncleavage !< number of cleavage systems per family
real(pReal), dimension(:), allocatable, private :: &
kinematics_cleavage_opening_sdot_0, &
kinematics_cleavage_opening_N
integer, dimension(:), allocatable :: &
kinematics_cleavage_opening_totalNcleavage !< total number of cleavage systems
integer, dimension(:,:), allocatable :: &
kinematics_cleavage_opening_Ncleavage !< number of cleavage systems per family
real(pReal), dimension(:), allocatable :: &
kinematics_cleavage_opening_sdot_0, &
kinematics_cleavage_opening_N
real(pReal), dimension(:,:), allocatable, private :: &
kinematics_cleavage_opening_critDisp, &
kinematics_cleavage_opening_critLoad
real(pReal), dimension(:,:), allocatable :: &
kinematics_cleavage_opening_critDisp, &
kinematics_cleavage_opening_critLoad
! End Deprecated
public :: &
kinematics_cleavage_opening_init, &
kinematics_cleavage_opening_LiAndItsTangent
public :: &
kinematics_cleavage_opening_init, &
kinematics_cleavage_opening_LiAndItsTangent
contains
@ -51,170 +58,144 @@ contains
!> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
subroutine kinematics_cleavage_opening_init()
use debug, only: &
debug_level,&
debug_constitutive,&
debug_levelBasic
use config, only: &
config_phase
use IO, only: &
IO_error
use material, only: &
phase_kinematics, &
KINEMATICS_cleavage_opening_label, &
KINEMATICS_cleavage_opening_ID
use lattice, only: &
lattice_maxNcleavageFamily, &
lattice_NcleavageSystem
subroutine kinematics_cleavage_opening_init
integer, allocatable, dimension(:) :: tempInt
real(pReal), allocatable, dimension(:) :: tempFloat
integer, allocatable, dimension(:) :: tempInt
real(pReal), allocatable, dimension(:) :: tempFloat
integer :: maxNinstance,p,instance,kinematics
integer :: maxNinstance,p,instance
write(6,'(/,a)') ' <<<+- kinematics_'//KINEMATICS_cleavage_opening_LABEL//' init -+>>>'
write(6,'(/,a)') ' <<<+- kinematics_'//KINEMATICS_cleavage_opening_LABEL//' init -+>>>'
maxNinstance = int(count(phase_kinematics == KINEMATICS_cleavage_opening_ID))
if (maxNinstance == 0) return
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) &
write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
allocate(kinematics_cleavage_opening_instance(size(config_phase)), source=0)
do p = 1, size(config_phase)
kinematics_cleavage_opening_instance(p) = count(phase_kinematics(:,1:p) == kinematics_cleavage_opening_ID) ! ToDo: count correct?
enddo
allocate(kinematics_cleavage_opening_critDisp(lattice_maxNcleavageFamily,maxNinstance), source=0.0_pReal)
allocate(kinematics_cleavage_opening_critLoad(lattice_maxNcleavageFamily,maxNinstance), source=0.0_pReal)
allocate(kinematics_cleavage_opening_Ncleavage(lattice_maxNcleavageFamily,maxNinstance), source=0)
allocate(kinematics_cleavage_opening_totalNcleavage(maxNinstance), source=0)
allocate(kinematics_cleavage_opening_sdot_0(maxNinstance), source=0.0_pReal)
allocate(kinematics_cleavage_opening_N(maxNinstance), source=0.0_pReal)
maxNinstance = count(phase_kinematics == KINEMATICS_cleavage_opening_ID)
if (maxNinstance == 0) return
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) &
write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
allocate(kinematics_cleavage_opening_instance(size(config_phase)), source=0)
do p = 1, size(config_phase)
kinematics_cleavage_opening_instance(p) = count(phase_kinematics(:,1:p) == kinematics_cleavage_opening_ID) ! ToDo: count correct?
enddo
allocate(kinematics_cleavage_opening_critDisp(lattice_maxNcleavageFamily,maxNinstance), source=0.0_pReal)
allocate(kinematics_cleavage_opening_critLoad(lattice_maxNcleavageFamily,maxNinstance), source=0.0_pReal)
allocate(kinematics_cleavage_opening_Ncleavage(lattice_maxNcleavageFamily,maxNinstance), source=0)
allocate(kinematics_cleavage_opening_totalNcleavage(maxNinstance), source=0)
allocate(kinematics_cleavage_opening_sdot_0(maxNinstance), source=0.0_pReal)
allocate(kinematics_cleavage_opening_N(maxNinstance), source=0.0_pReal)
do p = 1, size(config_phase)
if (all(phase_kinematics(:,p) /= KINEMATICS_cleavage_opening_ID)) cycle
instance = kinematics_cleavage_opening_instance(p)
kinematics_cleavage_opening_sdot_0(instance) = config_phase(p)%getFloat('anisobrittle_sdot0')
kinematics_cleavage_opening_N(instance) = config_phase(p)%getFloat('anisobrittle_ratesensitivity')
tempInt = config_phase(p)%getInts('ncleavage')
kinematics_cleavage_opening_Ncleavage(1:size(tempInt),instance) = tempInt
do p = 1, size(config_phase)
if (all(phase_kinematics(:,p) /= KINEMATICS_cleavage_opening_ID)) cycle
instance = kinematics_cleavage_opening_instance(p)
kinematics_cleavage_opening_sdot_0(instance) = config_phase(p)%getFloat('anisobrittle_sdot0')
kinematics_cleavage_opening_N(instance) = config_phase(p)%getFloat('anisobrittle_ratesensitivity')
tempInt = config_phase(p)%getInts('ncleavage')
kinematics_cleavage_opening_Ncleavage(1:size(tempInt),instance) = tempInt
tempFloat = config_phase(p)%getFloats('anisobrittle_criticaldisplacement',requiredSize=size(tempInt))
kinematics_cleavage_opening_critDisp(1:size(tempInt),instance) = tempFloat
tempFloat = config_phase(p)%getFloats('anisobrittle_criticaldisplacement',requiredSize=size(tempInt))
kinematics_cleavage_opening_critDisp(1:size(tempInt),instance) = tempFloat
tempFloat = config_phase(p)%getFloats('anisobrittle_criticalload',requiredSize=size(tempInt))
kinematics_cleavage_opening_critLoad(1:size(tempInt),instance) = tempFloat
tempFloat = config_phase(p)%getFloats('anisobrittle_criticalload',requiredSize=size(tempInt))
kinematics_cleavage_opening_critLoad(1:size(tempInt),instance) = tempFloat
kinematics_cleavage_opening_Ncleavage(1:lattice_maxNcleavageFamily,instance) = &
min(lattice_NcleavageSystem(1:lattice_maxNcleavageFamily,p),& ! limit active cleavage systems per family to min of available and requested
kinematics_cleavage_opening_Ncleavage(1:lattice_maxNcleavageFamily,instance))
kinematics_cleavage_opening_totalNcleavage(instance) = sum(kinematics_cleavage_opening_Ncleavage(:,instance)) ! how many cleavage systems altogether
if (kinematics_cleavage_opening_sdot_0(instance) <= 0.0_pReal) &
call IO_error(211,el=instance,ext_msg='sdot_0 ('//KINEMATICS_cleavage_opening_LABEL//')')
if (any(kinematics_cleavage_opening_critDisp(1:size(tempInt),instance) < 0.0_pReal)) &
call IO_error(211,el=instance,ext_msg='critical_displacement ('//KINEMATICS_cleavage_opening_LABEL//')')
if (any(kinematics_cleavage_opening_critLoad(1:size(tempInt),instance) < 0.0_pReal)) &
call IO_error(211,el=instance,ext_msg='critical_load ('//KINEMATICS_cleavage_opening_LABEL//')')
if (kinematics_cleavage_opening_N(instance) <= 0.0_pReal) &
call IO_error(211,el=instance,ext_msg='rate_sensitivity ('//KINEMATICS_cleavage_opening_LABEL//')')
enddo
kinematics_cleavage_opening_Ncleavage(1:lattice_maxNcleavageFamily,instance) = &
min(lattice_NcleavageSystem(1:lattice_maxNcleavageFamily,p),& ! limit active cleavage systems per family to min of available and requested
kinematics_cleavage_opening_Ncleavage(1:lattice_maxNcleavageFamily,instance))
kinematics_cleavage_opening_totalNcleavage(instance) = sum(kinematics_cleavage_opening_Ncleavage(:,instance)) ! how many cleavage systems altogether
if (kinematics_cleavage_opening_sdot_0(instance) <= 0.0_pReal) &
call IO_error(211,el=instance,ext_msg='sdot_0 ('//KINEMATICS_cleavage_opening_LABEL//')')
if (any(kinematics_cleavage_opening_critDisp(1:size(tempInt),instance) < 0.0_pReal)) &
call IO_error(211,el=instance,ext_msg='critical_displacement ('//KINEMATICS_cleavage_opening_LABEL//')')
if (any(kinematics_cleavage_opening_critLoad(1:size(tempInt),instance) < 0.0_pReal)) &
call IO_error(211,el=instance,ext_msg='critical_load ('//KINEMATICS_cleavage_opening_LABEL//')')
if (kinematics_cleavage_opening_N(instance) <= 0.0_pReal) &
call IO_error(211,el=instance,ext_msg='rate_sensitivity ('//KINEMATICS_cleavage_opening_LABEL//')')
enddo
end subroutine kinematics_cleavage_opening_init
!--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the velocity gradient
!--------------------------------------------------------------------------------------------------
subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, ip, el)
use math, only: &
math_mul33xx33
use material, only: &
material_phase, &
material_homogenizationAt, &
damage, &
damageMapping
use lattice, only: &
lattice_Scleavage, &
lattice_maxNcleavageFamily, &
lattice_NcleavageSystem
integer, intent(in) :: &
ipc, & !< grain number
ip, & !< integration point number
el !< element number
real(pReal), intent(in), dimension(3,3) :: &
S
real(pReal), intent(out), dimension(3,3) :: &
Ld !< damage velocity gradient
real(pReal), intent(out), dimension(3,3,3,3) :: &
dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor)
integer :: &
instance, phase, &
homog, damageOffset, &
f, i, index_myFamily, k, l, m, n
real(pReal) :: &
traction_d, traction_t, traction_n, traction_crit, &
udotd, dudotd_dt, udott, dudott_dt, udotn, dudotn_dt
integer, intent(in) :: &
ipc, & !< grain number
ip, & !< integration point number
el !< element number
real(pReal), intent(in), dimension(3,3) :: &
S
real(pReal), intent(out), dimension(3,3) :: &
Ld !< damage velocity gradient
real(pReal), intent(out), dimension(3,3,3,3) :: &
dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor)
integer :: &
instance, phase, &
homog, damageOffset, &
f, i, index_myFamily, k, l, m, n
real(pReal) :: &
traction_d, traction_t, traction_n, traction_crit, &
udotd, dudotd_dt, udott, dudott_dt, udotn, dudotn_dt
phase = material_phase(ipc,ip,el)
instance = kinematics_cleavage_opening_instance(phase)
homog = material_homogenizationAt(el)
damageOffset = damageMapping(homog)%p(ip,el)
Ld = 0.0_pReal
dLd_dTstar = 0.0_pReal
do f = 1,lattice_maxNcleavageFamily
index_myFamily = sum(lattice_NcleavageSystem(1:f-1,phase)) ! at which index starts my family
do i = 1,kinematics_cleavage_opening_Ncleavage(f,instance) ! process each (active) cleavage system in family
traction_d = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,1,index_myFamily+i,phase))
traction_t = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,2,index_myFamily+i,phase))
traction_n = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,3,index_myFamily+i,phase))
traction_crit = kinematics_cleavage_opening_critLoad(f,instance)* &
damage(homog)%p(damageOffset)*damage(homog)%p(damageOffset)
udotd = &
sign(1.0_pReal,traction_d)* &
kinematics_cleavage_opening_sdot_0(instance)* &
(max(0.0_pReal, abs(traction_d) - traction_crit)/traction_crit)**kinematics_cleavage_opening_N(instance)
if (abs(udotd) > tol_math_check) then
Ld = Ld + udotd*lattice_Scleavage(1:3,1:3,1,index_myFamily+i,phase)
dudotd_dt = sign(1.0_pReal,traction_d)*udotd*kinematics_cleavage_opening_N(instance)/ &
max(0.0_pReal, abs(traction_d) - traction_crit)
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + &
dudotd_dt*lattice_Scleavage(k,l,1,index_myFamily+i,phase)* &
lattice_Scleavage(m,n,1,index_myFamily+i,phase)
endif
phase = material_phase(ipc,ip,el)
instance = kinematics_cleavage_opening_instance(phase)
homog = material_homogenizationAt(el)
damageOffset = damageMapping(homog)%p(ip,el)
Ld = 0.0_pReal
dLd_dTstar = 0.0_pReal
do f = 1,lattice_maxNcleavageFamily
index_myFamily = sum(lattice_NcleavageSystem(1:f-1,phase)) ! at which index starts my family
do i = 1,kinematics_cleavage_opening_Ncleavage(f,instance) ! process each (active) cleavage system in family
traction_d = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,1,index_myFamily+i,phase))
traction_t = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,2,index_myFamily+i,phase))
traction_n = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,3,index_myFamily+i,phase))
traction_crit = kinematics_cleavage_opening_critLoad(f,instance)* &
damage(homog)%p(damageOffset)*damage(homog)%p(damageOffset)
udotd = &
sign(1.0_pReal,traction_d)* &
kinematics_cleavage_opening_sdot_0(instance)* &
(max(0.0_pReal, abs(traction_d) - traction_crit)/traction_crit)**kinematics_cleavage_opening_N(instance)
if (abs(udotd) > tol_math_check) then
Ld = Ld + udotd*lattice_Scleavage(1:3,1:3,1,index_myFamily+i,phase)
dudotd_dt = sign(1.0_pReal,traction_d)*udotd*kinematics_cleavage_opening_N(instance)/ &
max(0.0_pReal, abs(traction_d) - traction_crit)
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + &
dudotd_dt*lattice_Scleavage(k,l,1,index_myFamily+i,phase)* &
lattice_Scleavage(m,n,1,index_myFamily+i,phase)
endif
udott = &
sign(1.0_pReal,traction_t)* &
kinematics_cleavage_opening_sdot_0(instance)* &
(max(0.0_pReal, abs(traction_t) - traction_crit)/traction_crit)**kinematics_cleavage_opening_N(instance)
if (abs(udott) > tol_math_check) then
Ld = Ld + udott*lattice_Scleavage(1:3,1:3,2,index_myFamily+i,phase)
dudott_dt = sign(1.0_pReal,traction_t)*udott*kinematics_cleavage_opening_N(instance)/ &
max(0.0_pReal, abs(traction_t) - traction_crit)
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + &
dudott_dt*lattice_Scleavage(k,l,2,index_myFamily+i,phase)* &
lattice_Scleavage(m,n,2,index_myFamily+i,phase)
endif
udott = &
sign(1.0_pReal,traction_t)* &
kinematics_cleavage_opening_sdot_0(instance)* &
(max(0.0_pReal, abs(traction_t) - traction_crit)/traction_crit)**kinematics_cleavage_opening_N(instance)
if (abs(udott) > tol_math_check) then
Ld = Ld + udott*lattice_Scleavage(1:3,1:3,2,index_myFamily+i,phase)
dudott_dt = sign(1.0_pReal,traction_t)*udott*kinematics_cleavage_opening_N(instance)/ &
max(0.0_pReal, abs(traction_t) - traction_crit)
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + &
dudott_dt*lattice_Scleavage(k,l,2,index_myFamily+i,phase)* &
lattice_Scleavage(m,n,2,index_myFamily+i,phase)
endif
udotn = &
sign(1.0_pReal,traction_n)* &
kinematics_cleavage_opening_sdot_0(instance)* &
(max(0.0_pReal, abs(traction_n) - traction_crit)/traction_crit)**kinematics_cleavage_opening_N(instance)
if (abs(udotn) > tol_math_check) then
Ld = Ld + udotn*lattice_Scleavage(1:3,1:3,3,index_myFamily+i,phase)
dudotn_dt = sign(1.0_pReal,traction_n)*udotn*kinematics_cleavage_opening_N(instance)/ &
max(0.0_pReal, abs(traction_n) - traction_crit)
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + &
dudotn_dt*lattice_Scleavage(k,l,3,index_myFamily+i,phase)* &
lattice_Scleavage(m,n,3,index_myFamily+i,phase)
endif
enddo
enddo
udotn = &
sign(1.0_pReal,traction_n)* &
kinematics_cleavage_opening_sdot_0(instance)* &
(max(0.0_pReal, abs(traction_n) - traction_crit)/traction_crit)**kinematics_cleavage_opening_N(instance)
if (abs(udotn) > tol_math_check) then
Ld = Ld + udotn*lattice_Scleavage(1:3,1:3,3,index_myFamily+i,phase)
dudotn_dt = sign(1.0_pReal,traction_n)*udotn*kinematics_cleavage_opening_N(instance)/ &
max(0.0_pReal, abs(traction_n) - traction_crit)
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + &
dudotn_dt*lattice_Scleavage(k,l,3,index_myFamily+i,phase)* &
lattice_Scleavage(m,n,3,index_myFamily+i,phase)
endif
enddo
enddo
end subroutine kinematics_cleavage_opening_LiAndItsTangent

44
src/kinematics_slipplane_opening.f90
View File

@ -6,12 +6,19 @@
!--------------------------------------------------------------------------------------------------
module kinematics_slipplane_opening
use prec
use config
use IO
use debug
use math
use lattice
use material
implicit none
private
integer, dimension(:), allocatable, private :: kinematics_slipplane_opening_instance
integer, dimension(:), allocatable :: kinematics_slipplane_opening_instance
type, private :: tParameters !< container type for internal constitutive parameters
type :: tParameters !< container type for internal constitutive parameters
integer :: &
totalNslip
integer, dimension(:), allocatable :: &
@ -19,7 +26,7 @@ module kinematics_slipplane_opening
real(pReal) :: &
sdot0, &
n
real(pReal), dimension(:), allocatable :: &
real(pReal), dimension(:), allocatable :: &
critLoad
real(pReal), dimension(:,:), allocatable :: &
slip_direction, &
@ -27,7 +34,8 @@ module kinematics_slipplane_opening
slip_transverse
end type tParameters
type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance)
type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance)
public :: &
kinematics_slipplane_opening_init, &
kinematics_slipplane_opening_LiAndItsTangent
@ -39,25 +47,9 @@ contains
!> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
subroutine kinematics_slipplane_opening_init()
use debug, only: &
debug_level,&
debug_constitutive,&
debug_levelBasic
use config, only: &
config_phase
use IO, only: &
IO_error
use math, only: &
math_expand
use material, only: &
phase_kinematics, &
KINEMATICS_slipplane_opening_label, &
KINEMATICS_slipplane_opening_ID
use lattice
subroutine kinematics_slipplane_opening_init
integer :: maxNinstance,p,instance,kinematics
integer :: maxNinstance,p,instance
write(6,'(/,a)') ' <<<+- kinematics_'//KINEMATICS_slipplane_opening_LABEL//' init -+>>>'
@ -111,14 +103,6 @@ end subroutine kinematics_slipplane_opening_init
!> @brief contains the constitutive equation for calculating the velocity gradient
!--------------------------------------------------------------------------------------------------
subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, ip, el)
use math, only: &
math_mul33xx33, &
math_outer
use material, only: &
material_phase, &
material_homogenizationAt, &
damage, &
damageMapping
integer, intent(in) :: &
ipc, & !< grain number

50
src/kinematics_thermal_expansion.f90
View File

@ -5,11 +5,17 @@
!--------------------------------------------------------------------------------------------------
module kinematics_thermal_expansion
use prec
use IO
use config
use debug
use math
use lattice
use material
implicit none
private
type, private :: tParameters
type :: tParameters
real(pReal), allocatable, dimension(:,:,:) :: &
expansion
end type tParameters
@ -28,19 +34,9 @@ contains
!> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
subroutine kinematics_thermal_expansion_init()
use debug, only: &
debug_level,&
debug_constitutive,&
debug_levelBasic
use material, only: &
phase_kinematics, &
KINEMATICS_thermal_expansion_label, &
KINEMATICS_thermal_expansion_ID
use config, only: &
config_phase
subroutine kinematics_thermal_expansion_init
integer(pInt) :: &
integer :: &
Ninstance, &
p, i
real(pReal), dimension(:), allocatable :: &
@ -48,14 +44,14 @@ subroutine kinematics_thermal_expansion_init()
write(6,'(/,a)') ' <<<+- kinematics_'//KINEMATICS_thermal_expansion_LABEL//' init -+>>>'
Ninstance = int(count(phase_kinematics == KINEMATICS_thermal_expansion_ID),pInt)
Ninstance = count(phase_kinematics == KINEMATICS_thermal_expansion_ID)
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) &
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
allocate(param(Ninstance))
do p = 1_pInt, size(phase_kinematics)
do p = 1, size(phase_kinematics)
if (all(phase_kinematics(:,p) /= KINEMATICS_thermal_expansion_ID)) cycle
! ToDo: Here we need to decide how to extend the concept of instances to
@ -78,13 +74,8 @@ end subroutine kinematics_thermal_expansion_init
!> @brief report initial thermal strain based on current temperature deviation from reference
!--------------------------------------------------------------------------------------------------
pure function kinematics_thermal_expansion_initialStrain(homog,phase,offset)
use material, only: &
temperature
use lattice, only: &
lattice_thermalExpansion33, &
lattice_referenceTemperature
integer(pInt), intent(in) :: &
integer, intent(in) :: &
phase, &
homog, offset
real(pReal), dimension(3,3) :: &
@ -106,17 +97,8 @@ end function kinematics_thermal_expansion_initialStrain
!> @brief contains the constitutive equation for calculating the velocity gradient
!--------------------------------------------------------------------------------------------------
subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, ipc, ip, el)
use material, only: &
material_phase, &
material_homogenizationAt, &
temperature, &
temperatureRate, &
thermalMapping
use lattice, only: &
lattice_thermalExpansion33, &
lattice_referenceTemperature
integer(pInt), intent(in) :: &
integer, intent(in) :: &
ipc, & !< grain number
ip, & !< integration point number
el !< element number
@ -124,7 +106,7 @@ subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, ipc, ip,
Li !< thermal velocity gradient
real(pReal), intent(out), dimension(3,3,3,3) :: &
dLi_dTstar !< derivative of Li with respect to Tstar (4th-order tensor defined to be zero)
integer(pInt) :: &
integer :: &
phase, &
homog, offset
real(pReal) :: &

208
src/lattice.f90
View File

@ -7,8 +7,10 @@
! and cleavage as well as interaction among the various systems
!--------------------------------------------------------------------------------------------------
module lattice
use prec, only: &
pReal
use prec
use IO
use config
use math
use future
implicit none
@ -28,25 +30,25 @@ module lattice
!--------------------------------------------------------------------------------------------------
! face centered cubic
integer, dimension(2), parameter, private :: &
integer, dimension(2), parameter :: &
LATTICE_FCC_NSLIPSYSTEM = [12, 6] !< # of slip systems per family for fcc
integer, dimension(1), parameter, private :: &
integer, dimension(1), parameter :: &
LATTICE_FCC_NTWINSYSTEM = [12] !< # of twin systems per family for fcc
integer, dimension(1), parameter, private :: &
integer, dimension(1), parameter :: &
LATTICE_FCC_NTRANSSYSTEM = [12] !< # of transformation systems per family for fcc
integer, dimension(2), parameter, private :: &
integer, dimension(2), parameter :: &
LATTICE_FCC_NCLEAVAGESYSTEM = [3, 4] !< # of cleavage systems per family for fcc
integer, parameter, private :: &
integer, parameter :: &
LATTICE_FCC_NSLIP = sum(LATTICE_FCC_NSLIPSYSTEM), & !< total # of slip systems for fcc
LATTICE_FCC_NTWIN = sum(LATTICE_FCC_NTWINSYSTEM), & !< total # of twin systems for fcc
LATTICE_FCC_NTRANS = sum(LATTICE_FCC_NTRANSSYSTEM), & !< total # of transformation systems for fcc
LATTICE_FCC_NCLEAVAGE = sum(LATTICE_FCC_NCLEAVAGESYSTEM) !< total # of cleavage systems for fcc
real(pReal), dimension(3+3,LATTICE_FCC_NSLIP), parameter, private :: &
real(pReal), dimension(3+3,LATTICE_FCC_NSLIP), parameter :: &
LATTICE_FCC_SYSTEMSLIP = reshape(real([&
! Slip direction Plane normal ! SCHMID-BOAS notation
0, 1,-1, 1, 1, 1, & ! B2
@ -70,11 +72,11 @@ module lattice
0, 1,-1, 0, 1, 1 &
],pReal),shape(LATTICE_FCC_SYSTEMSLIP)) !< Slip system <110>{111} directions. Sorted according to Eisenlohr & Hantcherli
character(len=*), dimension(2), parameter, private :: LATTICE_FCC_SLIPFAMILY_NAME = &
character(len=*), dimension(2), parameter :: LATTICE_FCC_SLIPFAMILY_NAME = &
['<0 1 -1>{1 1 1}', &
'<0 1 -1>{0 1 1}']
real(pReal), dimension(3+3,LATTICE_FCC_NTWIN), parameter, private :: &
real(pReal), dimension(3+3,LATTICE_FCC_NTWIN), parameter :: &
LATTICE_FCC_SYSTEMTWIN = reshape(real( [&
-2, 1, 1, 1, 1, 1, &
1,-2, 1, 1, 1, 1, &
@ -90,7 +92,7 @@ module lattice
-1, 1, 2, -1, 1,-1 &
],pReal),shape(LATTICE_FCC_SYSTEMTWIN)) !< Twin system <112>{111} directions. Sorted according to Eisenlohr & Hantcherli
character(len=*), dimension(1), parameter, private :: LATTICE_FCC_TWINFAMILY_NAME = &
character(len=*), dimension(1), parameter :: LATTICE_FCC_TWINFAMILY_NAME = &
['<-2 1 1>{1 1 1}']
@ -110,7 +112,7 @@ module lattice
10,11 &
],shape(LATTICE_FCC_TWINNUCLEATIONSLIPPAIR))
real(pReal), dimension(3+3,LATTICE_FCC_NCLEAVAGE), parameter, private :: &
real(pReal), dimension(3+3,LATTICE_FCC_NCLEAVAGE), parameter :: &
LATTICE_FCC_SYSTEMCLEAVAGE = reshape(real([&
! Cleavage direction Plane normal
0, 1, 0, 1, 0, 0, &
@ -124,21 +126,21 @@ module lattice
!--------------------------------------------------------------------------------------------------
! body centered cubic
integer, dimension(2), parameter, private :: &
integer, dimension(2), parameter :: &
LATTICE_BCC_NSLIPSYSTEM = [12, 12] !< # of slip systems per family for bcc
integer, dimension(1), parameter, private :: &
integer, dimension(1), parameter :: &
LATTICE_BCC_NTWINSYSTEM = [12] !< # of twin systems per family for bcc
integer, dimension(2), parameter, private :: &
integer, dimension(2), parameter :: &
LATTICE_BCC_NCLEAVAGESYSTEM = [3, 6] !< # of cleavage systems per family for bcc
integer, parameter, private :: &
integer, parameter :: &
LATTICE_BCC_NSLIP = sum(LATTICE_BCC_NSLIPSYSTEM), & !< total # of slip systems for bcc
LATTICE_BCC_NTWIN = sum(LATTICE_BCC_NTWINSYSTEM), & !< total # of twin systems for bcc
LATTICE_BCC_NCLEAVAGE = sum(LATTICE_BCC_NCLEAVAGESYSTEM) !< total # of cleavage systems for bcc
real(pReal), dimension(3+3,LATTICE_BCC_NSLIP), parameter, private :: &
real(pReal), dimension(3+3,LATTICE_BCC_NSLIP), parameter :: &
LATTICE_BCC_SYSTEMSLIP = reshape(real([&
! Slip direction Plane normal
! Slip system <111>{110}
@ -169,11 +171,11 @@ module lattice
1, 1, 1, 1, 1,-2 &
],pReal),shape(LATTICE_BCC_SYSTEMSLIP))
character(len=*), dimension(2), parameter, private :: LATTICE_BCC_SLIPFAMILY_NAME = &
character(len=*), dimension(2), parameter :: LATTICE_BCC_SLIPFAMILY_NAME = &
['<1 -1 1>{0 1 1}', &
'<1 -1 1>{2 1 1}']
real(pReal), dimension(3+3,LATTICE_BCC_NTWIN), parameter, private :: &
real(pReal), dimension(3+3,LATTICE_BCC_NTWIN), parameter :: &
LATTICE_BCC_SYSTEMTWIN = reshape(real([&
! Twin system <111>{112}
-1, 1, 1, 2, 1, 1, &
@ -190,10 +192,10 @@ module lattice
1, 1, 1, 1, 1,-2 &
],pReal),shape(LATTICE_BCC_SYSTEMTWIN))
character(len=*), dimension(1), parameter, private :: LATTICE_BCC_TWINFAMILY_NAME = &
character(len=*), dimension(1), parameter :: LATTICE_BCC_TWINFAMILY_NAME = &
['<1 1 1>{2 1 1}']
real(pReal), dimension(3+3,LATTICE_BCC_NCLEAVAGE), parameter, private :: &
real(pReal), dimension(3+3,LATTICE_BCC_NCLEAVAGE), parameter :: &
LATTICE_BCC_SYSTEMCLEAVAGE = reshape(real([&
! Cleavage direction Plane normal
0, 1, 0, 1, 0, 0, &
@ -209,21 +211,21 @@ module lattice
!--------------------------------------------------------------------------------------------------
! hexagonal
integer, dimension(6), parameter, private :: &
integer, dimension(6), parameter :: &
LATTICE_HEX_NSLIPSYSTEM = [3, 3, 3, 6, 12, 6] !< # of slip systems per family for hex
integer, dimension(4), parameter, private :: &
integer, dimension(4), parameter :: &
LATTICE_HEX_NTWINSYSTEM = [6, 6, 6, 6] !< # of slip systems per family for hex
integer, dimension(1), parameter, private :: &
integer, dimension(1), parameter :: &
LATTICE_HEX_NCLEAVAGESYSTEM = [3] !< # of cleavage systems per family for hex
integer, parameter, private :: &
integer, parameter :: &
LATTICE_HEX_NSLIP = sum(LATTICE_HEX_NSLIPSYSTEM), & !< total # of slip systems for hex
LATTICE_HEX_NTWIN = sum(LATTICE_HEX_NTWINSYSTEM), & !< total # of twin systems for hex
LATTICE_HEX_NCLEAVAGE = sum(LATTICE_HEX_NCLEAVAGESYSTEM) !< total # of cleavage systems for hex
real(pReal), dimension(4+4,LATTICE_HEX_NSLIP), parameter, private :: &
real(pReal), dimension(4+4,LATTICE_HEX_NSLIP), parameter :: &
LATTICE_HEX_SYSTEMSLIP = reshape(real([&
! Slip direction Plane normal
! Basal systems <11.0>{00.1} (independent of c/a-ratio, Bravais notation (4 coordinate base))
@ -267,7 +269,7 @@ module lattice
1, 1, -2, 3, -1, -1, 2, 2 &
],pReal),shape(LATTICE_HEX_SYSTEMSLIP)) !< slip systems for hex sorted by A. Alankar & P. Eisenlohr
character(len=*), dimension(6), parameter, private :: LATTICE_HEX_SLIPFAMILY_NAME = &
character(len=*), dimension(6), parameter :: LATTICE_HEX_SLIPFAMILY_NAME = &
['<1 1 . 1>{0 0 . 1} ', &
'<1 1 . 1>{1 0 . 0} ', &
'<1 0 . 0>{1 1 . 0} ', &
@ -275,7 +277,7 @@ module lattice
'<1 1 . 3>{-1 0 . 1} ', &
'<1 1 . 3>{-1 -1 . 2}']
real(pReal), dimension(4+4,LATTICE_HEX_NTWIN), parameter, private :: &
real(pReal), dimension(4+4,LATTICE_HEX_NTWIN), parameter :: &
LATTICE_HEX_SYSTEMTWIN = reshape(real([&
! Compression or Tension =f(twinning shear=f(c/a)) for each metal ! (according to Yoo 1981)
1, -1, 0, 1, -1, 1, 0, 2, & ! <-10.1>{10.2} shear = (3-(c/a)^2)/(sqrt(3) c/a)
@ -307,13 +309,13 @@ module lattice
1, 1, -2, -3, 1, 1, -2, 2 &
],pReal),shape(LATTICE_HEX_SYSTEMTWIN)) !< twin systems for hex, order follows Prof. Tom Bieler's scheme
character(len=*), dimension(4), parameter, private :: LATTICE_HEX_TWINFAMILY_NAME = &
character(len=*), dimension(4), parameter :: LATTICE_HEX_TWINFAMILY_NAME = &
['<-1 0 . 1>{1 0 . 2} ', &
'<1 1 . 6>{-1 -1 . 1}', &
'<1 0 . -2>{1 0 . 1} ', &
'<1 1 . -3>{1 1 . 2} ']
real(pReal), dimension(4+4,LATTICE_HEX_NCLEAVAGE), parameter, private :: &
real(pReal), dimension(4+4,LATTICE_HEX_NCLEAVAGE), parameter :: &
LATTICE_HEX_SYSTEMCLEAVAGE = reshape(real([&
! Cleavage direction Plane normal
2,-1,-1, 0, 0, 0, 0, 1, &
@ -324,13 +326,13 @@ module lattice
!--------------------------------------------------------------------------------------------------
! body centered tetragonal
integer, dimension(13), parameter, private :: &
integer, dimension(13), parameter :: &
LATTICE_BCT_NSLIPSYSTEM = [2, 2, 2, 4, 2, 4, 2, 2, 4, 8, 4, 8, 8 ] !< # of slip systems per family for bct (Sn) Bieler J. Electr Mater 2009
integer, parameter, private :: &
integer, parameter :: &
LATTICE_BCT_NSLIP = sum(LATTICE_BCT_NSLIPSYSTEM) !< total # of slip systems for bct
real(pReal), dimension(3+3,LATTICE_BCT_NSLIP), parameter, private :: &
real(pReal), dimension(3+3,LATTICE_BCT_NSLIP), parameter :: &
LATTICE_BCT_SYSTEMSLIP = reshape(real([&
! Slip direction Plane normal
! Slip family 1 {100)<001] (Bravais notation {hkl)<uvw] for bct c/a = 0.5456)
@ -400,7 +402,7 @@ module lattice
1, 1, 1, 1,-2, 1 &
],pReal),[ 3 + 3,LATTICE_BCT_NSLIP]) !< slip systems for bct sorted by Bieler
character(len=*), dimension(13), parameter, private :: LATTICE_BCT_SLIPFAMILY_NAME = &
character(len=*), dimension(13), parameter :: LATTICE_BCT_SLIPFAMILY_NAME = &
['{1 0 0)<0 0 1] ', &
'{1 1 0)<0 0 1] ', &
'{1 0 0)<0 1 0] ', &
@ -418,13 +420,13 @@ module lattice
!--------------------------------------------------------------------------------------------------
! isotropic
integer, dimension(1), parameter, private :: &
integer, dimension(1), parameter :: &
LATTICE_ISO_NCLEAVAGESYSTEM = [3] !< # of cleavage systems per family for iso
integer, parameter, private :: &
integer, parameter :: &
LATTICE_ISO_NCLEAVAGE = sum(LATTICE_ISO_NCLEAVAGESYSTEM) !< total # of cleavage systems for iso
real(pReal), dimension(3+3,LATTICE_ISO_NCLEAVAGE), parameter, private :: &
real(pReal), dimension(3+3,LATTICE_ISO_NCLEAVAGE), parameter :: &
LATTICE_ISO_SYSTEMCLEAVAGE= reshape(real([&
! Cleavage direction Plane normal
0, 1, 0, 1, 0, 0, &
@ -435,13 +437,13 @@ module lattice
!--------------------------------------------------------------------------------------------------
! orthorhombic
integer, dimension(3), parameter, private :: &
integer, dimension(3), parameter :: &
LATTICE_ORT_NCLEAVAGESYSTEM = [1, 1, 1] !< # of cleavage systems per family for ortho
integer, parameter, private :: &
integer, parameter :: &
LATTICE_ORT_NCLEAVAGE = sum(LATTICE_ORT_NCLEAVAGESYSTEM) !< total # of cleavage systems for ortho
real(pReal), dimension(3+3,LATTICE_ORT_NCLEAVAGE), parameter, private :: &
real(pReal), dimension(3+3,LATTICE_ORT_NCLEAVAGE), parameter :: &
LATTICE_ORT_SYSTEMCLEAVAGE = reshape(real([&
! Cleavage direction Plane normal
0, 1, 0, 1, 0, 0, &
@ -541,10 +543,6 @@ module lattice
!> @brief Module initialization
!--------------------------------------------------------------------------------------------------
subroutine lattice_init
use IO, only: &
IO_error
use config, only: &
config_phase
integer :: Nphases
character(len=65536) :: &
@ -654,15 +652,7 @@ end subroutine lattice_init
!> @brief !!!!!!!DEPRECTATED!!!!!!
!--------------------------------------------------------------------------------------------------
subroutine lattice_initializeStructure(myPhase,CoverA)
use prec, only: &
tol_math_check
use math, only: &
math_sym3333to66, &
math_Voigt66to3333, &
math_cross
use IO, only: &
IO_error
integer, intent(in) :: myPhase
real(pReal), intent(in) :: &
CoverA
@ -690,9 +680,10 @@ subroutine lattice_initializeStructure(myPhase,CoverA)
call IO_error(135,el=i,ip=myPhase,ext_msg='matrix diagonal "el"ement of phase "ip"')
enddo
forall (i = 1:3) &
do i = 1,3
lattice_thermalExpansion33 (1:3,1:3,i,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),&
lattice_thermalExpansion33 (1:3,1:3,i,myPhase))
enddo
lattice_thermalConductivity33 (1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),&
lattice_thermalConductivity33 (1:3,1:3,myPhase))
@ -763,17 +754,17 @@ pure function lattice_symmetrizeC66(struct,C66)
select case(struct)
case (LATTICE_iso_ID)
forall(k=1:3)
do k=1,3
forall(j=1:3) lattice_symmetrizeC66(k,j) = C66(1,2)
lattice_symmetrizeC66(k,k) = C66(1,1)
lattice_symmetrizeC66(k+3,k+3) = 0.5_pReal*(C66(1,1)-C66(1,2))
end forall
enddo
case (LATTICE_fcc_ID,LATTICE_bcc_ID)
forall(k=1:3)
do k=1,3
forall(j=1:3) lattice_symmetrizeC66(k,j) = C66(1,2)
lattice_symmetrizeC66(k,k) = C66(1,1)
lattice_symmetrizeC66(k+3,k+3) = C66(4,4)
end forall
enddo
case (LATTICE_hex_ID)
lattice_symmetrizeC66(1,1) = C66(1,1)
lattice_symmetrizeC66(2,2) = C66(1,1)
@ -834,7 +825,9 @@ pure function lattice_symmetrize33(struct,T33)
select case(struct)
case (LATTICE_iso_ID,LATTICE_fcc_ID,LATTICE_bcc_ID)
forall(k=1:3) lattice_symmetrize33(k,k) = T33(1,1)
do k=1,3
lattice_symmetrize33(k,k) = T33(1,1)
enddo
case (LATTICE_hex_ID)
lattice_symmetrize33(1,1) = T33(1,1)
lattice_symmetrize33(2,2) = T33(1,1)
@ -854,10 +847,6 @@ end function lattice_symmetrize33
!> @brief figures whether unit quat falls into stereographic standard triangle
!--------------------------------------------------------------------------------------------------
logical pure function lattice_qInSST(Q, struct)
use, intrinsic :: &
IEEE_arithmetic
use math, only: &
math_qToRodrig
real(pReal), dimension(4), intent(in) :: Q ! orientation
integer(kind(LATTICE_undefined_ID)), intent(in) :: struct ! lattice structure
@ -888,11 +877,6 @@ end function lattice_qInSST
!> @brief calculates the disorientation for 2 unit quaternions
!--------------------------------------------------------------------------------------------------
pure function lattice_qDisorientation(Q1, Q2, struct)
use prec, only: &
tol_math_check
use math, only: &
math_qMul, &
math_qConj
real(pReal), dimension(4) :: lattice_qDisorientation
real(pReal), dimension(4), intent(in) :: &
@ -998,8 +982,6 @@ end function lattice_qDisorientation
!> @brief Characteristic shear for twinning
!--------------------------------------------------------------------------------------------------
function lattice_characteristicShear_Twin(Ntwin,structure,CoverA) result(characteristicShear)
use IO, only: &
IO_error
integer, dimension(:), intent(in) :: Ntwin !< number of active twin systems per family
character(len=*), intent(in) :: structure !< lattice structure
@ -1077,14 +1059,6 @@ end function lattice_characteristicShear_Twin
!> @brief Rotated elasticity matrices for twinning in 66-vector notation
!--------------------------------------------------------------------------------------------------
function lattice_C66_twin(Ntwin,C66,structure,CoverA)
use IO, only: &
IO_error
use math, only: &
PI, &
math_axisAngleToR, &
math_sym3333to66, &
math_66toSym3333, &
math_rotate_forward3333
integer, dimension(:), intent(in) :: Ntwin !< number of active twin systems per family
character(len=*), intent(in) :: structure !< lattice structure
@ -1125,17 +1099,6 @@ end function lattice_C66_twin
!--------------------------------------------------------------------------------------------------
function lattice_C66_trans(Ntrans,C_parent66,structure_target, &
CoverA_trans,a_bcc,a_fcc)
use prec, only: &
tol_math_check
use IO, only: &
IO_error
use math, only: &
INRAD, &
MATH_I3, &
math_axisAngleToR, &
math_sym3333to66, &
math_66toSym3333, &
math_rotate_forward3333
integer, dimension(:), intent(in) :: Ntrans !< number of active twin systems per family
character(len=*), intent(in) :: structure_target !< lattice structure
@ -1196,13 +1159,6 @@ function lattice_C66_trans(Ntrans,C_parent66,structure_target, &
! Gröger et al. 2008, Acta Materialia 56 (2008) 5412–5425, table 1
!--------------------------------------------------------------------------------------------------
function lattice_nonSchmidMatrix(Nslip,nonSchmidCoefficients,sense) result(nonSchmidMatrix)
use IO, only: &
IO_error
use math, only: &
INRAD, &
math_outer, &
math_cross, &
math_axisAngleToR
integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family
real(pReal), dimension(:), intent(in) :: nonSchmidCoefficients !< non-Schmid coefficients for projections
@ -1246,9 +1202,7 @@ end function lattice_nonSchmidMatrix
!> details only active slip systems are considered
!--------------------------------------------------------------------------------------------------
function lattice_interaction_SlipBySlip(Nslip,interactionValues,structure) result(interactionMatrix)
use IO, only: &
IO_error
integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family
real(pReal), dimension(:), intent(in) :: interactionValues !< values for slip-slip interaction
character(len=*), intent(in) :: structure !< lattice structure
@ -1468,8 +1422,6 @@ end function lattice_interaction_SlipBySlip
!> details only active twin systems are considered
!--------------------------------------------------------------------------------------------------
function lattice_interaction_TwinByTwin(Ntwin,interactionValues,structure) result(interactionMatrix)
use IO, only: &
IO_error
integer, dimension(:), intent(in) :: Ntwin !< number of active twin systems per family
real(pReal), dimension(:), intent(in) :: interactionValues !< values for twin-twin interaction
@ -1571,8 +1523,6 @@ end function lattice_interaction_TwinByTwin
!> details only active trans systems are considered
!--------------------------------------------------------------------------------------------------
function lattice_interaction_TransByTrans(Ntrans,interactionValues,structure) result(interactionMatrix)
use IO, only: &
IO_error
integer, dimension(:), intent(in) :: Ntrans !< number of active trans systems per family
real(pReal), dimension(:), intent(in) :: interactionValues !< values for trans-trans interaction
@ -1618,8 +1568,6 @@ end function lattice_interaction_TransByTrans
!> details only active slip and twin systems are considered
!--------------------------------------------------------------------------------------------------
function lattice_interaction_SlipByTwin(Nslip,Ntwin,interactionValues,structure) result(interactionMatrix)
use IO, only: &
IO_error
integer, dimension(:), intent(in) :: Nslip, & !< number of active slip systems per family
Ntwin !< number of active twin systems per family
@ -1760,8 +1708,6 @@ end function lattice_interaction_SlipByTwin
!> details only active slip and trans systems are considered
!--------------------------------------------------------------------------------------------------
function lattice_interaction_SlipByTrans(Nslip,Ntrans,interactionValues,structure) result(interactionMatrix)
use IO, only: &
IO_error
integer, dimension(:), intent(in) :: Nslip, & !< number of active slip systems per family
Ntrans !< number of active trans systems per family
@ -1818,8 +1764,6 @@ function lattice_interaction_SlipByTrans(Nslip,Ntrans,interactionValues,structur
!> details only active twin and slip systems are considered
!--------------------------------------------------------------------------------------------------
function lattice_interaction_TwinBySlip(Ntwin,Nslip,interactionValues,structure) result(interactionMatrix)
use IO, only: &
IO_error
integer, dimension(:), intent(in) :: Ntwin, & !< number of active twin systems per family
Nslip !< number of active slip systems per family
@ -1898,13 +1842,6 @@ end function lattice_interaction_TwinBySlip
!> details only active slip systems are considered
!--------------------------------------------------------------------------------------------------
function lattice_SchmidMatrix_slip(Nslip,structure,cOverA) result(SchmidMatrix)
use prec, only: &
tol_math_check
use IO, only: &
IO_error
use math, only: &
math_trace33, &
math_outer
integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family
character(len=*), intent(in) :: structure !< lattice structure
@ -1957,13 +1894,6 @@ end function lattice_SchmidMatrix_slip
!> details only active twin systems are considered
!--------------------------------------------------------------------------------------------------
function lattice_SchmidMatrix_twin(Ntwin,structure,cOverA) result(SchmidMatrix)
use prec, only: &
tol_math_check
use IO, only: &
IO_error
use math, only: &
math_trace33, &
math_outer
integer, dimension(:), intent(in) :: Ntwin !< number of active twin systems per family
character(len=*), intent(in) :: structure !< lattice structure
@ -2013,8 +1943,6 @@ function lattice_SchmidMatrix_twin(Ntwin,structure,cOverA) result(SchmidMatrix)
!> details only active twin systems are considered
!--------------------------------------------------------------------------------------------------
function lattice_SchmidMatrix_trans(Ntrans,structure_target,cOverA,a_bcc,a_fcc) result(SchmidMatrix)
use IO, only: &
IO_error
integer, dimension(:), intent(in) :: Ntrans !< number of active twin systems per family
real(pReal), intent(in) :: cOverA !< c/a ratio
@ -2041,11 +1969,7 @@ end function lattice_SchmidMatrix_trans
!> details only active cleavage systems are considered
!--------------------------------------------------------------------------------------------------
function lattice_SchmidMatrix_cleavage(Ncleavage,structure,cOverA) result(SchmidMatrix)
use math, only: &
math_outer
use IO, only: &
IO_error
integer, dimension(:), intent(in) :: Ncleavage !< number of active cleavage systems per family
character(len=*), intent(in) :: structure !< lattice structure
real(pReal), intent(in) :: cOverA !< c/a ratio
@ -2154,8 +2078,6 @@ end function lattice_slip_transverse
!> @details: This projection is used to calculate forest hardening for edge dislocations
!--------------------------------------------------------------------------------------------------
function slipProjection_transverse(Nslip,structure,cOverA) result(projection)
use math, only: &
math_inner
integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family
character(len=*), intent(in) :: structure !< lattice structure
@ -2179,8 +2101,6 @@ end function slipProjection_transverse
!> @details: This projection is used to calculate forest hardening for screw dislocations
!--------------------------------------------------------------------------------------------------
function slipProjection_direction(Nslip,structure,cOverA) result(projection)
use math, only: &
math_inner
integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family
character(len=*), intent(in) :: structure !< lattice structure
@ -2204,9 +2124,7 @@ end function slipProjection_direction
!> @details Order: Direction, plane (normal), and common perpendicular
!--------------------------------------------------------------------------------------------------
function coordinateSystem_slip(Nslip,structure,cOverA) result(coordinateSystem)
use IO, only: &
IO_error
integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family
character(len=*), intent(in) :: structure !< lattice structure
real(pReal), intent(in) :: cOverA !< c/a ratio
@ -2249,8 +2167,6 @@ end function coordinateSystem_slip
!> @brief Populates reduced interaction matrix
!--------------------------------------------------------------------------------------------------
function buildInteraction(reacting_used,acting_used,reacting_max,acting_max,values,matrix)
use IO, only: &
IO_error
integer, dimension(:), intent(in) :: &
reacting_used, & !< # of reacting systems per family as specified in material.config
@ -2295,10 +2211,6 @@ end function buildInteraction
!> @details Order: Direction, plane (normal), and common perpendicular
!--------------------------------------------------------------------------------------------------
function buildCoordinateSystem(active,complete,system,structure,cOverA)
use IO, only: &
IO_error
use math, only: &
math_cross
integer, dimension(:), intent(in) :: &
active, &
@ -2370,16 +2282,6 @@ end function buildCoordinateSystem
! set a_bcc = 0.0 for fcc -> hex transformation
!--------------------------------------------------------------------------------------------------
subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_fcc,a_bcc)
use prec, only: &
dEq0
use math, only: &
math_cross, &
math_outer, &
math_axisAngleToR, &
INRAD, &
MATH_I3
use IO, only: &
IO_error
integer, dimension(:), intent(in) :: &
Ntrans

34
src/list.f90
View File

@ -3,8 +3,8 @@
!> @brief linked list
!--------------------------------------------------------------------------------------------------
module list
use prec, only: &
pReal
use prec
use IO
implicit none
private
@ -65,10 +65,6 @@ contains
!! to lower case. The data is not stored in the new element but in the current.
!--------------------------------------------------------------------------------------------------
subroutine add(this,string)
use IO, only: &
IO_isBlank, &
IO_lc, &
IO_stringPos
class(tPartitionedStringList), target, intent(in) :: this
character(len=*), intent(in) :: string
@ -157,8 +153,6 @@ end subroutine finalizeArray
!> @brief reports wether a given key (string value at first position) exists in the list
!--------------------------------------------------------------------------------------------------
logical function keyExists(this,key)
use IO, only: &
IO_stringValue
class(tPartitionedStringList), target, intent(in) :: this
character(len=*), intent(in) :: key
@ -180,8 +174,6 @@ end function keyExists
!> @details traverses list and counts each occurrence of specified key
!--------------------------------------------------------------------------------------------------
integer function countKeys(this,key)
use IO, only: &
IO_stringValue
class(tPartitionedStringList), target, intent(in) :: this
character(len=*), intent(in) :: key
@ -205,10 +197,6 @@ end function countKeys
!! error unless default is given
!--------------------------------------------------------------------------------------------------
real(pReal) function getFloat(this,key,defaultVal)
use IO, only : &
IO_error, &
IO_stringValue, &
IO_FloatValue
class(tPartitionedStringList), target, intent(in) :: this
character(len=*), intent(in) :: key
@ -241,10 +229,6 @@ end function getFloat
!! error unless default is given
!--------------------------------------------------------------------------------------------------
integer function getInt(this,key,defaultVal)
use IO, only: &
IO_error, &
IO_stringValue, &
IO_IntValue
class(tPartitionedStringList), target, intent(in) :: this
character(len=*), intent(in) :: key
@ -278,9 +262,6 @@ end function getInt
!! the individual chunks are returned
!--------------------------------------------------------------------------------------------------
character(len=65536) function getString(this,key,defaultVal,raw)
use IO, only: &
IO_error, &
IO_stringValue
class(tPartitionedStringList), target, intent(in) :: this
character(len=*), intent(in) :: key
@ -327,10 +308,6 @@ end function getString
!! values from the last occurrence. If key is not found exits with error unless default is given.
!--------------------------------------------------------------------------------------------------
function getFloats(this,key,defaultVal,requiredSize)
use IO, only: &
IO_error, &
IO_stringValue, &
IO_FloatValue
real(pReal), dimension(:), allocatable :: getFloats
class(tPartitionedStringList), target, intent(in) :: this
@ -376,10 +353,6 @@ end function getFloats
!! values from the last occurrence. If key is not found exits with error unless default is given.
!--------------------------------------------------------------------------------------------------
function getInts(this,key,defaultVal,requiredSize)
use IO, only: &
IO_error, &
IO_stringValue, &
IO_IntValue
integer, dimension(:), allocatable :: getInts
class(tPartitionedStringList), target, intent(in) :: this
@ -426,9 +399,6 @@ end function getInts
!! If raw is true, the the complete string is returned, otherwise the individual chunks are returned
!--------------------------------------------------------------------------------------------------
function getStrings(this,key,defaultVal,raw)
use IO, only: &
IO_error, &
IO_StringValue
character(len=65536),dimension(:), allocatable :: getStrings
class(tPartitionedStringList),target, intent(in) :: this

29
src/math.f90
View File

@ -10,12 +10,20 @@ module math
use future
implicit none
real(pReal), parameter, public :: PI = acos(-1.0_pReal) !< ratio of a circle's circumference to its diameter
real(pReal), parameter, public :: INDEG = 180.0_pReal/PI !< conversion from radian into degree
real(pReal), parameter, public :: INRAD = PI/180.0_pReal !< conversion from degree into radian
complex(pReal), parameter, public :: TWOPIIMG = cmplx(0.0_pReal,2.0_pReal*PI) !< Re(0.0), Im(2xPi)
public
#if __INTEL_COMPILER >= 1900
! do not make use associated entities available to other modules
private :: &
prec, &
future
#endif
real(pReal), dimension(3,3), parameter, public :: &
real(pReal), parameter :: PI = acos(-1.0_pReal) !< ratio of a circle's circumference to its diameter
real(pReal), parameter :: INDEG = 180.0_pReal/PI !< conversion from radian into degree
real(pReal), parameter :: INRAD = PI/180.0_pReal !< conversion from degree into radian
complex(pReal), parameter :: TWOPIIMG = cmplx(0.0_pReal,2.0_pReal*PI) !< Re(0.0), Im(2xPi)
real(pReal), dimension(3,3), parameter :: &
MATH_I3 = reshape([&
1.0_pReal,0.0_pReal,0.0_pReal, &
0.0_pReal,1.0_pReal,0.0_pReal, &
@ -75,7 +83,7 @@ module math
!---------------------------------------------------------------------------------------------------
private :: &
math_check
unitTest
contains
@ -116,14 +124,15 @@ subroutine math_init
write(6,'(a,4(/,26x,f17.14),/)') ' start of random sequence: ', randTest
call random_seed(put = randInit)
call math_check
call unitTest
end subroutine math_init
!--------------------------------------------------------------------------------------------------
!> @brief check correctness of (some) math functions
!--------------------------------------------------------------------------------------------------
subroutine math_check
subroutine unitTest
use IO, only: IO_error
character(len=64) :: error_msg
@ -145,7 +154,7 @@ subroutine math_check
call IO_error(401,ext_msg=error_msg)
endif
end subroutine math_check
end subroutine unitTest
!--------------------------------------------------------------------------------------------------
@ -274,6 +283,7 @@ pure function math_identity2nd(dimen)
end function math_identity2nd
!--------------------------------------------------------------------------------------------------
!> @brief symmetric fourth rank identity tensor of specified dimension
! from http://en.wikipedia.org/wiki/Tensor_derivative_(continuum_mechanics)#Derivative_of_a_second-order_tensor_with_respect_to_itself
@ -626,6 +636,7 @@ pure function math_skew33(m)
end function math_skew33
!--------------------------------------------------------------------------------------------------
!> @brief hydrostatic part of a 33 matrix
!--------------------------------------------------------------------------------------------------

103
src/plastic_dislotwin.f90
View File

@ -8,17 +8,26 @@
!> @details to be done
!--------------------------------------------------------------------------------------------------
module plastic_dislotwin
use prec, only: &
pReal
use prec
use debug
use math
use IO
use material
use config
use lattice
#if defined(PETSc) || defined(DAMASK_HDF5)
use results
#endif
implicit none
private
integer, dimension(:,:), allocatable, target, public :: &
plastic_dislotwin_sizePostResult !< size of each post result output
character(len=64), dimension(:,:), allocatable, target, public :: &
plastic_dislotwin_output !< name of each post result output
real(pReal), parameter, private :: &
real(pReal), parameter :: &
kB = 1.38e-23_pReal !< Boltzmann constant in J/Kelvin
enum, bind(c)
@ -39,7 +48,7 @@ module plastic_dislotwin
f_tr_ID
end enum
type, private :: tParameters
type :: tParameters
real(pReal) :: &
mu, &
nu, &
@ -119,7 +128,7 @@ module plastic_dislotwin
dipoleFormation !< flag indicating consideration of dipole formation
end type !< container type for internal constitutive parameters
type, private :: tDislotwinState
type :: tDislotwinState
real(pReal), dimension(:,:), pointer :: &
rho_mob, &
rho_dip, &
@ -128,7 +137,7 @@ module plastic_dislotwin
f_tr
end type tDislotwinState
type, private :: tDislotwinMicrostructure
type :: tDislotwinMicrostructure
real(pReal), dimension(:,:), allocatable :: &
Lambda_sl, & !* mean free path between 2 obstacles seen by a moving dislocation
Lambda_tw, & !* mean free path between 2 obstacles seen by a growing twin
@ -144,11 +153,11 @@ module plastic_dislotwin
!--------------------------------------------------------------------------------------------------
! containers for parameters and state
type(tParameters), allocatable, dimension(:), private :: param
type(tDislotwinState), allocatable, dimension(:), private :: &
type(tParameters), allocatable, dimension(:) :: param
type(tDislotwinState), allocatable, dimension(:) :: &
dotState, &
state
type(tDislotwinMicrostructure), allocatable, dimension(:), private :: dependentState
type(tDislotwinMicrostructure), allocatable, dimension(:) :: dependentState
public :: &
plastic_dislotwin_init, &
@ -158,10 +167,6 @@ module plastic_dislotwin
plastic_dislotwin_dotState, &
plastic_dislotwin_postResults, &
plastic_dislotwin_results
private :: &
kinetics_slip, &
kinetics_twin, &
kinetics_trans
contains
@ -171,24 +176,6 @@ contains
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
subroutine plastic_dislotwin_init
use prec, only: &
pStringLen, &
dEq0, &
dNeq0, &
dNeq
use debug, only: &
debug_level,&
debug_constitutive,&
debug_levelBasic
use math, only: &
math_expand,&
PI
use IO, only: &
IO_error
use material
use config, only: &
config_phase
use lattice
integer :: &
Ninstance, &
@ -591,10 +578,6 @@ end subroutine plastic_dislotwin_init
!> @brief returns the homogenized elasticity matrix
!--------------------------------------------------------------------------------------------------
function plastic_dislotwin_homogenizedC(ipc,ip,el) result(homogenizedC)
use material, only: &
material_phase, &
phase_plasticityInstance, &
phasememberAt
real(pReal), dimension(6,6) :: &
homogenizedC
@ -634,14 +617,6 @@ end function plastic_dislotwin_homogenizedC
!> @brief calculates plastic velocity gradient and its tangent
!--------------------------------------------------------------------------------------------------
subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of)
use prec, only: &
tol_math_check, &
dNeq0
use math, only: &
math_eigenValuesVectorsSym, &
math_outer, &
math_symmetric33, &
math_mul33xx33
real(pReal), dimension(3,3), intent(out) :: Lp
real(pReal), dimension(3,3,3,3), intent(out) :: dLp_dMp
@ -757,13 +732,6 @@ end subroutine plastic_dislotwin_LpAndItsTangent
!> @brief calculates the rate of change of microstructure
!--------------------------------------------------------------------------------------------------
subroutine plastic_dislotwin_dotState(Mp,T,instance,of)
use prec, only: &
tol_math_check, &
dEq0
use math, only: &
math_clip, &
math_mul33xx33, &
PI
real(pReal), dimension(3,3), intent(in):: &
Mp !< Mandel stress
@ -854,8 +822,6 @@ end subroutine plastic_dislotwin_dotState
!> @brief calculates derived quantities from state
!--------------------------------------------------------------------------------------------------
subroutine plastic_dislotwin_dependentState(T,instance,of)
use math, only: &
PI
integer, intent(in) :: &
instance, &
@ -868,13 +834,13 @@ subroutine plastic_dislotwin_dependentState(T,instance,of)
real(pReal) :: &
sumf_twin,SFE,sumf_trans
real(pReal), dimension(param(instance)%sum_N_sl) :: &
inv_lambda_sl_sl, & !* 1/mean free distance between 2 forest dislocations seen by a moving dislocation
inv_lambda_sl_tw, & !* 1/mean free distance between 2 twin stacks from different systems seen by a moving dislocation
inv_lambda_sl_tr !* 1/mean free distance between 2 martensite lamellar from different systems seen by a moving dislocation
inv_lambda_sl_sl, & !< 1/mean free distance between 2 forest dislocations seen by a moving dislocation
inv_lambda_sl_tw, & !< 1/mean free distance between 2 twin stacks from different systems seen by a moving dislocation
inv_lambda_sl_tr !< 1/mean free distance between 2 martensite lamellar from different systems seen by a moving dislocation
real(pReal), dimension(param(instance)%sum_N_tw) :: &
inv_lambda_tw_tw !* 1/mean free distance between 2 twin stacks from different systems seen by a growing twin
inv_lambda_tw_tw !< 1/mean free distance between 2 twin stacks from different systems seen by a growing twin
real(pReal), dimension(param(instance)%sum_N_tr) :: &
inv_lambda_tr_tr !* 1/mean free distance between 2 martensite stacks from different systems seen by a growing martensite (1/lambda_trans)
inv_lambda_tr_tr !< 1/mean free distance between 2 martensite stacks from different systems seen by a growing martensite
real(pReal), dimension(:), allocatable :: &
x0, &
@ -967,12 +933,6 @@ end subroutine plastic_dislotwin_dependentState
!> @brief return array of constitutive results
!--------------------------------------------------------------------------------------------------
function plastic_dislotwin_postResults(Mp,T,instance,of) result(postResults)
use prec, only: &
tol_math_check, &
dEq0
use math, only: &
PI, &
math_mul33xx33
real(pReal), dimension(3,3),intent(in) :: &
Mp !< 2nd Piola Kirchhoff stress tensor in Mandel notation
@ -1050,8 +1010,6 @@ end function plastic_dislotwin_postResults
!--------------------------------------------------------------------------------------------------
subroutine plastic_dislotwin_results(instance,group)
#if defined(PETSc) || defined(DAMASK_HDF5)
use results, only: &
results_writeDataset
integer, intent(in) :: instance
character(len=*) :: group
@ -1112,11 +1070,6 @@ end subroutine plastic_dislotwin_results
!--------------------------------------------------------------------------------------------------
pure subroutine kinetics_slip(Mp,T,instance,of, &
dot_gamma_sl,ddot_gamma_dtau_slip,tau_slip)
use prec, only: &
tol_math_check, &
dNeq0
use math, only: &
math_mul33xx33
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
@ -1190,11 +1143,6 @@ end subroutine kinetics_slip
!--------------------------------------------------------------------------------------------------
pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,instance,of,&
dot_gamma_twin,ddot_gamma_dtau_twin)
use prec, only: &
tol_math_check, &
dNeq0
use math, only: &
math_mul33xx33
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
@ -1261,11 +1209,6 @@ end subroutine kinetics_twin
!--------------------------------------------------------------------------------------------------
pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,instance,of,&
dot_gamma_tr,ddot_gamma_dtau_trans)
use prec, only: &
tol_math_check, &
dNeq0
use math, only: &
math_mul33xx33
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress

63
src/plastic_isotropic.f90
View File

@ -8,11 +8,19 @@
!! untextured polycrystal
!--------------------------------------------------------------------------------------------------
module plastic_isotropic
use prec, only: &
pReal
use prec
use debug
use math
use IO
use material
use config
#if defined(PETSc) || defined(DAMASK_HDF5)
use results
#endif
implicit none
private
integer, dimension(:,:), allocatable, target, public :: &
plastic_isotropic_sizePostResult !< size of each post result output
character(len=64), dimension(:,:), allocatable, target, public :: &
@ -25,7 +33,7 @@ module plastic_isotropic
dot_gamma_ID
end enum
type, private :: tParameters
type :: tParameters
real(pReal) :: &
M, & !< Taylor factor
xi_0, & !< initial critical stress
@ -49,7 +57,7 @@ module plastic_isotropic
dilatation
end type tParameters
type, private :: tIsotropicState
type :: tIsotropicState
real(pReal), pointer, dimension(:) :: &
xi, &
gamma
@ -57,8 +65,8 @@ module plastic_isotropic
!--------------------------------------------------------------------------------------------------
! containers for parameters and state
type(tParameters), allocatable, dimension(:), private :: param
type(tIsotropicState), allocatable, dimension(:), private :: &
type(tParameters), allocatable, dimension(:) :: param
type(tIsotropicState), allocatable, dimension(:) :: &
dotState, &
state
@ -77,25 +85,7 @@ contains
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
subroutine plastic_isotropic_init
use prec, only: &
pStringLen
use debug, only: &
#ifdef DEBUG
debug_e, &
debug_i, &
debug_g, &
debug_levelExtensive, &
#endif
debug_level, &
debug_constitutive, &
debug_levelBasic
use IO, only: &
IO_error
use material
use config, only: &
config_phase
use lattice
integer :: &
Ninstance, &
p, i, &
@ -235,16 +225,6 @@ end subroutine plastic_isotropic_init
!> @brief calculates plastic velocity gradient and its tangent
!--------------------------------------------------------------------------------------------------
subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
#ifdef DEBUG
use debug, only: &
debug_level, &
debug_constitutive,&
debug_levelExtensive, &
debug_levelSelective
#endif
use math, only: &
math_deviatoric33, &
math_mul33xx33
real(pReal), dimension(3,3), intent(out) :: &
Lp !< plastic velocity gradient
@ -307,10 +287,6 @@ end subroutine plastic_isotropic_LpAndItsTangent
! ToDo: Rename Tstar to Mi?
!--------------------------------------------------------------------------------------------------
subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dTstar,Tstar,instance,of)
use math, only: &
math_I3, &
math_spherical33, &
math_mul33xx33
real(pReal), dimension(3,3), intent(out) :: &
Li !< inleastic velocity gradient
@ -362,11 +338,6 @@ subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dTstar,Tstar,instance,of)
!> @brief calculates the rate of change of microstructure
!--------------------------------------------------------------------------------------------------
subroutine plastic_isotropic_dotState(Mp,instance,of)
use prec, only: &
dEq0
use math, only: &
math_mul33xx33, &
math_deviatoric33
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
@ -416,9 +387,6 @@ end subroutine plastic_isotropic_dotState
!> @brief return array of constitutive results
!--------------------------------------------------------------------------------------------------
function plastic_isotropic_postResults(Mp,instance,of) result(postResults)
use math, only: &
math_mul33xx33, &
math_deviatoric33
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
@ -468,7 +436,6 @@ end function plastic_isotropic_postResults
!--------------------------------------------------------------------------------------------------
subroutine plastic_isotropic_results(instance,group)
#if defined(PETSc) || defined(DAMASKHDF5)
use results
integer, intent(in) :: instance
character(len=*), intent(in) :: group

62
src/plastic_kinematichardening.f90
View File

@ -6,11 +6,20 @@
!! and a Voce-type kinematic hardening rule
!--------------------------------------------------------------------------------------------------
module plastic_kinehardening
use prec, only: &
pReal
use prec
use debug
use math
use IO
use material
use config
use lattice
#if defined(PETSc) || defined(DAMASK_HDF5)
use results
#endif
implicit none
private
integer, dimension(:,:), allocatable, target, public :: &
plastic_kinehardening_sizePostResult !< size of each post result output
character(len=64), dimension(:,:), allocatable, target, public :: &
@ -29,7 +38,7 @@ module plastic_kinehardening
resolvedstress_ID
end enum
type, private :: tParameters
type :: tParameters
real(pReal) :: &
gdot0, & !< reference shear strain rate for slip
n, & !< stress exponent for slip
@ -59,7 +68,7 @@ module plastic_kinehardening
outputID !< ID of each post result output
end type tParameters
type, private :: tKinehardeningState
type :: tKinehardeningState
real(pReal), pointer, dimension(:,:) :: & !< vectors along NipcMyInstance
crss, & !< critical resolved stress
crss_back, & !< critical resolved back stress
@ -71,8 +80,8 @@ module plastic_kinehardening
!--------------------------------------------------------------------------------------------------
! containers for parameters and state
type(tParameters), allocatable, dimension(:), private :: param
type(tKinehardeningState), allocatable, dimension(:), private :: &
type(tParameters), allocatable, dimension(:) :: param
type(tKinehardeningState), allocatable, dimension(:) :: &
dotState, &
deltaState, &
state
@ -84,8 +93,6 @@ module plastic_kinehardening
plastic_kinehardening_deltaState, &
plastic_kinehardening_postResults, &
plastic_kinehardening_results
private :: &
kinetics
contains
@ -95,27 +102,6 @@ contains
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
subroutine plastic_kinehardening_init
use prec, only: &
dEq0, &
pStringLen
use debug, only: &
#ifdef DEBUG
debug_e, &
debug_i, &
debug_g, &
debug_levelExtensive, &
#endif
debug_level, &
debug_constitutive,&
debug_levelBasic
use math, only: &
math_expand
use IO, only: &
IO_error
use material
use config, only: &
config_phase
use lattice
integer :: &
Ninstance, &
@ -417,16 +403,6 @@ end subroutine plastic_kinehardening_dotState
!> @brief calculates (instantaneous) incremental change of microstructure
!--------------------------------------------------------------------------------------------------
subroutine plastic_kinehardening_deltaState(Mp,instance,of)
use prec, only: &
dNeq, &
dEq0
#ifdef DEBUG
use debug, only: &
debug_level, &
debug_constitutive,&
debug_levelExtensive, &
debug_levelSelective
#endif
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
@ -475,8 +451,6 @@ end subroutine plastic_kinehardening_deltaState
!> @brief return array of constitutive results
!--------------------------------------------------------------------------------------------------
function plastic_kinehardening_postResults(Mp,instance,of) result(postResults)
use math, only: &
math_mul33xx33
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
@ -535,8 +509,6 @@ end function plastic_kinehardening_postResults
!--------------------------------------------------------------------------------------------------
subroutine plastic_kinehardening_results(instance,group)
#if defined(PETSc) || defined(DAMASK_HDF5)
use results, only: &
results_writeDataset
integer, intent(in) :: instance
character(len=*) :: group
@ -585,10 +557,6 @@ end subroutine plastic_kinehardening_results
!--------------------------------------------------------------------------------------------------
pure subroutine kinetics(Mp,instance,of, &
gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg)
use prec, only: &
dNeq0
use math, only: &
math_mul33xx33
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress

7
src/plastic_none.f90
View File

@ -5,6 +5,8 @@
!> @brief Dummy plasticity for purely elastic material
!--------------------------------------------------------------------------------------------------
module plastic_none
use material
use debug
implicit none
private
@ -19,11 +21,6 @@ contains
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
subroutine plastic_none_init
use debug, only: &
debug_level, &
debug_constitutive, &
debug_levelBasic
use material
integer :: &
Ninstance, &

50
src/plastic_phenopowerlaw.f90
View File

@ -5,11 +5,20 @@
!> @brief phenomenological crystal plasticity formulation using a powerlaw fitting
!--------------------------------------------------------------------------------------------------
module plastic_phenopowerlaw
use prec, only: &
pReal
use prec
use debug
use math
use IO
use material
use config
use lattice
#if defined(PETSc) || defined(DAMASK_HDF5)
use results
#endif
implicit none
private
integer, dimension(:,:), allocatable, target, public :: &
plastic_phenopowerlaw_sizePostResult !< size of each post result output
character(len=64), dimension(:,:), allocatable, target, public :: &
@ -28,7 +37,7 @@ module plastic_phenopowerlaw
resolvedstress_twin_ID
end enum
type, private :: tParameters
type :: tParameters
real(pReal) :: &
gdot0_slip, & !< reference shear strain rate for slip
gdot0_twin, & !< reference shear strain rate for twin
@ -73,7 +82,7 @@ module plastic_phenopowerlaw
outputID !< ID of each post result output
end type tParameters
type, private :: tPhenopowerlawState
type :: tPhenopowerlawState
real(pReal), pointer, dimension(:,:) :: &
xi_slip, &
xi_twin, &
@ -83,8 +92,8 @@ module plastic_phenopowerlaw
!--------------------------------------------------------------------------------------------------
! containers for parameters and state
type(tParameters), allocatable, dimension(:), private :: param
type(tPhenopowerlawState), allocatable, dimension(:), private :: &
type(tParameters), allocatable, dimension(:) :: param
type(tPhenopowerlawState), allocatable, dimension(:) :: &
dotState, &
state
@ -94,9 +103,6 @@ module plastic_phenopowerlaw
plastic_phenopowerlaw_dotState, &
plastic_phenopowerlaw_postResults, &
plastic_phenopowerlaw_results
private :: &
kinetics_slip, &
kinetics_twin
contains
@ -106,20 +112,6 @@ contains
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
subroutine plastic_phenopowerlaw_init
use prec, only: &
pStringLen
use debug, only: &
debug_level, &
debug_constitutive,&
debug_levelBasic
use math, only: &
math_expand
use IO, only: &
IO_error
use material
use config, only: &
config_phase
use lattice
integer :: &
Ninstance, &
@ -484,8 +476,6 @@ end subroutine plastic_phenopowerlaw_dotState
!> @brief return array of constitutive results
!--------------------------------------------------------------------------------------------------
function plastic_phenopowerlaw_postResults(Mp,instance,of) result(postResults)
use math, only: &
math_mul33xx33
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
@ -552,8 +542,6 @@ end function plastic_phenopowerlaw_postResults
!--------------------------------------------------------------------------------------------------
subroutine plastic_phenopowerlaw_results(instance,group)
#if defined(PETSc) || defined(DAMASK_HDF5)
use results, only: &
results_writeDataset
integer, intent(in) :: instance
character(len=*), intent(in) :: group
@ -598,10 +586,6 @@ end subroutine plastic_phenopowerlaw_results
!--------------------------------------------------------------------------------------------------
pure subroutine kinetics_slip(Mp,instance,of, &
gdot_slip_pos,gdot_slip_neg,dgdot_dtau_slip_pos,dgdot_dtau_slip_neg)
use prec, only: &
dNeq0
use math, only: &
math_mul33xx33
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
@ -674,10 +658,6 @@ end subroutine kinetics_slip
!--------------------------------------------------------------------------------------------------
pure subroutine kinetics_twin(Mp,instance,of,&
gdot_twin,dgdot_dtau_twin)
use prec, only: &
dNeq0
use math, only: &
math_mul33xx33
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress

337
src/quaternions.f90
View File

@ -3,27 +3,27 @@
! Modified 2017-2019, Martin Diehl/Max-Planck-Institut für Eisenforschung GmbH
! All rights reserved.
!
! Redistribution and use in source and binary forms, with or without modification, are
! Redistribution and use in source and binary forms, with or without modification, are
! permitted provided that the following conditions are met:
!
! - Redistributions of source code must retain the above copyright notice, this list
! - Redistributions of source code must retain the above copyright notice, this list
! of conditions and the following disclaimer.
! - Redistributions in binary form must reproduce the above copyright notice, this
! list of conditions and the following disclaimer in the documentation and/or
! - Redistributions in binary form must reproduce the above copyright notice, this
! list of conditions and the following disclaimer in the documentation and/or
! other materials provided with the distribution.
! - Neither the names of Marc De Graef, Carnegie Mellon University nor the names
! of its contributors may be used to endorse or promote products derived from
! - Neither the names of Marc De Graef, Carnegie Mellon University nor the names
! of its contributors may be used to endorse or promote products derived from
! this software without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
! IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
! ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
! LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
! SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
! OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
! IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
! ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
! LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
! SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
! OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE
! USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
! ###################################################################
@ -34,58 +34,57 @@
!> @details w is the real part, (x, y, z) are the imaginary parts.
!---------------------------------------------------------------------------------------------------
module quaternions
use prec, only: &
pReal
use future
use prec
use future
implicit none
public
real(pReal), parameter, public :: P = -1.0_pReal !< parameter for orientation conversion.
type, public :: quaternion
real(pReal) :: w = 0.0_pReal
real(pReal) :: x = 0.0_pReal
real(pReal) :: y = 0.0_pReal
real(pReal) :: z = 0.0_pReal
implicit none
public
contains
procedure, private :: add__
procedure, private :: pos__
generic, public :: operator(+) => add__,pos__
real(pReal), parameter, public :: P = -1.0_pReal !< parameter for orientation conversion.
procedure, private :: sub__
procedure, private :: neg__
generic, public :: operator(-) => sub__,neg__
type, public :: quaternion
real(pReal) :: w = 0.0_pReal
real(pReal) :: x = 0.0_pReal
real(pReal) :: y = 0.0_pReal
real(pReal) :: z = 0.0_pReal
procedure, private :: mul_quat__
procedure, private :: mul_scal__
generic, public :: operator(*) => mul_quat__, mul_scal__
procedure, private :: div_quat__
procedure, private :: div_scal__
generic, public :: operator(/) => div_quat__, div_scal__
contains
procedure, private :: add__
procedure, private :: pos__
generic, public :: operator(+) => add__,pos__
procedure, private :: eq__
generic, public :: operator(==) => eq__
procedure, private :: sub__
procedure, private :: neg__
generic, public :: operator(-) => sub__,neg__
procedure, private :: neq__
generic, public :: operator(/=) => neq__
procedure, private :: mul_quat__
procedure, private :: mul_scal__
generic, public :: operator(*) => mul_quat__, mul_scal__
procedure, private :: pow_quat__
procedure, private :: pow_scal__
generic, public :: operator(**) => pow_quat__, pow_scal__
procedure, private :: div_quat__
procedure, private :: div_scal__
generic, public :: operator(/) => div_quat__, div_scal__
procedure, public :: abs__
procedure, public :: dot_product__
procedure, public :: conjg__
procedure, public :: exp__
procedure, public :: log__
procedure, private :: eq__
generic, public :: operator(==) => eq__
procedure, public :: homomorphed => quat_homomorphed
procedure, private :: neq__
generic, public :: operator(/=) => neq__
end type
procedure, private :: pow_quat__
procedure, private :: pow_scal__
generic, public :: operator(**) => pow_quat__, pow_scal__
procedure, public :: abs__
procedure, public :: dot_product__
procedure, public :: conjg__
procedure, public :: exp__
procedure, public :: log__
procedure, public :: homomorphed => quat_homomorphed
end type
interface assignment (=)
module procedure assign_quat__
@ -124,12 +123,12 @@ contains
!---------------------------------------------------------------------------------------------------
type(quaternion) pure function init__(array)
real(pReal), intent(in), dimension(4) :: array
init__%w=array(1)
init__%x=array(2)
init__%y=array(3)
init__%z=array(4)
real(pReal), intent(in), dimension(4) :: array
init__%w=array(1)
init__%x=array(2)
init__%y=array(3)
init__%z=array(4)
end function init__
@ -139,14 +138,14 @@ end function init__
!---------------------------------------------------------------------------------------------------
elemental subroutine assign_quat__(self,other)
type(quaternion), intent(out) :: self
type(quaternion), intent(in) :: other
self%w = other%w
self%x = other%x
self%y = other%y
self%z = other%z
type(quaternion), intent(out) :: self
type(quaternion), intent(in) :: other
self%w = other%w
self%x = other%x
self%y = other%y
self%z = other%z
end subroutine assign_quat__
@ -155,14 +154,14 @@ end subroutine assign_quat__
!---------------------------------------------------------------------------------------------------
pure subroutine assign_vec__(self,other)
type(quaternion), intent(out) :: self
real(pReal), intent(in), dimension(4) :: other
self%w = other(1)
self%x = other(2)
self%y = other(3)
self%z = other(4)
type(quaternion), intent(out) :: self
real(pReal), intent(in), dimension(4) :: other
self%w = other(1)
self%x = other(2)
self%y = other(3)
self%z = other(4)
end subroutine assign_vec__
@ -171,13 +170,13 @@ end subroutine assign_vec__
!---------------------------------------------------------------------------------------------------
type(quaternion) elemental function add__(self,other)
class(quaternion), intent(in) :: self,other
add__%w = self%w + other%w
add__%x = self%x + other%x
add__%y = self%y + other%y
add__%z = self%z + other%z
class(quaternion), intent(in) :: self,other
add__%w = self%w + other%w
add__%x = self%x + other%x
add__%y = self%y + other%y
add__%z = self%z + other%z
end function add__
@ -186,13 +185,13 @@ end function add__
!---------------------------------------------------------------------------------------------------
type(quaternion) elemental function pos__(self)
class(quaternion), intent(in) :: self
pos__%w = self%w
pos__%x = self%x
pos__%y = self%y
pos__%z = self%z
class(quaternion), intent(in) :: self
pos__%w = self%w
pos__%x = self%x
pos__%y = self%y
pos__%z = self%z
end function pos__
@ -201,13 +200,13 @@ end function pos__
!---------------------------------------------------------------------------------------------------
type(quaternion) elemental function sub__(self,other)
class(quaternion), intent(in) :: self,other
sub__%w = self%w - other%w
sub__%x = self%x - other%x
sub__%y = self%y - other%y
sub__%z = self%z - other%z
class(quaternion), intent(in) :: self,other
sub__%w = self%w - other%w
sub__%x = self%x - other%x
sub__%y = self%y - other%y
sub__%z = self%z - other%z
end function sub__
@ -216,13 +215,13 @@ end function sub__
!---------------------------------------------------------------------------------------------------
type(quaternion) elemental function neg__(self)
class(quaternion), intent(in) :: self
neg__%w = -self%w
neg__%x = -self%x
neg__%y = -self%y
neg__%z = -self%z
class(quaternion), intent(in) :: self
neg__%w = -self%w
neg__%x = -self%x
neg__%y = -self%y
neg__%z = -self%z
end function neg__
@ -231,13 +230,13 @@ end function neg__
!---------------------------------------------------------------------------------------------------
type(quaternion) elemental function mul_quat__(self,other)
class(quaternion), intent(in) :: self, other
class(quaternion), intent(in) :: self, other
mul_quat__%w = self%w*other%w - self%x*other%x - self%y*other%y - self%z*other%z
mul_quat__%x = self%w*other%x + self%x*other%w + P * (self%y*other%z - self%z*other%y)
mul_quat__%y = self%w*other%y + self%y*other%w + P * (self%z*other%x - self%x*other%z)
mul_quat__%z = self%w*other%z + self%z*other%w + P * (self%x*other%y - self%y*other%x)
mul_quat__%w = self%w*other%w - self%x*other%x - self%y*other%y - self%z*other%z
mul_quat__%x = self%w*other%x + self%x*other%w + P * (self%y*other%z - self%z*other%y)
mul_quat__%y = self%w*other%y + self%y*other%w + P * (self%z*other%x - self%x*other%z)
mul_quat__%z = self%w*other%z + self%z*other%w + P * (self%x*other%y - self%y*other%x)
end function mul_quat__
@ -246,14 +245,14 @@ end function mul_quat__
!---------------------------------------------------------------------------------------------------
type(quaternion) elemental function mul_scal__(self,scal)
class(quaternion), intent(in) :: self
real(pReal), intent(in) :: scal
class(quaternion), intent(in) :: self
real(pReal), intent(in) :: scal
mul_scal__%w = self%w*scal
mul_scal__%x = self%x*scal
mul_scal__%y = self%y*scal
mul_scal__%z = self%z*scal
mul_scal__%w = self%w*scal
mul_scal__%x = self%x*scal
mul_scal__%y = self%y*scal
mul_scal__%z = self%z*scal
end function mul_scal__
@ -262,9 +261,9 @@ end function mul_scal__
!---------------------------------------------------------------------------------------------------
type(quaternion) elemental function div_quat__(self,other)
class(quaternion), intent(in) :: self, other
class(quaternion), intent(in) :: self, other
div_quat__ = self * (conjg(other)/(abs(other)**2.0_pReal))
div_quat__ = self * (conjg(other)/(abs(other)**2.0_pReal))
end function div_quat__
@ -274,10 +273,10 @@ end function div_quat__
!---------------------------------------------------------------------------------------------------
type(quaternion) elemental function div_scal__(self,scal)
class(quaternion), intent(in) :: self
real(pReal), intent(in) :: scal
class(quaternion), intent(in) :: self
real(pReal), intent(in) :: scal
div_scal__ = [self%w,self%x,self%y,self%z]/scal
div_scal__ = [self%w,self%x,self%y,self%z]/scal
end function div_scal__
@ -286,14 +285,12 @@ end function div_scal__
!> equality of two quaternions
!---------------------------------------------------------------------------------------------------
logical elemental function eq__(self,other)
use prec, only: &
dEq
class(quaternion), intent(in) :: self,other
class(quaternion), intent(in) :: self,other
eq__ = all(dEq([ self%w, self%x, self%y, self%z], &
[other%w,other%x,other%y,other%z]))
eq__ = all(dEq([ self%w, self%x, self%y, self%z], &
[other%w,other%x,other%y,other%z]))
end function eq__
@ -302,10 +299,10 @@ end function eq__
!---------------------------------------------------------------------------------------------------
logical elemental function neq__(self,other)
class(quaternion), intent(in) :: self,other
class(quaternion), intent(in) :: self,other
neq__ = .not. self%eq__(other)
neq__ = .not. self%eq__(other)
end function neq__
@ -314,11 +311,11 @@ end function neq__
!---------------------------------------------------------------------------------------------------
type(quaternion) elemental function pow_scal__(self,expon)
class(quaternion), intent(in) :: self
real(pReal), intent(in) :: expon
pow_scal__ = exp(log(self)*expon)
class(quaternion), intent(in) :: self
real(pReal), intent(in) :: expon
pow_scal__ = exp(log(self)*expon)
end function pow_scal__
@ -327,11 +324,11 @@ end function pow_scal__
!---------------------------------------------------------------------------------------------------
type(quaternion) elemental function pow_quat__(self,expon)
class(quaternion), intent(in) :: self
type(quaternion), intent(in) :: expon
pow_quat__ = exp(log(self)*expon)
class(quaternion), intent(in) :: self
type(quaternion), intent(in) :: expon
pow_quat__ = exp(log(self)*expon)
end function pow_quat__
@ -341,15 +338,15 @@ end function pow_quat__
!---------------------------------------------------------------------------------------------------
type(quaternion) elemental function exp__(self)
class(quaternion), intent(in) :: self
real(pReal) :: absImag
class(quaternion), intent(in) :: self
real(pReal) :: absImag
absImag = norm2([self%x, self%y, self%z])
absImag = norm2([self%x, self%y, self%z])
exp__ = exp(self%w) * [ cos(absImag), &
self%x/absImag * sin(absImag), &
self%y/absImag * sin(absImag), &
self%z/absImag * sin(absImag)]
exp__ = exp(self%w) * [ cos(absImag), &
self%x/absImag * sin(absImag), &
self%y/absImag * sin(absImag), &
self%z/absImag * sin(absImag)]
end function exp__
@ -360,16 +357,16 @@ end function exp__
!---------------------------------------------------------------------------------------------------
type(quaternion) elemental function log__(self)
class(quaternion), intent(in) :: self
real(pReal) :: absImag
class(quaternion), intent(in) :: self
real(pReal) :: absImag
absImag = norm2([self%x, self%y, self%z])
absImag = norm2([self%x, self%y, self%z])
log__ = [log(abs(self)), &
self%x/absImag * acos(self%w/abs(self)), &
self%y/absImag * acos(self%w/abs(self)), &
self%z/absImag * acos(self%w/abs(self))]
log__ = [log(abs(self)), &
self%x/absImag * acos(self%w/abs(self)), &
self%y/absImag * acos(self%w/abs(self)), &
self%z/absImag * acos(self%w/abs(self))]
end function log__
@ -378,10 +375,10 @@ end function log__
!---------------------------------------------------------------------------------------------------
real(pReal) elemental function abs__(a)
class(quaternion), intent(in) :: a
class(quaternion), intent(in) :: a
abs__ = norm2([a%w,a%x,a%y,a%z])
abs__ = norm2([a%w,a%x,a%y,a%z])
end function abs__
@ -390,10 +387,10 @@ end function abs__
!---------------------------------------------------------------------------------------------------
real(pReal) elemental function dot_product__(a,b)
class(quaternion), intent(in) :: a,b
class(quaternion), intent(in) :: a,b
dot_product__ = a%w*b%w + a%x*b%x + a%y*b%y + a%z*b%z
dot_product__ = a%w*b%w + a%x*b%x + a%y*b%y + a%z*b%z
end function dot_product__
@ -402,10 +399,10 @@ end function dot_product__
!---------------------------------------------------------------------------------------------------
type(quaternion) elemental function conjg__(a)
class(quaternion), intent(in) :: a
class(quaternion), intent(in) :: a
conjg__ = quaternion([a%w, -a%x, -a%y, -a%z])
conjg__ = quaternion([a%w, -a%x, -a%y, -a%z])
end function conjg__
@ -414,10 +411,10 @@ end function conjg__
!---------------------------------------------------------------------------------------------------
type(quaternion) elemental function quat_homomorphed(a)
class(quaternion), intent(in) :: a
class(quaternion), intent(in) :: a
quat_homomorphed = quaternion(-[a%w,a%x,a%y,a%z])
quat_homomorphed = quaternion(-[a%w,a%x,a%y,a%z])
end function quat_homomorphed
end module quaternions

16
src/results.f90
View File

@ -5,6 +5,9 @@
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!--------------------------------------------------------------------------------------------------
module results
use DAMASK_interface
use rotations
use numerics
use HDF5_utilities
#ifdef PETSc
use PETSC
@ -55,8 +58,6 @@ module results
contains
subroutine results_init
use DAMASK_interface, only: &
getSolverJobName
character(len=pStringLen) :: commandLine
@ -83,9 +84,6 @@ end subroutine results_init
!> @brief opens the results file to append data
!--------------------------------------------------------------------------------------------------
subroutine results_openJobFile
use DAMASK_interface, only: &
getSolverJobName
resultsFile = HDF5_openFile(trim(getSolverJobName())//'.hdf5','a',.true.)
@ -396,8 +394,6 @@ end subroutine results_writeTensorDataset_int
!> @brief stores a scalar dataset in a group
!--------------------------------------------------------------------------------------------------
subroutine results_writeScalarDataset_rotation(group,dataset,label,description,lattice_structure)
use rotations, only: &
rotation
character(len=*), intent(in) :: label,group,description
character(len=*), intent(in), optional :: lattice_structure
@ -428,9 +424,6 @@ end subroutine results_writeScalarDataset_rotation
!> @brief adds the unique mapping from spatial position and constituent ID to results
!--------------------------------------------------------------------------------------------------
subroutine results_mapping_constituent(phaseAt,memberAt,label)
use numerics, only: &
worldrank, &
worldsize
integer, dimension(:,:), intent(in) :: phaseAt !< phase section at (constituent,element)
integer, dimension(:,:,:), intent(in) :: memberAt !< phase member at (constituent,IP,element)
@ -566,9 +559,6 @@ end subroutine results_mapping_constituent
!> @brief adds the unique mapping from spatial position and constituent ID to results
!--------------------------------------------------------------------------------------------------
subroutine results_mapping_materialpoint(homogenizationAt,memberAt,label)
use numerics, only: &
worldrank, &
worldsize
integer, dimension(:), intent(in) :: homogenizationAt !< homogenization section at (element)
integer, dimension(:,:), intent(in) :: memberAt !< homogenization member at (IP,element)

86
src/rotations.f90
View File

@ -46,12 +46,15 @@
!---------------------------------------------------------------------------------------------------
module rotations
use prec, only: &
pReal
use prec
use IO
use math
use Lambert
use quaternions
implicit none
private
type, public :: rotation
type(quaternion), private :: q
contains
@ -148,8 +151,6 @@ end subroutine
!> @details: rotation is based on unit quaternion or rotation matrix (fallback)
!---------------------------------------------------------------------------------------------------
function rotVector(self,v,active)
use prec, only: &
dEq
real(pReal), dimension(3) :: rotVector
class(rotation), intent(in) :: self
@ -260,10 +261,6 @@ end function qu2om
!> @brief convert unit quaternion to Euler angles
!---------------------------------------------------------------------------------------------------
pure function qu2eu(qu) result(eu)
use prec, only: &
dEq0
use math, only: &
PI
type(quaternion), intent(in) :: qu
real(pReal), dimension(3) :: eu
@ -294,12 +291,6 @@ end function qu2eu
!> @brief convert unit quaternion to axis angle pair
!---------------------------------------------------------------------------------------------------
pure function qu2ax(qu) result(ax)
use prec, only: &
dEq0, &
dNeq0
use math, only: &
PI, &
math_clip
type(quaternion), intent(in) :: qu
real(pReal), dimension(4) :: ax
@ -324,13 +315,6 @@ end function qu2ax
!> @brief convert unit quaternion to Rodrigues vector
!---------------------------------------------------------------------------------------------------
pure function qu2ro(qu) result(ro)
use, intrinsic :: IEEE_ARITHMETIC, only: &
IEEE_value, &
IEEE_positive_inf
use prec, only: &
dEq0
use math, only: &
math_clip
type(quaternion), intent(in) :: qu
real(pReal), dimension(4) :: ro
@ -358,10 +342,6 @@ end function qu2ro
!> @brief convert unit quaternion to homochoric
!---------------------------------------------------------------------------------------------------
pure function qu2ho(qu) result(ho)
use prec, only: &
dEq0
use math, only: &
math_clip
type(quaternion), intent(in) :: qu
real(pReal), dimension(3) :: ho
@ -415,8 +395,6 @@ end function om2qu
!> @brief orientation matrix to Euler angles
!---------------------------------------------------------------------------------------------------
pure function om2eu(om) result(eu)
use math, only: &
PI
real(pReal), intent(in), dimension(3,3) :: om
real(pReal), dimension(3) :: eu
@ -441,15 +419,6 @@ end function om2eu
!> @brief convert orientation matrix to axis angle pair
!---------------------------------------------------------------------------------------------------
function om2ax(om) result(ax)
use prec, only: &
dEq0, &
cEq, &
dNeq0
use IO, only: &
IO_error
use math, only: &
math_clip, &
math_trace33
real(pReal), intent(in) :: om(3,3)
real(pReal) :: ax(4)
@ -560,8 +529,6 @@ end function eu2qu
!> @brief Euler angles to orientation matrix
!---------------------------------------------------------------------------------------------------
pure function eu2om(eu) result(om)
use prec, only: &
dEq0
real(pReal), intent(in), dimension(3) :: eu
real(pReal), dimension(3,3) :: om
@ -591,11 +558,6 @@ end function eu2om
!> @brief convert euler to axis angle
!---------------------------------------------------------------------------------------------------
pure function eu2ax(eu) result(ax)
use prec, only: &
dEq0, &
dEq
use math, only: &
PI
real(pReal), intent(in), dimension(3) :: eu
real(pReal), dimension(4) :: ax
@ -625,13 +587,6 @@ end function eu2ax
!> @brief Euler angles to Rodrigues vector
!---------------------------------------------------------------------------------------------------
pure function eu2ro(eu) result(ro)
use prec, only: &
dEq0
use, intrinsic :: IEEE_ARITHMETIC, only: &
IEEE_value, &
IEEE_positive_inf
use math, only: &
PI
real(pReal), intent(in), dimension(3) :: eu
real(pReal), dimension(4) :: ro
@ -681,8 +636,6 @@ end function eu2cu
!> @brief convert axis angle pair to quaternion
!---------------------------------------------------------------------------------------------------
pure function ax2qu(ax) result(qu)
use prec, only: &
dEq0
real(pReal), intent(in), dimension(4) :: ax
type(quaternion) :: qu
@ -755,13 +708,6 @@ end function ax2eu
!> @brief convert axis angle pair to Rodrigues vector
!---------------------------------------------------------------------------------------------------
pure function ax2ro(ax) result(ro)
use, intrinsic :: IEEE_ARITHMETIC, only: &
IEEE_value, &
IEEE_positive_inf
use prec, only: &
dEq0
use math, only: &
PI
real(pReal), intent(in), dimension(4) :: ax
real(pReal), dimension(4) :: ro
@ -858,12 +804,6 @@ end function ro2eu
!> @brief convert Rodrigues vector to axis angle pair
!---------------------------------------------------------------------------------------------------
pure function ro2ax(ro) result(ax)
use, intrinsic :: IEEE_ARITHMETIC, only: &
IEEE_is_finite
use prec, only: &
dEq0
use math, only: &
PI
real(pReal), intent(in), dimension(4) :: ro
real(pReal), dimension(4) :: ax
@ -890,12 +830,6 @@ end function ro2ax
!> @brief convert Rodrigues vector to homochoric
!---------------------------------------------------------------------------------------------------
pure function ro2ho(ro) result(ho)
use, intrinsic :: IEEE_ARITHMETIC, only: &
IEEE_is_finite
use prec, only: &
dEq0
use math, only: &
PI
real(pReal), intent(in), dimension(4) :: ro
real(pReal), dimension(3) :: ho
@ -973,8 +907,6 @@ end function ho2eu
!> @brief convert homochoric to axis angle pair
!---------------------------------------------------------------------------------------------------
pure function ho2ax(ho) result(ax)
use prec, only: &
dEq0
real(pReal), intent(in), dimension(3) :: ho
real(pReal), dimension(4) :: ax
@ -1029,13 +961,11 @@ end function ho2ro
!> @brief convert homochoric to cubochoric
!---------------------------------------------------------------------------------------------------
function ho2cu(ho) result(cu)
use Lambert, only: &
LambertBallToCube
real(pReal), intent(in), dimension(3) :: ho
real(pReal), dimension(3) :: cu
cu = LambertBallToCube(ho)
cu = Lambert_BallToCube(ho)
end function ho2cu
@ -1115,13 +1045,11 @@ end function cu2ro
!> @brief convert cubochoric to homochoric
!---------------------------------------------------------------------------------------------------
function cu2ho(cu) result(ho)
use Lambert, only: &
LambertCubeToBall
real(pReal), intent(in), dimension(3) :: cu
real(pReal), dimension(3) :: ho
ho = LambertCubeToBall(cu)
ho = Lambert_CubeToBall(cu)
end function cu2ho

465
src/source_damage_anisoBrittle.f90
View File

@ -5,55 +5,62 @@
!> @details to be done
!--------------------------------------------------------------------------------------------------
module source_damage_anisoBrittle
use prec
use prec
use debug
use IO
use math
use material
use config
use lattice
implicit none
private
integer, dimension(:), allocatable, public, protected :: &
source_damage_anisoBrittle_offset, & !< which source is my current source mechanism?
source_damage_anisoBrittle_instance !< instance of source mechanism
implicit none
private
integer, dimension(:,:), allocatable, target, public :: &
source_damage_anisoBrittle_sizePostResult !< size of each post result output
integer, dimension(:), allocatable, public, protected :: &
source_damage_anisoBrittle_offset, & !< which source is my current source mechanism?
source_damage_anisoBrittle_instance !< instance of source mechanism
character(len=64), dimension(:,:), allocatable, target, public :: &
source_damage_anisoBrittle_output !< name of each post result output
integer, dimension(:,:), allocatable, private :: &
source_damage_anisoBrittle_Ncleavage !< number of cleavage systems per family
integer, dimension(:,:), allocatable, target, public :: &
source_damage_anisoBrittle_sizePostResult !< size of each post result output
enum, bind(c)
enumerator :: undefined_ID, &
damage_drivingforce_ID
end enum
character(len=64), dimension(:,:), allocatable, target, public :: &
source_damage_anisoBrittle_output !< name of each post result output
integer, dimension(:,:), allocatable :: &
source_damage_anisoBrittle_Ncleavage !< number of cleavage systems per family
enum, bind(c)
enumerator :: undefined_ID, &
damage_drivingforce_ID
end enum
type, private :: tParameters !< container type for internal constitutive parameters
real(pReal) :: &
aTol, &
sdot_0, &
N
real(pReal), dimension(:), allocatable :: &
critDisp, &
critLoad
real(pReal), dimension(:,:,:,:), allocatable :: &
cleavage_systems
integer :: &
totalNcleavage
integer, dimension(:), allocatable :: &
Ncleavage
integer(kind(undefined_ID)), allocatable, dimension(:) :: &
outputID !< ID of each post result output
end type tParameters
type :: tParameters !< container type for internal constitutive parameters
real(pReal) :: &
aTol, &
sdot_0, &
N
real(pReal), dimension(:), allocatable :: &
critDisp, &
critLoad
real(pReal), dimension(:,:,:,:), allocatable :: &
cleavage_systems
integer :: &
totalNcleavage
integer, dimension(:), allocatable :: &
Ncleavage
integer(kind(undefined_ID)), allocatable, dimension(:) :: &
outputID !< ID of each post result output
end type tParameters
type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance)
type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance)
public :: &
source_damage_anisoBrittle_init, &
source_damage_anisoBrittle_dotState, &
source_damage_anisobrittle_getRateAndItsTangent, &
source_damage_anisoBrittle_postResults
public :: &
source_damage_anisoBrittle_init, &
source_damage_anisoBrittle_dotState, &
source_damage_anisobrittle_getRateAndItsTangent, &
source_damage_anisoBrittle_postResults
contains
@ -63,266 +70,230 @@ contains
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
subroutine source_damage_anisoBrittle_init
use debug, only: &
debug_level,&
debug_constitutive,&
debug_levelBasic
use IO, only: &
IO_error
use math, only: &
math_expand
use material, only: &
material_allocateSourceState, &
phase_source, &
phase_Nsources, &
phase_Noutput, &
SOURCE_damage_anisoBrittle_label, &
SOURCE_damage_anisoBrittle_ID, &
material_phase, &
sourceState
use config, only: &
config_phase, &
material_Nphase
use lattice, only: &
lattice_SchmidMatrix_cleavage, &
lattice_maxNcleavageFamily
integer :: Ninstance,phase,instance,source,sourceOffset
integer :: NofMyPhase,p ,i
integer, dimension(0), parameter :: emptyIntArray = [integer::]
character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::]
integer(kind(undefined_ID)) :: &
outputID
integer :: Ninstance,phase,instance,source,sourceOffset
integer :: NofMyPhase,p ,i
integer, dimension(0), parameter :: emptyIntArray = [integer::]
character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::]
integer(kind(undefined_ID)) :: &
outputID
character(len=pStringLen) :: &
extmsg = ''
character(len=65536), dimension(:), allocatable :: &
outputs
character(len=pStringLen) :: &
extmsg = ''
character(len=65536), dimension(:), allocatable :: &
outputs
write(6,'(/,a)') ' <<<+- source_'//SOURCE_DAMAGE_ANISOBRITTLE_LABEL//' init -+>>>'
write(6,'(/,a)') ' <<<+- source_'//SOURCE_DAMAGE_ANISOBRITTLE_LABEL//' init -+>>>'
Ninstance = int(count(phase_source == SOURCE_damage_anisoBrittle_ID))
if (Ninstance == 0) return
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) &
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
allocate(source_damage_anisoBrittle_offset(material_Nphase), source=0)
allocate(source_damage_anisoBrittle_instance(material_Nphase), source=0)
do phase = 1, material_Nphase
source_damage_anisoBrittle_instance(phase) = count(phase_source(:,1:phase) == source_damage_anisoBrittle_ID)
do source = 1, phase_Nsources(phase)
if (phase_source(source,phase) == source_damage_anisoBrittle_ID) &
source_damage_anisoBrittle_offset(phase) = source
enddo
enddo
allocate(source_damage_anisoBrittle_sizePostResult(maxval(phase_Noutput),Ninstance), source=0)
allocate(source_damage_anisoBrittle_output(maxval(phase_Noutput),Ninstance))
source_damage_anisoBrittle_output = ''
Ninstance = count(phase_source == SOURCE_damage_anisoBrittle_ID)
if (Ninstance == 0) return
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) &
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
allocate(source_damage_anisoBrittle_offset(material_Nphase), source=0)
allocate(source_damage_anisoBrittle_instance(material_Nphase), source=0)
do phase = 1, material_Nphase
source_damage_anisoBrittle_instance(phase) = count(phase_source(:,1:phase) == source_damage_anisoBrittle_ID)
do source = 1, phase_Nsources(phase)
if (phase_source(source,phase) == source_damage_anisoBrittle_ID) &
source_damage_anisoBrittle_offset(phase) = source
enddo
enddo
allocate(source_damage_anisoBrittle_sizePostResult(maxval(phase_Noutput),Ninstance), source=0)
allocate(source_damage_anisoBrittle_output(maxval(phase_Noutput),Ninstance))
source_damage_anisoBrittle_output = ''
allocate(source_damage_anisoBrittle_Ncleavage(lattice_maxNcleavageFamily,Ninstance), source=0)
allocate(source_damage_anisoBrittle_Ncleavage(lattice_maxNcleavageFamily,Ninstance), source=0)
allocate(param(Ninstance))
do p=1, size(config_phase)
if (all(phase_source(:,p) /= SOURCE_DAMAGE_ANISOBRITTLE_ID)) cycle
associate(prm => param(source_damage_anisoBrittle_instance(p)), &
config => config_phase(p))
prm%aTol = config%getFloat('anisobrittle_atol',defaultVal = 1.0e-3_pReal)
allocate(param(Ninstance))
do p=1, size(config_phase)
if (all(phase_source(:,p) /= SOURCE_DAMAGE_ANISOBRITTLE_ID)) cycle
associate(prm => param(source_damage_anisoBrittle_instance(p)), &
config => config_phase(p))
prm%aTol = config%getFloat('anisobrittle_atol',defaultVal = 1.0e-3_pReal)
prm%N = config%getFloat('anisobrittle_ratesensitivity')
prm%sdot_0 = config%getFloat('anisobrittle_sdot0')
! sanity checks
if (prm%aTol < 0.0_pReal) extmsg = trim(extmsg)//' anisobrittle_atol'
if (prm%N <= 0.0_pReal) extmsg = trim(extmsg)//' anisobrittle_ratesensitivity'
if (prm%sdot_0 <= 0.0_pReal) extmsg = trim(extmsg)//' anisobrittle_sdot0'
prm%Ncleavage = config%getInts('ncleavage',defaultVal=emptyIntArray)
prm%N = config%getFloat('anisobrittle_ratesensitivity')
prm%sdot_0 = config%getFloat('anisobrittle_sdot0')
! sanity checks
if (prm%aTol < 0.0_pReal) extmsg = trim(extmsg)//' anisobrittle_atol'
if (prm%N <= 0.0_pReal) extmsg = trim(extmsg)//' anisobrittle_ratesensitivity'
if (prm%sdot_0 <= 0.0_pReal) extmsg = trim(extmsg)//' anisobrittle_sdot0'
prm%Ncleavage = config%getInts('ncleavage',defaultVal=emptyIntArray)
prm%critDisp = config%getFloats('anisobrittle_criticaldisplacement',requiredSize=size(prm%Ncleavage))
prm%critLoad = config%getFloats('anisobrittle_criticalload', requiredSize=size(prm%Ncleavage))
prm%cleavage_systems = lattice_SchmidMatrix_cleavage (prm%Ncleavage,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal))
prm%critDisp = config%getFloats('anisobrittle_criticaldisplacement',requiredSize=size(prm%Ncleavage))
prm%critLoad = config%getFloats('anisobrittle_criticalload', requiredSize=size(prm%Ncleavage))
prm%cleavage_systems = lattice_SchmidMatrix_cleavage (prm%Ncleavage,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal))
! expand: family => system
prm%critDisp = math_expand(prm%critDisp, prm%Ncleavage)
prm%critLoad = math_expand(prm%critLoad, prm%Ncleavage)
if (any(prm%critLoad < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_criticalload'
if (any(prm%critDisp < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_criticaldisplacement'
! expand: family => system
prm%critDisp = math_expand(prm%critDisp, prm%Ncleavage)
prm%critLoad = math_expand(prm%critLoad, prm%Ncleavage)
if (any(prm%critLoad < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_criticalload'
if (any(prm%critDisp < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_criticaldisplacement'
!--------------------------------------------------------------------------------------------------
! exit if any parameter is out of range
if (extmsg /= '') &
call IO_error(211,ext_msg=trim(extmsg)//'('//SOURCE_DAMAGE_ANISOBRITTLE_LABEL//')')
if (extmsg /= '') &
call IO_error(211,ext_msg=trim(extmsg)//'('//SOURCE_DAMAGE_ANISOBRITTLE_LABEL//')')
!--------------------------------------------------------------------------------------------------
! output pararameters
outputs = config%getStrings('(output)',defaultVal=emptyStringArray)
allocate(prm%outputID(0))
do i=1, size(outputs)
outputID = undefined_ID
select case(outputs(i))
case ('anisobrittle_drivingforce')
source_damage_anisoBrittle_sizePostResult(i,source_damage_anisoBrittle_instance(p)) = 1
source_damage_anisoBrittle_output(i,source_damage_anisoBrittle_instance(p)) = outputs(i)
prm%outputID = [prm%outputID, damage_drivingforce_ID]
outputs = config%getStrings('(output)',defaultVal=emptyStringArray)
allocate(prm%outputID(0))
do i=1, size(outputs)
outputID = undefined_ID
select case(outputs(i))
case ('anisobrittle_drivingforce')
source_damage_anisoBrittle_sizePostResult(i,source_damage_anisoBrittle_instance(p)) = 1
source_damage_anisoBrittle_output(i,source_damage_anisoBrittle_instance(p)) = outputs(i)
prm%outputID = [prm%outputID, damage_drivingforce_ID]
end select
end select
enddo
enddo
end associate
phase = p
NofMyPhase=count(material_phase==phase)
instance = source_damage_anisoBrittle_instance(phase)
sourceOffset = source_damage_anisoBrittle_offset(phase)
end associate
phase = p
NofMyPhase=count(material_phase==phase)
instance = source_damage_anisoBrittle_instance(phase)
sourceOffset = source_damage_anisoBrittle_offset(phase)
call material_allocateSourceState(phase,sourceOffset,NofMyPhase,1,1,0)
sourceState(phase)%p(sourceOffset)%sizePostResults = sum(source_damage_anisoBrittle_sizePostResult(:,instance))
sourceState(phase)%p(sourceOffset)%aTolState=param(instance)%aTol
call material_allocateSourceState(phase,sourceOffset,NofMyPhase,1,1,0)
sourceState(phase)%p(sourceOffset)%sizePostResults = sum(source_damage_anisoBrittle_sizePostResult(:,instance))
sourceState(phase)%p(sourceOffset)%aTolState=param(instance)%aTol
source_damage_anisoBrittle_Ncleavage(1:size(param(instance)%Ncleavage),instance) = param(instance)%Ncleavage
enddo
source_damage_anisoBrittle_Ncleavage(1:size(param(instance)%Ncleavage),instance) = param(instance)%Ncleavage
enddo
end subroutine source_damage_anisoBrittle_init
!--------------------------------------------------------------------------------------------------
!> @brief calculates derived quantities from state
!--------------------------------------------------------------------------------------------------
subroutine source_damage_anisoBrittle_dotState(S, ipc, ip, el)
use math, only: &
math_mul33xx33
use material, only: &
phaseAt, phasememberAt, &
sourceState, &
material_homogenizationAt, &
damage, &
damageMapping
use lattice, only: &
lattice_Scleavage, &
lattice_maxNcleavageFamily, &
lattice_NcleavageSystem
integer, intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
real(pReal), intent(in), dimension(3,3) :: &
S
integer :: &
phase, &
constituent, &
instance, &
sourceOffset, &
damageOffset, &
homog, &
f, i, index_myFamily, index
real(pReal) :: &
traction_d, traction_t, traction_n, traction_crit
integer, intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
real(pReal), intent(in), dimension(3,3) :: &
S
integer :: &
phase, &
constituent, &
instance, &
sourceOffset, &
damageOffset, &
homog, &
f, i, index_myFamily, index
real(pReal) :: &
traction_d, traction_t, traction_n, traction_crit
phase = phaseAt(ipc,ip,el)
constituent = phasememberAt(ipc,ip,el)
instance = source_damage_anisoBrittle_instance(phase)
sourceOffset = source_damage_anisoBrittle_offset(phase)
homog = material_homogenizationAt(el)
damageOffset = damageMapping(homog)%p(ip,el)
sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = 0.0_pReal
index = 1
do f = 1,lattice_maxNcleavageFamily
index_myFamily = sum(lattice_NcleavageSystem(1:f-1,phase)) ! at which index starts my family
do i = 1,source_damage_anisoBrittle_Ncleavage(f,instance) ! process each (active) cleavage system in family
phase = phaseAt(ipc,ip,el)
constituent = phasememberAt(ipc,ip,el)
instance = source_damage_anisoBrittle_instance(phase)
sourceOffset = source_damage_anisoBrittle_offset(phase)
homog = material_homogenizationAt(el)
damageOffset = damageMapping(homog)%p(ip,el)
sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = 0.0_pReal
index = 1
do f = 1,lattice_maxNcleavageFamily
index_myFamily = sum(lattice_NcleavageSystem(1:f-1,phase)) ! at which index starts my family
do i = 1,source_damage_anisoBrittle_Ncleavage(f,instance) ! process each (active) cleavage system in family
traction_d = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,1,index_myFamily+i,phase))
traction_t = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,2,index_myFamily+i,phase))
traction_n = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,3,index_myFamily+i,phase))
traction_crit = param(instance)%critLoad(index)* &
damage(homog)%p(damageOffset)*damage(homog)%p(damageOffset)
traction_d = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,1,index_myFamily+i,phase))
traction_t = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,2,index_myFamily+i,phase))
traction_n = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,3,index_myFamily+i,phase))
traction_crit = param(instance)%critLoad(index)* &
damage(homog)%p(damageOffset)*damage(homog)%p(damageOffset)
sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = &
sourceState(phase)%p(sourceOffset)%dotState(1,constituent) + &
param(instance)%sdot_0* &
((max(0.0_pReal, abs(traction_d) - traction_crit)/traction_crit)**param(instance)%N + &
(max(0.0_pReal, abs(traction_t) - traction_crit)/traction_crit)**param(instance)%N + &
(max(0.0_pReal, abs(traction_n) - traction_crit)/traction_crit)**param(instance)%N)/ &
param(instance)%critDisp(index)
sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = &
sourceState(phase)%p(sourceOffset)%dotState(1,constituent) + &
param(instance)%sdot_0* &
((max(0.0_pReal, abs(traction_d) - traction_crit)/traction_crit)**param(instance)%N + &
(max(0.0_pReal, abs(traction_t) - traction_crit)/traction_crit)**param(instance)%N + &
(max(0.0_pReal, abs(traction_n) - traction_crit)/traction_crit)**param(instance)%N)/ &
param(instance)%critDisp(index)
index = index + 1
enddo
enddo
index = index + 1
enddo
enddo
end subroutine source_damage_anisoBrittle_dotState
!--------------------------------------------------------------------------------------------------
!> @brief returns local part of nonlocal damage driving force
!--------------------------------------------------------------------------------------------------
subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
use material, only: &
sourceState
integer, intent(in) :: &
phase, &
constituent
real(pReal), intent(in) :: &
phi
real(pReal), intent(out) :: &
localphiDot, &
dLocalphiDot_dPhi
integer :: &
sourceOffset
integer, intent(in) :: &
phase, &
constituent
real(pReal), intent(in) :: &
phi
real(pReal), intent(out) :: &
localphiDot, &
dLocalphiDot_dPhi
integer :: &
sourceOffset
sourceOffset = source_damage_anisoBrittle_offset(phase)
localphiDot = 1.0_pReal &
- sourceState(phase)%p(sourceOffset)%state(1,constituent)*phi
dLocalphiDot_dPhi = -sourceState(phase)%p(sourceOffset)%state(1,constituent)
sourceOffset = source_damage_anisoBrittle_offset(phase)
localphiDot = 1.0_pReal &
- sourceState(phase)%p(sourceOffset)%state(1,constituent)*phi
dLocalphiDot_dPhi = -sourceState(phase)%p(sourceOffset)%state(1,constituent)
end subroutine source_damage_anisobrittle_getRateAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief return array of local damage results
!--------------------------------------------------------------------------------------------------
function source_damage_anisoBrittle_postResults(phase, constituent)
use material, only: &
sourceState
integer, intent(in) :: &
phase, &
constituent
real(pReal), dimension(sum(source_damage_anisoBrittle_sizePostResult(:, &
source_damage_anisoBrittle_instance(phase)))) :: &
source_damage_anisoBrittle_postResults
integer, intent(in) :: &
phase, &
constituent
integer :: &
instance, sourceOffset, o, c
instance = source_damage_anisoBrittle_instance(phase)
sourceOffset = source_damage_anisoBrittle_offset(phase)
real(pReal), dimension(sum(source_damage_anisoBrittle_sizePostResult(:, &
source_damage_anisoBrittle_instance(phase)))) :: &
source_damage_anisoBrittle_postResults
c = 0
integer :: &
instance, sourceOffset, o, c
instance = source_damage_anisoBrittle_instance(phase)
sourceOffset = source_damage_anisoBrittle_offset(phase)
do o = 1,size(param(instance)%outputID)
select case(param(instance)%outputID(o))
case (damage_drivingforce_ID)
source_damage_anisoBrittle_postResults(c+1) = &
sourceState(phase)%p(sourceOffset)%state(1,constituent)
c = c + 1
c = 0
end select
enddo
do o = 1,size(param(instance)%outputID)
select case(param(instance)%outputID(o))
case (damage_drivingforce_ID)
source_damage_anisoBrittle_postResults(c+1) = &
sourceState(phase)%p(sourceOffset)%state(1,constituent)
c = c + 1
end select
enddo
end function source_damage_anisoBrittle_postResults
end module source_damage_anisoBrittle

2
src/source_damage_isoBrittle.f90
View File

@ -84,7 +84,7 @@ subroutine source_damage_isoBrittle_init
write(6,'(/,a)') ' <<<+- source_'//SOURCE_DAMAGE_ISOBRITTLE_LABEL//' init -+>>>'
Ninstance = int(count(phase_source == SOURCE_damage_isoBrittle_ID))
Ninstance = count(phase_source == SOURCE_damage_isoBrittle_ID)
if (Ninstance == 0) return
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) &

34
src/source_thermal_dissipation.f90
View File

@ -5,27 +5,30 @@
!> @details to be done
!--------------------------------------------------------------------------------------------------
module source_thermal_dissipation
use prec, only: &
pReal
use prec
use debug
use material
use config
implicit none
private
integer, dimension(:), allocatable, public, protected :: &
source_thermal_dissipation_offset, & !< which source is my current thermal dissipation mechanism?
source_thermal_dissipation_instance !< instance of thermal dissipation source mechanism
source_thermal_dissipation_offset, & !< which source is my current thermal dissipation mechanism?
source_thermal_dissipation_instance !< instance of thermal dissipation source mechanism
integer, dimension(:,:), allocatable, target, public :: &
source_thermal_dissipation_sizePostResult !< size of each post result output
source_thermal_dissipation_sizePostResult !< size of each post result output
character(len=64), dimension(:,:), allocatable, target, public :: &
source_thermal_dissipation_output !< name of each post result output
source_thermal_dissipation_output !< name of each post result output
type, private :: tParameters !< container type for internal constitutive parameters
type :: tParameters !< container type for internal constitutive parameters
real(pReal) :: &
kappa
end type tParameters
type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance)
type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance)
public :: &
@ -40,21 +43,6 @@ contains
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
subroutine source_thermal_dissipation_init
use debug, only: &
debug_level,&
debug_constitutive,&
debug_levelBasic
use material, only: &
material_allocateSourceState, &
phase_source, &
phase_Nsources, &
phase_Noutput, &
SOURCE_thermal_dissipation_label, &
SOURCE_thermal_dissipation_ID, &
material_phase
use config, only: &
config_phase, &
material_Nphase
integer :: Ninstance,instance,source,sourceOffset
integer :: NofMyPhase,p

31
src/source_thermal_externalheat.f90
View File

@ -5,11 +5,14 @@
!> @brief material subroutine for variable heat source
!--------------------------------------------------------------------------------------------------
module source_thermal_externalheat
use prec, only: &
pReal
use prec
use debug
use material
use config
implicit none
private
integer, dimension(:), allocatable, public, protected :: &
source_thermal_externalheat_offset, & !< which source is my current thermal dissipation mechanism?
source_thermal_externalheat_instance !< instance of thermal dissipation source mechanism
@ -23,7 +26,7 @@ module source_thermal_externalheat
integer, dimension(:), allocatable, target, public :: &
source_thermal_externalheat_Noutput !< number of outputs per instance of this source
type, private :: tParameters !< container type for internal constitutive parameters
type :: tParameters !< container type for internal constitutive parameters
real(pReal), dimension(:), allocatable :: &
time, &
heat_rate
@ -31,7 +34,7 @@ module source_thermal_externalheat
nIntervals
end type tParameters
type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance)
type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance)
public :: &
@ -47,22 +50,6 @@ contains
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
subroutine source_thermal_externalheat_init
use debug, only: &
debug_level,&
debug_constitutive,&
debug_levelBasic
use material, only: &
material_allocateSourceState, &
material_phase, &
phase_source, &
phase_Nsources, &
phase_Noutput, &
SOURCE_thermal_externalheat_label, &
SOURCE_thermal_externalheat_ID
use config, only: &
config_phase, &
material_Nphase
integer :: maxNinstance,instance,source,sourceOffset,NofMyPhase,p
@ -116,8 +103,6 @@ end subroutine source_thermal_externalheat_init
!> @details state only contains current time to linearly interpolate given heat powers
!--------------------------------------------------------------------------------------------------
subroutine source_thermal_externalheat_dotState(phase, of)
use material, only: &
sourceState
integer, intent(in) :: &
phase, &
@ -135,8 +120,6 @@ end subroutine source_thermal_externalheat_dotState
!> @brief returns local heat generation rate
!--------------------------------------------------------------------------------------------------
subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, dTDot_dT, phase, of)
use material, only: &
sourceState
integer, intent(in) :: &
phase, &

78
src/thermal_adiabatic.f90
View File

@ -3,9 +3,16 @@
!> @brief material subroutine for adiabatic temperature evolution
!--------------------------------------------------------------------------------------------------
module thermal_adiabatic
use prec, only: &
pReal
use prec
use config
use numerics
use material
use source_thermal_dissipation
use source_thermal_externalheat
use crystallite
use lattice
use mesh
implicit none
private
@ -21,7 +28,7 @@ module thermal_adiabatic
enumerator :: undefined_ID, &
temperature_ID
end enum
integer(kind(undefined_ID)), dimension(:,:), allocatable, private :: &
integer(kind(undefined_ID)), dimension(:,:), allocatable :: &
thermal_adiabatic_outputID !< ID of each post result output
@ -41,21 +48,6 @@ contains
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
subroutine thermal_adiabatic_init
use material, only: &
thermal_type, &
thermal_typeInstance, &
homogenization_Noutput, &
THERMAL_ADIABATIC_label, &
THERMAL_adiabatic_ID, &
material_homogenizationAt, &
mappingHomogenization, &
thermalState, &
thermalMapping, &
thermal_initialT, &
temperature, &
temperatureRate
use config, only: &
config_homogenization
integer :: maxNinstance,section,instance,i,sizeState,NofMyHomog
character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::]
@ -112,16 +104,6 @@ end subroutine thermal_adiabatic_init
!> @brief calculates adiabatic change in temperature based on local heat generation model
!--------------------------------------------------------------------------------------------------
function thermal_adiabatic_updateState(subdt, ip, el)
use numerics, only: &
err_thermal_tolAbs, &
err_thermal_tolRel
use material, only: &
material_homogenizationAt, &
mappingHomogenization, &
thermalState, &
temperature, &
temperatureRate, &
thermalMapping
integer, intent(in) :: &
ip, & !< integration point number
@ -156,28 +138,11 @@ function thermal_adiabatic_updateState(subdt, ip, el)
end function thermal_adiabatic_updateState
!--------------------------------------------------------------------------------------------------
!> @brief returns heat generation rate
!--------------------------------------------------------------------------------------------------
subroutine thermal_adiabatic_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el)
use material, only: &
homogenization_Ngrains, &
material_homogenizationAt, &
mappingHomogenization, &
phaseAt, &
phasememberAt, &
thermal_typeInstance, &
phase_Nsources, &
phase_source, &
SOURCE_thermal_dissipation_ID, &
SOURCE_thermal_externalheat_ID
use source_thermal_dissipation, only: &
source_thermal_dissipation_getRateAndItsTangent
use source_thermal_externalheat, only: &
source_thermal_externalheat_getRateAndItsTangent
use crystallite, only: &
crystallite_S, &
crystallite_Lp
integer, intent(in) :: &
ip, & !< integration point number
@ -230,18 +195,12 @@ subroutine thermal_adiabatic_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el)
dTdot_dT = dTdot_dT/real(homogenization_Ngrains(homog),pReal)
end subroutine thermal_adiabatic_getSourceAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief returns homogenized specific heat capacity
!--------------------------------------------------------------------------------------------------
function thermal_adiabatic_getSpecificHeat(ip,el)
use lattice, only: &
lattice_specificHeat
use material, only: &
homogenization_Ngrains, &
material_phase
use mesh, only: &
mesh_element
integer, intent(in) :: &
ip, & !< integration point number
@ -270,13 +229,6 @@ end function thermal_adiabatic_getSpecificHeat
!> @brief returns homogenized mass density
!--------------------------------------------------------------------------------------------------
function thermal_adiabatic_getMassDensity(ip,el)
use lattice, only: &
lattice_massDensity
use material, only: &
homogenization_Ngrains, &
material_phase
use mesh, only: &
mesh_element
integer, intent(in) :: &
ip, & !< integration point number
@ -304,8 +256,6 @@ end function thermal_adiabatic_getMassDensity
!> @brief return array of thermal results
!--------------------------------------------------------------------------------------------------
function thermal_adiabatic_postResults(homog,instance,of) result(postResults)
use material, only: &
temperature
integer, intent(in) :: &
homog, &