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

Merge branch 'development' into geometry-class

This commit is contained in:
Martin Diehl 2019-05-31 16:58:16 +02:00
parent bf27505faa 236b5c8bde
commit e0971a9b05
40 changed files with 1228 additions and 2568 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
VERSION
View File

@ -1 +1 @@
v2.0.3-332-g5abcca50 v2.0.3-367-g70428155

2
cmake/Compiler-Intel.cmake
View File

@ -32,6 +32,8 @@
# disables warnings ... # disables warnings ...
set (COMPILE_FLAGS "${COMPILE_FLAGS} 5268") set (COMPILE_FLAGS "${COMPILE_FLAGS} 5268")
# ... the text exceeds right hand column allowed on the line (we have only comments there) # ... 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") set (COMPILE_FLAGS "${COMPILE_FLAGS} -warn")
# enables warnings ... # enables warnings ...

1
examples/SpectralMethod/Polycrystal/material.config
View File

@ -11,7 +11,6 @@ mech none
[almostAll] [almostAll]
(output) phase (output) phase
(output) texture (output) texture
(output) volume
(output) orientation # quaternion (output) orientation # quaternion
(output) grainrotation # deviation from initial orientation as axis (1-3) and angle in degree (4) (output) grainrotation # deviation from initial orientation as axis (1-3) and angle in degree (4)
(output) f # deformation gradient tensor; synonyms: "defgrad" (output) f # deformation gradient tensor; synonyms: "defgrad"

7
src/DAMASK_interface.f90
View File

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

10
src/HDF5_utilities.f90
View File

@ -8,6 +8,8 @@ module HDF5_utilities
use prec use prec
use IO use IO
use HDF5 use HDF5
use rotations
use numerics
#ifdef PETSc #ifdef PETSc
use PETSC use PETSC
#endif #endif
@ -1676,8 +1678,6 @@ end subroutine HDF5_write_int7
! ToDo: We could optionally write out other representations (axis angle, euler, ...) ! ToDo: We could optionally write out other representations (axis angle, euler, ...)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine HDF5_write_rotation(loc_id,dataset,datasetName,parallel) subroutine HDF5_write_rotation(loc_id,dataset,datasetName,parallel)
use rotations, only: &
rotation
type(rotation), intent(in), dimension(:) :: dataset type(rotation), intent(in), dimension(:) :: dataset
integer(HID_T), intent(in) :: loc_id !< file or group handle 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, & subroutine initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, &
myStart, globalShape, & myStart, globalShape, &
loc_id,localShape,datasetName,parallel) loc_id,localShape,datasetName,parallel)
use numerics, only: &
worldrank, &
worldsize
integer(HID_T), intent(in) :: loc_id !< file or group handle integer(HID_T), intent(in) :: loc_id !< file or group handle
character(len=*), intent(in) :: datasetName !< name of the dataset in the file 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, & subroutine initialize_write(dset_id, filespace_id, memspace_id, plist_id, &
myStart, totalShape, & myStart, totalShape, &
loc_id,myShape,datasetName,datatype,parallel) loc_id,myShape,datasetName,datatype,parallel)
use numerics, only: &
worldrank, &
worldsize
integer(HID_T), intent(in) :: loc_id !< file or group handle integer(HID_T), intent(in) :: loc_id !< file or group handle
character(len=*), intent(in) :: datasetName !< name of the dataset in the file 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). !> Modeling and Simulations in Materials Science and Engineering 22, 075013 (2014).
!-------------------------------------------------------------------------- !--------------------------------------------------------------------------
module Lambert module Lambert
use prec
use math use math
implicit none implicit none
private private
real(pReal), parameter, private :: &
real(pReal), parameter :: &
SPI = sqrt(PI), & SPI = sqrt(PI), &
PREF = sqrt(6.0_pReal/PI), & PREF = sqrt(6.0_pReal/PI), &
A = PI**(5.0_pReal/6.0_pReal)/6.0_pReal**(1.0_pReal/6.0_pReal), & 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 PREK = R1 * 2.0_pReal**(1.0_pReal/4.0_pReal)/BETA
public :: & public :: &
LambertCubeToBall, & Lambert_CubeToBall, &
LambertBallToCube Lambert_BallToCube
private :: &
GetPyramidOrder
contains contains
@ -68,7 +68,7 @@ contains
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief map from 3D cubic grid to 3D ball !> @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), intent(in), dimension(3) :: cube
real(pReal), dimension(3) :: ball, LamXYZ, XYZ real(pReal), dimension(3) :: ball, LamXYZ, XYZ
@ -116,7 +116,7 @@ function LambertCubeToBall(cube) result(ball)
endif center 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 !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief map from 3D ball to 3D cubic grid !> @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), intent(in), dimension(3) :: xyz
real(pReal), dimension(3) :: cube, xyz1, xyz3 real(pReal), dimension(3) :: cube, xyz1, xyz3
@ -170,7 +170,7 @@ pure function LambertBallToCube(xyz) result(cube)
endif center endif center
end function LambertBallToCube end function Lambert_BallToCube
!-------------------------------------------------------------------------- !--------------------------------------------------------------------------

1
src/commercialFEM_fileList.f90
View File

@ -14,6 +14,7 @@
#include "Lambert.f90" #include "Lambert.f90"
#include "rotations.f90" #include "rotations.f90"
#include "FEsolving.f90" #include "FEsolving.f90"
#include "geometry_plastic_nonlocal.f90"
#include "element.f90" #include "element.f90"
#include "mesh_base.f90" #include "mesh_base.f90"
#ifdef Abaqus #ifdef Abaqus

40
src/config.f90
View File

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

338
src/crystallite.f90
View File

@ -9,36 +9,43 @@
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module crystallite module crystallite
use prec, only: & use prec
pReal, & use IO
pStringLen use config
use rotations, only: & use debug
rotation use numerics
use FEsolving, only: & use rotations
FEsolving_execElem, & use math
FEsolving_execIP use mesh
use material, only: & use FEsolving
homogenization_Ngrains use material
use constitutive
use lattice
use future use future
use plastic_nonlocal
#if defined(PETSc) || defined(DAMASK_HDF5)
use HDF5_utilities
use results
#endif
implicit none implicit none
private private
character(len=64), dimension(:,:), allocatable, private :: & character(len=64), dimension(:,:), allocatable :: &
crystallite_output !< name of each post result output crystallite_output !< name of each post result output
integer, public, protected :: & integer, public, protected :: &
crystallite_maxSizePostResults !< description not available crystallite_maxSizePostResults !< description not available
integer, dimension(:), allocatable, public, protected :: & integer, dimension(:), allocatable, public, protected :: &
crystallite_sizePostResults !< description not available crystallite_sizePostResults !< description not available
integer, dimension(:,:), allocatable, private :: & integer, dimension(:,:), allocatable :: &
crystallite_sizePostResult !< description not available crystallite_sizePostResult !< description not available
real(pReal), dimension(:,:,:), allocatable, public :: & real(pReal), dimension(:,:,:), allocatable, public :: &
crystallite_dt !< requested time increment of each grain 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_subdt, & !< substepped time increment of each grain
crystallite_subFrac, & !< already calculated fraction of increment crystallite_subFrac, & !< already calculated fraction of increment
crystallite_subStep !< size of next integration step crystallite_subStep !< size of next integration step
type(rotation), dimension(:,:,:), allocatable, private :: & type(rotation), dimension(:,:,:), allocatable :: &
crystallite_orientation, & !< orientation crystallite_orientation, & !< orientation
crystallite_orientation0 !< initial orientation crystallite_orientation0 !< initial orientation
real(pReal), dimension(:,:,:,:,:), allocatable, public, protected :: & real(pReal), dimension(:,:,:,:,:), allocatable, public, protected :: &
@ -63,7 +70,7 @@ module crystallite
crystallite_Li, & !< current intermediate velocitiy grad (end of converged time step) crystallite_Li, & !< current intermediate velocitiy grad (end of converged time step)
crystallite_Li0, & !< intermediate velocitiy grad at start of FE inc crystallite_Li0, & !< intermediate velocitiy grad at start of FE inc
crystallite_partionedLi0 !< intermediate velocity grad at start of homog 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_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_invFp, & !< inverse of current plastic def grad (end of converged time step)
crystallite_subFp0,& !< plastic def grad at start of crystallite inc 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) crystallite_dPdF !< current individual dPdF per grain (end of converged time step)
logical, dimension(:,:,:), allocatable, public :: & logical, dimension(:,:,:), allocatable, public :: &
crystallite_requested !< used by upper level (homogenization) to request crystallite calculation crystallite_requested !< used by upper level (homogenization) to request crystallite calculation
logical, dimension(:,:,:), allocatable, private :: & logical, dimension(:,:,:), allocatable :: &
crystallite_converged, & !< convergence flag crystallite_converged, & !< convergence flag
crystallite_todo, & !< flag to indicate need for further computation crystallite_todo, & !< flag to indicate need for further computation
crystallite_localPlasticity !< indicates this grain to have purely local constitutive law crystallite_localPlasticity !< indicates this grain to have purely local constitutive law
@ -86,7 +93,6 @@ module crystallite
enumerator :: undefined_ID, & enumerator :: undefined_ID, &
phase_ID, & phase_ID, &
texture_ID, & texture_ID, &
volume_ID, &
orientation_ID, & orientation_ID, &
grainrotation_ID, & grainrotation_ID, &
defgrad_ID, & defgrad_ID, &
@ -101,16 +107,16 @@ module crystallite
neighboringip_ID, & neighboringip_ID, &
neighboringelement_ID neighboringelement_ID
end enum end enum
integer(kind(undefined_ID)),dimension(:,:), allocatable, private :: & integer(kind(undefined_ID)),dimension(:,:), allocatable :: &
crystallite_outputID !< ID of each post result output 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(:) :: & character(len=65536), allocatable, dimension(:) :: &
label label
end type tOutput end type tOutput
type(tOutput), allocatable, dimension(:), private :: output_constituent type(tOutput), allocatable, dimension(:) :: output_constituent
type, private :: tNumerics type :: tNumerics
integer :: & integer :: &
iJacoLpresiduum, & !< frequency of Jacobian update of residuum in Lp iJacoLpresiduum, & !< frequency of Jacobian update of residuum in Lp
nState, & !< state loop limit nState, & !< state loop limit
@ -138,15 +144,6 @@ module crystallite
crystallite_push33ToRef, & crystallite_push33ToRef, &
crystallite_postResults, & crystallite_postResults, &
crystallite_results crystallite_results
private :: &
integrateStress, &
integrateState, &
integrateStateFPI, &
integrateStateEuler, &
integrateStateAdaptiveEuler, &
integrateStateRK4, &
integrateStateRKCK45, &
stateJump
contains contains
@ -155,39 +152,6 @@ contains
!> @brief allocates and initialize per grain variables !> @brief allocates and initialize per grain variables
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine crystallite_init 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 integer, parameter :: FILEUNIT=434
logical, dimension(:,:), allocatable :: devNull logical, dimension(:,:), allocatable :: devNull
@ -321,8 +285,6 @@ subroutine crystallite_init
crystallite_outputID(o,c) = phase_ID crystallite_outputID(o,c) = phase_ID
case ('texture') outputName case ('texture') outputName
crystallite_outputID(o,c) = texture_ID crystallite_outputID(o,c) = texture_ID
case ('volume') outputName
crystallite_outputID(o,c) = volume_ID
case ('orientation') outputName case ('orientation') outputName
crystallite_outputID(o,c) = orientation_ID crystallite_outputID(o,c) = orientation_ID
case ('grainrotation') outputName case ('grainrotation') outputName
@ -371,7 +333,7 @@ subroutine crystallite_init
do r = 1,size(config_crystallite) do r = 1,size(config_crystallite)
do o = 1,crystallite_Noutput(r) do o = 1,crystallite_Noutput(r)
select case(crystallite_outputID(o,r)) select case(crystallite_outputID(o,r))
case(phase_ID,texture_ID,volume_ID) case(phase_ID,texture_ID)
mySize = 1 mySize = 1
case(orientation_ID,grainrotation_ID) case(orientation_ID,grainrotation_ID)
mySize = 4 mySize = 4
@ -478,34 +440,6 @@ end subroutine crystallite_init
!> @brief calculate stress (P) !> @brief calculate stress (P)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC) 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 logical, dimension(theMesh%elem%nIPs,theMesh%Nelems) :: crystallite_stress
real(pReal), intent(in), optional :: & real(pReal), intent(in), optional :: &
@ -746,30 +680,6 @@ end function crystallite_stress
!> @brief calculate tangent (dPdF) !> @brief calculate tangent (dPdF)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine crystallite_stressTangent 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 :: & integer :: &
c, & !< counter in integration point component loop c, & !< counter in integration point component loop
@ -910,19 +820,6 @@ end subroutine crystallite_stressTangent
!> @brief calculates orientations !> @brief calculates orientations
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine crystallite_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 & integer &
c, & !< counter in integration point component loop c, & !< counter in integration point component loop
@ -979,28 +876,6 @@ end function crystallite_push33ToRef
!> @brief return results of particular grain !> @brief return results of particular grain
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function crystallite_postResults(ipc, ip, el) 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):: & integer, intent(in):: &
el, & !< element index el, & !< element index
@ -1036,11 +911,6 @@ function crystallite_postResults(ipc, ip, el)
case (texture_ID) case (texture_ID)
mySize = 1 mySize = 1
crystallite_postResults(c+1) = real(material_texture(ipc,ip,el),pReal) ! textureID of grain crystallite_postResults(c+1) = real(material_texture(ipc,ip,el),pReal) ! textureID of grain
case (volume_ID)
mySize = 1
detF = math_det33(crystallite_partionedF(1:3,1:3,ipc,ip,el)) ! V_current = det(F) * V_reference
crystallite_postResults(c+1) = detF * mesh_ipVolume(ip,el) &
/ real(homogenization_Ngrains(mesh_element(3,el)),pReal) ! grain volume (not fraction but absolute)
case (orientation_ID) case (orientation_ID)
mySize = 4 mySize = 4
crystallite_postResults(c+1:c+mySize) = crystallite_orientation(ipc,ip,el)%asQuaternion() crystallite_postResults(c+1:c+mySize) = crystallite_orientation(ipc,ip,el)%asQuaternion()
@ -1118,16 +988,9 @@ end function crystallite_postResults
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine crystallite_results subroutine crystallite_results
#if defined(PETSc) || defined(DAMASK_HDF5) #if defined(PETSc) || defined(DAMASK_HDF5)
use lattice
use results
use HDF5_utilities
use rotations
use config, only: & use config, only: &
config_name_phase => phase_name ! anticipate logical name config_name_phase => phase_name ! anticipate logical name
use material, only: &
material_phase_plasticity_type => phase_plasticity
integer :: p,o integer :: p,o
real(pReal), allocatable, dimension(:,:,:) :: selected_tensors real(pReal), allocatable, dimension(:,:,:) :: selected_tensors
type(rotation), allocatable, dimension(:) :: selected_rotations type(rotation), allocatable, dimension(:) :: selected_rotations
@ -1267,33 +1130,6 @@ end subroutine crystallite_results
!> intermediate acceleration of the Newton-Raphson correction !> intermediate acceleration of the Newton-Raphson correction
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
logical function integrateStress(ipc,ip,el,timeFraction) 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 integer, intent(in):: el, & ! element index
ip, & ! integration point index ip, & ! integration point index
@ -1693,27 +1529,6 @@ end function integrateStress
!> using Fixed Point Iteration to adapt the stepsize !> using Fixed Point Iteration to adapt the stepsize
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine integrateStateFPI 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 :: & integer :: &
NiterationState, & !< number of iterations in state loop NiterationState, & !< number of iterations in state loop
@ -1901,8 +1716,6 @@ end subroutine integrateStateFPI
!> @brief integrate state with 1st order explicit Euler method !> @brief integrate state with 1st order explicit Euler method
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine integrateStateEuler subroutine integrateStateEuler
use material, only: &
plasticState
call update_dotState(1.0_pReal) call update_dotState(1.0_pReal)
call update_state(1.0_pReal) call update_state(1.0_pReal)
@ -1919,19 +1732,6 @@ end subroutine integrateStateEuler
!> @brief integrate stress, state with 1st order Euler method with adaptive step size !> @brief integrate stress, state with 1st order Euler method with adaptive step size
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine integrateStateAdaptiveEuler 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 :: & integer :: &
e, & ! element index in element loop e, & ! element index in element loop
@ -2025,14 +1825,6 @@ end subroutine integrateStateAdaptiveEuler
! ToDo: This is totally BROKEN: RK4dotState is never used!!! ! ToDo: This is totally BROKEN: RK4dotState is never used!!!
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine integrateStateRK4 subroutine integrateStateRK4
use mesh, only: &
mesh_element
use material, only: &
homogenization_Ngrains, &
plasticState, &
sourceState, &
phase_Nsources, &
phaseAt, phasememberAt
real(pReal), dimension(4), parameter :: & 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 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 +1884,6 @@ end subroutine integrateStateRK4
!> adaptive step size (use 5th order solution to advance = "local extrapolation") !> adaptive step size (use 5th order solution to advance = "local extrapolation")
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine integrateStateRKCK45 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 :: & real(pReal), dimension(5,5), parameter :: &
A = reshape([& A = reshape([&
@ -2287,8 +2066,6 @@ end subroutine nonlocalConvergenceCheck
!> @details: For explicitEuler, RK4 and RKCK45, adaptive Euler and FPI have their on criteria !> @details: For explicitEuler, RK4 and RKCK45, adaptive Euler and FPI have their on criteria
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine setConvergenceFlag subroutine setConvergenceFlag
use mesh, only: &
mesh_element
integer :: & integer :: &
e, & !< element index in element loop e, & !< element index in element loop
@ -2327,8 +2104,6 @@ end subroutine setConvergenceFlag
!> @brief Standard forwarding of state as state = state0 + dotState * (delta t) !> @brief Standard forwarding of state as state = state0 + dotState * (delta t)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine update_stress(timeFraction) subroutine update_stress(timeFraction)
use mesh, only: &
mesh_element
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
timeFraction timeFraction
@ -2360,8 +2135,6 @@ end subroutine update_stress
!> @brief tbd !> @brief tbd
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine update_dependentState subroutine update_dependentState
use mesh, only: &
mesh_element
use constitutive, only: & use constitutive, only: &
constitutive_dependentState => constitutive_microstructure constitutive_dependentState => constitutive_microstructure
@ -2387,13 +2160,6 @@ end subroutine update_dependentState
!> @brief Standard forwarding of state as state = state0 + dotState * (delta t) !> @brief Standard forwarding of state as state = state0 + dotState * (delta t)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine update_state(timeFraction) subroutine update_state(timeFraction)
use material, only: &
plasticState, &
sourceState, &
phase_Nsources, &
phaseAt, phasememberAt
use mesh, only: &
mesh_element
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
timeFraction timeFraction
@ -2435,17 +2201,6 @@ end subroutine update_state
!> if NaN occurs, crystallite_todo is set to FALSE. Any NaN in a nonlocal propagates to all others !> if NaN occurs, crystallite_todo is set to FALSE. Any NaN in a nonlocal propagates to all others
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine update_dotState(timeFraction) 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) :: & real(pReal), intent(in) :: &
timeFraction timeFraction
@ -2492,19 +2247,7 @@ end subroutine update_DotState
subroutine update_deltaState 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 :: & integer :: &
e, & !< element index in element loop e, & !< element index in element loop
i, & !< integration point index in ip loop i, & !< integration point index in ip loop
@ -2569,29 +2312,6 @@ end subroutine update_deltaState
!> returns true, if state jump was successfull or not needed. false indicates NaN in delta state !> returns true, if state jump was successfull or not needed. false indicates NaN in delta state
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
logical function stateJump(ipc,ip,el) 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):: & integer, intent(in):: &
el, & ! element index el, & ! element index

57
src/damage_local.f90
View File

@ -4,9 +4,13 @@
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module damage_local module damage_local
use prec use prec
use material
use numerics
use config
implicit none implicit none
private private
integer, dimension(:,:), allocatable, target, public :: & integer, dimension(:,:), allocatable, target, public :: &
damage_local_sizePostResult !< size of each post result output damage_local_sizePostResult !< size of each post result output
@ -20,23 +24,22 @@ module damage_local
enumerator :: undefined_ID, & enumerator :: undefined_ID, &
damage_ID damage_ID
end enum end enum
integer(kind(undefined_ID)), dimension(:,:), allocatable, private :: & integer(kind(undefined_ID)), dimension(:,:), allocatable :: &
damage_local_outputID !< ID of each post result output damage_local_outputID !< ID of each post result output
type, private :: tParameters type :: tParameters
integer(kind(undefined_ID)), dimension(:), allocatable :: & integer(kind(undefined_ID)), dimension(:), allocatable :: &
outputID outputID
end type tParameters end type tParameters
type(tparameters), dimension(:), allocatable, private :: & type(tparameters), dimension(:), allocatable :: &
param param
public :: & public :: &
damage_local_init, & damage_local_init, &
damage_local_updateState, & damage_local_updateState, &
damage_local_postResults damage_local_postResults
private :: &
damage_local_getSourceAndItsTangent
contains contains
@ -45,23 +48,8 @@ contains
!> @details reads in material parameters, allocates arrays, and does sanity checks !> @details reads in material parameters, allocates arrays, and does sanity checks
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine damage_local_init 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 :: sizeState
integer :: NofMyHomog, h integer :: NofMyHomog, h
integer(kind(undefined_ID)) :: & integer(kind(undefined_ID)) :: &
@ -72,7 +60,7 @@ subroutine damage_local_init
write(6,'(/,a)') ' <<<+- damage_'//DAMAGE_local_label//' 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 if (maxNinstance == 0) return
allocate(damage_local_sizePostResult (maxval(homogenization_Noutput),maxNinstance),source=0) 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 !> @brief calculates local change in damage field
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function damage_local_updateState(subdt, ip, el) 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) :: & integer, intent(in) :: &
ip, & !< integration point number ip, & !< integration point number
@ -177,17 +157,6 @@ end function damage_local_updateState
!> @brief calculates homogenized local damage driving forces !> @brief calculates homogenized local damage driving forces
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine damage_local_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el) 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: & use source_damage_isoBrittle, only: &
source_damage_isobrittle_getRateAndItsTangent source_damage_isobrittle_getRateAndItsTangent
use source_damage_isoDuctile, only: & use source_damage_isoDuctile, only: &
@ -244,15 +213,11 @@ subroutine damage_local_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el
end subroutine damage_local_getSourceAndItsTangent end subroutine damage_local_getSourceAndItsTangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief return array of damage results !> @brief return array of damage results
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function damage_local_postResults(ip,el) function damage_local_postResults(ip,el)
use material, only: &
material_homogenizationAt, &
damage_typeInstance, &
damageMapping, &
damage
integer, intent(in) :: & integer, intent(in) :: &
ip, & !< integration point ip, & !< integration point

16
src/damage_none.f90
View File

@ -3,6 +3,8 @@
!> @brief material subroutine for constant damage field !> @brief material subroutine for constant damage field
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module damage_none module damage_none
use config
use material
implicit none implicit none
private private
@ -15,18 +17,8 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief allocates all neccessary fields, reads information from material configuration file !> @brief allocates all neccessary fields, reads information from material configuration file
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine damage_none_init() 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
integer :: & integer :: &
homog, & homog, &
NofMyHomog NofMyHomog

450
src/damage_nonlocal.f90
View File

@ -4,39 +4,50 @@
!> @details to be done !> @details to be done
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module damage_nonlocal 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 implicit none
private private
integer, dimension(:,:), allocatable, target, public :: &
damage_nonlocal_sizePostResult !< size of each post result output integer, dimension(:,:), allocatable, target, public :: &
damage_nonlocal_sizePostResult !< size of each post result output
character(len=64), dimension(:,:), allocatable, target, public :: & character(len=64), dimension(:,:), allocatable, target, public :: &
damage_nonlocal_output !< name of each post result output damage_nonlocal_output !< name of each post result output
integer, dimension(:), allocatable, target, public :: & integer, dimension(:), allocatable, target, public :: &
damage_nonlocal_Noutput !< number of outputs per instance of this damage damage_nonlocal_Noutput !< number of outputs per instance of this damage
enum, bind(c) enum, bind(c)
enumerator :: undefined_ID, & enumerator :: undefined_ID, &
damage_ID damage_ID
end enum end enum
type, private :: tParameters type :: tParameters
integer(kind(undefined_ID)), dimension(:), allocatable :: & integer(kind(undefined_ID)), dimension(:), allocatable :: &
outputID outputID
end type tParameters end type tParameters
type(tparameters), dimension(:), allocatable, private :: & type(tparameters), dimension(:), allocatable :: &
param param
public :: & public :: &
damage_nonlocal_init, & damage_nonlocal_init, &
damage_nonlocal_getSourceAndItsTangent, & damage_nonlocal_getSourceAndItsTangent, &
damage_nonlocal_getDiffusion33, & damage_nonlocal_getDiffusion33, &
damage_nonlocal_getMobility, & damage_nonlocal_getMobility, &
damage_nonlocal_putNonLocalDamage, & damage_nonlocal_putNonLocalDamage, &
damage_nonlocal_postResults damage_nonlocal_postResults
contains contains
@ -45,283 +56,228 @@ contains
!> @details reads in material parameters, allocates arrays, and does sanity checks !> @details reads in material parameters, allocates arrays, and does sanity checks
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine damage_nonlocal_init 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 write(6,'(/,a)') ' <<<+- damage_'//DAMAGE_nonlocal_label//' init -+>>>'
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))
do h = 1, size(damage_type) maxNinstance = count(damage_type == DAMAGE_nonlocal_ID)
if (damage_type(h) /= DAMAGE_NONLOCAL_ID) cycle if (maxNinstance == 0) return
associate(prm => param(damage_typeInstance(h)), &
config => config_homogenization(h)) allocate(damage_nonlocal_sizePostResult (maxval(homogenization_Noutput),maxNinstance),source=0)
allocate(damage_nonlocal_output (maxval(homogenization_Noutput),maxNinstance))
instance = damage_typeInstance(h) damage_nonlocal_output = ''
outputs = config%getStrings('(output)',defaultVal=emptyStringArray) allocate(damage_nonlocal_Noutput (maxNinstance), source=0)
allocate(prm%outputID(0))
allocate(param(maxNinstance))
do i=1, size(outputs) do h = 1, size(damage_type)
outputID = undefined_ID if (damage_type(h) /= DAMAGE_NONLOCAL_ID) cycle
select case(outputs(i)) associate(prm => param(damage_typeInstance(h)), &
config => config_homogenization(h))
case ('damage')
damage_nonlocal_output(i,damage_typeInstance(h)) = outputs(i) instance = damage_typeInstance(h)
damage_nonlocal_Noutput(instance) = damage_nonlocal_Noutput(instance) + 1 outputs = config%getStrings('(output)',defaultVal=emptyStringArray)
damage_nonlocal_sizePostResult(i,damage_typeInstance(h)) = 1 allocate(prm%outputID(0))
prm%outputID = [prm%outputID , damage_ID]
end select do i=1, size(outputs)
outputID = undefined_ID
enddo 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) NofMyHomog = count(material_homogenizationAt == homog)
instance = damage_typeInstance(homog) instance = damage_typeInstance(homog)
! allocate state arrays ! allocate state arrays
sizeState = 1 sizeState = 1
damageState(homog)%sizeState = sizeState damageState(homog)%sizeState = sizeState
damageState(homog)%sizePostResults = sum(damage_nonlocal_sizePostResult(:,instance)) damageState(homog)%sizePostResults = sum(damage_nonlocal_sizePostResult(:,instance))
allocate(damageState(homog)%state0 (sizeState,NofMyHomog), source=damage_initialPhi(homog)) allocate(damageState(homog)%state0 (sizeState,NofMyHomog), source=damage_initialPhi(homog))
allocate(damageState(homog)%subState0(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(damageState(homog)%state (sizeState,NofMyHomog), source=damage_initialPhi(homog))
nullify(damageMapping(homog)%p) nullify(damageMapping(homog)%p)
damageMapping(homog)%p => mappingHomogenization(1,:,:) damageMapping(homog)%p => mappingHomogenization(1,:,:)
deallocate(damage(homog)%p) deallocate(damage(homog)%p)
damage(homog)%p => damageState(homog)%state(1,:) damage(homog)%p => damageState(homog)%state(1,:)
end associate end associate
enddo enddo
end subroutine damage_nonlocal_init end subroutine damage_nonlocal_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculates homogenized damage driving forces !> @brief calculates homogenized damage driving forces
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el) 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) :: & integer, intent(in) :: &
ip, & !< integration point number ip, & !< integration point number
el !< element number el !< element number
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
phi phi
integer :: & integer :: &
phase, & phase, &
grain, & grain, &
source, & source, &
constituent constituent
real(pReal) :: & real(pReal) :: &
phiDot, dPhiDot_dPhi, localphiDot, dLocalphiDot_dPhi phiDot, dPhiDot_dPhi, localphiDot, dLocalphiDot_dPhi
phiDot = 0.0_pReal phiDot = 0.0_pReal
dPhiDot_dPhi = 0.0_pReal dPhiDot_dPhi = 0.0_pReal
do grain = 1, homogenization_Ngrains(material_homogenizationAt(el)) do grain = 1, homogenization_Ngrains(material_homogenizationAt(el))
phase = phaseAt(grain,ip,el) phase = phaseAt(grain,ip,el)
constituent = phasememberAt(grain,ip,el) constituent = phasememberAt(grain,ip,el)
do source = 1, phase_Nsources(phase) do source = 1, phase_Nsources(phase)
select case(phase_source(source,phase)) select case(phase_source(source,phase))
case (SOURCE_damage_isoBrittle_ID) case (SOURCE_damage_isoBrittle_ID)
call source_damage_isobrittle_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) call source_damage_isobrittle_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
case (SOURCE_damage_isoDuctile_ID) case (SOURCE_damage_isoDuctile_ID)
call source_damage_isoductile_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) call source_damage_isoductile_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
case (SOURCE_damage_anisoBrittle_ID) case (SOURCE_damage_anisoBrittle_ID)
call source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) call source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
case (SOURCE_damage_anisoDuctile_ID) case (SOURCE_damage_anisoDuctile_ID)
call source_damage_anisoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) call source_damage_anisoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
case default case default
localphiDot = 0.0_pReal localphiDot = 0.0_pReal
dLocalphiDot_dPhi = 0.0_pReal dLocalphiDot_dPhi = 0.0_pReal
end select end select
phiDot = phiDot + localphiDot phiDot = phiDot + localphiDot
dPhiDot_dPhi = dPhiDot_dPhi + dLocalphiDot_dPhi dPhiDot_dPhi = dPhiDot_dPhi + dLocalphiDot_dPhi
enddo enddo
enddo enddo
phiDot = phiDot/real(homogenization_Ngrains(material_homogenizationAt(el)),pReal) phiDot = phiDot/real(homogenization_Ngrains(material_homogenizationAt(el)),pReal)
dPhiDot_dPhi = dPhiDot_dPhi/real(homogenization_Ngrains(material_homogenizationAt(el)),pReal) dPhiDot_dPhi = dPhiDot_dPhi/real(homogenization_Ngrains(material_homogenizationAt(el)),pReal)
end subroutine damage_nonlocal_getSourceAndItsTangent end subroutine damage_nonlocal_getSourceAndItsTangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief returns homogenized non local damage diffusion tensor in reference configuration !> @brief returns homogenized non local damage diffusion tensor in reference configuration
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function damage_nonlocal_getDiffusion33(ip,el) 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) :: & integer, intent(in) :: &
ip, & !< integration point number ip, & !< integration point number
el !< element number el !< element number
real(pReal), dimension(3,3) :: & real(pReal), dimension(3,3) :: &
damage_nonlocal_getDiffusion33 damage_nonlocal_getDiffusion33
integer :: & integer :: &
homog, & homog, &
grain grain
homog = material_homogenizationAt(el) homog = material_homogenizationAt(el)
damage_nonlocal_getDiffusion33 = 0.0_pReal damage_nonlocal_getDiffusion33 = 0.0_pReal
do grain = 1, homogenization_Ngrains(homog) do grain = 1, homogenization_Ngrains(homog)
damage_nonlocal_getDiffusion33 = damage_nonlocal_getDiffusion33 + & damage_nonlocal_getDiffusion33 = damage_nonlocal_getDiffusion33 + &
crystallite_push33ToRef(grain,ip,el,lattice_DamageDiffusion33(1:3,1:3,material_phase(grain,ip,el))) crystallite_push33ToRef(grain,ip,el,lattice_DamageDiffusion33(1:3,1:3,material_phase(grain,ip,el)))
enddo enddo
damage_nonlocal_getDiffusion33 = & damage_nonlocal_getDiffusion33 = &
charLength**2*damage_nonlocal_getDiffusion33/real(homogenization_Ngrains(homog),pReal) charLength**2*damage_nonlocal_getDiffusion33/real(homogenization_Ngrains(homog),pReal)
end function damage_nonlocal_getDiffusion33 end function damage_nonlocal_getDiffusion33
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Returns homogenized nonlocal damage mobility !> @brief Returns homogenized nonlocal damage mobility
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
real(pReal) function damage_nonlocal_getMobility(ip,el) 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) :: & integer, intent(in) :: &
ip, & !< integration point number ip, & !< integration point number
el !< element number el !< element number
integer :: & integer :: &
ipc ipc
damage_nonlocal_getMobility = 0.0_pReal damage_nonlocal_getMobility = 0.0_pReal
do ipc = 1, homogenization_Ngrains(mesh_element(3,el)) do ipc = 1, homogenization_Ngrains(mesh_element(3,el))
damage_nonlocal_getMobility = damage_nonlocal_getMobility + lattice_DamageMobility(material_phase(ipc,ip,el)) damage_nonlocal_getMobility = damage_nonlocal_getMobility + lattice_DamageMobility(material_phase(ipc,ip,el))
enddo enddo
damage_nonlocal_getMobility = damage_nonlocal_getMobility/& damage_nonlocal_getMobility = damage_nonlocal_getMobility/&
real(homogenization_Ngrains(mesh_element(3,el)),pReal) real(homogenization_Ngrains(mesh_element(3,el)),pReal)
end function damage_nonlocal_getMobility end function damage_nonlocal_getMobility
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief updated nonlocal damage field with solution from damage phase field PDE !> @brief updated nonlocal damage field with solution from damage phase field PDE
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine damage_nonlocal_putNonLocalDamage(phi,ip,el) subroutine damage_nonlocal_putNonLocalDamage(phi,ip,el)
use material, only: &
material_homogenizationAt, &
damageMapping, &
damage
integer, intent(in) :: & integer, intent(in) :: &
ip, & !< integration point number ip, & !< integration point number
el !< element number el !< element number
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
phi phi
integer :: & integer :: &
homog, & homog, &
offset offset
homog = material_homogenizationAt(el) homog = material_homogenizationAt(el)
offset = damageMapping(homog)%p(ip,el) offset = damageMapping(homog)%p(ip,el)
damage(homog)%p(offset) = phi damage(homog)%p(offset) = phi
end subroutine damage_nonlocal_putNonLocalDamage end subroutine damage_nonlocal_putNonLocalDamage
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief return array of damage results !> @brief return array of damage results
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function damage_nonlocal_postResults(ip,el) function damage_nonlocal_postResults(ip,el)
use material, only: &
material_homogenizationAt, &
damage_typeInstance, &
damageMapping, &
damage
integer, intent(in) :: & integer, intent(in) :: &
ip, & !< integration point ip, & !< integration point
el !< element el !< element
real(pReal), dimension(sum(damage_nonlocal_sizePostResult(:,damage_typeInstance(material_homogenizationAt(el))))) :: & real(pReal), dimension(sum(damage_nonlocal_sizePostResult(:,damage_typeInstance(material_homogenizationAt(el))))) :: &
damage_nonlocal_postResults damage_nonlocal_postResults
integer :: & integer :: &
instance, homog, offset, o, c instance, homog, offset, o, c
homog = material_homogenizationAt(el) homog = material_homogenizationAt(el)
offset = damageMapping(homog)%p(ip,el) offset = damageMapping(homog)%p(ip,el)
instance = damage_typeInstance(homog) instance = damage_typeInstance(homog)
associate(prm => param(instance)) associate(prm => param(instance))
c = 0 c = 0
outputsLoop: do o = 1,size(prm%outputID) outputsLoop: do o = 1,size(prm%outputID)
select case(prm%outputID(o)) select case(prm%outputID(o))
case (damage_ID) case (damage_ID)
damage_nonlocal_postResults(c+1) = damage(homog)%p(offset) damage_nonlocal_postResults(c+1) = damage(homog)%p(offset)
c = c + 1 c = c + 1
end select end select
enddo outputsLoop enddo outputsLoop
end associate end associate
end function damage_nonlocal_postResults end function damage_nonlocal_postResults
end module damage_nonlocal 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 !> @brief Reading in and interpretating the debugging settings for the various modules
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module debug module debug
use prec, only: & use prec
pInt, & use IO
pReal
implicit none implicit none
private private
integer(pInt), parameter, public :: & integer(pInt), parameter, public :: &
debug_LEVELSELECTIVE = 2_pInt**0_pInt, & debug_LEVELSELECTIVE = 2_pInt**0_pInt, &
debug_LEVELBASIC = 2_pInt**1_pInt, & debug_LEVELBASIC = 2_pInt**1_pInt, &
@ -78,19 +78,7 @@ contains
!> @brief reads in parameters from debug.config and allocates arrays !> @brief reads in parameters from debug.config and allocates arrays
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine debug_init 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 character(len=pStringLen), dimension(:), allocatable :: fileContent
integer :: i, what, j integer :: i, what, j
@ -253,8 +241,6 @@ end subroutine debug_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine debug_reset subroutine debug_reset
implicit none
debug_stressMaxLocation = 0_pInt debug_stressMaxLocation = 0_pInt
debug_stressMinLocation = 0_pInt debug_stressMinLocation = 0_pInt
debug_jacobianMaxLocation = 0_pInt debug_jacobianMaxLocation = 0_pInt
@ -272,8 +258,6 @@ end subroutine debug_reset
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine debug_info subroutine debug_info
implicit none
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
debugOutputCPFEM: if (iand(debug_level(debug_CPFEM),debug_LEVELBASIC) /= 0 & debugOutputCPFEM: if (iand(debug_level(debug_CPFEM),debug_LEVELBASIC) /= 0 &
.and. any(debug_stressMinLocation /= 0_pInt) & .and. any(debug_stressMinLocation /= 0_pInt) &

52
src/geometry_plastic_nonlocal.f90 Normal file
View File

@ -0,0 +1,52 @@
!--------------------------------------------------------------------------------------------------
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @author Christoph Koords, Max-Planck-Institut für Eisenforschung GmbH
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Geometric information about the IP cells needed for the nonlocal
! plasticity model
!--------------------------------------------------------------------------------------------------
module geometry_plastic_nonlocal
use prec
implicit none
private
logical, dimension(3), public, parameter :: &
geometry_plastic_nonlocal_periodicSurface = .true. !< flag indicating periodic outer surfaces (used for fluxes) NEEDED?
integer, dimension(:,:,:,:), allocatable, public, protected :: &
geometry_plastic_nonlocal_IPneighborhood !< 6 or less neighboring IPs as [element_num, IP_index, neighbor_index that points to me]
real(pReal), dimension(:,:), allocatable, public, protected :: &
geometry_plastic_nonlocal_IPvolume !< volume associated with IP (initially!)
real(pReal), dimension(:,:,:), allocatable, public, protected :: &
geometry_plastic_nonlocal_IParea !< area of interface to neighboring IP (initially!)
real(pReal),dimension(:,:,:,:), allocatable, public, protected :: &
geometry_plastic_nonlocal_IPareaNormal !< area normal of interface to neighboring IP (initially!)
public :: &
geometry_plastic_nonlocal_set_IPneighborhood, &
geometry_plastic_nonlocal_set_IPvolume
contains
subroutine geometry_plastic_nonlocal_set_IPneighborhood(IPneighborhood)
integer, dimension(:,:,:,:), intent(in) :: IPneighborhood
geometry_plastic_nonlocal_IPneighborhood = IPneighborhood
end subroutine geometry_plastic_nonlocal_set_IPneighborhood
subroutine geometry_plastic_nonlocal_set_IPvolume(IPvolume)
real(pReal), dimension(:,:), intent(in) :: IPvolume
geometry_plastic_nonlocal_IPvolume = IPvolume
end subroutine geometry_plastic_nonlocal_set_IPvolume
end module geometry_plastic_nonlocal

29
src/grid/spectral_utilities.f90
View File

@ -196,7 +196,6 @@ subroutine utilities_init
grid3Offset, & grid3Offset, &
geomSize geomSize
implicit none
PetscErrorCode :: ierr PetscErrorCode :: ierr
integer :: i, j, k, & integer :: i, j, k, &
FFTW_planner_flag FFTW_planner_flag
@ -425,7 +424,6 @@ subroutine utilities_updateGamma(C,saveReference)
math_det33, & math_det33, &
math_invert2 math_invert2
implicit none
real(pReal), intent(in), dimension(3,3,3,3) :: C !< input stiffness to store as reference stiffness real(pReal), intent(in), dimension(3,3,3,3) :: C !< input stiffness to store as reference stiffness
logical , intent(in) :: saveReference !< save reference stiffness to file for restart logical , intent(in) :: saveReference !< save reference stiffness to file for restart
complex(pReal), dimension(3,3) :: temp33_complex, xiDyad_cmplx complex(pReal), dimension(3,3) :: temp33_complex, xiDyad_cmplx
@ -473,7 +471,6 @@ end subroutine utilities_updateGamma
!> @details Does an unweighted filtered FFT transform from real to complex !> @details Does an unweighted filtered FFT transform from real to complex
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine utilities_FFTtensorForward subroutine utilities_FFTtensorForward
implicit none
call fftw_mpi_execute_dft_r2c(planTensorForth,tensorField_real,tensorField_fourier) call fftw_mpi_execute_dft_r2c(planTensorForth,tensorField_real,tensorField_fourier)
@ -485,7 +482,6 @@ end subroutine utilities_FFTtensorForward
!> @details Does an weighted inverse FFT transform from complex to real !> @details Does an weighted inverse FFT transform from complex to real
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine utilities_FFTtensorBackward subroutine utilities_FFTtensorBackward
implicit none
call fftw_mpi_execute_dft_c2r(planTensorBack,tensorField_fourier,tensorField_real) call fftw_mpi_execute_dft_c2r(planTensorBack,tensorField_fourier,tensorField_real)
tensorField_real = tensorField_real * wgt ! normalize the result by number of elements tensorField_real = tensorField_real * wgt ! normalize the result by number of elements
@ -497,7 +493,6 @@ end subroutine utilities_FFTtensorBackward
!> @details Does an unweighted filtered FFT transform from real to complex !> @details Does an unweighted filtered FFT transform from real to complex
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine utilities_FFTscalarForward subroutine utilities_FFTscalarForward
implicit none
call fftw_mpi_execute_dft_r2c(planScalarForth,scalarField_real,scalarField_fourier) call fftw_mpi_execute_dft_r2c(planScalarForth,scalarField_real,scalarField_fourier)
@ -509,7 +504,6 @@ end subroutine utilities_FFTscalarForward
!> @details Does an weighted inverse FFT transform from complex to real !> @details Does an weighted inverse FFT transform from complex to real
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine utilities_FFTscalarBackward subroutine utilities_FFTscalarBackward
implicit none
call fftw_mpi_execute_dft_c2r(planScalarBack,scalarField_fourier,scalarField_real) call fftw_mpi_execute_dft_c2r(planScalarBack,scalarField_fourier,scalarField_real)
scalarField_real = scalarField_real * wgt ! normalize the result by number of elements scalarField_real = scalarField_real * wgt ! normalize the result by number of elements
@ -522,7 +516,6 @@ end subroutine utilities_FFTscalarBackward
!> @details Does an unweighted filtered FFT transform from real to complex. !> @details Does an unweighted filtered FFT transform from real to complex.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine utilities_FFTvectorForward subroutine utilities_FFTvectorForward
implicit none
call fftw_mpi_execute_dft_r2c(planVectorForth,vectorField_real,vectorField_fourier) call fftw_mpi_execute_dft_r2c(planVectorForth,vectorField_real,vectorField_fourier)
@ -534,7 +527,6 @@ end subroutine utilities_FFTvectorForward
!> @details Does an weighted inverse FFT transform from complex to real !> @details Does an weighted inverse FFT transform from complex to real
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine utilities_FFTvectorBackward subroutine utilities_FFTvectorBackward
implicit none
call fftw_mpi_execute_dft_c2r(planVectorBack,vectorField_fourier,vectorField_real) call fftw_mpi_execute_dft_c2r(planVectorBack,vectorField_fourier,vectorField_real)
vectorField_real = vectorField_real * wgt ! normalize the result by number of elements vectorField_real = vectorField_real * wgt ! normalize the result by number of elements
@ -554,7 +546,6 @@ subroutine utilities_fourierGammaConvolution(fieldAim)
grid, & grid, &
grid3Offset grid3Offset
implicit none
real(pReal), intent(in), dimension(3,3) :: fieldAim !< desired average value of the field after convolution real(pReal), intent(in), dimension(3,3) :: fieldAim !< desired average value of the field after convolution
complex(pReal), dimension(3,3) :: temp33_complex, xiDyad_cmplx complex(pReal), dimension(3,3) :: temp33_complex, xiDyad_cmplx
real(pReal), dimension(6,6) :: A, A_inv real(pReal), dimension(6,6) :: A, A_inv
@ -615,7 +606,6 @@ subroutine utilities_fourierGreenConvolution(D_ref, mobility_ref, deltaT)
grid, & grid, &
grid3 grid3
implicit none
real(pReal), dimension(3,3), intent(in) :: D_ref real(pReal), dimension(3,3), intent(in) :: D_ref
real(pReal), intent(in) :: mobility_ref, deltaT real(pReal), intent(in) :: mobility_ref, deltaT
complex(pReal) :: GreenOp_hat complex(pReal) :: GreenOp_hat
@ -644,7 +634,6 @@ real(pReal) function utilities_divergenceRMS()
grid, & grid, &
grid3 grid3
implicit none
integer :: i, j, k, ierr integer :: i, j, k, ierr
complex(pReal), dimension(3) :: rescaledGeom complex(pReal), dimension(3) :: rescaledGeom
@ -694,7 +683,6 @@ real(pReal) function utilities_curlRMS()
grid, & grid, &
grid3 grid3
implicit none
integer :: i, j, k, l, ierr integer :: i, j, k, l, ierr
complex(pReal), dimension(3,3) :: curl_fourier complex(pReal), dimension(3,3) :: curl_fourier
complex(pReal), dimension(3) :: rescaledGeom complex(pReal), dimension(3) :: rescaledGeom
@ -766,7 +754,6 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
math_rotate_forward33, & math_rotate_forward33, &
math_invert2 math_invert2
implicit none
real(pReal), dimension(3,3,3,3) :: utilities_maskedCompliance !< masked compliance real(pReal), dimension(3,3,3,3) :: utilities_maskedCompliance !< masked compliance
real(pReal), intent(in) , dimension(3,3,3,3) :: C !< current average stiffness real(pReal), intent(in) , dimension(3,3,3,3) :: C !< current average stiffness
real(pReal), intent(in) , dimension(3,3) :: rot_BC !< rotation of load frame real(pReal), intent(in) , dimension(3,3) :: rot_BC !< rotation of load frame
@ -861,7 +848,6 @@ subroutine utilities_fourierScalarGradient()
grid3, & grid3, &
grid grid
implicit none
integer :: i, j, k integer :: i, j, k
vectorField_fourier = cmplx(0.0_pReal,0.0_pReal,pReal) vectorField_fourier = cmplx(0.0_pReal,0.0_pReal,pReal)
@ -879,7 +865,6 @@ subroutine utilities_fourierVectorDivergence()
grid3, & grid3, &
grid grid
implicit none
integer :: i, j, k integer :: i, j, k
scalarField_fourier = cmplx(0.0_pReal,0.0_pReal,pReal) scalarField_fourier = cmplx(0.0_pReal,0.0_pReal,pReal)
@ -898,7 +883,6 @@ subroutine utilities_fourierVectorGradient()
grid3, & grid3, &
grid grid
implicit none
integer :: i, j, k, m, n integer :: i, j, k, m, n
tensorField_fourier = cmplx(0.0_pReal,0.0_pReal,pReal) tensorField_fourier = cmplx(0.0_pReal,0.0_pReal,pReal)
@ -919,7 +903,6 @@ subroutine utilities_fourierTensorDivergence()
grid3, & grid3, &
grid grid
implicit none
integer :: i, j, k, m, n integer :: i, j, k, m, n
vectorField_fourier = cmplx(0.0_pReal,0.0_pReal,pReal) vectorField_fourier = cmplx(0.0_pReal,0.0_pReal,pReal)
@ -942,9 +925,6 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
IO_error IO_error
use numerics, only: & use numerics, only: &
worldrank worldrank
use debug, only: &
debug_reset, &
debug_info
use math, only: & use math, only: &
math_rotate_forward33, & math_rotate_forward33, &
math_det33 math_det33
@ -957,7 +937,6 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
materialpoint_dPdF, & materialpoint_dPdF, &
materialpoint_stressAndItsTangent materialpoint_stressAndItsTangent
implicit none
real(pReal),intent(out), dimension(3,3,3,3) :: C_volAvg, C_minmaxAvg !< average stiffness real(pReal),intent(out), dimension(3,3,3,3) :: C_volAvg, C_minmaxAvg !< average stiffness
real(pReal),intent(out), dimension(3,3) :: P_av !< average PK stress real(pReal),intent(out), dimension(3,3) :: P_av !< average PK stress
real(pReal),intent(out), dimension(3,3,grid(1),grid(2),grid3) :: P !< PK stress real(pReal),intent(out), dimension(3,3,grid(1),grid(2),grid3) :: P !< PK stress
@ -977,7 +956,6 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
materialpoint_F = reshape(F,[3,3,1,product(grid(1:2))*grid3]) ! set materialpoint target F to estimated field materialpoint_F = reshape(F,[3,3,1,product(grid(1:2))*grid3]) ! set materialpoint target F to estimated field
call debug_reset() ! this has no effect on rank >0
call materialpoint_stressAndItsTangent(.true.,timeinc) ! calculate P field call materialpoint_stressAndItsTangent(.true.,timeinc) ! calculate P field
P = reshape(materialpoint_P, [3,3,grid(1),grid(2),grid3]) P = reshape(materialpoint_P, [3,3,grid(1),grid(2),grid3])
@ -1023,8 +1001,7 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
C_volAvg = sum(sum(materialpoint_dPdF,dim=6),dim=5) C_volAvg = sum(sum(materialpoint_dPdF,dim=6),dim=5)
call MPI_Allreduce(MPI_IN_PLACE,C_volAvg,81,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,C_volAvg,81,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
C_volAvg = C_volAvg * wgt C_volAvg = C_volAvg * wgt
call debug_info() ! this has no effect on rank >0
end subroutine utilities_constitutiveResponse end subroutine utilities_constitutiveResponse
@ -1037,7 +1014,6 @@ pure function utilities_calculateRate(heterogeneous,field0,field,dt,avRate)
grid3, & grid3, &
grid grid
implicit none
real(pReal), intent(in), dimension(3,3) :: & real(pReal), intent(in), dimension(3,3) :: &
avRate !< homogeneous addon avRate !< homogeneous addon
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
@ -1068,7 +1044,6 @@ function utilities_forwardField(timeinc,field_lastInc,rate,aim)
grid3, & grid3, &
grid grid
implicit none
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
timeinc !< timeinc of current step timeinc !< timeinc of current step
real(pReal), intent(in), dimension(3,3,grid(1),grid(2),grid3) :: & real(pReal), intent(in), dimension(3,3,grid(1),grid(2),grid3) :: &
@ -1105,7 +1080,6 @@ pure function utilities_getFreqDerivative(k_s)
geomSize, & geomSize, &
grid grid
implicit none
integer, intent(in), dimension(3) :: k_s !< indices of frequency integer, intent(in), dimension(3) :: k_s !< indices of frequency
complex(pReal), dimension(3) :: utilities_getFreqDerivative complex(pReal), dimension(3) :: utilities_getFreqDerivative
@ -1163,7 +1137,6 @@ subroutine utilities_updateIPcoords(F)
grid3Offset, & grid3Offset, &
geomSize, & geomSize, &
mesh_ipCoordinates mesh_ipCoordinates
implicit none
real(pReal), dimension(3,3,grid(1),grid(2),grid3), intent(in) :: F real(pReal), dimension(3,3,grid(1),grid(2),grid3), intent(in) :: F
integer :: i, j, k, m, ierr integer :: i, j, k, m, ierr

24
src/homogenization.f90
View File

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

335
src/kinematics_cleavage_opening.f90
View File

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

44
src/kinematics_slipplane_opening.f90
View File

@ -6,12 +6,19 @@
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module kinematics_slipplane_opening module kinematics_slipplane_opening
use prec use prec
use config
use IO
use debug
use math
use lattice
use material
implicit none implicit none
private 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 :: & integer :: &
totalNslip totalNslip
integer, dimension(:), allocatable :: & integer, dimension(:), allocatable :: &
@ -19,7 +26,7 @@ module kinematics_slipplane_opening
real(pReal) :: & real(pReal) :: &
sdot0, & sdot0, &
n n
real(pReal), dimension(:), allocatable :: & real(pReal), dimension(:), allocatable :: &
critLoad critLoad
real(pReal), dimension(:,:), allocatable :: & real(pReal), dimension(:,:), allocatable :: &
slip_direction, & slip_direction, &
@ -27,7 +34,8 @@ module kinematics_slipplane_opening
slip_transverse slip_transverse
end type tParameters 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 :: & public :: &
kinematics_slipplane_opening_init, & kinematics_slipplane_opening_init, &
kinematics_slipplane_opening_LiAndItsTangent kinematics_slipplane_opening_LiAndItsTangent
@ -39,25 +47,9 @@ contains
!> @brief module initialization !> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks !> @details reads in material parameters, allocates arrays, and does sanity checks
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine kinematics_slipplane_opening_init() 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
integer :: maxNinstance,p,instance
integer :: maxNinstance,p,instance,kinematics
write(6,'(/,a)') ' <<<+- kinematics_'//KINEMATICS_slipplane_opening_LABEL//' init -+>>>' 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 !> @brief contains the constitutive equation for calculating the velocity gradient
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, ip, el) 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) :: & integer, intent(in) :: &
ipc, & !< grain number ipc, & !< grain number

50
src/kinematics_thermal_expansion.f90
View File

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

208
src/lattice.f90
View File

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

34
src/list.f90
View File

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

29
src/math.f90
View File

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

14
src/mesh/DAMASK_FEM.f90
View File

@ -28,8 +28,7 @@ program DAMASK_FEM
IO_intOut, & IO_intOut, &
IO_warning IO_warning
use math ! need to include the whole module for FFTW use math ! need to include the whole module for FFTW
use CPFEM2, only: & use CPFEM2
CPFEM_initAll
use FEsolving, only: & use FEsolving, only: &
restartWrite, & restartWrite, &
restartInc restartInc
@ -114,7 +113,7 @@ program DAMASK_FEM
write(6,'(/,a)') ' <<<+- DAMASK_FEM init -+>>>' write(6,'(/,a)') ' <<<+- DAMASK_FEM init -+>>>'
! reading basic information from load case file and allocate data structure containing load cases ! reading basic information from load case file and allocate data structure containing load cases
call DMGetDimension(geomMesh,dimPlex,ierr)! CHKERRQ(ierr) !< dimension of mesh (2D or 3D) call DMGetDimension(geomMesh,dimPlex,ierr); CHKERRA(ierr) !< dimension of mesh (2D or 3D)
nActiveFields = 1 nActiveFields = 1
allocate(solres(nActiveFields)) allocate(solres(nActiveFields))
@ -394,8 +393,7 @@ program DAMASK_FEM
cutBack = .False. cutBack = .False.
if(.not. all(solres(:)%converged .and. solres(:)%stagConverged)) then ! no solution found if(.not. all(solres(:)%converged .and. solres(:)%stagConverged)) then ! no solution found
if (cutBackLevel < maxCutBack) then ! do cut back if (cutBackLevel < maxCutBack) then ! do cut back
if (worldrank == 0) & write(6,'(/,a)') ' cut back detected'
write(6,'(/,a)') ' cut back detected'
cutBack = .True. cutBack = .True.
stepFraction = (stepFraction - 1_pInt) * subStepFactor ! adjust to new denominator stepFraction = (stepFraction - 1_pInt) * subStepFactor ! adjust to new denominator
cutBackLevel = cutBackLevel + 1_pInt cutBackLevel = cutBackLevel + 1_pInt
@ -403,7 +401,7 @@ program DAMASK_FEM
timeinc = timeinc/2.0_pReal timeinc = timeinc/2.0_pReal
else ! default behavior, exit if spectral solver does not converge else ! default behavior, exit if spectral solver does not converge
call IO_warning(850_pInt) call IO_warning(850_pInt)
call quit(-1_pInt*(lastRestartWritten+1_pInt)) ! quit and provide information about last restart inc written (e.g. for regridding) ! continue from non-converged solution and start guessing after accepted (sub)inc call quit(-1_pInt*(lastRestartWritten+1_pInt)) ! quit and provide information about last restart inc written
endif endif
else else
guess = .true. ! start guessing after first converged (sub)inc guess = .true. ! start guessing after first converged (sub)inc
@ -428,7 +426,8 @@ program DAMASK_FEM
endif; flush(6) endif; flush(6)
if (mod(inc,loadCases(currentLoadCase)%outputFrequency) == 0_pInt) then ! at output frequency if (mod(inc,loadCases(currentLoadCase)%outputFrequency) == 0_pInt) then ! at output frequency
write(6,'(1/,a)') ' ToDo: ... writing results to file ......................................' write(6,'(1/,a)') ' ... writing results to file ......................................'
call CPFEM_results(totalIncsCounter,time)
endif endif
if ( loadCases(currentLoadCase)%restartFrequency > 0_pInt & ! writing of restart info requested ... if ( loadCases(currentLoadCase)%restartFrequency > 0_pInt & ! writing of restart info requested ...
.and. mod(inc,loadCases(currentLoadCase)%restartFrequency) == 0_pInt) then ! ... and at frequency of writing restart information .and. mod(inc,loadCases(currentLoadCase)%restartFrequency) == 0_pInt) then ! ... and at frequency of writing restart information
@ -452,7 +451,6 @@ program DAMASK_FEM
real(convergedCounter, pReal)/& real(convergedCounter, pReal)/&
real(notConvergedCounter + convergedCounter,pReal)*100.0_pReal, ' %) increments converged!' real(notConvergedCounter + convergedCounter,pReal)*100.0_pReal, ' %) increments converged!'
flush(6) flush(6)
call MPI_file_close(fileUnit,ierr)
close(statUnit) close(statUnit)
if (notConvergedCounter > 0_pInt) call quit(2_pInt) ! error if some are not converged if (notConvergedCounter > 0_pInt) call quit(2_pInt) ! error if some are not converged

66
src/mesh/FEM_mech.f90
View File

@ -84,11 +84,9 @@ subroutine FEM_mech_init(fieldBC)
PetscDS :: mechDS PetscDS :: mechDS
PetscDualSpace :: mechDualSpace PetscDualSpace :: mechDualSpace
DMLabel :: BCLabel DMLabel :: BCLabel
PetscInt, dimension(:), allocatable, target :: numComp, numDoF, bcField
PetscInt, dimension(:), pointer :: pNumComp, pNumDof, pBcField, pBcPoint PetscInt, dimension(:), pointer :: pNumComp, pNumDof, pBcField, pBcPoint
PetscInt :: numBC, bcSize, nc PetscInt :: numBC, bcSize, nc
IS :: bcPoint IS :: bcPoint
IS, allocatable, target :: bcComps(:), bcPoints(:)
IS, pointer :: pBcComps(:), pBcPoints(:) IS, pointer :: pBcComps(:), pBcPoints(:)
PetscSection :: section PetscSection :: section
PetscInt :: field, faceSet, topologDim, nNodalPoints PetscInt :: field, faceSet, topologDim, nNodalPoints
@ -98,7 +96,7 @@ subroutine FEM_mech_init(fieldBC)
PetscScalar, pointer :: px_scal(:) PetscScalar, pointer :: px_scal(:)
PetscScalar, allocatable, target :: x_scal(:) PetscScalar, allocatable, target :: x_scal(:)
PetscReal :: detJ PetscReal :: detJ
PetscReal, allocatable, target :: v0(:), cellJ(:), invcellJ(:), cellJMat(:,:) PetscReal, allocatable, target :: cellJMat(:,:)
PetscReal, pointer :: pV0(:), pCellJ(:), pInvcellJ(:) PetscReal, pointer :: pV0(:), pCellJ(:), pInvcellJ(:)
PetscInt :: cellStart, cellEnd, cell, basis PetscInt :: cellStart, cellEnd, cell, basis
character(len=7) :: prefix = 'mechFE_' character(len=7) :: prefix = 'mechFE_'
@ -139,26 +137,26 @@ subroutine FEM_mech_init(fieldBC)
call DMGetLabel(mech_mesh,'Face Sets',BCLabel,ierr); CHKERRQ(ierr) call DMGetLabel(mech_mesh,'Face Sets',BCLabel,ierr); CHKERRQ(ierr)
call DMPlexLabelComplete(mech_mesh,BCLabel,ierr); CHKERRQ(ierr) call DMPlexLabelComplete(mech_mesh,BCLabel,ierr); CHKERRQ(ierr)
call DMGetSection(mech_mesh,section,ierr); CHKERRQ(ierr) call DMGetSection(mech_mesh,section,ierr); CHKERRQ(ierr)
allocate(numComp(1), source=dimPlex); pNumComp => numComp allocate(pnumComp(1), source=dimPlex)
allocate(numDof(dimPlex+1), source = 0); pNumDof => numDof allocate(pnumDof(dimPlex+1), source = 0)
do topologDim = 0, dimPlex do topologDim = 0, dimPlex
call DMPlexGetDepthStratum(mech_mesh,topologDim,cellStart,cellEnd,ierr) call DMPlexGetDepthStratum(mech_mesh,topologDim,cellStart,cellEnd,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call PetscSectionGetDof(section,cellStart,numDof(topologDim+1),ierr) call PetscSectionGetDof(section,cellStart,pnumDof(topologDim+1),ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
enddo enddo
numBC = 0 numBC = 0
do field = 1, dimPlex; do faceSet = 1, mesh_Nboundaries do field = 1, dimPlex; do faceSet = 1, mesh_Nboundaries
if (fieldBC%componentBC(field)%Mask(faceSet)) numBC = numBC + 1 if (fieldBC%componentBC(field)%Mask(faceSet)) numBC = numBC + 1
enddo; enddo enddo; enddo
allocate(bcField(numBC), source=0); pBcField => bcField allocate(pbcField(numBC), source=0)
allocate(bcComps(numBC)); pBcComps => bcComps allocate(pbcComps(numBC))
allocate(bcPoints(numBC)); pBcPoints => bcPoints allocate(pbcPoints(numBC))
numBC = 0 numBC = 0
do field = 1, dimPlex; do faceSet = 1, mesh_Nboundaries do field = 1, dimPlex; do faceSet = 1, mesh_Nboundaries
if (fieldBC%componentBC(field)%Mask(faceSet)) then if (fieldBC%componentBC(field)%Mask(faceSet)) then
numBC = numBC + 1 numBC = numBC + 1
call ISCreateGeneral(PETSC_COMM_WORLD,1,[field-1],PETSC_COPY_VALUES,bcComps(numBC),ierr) call ISCreateGeneral(PETSC_COMM_WORLD,1,[field-1],PETSC_COPY_VALUES,pbcComps(numBC),ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call DMGetStratumSize(mech_mesh,'Face Sets',mesh_boundaries(faceSet),bcSize,ierr) call DMGetStratumSize(mech_mesh,'Face Sets',mesh_boundaries(faceSet),bcSize,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
@ -166,12 +164,12 @@ subroutine FEM_mech_init(fieldBC)
call DMGetStratumIS(mech_mesh,'Face Sets',mesh_boundaries(faceSet),bcPoint,ierr) call DMGetStratumIS(mech_mesh,'Face Sets',mesh_boundaries(faceSet),bcPoint,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call ISGetIndicesF90(bcPoint,pBcPoint,ierr); CHKERRQ(ierr) call ISGetIndicesF90(bcPoint,pBcPoint,ierr); CHKERRQ(ierr)
call ISCreateGeneral(PETSC_COMM_WORLD,bcSize,pBcPoint,PETSC_COPY_VALUES,bcPoints(numBC),ierr) call ISCreateGeneral(PETSC_COMM_WORLD,bcSize,pBcPoint,PETSC_COPY_VALUES,pbcPoints(numBC),ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call ISRestoreIndicesF90(bcPoint,pBcPoint,ierr); CHKERRQ(ierr) call ISRestoreIndicesF90(bcPoint,pBcPoint,ierr); CHKERRQ(ierr)
call ISDestroy(bcPoint,ierr); CHKERRQ(ierr) call ISDestroy(bcPoint,ierr); CHKERRQ(ierr)
else else
call ISCreateGeneral(PETSC_COMM_WORLD,0,[0],PETSC_COPY_VALUES,bcPoints(numBC),ierr) call ISCreateGeneral(PETSC_COMM_WORLD,0,[0],PETSC_COPY_VALUES,pbcPoints(numBC),ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
endif endif
endif endif
@ -182,7 +180,7 @@ subroutine FEM_mech_init(fieldBC)
CHKERRQ(ierr) CHKERRQ(ierr)
call DMSetSection(mech_mesh,section,ierr); CHKERRQ(ierr) call DMSetSection(mech_mesh,section,ierr); CHKERRQ(ierr)
do faceSet = 1, numBC do faceSet = 1, numBC
call ISDestroy(bcPoints(faceSet),ierr); CHKERRQ(ierr) call ISDestroy(pbcPoints(faceSet),ierr); CHKERRQ(ierr)
enddo enddo
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -213,13 +211,10 @@ subroutine FEM_mech_init(fieldBC)
allocate(nodalWeights(1)) allocate(nodalWeights(1))
nodalPointsP => nodalPoints nodalPointsP => nodalPoints
nodalWeightsP => nodalWeights nodalWeightsP => nodalWeights
allocate(v0(dimPlex)) allocate(pv0(dimPlex))
allocate(cellJ(dimPlex*dimPlex)) allocate(pcellJ(dimPlex*dimPlex))
allocate(invcellJ(dimPlex*dimPlex)) allocate(pinvcellJ(dimPlex*dimPlex))
allocate(cellJMat(dimPlex,dimPlex)) allocate(cellJMat(dimPlex,dimPlex))
pV0 => v0
pCellJ => cellJ
pInvcellJ => invcellJ
call DMGetSection(mech_mesh,section,ierr); CHKERRQ(ierr) call DMGetSection(mech_mesh,section,ierr); CHKERRQ(ierr)
call DMGetDS(mech_mesh,mechDS,ierr); CHKERRQ(ierr) call DMGetDS(mech_mesh,mechDS,ierr); CHKERRQ(ierr)
call PetscDSGetDiscretization(mechDS,0,mechFE,ierr) call PetscDSGetDiscretization(mechDS,0,mechFE,ierr)
@ -325,22 +320,19 @@ subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr)
PetscScalar, dimension(:), pointer :: x_scal, pf_scal PetscScalar, dimension(:), pointer :: x_scal, pf_scal
PetscScalar, target :: f_scal(cellDof) PetscScalar, target :: f_scal(cellDof)
PetscReal :: detJ, IcellJMat(dimPlex,dimPlex) PetscReal :: detJ, IcellJMat(dimPlex,dimPlex)
PetscReal, target :: v0(dimPlex), cellJ(dimPlex*dimPlex), & PetscReal, pointer,dimension(:) :: pV0, pCellJ, pInvcellJ, basisField, basisFieldDer
invcellJ(dimPlex*dimPlex)
PetscReal, pointer :: pV0(:), pCellJ(:), pInvcellJ(:)
PetscReal, pointer :: basisField(:), basisFieldDer(:)
PetscInt :: cellStart, cellEnd, cell, field, face, & PetscInt :: cellStart, cellEnd, cell, field, face, &
qPt, basis, comp, cidx qPt, basis, comp, cidx
PetscReal :: detFAvg PetscReal :: detFAvg
PetscReal :: BMat(dimPlex*dimPlex,cellDof) PetscReal :: BMat(dimPlex*dimPlex,cellDof)
PetscObject :: dummy PetscObject,intent(in) :: dummy
PetscInt :: bcSize PetscInt :: bcSize
IS :: bcPoints IS :: bcPoints
PetscErrorCode :: ierr PetscErrorCode :: ierr
pV0 => v0 allocate(pV0(dimPlex))
pCellJ => cellJ allocate(pcellJ(dimPlex**2))
pInvcellJ => invcellJ allocate(pinvcellJ(dimPlex**2))
call DMGetSection(dm_local,section,ierr); CHKERRQ(ierr) call DMGetSection(dm_local,section,ierr); CHKERRQ(ierr)
call DMGetDS(dm_local,prob,ierr); CHKERRQ(ierr) call DMGetDS(dm_local,prob,ierr); CHKERRQ(ierr)
call PetscDSGetTabulation(prob,0,basisField,basisFieldDer,ierr) call PetscDSGetTabulation(prob,0,basisField,basisFieldDer,ierr)
@ -460,13 +452,11 @@ subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr)
Vec :: x_local, xx_local Vec :: x_local, xx_local
Mat :: Jac_pre, Jac Mat :: Jac_pre, Jac
PetscSection :: section, gSection PetscSection :: section, gSection
PetscReal :: detJ, IcellJMat(dimPlex,dimPlex) PetscReal :: detJ
PetscReal, target :: v0(dimPlex), cellJ(dimPlex*dimPlex), &
invcellJ(dimPlex*dimPlex)
PetscReal, dimension(:), pointer :: basisField, basisFieldDer, & PetscReal, dimension(:), pointer :: basisField, basisFieldDer, &
pV0, pCellJ, pInvcellJ pV0, pCellJ, pInvcellJ
PetscInt :: cellStart, cellEnd, cell, field, face, & PetscInt :: cellStart, cellEnd, cell, field, face, &
qPt, basis, comp, cidx qPt, basis, comp, cidx,bcSize
PetscScalar,dimension(cellDOF,cellDOF), target :: K_e, & PetscScalar,dimension(cellDOF,cellDOF), target :: K_e, &
K_eA , & K_eA , &
K_eB K_eB
@ -477,14 +467,14 @@ subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr)
MatB (1 ,cellDof) MatB (1 ,cellDof)
PetscScalar, dimension(:), pointer :: pK_e, x_scal PetscScalar, dimension(:), pointer :: pK_e, x_scal
PetscReal, dimension(3,3) :: F, FAvg, FInv PetscReal, dimension(3,3) :: F, FAvg, FInv
PetscObject :: dummy PetscObject, intent(in) :: dummy
PetscInt :: bcSize
IS :: bcPoints IS :: bcPoints
PetscErrorCode :: ierr PetscErrorCode :: ierr
pV0 => v0 allocate(pV0(dimPlex))
pCellJ => cellJ allocate(pcellJ(dimPlex**2))
pInvcellJ => invcellJ allocate(pinvcellJ(dimPlex**2))
call MatSetOption(Jac,MAT_KEEP_NONZERO_PATTERN,PETSC_TRUE,ierr); CHKERRQ(ierr) call MatSetOption(Jac,MAT_KEEP_NONZERO_PATTERN,PETSC_TRUE,ierr); CHKERRQ(ierr)
call MatSetOption(Jac,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE,ierr); CHKERRQ(ierr) call MatSetOption(Jac,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE,ierr); CHKERRQ(ierr)
call MatZeroEntries(Jac,ierr); CHKERRQ(ierr) call MatZeroEntries(Jac,ierr); CHKERRQ(ierr)
@ -513,7 +503,6 @@ subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call DMPlexComputeCellGeometryAffineFEM(dm_local,cell,pV0,pCellJ,pInvcellJ,detJ,ierr) call DMPlexComputeCellGeometryAffineFEM(dm_local,cell,pV0,pCellJ,pInvcellJ,detJ,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
IcellJMat = reshape(pInvcellJ, shape = [dimPlex,dimPlex])
K_eA = 0.0 K_eA = 0.0
K_eB = 0.0 K_eB = 0.0
MatB = 0.0 MatB = 0.0
@ -525,7 +514,8 @@ subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr)
do comp = 0, dimPlex-1 do comp = 0, dimPlex-1
cidx = basis*dimPlex+comp cidx = basis*dimPlex+comp
BMat(comp*dimPlex+1:(comp+1)*dimPlex,basis*dimPlex+comp+1) = & BMat(comp*dimPlex+1:(comp+1)*dimPlex,basis*dimPlex+comp+1) = &
matmul(IcellJMat,basisFieldDer((((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp )*dimPlex+1: & matmul( reshape(pInvcellJ, shape = [dimPlex,dimPlex]),&
basisFieldDer((((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp )*dimPlex+1: &
(((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp+1)*dimPlex)) (((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp+1)*dimPlex))
enddo enddo
enddo enddo

40
src/mesh/FEM_utilities.f90
View File

@ -24,9 +24,7 @@ use PETScis
! grid related information information ! grid related information information
real(pReal), public :: wgt !< weighting factor 1/Nelems real(pReal), public :: wgt !< weighting factor 1/Nelems
!--------------------------------------------------------------------------------------------------
! output data
Vec, public :: coordinatesVec
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! field labels information ! field labels information
character(len=*), parameter, public :: & character(len=*), parameter, public :: &
@ -53,7 +51,6 @@ use PETScis
type, public :: tSolutionState !< return type of solution from FEM solver variants type, public :: tSolutionState !< return type of solution from FEM solver variants
logical :: converged = .true. logical :: converged = .true.
logical :: stagConverged = .true. logical :: stagConverged = .true.
logical :: regrid = .false.
integer(pInt) :: iterationsNeeded = 0_pInt integer(pInt) :: iterationsNeeded = 0_pInt
end type tSolutionState end type tSolutionState
@ -79,18 +76,6 @@ use PETScis
integer(pInt), allocatable :: faceID(:) integer(pInt), allocatable :: faceID(:)
type(tFieldBC), allocatable :: fieldBC(:) type(tFieldBC), allocatable :: fieldBC(:)
end type tLoadCase end type tLoadCase
type, public :: tFEMInterpolation
integer(pInt) :: n
real(pReal), dimension(:,:) , allocatable :: shapeFunc, shapeDerivReal, geomShapeDerivIso
real(pReal), dimension(:,:,:), allocatable :: shapeDerivIso
end type tFEMInterpolation
type, public :: tQuadrature
integer(pInt) :: n
real(pReal), dimension(:) , allocatable :: Weights
real(pReal), dimension(:,:), allocatable :: Points
end type tQuadrature
public :: & public :: &
utilities_init, & utilities_init, &
@ -119,11 +104,8 @@ subroutine utilities_init
use math ! must use the whole module for use of FFTW use math ! must use the whole module for use of FFTW
use mesh, only: & use mesh, only: &
mesh_NcpElemsGlobal, & mesh_NcpElemsGlobal, &
mesh_maxNips, & mesh_maxNips
geomMesh
implicit none
character(len=1024) :: petsc_optionsPhysics character(len=1024) :: petsc_optionsPhysics
PetscErrorCode :: ierr PetscErrorCode :: ierr
@ -157,35 +139,21 @@ end subroutine utilities_init
!> @brief calculates constitutive response !> @brief calculates constitutive response
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData) subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData)
use math, only: &
math_det33
use FEsolving, only: & use FEsolving, only: &
restartWrite restartWrite
use homogenization, only: & use homogenization, only: &
materialpoint_P, & materialpoint_P, &
materialpoint_stressAndItsTangent materialpoint_stressAndItsTangent
implicit none
real(pReal), intent(in) :: timeinc !< loading time real(pReal), intent(in) :: timeinc !< loading time
logical, intent(in) :: forwardData !< age results logical, intent(in) :: forwardData !< age results
real(pReal),intent(out), dimension(3,3) :: P_av !< average PK stress real(pReal),intent(out), dimension(3,3) :: P_av !< average PK stress
logical :: &
age
PetscErrorCode :: ierr PetscErrorCode :: ierr
write(6,'(/,a)') ' ... evaluating constitutive response ......................................' write(6,'(/,a)') ' ... evaluating constitutive response ......................................'
age = .False.
if (forwardData) then ! aging results
age = .True.
endif
if (cutBack) then ! restore saved variables
age = .False.
endif
call materialpoint_stressAndItsTangent(.true.,timeinc) ! calculate P field call materialpoint_stressAndItsTangent(.true.,timeinc) ! calculate P field
restartWrite = .false. ! reset restartWrite status restartWrite = .false. ! reset restartWrite status
@ -202,8 +170,6 @@ end subroutine utilities_constitutiveResponse
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine utilities_projectBCValues(localVec,section,field,comp,bcPointsIS,BCValue,BCDotValue,timeinc) subroutine utilities_projectBCValues(localVec,section,field,comp,bcPointsIS,BCValue,BCDotValue,timeinc)
implicit none
Vec :: localVec Vec :: localVec
PetscInt :: field, comp, nBcPoints, point, dof, numDof, numComp, offset PetscInt :: field, comp, nBcPoints, point, dof, numDof, numComp, offset
PetscSection :: section PetscSection :: section

562
src/mesh_grid.f90
View File

@ -7,18 +7,14 @@
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module mesh module mesh
use, intrinsic :: iso_c_binding use, intrinsic :: iso_c_binding
use prec, only: pReal, pInt use prec
use geometry_plastic_nonlocal
use mesh_base use mesh_base
implicit none implicit none
private private
integer(pInt), public, protected :: & integer(pInt), public, protected :: &
mesh_Nnodes, & !< total number of nodes in mesh mesh_Nnodes
mesh_Ncellnodes, & !< total number of cell nodes in mesh (including duplicates)
mesh_Ncells, & !< total number of cells in mesh
mesh_maxNipNeighbors, & !< max number of IP neighbors in any CP element
mesh_maxNsharedElems !< max number of CP elements sharing a node
integer(pInt), dimension(:), allocatable, private :: & integer(pInt), dimension(:), allocatable, private :: &
microGlobal microGlobal
@ -34,9 +30,9 @@ module mesh
real(pReal), public, protected :: & real(pReal), public, protected :: &
mesh_unitlength !< physical length of one unit in mesh mesh_unitlength !< physical length of one unit in mesh
real(pReal), dimension(:,:), allocatable, public :: & real(pReal), dimension(:,:), allocatable, private :: &
mesh_node, & !< node x,y,z coordinates (after deformation! ONLY FOR MARC!!!) mesh_node !< node x,y,z coordinates (after deformation! ONLY FOR MARC!!!)
mesh_cellnode !< cell node x,y,z coordinates (after deformation! ONLY FOR MARC!!!)
real(pReal), dimension(:,:), allocatable, public, protected :: & real(pReal), dimension(:,:), allocatable, public, protected :: &
mesh_ipVolume, & !< volume associated with IP (initially!) mesh_ipVolume, & !< volume associated with IP (initially!)
@ -53,56 +49,8 @@ module mesh
logical, dimension(3), public, parameter :: mesh_periodicSurface = .true. !< flag indicating periodic outer surfaces (used for fluxes) logical, dimension(3), public, parameter :: mesh_periodicSurface = .true. !< flag indicating periodic outer surfaces (used for fluxes)
integer(pInt), dimension(:,:), allocatable, private :: &
mesh_cellnodeParent !< cellnode's parent element ID, cellnode's intra-element ID
integer(pInt),dimension(:,:,:), allocatable, private :: &
mesh_cell !< cell connectivity for each element,ip/cell
integer(pInt), dimension(:,:,:), allocatable, private :: &
FE_cellface !< list of intra-cell cell node IDs that constitute the cell faces of a specific type of cell
! These definitions should actually reside in the FE-solver specific part (different for MARC/ABAQUS)
! Hence, I suggest to prefix with "FE_"
integer(pInt), parameter, private :: &
FE_Ngeomtypes = 10_pInt, &
FE_Ncelltypes = 4_pInt, &
FE_maxNmatchingNodesPerFace = 4_pInt, &
FE_maxNfaces = 6_pInt, &
FE_maxNcellnodesPerCell = 8_pInt, &
FE_maxNcellfaces = 6_pInt, &
FE_maxNcellnodesPerCellface = 4_pInt
integer(pInt), dimension(FE_Ncelltypes), parameter, private :: FE_NcellnodesPerCell = & !< number of cell nodes in a specific cell type
int([ &
3, & ! (2D 3node)
4, & ! (2D 4node)
4, & ! (3D 4node)
8 & ! (3D 8node)
],pInt)
integer(pInt), dimension(FE_Ncelltypes), parameter, private :: FE_NcellnodesPerCellface = & !< number of cell nodes per cell face in a specific cell type
int([&
2, & ! (2D 3node)
2, & ! (2D 4node)
3, & ! (3D 4node)
4 & ! (3D 8node)
],pInt)
integer(pInt), dimension(FE_Ncelltypes), parameter, private :: FE_NipNeighbors = & !< number of ip neighbors / cell faces in a specific cell type
int([&
3, & ! (2D 3node)
4, & ! (2D 4node)
4, & ! (3D 4node)
6 & ! (3D 8node)
],pInt)
! grid specific
integer(pInt), dimension(3), public, protected :: & integer(pInt), dimension(3), public, protected :: &
grid !< (global) grid grid !< (global) grid
integer(pInt), public, protected :: & integer(pInt), public, protected :: &
@ -116,18 +64,14 @@ integer(pInt), dimension(:,:), allocatable, private :: &
size3offset !< (local) size offset in 3rd direction size3offset !< (local) size offset in 3rd direction
public :: & public :: &
mesh_init, & mesh_init
mesh_cellCenterCoordinates
private :: & private :: &
mesh_build_cellconnectivity, &
mesh_build_ipAreas, & mesh_build_ipAreas, &
mesh_build_FEdata, & mesh_build_ipNormals, &
mesh_spectral_build_nodes, & mesh_spectral_build_nodes, &
mesh_spectral_build_elements, & mesh_spectral_build_elements, &
mesh_spectral_build_ipNeighborhood, & mesh_spectral_build_ipNeighborhood, &
mesh_build_cellnodes, &
mesh_build_ipVolumes, &
mesh_build_ipCoordinates mesh_build_ipCoordinates
type, public, extends(tMesh) :: tMesh_grid type, public, extends(tMesh) :: tMesh_grid
@ -190,9 +134,8 @@ subroutine mesh_init(ip,el)
implicit none implicit none
include 'fftw3-mpi.f03' include 'fftw3-mpi.f03'
integer(C_INTPTR_T) :: devNull, local_K, local_K_offset integer(C_INTPTR_T) :: devNull, local_K, local_K_offset
integer :: ierr, worldsize integer :: ierr, worldsize, j
integer(pInt), intent(in), optional :: el, ip integer(pInt), intent(in), optional :: el, ip
integer(pInt) :: j
logical :: myDebug logical :: myDebug
write(6,'(/,a)') ' <<<+- mesh init -+>>>' write(6,'(/,a)') ' <<<+- mesh init -+>>>'
@ -225,31 +168,31 @@ subroutine mesh_init(ip,el)
mesh_Nnodes = product(grid(1:2) + 1_pInt)*(grid3 + 1_pInt) mesh_Nnodes = product(grid(1:2) + 1_pInt)*(grid3 + 1_pInt)
call mesh_spectral_build_nodes() mesh_node0 = mesh_spectral_build_nodes()
mesh_node = mesh_node0
if (myDebug) write(6,'(a)') ' Built nodes'; flush(6) if (myDebug) write(6,'(a)') ' Built nodes'; flush(6)
call theMesh%init(mesh_node) call theMesh%init(mesh_node)
call theMesh%setNelems(product(grid(1:2))*grid3) call theMesh%setNelems(product(grid(1:2))*grid3)
mesh_homogenizationAt = mesh_homogenizationAt(product(grid(1:2))*grid3) ! reallocate/shrink in case of MPI
mesh_maxNipNeighbors = theMesh%elem%nIPneighbors
call mesh_spectral_build_elements() call mesh_spectral_build_elements()
mesh_homogenizationAt = mesh_homogenizationAt(product(grid(1:2))*grid3Offset+1: &
product(grid(1:2))*(grid3Offset+grid3)) ! reallocate/shrink in case of MPI
if (myDebug) write(6,'(a)') ' Built elements'; flush(6) if (myDebug) write(6,'(a)') ' Built elements'; flush(6)
call mesh_build_FEdata ! get properties of the different types of elements
call mesh_build_cellconnectivity
if (myDebug) write(6,'(a)') ' Built cell connectivity'; flush(6)
mesh_cellnode = mesh_build_cellnodes(mesh_node,mesh_Ncellnodes)
if (myDebug) write(6,'(a)') ' Built cell nodes'; flush(6) if (myDebug) write(6,'(a)') ' Built cell nodes'; flush(6)
call mesh_build_ipCoordinates mesh_ipCoordinates = mesh_build_ipCoordinates()
if (myDebug) write(6,'(a)') ' Built IP coordinates'; flush(6) if (myDebug) write(6,'(a)') ' Built IP coordinates'; flush(6)
call mesh_build_ipVolumes allocate(mesh_ipVolume(1,theMesh%nElems),source=product([geomSize(1:2),size3]/real([grid(1:2),grid3])))
if (myDebug) write(6,'(a)') ' Built IP volumes'; flush(6) if (myDebug) write(6,'(a)') ' Built IP volumes'; flush(6)
call mesh_build_ipAreas mesh_ipArea = mesh_build_ipAreas()
mesh_ipAreaNormal = mesh_build_ipNormals()
if (myDebug) write(6,'(a)') ' Built IP areas'; flush(6) if (myDebug) write(6,'(a)') ' Built IP areas'; flush(6)
call mesh_spectral_build_ipNeighborhood call mesh_spectral_build_ipNeighborhood
call geometry_plastic_nonlocal_set_IPneighborhood(mesh_ipNeighborhood)
if (myDebug) write(6,'(a)') ' Built IP neighborhood'; flush(6) if (myDebug) write(6,'(a)') ' Built IP neighborhood'; flush(6)
@ -264,13 +207,10 @@ subroutine mesh_init(ip,el)
!!!! COMPATIBILITY HACK !!!! !!!! COMPATIBILITY HACK !!!!
! for a homogeneous mesh, all elements have the same number of IPs and and cell nodes.
! hence, xxPerElem instead of maxXX
! better name
theMesh%homogenizationAt = mesh_element(3,:) theMesh%homogenizationAt = mesh_element(3,:)
theMesh%microstructureAt = mesh_element(4,:) theMesh%microstructureAt = mesh_element(4,:)
!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!
deallocate(mesh_cell)
end subroutine mesh_init end subroutine mesh_init
@ -394,7 +334,7 @@ subroutine mesh_spectral_read_grid()
allocate(mesh_homogenizationAt(product(grid)), source = h) ! too large in case of MPI (shrink later, not very elegant) allocate(mesh_homogenizationAt(product(grid)), source = h) ! too large in case of MPI (shrink later, not very elegant)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! read and interprete content ! read and interpret content
e = 1_pInt e = 1_pInt
do while (startPos < len(rawData)) do while (startPos < len(rawData))
endPos = startPos + index(rawData(startPos:),new_line('')) - 1_pInt endPos = startPos + index(rawData(startPos:),new_line('')) - 1_pInt
@ -429,44 +369,53 @@ subroutine mesh_spectral_read_grid()
end subroutine mesh_spectral_read_grid end subroutine mesh_spectral_read_grid
!-------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
!> @brief Store x,y,z coordinates of all nodes in mesh. !> @brief Calculates position of nodes (pretend to be an element)
!! Allocates global arrays 'mesh_node0' and 'mesh_node' !---------------------------------------------------------------------------------------------------
!-------------------------------------------------------------------------------------------------- pure function mesh_spectral_build_nodes()
subroutine mesh_spectral_build_nodes()
implicit none real(pReal), dimension(3,mesh_Nnodes) :: mesh_spectral_build_nodes
integer(pInt) :: n integer :: n,a,b,c
allocate (mesh_node0 (3,mesh_Nnodes), source = 0.0_pReal) n = 0
do c = 0, grid3
do b = 0, grid(2)
do a = 0, grid(1)
n = n + 1
mesh_spectral_build_nodes(1:3,n) = geomSize/real(grid,pReal) * real([a,b,grid3Offset+c],pReal)
enddo
enddo
enddo
forall (n = 0_pInt:mesh_Nnodes-1_pInt) end function mesh_spectral_build_nodes
mesh_node0(1,n+1_pInt) = mesh_unitlength * &
geomSize(1)*real(mod(n,(grid(1)+1_pInt) ),pReal) &
/ real(grid(1),pReal)
mesh_node0(2,n+1_pInt) = mesh_unitlength * &
geomSize(2)*real(mod(n/(grid(1)+1_pInt),(grid(2)+1_pInt)),pReal) &
/ real(grid(2),pReal)
mesh_node0(3,n+1_pInt) = mesh_unitlength * &
size3*real(mod(n/(grid(1)+1_pInt)/(grid(2)+1_pInt),(grid3+1_pInt)),pReal) &
/ real(grid3,pReal) + &
size3offset
end forall
mesh_node = mesh_node0
end subroutine mesh_spectral_build_nodes !---------------------------------------------------------------------------------------------------
!> @brief Calculates position of IPs/cell centres (pretend to be an element)
!---------------------------------------------------------------------------------------------------
function mesh_build_ipCoordinates()
real(pReal), dimension(3,1,theMesh%nElems) :: mesh_build_ipCoordinates
integer :: n,a,b,c
n = 0
do c = 1, grid3
do b = 1, grid(2)
do a = 1, grid(1)
n = n + 1
mesh_build_ipCoordinates(1:3,1,n) = geomSize/real(grid,pReal) * (real([a,b,grid3Offset+c],pReal) -0.5_pReal)
enddo
enddo
enddo
end function mesh_build_ipCoordinates
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Store FEid, type, material, texture, and node list per element. !> @brief Store FEid, type, material, texture, and node list per element.
!! Allocates global array 'mesh_element' !! Allocates global array 'mesh_element'
!> @todo does the IO_error makes sense?
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine mesh_spectral_build_elements() subroutine mesh_spectral_build_elements()
use IO, only: &
IO_error
implicit none
integer(pInt) :: & integer(pInt) :: &
e, & e, &
elemOffset elemOffset
@ -475,11 +424,9 @@ subroutine mesh_spectral_build_elements()
allocate(mesh_element (4_pInt+8_pInt,theMesh%nElems), source = 0_pInt) allocate(mesh_element (4_pInt+8_pInt,theMesh%nElems), source = 0_pInt)
elemOffset = product(grid(1:2))*grid3Offset elemOffset = product(grid(1:2))*grid3Offset
e = 0_pInt do e=1, theMesh%nElems
do while (e < theMesh%nElems) ! fill expected number of elements, stop at end of data
e = e+1_pInt ! valid element entry
mesh_element( 1,e) = -1_pInt ! DEPRECATED mesh_element( 1,e) = -1_pInt ! DEPRECATED
mesh_element( 2,e) = 10_pInt mesh_element( 2,e) = -1_pInt ! DEPRECATED
mesh_element( 3,e) = mesh_homogenizationAt(e) mesh_element( 3,e) = mesh_homogenizationAt(e)
mesh_element( 4,e) = microGlobal(e+elemOffset) ! microstructure mesh_element( 4,e) = microGlobal(e+elemOffset) ! microstructure
mesh_element( 5,e) = e + (e-1_pInt)/grid(1) + & mesh_element( 5,e) = e + (e-1_pInt)/grid(1) + &
@ -493,8 +440,6 @@ subroutine mesh_spectral_build_elements()
mesh_element(12,e) = mesh_element(9,e) + grid(1) + 1_pInt mesh_element(12,e) = mesh_element(9,e) + grid(1) + 1_pInt
enddo enddo
if (e /= theMesh%nElems) call IO_error(880_pInt,e)
end subroutine mesh_spectral_build_elements end subroutine mesh_spectral_build_elements
@ -508,7 +453,7 @@ subroutine mesh_spectral_build_ipNeighborhood
integer(pInt) :: & integer(pInt) :: &
x,y,z, & x,y,z, &
e e
allocate(mesh_ipNeighborhood(3,theMesh%elem%nIPneighbors,theMesh%elem%nIPs,theMesh%nElems),source=0_pInt) allocate(mesh_ipNeighborhood(3,6,1,theMesh%nElems),source=0_pInt)
e = 0_pInt e = 0_pInt
do z = 0_pInt,grid3-1_pInt do z = 0_pInt,grid3-1_pInt
@ -562,7 +507,6 @@ function mesh_nodesAroundCentres(gDim,Favg,centres) result(nodes)
debug_level, & debug_level, &
debug_levelBasic debug_levelBasic
implicit none
real(pReal), intent(in), dimension(:,:,:,:) :: & real(pReal), intent(in), dimension(:,:,:,:) :: &
centres centres
real(pReal), dimension(3,size(centres,2)+1,size(centres,3)+1,size(centres,4)+1) :: & real(pReal), dimension(3,size(centres,2)+1,size(centres,3)+1,size(centres,4)+1) :: &
@ -641,385 +585,35 @@ function mesh_nodesAroundCentres(gDim,Favg,centres) result(nodes)
end function mesh_nodesAroundCentres end function mesh_nodesAroundCentres
!#################################################################################################################
!#################################################################################################################
!#################################################################################################################
! The following routines are not solver specific and should be included in mesh_base (most likely in modified form)
!#################################################################################################################
!#################################################################################################################
!#################################################################################################################
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Split CP elements into cells. !> @brief calculation of IP interface areas, allocate globals '_ipArea', and '_ipAreaNormal'
!> @details Build a mapping between cells and the corresponding cell nodes ('mesh_cell').
!> Cell nodes that are also matching nodes are unique in the list of cell nodes,
!> all others (currently) might be stored more than once.
!> Also allocates the 'mesh_node' array.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine mesh_build_cellconnectivity pure function mesh_build_ipAreas()
implicit none real(pReal), dimension(6,1,theMesh%nElems) :: mesh_build_ipAreas
integer(pInt), dimension(:), allocatable :: &
matchingNode2cellnode
integer(pInt), dimension(:,:), allocatable :: &
cellnodeParent
integer(pInt), dimension(theMesh%elem%Ncellnodes) :: &
localCellnode2globalCellnode
integer(pInt) :: &
e,n,i, &
matchingNodeID, &
localCellnodeID
integer(pInt), dimension(FE_Ngeomtypes), parameter :: FE_NmatchingNodes = & !< number of nodes that are needed for face matching in a specific type of element geometry
int([ &
3, & ! element 6 (2D 3node 1ip)
3, & ! element 125 (2D 6node 3ip)
4, & ! element 11 (2D 4node 4ip)
4, & ! element 27 (2D 8node 9ip)
4, & ! element 134 (3D 4node 1ip)
4, & ! element 127 (3D 10node 4ip)
6, & ! element 136 (3D 6node 6ip)
8, & ! element 117 (3D 8node 1ip)
8, & ! element 7 (3D 8node 8ip)
8 & ! element 21 (3D 20node 27ip)
],pInt)
allocate(mesh_cell(FE_maxNcellnodesPerCell,theMesh%elem%nIPs,theMesh%nElems), source=0_pInt) mesh_build_ipAreas(1:2,1,:) = geomSize(2)/real(grid(2)) * geomSize(3)/real(grid(3))
allocate(matchingNode2cellnode(theMesh%nNodes), source=0_pInt) mesh_build_ipAreas(3:4,1,:) = geomSize(3)/real(grid(3)) * geomSize(1)/real(grid(1))
allocate(cellnodeParent(2_pInt,theMesh%elem%Ncellnodes*theMesh%nElems), source=0_pInt) mesh_build_ipAreas(5:6,1,:) = geomSize(1)/real(grid(1)) * geomSize(2)/real(grid(2))
mesh_Ncells = theMesh%nElems*theMesh%elem%nIPs end function mesh_build_ipAreas
!--------------------------------------------------------------------------------------------------
! Count cell nodes (including duplicates) and generate cell connectivity list
mesh_Ncellnodes = 0_pInt
do e = 1_pInt,theMesh%nElems
localCellnode2globalCellnode = 0_pInt
do i = 1_pInt,theMesh%elem%nIPs
do n = 1_pInt,theMesh%elem%NcellnodesPerCell
localCellnodeID = theMesh%elem%cell(n,i)
if (localCellnodeID <= FE_NmatchingNodes(theMesh%elem%geomType)) then ! this cell node is a matching node
matchingNodeID = mesh_element(4_pInt+localCellnodeID,e)
if (matchingNode2cellnode(matchingNodeID) == 0_pInt) then ! if this matching node does not yet exist in the glbal cell node list ...
mesh_Ncellnodes = mesh_Ncellnodes + 1_pInt ! ... count it as cell node ...
matchingNode2cellnode(matchingNodeID) = mesh_Ncellnodes ! ... and remember its global ID
cellnodeParent(1_pInt,mesh_Ncellnodes) = e ! ... and where it belongs to
cellnodeParent(2_pInt,mesh_Ncellnodes) = localCellnodeID
endif
mesh_cell(n,i,e) = matchingNode2cellnode(matchingNodeID)
else ! this cell node is no matching node
if (localCellnode2globalCellnode(localCellnodeID) == 0_pInt) then ! if this local cell node does not yet exist in the global cell node list ...
mesh_Ncellnodes = mesh_Ncellnodes + 1_pInt ! ... count it as cell node ...
localCellnode2globalCellnode(localCellnodeID) = mesh_Ncellnodes ! ... and remember its global ID ...
cellnodeParent(1_pInt,mesh_Ncellnodes) = e ! ... and it belongs to
cellnodeParent(2_pInt,mesh_Ncellnodes) = localCellnodeID
endif
mesh_cell(n,i,e) = localCellnode2globalCellnode(localCellnodeID)
endif
enddo
enddo
enddo
allocate(mesh_cellnodeParent(2_pInt,mesh_Ncellnodes))
allocate(mesh_cellnode(3_pInt,mesh_Ncellnodes))
forall(n = 1_pInt:mesh_Ncellnodes)
mesh_cellnodeParent(1,n) = cellnodeParent(1,n)
mesh_cellnodeParent(2,n) = cellnodeParent(2,n)
endforall
end subroutine mesh_build_cellconnectivity
!--------------------------------------------------------------------------------------------------
!> @brief Calculate position of cellnodes from the given position of nodes
!> Build list of cellnodes' coordinates.
!> Cellnode coordinates are calculated from a weighted sum of node coordinates.
!--------------------------------------------------------------------------------------------------
function mesh_build_cellnodes(nodes,Ncellnodes)
implicit none
integer(pInt), intent(in) :: Ncellnodes !< requested number of cellnodes
real(pReal), dimension(3,mesh_Nnodes), intent(in) :: nodes
real(pReal), dimension(3,Ncellnodes) :: mesh_build_cellnodes
integer(pInt) :: &
e,n,m, &
localCellnodeID
real(pReal), dimension(3) :: &
myCoords
mesh_build_cellnodes = 0.0_pReal
!$OMP PARALLEL DO PRIVATE(e,localCellnodeID,myCoords)
do n = 1_pInt,Ncellnodes ! loop over cell nodes
e = mesh_cellnodeParent(1,n)
localCellnodeID = mesh_cellnodeParent(2,n)
myCoords = 0.0_pReal
do m = 1_pInt,theMesh%elem%nNodes
myCoords = myCoords + nodes(1:3,mesh_element(4_pInt+m,e)) &
* theMesh%elem%cellNodeParentNodeWeights(m,localCellnodeID)
enddo
mesh_build_cellnodes(1:3,n) = myCoords / sum(theMesh%elem%cellNodeParentNodeWeights(:,localCellnodeID))
enddo
!$OMP END PARALLEL DO
end function mesh_build_cellnodes
!--------------------------------------------------------------------------------------------------
!> @brief Calculates IP volume. Allocates global array 'mesh_ipVolume'
!> @details The IP volume is calculated differently depending on the cell type.
!> 2D cells assume an element depth of one in order to calculate the volume.
!> For the hexahedral cell we subdivide the cell into subvolumes of pyramidal
!> shape with a cell face as basis and the central ip at the tip. This subvolume is
!> calculated as an average of four tetrahedals with three corners on the cell face
!> and one corner at the central ip.
!--------------------------------------------------------------------------------------------------
subroutine mesh_build_ipVolumes
use math, only: &
math_volTetrahedron, &
math_areaTriangle
implicit none
integer(pInt) :: e,t,g,c,i,m,f,n
real(pReal), dimension(FE_maxNcellnodesPerCellface,FE_maxNcellfaces) :: subvolume
allocate(mesh_ipVolume(theMesh%elem%nIPs,theMesh%nElems),source=0.0_pReal)
!$OMP PARALLEL DO PRIVATE(t,g,c,m,subvolume)
do e = 1_pInt,theMesh%nElems ! loop over cpElems
select case (theMesh%elem%cellType)
case (1_pInt) ! 2D 3node
forall (i = 1_pInt:theMesh%elem%nIPs) & ! loop over ips=cells in this element
mesh_ipVolume(i,e) = math_areaTriangle(mesh_cellnode(1:3,mesh_cell(1,i,e)), &
mesh_cellnode(1:3,mesh_cell(2,i,e)), &
mesh_cellnode(1:3,mesh_cell(3,i,e)))
case (2_pInt) ! 2D 4node
forall (i = 1_pInt:theMesh%elem%nIPs) & ! loop over ips=cells in this element
mesh_ipVolume(i,e) = math_areaTriangle(mesh_cellnode(1:3,mesh_cell(1,i,e)), & ! here we assume a planar shape, so division in two triangles suffices
mesh_cellnode(1:3,mesh_cell(2,i,e)), &
mesh_cellnode(1:3,mesh_cell(3,i,e))) &
+ math_areaTriangle(mesh_cellnode(1:3,mesh_cell(3,i,e)), &
mesh_cellnode(1:3,mesh_cell(4,i,e)), &
mesh_cellnode(1:3,mesh_cell(1,i,e)))
case (3_pInt) ! 3D 4node
forall (i = 1_pInt:theMesh%elem%nIPs) & ! loop over ips=cells in this element
mesh_ipVolume(i,e) = math_volTetrahedron(mesh_cellnode(1:3,mesh_cell(1,i,e)), &
mesh_cellnode(1:3,mesh_cell(2,i,e)), &
mesh_cellnode(1:3,mesh_cell(3,i,e)), &
mesh_cellnode(1:3,mesh_cell(4,i,e)))
case (4_pInt)
c = theMesh%elem%cellType ! 3D 8node
m = FE_NcellnodesPerCellface(c)
do i = 1_pInt,theMesh%elem%nIPs ! loop over ips=cells in this element
subvolume = 0.0_pReal
forall(f = 1_pInt:FE_NipNeighbors(c), n = 1_pInt:FE_NcellnodesPerCellface(c)) &
subvolume(n,f) = math_volTetrahedron(&
mesh_cellnode(1:3,mesh_cell(FE_cellface( n ,f,c),i,e)), &
mesh_cellnode(1:3,mesh_cell(FE_cellface(1+mod(n ,m),f,c),i,e)), &
mesh_cellnode(1:3,mesh_cell(FE_cellface(1+mod(n+1,m),f,c),i,e)), &
mesh_ipCoordinates(1:3,i,e))
mesh_ipVolume(i,e) = 0.5_pReal * sum(subvolume) ! each subvolume is based on four tetrahedrons, altough the face consists of only two triangles -> averaging factor two
enddo
end select
enddo
!$OMP END PARALLEL DO
end subroutine mesh_build_ipVolumes
!--------------------------------------------------------------------------------------------------
!> @brief Calculates IP Coordinates. Allocates global array 'mesh_ipCoordinates'
! Called by all solvers in mesh_init in order to initialize the ip coordinates.
! Later on the current ip coordinates are directly prvided by the spectral solver and by Abaqus,
! so no need to use this subroutine anymore; Marc however only provides nodal displacements,
! so in this case the ip coordinates are always calculated on the basis of this subroutine.
! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! FOR THE MOMENT THIS SUBROUTINE ACTUALLY CALCULATES THE CELL CENTER AND NOT THE IP COORDINATES,
! AS THE IP IS NOT (ALWAYS) LOCATED IN THE CENTER OF THE IP VOLUME.
! HAS TO BE CHANGED IN A LATER VERSION.
! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!--------------------------------------------------------------------------------------------------
subroutine mesh_build_ipCoordinates
implicit none
integer(pInt) :: e,c,i,n
real(pReal), dimension(3) :: myCoords
if (.not. allocated(mesh_ipCoordinates)) &
allocate(mesh_ipCoordinates(3,theMesh%elem%nIPs,theMesh%nElems),source=0.0_pReal)
!$OMP PARALLEL DO PRIVATE(c,myCoords)
do e = 1_pInt,theMesh%nElems ! loop over cpElems
c = theMesh%elem%cellType
do i = 1_pInt,theMesh%elem%nIPs ! loop over ips=cells in this element
myCoords = 0.0_pReal
do n = 1_pInt,FE_NcellnodesPerCell(c) ! loop over cell nodes in this cell
myCoords = myCoords + mesh_cellnode(1:3,mesh_cell(n,i,e))
enddo
mesh_ipCoordinates(1:3,i,e) = myCoords / real(FE_NcellnodesPerCell(c),pReal)
enddo
enddo
!$OMP END PARALLEL DO
end subroutine mesh_build_ipCoordinates
!--------------------------------------------------------------------------------------------------
!> @brief Calculates cell center coordinates.
!--------------------------------------------------------------------------------------------------
pure function mesh_cellCenterCoordinates(ip,el)
implicit none
integer(pInt), intent(in) :: el, & !< element number
ip !< integration point number
real(pReal), dimension(3) :: mesh_cellCenterCoordinates !< x,y,z coordinates of the cell center of the requested IP cell
integer(pInt) :: c,n
c = theMesh%elem%cellType
mesh_cellCenterCoordinates = 0.0_pReal
do n = 1_pInt,FE_NcellnodesPerCell(c) ! loop over cell nodes in this cell
mesh_cellCenterCoordinates = mesh_cellCenterCoordinates + mesh_cellnode(1:3,mesh_cell(n,ip,el))
enddo
mesh_cellCenterCoordinates = mesh_cellCenterCoordinates / real(FE_NcellnodesPerCell(c),pReal)
end function mesh_cellCenterCoordinates
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculation of IP interface areas, allocate globals '_ipArea', and '_ipAreaNormal' !> @brief calculation of IP interface areas, allocate globals '_ipArea', and '_ipAreaNormal'
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine mesh_build_ipAreas pure function mesh_build_ipNormals()
use math, only: &
math_cross
implicit none real, dimension(3,6,1,theMesh%nElems) :: mesh_build_ipNormals
integer(pInt) :: e,t,g,c,i,f,n,m
real(pReal), dimension (3,FE_maxNcellnodesPerCellface) :: nodePos, normals
real(pReal), dimension(3) :: normal
allocate(mesh_ipArea(theMesh%elem%nIPneighbors,theMesh%elem%nIPs,theMesh%nElems), source=0.0_pReal) mesh_build_ipNormals(1:3,1,1,:) = spread([+1.0_pReal, 0.0_pReal, 0.0_pReal],2,theMesh%nElems)
allocate(mesh_ipAreaNormal(3_pInt,theMesh%elem%nIPneighbors,theMesh%elem%nIPs,theMesh%nElems), source=0.0_pReal) mesh_build_ipNormals(1:3,2,1,:) = spread([-1.0_pReal, 0.0_pReal, 0.0_pReal],2,theMesh%nElems)
mesh_build_ipNormals(1:3,3,1,:) = spread([ 0.0_pReal,+1.0_pReal, 0.0_pReal],2,theMesh%nElems)
!$OMP PARALLEL DO PRIVATE(t,g,c,nodePos,normal,normals) mesh_build_ipNormals(1:3,4,1,:) = spread([ 0.0_pReal,-1.0_pReal, 0.0_pReal],2,theMesh%nElems)
do e = 1_pInt,theMesh%nElems ! loop over cpElems mesh_build_ipNormals(1:3,5,1,:) = spread([ 0.0_pReal, 0.0_pReal,+1.0_pReal],2,theMesh%nElems)
c = theMesh%elem%cellType mesh_build_ipNormals(1:3,6,1,:) = spread([ 0.0_pReal, 0.0_pReal,-1.0_pReal],2,theMesh%nElems)
select case (c)
case (1_pInt,2_pInt) ! 2D 3 or 4 node
do i = 1_pInt,theMesh%elem%nIPs ! loop over ips=cells in this element
do f = 1_pInt,FE_NipNeighbors(c) ! loop over cell faces
forall(n = 1_pInt:FE_NcellnodesPerCellface(c)) &
nodePos(1:3,n) = mesh_cellnode(1:3,mesh_cell(FE_cellface(n,f,c),i,e))
normal(1) = nodePos(2,2) - nodePos(2,1) ! x_normal = y_connectingVector
normal(2) = -(nodePos(1,2) - nodePos(1,1)) ! y_normal = -x_connectingVector
normal(3) = 0.0_pReal
mesh_ipArea(f,i,e) = norm2(normal)
mesh_ipAreaNormal(1:3,f,i,e) = normal / norm2(normal) ! ensure unit length of area normal
enddo
enddo
case (3_pInt) ! 3D 4node
do i = 1_pInt,theMesh%elem%nIPs ! loop over ips=cells in this element
do f = 1_pInt,FE_NipNeighbors(c) ! loop over cell faces
forall(n = 1_pInt:FE_NcellnodesPerCellface(c)) &
nodePos(1:3,n) = mesh_cellnode(1:3,mesh_cell(FE_cellface(n,f,c),i,e))
normal = math_cross(nodePos(1:3,2) - nodePos(1:3,1), &
nodePos(1:3,3) - nodePos(1:3,1))
mesh_ipArea(f,i,e) = norm2(normal)
mesh_ipAreaNormal(1:3,f,i,e) = normal / norm2(normal) ! ensure unit length of area normal
enddo
enddo
case (4_pInt) ! 3D 8node
! for this cell type we get the normal of the quadrilateral face as an average of
! four normals of triangular subfaces; since the face consists only of two triangles,
! the sum has to be divided by two; this whole prcedure tries to compensate for
! probable non-planar cell surfaces
m = FE_NcellnodesPerCellface(c)
do i = 1_pInt,theMesh%elem%nIPs ! loop over ips=cells in this element
do f = 1_pInt,FE_NipNeighbors(c) ! loop over cell faces
forall(n = 1_pInt:FE_NcellnodesPerCellface(c)) &
nodePos(1:3,n) = mesh_cellnode(1:3,mesh_cell(FE_cellface(n,f,c),i,e))
forall(n = 1_pInt:FE_NcellnodesPerCellface(c)) &
normals(1:3,n) = 0.5_pReal &
* math_cross(nodePos(1:3,1+mod(n ,m)) - nodePos(1:3,n), &
nodePos(1:3,1+mod(n+1,m)) - nodePos(1:3,n))
normal = 0.5_pReal * sum(normals,2)
mesh_ipArea(f,i,e) = norm2(normal)
mesh_ipAreaNormal(1:3,f,i,e) = normal / norm2(normal)
enddo
enddo
end select
enddo
!$OMP END PARALLEL DO
end subroutine mesh_build_ipAreas
!--------------------------------------------------------------------------------------------------
!> @brief get properties of different types of finite elements
!> @details assign globals: FE_nodesAtIP, FE_ipNeighbor, FE_subNodeOnIPFace
!--------------------------------------------------------------------------------------------------
subroutine mesh_build_FEdata
implicit none
integer(pInt) :: me
allocate(FE_cellface(FE_maxNcellnodesPerCellface,FE_maxNcellfaces,FE_Ncelltypes), source=0_pInt)
! *** FE_cellface ***
me = 0_pInt
me = me + 1_pInt
FE_cellface(1:FE_NcellnodesPerCellface(me),1:FE_NipNeighbors(me),me) = & ! 2D 3node, VTK_TRIANGLE (5)
reshape(int([&
2,3, &
3,1, &
1,2 &
],pInt),[FE_NcellnodesPerCellface(me),FE_NipNeighbors(me)])
me = me + 1_pInt
FE_cellface(1:FE_NcellnodesPerCellface(me),1:FE_NipNeighbors(me),me) = & ! 2D 4node, VTK_QUAD (9)
reshape(int([&
2,3, &
4,1, &
3,4, &
1,2 &
],pInt),[FE_NcellnodesPerCellface(me),FE_NipNeighbors(me)])
me = me + 1_pInt
FE_cellface(1:FE_NcellnodesPerCellface(me),1:FE_NipNeighbors(me),me) = & ! 3D 4node, VTK_TETRA (10)
reshape(int([&
1,3,2, &
1,2,4, &
2,3,4, &
1,4,3 &
],pInt),[FE_NcellnodesPerCellface(me),FE_NipNeighbors(me)])
me = me + 1_pInt
FE_cellface(1:FE_NcellnodesPerCellface(me),1:FE_NipNeighbors(me),me) = & ! 3D 8node, VTK_HEXAHEDRON (12)
reshape(int([&
2,3,7,6, &
4,1,5,8, &
3,4,8,7, &
1,2,6,5, &
5,6,7,8, &
1,4,3,2 &
],pInt),[FE_NcellnodesPerCellface(me),FE_NipNeighbors(me)])
end function mesh_build_ipNormals
end subroutine mesh_build_FEdata
end module mesh end module mesh

103
src/plastic_dislotwin.f90
View File

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

63
src/plastic_isotropic.f90
View File

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

62
src/plastic_kinematichardening.f90
View File

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

7
src/plastic_none.f90
View File

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

11
src/plastic_nonlocal.f90
View File

@ -5,10 +5,15 @@
!> @brief material subroutine for plasticity including dislocation flux !> @brief material subroutine for plasticity including dislocation flux
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module plastic_nonlocal module plastic_nonlocal
use prec, only: & use prec
pReal
use future use future
use geometry_plastic_nonlocal, only: &
periodicSurface => geometry_plastic_nonlocal_periodicSurface, &
IPneighborhood => geometry_plastic_nonlocal_IPneighborhood, &
IPvolume => geometry_plastic_nonlocal_IPvolume, &
IParea => geometry_plastic_nonlocal_IParea, &
IPareaNormal => geometry_plastic_nonlocal_IPareaNormal
implicit none implicit none
private private
real(pReal), parameter, private :: & real(pReal), parameter, private :: &

50
src/plastic_phenopowerlaw.f90
View File

@ -5,11 +5,20 @@
!> @brief phenomenological crystal plasticity formulation using a powerlaw fitting !> @brief phenomenological crystal plasticity formulation using a powerlaw fitting
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module plastic_phenopowerlaw module plastic_phenopowerlaw
use prec, only: & use prec
pReal use debug
use math
use IO
use material
use config
use lattice
#if defined(PETSc) || defined(DAMASK_HDF5)
use results
#endif
implicit none implicit none
private private
integer, dimension(:,:), allocatable, target, public :: & integer, dimension(:,:), allocatable, target, public :: &
plastic_phenopowerlaw_sizePostResult !< size of each post result output plastic_phenopowerlaw_sizePostResult !< size of each post result output
character(len=64), dimension(:,:), allocatable, target, public :: & character(len=64), dimension(:,:), allocatable, target, public :: &
@ -28,7 +37,7 @@ module plastic_phenopowerlaw
resolvedstress_twin_ID resolvedstress_twin_ID
end enum end enum
type, private :: tParameters type :: tParameters
real(pReal) :: & real(pReal) :: &
gdot0_slip, & !< reference shear strain rate for slip gdot0_slip, & !< reference shear strain rate for slip
gdot0_twin, & !< reference shear strain rate for twin gdot0_twin, & !< reference shear strain rate for twin
@ -73,7 +82,7 @@ module plastic_phenopowerlaw
outputID !< ID of each post result output outputID !< ID of each post result output
end type tParameters end type tParameters
type, private :: tPhenopowerlawState type :: tPhenopowerlawState
real(pReal), pointer, dimension(:,:) :: & real(pReal), pointer, dimension(:,:) :: &
xi_slip, & xi_slip, &
xi_twin, & xi_twin, &
@ -83,8 +92,8 @@ module plastic_phenopowerlaw
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! containers for parameters and state ! containers for parameters and state
type(tParameters), allocatable, dimension(:), private :: param type(tParameters), allocatable, dimension(:) :: param
type(tPhenopowerlawState), allocatable, dimension(:), private :: & type(tPhenopowerlawState), allocatable, dimension(:) :: &
dotState, & dotState, &
state state
@ -94,9 +103,6 @@ module plastic_phenopowerlaw
plastic_phenopowerlaw_dotState, & plastic_phenopowerlaw_dotState, &
plastic_phenopowerlaw_postResults, & plastic_phenopowerlaw_postResults, &
plastic_phenopowerlaw_results plastic_phenopowerlaw_results
private :: &
kinetics_slip, &
kinetics_twin
contains contains
@ -106,20 +112,6 @@ contains
!> @details reads in material parameters, allocates arrays, and does sanity checks !> @details reads in material parameters, allocates arrays, and does sanity checks
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine plastic_phenopowerlaw_init 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 :: & integer :: &
Ninstance, & Ninstance, &
@ -484,8 +476,6 @@ end subroutine plastic_phenopowerlaw_dotState
!> @brief return array of constitutive results !> @brief return array of constitutive results
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function plastic_phenopowerlaw_postResults(Mp,instance,of) result(postResults) function plastic_phenopowerlaw_postResults(Mp,instance,of) result(postResults)
use math, only: &
math_mul33xx33
real(pReal), dimension(3,3), intent(in) :: & real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress Mp !< Mandel stress
@ -552,8 +542,6 @@ end function plastic_phenopowerlaw_postResults
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine plastic_phenopowerlaw_results(instance,group) subroutine plastic_phenopowerlaw_results(instance,group)
#if defined(PETSc) || defined(DAMASK_HDF5) #if defined(PETSc) || defined(DAMASK_HDF5)
use results, only: &
results_writeDataset
integer, intent(in) :: instance integer, intent(in) :: instance
character(len=*), intent(in) :: group character(len=*), intent(in) :: group
@ -598,10 +586,6 @@ end subroutine plastic_phenopowerlaw_results
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure subroutine kinetics_slip(Mp,instance,of, & pure subroutine kinetics_slip(Mp,instance,of, &
gdot_slip_pos,gdot_slip_neg,dgdot_dtau_slip_pos,dgdot_dtau_slip_neg) 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) :: & real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress Mp !< Mandel stress
@ -674,10 +658,6 @@ end subroutine kinetics_slip
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure subroutine kinetics_twin(Mp,instance,of,& pure subroutine kinetics_twin(Mp,instance,of,&
gdot_twin,dgdot_dtau_twin) gdot_twin,dgdot_dtau_twin)
use prec, only: &
dNeq0
use math, only: &
math_mul33xx33
real(pReal), dimension(3,3), intent(in) :: & real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress 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 ! Modified 2017-2019, Martin Diehl/Max-Planck-Institut für Eisenforschung GmbH
! All rights reserved. ! 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: ! 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. ! of conditions and the following disclaimer.
! - Redistributions in binary form must reproduce the above copyright notice, this ! - Redistributions in binary form must reproduce the above copyright notice, this
! list of conditions and the following disclaimer in the documentation and/or ! list of conditions and the following disclaimer in the documentation and/or
! other materials provided with the distribution. ! other materials provided with the distribution.
! - Neither the names of Marc De Graef, Carnegie Mellon University nor the names ! - 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 ! of its contributors may be used to endorse or promote products derived from
! this software without specific prior written permission. ! this software without specific prior written permission.
! !
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE ! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
! IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ! IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
! ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE ! ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
! LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL ! LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR ! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
! SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER ! SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, ! 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 ! 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. ! 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. !> @details w is the real part, (x, y, z) are the imaginary parts.
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
module quaternions module quaternions
use prec, only: & use prec
pReal use future
use future
implicit none implicit none
public 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
contains real(pReal), parameter, public :: P = -1.0_pReal !< parameter for orientation conversion.
procedure, private :: add__
procedure, private :: pos__
generic, public :: operator(+) => add__,pos__
procedure, private :: sub__ type, public :: quaternion
procedure, private :: neg__ real(pReal) :: w = 0.0_pReal
generic, public :: operator(-) => sub__,neg__ 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__ contains
procedure, private :: div_scal__ procedure, private :: add__
generic, public :: operator(/) => div_quat__, div_scal__ procedure, private :: pos__
generic, public :: operator(+) => add__,pos__
procedure, private :: eq__ procedure, private :: sub__
generic, public :: operator(==) => eq__ procedure, private :: neg__
generic, public :: operator(-) => sub__,neg__
procedure, private :: neq__ procedure, private :: mul_quat__
generic, public :: operator(/=) => neq__ procedure, private :: mul_scal__
generic, public :: operator(*) => mul_quat__, mul_scal__
procedure, private :: pow_quat__ procedure, private :: div_quat__
procedure, private :: pow_scal__ procedure, private :: div_scal__
generic, public :: operator(**) => pow_quat__, pow_scal__ generic, public :: operator(/) => div_quat__, div_scal__
procedure, public :: abs__ procedure, private :: eq__
procedure, public :: dot_product__ generic, public :: operator(==) => eq__
procedure, public :: conjg__
procedure, public :: exp__
procedure, public :: log__
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 (=) interface assignment (=)
module procedure assign_quat__ module procedure assign_quat__
@ -124,12 +123,12 @@ contains
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
type(quaternion) pure function init__(array) type(quaternion) pure function init__(array)
real(pReal), intent(in), dimension(4) :: array real(pReal), intent(in), dimension(4) :: array
init__%w=array(1) init__%w=array(1)
init__%x=array(2) init__%x=array(2)
init__%y=array(3) init__%y=array(3)
init__%z=array(4) init__%z=array(4)
end function init__ end function init__
@ -139,14 +138,14 @@ end function init__
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
elemental subroutine assign_quat__(self,other) elemental subroutine assign_quat__(self,other)
type(quaternion), intent(out) :: self type(quaternion), intent(out) :: self
type(quaternion), intent(in) :: other type(quaternion), intent(in) :: other
self%w = other%w self%w = other%w
self%x = other%x self%x = other%x
self%y = other%y self%y = other%y
self%z = other%z self%z = other%z
end subroutine assign_quat__ end subroutine assign_quat__
@ -155,14 +154,14 @@ end subroutine assign_quat__
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
pure subroutine assign_vec__(self,other) pure subroutine assign_vec__(self,other)
type(quaternion), intent(out) :: self type(quaternion), intent(out) :: self
real(pReal), intent(in), dimension(4) :: other real(pReal), intent(in), dimension(4) :: other
self%w = other(1) self%w = other(1)
self%x = other(2) self%x = other(2)
self%y = other(3) self%y = other(3)
self%z = other(4) self%z = other(4)
end subroutine assign_vec__ end subroutine assign_vec__
@ -171,13 +170,13 @@ end subroutine assign_vec__
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
type(quaternion) elemental function add__(self,other) type(quaternion) elemental function add__(self,other)
class(quaternion), intent(in) :: self,other class(quaternion), intent(in) :: self,other
add__%w = self%w + other%w add__%w = self%w + other%w
add__%x = self%x + other%x add__%x = self%x + other%x
add__%y = self%y + other%y add__%y = self%y + other%y
add__%z = self%z + other%z add__%z = self%z + other%z
end function add__ end function add__
@ -186,13 +185,13 @@ end function add__
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
type(quaternion) elemental function pos__(self) type(quaternion) elemental function pos__(self)
class(quaternion), intent(in) :: self class(quaternion), intent(in) :: self
pos__%w = self%w pos__%w = self%w
pos__%x = self%x pos__%x = self%x
pos__%y = self%y pos__%y = self%y
pos__%z = self%z pos__%z = self%z
end function pos__ end function pos__
@ -201,13 +200,13 @@ end function pos__
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
type(quaternion) elemental function sub__(self,other) type(quaternion) elemental function sub__(self,other)
class(quaternion), intent(in) :: self,other class(quaternion), intent(in) :: self,other
sub__%w = self%w - other%w sub__%w = self%w - other%w
sub__%x = self%x - other%x sub__%x = self%x - other%x
sub__%y = self%y - other%y sub__%y = self%y - other%y
sub__%z = self%z - other%z sub__%z = self%z - other%z
end function sub__ end function sub__
@ -216,13 +215,13 @@ end function sub__
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
type(quaternion) elemental function neg__(self) type(quaternion) elemental function neg__(self)
class(quaternion), intent(in) :: self class(quaternion), intent(in) :: self
neg__%w = -self%w neg__%w = -self%w
neg__%x = -self%x neg__%x = -self%x
neg__%y = -self%y neg__%y = -self%y
neg__%z = -self%z neg__%z = -self%z
end function neg__ end function neg__
@ -231,13 +230,13 @@ end function neg__
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
type(quaternion) elemental function mul_quat__(self,other) 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__ end function mul_quat__
@ -246,14 +245,14 @@ end function mul_quat__
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
type(quaternion) elemental function mul_scal__(self,scal) type(quaternion) elemental function mul_scal__(self,scal)
class(quaternion), intent(in) :: self class(quaternion), intent(in) :: self
real(pReal), intent(in) :: scal 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__ end function mul_scal__
@ -262,9 +261,9 @@ end function mul_scal__
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
type(quaternion) elemental function div_quat__(self,other) 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__ end function div_quat__
@ -274,10 +273,10 @@ end function div_quat__
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
type(quaternion) elemental function div_scal__(self,scal) type(quaternion) elemental function div_scal__(self,scal)
class(quaternion), intent(in) :: self class(quaternion), intent(in) :: self
real(pReal), intent(in) :: scal 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__ end function div_scal__
@ -286,14 +285,12 @@ end function div_scal__
!> equality of two quaternions !> equality of two quaternions
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
logical elemental function eq__(self,other) 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__ end function eq__
@ -302,10 +299,10 @@ end function eq__
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
logical elemental function neq__(self,other) 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__ end function neq__
@ -314,11 +311,11 @@ end function neq__
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
type(quaternion) elemental function pow_scal__(self,expon) type(quaternion) elemental function pow_scal__(self,expon)
class(quaternion), intent(in) :: self class(quaternion), intent(in) :: self
real(pReal), intent(in) :: expon real(pReal), intent(in) :: expon
pow_scal__ = exp(log(self)*expon) pow_scal__ = exp(log(self)*expon)
end function pow_scal__ end function pow_scal__
@ -327,11 +324,11 @@ end function pow_scal__
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
type(quaternion) elemental function pow_quat__(self,expon) type(quaternion) elemental function pow_quat__(self,expon)
class(quaternion), intent(in) :: self class(quaternion), intent(in) :: self
type(quaternion), intent(in) :: expon type(quaternion), intent(in) :: expon
pow_quat__ = exp(log(self)*expon) pow_quat__ = exp(log(self)*expon)
end function pow_quat__ end function pow_quat__
@ -341,15 +338,15 @@ end function pow_quat__
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
type(quaternion) elemental function exp__(self) type(quaternion) elemental function exp__(self)
class(quaternion), intent(in) :: self class(quaternion), intent(in) :: self
real(pReal) :: absImag 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), & exp__ = exp(self%w) * [ cos(absImag), &
self%x/absImag * sin(absImag), & self%x/absImag * sin(absImag), &
self%y/absImag * sin(absImag), & self%y/absImag * sin(absImag), &
self%z/absImag * sin(absImag)] self%z/absImag * sin(absImag)]
end function exp__ end function exp__
@ -360,16 +357,16 @@ end function exp__
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
type(quaternion) elemental function log__(self) type(quaternion) elemental function log__(self)
class(quaternion), intent(in) :: self class(quaternion), intent(in) :: self
real(pReal) :: absImag 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__ end function log__
@ -378,10 +375,10 @@ end function log__
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
real(pReal) elemental function abs__(a) 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__ end function abs__
@ -390,10 +387,10 @@ end function abs__
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
real(pReal) elemental function dot_product__(a,b) 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__ end function dot_product__
@ -402,10 +399,10 @@ end function dot_product__
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
type(quaternion) elemental function conjg__(a) 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__ end function conjg__
@ -414,10 +411,10 @@ end function conjg__
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
type(quaternion) elemental function quat_homomorphed(a) 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 function quat_homomorphed
end module quaternions end module quaternions

16
src/results.f90
View File

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

86
src/rotations.f90
View File

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

465
src/source_damage_anisoBrittle.f90
View File

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

31
src/source_thermal_externalheat.f90
View File

@ -5,11 +5,14 @@
!> @brief material subroutine for variable heat source !> @brief material subroutine for variable heat source
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module source_thermal_externalheat module source_thermal_externalheat
use prec, only: & use prec
pReal use debug
use material
use config
implicit none implicit none
private private
integer, dimension(:), allocatable, public, protected :: & integer, dimension(:), allocatable, public, protected :: &
source_thermal_externalheat_offset, & !< which source is my current thermal dissipation mechanism? source_thermal_externalheat_offset, & !< which source is my current thermal dissipation mechanism?
source_thermal_externalheat_instance !< instance of thermal dissipation source mechanism source_thermal_externalheat_instance !< instance of thermal dissipation source mechanism
@ -23,7 +26,7 @@ module source_thermal_externalheat
integer, dimension(:), allocatable, target, public :: & integer, dimension(:), allocatable, target, public :: &
source_thermal_externalheat_Noutput !< number of outputs per instance of this source 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 :: & real(pReal), dimension(:), allocatable :: &
time, & time, &
heat_rate heat_rate
@ -31,7 +34,7 @@ module source_thermal_externalheat
nIntervals nIntervals
end type tParameters 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 :: & public :: &
@ -47,22 +50,6 @@ contains
!> @details reads in material parameters, allocates arrays, and does sanity checks !> @details reads in material parameters, allocates arrays, and does sanity checks
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine source_thermal_externalheat_init 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 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 !> @details state only contains current time to linearly interpolate given heat powers
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine source_thermal_externalheat_dotState(phase, of) subroutine source_thermal_externalheat_dotState(phase, of)
use material, only: &
sourceState
integer, intent(in) :: & integer, intent(in) :: &
phase, & phase, &
@ -135,8 +120,6 @@ end subroutine source_thermal_externalheat_dotState
!> @brief returns local heat generation rate !> @brief returns local heat generation rate
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, dTDot_dT, phase, of) subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, dTDot_dT, phase, of)
use material, only: &
sourceState
integer, intent(in) :: & integer, intent(in) :: &
phase, & phase, &

78
src/thermal_adiabatic.f90
View File

@ -3,9 +3,16 @@
!> @brief material subroutine for adiabatic temperature evolution !> @brief material subroutine for adiabatic temperature evolution
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module thermal_adiabatic module thermal_adiabatic
use prec, only: & use prec
pReal use config
use numerics
use material
use source_thermal_dissipation
use source_thermal_externalheat
use crystallite
use lattice
use mesh
implicit none implicit none
private private
@ -21,7 +28,7 @@ module thermal_adiabatic
enumerator :: undefined_ID, & enumerator :: undefined_ID, &
temperature_ID temperature_ID
end enum end enum
integer(kind(undefined_ID)), dimension(:,:), allocatable, private :: & integer(kind(undefined_ID)), dimension(:,:), allocatable :: &
thermal_adiabatic_outputID !< ID of each post result output 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 !> @details reads in material parameters, allocates arrays, and does sanity checks
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine thermal_adiabatic_init 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 integer :: maxNinstance,section,instance,i,sizeState,NofMyHomog
character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::] 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 !> @brief calculates adiabatic change in temperature based on local heat generation model
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function thermal_adiabatic_updateState(subdt, ip, el) 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) :: & integer, intent(in) :: &
ip, & !< integration point number ip, & !< integration point number
@ -156,28 +138,11 @@ function thermal_adiabatic_updateState(subdt, ip, el)
end function thermal_adiabatic_updateState end function thermal_adiabatic_updateState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief returns heat generation rate !> @brief returns heat generation rate
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine thermal_adiabatic_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el) 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) :: & integer, intent(in) :: &
ip, & !< integration point number 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) dTdot_dT = dTdot_dT/real(homogenization_Ngrains(homog),pReal)
end subroutine thermal_adiabatic_getSourceAndItsTangent end subroutine thermal_adiabatic_getSourceAndItsTangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief returns homogenized specific heat capacity !> @brief returns homogenized specific heat capacity
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function thermal_adiabatic_getSpecificHeat(ip,el) 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) :: & integer, intent(in) :: &
ip, & !< integration point number ip, & !< integration point number
@ -270,13 +229,6 @@ end function thermal_adiabatic_getSpecificHeat
!> @brief returns homogenized mass density !> @brief returns homogenized mass density
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function thermal_adiabatic_getMassDensity(ip,el) 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) :: & integer, intent(in) :: &
ip, & !< integration point number ip, & !< integration point number
@ -304,8 +256,6 @@ end function thermal_adiabatic_getMassDensity
!> @brief return array of thermal results !> @brief return array of thermal results
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function thermal_adiabatic_postResults(homog,instance,of) result(postResults) function thermal_adiabatic_postResults(homog,instance,of) result(postResults)
use material, only: &
temperature
integer, intent(in) :: & integer, intent(in) :: &
homog, & homog, &