Merge branch 'simplify-cmake' into fix-thermalexpansion

This commit is contained in:
Martin Diehl 2021-05-23 21:59:15 +02:00
commit b0f4ea45a4
55 changed files with 705 additions and 671 deletions

View File

@ -1,6 +1,12 @@
cmake_minimum_required (VERSION 3.10.0)
include (FindPkgConfig REQUIRED)
if (DEFINED ENV{PETSC_DIR})
message ("PETSC_DIR:\n$ENV{PETSC_DIR}\n")
else ()
message (FATAL_ERROR "PETSc location (PETSC_DIR) is not defined")
endif ()
# Dummy project to determine compiler names and version
project (Prerequisites LANGUAGES)
set(ENV{PKG_CONFIG_PATH} "$ENV{PETSC_DIR}/$ENV{PETSC_ARCH}/lib/pkgconfig")
@ -8,81 +14,8 @@ pkg_check_modules (PETSC REQUIRED PETSc>=3.12.0 PETSc<3.16.0)
pkg_get_variable (CMAKE_Fortran_COMPILER PETSc fcompiler)
pkg_get_variable (CMAKE_C_COMPILER PETSc ccompiler)
find_program (CAT_EXECUTABLE NAMES cat)
execute_process (COMMAND ${CAT_EXECUTABLE} ${PROJECT_SOURCE_DIR}/VERSION
RESULT_VARIABLE DAMASK_VERSION_RETURN
OUTPUT_VARIABLE DAMASK_VERSION
OUTPUT_STRIP_TRAILING_WHITESPACE)
#---------------------------------------------------------------------------------------
# Find PETSc from system environment
set(PETSC_DIR $ENV{PETSC_DIR})
if (PETSC_DIR STREQUAL "")
message (FATAL_ERROR "PETSc location (PETSC_DIR) is not defined")
endif ()
set (petsc_conf_variables "${PETSC_DIR}/lib/petsc/conf/variables")
set (petsc_conf_rules "${PETSC_DIR}/lib/petsc/conf/rules" )
# Use existing variables from PETSc
# https://github.com/jedbrown/cmake-modules/blob/master/FindPETSc.cmake
# Generate a temporary makefile to probe the PETSc configuration
# This file will be deleted once the settings from PETSc are parsed into CMake
exec_program (mktemp ARGS -d OUTPUT_VARIABLE TEMPDIR)
set (petsc_config_makefile "${TEMPDIR}/Makefile.petsc")
file (WRITE
"${petsc_config_makefile}"
"## This file was auto generated by CMake
# PETSC_DIR = ${PETSC_DIR}
SHELL = /bin/sh
include ${petsc_conf_rules}
include ${petsc_conf_variables}
INCLUDE_DIRS := \${PETSC_FC_INCLUDES}
LIBRARIES := \${PETSC_WITH_EXTERNAL_LIB}
includes:
\t@echo \${INCLUDE_DIRS}
extlibs:
\t@echo \${LIBRARIES}
")
# CMake will execute each target in the ${petsc_config_makefile}
# to acquire corresponding PETSc Variables.
find_program (MAKE_EXECUTABLE NAMES gmake make)
# Find the PETSc includes directory settings
execute_process (COMMAND ${MAKE_EXECUTABLE} --no-print-directory -f ${petsc_config_makefile} "includes"
RESULT_VARIABLE PETSC_INCLUDES_RETURN
OUTPUT_VARIABLE petsc_includes
OUTPUT_STRIP_TRAILING_WHITESPACE)
# Find the PETSc external linking directory settings
execute_process (COMMAND ${MAKE_EXECUTABLE} --no-print-directory -f ${petsc_config_makefile} "extlibs"
RESULT_VARIABLE PETSC_EXTERNAL_LIB_RETURN
OUTPUT_VARIABLE petsc_external_lib
OUTPUT_STRIP_TRAILING_WHITESPACE)
# Remove temporary makefile, no need to keep it anymore.
file (REMOVE_RECURSE ${TEMPDIR})
# Remove duplicate compiler and linker flags
string (REGEX MATCHALL "-I([^\" ]+)" TMP_LIST "${petsc_includes}")
list (REMOVE_DUPLICATES TMP_LIST)
foreach (dir ${TMP_LIST})
set (PETSC_INCLUDES "${PETSC_INCLUDES} ${dir}")
endforeach (dir)
string (REGEX MATCHALL "-[lLW]([^\" ]+)" TMP_LIST "${petsc_external_lib}")
list (REMOVE_DUPLICATES TMP_LIST)
foreach (exlib ${TMP_LIST})
set (PETSC_EXTERNAL_LIB "${PETSC_EXTERNAL_LIB} ${exlib}")
endforeach (exlib)
message ("Found PETSC_DIR:\n${PETSC_DIR}\n" )
message ("Found PETSC_INCLUDES:\n${PETSC_INCLUDES}\n" )
message ("Found PETSC_EXTERNAL_LIB:\n${PETSC_EXTERNAL_LIB}\n")
#---------------------------------------------------------------------------------------
# Now start to care about DAMASK
# DAMASK solver defines project to build
string(TOLOWER ${DAMASK_SOLVER} DAMASK_SOLVER)
# Solver determines name of project
string(TOLOWER "${DAMASK_SOLVER}" DAMASK_SOLVER)
if (DAMASK_SOLVER STREQUAL "grid")
project (damask-grid HOMEPAGE_URL https://damask.mpie.de LANGUAGES Fortran C)
add_definitions (-DGrid)
@ -92,10 +25,13 @@ elseif (DAMASK_SOLVER STREQUAL "mesh")
else ()
message (FATAL_ERROR "Build target (DAMASK_SOLVER) is not defined")
endif ()
add_definitions (-DDAMASKVERSION="${DAMASK_VERSION}")
add_definitions (-DPETSc)
message ("\nBuilding ${CMAKE_PROJECT_NAME}\n")
file (STRINGS ${PROJECT_SOURCE_DIR}/VERSION DAMASK_VERSION)
message ("\nBuilding ${CMAKE_PROJECT_NAME} ${DAMASK_VERSION}\n")
add_definitions (-DPETSc)
add_definitions (-DDAMASKVERSION="${DAMASK_VERSION}")
if (CMAKE_BUILD_TYPE STREQUAL "")
set (CMAKE_BUILD_TYPE "RELEASE")
@ -139,12 +75,22 @@ if (CMAKE_Fortran_COMPILER_ID STREQUAL "Intel")
include (Compiler-Intel)
elseif(CMAKE_Fortran_COMPILER_ID STREQUAL "GNU")
include (Compiler-GNU)
elseif(CMAKE_Fortran_COMPILER_ID STREQUAL "PGI")
include (Compiler-PGI)
else ()
message (FATAL_ERROR "Compiler type (CMAKE_Fortran_COMPILER_ID) not recognized")
endif ()
file (STRINGS "$ENV{PETSC_DIR}/$ENV{PETSC_ARCH}/lib/petsc/conf/petscvariables" PETSC_EXTERNAL_LIB REGEX "PETSC_WITH_EXTERNAL_LIB = .*$?")
string (REGEX MATCHALL "-[lLW]([^\" ]+)" PETSC_EXTERNAL_LIB "${PETSC_EXTERNAL_LIB}")
list (REMOVE_DUPLICATES PETSC_EXTERNAL_LIB)
string (REPLACE ";" " " PETSC_EXTERNAL_LIB "${PETSC_EXTERNAL_LIB}")
message ("PETSC_EXTERNAL_LIB:\n${PETSC_EXTERNAL_LIB}\n")
file (STRINGS "$ENV{PETSC_DIR}/$ENV{PETSC_ARCH}/lib/petsc/conf/petscvariables" PETSC_INCLUDES REGEX "PETSC_FC_INCLUDES = .*$?")
string (REGEX MATCHALL "-I([^\" ]+)" PETSC_INCLUDES "${PETSC_INCLUDES}")
list (REMOVE_DUPLICATES PETSC_INCLUDES)
string (REPLACE ";" " " PETSC_INCLUDES "${PETSC_INCLUDES}")
message ("PETSC_INCLUDES:\n${PETSC_INCLUDES}\n")
set (CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE} "${BUILDCMD_PRE} ${OPENMP_FLAGS} ${STANDARD_CHECK} ${OPTIMIZATION_FLAGS} ${COMPILE_FLAGS} ${PRECISION_FLAGS}")
set (CMAKE_Fortran_LINK_EXECUTABLE "${BUILDCMD_PRE} ${CMAKE_Fortran_COMPILER} ${OPENMP_FLAGS} ${OPTIMIZATION_FLAGS} ${LINKER_FLAGS}")

@ -1 +1 @@
Subproject commit 95f7faea920dd6956884e4a55f72e5d5b1ffcdc8
Subproject commit 0744cf7f93dbd06423baaae97a67959f11c3e6a5

View File

@ -1 +1 @@
v3.0.0-alpha3-51-g38ca85b58
v3.0.0-alpha3-137-g0b80252d9

View File

@ -1,52 +0,0 @@
###################################################################################################
# PGI Compiler
###################################################################################################
if (OPENMP)
set (OPENMP_FLAGS "-mp")
else ()
set (OPENMP_FLAGS "-nomp")
endif ()
if (OPTIMIZATION STREQUAL "OFF")
set (OPTIMIZATION_FLAGS "-O0" )
elseif (OPTIMIZATION STREQUAL "DEFENSIVE")
set (OPTIMIZATION_FLAGS "-O2 -fast")
elseif (OPTIMIZATION STREQUAL "AGGRESSIVE")
set (OPTIMIZATION_FLAGS "-O4 -fast -Mvect=sse")
endif ()
set (STANDARD_CHECK "-Mallocatable=03 -Mstandard")
#------------------------------------------------------------------------------------------------
# Fine tuning compilation options
set (COMPILE_FLAGS "${COMPILE_FLAGS} -Mpreprocess")
# preprocessor
set (COMPILE_FLAGS "${COMPILE_FLAGS} -Minfo=all")
# instructs the compiler to produce information on standard error
set (COMPILE_FLAGS "${COMPILE_FLAGS} -Minform=warn")
# instructs the compiler to display error messages at the specified and higher levels
set (COMPILE_FLAGS "${COMPILE_FLAGS} -Mdclchk")
# instructs the compiler to require that all program variables be declared
#------------------------------------------------------------------------------------------------O
# Runtime debugging
set (DEBUG_FLAGS "${DEBUG_FLAGS} -g")
# Includes debugging information in the object module; sets the optimization level to zero unless a -O option is present on the command line
set (DEBUG_FLAGS "${DEBUG_FLAGS} -C")
# Generates code to check array bounds
set (DEBUG_FLAGS "${DEBUG_FLAGS} -Mchkptr")
# Check for NULL pointers (pgf95, pgfortran only)
set (DEBUG_FLAGS "${DEBUG_FLAGS} -Mchkstk")
# Check the stack for available space upon entry to and before the start of a parallel region. Useful when many private variables are declared
set (DEBUG_FLAGS "${DEBUG_FLAGS} -Mbounds")
# Specifies whether array bounds checking is enabled or disabled
#------------------------------------------------------------------------------------------------
# precision settings
set (PRECISION_FLAGS "${PRECISION_FLAGS} -r8")
# Determines whether the compiler promotes REAL variables and constants to DOUBLE PRECISION

View File

@ -17,7 +17,7 @@ phase:
atol_xi: 1.0
dot_gamma_0_sl: 0.001
h_0_sl_sl: 75e6
h_sl_sl: [1, 1, 1.4, 1.4, 1.4, 1.4]
h_sl_sl: [1, 1, 1.4, 1.4, 1.4, 1.4, 1.4]
n_sl: 20
output: [xi_sl]
xi_0_sl: [31e6]

View File

@ -51,5 +51,5 @@ D0 4.0e-5 # Vacancy diffusion prefactor [m**2/
Qsd 4.5e-19 # Activation energy for climb [J]
Catomicvolume 1.0 # Adj. parameter controlling the atomic volume [in b]
Cedgedipmindistance 1.0 # Adj. parameter controlling the minimum dipole distance [in b]
interaction_slipslip 0.009 0.009 0.72 0.05 0.09 0.06
interaction_slipslip 0.009 0.72 0.009 0.05 0.05 0.06 0.09
nonschmid_coefficients 0.938 0.71 4.43 0.0 0.0 0.0

View File

@ -19,7 +19,7 @@ TWIP_Steel_FeMnC:
D_0: 4.0e-5 # Vacancy diffusion prefactor / m^2/s
D_a: 1.0 # minimum dipole distance / b
Q_cl: 4.5e-19 # Activation energy for climb / J
h_sl_sl: [0.122, 0.122, 0.625, 0.07, 0.137, 0.122] # Interaction coefficients (Kubin et al. 2008)
h_sl_sl: [0.122, 0.122, 0.625, 0.07, 0.137, 0.137, 0.122] # Interaction coefficients (Kubin et al. 2008)
# shear band parameters
xi_sb: 180.0e6
Q_sb: 3.7e-19

View File

@ -18,4 +18,4 @@ Tungsten:
D_0: 4.0e-5 # Vacancy diffusion prefactor / m^2/s
D_a: 1.0 # minimum dipole distance / b
Q_cl: 4.5e-19 # Activation energy for climb / J
h_sl_sl: [1, 1, 1.4, 1.4, 1.4, 1.4]
h_sl_sl: [1, 1.4, 1, 1.4, 1.4, 1.4, 1.4]

View File

@ -52,7 +52,7 @@ q 1 # exponent for thermal barrier
attackFrequency 50e9 # attack frequency in Hz
surfaceTransmissivity 1.0 # transmissivity of free surfaces for dislocation flux
grainboundaryTransmissivity 0.0 # transmissivity of grain boundaries for dislocation flux (grain bundaries are identified as interfaces with different textures on both sides); if not set or set to negative number, the subroutine automatically determines the transmissivity at the grain boundary
interaction_SlipSlip 0 0 0.625 0.07 0.137 0.122 # Dislocation interaction coefficient
interaction_SlipSlip 0 0 0.625 0.07 0.137 0.137 0.122 # Dislocation interaction coefficient
linetension 0.8 # constant indicating the effect of the line tension on the hardening coefficients (0 to 1)
edgejog 1.0 # fraction of annihilated screw dipoles that forms edge jogs (0 to 1)
shortRangeStressCorrection 0 # switch for use of short range correction stress

View File

@ -57,6 +57,6 @@ significantN 1
shortRangeStressCorrection 0
CFLfactor 1.1 # safety factor for CFL flux check (numerical parameter)
r 1
interaction_SlipSlip 0 0 0.625 0.07 0.137 0.122 # Dislocation interaction coefficient
interaction_SlipSlip 0 0 0.625 0.07 0.137 0.137 0.122 # Dislocation interaction coefficient
linetension 0.8
edgejog 0.01 # 0.2

View File

@ -8,7 +8,7 @@ Aluminum:
a_sl: 2.25
dot_gamma_0_sl: 0.001
h_0_sl_sl: 75e6
h_sl_sl: [1, 1, 1.4, 1.4, 1.4, 1.4]
h_sl_sl: [1, 1, 1.4, 1.4, 1.4, 1.4, 1.4]
n_sl: 20
output: [xi_sl, gamma_sl]
type: phenopowerlaw

View File

@ -10,7 +10,7 @@ Ferrite:
a_sl: 2.0
dot_gamma_0_sl: 0.001
h_0_sl_sl: 1000.0e6
h_sl_sl: [1, 1, 1.4, 1.4, 1.4, 1.4]
h_sl_sl: [1, 1.4, 1, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4]
n_sl: 20
type: phenopowerlaw
xi_0_sl: [95.e6, 96.e6]

View File

@ -10,7 +10,7 @@ Martensite:
a_sl: 2.0
dot_gamma_0_sl: 0.001
h_0_sl_sl: 563.0e9
h_sl_sl: [1, 1, 1.4, 1.4, 1.4, 1.4]
h_sl_sl: [1, 1.4, 1, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4]
n_sl: 20
type: phenopowerlaw
xi_0_sl: [405.8e6, 456.7e6]

View File

@ -1,8 +1,11 @@
N_cl: [3]
dot_o: 1e-3
g_crit: [0.50e7]
q: 20
s_crit: [0.006666]
type: anisobrittle
N_cl: [3]
g_crit: [0.50e7]
s_crit: [0.006666]
dot_o: 1e-3
q: 20
output: [f_phi]
D_11: 1.0
M: 0.001

View File

@ -0,0 +1,9 @@
type: isobrittle
W_crit: 1400000.0
m: 1.0
isoBrittle_atol: 0.01
output: [f_phi]
D_11: 1.0
M: 0.001

View File

@ -18,7 +18,7 @@ q_sl: [1.55, 1.55]
i_sl: [23.3, 23.3]
D_a: 7.4 # C_anni
B: [0.001, 0.001]
h_sl_sl: [0.1, 0.1, 0.72, 0.053, 0.137, 0.073]
h_sl_sl: [0.1, 0.72, 0.1, 0.053, 0.053, 0.073, 0.137, 0.72, 0.72, 0.053, 0.053, 0.053, 0.053, 0.073, 0.073, 0.073, 0.073, 0.073, 0.073, 0.137, 0.073, 0.073, 0.137, 0.073]
D_0: 4.000E-05
Q_cl: 5.400E-19 # no recovery!
D: 40e-6 # estimated

View File

@ -8,7 +8,7 @@ N_sl: [12]
n_sl: 83
dot_gamma_0_sl: 0.001
h_0_sl_sl: 75e6
h_sl_sl: [1, 1, 1.4, 1.4, 1.4, 1.4]
h_sl_sl: [1, 1, 1.4, 1.4, 1.4, 1.4, 1.4]
a_sl: 1.0
xi_0_sl: [26e6]
xi_inf_sl: [53e6]

View File

@ -17,7 +17,7 @@ phase:
atol_xi: 1.0
dot_gamma_0_sl: 0.001
h_0_sl_sl: 75e6
h_sl_sl: [1, 1, 1.4, 1.4, 1.4, 1.4]
h_sl_sl: [1, 1, 1.4, 1.4, 1.4, 1.4, 1.4]
n_sl: 20
output: [xi_sl]
xi_0_sl: [31e6]

View File

@ -1,22 +1,22 @@
homogenization:
SX:
N_constituents: 1
mechanics: {type: pass}
mechanical: {type: pass}
phase:
Aluminum:
lattice: cF
mechanics:
mechanical:
output: [F, P, F_e, F_p, L_p]
elasticity: {type: Hooke, C_11: 106.75e9, C_12: 60.41e9, C_44: 28.34e9}
plasticity:
elastic: {type: Hooke, C_11: 106.75e9, C_12: 60.41e9, C_44: 28.34e9}
plastic:
type: phenopowerlaw
N_sl: [12]
a_sl: 2.25
atol_xi: 1.0
dot_gamma_0_sl: 0.001
h_0_sl_sl: 75e6
h_sl_sl: [1, 1, 1.4, 1.4, 1.4, 1.4]
h_sl_sl: [1, 1, 1.4, 1.4, 1.4, 1.4, 1.4]
n_sl: 20
output: [xi_sl]
xi_0_sl: [31e6]

View File

@ -503,7 +503,7 @@ then
PROFILE=" $PROFILE -pg"
fi
FORT_OPT="-c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr -mp1 -WB -fp-model source"
FORT_OPT="-c -implicitnone -stand f18 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr -mp1 -WB -fp-model source"
if test "$MTHREAD" = "OPENMP"
then
FORT_OPT=" $FORT_OPT -qopenmp"
@ -541,15 +541,15 @@ fi
# DAMASK compiler calls: additional flags are in line 2 OpenMP flags in line 3
DFORTLOWMP="$FCOMP -c -O0 -qno-offload -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB $I8FFLAGS -I$MARC_SOURCE/common \
DFORTLOWMP="$FCOMP -c -O0 -qno-offload -implicitnone -stand f18 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2020 -DDAMASKVERSION=$DAMASKVERSION \
-qopenmp -qopenmp-threadprivate=compat\
$MUMPS_INCLUDE $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS -I$KDTREE2_MOD"
DFORTRANMP="$FCOMP -c -O1 -qno-offload -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB $I8FFLAGS -I$MARC_SOURCE/common \
DFORTRANMP="$FCOMP -c -O1 -qno-offload -implicitnone -stand f18 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2020 -DDAMASKVERSION=$DAMASKVERSION \
-qopenmp -qopenmp-threadprivate=compat\
$MUMPS_INCLUDE $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS -I$KDTREE2_MOD"
DFORTHIGHMP="$FCOMP -c -O3 -qno-offload -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB $I8FFLAGS -I$MARC_SOURCE/common \
DFORTHIGHMP="$FCOMP -c -O3 -qno-offload -implicitnone -stand f18 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2020 -DDAMASKVERSION=$DAMASKVERSION \
-qopenmp -qopenmp-threadprivate=compat\
$MUMPS_INCLUDE $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS -I$KDTREE2_MOD"

View File

@ -14,7 +14,6 @@ from . import tensor # noqa
from . import mechanics # noqa
from . import solver # noqa
from . import grid_filters # noqa
from . import lattice # noqa
#Modules that contain only one class (of the same name), are prefixed by a '_'.
#For example, '_colormap' containsa class called 'Colormap' which is imported as 'damask.Colormap'.
from ._rotation import Rotation # noqa

View File

@ -110,13 +110,13 @@ class Colormap(mpl.colors.ListedColormap):
low_,high_ = map(Colormap._rgb2msh,low_high)
elif model.lower() == 'hsv':
if np.any(low_high<0) or np.any(low_high[:,1:3]>1) or np.any(low_high[:,0]>360):
if np.any(low_high<0) or np.any(low_high>[360,1,1]):
raise ValueError(f'HSV color {low} | {high} are out of range.')
low_,high_ = map(Colormap._hsv2msh,low_high)
elif model.lower() == 'hsl':
if np.any(low_high<0) or np.any(low_high[:,1:3]>1) or np.any(low_high[:,0]>360):
if np.any(low_high<0) or np.any(low_high>[360,1,1]):
raise ValueError(f'HSL color {low} | {high} are out of range.')
low_,high_ = map(Colormap._hsl2msh,low_high)

View File

@ -63,6 +63,30 @@ kinematics = {
[+1,-1,+1 , -1,+1,+2],
[-1,+1,+1 , +1,-1,+2],
[+1,+1,+1 , +1,+1,-2],
[+1,+1,-1 , +1,+2,+3],
[+1,-1,+1 , -1,+2,+3],
[-1,+1,+1 , +1,-2,+3],
[+1,+1,+1 , +1,+2,-3],
[+1,-1,+1 , +1,+3,+2],
[+1,+1,-1 , -1,+3,+2],
[+1,+1,+1 , +1,-3,+2],
[-1,+1,+1 , +1,+3,-2],
[+1,+1,-1 , +2,+1,+3],
[+1,-1,+1 , -2,+1,+3],
[-1,+1,+1 , +2,-1,+3],
[+1,+1,+1 , +2,+1,-3],
[+1,-1,+1 , +2,+3,+1],
[+1,+1,-1 , -2,+3,+1],
[+1,+1,+1 , +2,-3,+1],
[-1,+1,+1 , +2,+3,-1],
[-1,+1,+1 , +3,+1,+2],
[+1,+1,+1 , -3,+1,+2],
[+1,+1,-1 , +3,-1,+2],
[+1,-1,+1 , +3,+1,-2],
[-1,+1,+1 , +3,+2,+1],
[+1,+1,+1 , -3,+2,+1],
[+1,+1,-1 , +3,-2,+1],
[+1,-1,+1 , +3,+2,-1],
],'d'),
'twin' : _np.array([
[-1, 1, 1, 2, 1, 1],

View File

@ -5,6 +5,35 @@ import numpy as np
from . import Rotation
from . import util
from . import tensor
from . import _lattice
_crystal_families = ['triclinic',
'monoclinic',
'orthorhombic',
'tetragonal',
'hexagonal',
'cubic']
_lattice_symmetries = {
'aP': 'triclinic',
'mP': 'monoclinic',
'mS': 'monoclinic',
'oP': 'orthorhombic',
'oS': 'orthorhombic',
'oI': 'orthorhombic',
'oF': 'orthorhombic',
'tP': 'tetragonal',
'tI': 'tetragonal',
'hP': 'hexagonal',
'cP': 'cubic',
'cI': 'cubic',
'cF': 'cubic',
}
_parameter_doc = \
"""lattice : str
@ -33,7 +62,7 @@ class Orientation(Rotation):
"""
Representation of crystallographic orientation as combination of rotation and either crystal family or Bravais lattice.
The crystal family is one of Orientation.crystal_families:
The crystal family is one of:
- triclinic
- monoclinic
@ -45,7 +74,7 @@ class Orientation(Rotation):
and enables symmetry-related operations such as
"equivalent", "reduced", "disorientation", "IPF_color", or "to_SST".
The Bravais lattice is one of Orientation.lattice_symmetries:
The Bravais lattice is given in the Pearson notation:
- triclinic
- aP : primitive
@ -85,35 +114,6 @@ class Orientation(Rotation):
"""
crystal_families = ['triclinic',
'monoclinic',
'orthorhombic',
'tetragonal',
'hexagonal',
'cubic']
lattice_symmetries = {
'aP': 'triclinic',
'mP': 'monoclinic',
'mS': 'monoclinic',
'oP': 'orthorhombic',
'oS': 'orthorhombic',
'oI': 'orthorhombic',
'oF': 'orthorhombic',
'tP': 'tetragonal',
'tI': 'tetragonal',
'hP': 'hexagonal',
'cP': 'cubic',
'cI': 'cubic',
'cF': 'cubic',
}
@util.extend_docstring(_parameter_doc)
def __init__(self,
rotation = None,
@ -132,34 +132,17 @@ class Orientation(Rotation):
Defaults to no rotation.
"""
from damask.lattice import kinematics
Rotation.__init__(self) if rotation is None else Rotation.__init__(self,rotation=rotation)
if ( lattice not in self.lattice_symmetries
and lattice not in self.crystal_families):
raise KeyError(f'Lattice "{lattice}" is unknown')
self.family = None
self.lattice = None
self.a = None
self.b = None
self.c = None
self.alpha = None
self.beta = None
self.gamma = None
self.kinematics = None
if lattice in self.lattice_symmetries:
self.family = self.lattice_symmetries[lattice]
if lattice in _lattice_symmetries:
self.family = _lattice_symmetries[lattice]
self.lattice = lattice
self.a = 1 if a is None else a
self.b = b
self.c = c
self.alpha = (np.radians(alpha) if degrees else alpha) if alpha is not None else None
self.beta = (np.radians(beta) if degrees else beta) if beta is not None else None
self.gamma = (np.radians(gamma) if degrees else gamma) if gamma is not None else None
self.a = float(self.a) if self.a is not None else \
(self.b / self.ratio['b'] if self.b is not None and self.ratio['b'] is not None else
self.c / self.ratio['c'] if self.c is not None and self.ratio['c'] is not None else None)
@ -171,9 +154,13 @@ class Orientation(Rotation):
(self.a * self.ratio['c'] if self.a is not None and self.ratio['c'] is not None else
self.b / self.ratio['b'] * self.ratio['c']
if self.c is not None and self.ratio['b'] is not None and self.ratio['c'] is not None else None)
self.alpha = self.alpha if self.alpha is not None else self.immutable['alpha'] if 'alpha' in self.immutable else None
self.beta = self.beta if self.beta is not None else self.immutable['beta'] if 'beta' in self.immutable else None
self.gamma = self.gamma if self.gamma is not None else self.immutable['gamma'] if 'gamma' in self.immutable else None
self.alpha = np.radians(alpha) if degrees and alpha is not None else alpha
self.beta = np.radians(beta) if degrees and beta is not None else beta
self.gamma = np.radians(gamma) if degrees and gamma is not None else gamma
if self.alpha is None and 'alpha' in self.immutable: self.alpha = self.immutable['alpha']
if self.beta is None and 'beta' in self.immutable: self.beta = self.immutable['beta']
if self.gamma is None and 'gamma' in self.immutable: self.gamma = self.immutable['gamma']
if \
(self.a is None) \
@ -190,16 +177,22 @@ class Orientation(Rotation):
> np.sum(np.roll([self.alpha,self.beta,self.gamma],r)[1:]) for r in range(3)]):
raise ValueError ('Each lattice angle must be less than sum of others')
if self.lattice in kinematics:
master = kinematics[self.lattice]
if self.lattice in _lattice.kinematics:
master = _lattice.kinematics[self.lattice]
self.kinematics = {}
for m in master:
self.kinematics[m] = {'direction':master[m][:,0:3],'plane':master[m][:,3:6]} \
if master[m].shape[-1] == 6 else \
{'direction':self.Bravais_to_Miller(uvtw=master[m][:,0:4]),
'plane': self.Bravais_to_Miller(hkil=master[m][:,4:8])}
elif lattice in self.crystal_families:
self.family = lattice
elif lattice in _crystal_families:
self.family = lattice
self.lattice = None
self.a = self.b = self.c = None
self.alpha = self.beta = self.gamma = None
else:
raise KeyError(f'Lattice "{lattice}" is unknown')
def __repr__(self):
@ -676,11 +669,9 @@ class Orientation(Rotation):
https://doi.org/10.1016/j.actamat.2004.11.021
"""
from damask.lattice import relations
if model not in relations:
if model not in _lattice.relations:
raise KeyError(f'Orientation relationship "{model}" is unknown')
r = relations[model]
r = _lattice.relations[model]
if self.lattice not in r:
raise KeyError(f'Relationship "{model}" not supported for lattice "{self.lattice}"')

View File

@ -629,7 +629,7 @@ class Rotation:
else:
if np.any(qu[...,0] < 0.0):
raise ValueError('Quaternion with negative first (real) component.')
if not np.all(np.isclose(np.linalg.norm(qu,axis=-1), 1.0)):
if not np.all(np.isclose(np.linalg.norm(qu,axis=-1), 1.0,rtol=0.0)):
raise ValueError('Quaternion is not of unit length.')
return Rotation(qu)

View File

@ -384,7 +384,7 @@ class Table:
Returns
-------
udated : damask.Table
updated : damask.Table
Updated table.
"""
@ -412,7 +412,7 @@ class Table:
Returns
-------
udated : damask.Table
updated : damask.Table
Updated table.
"""
@ -435,7 +435,7 @@ class Table:
Returns
-------
udated : damask.Table
updated : damask.Table
Updated table.
"""
@ -461,7 +461,7 @@ class Table:
Returns
-------
udated : damask.Table
updated : damask.Table
Updated table.
"""
@ -494,7 +494,7 @@ class Table:
Returns
-------
udated : damask.Table
updated : damask.Table
Updated table.
"""
@ -519,7 +519,7 @@ class Table:
Returns
-------
udated : damask.Table
updated : damask.Table
Updated table.
"""

View File

@ -150,9 +150,8 @@ class VTK:
----------
fname : str or pathlib.Path
Filename for reading. Valid extensions are .vtr, .vtu, .vtp, and .vtk.
dataset_type : str, optional
dataset_type : {'vtkRectilinearGrid', 'vtkUnstructuredGrid', 'vtkPolyData'}, optional
Name of the vtk.vtkDataSet subclass when opening a .vtk file.
Valid types are vtkRectilinearGrid, vtkUnstructuredGrid, and vtkPolyData.
Returns
-------

View File

@ -124,9 +124,6 @@ def strain(F,t,m):
"""
Calculate strain tensor (SethHill family).
For details refer to https://en.wikipedia.org/wiki/Finite_strain_theory and
https://de.wikipedia.org/wiki/Verzerrungstensor
Parameters
----------
F : numpy.ndarray of shape (...,3,3)
@ -142,6 +139,11 @@ def strain(F,t,m):
epsilon : numpy.ndarray of shape (...,3,3)
Strain of F.
References
----------
https://en.wikipedia.org/wiki/Finite_strain_theory
https://de.wikipedia.org/wiki/Verzerrungstensor
"""
if t == 'V':
w,n = _np.linalg.eigh(deformation_Cauchy_Green_left(F))
@ -150,7 +152,6 @@ def strain(F,t,m):
if m > 0.0:
eps = 1.0/(2.0*abs(m)) * (+ _np.einsum('...j,...kj,...lj',w**m,n,n) - _np.eye(3))
elif m < 0.0:
eps = 1.0/(2.0*abs(m)) * (- _np.einsum('...j,...kj,...lj',w**m,n,n) + _np.eye(3))
else:

View File

@ -3,8 +3,8 @@
from scipy import spatial as _spatial
import numpy as _np
from . import util
from . import grid_filters
from . import util as _util
from . import grid_filters as _grid_filters
def from_random(size,N_seeds,cells=None,rng_seed=None):
@ -34,7 +34,7 @@ def from_random(size,N_seeds,cells=None,rng_seed=None):
if cells is None:
coords = rng.random((N_seeds,3)) * size
else:
grid_coords = grid_filters.coordinates0_point(cells,size).reshape(-1,3,order='F')
grid_coords = _grid_filters.coordinates0_point(cells,size).reshape(-1,3,order='F')
coords = grid_coords[rng.choice(_np.prod(cells),N_seeds, replace=False)] \
+ _np.broadcast_to(size/cells,(N_seeds,3))*(rng.random((N_seeds,3))*.5-.25) # wobble without leaving cells
@ -73,7 +73,7 @@ def from_Poisson_disc(size,N_seeds,N_candidates,distance,periodic=True,rng_seed=
s = 1
i = 0
progress = util._ProgressBar(N_seeds+1,'',50)
progress = _util._ProgressBar(N_seeds+1,'',50)
while s < N_seeds:
candidates = rng.random((N_candidates,3))*_np.broadcast_to(size,(N_candidates,3))
tree = _spatial.cKDTree(coords[:s],boxsize=size) if periodic else \
@ -120,7 +120,7 @@ def from_grid(grid,selection=None,invert=False,average=False,periodic=True):
material = grid.material.reshape((-1,1),order='F')
mask = _np.full(grid.cells.prod(),True,dtype=bool) if selection is None else \
_np.isin(material,selection,invert=invert).flatten()
coords = grid_filters.coordinates0_point(grid.cells,grid.size).reshape(-1,3,order='F')
coords = _grid_filters.coordinates0_point(grid.cells,grid.size).reshape(-1,3,order='F')
if not average:
return (coords[mask],material[mask])

View File

@ -1,3 +1,5 @@
"""Miscellaneous helper functionality."""
import sys
import datetime
import os
@ -177,26 +179,36 @@ def execute(cmd,wd='./',env=None):
def natural_sort(key):
"""
Natural sort.
For use in python's 'sorted'.
References
----------
https://en.wikipedia.org/wiki/Natural_sort_order
"""
convert = lambda text: int(text) if text.isdigit() else text
return [ convert(c) for c in re.split('([0-9]+)', key) ]
def show_progress(iterable,N_iter=None,prefix='',bar_length=50):
"""
Decorate a loop with a status bar.
Decorate a loop with a progress bar.
Use similar like enumerate.
Parameters
----------
iterable : iterable/function with yield statement
iterable : iterable or function with yield statement
Iterable (or function with yield statement) to be decorated.
N_iter : int
Total # of iterations. Needed if number of iterations can not be obtained as len(iterable).
prefix : str, optional.
N_iter : int, optional
Total number of iterations. Required unless obtainable as len(iterable).
prefix : str, optional
Prefix string.
bar_length : int, optional
Character length of bar. Defaults to 50.
Length of progress bar in characters. Defaults to 50.
"""
if N_iter in [0,1] or (hasattr(iterable,'__len__') and len(iterable) <= 1):
@ -509,6 +521,7 @@ def dict_prune(d):
v = dict_prune(v)
if not isinstance(v,dict) or v != {}:
new[k] = v
return new

View File

@ -12,7 +12,7 @@ setuptools.setup(
author='The DAMASK team',
author_email='damask@mpie.de',
description='DAMASK library',
long_description='Python library for pre and post processing of DAMASK simulations',
long_description='Python library for managing DAMASK simulations',
url='https://damask.mpie.de',
packages=setuptools.find_packages(),
include_package_data=True,

View File

@ -23,3 +23,27 @@
-0.2357022603955159 0.2357022603955159 0.4714045207910318 0.23570226039551587 -0.23570226039551587 -0.47140452079103173 -0.2357022603955159 0.2357022603955159 0.4714045207910318
-0.2357022603955158 0.2357022603955158 -0.4714045207910316 0.2357022603955159 -0.2357022603955159 0.4714045207910318 0.2357022603955159 -0.2357022603955159 0.4714045207910318
0.2357022603955159 0.23570226039551587 -0.4714045207910318 0.2357022603955159 0.23570226039551587 -0.4714045207910318 0.2357022603955159 0.23570226039551587 -0.4714045207910318
0.1543033499620919 0.3086066999241838 0.4629100498862757 0.1543033499620919 0.3086066999241838 0.4629100498862757 -0.15430334996209194 -0.3086066999241839 -0.4629100498862758
-0.15430334996209194 0.3086066999241839 0.4629100498862758 0.1543033499620919 -0.3086066999241838 -0.4629100498862757 -0.15430334996209194 0.3086066999241839 0.4629100498862758
-0.15430334996209188 0.30860669992418377 -0.46291004988627565 0.15430334996209194 -0.3086066999241839 0.4629100498862758 0.15430334996209194 -0.3086066999241839 0.4629100498862758
0.15430334996209194 0.3086066999241839 -0.4629100498862758 0.15430334996209194 0.3086066999241839 -0.4629100498862758 0.15430334996209194 0.3086066999241839 -0.4629100498862758
0.15430334996209194 0.4629100498862758 0.3086066999241838 -0.1543033499620919 -0.4629100498862757 -0.30860669992418377 0.15430334996209194 0.4629100498862758 0.3086066999241838
-0.1543033499620919 0.4629100498862757 0.30860669992418377 -0.1543033499620919 0.4629100498862757 0.30860669992418377 0.15430334996209194 -0.4629100498862758 -0.3086066999241838
0.15430334996209194 -0.4629100498862758 0.3086066999241839 0.15430334996209194 -0.4629100498862758 0.3086066999241839 0.15430334996209194 -0.4629100498862758 0.3086066999241839
-0.15430334996209188 -0.46291004988627565 0.30860669992418377 0.15430334996209194 0.4629100498862758 -0.3086066999241839 0.15430334996209194 0.4629100498862758 -0.3086066999241839
0.3086066999241838 0.15430334996209188 0.4629100498862757 0.3086066999241838 0.15430334996209188 0.4629100498862757 -0.3086066999241839 -0.1543033499620919 -0.4629100498862758
-0.3086066999241839 0.15430334996209197 0.4629100498862758 0.3086066999241838 -0.15430334996209194 -0.4629100498862757 -0.3086066999241839 0.15430334996209197 0.4629100498862758
-0.30860669992418377 0.1543033499620919 -0.46291004988627565 0.3086066999241839 -0.15430334996209197 0.4629100498862758 0.3086066999241839 -0.15430334996209197 0.4629100498862758
0.3086066999241839 0.1543033499620919 -0.4629100498862758 0.3086066999241839 0.1543033499620919 -0.4629100498862758 0.3086066999241839 0.1543033499620919 -0.4629100498862758
0.3086066999241839 0.4629100498862758 0.15430334996209188 -0.3086066999241838 -0.4629100498862757 -0.15430334996209186 0.3086066999241839 0.4629100498862758 0.15430334996209188
-0.3086066999241838 0.4629100498862757 0.15430334996209188 -0.3086066999241838 0.4629100498862757 0.15430334996209188 0.3086066999241839 -0.4629100498862758 -0.1543033499620919
0.3086066999241839 -0.4629100498862758 0.15430334996209194 0.3086066999241839 -0.4629100498862758 0.15430334996209194 0.3086066999241839 -0.4629100498862758 0.15430334996209194
-0.30860669992418377 -0.46291004988627565 0.15430334996209194 0.3086066999241839 0.4629100498862758 -0.154303349962092 0.3086066999241839 0.4629100498862758 -0.154303349962092
-0.46291004988627565 -0.15430334996209186 -0.3086066999241837 0.4629100498862758 0.1543033499620919 0.3086066999241838 0.4629100498862758 0.1543033499620919 0.3086066999241838
-0.4629100498862758 0.15430334996209197 0.3086066999241839 -0.4629100498862758 0.15430334996209197 0.3086066999241839 -0.4629100498862758 0.15430334996209197 0.3086066999241839
0.4629100498862757 -0.15430334996209194 0.30860669992418377 0.4629100498862757 -0.15430334996209194 0.30860669992418377 -0.4629100498862758 0.15430334996209197 -0.3086066999241838
0.4629100498862758 0.1543033499620919 -0.3086066999241839 -0.4629100498862757 -0.15430334996209188 0.3086066999241838 0.4629100498862758 0.1543033499620919 -0.3086066999241839
-0.46291004988627565 -0.3086066999241837 -0.1543033499620918 0.4629100498862758 0.3086066999241838 0.15430334996209188 0.4629100498862758 0.3086066999241838 0.15430334996209188
-0.4629100498862758 0.3086066999241839 0.15430334996209194 -0.4629100498862758 0.3086066999241839 0.15430334996209194 -0.4629100498862758 0.3086066999241839 0.15430334996209194
0.4629100498862757 -0.3086066999241838 0.1543033499620919 0.4629100498862757 -0.3086066999241838 0.1543033499620919 -0.4629100498862758 0.3086066999241839 -0.15430334996209194
0.4629100498862758 0.3086066999241838 -0.154303349962092 -0.4629100498862757 -0.30860669992418377 0.15430334996209197 0.4629100498862758 0.3086066999241838 -0.154303349962092

View File

@ -16,7 +16,7 @@ phase:
atol_xi: 1.0
dot_gamma_0_sl: 0.001
h_0_sl_sl: 75e6
h_sl_sl: [1, 1, 1.4, 1.4, 1.4, 1.4]
h_sl_sl: [1, 1, 1.4, 1.4, 1.4, 1.4, 1.4]
n_sl: 20
output: [xi_sl]
type: phenopowerlaw
@ -33,7 +33,7 @@ phase:
atol_xi: 1.0
dot_gamma_0_sl: 0.001
h_0_sl_sl: 75e6
h_sl_sl: [1, 1, 1.4, 1.4, 1.4, 1.4]
h_sl_sl: [1, 1.4, 1, 1.4, 1.4, 1.4, 1.4]
n_sl: 20
output: [xi_sl]
type: phenopowerlaw

View File

@ -665,7 +665,7 @@ phase:
atol_xi: 1.0
dot_gamma_0_sl: 0.001
h_0_sl_sl: 75e6
h_sl_sl: [1, 1, 1.4, 1.4, 1.4, 1.4]
h_sl_sl: [1, 1.4, 1, 1.4, 1.4, 1.4, 1.4]
n_sl: 20
output: [xi_sl]
type: phenopowerlaw

View File

@ -16,7 +16,7 @@ phase:
atol_xi: 1.0
dot_gamma_0_sl: 0.001
h_0_sl_sl: 75e6
h_sl_sl: [1, 1, 1.4, 1.4, 1.4, 1.4]
h_sl_sl: [1, 1, 1.4, 1.4, 1.4, 1.4, 1.4]
n_sl: 20
output: [xi_sl]
type: phenopowerlaw

View File

@ -5,9 +5,10 @@ from itertools import permutations
from damask import Rotation
from damask import Orientation
from damask import Table
from damask import lattice
from damask import util
from damask import grid_filters
from damask import _lattice as lattice
from damask._orientation import _crystal_families as crystal_families
@pytest.fixture
@ -22,7 +23,7 @@ def set_of_rodrigues(set_of_quaternions):
class TestOrientation:
@pytest.mark.parametrize('lattice',Orientation.crystal_families)
@pytest.mark.parametrize('lattice',crystal_families)
@pytest.mark.parametrize('shape',[None,5,(4,6)])
def test_equal(self,lattice,shape):
R = Rotation.from_random(shape)
@ -30,14 +31,14 @@ class TestOrientation:
(Orientation(R,lattice) == Orientation(R,lattice)).all()
@pytest.mark.parametrize('lattice',Orientation.crystal_families)
@pytest.mark.parametrize('lattice',crystal_families)
@pytest.mark.parametrize('shape',[None,5,(4,6)])
def test_unequal(self,lattice,shape):
R = Rotation.from_random(shape)
assert not ( Orientation(R,lattice) != Orientation(R,lattice) if shape is None else \
(Orientation(R,lattice) != Orientation(R,lattice)).any())
@pytest.mark.parametrize('lattice',Orientation.crystal_families)
@pytest.mark.parametrize('lattice',crystal_families)
@pytest.mark.parametrize('shape',[None,5,(4,6)])
def test_close(self,lattice,shape):
R = Orientation.from_random(lattice=lattice,shape=shape)
@ -182,14 +183,14 @@ class TestOrientation:
with pytest.raises(ValueError):
Orientation(lattice='aP',a=1,b=2,c=3,alpha=45,beta=45,gamma=90.0001,degrees=True) # noqa
@pytest.mark.parametrize('lattice',Orientation.crystal_families)
@pytest.mark.parametrize('lattice',crystal_families)
@pytest.mark.parametrize('angle',[10,20,30,40])
def test_average(self,angle,lattice):
o = Orientation.from_axis_angle(lattice=lattice,axis_angle=[[0,0,1,10],[0,0,1,angle]],degrees=True)
avg_angle = o.average().as_axis_angle(degrees=True,pair=True)[1]
assert np.isclose(avg_angle,10+(angle-10)/2.)
@pytest.mark.parametrize('lattice',Orientation.crystal_families)
@pytest.mark.parametrize('lattice',crystal_families)
def test_reduced_equivalent(self,lattice):
i = Orientation(lattice=lattice)
o = Orientation.from_random(lattice=lattice)
@ -197,7 +198,7 @@ class TestOrientation:
FZ = np.argmin(abs(eq.misorientation(i.broadcast_to(len(eq))).as_axis_angle(pair=True)[1]))
assert o.reduced == eq[FZ]
@pytest.mark.parametrize('lattice',Orientation.crystal_families)
@pytest.mark.parametrize('lattice',crystal_families)
@pytest.mark.parametrize('N',[1,8,32])
def test_disorientation(self,lattice,N):
o = Orientation.from_random(lattice=lattice,shape=N)
@ -215,7 +216,7 @@ class TestOrientation:
.misorientation(p[n].equivalent[ops[n][1]])
.as_quaternion())
@pytest.mark.parametrize('lattice',Orientation.crystal_families)
@pytest.mark.parametrize('lattice',crystal_families)
@pytest.mark.parametrize('a,b',[
((2,3,2),(2,3,2)),
((2,2),(4,4)),
@ -230,20 +231,20 @@ class TestOrientation:
assert o[tuple(loc[:len(o.shape)])].disorientation(p[tuple(loc[-len(p.shape):])]) \
.isclose(o.disorientation(p)[tuple(loc)])
@pytest.mark.parametrize('lattice',Orientation.crystal_families)
@pytest.mark.parametrize('lattice',crystal_families)
def test_disorientation360(self,lattice):
o_1 = Orientation(Rotation(),lattice)
o_2 = Orientation.from_Euler_angles(lattice=lattice,phi=[360,0,0],degrees=True)
assert np.allclose((o_1.disorientation(o_2)).as_matrix(),np.eye(3))
@pytest.mark.parametrize('lattice',Orientation.crystal_families)
@pytest.mark.parametrize('lattice',crystal_families)
@pytest.mark.parametrize('shape',[(1),(2,3),(4,3,2)])
def test_reduced_vectorization(self,lattice,shape):
o = Orientation.from_random(lattice=lattice,shape=shape)
for r, theO in zip(o.reduced.flatten(),o.flatten()):
assert r == theO.reduced
@pytest.mark.parametrize('lattice',Orientation.crystal_families)
@pytest.mark.parametrize('lattice',crystal_families)
def test_reduced_corner_cases(self,lattice):
# test whether there is always a sym-eq rotation that falls into the FZ
N = np.random.randint(10,40)
@ -253,7 +254,7 @@ class TestOrientation:
assert evenly_distributed.shape == evenly_distributed.reduced.shape
@pytest.mark.parametrize('lattice',Orientation.crystal_families)
@pytest.mark.parametrize('lattice',crystal_families)
@pytest.mark.parametrize('shape',[(1),(2,3),(4,3,2)])
@pytest.mark.parametrize('vector',np.array([[1,0,0],[1,2,3],[-1,1,-1]]))
@pytest.mark.parametrize('proper',[True,False])
@ -262,7 +263,7 @@ class TestOrientation:
for r, theO in zip(o.to_SST(vector=vector,proper=proper).reshape((-1,3)),o.flatten()):
assert np.allclose(r,theO.to_SST(vector=vector,proper=proper))
@pytest.mark.parametrize('lattice',Orientation.crystal_families)
@pytest.mark.parametrize('lattice',crystal_families)
@pytest.mark.parametrize('shape',[(1),(2,3),(4,3,2)])
@pytest.mark.parametrize('vector',np.array([[1,0,0],[1,2,3],[-1,1,-1]]))
@pytest.mark.parametrize('proper',[True,False])
@ -272,7 +273,7 @@ class TestOrientation:
for r, theO in zip(o.IPF_color(vector,in_SST=in_SST,proper=proper).reshape((-1,3)),o.flatten()):
assert np.allclose(r,theO.IPF_color(vector,in_SST=in_SST,proper=proper))
@pytest.mark.parametrize('lattice',Orientation.crystal_families)
@pytest.mark.parametrize('lattice',crystal_families)
@pytest.mark.parametrize('a,b',[
((2,3,2),(2,3,2)),
((2,2),(4,4)),
@ -300,7 +301,7 @@ class TestOrientation:
assert np.allclose(np.array(color['RGB']),
cube.IPF_color(vector=np.array(direction),proper=proper))
@pytest.mark.parametrize('lattice',Orientation.crystal_families)
@pytest.mark.parametrize('lattice',crystal_families)
@pytest.mark.parametrize('proper',[True,False])
def test_IPF_equivalent(self,set_of_quaternions,lattice,proper):
direction = np.random.random(3)*2.0-1.0
@ -308,13 +309,13 @@ class TestOrientation:
color = o.IPF_color(vector=direction,proper=proper)
assert np.allclose(np.broadcast_to(color[0,...],color.shape),color)
@pytest.mark.parametrize('lattice',Orientation.crystal_families)
@pytest.mark.parametrize('lattice',crystal_families)
def test_in_FZ_vectorization(self,set_of_rodrigues,lattice):
result = Orientation.from_Rodrigues_vector(rho=set_of_rodrigues.reshape((-1,4,4)),lattice=lattice).in_FZ.reshape(-1)
for r,rho in zip(result,set_of_rodrigues[:len(result)]):
assert r == Orientation.from_Rodrigues_vector(rho=rho,lattice=lattice).in_FZ
@pytest.mark.parametrize('lattice',Orientation.crystal_families)
@pytest.mark.parametrize('lattice',crystal_families)
def test_in_disorientation_FZ_vectorization(self,set_of_rodrigues,lattice):
result = Orientation.from_Rodrigues_vector(rho=set_of_rodrigues.reshape((-1,4,4)),
lattice=lattice).in_disorientation_FZ.reshape(-1)
@ -322,7 +323,7 @@ class TestOrientation:
assert r == Orientation.from_Rodrigues_vector(rho=rho,lattice=lattice).in_disorientation_FZ
@pytest.mark.parametrize('proper',[True,False])
@pytest.mark.parametrize('lattice',Orientation.crystal_families)
@pytest.mark.parametrize('lattice',crystal_families)
def test_in_SST_vectorization(self,lattice,proper):
vecs = np.random.rand(20,4,3)
result = Orientation(lattice=lattice).in_SST(vecs,proper).flatten()
@ -393,7 +394,7 @@ class TestOrientation:
a=a,b=b,c=c,
alpha=alpha,beta=beta,gamma=gamma).related(relation) # noqa
@pytest.mark.parametrize('lattice',Orientation.crystal_families)
@pytest.mark.parametrize('lattice',crystal_families)
@pytest.mark.parametrize('proper',[True,False])
def test_in_SST(self,lattice,proper):
assert Orientation(lattice=lattice).in_SST(np.zeros(3),proper)

View File

@ -14,23 +14,23 @@ list(FILTER damask-sources EXCLUDE REGEX ".*commercialFEM_fileList.*.f90")
if (PROJECT_NAME STREQUAL "damask-grid")
file(GLOB grid-sources grid/*.f90)
file(GLOB grid-sources grid/*.f90)
if(NOT CMAKE_BUILD_TYPE STREQUAL "SYNTAXONLY")
if (NOT CMAKE_BUILD_TYPE STREQUAL "SYNTAXONLY")
add_executable(DAMASK_grid ${damask-sources} ${grid-sources})
install (TARGETS DAMASK_grid RUNTIME DESTINATION bin)
else()
else ()
add_library(DAMASK_grid OBJECT ${damask-sources} ${grid-sources})
exec_program (mktemp OUTPUT_VARIABLE nothing)
exec_program (mktemp ARGS -d OUTPUT_VARIABLE black_hole)
install (PROGRAMS ${nothing} DESTINATION ${black_hole})
endif()
endif ()
elseif (PROJECT_NAME STREQUAL "damask-mesh")
file(GLOB mesh-sources mesh/*.f90)
file(GLOB mesh-sources mesh/*.f90)
add_executable(DAMASK_mesh ${damask-sources} ${mesh-sources})
install (TARGETS DAMASK_mesh RUNTIME DESTINATION bin)
add_executable(DAMASK_mesh ${damask-sources} ${mesh-sources})
install (TARGETS DAMASK_mesh RUNTIME DESTINATION bin)
endif()
endif ()

View File

@ -96,15 +96,10 @@ subroutine DAMASK_interface_init
print'(/,a)', ' Version: '//DAMASKVERSION
! https://github.com/jeffhammond/HPCInfo/blob/master/docs/Preprocessor-Macros.md
#if defined(__PGI)
print'(/,a,i4.4,a,i8.8)', ' Compiled with PGI fortran version :', __PGIC__,&
'.', __PGIC_MINOR__
#else
print'(/,a)', ' Compiled with: '//compiler_version()
print'(a)', ' Compiler options: '//compiler_options()
#endif
! https://github.com/jeffhammond/HPCInfo/blob/master/docs/Preprocessor-Macros.md
print'(/,a)', ' Compiled on: '//__DATE__//' at '//__TIME__
print'(/,a,i0,a,i0,a,i0)', &

View File

@ -113,8 +113,7 @@ end function IO_readlines
!--------------------------------------------------------------------------------------------------
!> @brief Read whole file.
!> @details ensures that the string ends with a new line (expected UNIX behavior) and rejects
! windows (CRLF) line endings
!> @details ensures that the string ends with a new line (expected UNIX behavior)
!--------------------------------------------------------------------------------------------------
function IO_read(fileName) result(fileContent)
@ -124,8 +123,7 @@ function IO_read(fileName) result(fileContent)
integer :: &
fileLength, &
fileUnit, &
myStat, &
firstEOL
myStat
character, parameter :: CR = achar(13)
@ -143,12 +141,17 @@ function IO_read(fileName) result(fileContent)
if(myStat /= 0) call IO_error(102,ext_msg=trim(fileName))
close(fileUnit)
foundCRLF: if (scan(fileContent(:index(fileContent,IO_EOL)),CR) /= 0) then
CRLF2LF: block
integer :: c
do c=1, len(fileContent)
if (fileContent(c:c) == CR) fileContent(c:c) = ' '
enddo
end block CRLF2LF
endif foundCRLF
if(fileContent(fileLength:fileLength) /= IO_EOL) fileContent = fileContent//IO_EOL ! ensure EOL@EOF
firstEOL = index(fileContent,IO_EOL)
if(scan(fileContent(firstEOL:firstEOL),CR) /= 0) call IO_error(115)
end function IO_read
@ -400,9 +403,6 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg)
msg = 'invalid character for logical:'
case (114)
msg = 'cannot decode base64 string:'
case (115)
msg = 'found CR. Windows file endings (CRLF) are not supported.'
!--------------------------------------------------------------------------------------------------
! lattice error messages
@ -644,11 +644,16 @@ subroutine selfTest
if(dNeq(1.0_pReal, IO_stringAsFloat('1.0'))) error stop 'IO_stringAsFloat'
if(dNeq(1.0_pReal, IO_stringAsFloat('1e0'))) error stop 'IO_stringAsFloat'
if(dNeq(0.1_pReal, IO_stringAsFloat('1e-1'))) error stop 'IO_stringAsFloat'
if(dNeq(0.1_pReal, IO_stringAsFloat('1.0e-1'))) error stop 'IO_stringAsFloat'
if(dNeq(0.1_pReal, IO_stringAsFloat('1.00e-1'))) error stop 'IO_stringAsFloat'
if(dNeq(10._pReal, IO_stringAsFloat(' 1.0e+1 '))) error stop 'IO_stringAsFloat'
if(3112019 /= IO_stringAsInt( '3112019')) error stop 'IO_stringAsInt'
if(3112019 /= IO_stringAsInt(' 3112019')) error stop 'IO_stringAsInt'
if(-3112019 /= IO_stringAsInt('-3112019')) error stop 'IO_stringAsInt'
if(3112019 /= IO_stringAsInt('+3112019 ')) error stop 'IO_stringAsInt'
if(3112019 /= IO_stringAsInt('03112019 ')) error stop 'IO_stringAsInt'
if(3112019 /= IO_stringAsInt('+03112019')) error stop 'IO_stringAsInt'
if(.not. IO_stringAsBool(' true')) error stop 'IO_stringAsBool'
if(.not. IO_stringAsBool(' True ')) error stop 'IO_stringAsBool'

View File

@ -21,6 +21,7 @@
#include "lattice.f90"
#include "phase.f90"
#include "phase_mechanical.f90"
#include "phase_mechanical_elastic.f90"
#include "phase_mechanical_plastic.f90"
#include "phase_mechanical_plastic_none.f90"
#include "phase_mechanical_plastic_isotropic.f90"

View File

@ -31,44 +31,38 @@ module lattice
FCC_NCLEAVAGESYSTEM = [3] !< # of cleavage systems per family for fcc
integer, parameter :: &
#ifndef __PGI
FCC_NSLIP = sum(FCC_NSLIPSYSTEM), & !< total # of slip systems for fcc
FCC_NTWIN = sum(FCC_NTWINSYSTEM), & !< total # of twin systems for fcc
FCC_NTRANS = sum(FCC_NTRANSSYSTEM), & !< total # of transformation systems for fcc
FCC_NCLEAVAGE = sum(FCC_NCLEAVAGESYSTEM) !< total # of cleavage systems for fcc
#else
FCC_NSLIP = 18, &
FCC_NTWIN = 12, &
FCC_NTRANS = 12, &
FCC_NCLEAVAGE = 3
#endif
real(pReal), dimension(3+3,FCC_NSLIP), parameter :: &
FCC_SYSTEMSLIP = reshape(real([&
! Slip direction Plane normal ! SCHMID-BOAS notation
0, 1,-1, 1, 1, 1, & ! B2
-1, 0, 1, 1, 1, 1, & ! B4
1,-1, 0, 1, 1, 1, & ! B5
0,-1,-1, -1,-1, 1, & ! C1
1, 0, 1, -1,-1, 1, & ! C3
-1, 1, 0, -1,-1, 1, & ! C5
0,-1, 1, 1,-1,-1, & ! A2
-1, 0,-1, 1,-1,-1, & ! A3
1, 1, 0, 1,-1,-1, & ! A6
0, 1, 1, -1, 1,-1, & ! D1
1, 0,-1, -1, 1,-1, & ! D4
-1,-1, 0, -1, 1,-1, & ! D6
! Slip system <110>{110}
! <110>{111} systems
0, 1,-1, 1, 1, 1, & ! B2
-1, 0, 1, 1, 1, 1, & ! B4
1,-1, 0, 1, 1, 1, & ! B5
0,-1,-1, -1,-1, 1, & ! C1
1, 0, 1, -1,-1, 1, & ! C3
-1, 1, 0, -1,-1, 1, & ! C5
0,-1, 1, 1,-1,-1, & ! A2
-1, 0,-1, 1,-1,-1, & ! A3
1, 1, 0, 1,-1,-1, & ! A6
0, 1, 1, -1, 1,-1, & ! D1
1, 0,-1, -1, 1,-1, & ! D4
-1,-1, 0, -1, 1,-1, & ! D6
! <110>{110}/non-octahedral systems
1, 1, 0, 1,-1, 0, &
1,-1, 0, 1, 1, 0, &
1, 0, 1, 1, 0,-1, &
1, 0,-1, 1, 0, 1, &
0, 1, 1, 0, 1,-1, &
0, 1,-1, 0, 1, 1 &
],pReal),shape(FCC_SYSTEMSLIP)) !< Slip system <110>{111} directions. Sorted according to Eisenlohr & Hantcherli
],pReal),shape(FCC_SYSTEMSLIP)) !< fcc slip systems
real(pReal), dimension(3+3,FCC_NTWIN), parameter :: &
FCC_SYSTEMTWIN = reshape(real( [&
! <112>{111} systems
-2, 1, 1, 1, 1, 1, &
1,-2, 1, 1, 1, 1, &
1, 1,-2, 1, 1, 1, &
@ -81,7 +75,7 @@ module lattice
2, 1,-1, -1, 1,-1, &
-1,-2,-1, -1, 1,-1, &
-1, 1, 2, -1, 1,-1 &
],pReal),shape(FCC_SYSTEMTWIN)) !< Twin system <112>{111} directions. Sorted according to Eisenlohr & Hantcherli
],pReal),shape(FCC_SYSTEMTWIN)) !< fcc twin systems
integer, dimension(2,FCC_NTWIN), parameter, public :: &
lattice_FCC_TWINNUCLEATIONSLIPPAIR = reshape( [&
@ -101,16 +95,16 @@ module lattice
real(pReal), dimension(3+3,FCC_NCLEAVAGE), parameter :: &
FCC_SYSTEMCLEAVAGE = reshape(real([&
! Cleavage direction Plane normal
! <001>{001} systems
0, 1, 0, 1, 0, 0, &
0, 0, 1, 0, 1, 0, &
1, 0, 0, 0, 0, 1 &
],pReal),shape(FCC_SYSTEMCLEAVAGE))
],pReal),shape(FCC_SYSTEMCLEAVAGE)) !< fcc cleavage systems
!--------------------------------------------------------------------------------------------------
! body centered cubic (cI)
integer, dimension(*), parameter :: &
BCC_NSLIPSYSTEM = [12, 12] !< # of slip systems per family for bcc
BCC_NSLIPSYSTEM = [12, 12, 24] !< # of slip systems per family for bcc
integer, dimension(*), parameter :: &
BCC_NTWINSYSTEM = [12] !< # of twin systems per family for bcc
@ -119,50 +113,68 @@ module lattice
BCC_NCLEAVAGESYSTEM = [3] !< # of cleavage systems per family for bcc
integer, parameter :: &
#ifndef __PGI
BCC_NSLIP = sum(BCC_NSLIPSYSTEM), & !< total # of slip systems for bcc
BCC_NTWIN = sum(BCC_NTWINSYSTEM), & !< total # of twin systems for bcc
BCC_NCLEAVAGE = sum(BCC_NCLEAVAGESYSTEM) !< total # of cleavage systems for bcc
#else
BCC_NSLIP = 24, &
BCC_NTWIN = 12, &
BCC_NCLEAVAGE = 3
#endif
real(pReal), dimension(3+3,BCC_NSLIP), parameter :: &
BCC_SYSTEMSLIP = reshape(real([&
! Slip direction Plane normal
! Slip system <111>{110}
1,-1, 1, 0, 1, 1, &
-1,-1, 1, 0, 1, 1, &
1, 1, 1, 0,-1, 1, &
-1, 1, 1, 0,-1, 1, &
-1, 1, 1, 1, 0, 1, &
-1,-1, 1, 1, 0, 1, &
1, 1, 1, -1, 0, 1, &
1,-1, 1, -1, 0, 1, &
-1, 1, 1, 1, 1, 0, &
-1, 1,-1, 1, 1, 0, &
1, 1, 1, -1, 1, 0, &
1, 1,-1, -1, 1, 0, &
! Slip system <111>{112}
-1, 1, 1, 2, 1, 1, &
1, 1, 1, -2, 1, 1, &
1, 1,-1, 2,-1, 1, &
1,-1, 1, 2, 1,-1, &
1,-1, 1, 1, 2, 1, &
1, 1,-1, -1, 2, 1, &
1, 1, 1, 1,-2, 1, &
-1, 1, 1, 1, 2,-1, &
1, 1,-1, 1, 1, 2, &
1,-1, 1, -1, 1, 2, &
-1, 1, 1, 1,-1, 2, &
1, 1, 1, 1, 1,-2 &
],pReal),shape(BCC_SYSTEMSLIP))
! <111>{110} systems
1,-1, 1, 0, 1, 1, & ! D1
-1,-1, 1, 0, 1, 1, & ! C1
1, 1, 1, 0,-1, 1, & ! B2
-1, 1, 1, 0,-1, 1, & ! A2
-1, 1, 1, 1, 0, 1, & ! A3
-1,-1, 1, 1, 0, 1, & ! C3
1, 1, 1, -1, 0, 1, & ! B4
1,-1, 1, -1, 0, 1, & ! D4
-1, 1, 1, 1, 1, 0, & ! A6
-1, 1,-1, 1, 1, 0, & ! D6
1, 1, 1, -1, 1, 0, & ! B5
1, 1,-1, -1, 1, 0, & ! C5
! <111>{112} systems
-1, 1, 1, 2, 1, 1, & ! A-4
1, 1, 1, -2, 1, 1, & ! B-3
1, 1,-1, 2,-1, 1, & ! C-10
1,-1, 1, 2, 1,-1, & ! D-9
1,-1, 1, 1, 2, 1, & ! D-6
1, 1,-1, -1, 2, 1, & ! C-5
1, 1, 1, 1,-2, 1, & ! B-12
-1, 1, 1, 1, 2,-1, & ! A-11
1, 1,-1, 1, 1, 2, & ! C-2
1,-1, 1, -1, 1, 2, & ! D-1
-1, 1, 1, 1,-1, 2, & ! A-8
1, 1, 1, 1, 1,-2, & ! B-7
! Slip system <111>{123}
1, 1,-1, 1, 2, 3, &
1,-1, 1, -1, 2, 3, &
-1, 1, 1, 1,-2, 3, &
1, 1, 1, 1, 2,-3, &
1,-1, 1, 1, 3, 2, &
1, 1,-1, -1, 3, 2, &
1, 1, 1, 1,-3, 2, &
-1, 1, 1, 1, 3,-2, &
1, 1,-1, 2, 1, 3, &
1,-1, 1, -2, 1, 3, &
-1, 1, 1, 2,-1, 3, &
1, 1, 1, 2, 1,-3, &
1,-1, 1, 2, 3, 1, &
1, 1,-1, -2, 3, 1, &
1, 1, 1, 2,-3, 1, &
-1, 1, 1, 2, 3,-1, &
-1, 1, 1, 3, 1, 2, &
1, 1, 1, -3, 1, 2, &
1, 1,-1, 3,-1, 2, &
1,-1, 1, 3, 1,-2, &
-1, 1, 1, 3, 2, 1, &
1, 1, 1, -3, 2, 1, &
1, 1,-1, 3,-2, 1, &
1,-1, 1, 3, 2,-1 &
],pReal),shape(BCC_SYSTEMSLIP)) !< bcc slip systems
real(pReal), dimension(3+3,BCC_NTWIN), parameter :: &
BCC_SYSTEMTWIN = reshape(real([&
! Twin system <111>{112}
! <111>{112} systems
-1, 1, 1, 2, 1, 1, &
1, 1, 1, -2, 1, 1, &
1, 1,-1, 2,-1, 1, &
@ -175,15 +187,15 @@ module lattice
1,-1, 1, -1, 1, 2, &
-1, 1, 1, 1,-1, 2, &
1, 1, 1, 1, 1,-2 &
],pReal),shape(BCC_SYSTEMTWIN))
],pReal),shape(BCC_SYSTEMTWIN)) !< bcc twin systems
real(pReal), dimension(3+3,BCC_NCLEAVAGE), parameter :: &
BCC_SYSTEMCLEAVAGE = reshape(real([&
! Cleavage direction Plane normal
! <001>{001} systems
0, 1, 0, 1, 0, 0, &
0, 0, 1, 0, 1, 0, &
1, 0, 0, 0, 0, 1 &
],pReal),shape(BCC_SYSTEMCLEAVAGE))
],pReal),shape(BCC_SYSTEMCLEAVAGE)) !< bcc cleavage systems
!--------------------------------------------------------------------------------------------------
! hexagonal (hP)
@ -194,37 +206,31 @@ module lattice
HEX_NTWINSYSTEM = [6, 6, 6, 6] !< # of slip systems per family for hex
integer, parameter :: &
#ifndef __PGI
HEX_NSLIP = sum(HEX_NSLIPSYSTEM), & !< total # of slip systems for hex
HEX_NTWIN = sum(HEX_NTWINSYSTEM) !< total # of twin systems for hex
#else
HEX_NSLIP = 33, &
HEX_NTWIN = 24
#endif
real(pReal), dimension(4+4,HEX_NSLIP), parameter :: &
HEX_SYSTEMSLIP = reshape(real([&
! Slip direction Plane normal
! Basal systems <-1-1.0>{00.1} (independent of c/a-ratio, Bravais notation (4 coordinate base))
! <-1-1.0>{00.1}/basal systems (independent of c/a-ratio)
2, -1, -1, 0, 0, 0, 0, 1, &
-1, 2, -1, 0, 0, 0, 0, 1, &
-1, -1, 2, 0, 0, 0, 0, 1, &
! 1st type prismatic systems <-1-1.0>{1-1.0} (independent of c/a-ratio)
! <-1-1.0>{1-1.0}/prismatic systems (independent of c/a-ratio)
2, -1, -1, 0, 0, 1, -1, 0, &
-1, 2, -1, 0, -1, 0, 1, 0, &
-1, -1, 2, 0, 1, -1, 0, 0, &
! 2nd type prismatic systems <-11.0>{11.0} -- a slip; plane normals independent of c/a-ratio
! <-11.0>{11.0}/2nd order prismatic compound systems (plane normal independent of c/a-ratio)
-1, 1, 0, 0, 1, 1, -2, 0, &
0, -1, 1, 0, -2, 1, 1, 0, &
1, 0, -1, 0, 1, -2, 1, 0, &
! 1st type 1st order pyramidal systems <-1-1.0>{-11.1} -- plane normals depend on the c/a-ratio
! <-1-1.0>{-11.1}/1st order pyramidal <a> systems (direction independent of c/a-ratio)
-1, 2, -1, 0, 1, 0, -1, 1, &
-2, 1, 1, 0, 0, 1, -1, 1, &
-1, -1, 2, 0, -1, 1, 0, 1, &
1, -2, 1, 0, -1, 0, 1, 1, &
2, -1, -1, 0, 0, -1, 1, 1, &
1, 1, -2, 0, 1, -1, 0, 1, &
! pyramidal system: c+a slip <11.3>{-10.1} -- plane normals depend on the c/a-ratio
! <11.3>{-10.1}/1st order pyramidal <c+a> systems (direction independent of c/a-ratio)
-2, 1, 1, 3, 1, 0, -1, 1, &
-1, -1, 2, 3, 1, 0, -1, 1, &
-1, -1, 2, 3, 0, 1, -1, 1, &
@ -237,96 +243,96 @@ module lattice
-1, 2, -1, 3, 0, -1, 1, 1, &
-1, 2, -1, 3, 1, -1, 0, 1, &
-2, 1, 1, 3, 1, -1, 0, 1, &
! pyramidal system: c+a slip <11.3>{-1-1.2} -- as for hexagonal ice (Castelnau et al. 1996, similar to twin system found below)
-1, -1, 2, 3, 1, 1, -2, 2, & ! <11.3>{-1-1.2} shear = 2((c/a)^2-2)/(3 c/a)
! <11.3>{-1-1.2}/2nd order pyramidal <c+a> systems
-1, -1, 2, 3, 1, 1, -2, 2, &
1, -2, 1, 3, -1, 2, -1, 2, &
2, -1, -1, 3, -2, 1, 1, 2, &
1, 1, -2, 3, -1, -1, 2, 2, &
-1, 2, -1, 3, 1, -2, 1, 2, &
-2, 1, 1, 3, 2, -1, -1, 2 &
],pReal),shape(HEX_SYSTEMSLIP)) !< slip systems for hex, sorted by P. Eisenlohr CCW around <c> starting next to a_1 axis
],pReal),shape(HEX_SYSTEMSLIP)) !< hex slip systems, sorted by P. Eisenlohr CCW around <c> starting next to a_1 axis
real(pReal), dimension(4+4,HEX_NTWIN), parameter :: &
HEX_SYSTEMTWIN = reshape(real([&
! Compression or Tension = f(twinning shear=f(c/a)) for each metal ! (according to Yoo 1981)
-1, 0, 1, 1, 1, 0, -1, 2, & ! <-10.1>{10.2} shear = (3-(c/a)^2)/(sqrt(3) c/a)
! <-10.1>{10.2} systems, shear = (3-(c/a)^2)/(sqrt(3) c/a)
! tension in Co, Mg, Zr, Ti, and Be; compression in Cd and Zn
-1, 0, 1, 1, 1, 0, -1, 2, & !
0, -1, 1, 1, 0, 1, -1, 2, &
1, -1, 0, 1, -1, 1, 0, 2, &
1, 0, -1, 1, -1, 0, 1, 2, &
0, 1, -1, 1, 0, -1, 1, 2, &
-1, 1, 0, 1, 1, -1, 0, 2, &
!
-1, -1, 2, 6, 1, 1, -2, 1, & ! <11.6>{-1-1.1} shear = 1/(c/a)
! <11.6>{-1-1.1} systems, shear = 1/(c/a)
! tension in Co, Re, and Zr
-1, -1, 2, 6, 1, 1, -2, 1, &
1, -2, 1, 6, -1, 2, -1, 1, &
2, -1, -1, 6, -2, 1, 1, 1, &
1, 1, -2, 6, -1, -1, 2, 1, &
-1, 2, -1, 6, 1, -2, 1, 1, &
-2, 1, 1, 6, 2, -1, -1, 1, &
!
1, 0, -1, -2, 1, 0, -1, 1, & ! <10.-2>{10.1} shear = (4(c/a)^2-9)/(4 sqrt(3) c/a)
! <10.-2>{10.1} systems, shear = (4(c/a)^2-9)/(4 sqrt(3) c/a)
! compression in Mg
1, 0, -1, -2, 1, 0, -1, 1, &
0, 1, -1, -2, 0, 1, -1, 1, &
-1, 1, 0, -2, -1, 1, 0, 1, &
-1, 0, 1, -2, -1, 0, 1, 1, &
0, -1, 1, -2, 0, -1, 1, 1, &
1, -1, 0, -2, 1, -1, 0, 1, &
!
1, 1, -2, -3, 1, 1, -2, 2, & ! <11.-3>{11.2} shear = 2((c/a)^2-2)/(3 c/a)
! <11.-3>{11.2} systems, shear = 2((c/a)^2-2)/(3 c/a)
! compression in Ti and Zr
1, 1, -2, -3, 1, 1, -2, 2, &
-1, 2, -1, -3, -1, 2, -1, 2, &
-2, 1, 1, -3, -2, 1, 1, 2, &
-1, -1, 2, -3, -1, -1, 2, 2, &
1, -2, 1, -3, 1, -2, 1, 2, &
2, -1, -1, -3, 2, -1, -1, 2 &
],pReal),shape(HEX_SYSTEMTWIN)) !< twin systems for hex, sorted by P. Eisenlohr CCW around <c> starting next to a_1 axis
],pReal),shape(HEX_SYSTEMTWIN)) !< hex twin systems, sorted by P. Eisenlohr CCW around <c> starting next to a_1 axis
!--------------------------------------------------------------------------------------------------
! body centered tetragonal (tI)
integer, dimension(*), parameter :: &
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
BCT_NSLIPSYSTEM = [2, 2, 2, 4, 2, 4, 2, 2, 4, 8, 4, 8, 8 ] !< # of slip systems per family for bct
integer, parameter :: &
#ifndef __PGI
BCT_NSLIP = sum(BCT_NSLIPSYSTEM) !< total # of slip systems for bct
#else
BCT_NSLIP = 52
#endif
real(pReal), dimension(3+3,BCT_NSLIP), parameter :: &
BCT_SYSTEMSLIP = reshape(real([&
! Slip direction Plane normal
! Slip family 1 {100)<001] (Bravais notation {hkl)<uvw] for bct c/a = 0.5456)
! {100)<001] systems
0, 0, 1, 1, 0, 0, &
0, 0, 1, 0, 1, 0, &
! Slip family 2 {110)<001]
! {110)<001] systems
0, 0, 1, 1, 1, 0, &
0, 0, 1, -1, 1, 0, &
! slip family 3 {100)<010]
! {100)<010] systems
0, 1, 0, 1, 0, 0, &
1, 0, 0, 0, 1, 0, &
! Slip family 4 {110)<1-11]/2
! {110)<1-11]/2 systems
1,-1, 1, 1, 1, 0, &
1,-1,-1, 1, 1, 0, &
-1,-1,-1, -1, 1, 0, &
-1,-1, 1, -1, 1, 0, &
! Slip family 5 {110)<1-10]
! {110)<1-10] systems
1, -1, 0, 1, 1, 0, &
1, 1, 0, 1,-1, 0, &
! Slip family 6 {100)<011]
! {100)<011] systems
0, 1, 1, 1, 0, 0, &
0,-1, 1, 1, 0, 0, &
-1, 0, 1, 0, 1, 0, &
1, 0, 1, 0, 1, 0, &
! Slip family 7 {001)<010]
! {001)<010] systems
0, 1, 0, 0, 0, 1, &
1, 0, 0, 0, 0, 1, &
! Slip family 8 {001)<110]
! {001)<110] systems
1, 1, 0, 0, 0, 1, &
-1, 1, 0, 0, 0, 1, &
! Slip family 9 {011)<01-1]
! {011)<01-1] systems
0, 1,-1, 0, 1, 1, &
0,-1,-1, 0,-1, 1, &
-1, 0,-1, -1, 0, 1, &
1, 0,-1, 1, 0, 1, &
! Slip family 10 {011)<1-11]/2
! {011)<1-11]/2 systems
1,-1, 1, 0, 1, 1, &
1, 1,-1, 0, 1, 1, &
1, 1, 1, 0, 1,-1, &
@ -335,12 +341,12 @@ module lattice
-1,-1, 1, 1, 0, 1, &
1, 1, 1, 1, 0,-1, &
1,-1, 1, 1, 0,-1, &
! Slip family 11 {011)<100]
! {011)<100] systems
1, 0, 0, 0, 1, 1, &
1, 0, 0, 0, 1,-1, &
0, 1, 0, 1, 0, 1, &
0, 1, 0, 1, 0,-1, &
! Slip family 12 {211)<01-1]
! {211)<01-1] systems
0, 1,-1, 2, 1, 1, &
0,-1,-1, 2,-1, 1, &
1, 0,-1, 1, 2, 1, &
@ -349,7 +355,7 @@ module lattice
0,-1,-1, -2,-1, 1, &
-1, 0,-1, -1,-2, 1, &
1, 0,-1, 1,-2, 1, &
! Slip family 13 {211)<-111]/2
! {211)<-111]/2 systems
-1, 1, 1, 2, 1, 1, &
-1,-1, 1, 2,-1, 1, &
1,-1, 1, 1, 2, 1, &
@ -358,27 +364,22 @@ module lattice
1,-1, 1, -2,-1, 1, &
-1, 1, 1, -1,-2, 1, &
1, 1, 1, 1,-2, 1 &
],pReal),shape(BCT_SYSTEMSLIP)) !< slip systems for bct sorted by Bieler
],pReal),shape(BCT_SYSTEMSLIP)) !< bct slip systems for c/a = 0.5456 (Sn), sorted by Bieler 2009 (https://doi.org/10.1007/s11664-009-0909-x)
!--------------------------------------------------------------------------------------------------
! orthorhombic primitive (oP)
integer, dimension(*), parameter :: &
ORT_NCLEAVAGESYSTEM = [1, 1, 1] !< # of cleavage systems per family for ortho
ORT_NCLEAVAGESYSTEM = [1, 1, 1] !< # of cleavage systems per family for orthorhombic primitive
integer, parameter :: &
#ifndef __PGI
ORT_NCLEAVAGE = sum(ORT_NCLEAVAGESYSTEM) !< total # of cleavage systems for ortho
#else
ORT_NCLEAVAGE = 3
#endif
ORT_NCLEAVAGE = sum(ORT_NCLEAVAGESYSTEM) !< total # of cleavage systems for orthorhombic primitive
real(pReal), dimension(3+3,ORT_NCLEAVAGE), parameter :: &
ORT_SYSTEMCLEAVAGE = reshape(real([&
! Cleavage direction Plane normal
0, 1, 0, 1, 0, 0, &
0, 0, 1, 0, 1, 0, &
1, 0, 0, 0, 0, 1 &
],pReal),shape(ORT_SYSTEMCLEAVAGE))
],pReal),shape(ORT_SYSTEMCLEAVAGE)) !< orthorhombic primitive cleavage systems
enum, bind(c); enumerator :: &
@ -714,10 +715,9 @@ function lattice_nonSchmidMatrix(Nslip,nonSchmidCoefficients,sense) result(nonSc
if (abs(sense) /= 1) error stop 'Sense in lattice_nonSchmidMatrix'
coordinateSystem = buildCoordinateSystem(Nslip,BCC_NSLIPSYSTEM,BCC_SYSTEMSLIP,&
'cI',0.0_pReal)
coordinateSystem = buildCoordinateSystem(Nslip,BCC_NSLIPSYSTEM,BCC_SYSTEMSLIP,'cI',0.0_pReal)
coordinateSystem(1:3,1,1:sum(Nslip)) = coordinateSystem(1:3,1,1:sum(Nslip))*real(sense,pReal) ! convert unidirectional coordinate system
nonSchmidMatrix = lattice_SchmidMatrix_slip(Nslip,'cI',0.0_pReal) ! Schmid contribution
nonSchmidMatrix = lattice_SchmidMatrix_slip(Nslip,'cI',0.0_pReal) ! Schmid contribution
do i = 1,sum(Nslip)
direction = coordinateSystem(1:3,1,i)
@ -759,81 +759,112 @@ function lattice_interaction_SlipBySlip(Nslip,interactionValues,structure) resul
integer, dimension(FCC_NSLIP,FCC_NSLIP), parameter :: &
FCC_INTERACTIONSLIPSLIP = reshape( [&
1, 2, 2, 4, 6, 5, 3, 5, 5, 4, 5, 6, 9,10, 9,10,11,12, & ! -----> acting
2, 1, 2, 6, 4, 5, 5, 4, 6, 5, 3, 5, 9,10,11,12, 9,10, & ! |
2, 2, 1, 5, 5, 3, 5, 6, 4, 6, 5, 4, 11,12, 9,10, 9,10, & ! |
4, 6, 5, 1, 2, 2, 4, 5, 6, 3, 5, 5, 9,10,10, 9,12,11, & ! v
6, 4, 5, 2, 1, 2, 5, 3, 5, 5, 4, 6, 9,10,12,11,10, 9, & ! reacting
5, 5, 3, 2, 2, 1, 6, 5, 4, 5, 6, 4, 11,12,10, 9,10, 9, &
3, 5, 5, 4, 5, 6, 1, 2, 2, 4, 6, 5, 10, 9,10, 9,11,12, &
5, 4, 6, 5, 3, 5, 2, 1, 2, 6, 4, 5, 10, 9,12,11, 9,10, &
5, 6, 4, 6, 5, 4, 2, 2, 1, 5, 5, 3, 12,11,10, 9, 9,10, &
4, 5, 6, 3, 5, 5, 4, 6, 5, 1, 2, 2, 10, 9, 9,10,12,11, &
5, 3, 5, 5, 4, 6, 6, 4, 5, 2, 1, 2, 10, 9,11,12,10, 9, &
6, 5, 4, 5, 6, 4, 5, 5, 3, 2, 2, 1, 12,11, 9,10,10, 9, &
1, 2, 2, 4, 7, 5, 3, 5, 5, 4, 6, 7, 10,11,10,11,12,13, & ! -----> acting (forest)
2, 1, 2, 7, 4, 5, 6, 4, 7, 5, 3, 5, 10,11,12,13,10,11, & ! |
2, 2, 1, 5, 5, 3, 6, 7, 4, 7, 6, 4, 12,13,10,11,10,11, & ! |
4, 7, 6, 1, 2, 2, 4, 6, 7, 3, 5, 5, 10,11,11,10,13,12, & ! v
7, 4, 6, 2, 1, 2, 5, 3, 5, 6, 4, 7, 10,11,13,12,11,10, & ! reacting (primary)
5, 5, 3, 2, 2, 1, 7, 6, 4, 6, 7, 4, 12,13,11,10,11,10, &
3, 5, 5, 4, 6, 7, 1, 2, 2, 4, 7, 6, 11,10,11,10,12,13, &
6, 4, 7, 5, 3, 5, 2, 1, 2, 7, 4, 6, 11,10,13,12,10,11, &
6, 7, 4, 7, 6, 4, 2, 2, 1, 5, 5, 3, 13,12,11,10,10,11, &
4, 6, 7, 3, 5, 5, 4, 7, 6, 1, 2, 2, 11,10,10,11,13,12, &
5, 3, 5, 6, 4, 7, 7, 4, 6, 2, 1, 2, 11,10,12,13,11,10, &
7, 6, 4, 6, 7, 4, 5, 5, 3, 2, 2, 1, 13,12,10,11,11,10, &
9, 9,11, 9, 9,11,10,10,12,10,10,12, 1, 7, 8, 8, 8, 8, &
10,10,12,10,10,12, 9, 9,11, 9, 9,11, 7, 1, 8, 8, 8, 8, &
9,11, 9,10,12,10,10,12,10, 9,11, 9, 8, 8, 1, 7, 8, 8, &
10,12,10, 9,11, 9, 9,11, 9,10,12,10, 8, 8, 7, 1, 8, 8, &
11, 9, 9,12,10,10,11, 9, 9,12,10,10, 8, 8, 8, 8, 1, 7, &
12,10,10,11, 9, 9,12,10,10,11, 9, 9, 8, 8, 8, 8, 7, 1 &
],shape(FCC_INTERACTIONSLIPSLIP)) !< Slip--slip interaction types for fcc
!< 1: self interaction
!< 2: coplanar interaction
!< 3: collinear interaction
!< 4: Hirth locks
!< 5: glissile junctions
!< 6: Lomer locks
!< 7: crossing (similar to Hirth locks in <110>{111} for two {110} planes)
!< 8: similar to Lomer locks in <110>{111} for two {110} planes
!< 9: similar to Lomer locks in <110>{111} btw one {110} and one {111} plane
!<10: similar to glissile junctions in <110>{111} btw one {110} and one {111} plane
!<11: crossing btw one {110} and one {111} plane
!<12: collinear btw one {110} and one {111} plane
10,10,12,10,10,12,11,11,13,11,11,13, 1, 8, 9, 9, 9, 9, &
11,11,13,11,11,13,10,10,12,10,10,12, 8, 1, 9, 9, 9, 9, &
10,12,10,11,13,11,11,13,11,10,12,10, 9, 9, 1, 8, 9, 9, &
11,13,11,10,12,10,10,12,10,11,13,11, 9, 9, 8, 1, 9, 9, &
12,10,10,13,11,11,12,10,10,13,11,11, 9, 9, 9, 9, 1, 8, &
13,11,11,12,10,10,13,11,11,12,10,10, 9, 9, 9, 9, 8, 1 &
],shape(FCC_INTERACTIONSLIPSLIP)) !< Slip-slip interaction types for fcc / Madec 2017 (https://doi.org/10.1016/j.actamat.2016.12.040)
!< 1: self interaction --> alpha 0
!< 2: coplanar interaction --> alpha copla
!< 3: collinear interaction --> alpha coli
!< 4: Hirth locks --> alpha 1
!< 5: glissile junctions I --> alpha 2
!< 6: glissile junctions II --> alpha 2*
!< 7: Lomer locks --> alpha 3
!< 8: crossing (similar to Hirth locks in <110>{111} for two {110} planes)
!< 9: similar to Lomer locks in <110>{111} for two {110} planes
!<10: similar to Lomer locks in <110>{111} btw one {110} and one {111} plane
!<11: similar to glissile junctions in <110>{111} btw one {110} and one {111} plane
!<12: crossing btw one {110} and one {111} plane
!<13: collinear btw one {110} and one {111} plane
integer, dimension(BCC_NSLIP,BCC_NSLIP), parameter :: &
BCC_INTERACTIONSLIPSLIP = reshape( [&
1,2,6,6,5,4,4,3,4,3,5,4, 6,6,4,3,3,4,6,6,4,3,6,6, & ! -----> acting
2,1,6,6,4,3,5,4,5,4,4,3, 6,6,3,4,4,3,6,6,3,4,6,6, & ! |
6,6,1,2,4,5,3,4,4,5,3,4, 4,3,6,6,6,6,3,4,6,6,4,3, & ! |
6,6,2,1,3,4,4,5,3,4,4,5, 3,4,6,6,6,6,4,3,6,6,3,4, & ! v
5,4,4,3,1,2,6,6,3,4,5,4, 3,6,4,6,6,4,6,3,4,6,3,6, & ! reacting
4,3,5,4,2,1,6,6,4,5,4,3, 4,6,3,6,6,3,6,4,3,6,4,6, &
4,5,3,4,6,6,1,2,5,4,3,4, 6,3,6,4,4,6,3,6,6,4,6,3, &
3,4,4,5,6,6,2,1,4,3,4,5, 6,4,6,3,3,6,4,6,6,3,6,4, &
4,5,4,3,3,4,5,4,1,2,6,6, 3,6,6,4,4,6,6,3,6,4,3,6, &
3,4,5,4,4,5,4,3,2,1,6,6, 4,6,6,3,3,6,6,4,6,3,4,6, &
5,4,3,4,5,4,3,4,6,6,1,2, 6,3,4,6,6,4,3,6,4,6,6,3, &
4,3,4,5,4,3,4,5,6,6,2,1, 6,4,3,6,6,3,4,6,3,6,6,4, &
!
6,6,4,3,3,4,6,6,3,4,6,6, 1,5,6,6,5,6,6,3,5,6,3,6, &
6,6,3,4,6,6,3,4,6,6,3,4, 5,1,6,6,6,5,3,6,6,5,6,3, &
4,3,6,6,4,3,6,6,6,6,4,3, 6,6,1,5,6,3,5,6,3,6,5,6, &
3,4,6,6,6,6,4,3,4,3,6,6, 6,6,5,1,3,6,6,5,6,3,6,5, &
3,4,6,6,6,6,4,3,4,3,6,6, 5,6,6,3,1,6,5,6,5,3,6,6, &
4,3,6,6,4,3,6,6,6,6,4,3, 6,5,3,6,6,1,6,5,3,5,6,6, &
6,6,3,4,6,6,3,4,6,6,3,4, 6,3,5,6,5,6,1,6,6,6,5,3, &
6,6,4,3,3,4,6,6,3,4,6,6, 3,6,6,5,6,5,6,1,6,6,3,5, &
4,3,6,6,4,3,6,6,6,6,4,3, 5,6,3,6,5,3,6,6,1,6,6,5, &
3,4,6,6,6,6,4,3,4,3,6,6, 6,5,6,3,3,5,6,6,6,1,5,6, &
6,6,4,3,3,4,6,6,3,4,6,6, 3,6,5,6,6,6,5,3,6,5,1,6, &
6,6,3,4,6,6,3,4,6,6,3,4, 6,3,6,5,6,6,3,5,5,6,6,1 &
],shape(BCC_INTERACTIONSLIPSLIP)) !< Slip--slip interaction types for bcc from Queyreau et al. Int J Plast 25 (2009) 361377
!< 1: self interaction
!< 2: coplanar interaction
!< 3: collinear interaction
!< 4: mixed-asymmetrical junction
!< 5: mixed-symmetrical junction
!< 6: edge junction
1, 3, 6, 6, 7, 5, 4, 2, 4, 2, 7, 5, 18, 18, 11, 8, 9, 13, 17, 14, 13, 9, 17, 14, 28, 25, 28, 28, 25, 28, 28, 28, 28, 25, 28, 28, 25, 28, 28, 28, 28, 28, 28, 25, 28, 28, 28, 25, &! -----> acting (forest)
3, 1, 6, 6, 4, 2, 7, 5, 7, 5, 4, 2, 18, 18, 8, 11, 13, 9, 14, 17, 9, 13, 14, 17, 25, 28, 28, 28, 28, 25, 28, 28, 25, 28, 28, 28, 28, 25, 28, 28, 28, 28, 25, 28, 28, 28, 25, 28, &! |
6, 6, 1, 3, 5, 7, 2, 4, 5, 7, 2, 4, 11, 8, 18, 18, 17, 14, 9, 13, 17, 14, 13, 9, 28, 28, 28, 25, 28, 28, 25, 28, 28, 28, 28, 25, 28, 28, 25, 28, 28, 25, 28, 28, 28, 25, 28, 28, &! |
6, 6, 3, 1, 2, 4, 5, 7, 2, 4, 5, 7, 8, 11, 18, 18, 14, 17, 13, 9, 14, 17, 9, 13, 28, 28, 25, 28, 28, 28, 28, 25, 28, 28, 25, 28, 28, 28, 28, 25, 25, 28, 28, 28, 25, 28, 28, 28, &! v
7, 5, 4, 2, 1, 3, 6, 6, 2, 4, 7, 5, 9, 17, 13, 14, 18, 11, 18, 8, 13, 17, 9, 14, 28, 28, 25, 28, 28, 28, 28, 25, 28, 28, 25, 28, 28, 28, 28, 25, 25, 28, 28, 28, 25, 28, 28, 28, &! reacting (primary)
4, 2, 7, 5, 3, 1, 6, 6, 5, 7, 4, 2, 13, 14, 9, 17, 18, 8, 18, 11, 9, 14, 13, 17, 25, 28, 28, 28, 28, 25, 28, 28, 25, 28, 28, 28, 28, 25, 28, 28, 28, 28, 25, 28, 28, 28, 25, 28, &
5, 7, 2, 4, 6, 6, 1, 3, 7, 5, 2, 4, 17, 9, 14, 13, 11, 18, 8, 18, 17, 13, 14, 9, 28, 28, 28, 25, 28, 28, 25, 28, 28, 28, 28, 25, 28, 28, 25, 28, 28, 25, 28, 28, 28, 25, 28, 28, &
2, 4, 5, 7, 6, 6, 3, 1, 4, 2, 5, 7, 14, 13, 17, 9, 8, 18, 11, 18, 14, 9, 17, 13, 28, 25, 28, 28, 25, 28, 28, 28, 28, 25, 28, 28, 25, 28, 28, 28, 28, 28, 28, 25, 28, 28, 28, 25, &
5, 7, 4, 2, 2, 4, 7, 5, 1, 3, 6, 6, 9, 17, 14, 13, 13, 17, 14, 9, 18, 11, 8, 18, 28, 28, 25, 28, 28, 28, 28, 25, 28, 28, 25, 28, 28, 28, 28, 25, 25, 28, 28, 28, 25, 28, 28, 28, &
2, 4, 7, 5, 5, 7, 4, 2, 3, 1, 6, 6, 13, 14, 17, 9, 9, 14, 17, 13, 18, 8, 11, 18, 28, 25, 28, 28, 25, 28, 28, 28, 28, 25, 28, 28, 25, 28, 28, 28, 28, 28, 28, 25, 28, 28, 28, 25, &
7, 5, 2, 4, 7, 5, 2, 4, 6, 6, 1, 3, 17, 9, 13, 14, 17, 13, 9, 14, 11, 18, 18, 8, 28, 28, 28, 25, 28, 28, 25, 28, 28, 28, 28, 25, 28, 28, 25, 28, 28, 25, 28, 28, 28, 25, 28, 28, &
4, 2, 5, 7, 4, 2, 5, 7, 6, 6, 3, 1, 14, 13, 9, 17, 14, 9, 13, 17, 8, 18, 18, 11, 25, 28, 28, 28, 28, 25, 28, 28, 25, 28, 28, 28, 28, 25, 28, 28, 28, 28, 25, 28, 28, 28, 25, 28, &
!
19, 19, 10, 8, 9, 12, 16, 15, 9, 12, 16, 15, 1, 20, 24, 24, 23, 22, 21, 2, 23, 22, 2, 21, 28, 28, 26, 28, 28, 28, 28, 26, 28, 28, 26, 28, 28, 28, 28, 26, 26, 28, 28, 28, 26, 28, 28, 28, &
19, 19, 8, 10, 16, 15, 9, 12, 16, 15, 9, 12, 20, 1, 24, 24, 22, 23, 2, 21, 22, 23, 21, 2, 28, 28, 28, 26, 28, 28, 26, 28, 28, 28, 28, 26, 28, 28, 26, 28, 28, 26, 28, 28, 28, 26, 28, 28, &
10, 8, 19, 19, 12, 9, 15, 16, 15, 16, 12, 9, 24, 24, 1, 20, 21, 2, 23, 22, 2, 21, 23, 22, 26, 28, 28, 28, 28, 26, 28, 28, 26, 28, 28, 28, 28, 26, 28, 28, 28, 28, 26, 28, 28, 28, 26, 28, &
8, 10, 19, 19, 15, 16, 12, 9, 12, 9, 15, 16, 24, 24, 20, 1, 2, 21, 22, 23, 21, 2, 22, 23, 28, 26, 28, 28, 26, 28, 28, 28, 28, 26, 28, 28, 26, 28, 28, 28, 28, 28, 28, 26, 28, 28, 28, 26, &
9, 12, 16, 15, 19, 19, 10, 8, 12, 9, 16, 15, 23, 21, 22, 2, 1, 24, 20, 24, 23, 2, 22, 21, 28, 26, 28, 28, 26, 28, 28, 28, 28, 26, 28, 28, 26, 28, 28, 28, 28, 28, 28, 26, 28, 28, 28, 26, &
12, 9, 15, 16, 10, 8, 19, 19, 16, 15, 12, 9, 21, 23, 2, 21, 24, 1, 24, 20, 2, 23, 21, 22, 26, 28, 28, 28, 28, 26, 28, 28, 26, 28, 28, 28, 28, 26, 28, 28, 28, 28, 26, 28, 28, 28, 26, 28, &
16, 15, 9, 12, 19, 19, 8, 10, 15, 16, 9, 12, 22, 2, 23, 22, 20, 24, 1, 24, 22, 21, 23, 2, 28, 28, 28, 26, 28, 28, 26, 28, 28, 28, 28, 26, 28, 28, 26, 28, 28, 26, 28, 28, 28, 26, 28, 28, &
15, 16, 12, 9, 8, 10, 19, 19, 9, 12, 15, 16, 2, 22, 21, 23, 24, 20, 24, 1, 21, 22, 2, 23, 28, 28, 26, 28, 28, 28, 28, 26, 28, 28, 26, 28, 28, 28, 28, 26, 26, 28, 28, 28, 26, 28, 28, 28, &
12, 9, 16, 15, 12, 9, 16, 15, 19, 19, 10, 8, 23, 21, 2, 22, 23, 2, 21, 22, 1, 24, 24, 20, 26, 28, 28, 28, 28, 26, 28, 28, 26, 28, 28, 28, 28, 26, 28, 28, 28, 28, 26, 28, 28, 28, 26, 28, &
9, 12, 15, 16, 16, 15, 12, 9, 10, 8, 19, 19, 21, 23, 22, 2, 2, 23, 22, 21, 24, 1, 20, 24, 28, 26, 28, 28, 26, 28, 28, 28, 28, 26, 28, 28, 26, 28, 28, 28, 28, 28, 28, 26, 28, 28, 28, 26, &
16, 15, 12, 9, 9, 12, 15, 16, 8, 10, 19, 19, 2, 22, 23, 21, 21, 22, 23, 2, 24, 20, 1, 24, 28, 28, 26, 28, 28, 28, 28, 26, 28, 28, 26, 28, 28, 28, 28, 26, 26, 28, 28, 28, 26, 28, 28, 28, &
15, 16, 9, 12, 15, 16, 9, 12, 19, 19, 8, 10, 22, 2, 21, 23, 22, 21, 2, 23, 20, 24, 24, 1, 28, 28, 28, 26, 28, 28, 26, 28, 28, 28, 28, 26, 28, 28, 26, 28, 28, 26, 28, 28, 28, 26, 28, 28, &
!
28, 25, 28, 28, 28, 25, 28, 28, 28, 28, 28, 25, 28, 28, 26, 28, 28, 26, 28, 28, 26, 28, 28, 28, 1, 28, 28, 28, 28, 27, 28, 28, 27, 28, 28, 28, 28, 27, 28, 28, 28, 28, 27, 28, 28, 28, 27, 28, &
25, 28, 28, 28, 28, 28, 28, 25, 28, 25, 28, 28, 28, 28, 28, 26, 26, 28, 28, 28, 28, 26, 28, 28, 28, 1, 28, 28, 27, 28, 28, 28, 28, 27, 28, 28, 27, 28, 28, 28, 28, 28, 28, 27, 28, 28, 28, 27, &
28, 28, 28, 25, 25, 28, 28, 28, 25, 28, 28, 28, 26, 28, 28, 28, 28, 28, 28, 26, 28, 28, 26, 28, 28, 28, 1, 28, 28, 28, 28, 27, 28, 28, 27, 28, 28, 28, 28, 27, 27, 28, 28, 28, 27, 28, 28, 28, &
28, 28, 25, 28, 28, 28, 25, 28, 28, 28, 25, 28, 28, 26, 28, 28, 28, 28, 26, 28, 28, 28, 28, 26, 28, 28, 28, 1, 28, 28, 27, 28, 28, 28, 28, 27, 28, 28, 27, 28, 28, 27, 28, 28, 28, 27, 28, 28, &
25, 28, 28, 28, 28, 28, 28, 25, 28, 25, 28, 28, 28, 28, 28, 26, 26, 28, 28, 28, 28, 26, 28, 28, 28, 27, 28, 28, 1, 28, 28, 28, 28, 27, 28, 28, 27, 28, 28, 28, 28, 28, 28, 27, 28, 28, 28, 27, &
28, 25, 28, 28, 28, 25, 28, 28, 28, 28, 28, 25, 28, 28, 26, 28, 28, 26, 28, 28, 26, 28, 28, 28, 27, 28, 28, 28, 28, 1, 28, 28, 27, 28, 28, 28, 28, 27, 28, 28, 28, 28, 27, 28, 28, 28, 27, 28, &
28, 28, 25, 28, 28, 28, 25, 28, 28, 28, 25, 28, 28, 26, 28, 28, 28, 28, 26, 28, 28, 28, 28, 26, 28, 28, 28, 27, 28, 28, 1, 28, 28, 28, 28, 27, 28, 28, 27, 28, 28, 27, 28, 28, 28, 27, 28, 28, &
28, 28, 28, 25, 25, 28, 28, 28, 25, 28, 28, 28, 26, 28, 28, 28, 28, 28, 28, 26, 28, 28, 26, 28, 28, 28, 27, 28, 28, 28, 28, 1, 28, 28, 27, 28, 28, 28, 28, 27, 27, 28, 28, 28, 27, 28, 28, 28, &
28, 25, 28, 28, 28, 25, 28, 28, 28, 28, 28, 25, 28, 28, 26, 28, 28, 26, 28, 28, 26, 28, 28, 28, 27, 28, 28, 28, 28, 27, 28, 28, 1, 28, 28, 28, 28, 27, 28, 28, 28, 28, 27, 28, 28, 28, 27, 28, &
25, 28, 28, 28, 28, 28, 28, 25, 28, 25, 28, 28, 28, 28, 28, 26, 26, 28, 28, 28, 28, 26, 28, 28, 28, 27, 28, 28, 27, 28, 28, 28, 28, 1, 28, 28, 27, 28, 28, 28, 28, 28, 28, 27, 28, 28, 28, 27, &
28, 28, 28, 25, 25, 28, 28, 28, 25, 28, 28, 28, 26, 28, 28, 28, 28, 28, 28, 26, 28, 28, 26, 28, 28, 28, 27, 28, 28, 28, 28, 27, 28, 28, 1, 28, 28, 28, 28, 27, 27, 28, 28, 28, 27, 28, 28, 28, &
28, 28, 25, 28, 28, 28, 25, 28, 28, 28, 25, 28, 28, 26, 28, 28, 28, 28, 26, 28, 28, 28, 28, 26, 28, 28, 28, 27, 28, 28, 27, 28, 28, 28, 28, 1, 28, 28, 27, 28, 28, 27, 28, 28, 28, 27, 28, 28, &
25, 28, 28, 28, 28, 28, 28, 25, 28, 25, 28, 28, 28, 28, 28, 26, 26, 28, 28, 28, 28, 26, 28, 28, 28, 27, 28, 28, 27, 28, 28, 28, 28, 27, 28, 28, 1, 28, 28, 28, 28, 28, 28, 27, 28, 28, 28, 27, &
28, 25, 28, 28, 28, 25, 28, 28, 28, 28, 28, 25, 28, 28, 26, 28, 28, 26, 28, 28, 26, 28, 28, 28, 27, 28, 28, 28, 28, 27, 28, 28, 27, 28, 28, 28, 28, 1, 28, 28, 28, 28, 27, 28, 28, 28, 27, 28, &
28, 28, 25, 28, 28, 28, 25, 28, 28, 28, 25, 28, 28, 26, 28, 28, 28, 28, 26, 28, 28, 28, 28, 26, 28, 28, 28, 27, 28, 28, 27, 28, 28, 28, 28, 27, 28, 28, 1, 28, 28, 27, 28, 28, 28, 27, 28, 28, &
28, 28, 28, 25, 25, 28, 28, 28, 25, 28, 28, 28, 26, 28, 28, 28, 28, 28, 28, 26, 28, 28, 26, 28, 28, 28, 27, 28, 28, 28, 28, 27, 28, 28, 27, 28, 28, 28, 28, 1, 27, 28, 28, 28, 27, 28, 28, 28, &
28, 28, 28, 25, 25, 28, 28, 28, 25, 28, 28, 28, 26, 28, 28, 28, 28, 28, 28, 26, 28, 28, 26, 28, 28, 28, 27, 28, 28, 28, 28, 27, 28, 28, 27, 28, 28, 28, 28, 27, 1, 28, 28, 28, 27, 28, 28, 28, &
28, 28, 25, 28, 28, 28, 25, 28, 28, 28, 25, 28, 28, 26, 28, 28, 28, 28, 26, 28, 28, 28, 28, 26, 28, 28, 28, 27, 28, 28, 27, 28, 28, 28, 28, 27, 28, 28, 27, 28, 28, 1, 28, 28, 28, 27, 28, 28, &
28, 25, 28, 28, 28, 25, 28, 28, 28, 28, 28, 25, 28, 28, 26, 28, 28, 26, 28, 28, 26, 28, 28, 28, 27, 28, 28, 28, 28, 27, 28, 28, 27, 28, 28, 28, 28, 27, 28, 28, 28, 28, 1, 28, 28, 28, 27, 28, &
25, 28, 28, 28, 28, 28, 28, 25, 28, 25, 28, 28, 28, 28, 28, 26, 26, 28, 28, 28, 28, 26, 28, 28, 28, 27, 28, 28, 27, 28, 28, 28, 28, 27, 28, 28, 27, 28, 28, 28, 28, 28, 28, 1, 28, 28, 28, 27, &
28, 28, 28, 25, 25, 28, 28, 28, 25, 28, 28, 28, 26, 28, 28, 28, 28, 28, 28, 26, 28, 28, 26, 28, 28, 28, 27, 28, 28, 28, 28, 27, 28, 28, 27, 28, 28, 28, 28, 27, 27, 28, 28, 28, 1, 28, 28, 28, &
28, 28, 25, 28, 28, 28, 25, 28, 28, 28, 25, 28, 28, 26, 28, 28, 28, 28, 26, 28, 28, 28, 28, 26, 28, 28, 28, 27, 28, 28, 27, 28, 28, 28, 28, 27, 28, 28, 27, 28, 28, 27, 28, 28, 28, 1, 28, 28, &
28, 25, 28, 28, 28, 25, 28, 28, 28, 28, 28, 25, 28, 28, 26, 28, 28, 26, 28, 28, 26, 28, 28, 28, 27, 28, 28, 28, 28, 27, 28, 28, 27, 28, 28, 28, 28, 27, 28, 28, 28, 28, 27, 28, 28, 28, 1, 28, &
25, 28, 28, 28, 28, 28, 28, 25, 28, 25, 28, 28, 28, 28, 28, 26, 26, 28, 28, 28, 28, 26, 28, 28, 28, 27, 28, 28, 27, 28, 28, 28, 28, 27, 28, 28, 27, 28, 28, 28, 28, 28, 28, 27, 28, 28, 28, 1 &
],shape(BCC_INTERACTIONSLIPSLIP)) !< Slip-slip interaction types for bcc / Madec 2017 (https://doi.org/10.1016/j.actamat.2016.12.040)
!< 1: self interaction --> alpha 0
!< 2: collinear interaction --> alpha 1
!< 3: coplanar interaction --> alpha 2
!< 4-7: other coefficients
!< 8: {110}-{112}, collinear and perpendicular planes --> alpha 6
!< 9: {110}-{112}, just collinear --> alpha 7
!< 10-24: other coefficients
!< 25: {110}-{123}, just collinear
!< 26: {112}-{123}, just collinear
!< 27: {123}-{123}, just collinear
!< 28: other interaction
integer, dimension(HEX_NSLIP,HEX_NSLIP), parameter :: &
HEX_INTERACTIONSLIPSLIP = reshape( [&
1, 2, 2, 3, 3, 3, 7, 7, 7, 13,13,13,13,13,13, 21,21,21,21,21,21,21,21,21,21,21,21, 31,31,31,31,31,31, & ! -----> acting
1, 2, 2, 3, 3, 3, 7, 7, 7, 13,13,13,13,13,13, 21,21,21,21,21,21,21,21,21,21,21,21, 31,31,31,31,31,31, & ! -----> acting (forest)
2, 1, 2, 3, 3, 3, 7, 7, 7, 13,13,13,13,13,13, 21,21,21,21,21,21,21,21,21,21,21,21, 31,31,31,31,31,31, & ! |
2, 2, 1, 3, 3, 3, 7, 7, 7, 13,13,13,13,13,13, 21,21,21,21,21,21,21,21,21,21,21,21, 31,31,31,31,31,31, & ! |
! ! v
6, 6, 6, 4, 5, 5, 8, 8, 8, 14,14,14,14,14,14, 22,22,22,22,22,22,22,22,22,22,22,22, 32,32,32,32,32,32, & ! reacting
6, 6, 6, 4, 5, 5, 8, 8, 8, 14,14,14,14,14,14, 22,22,22,22,22,22,22,22,22,22,22,22, 32,32,32,32,32,32, & ! reacting (primary)
6, 6, 6, 5, 4, 5, 8, 8, 8, 14,14,14,14,14,14, 22,22,22,22,22,22,22,22,22,22,22,22, 32,32,32,32,32,32, &
6, 6, 6, 5, 5, 4, 8, 8, 8, 14,14,14,14,14,14, 22,22,22,22,22,22,22,22,22,22,22,22, 32,32,32,32,32,32, &
!
@ -867,7 +898,7 @@ function lattice_interaction_SlipBySlip(Nslip,interactionValues,structure) resul
42,42,42, 41,41,41, 40,40,40, 39,39,39,39,39,39, 38,38,38,38,38,38,38,38,38,38,38,38, 37,37,37,36,37,37, &
42,42,42, 41,41,41, 40,40,40, 39,39,39,39,39,39, 38,38,38,38,38,38,38,38,38,38,38,38, 37,37,37,37,36,37, &
42,42,42, 41,41,41, 40,40,40, 39,39,39,39,39,39, 38,38,38,38,38,38,38,38,38,38,38,38, 37,37,37,37,37,36 &
],shape(HEX_INTERACTIONSLIPSLIP)) !< Slip--slip interaction types for hex (onion peel naming scheme)
],shape(HEX_INTERACTIONSLIPSLIP)) !< Slip-slip interaction types for hex (onion peel naming scheme)
integer, dimension(BCT_NSLIP,BCT_NSLIP), parameter :: &
BCT_INTERACTIONSLIPSLIP = reshape( [&
@ -1167,11 +1198,38 @@ function lattice_interaction_SlipByTwin(Nslip,Ntwin,interactionValues,structure)
3,3,2,3,3,2,3,3,1,3,3,3, &
3,3,3,2,2,3,3,3,3,1,3,3, &
2,3,3,3,3,3,3,2,3,3,1,3, &
3,2,3,3,3,3,2,3,3,3,3,1 &
3,2,3,3,3,3,2,3,3,3,3,1, &
!
4,4,4,4,4,4,4,4,4,4,4,4, &
4,4,4,4,4,4,4,4,4,4,4,4, &
4,4,4,4,4,4,4,4,4,4,4,4, &
4,4,4,4,4,4,4,4,4,4,4,4, &
4,4,4,4,4,4,4,4,4,4,4,4, &
4,4,4,4,4,4,4,4,4,4,4,4, &
4,4,4,4,4,4,4,4,4,4,4,4, &
4,4,4,4,4,4,4,4,4,4,4,4, &
4,4,4,4,4,4,4,4,4,4,4,4, &
4,4,4,4,4,4,4,4,4,4,4,4, &
4,4,4,4,4,4,4,4,4,4,4,4, &
4,4,4,4,4,4,4,4,4,4,4,4, &
4,4,4,4,4,4,4,4,4,4,4,4, &
4,4,4,4,4,4,4,4,4,4,4,4, &
4,4,4,4,4,4,4,4,4,4,4,4, &
4,4,4,4,4,4,4,4,4,4,4,4, &
4,4,4,4,4,4,4,4,4,4,4,4, &
4,4,4,4,4,4,4,4,4,4,4,4, &
4,4,4,4,4,4,4,4,4,4,4,4, &
4,4,4,4,4,4,4,4,4,4,4,4, &
4,4,4,4,4,4,4,4,4,4,4,4, &
4,4,4,4,4,4,4,4,4,4,4,4, &
4,4,4,4,4,4,4,4,4,4,4,4, &
4,4,4,4,4,4,4,4,4,4,4,4 &
],shape(BCC_INTERACTIONSLIPTWIN)) !< Slip-twin interaction types for bcc
!< 1: coplanar interaction
!< 2: screw trace between slip system and twin habit plane (easy cross slip)
!< 3: other interaction
!< 4: other interaction with slip family {123}
integer, dimension(HEX_NTWIN,HEX_NSLIP), parameter :: &
HEX_INTERACTIONSLIPTWIN = reshape( [&
1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, & ! ----> twin (acting)
@ -1574,7 +1632,7 @@ end function lattice_slip_normal
!--------------------------------------------------------------------------------------------------
!> @brief Transverse direction of slip systems ( || t = b x n)
!> @brief Transverse direction of slip systems (|| t = b x n)
!--------------------------------------------------------------------------------------------------
function lattice_slip_transverse(Nslip,structure,cOverA) result(t)

View File

@ -58,10 +58,6 @@ module phase
grain
end type tDebugOptions
integer, dimension(:), allocatable, public :: & !< ToDo: should be protected (bug in Intel compiler)
phase_elasticityInstance, &
phase_NstiffnessDegradations
logical, dimension(:), allocatable, public :: & ! ToDo: should be protected (bug in Intel Compiler)
phase_localPlasticity !< flags phases with local constitutive law
@ -71,10 +67,6 @@ module phase
damageState
integer, public, protected :: &
phase_plasticity_maxSizeDotState, &
phase_source_maxSizeDotState
interface
! == cleaned:begin =================================================================================
@ -239,11 +231,16 @@ module phase
logical :: converged_
end function crystallite_stress
!ToDo: Try to merge the all stiffness functions
module function phase_homogenizedC(ph,en) result(C)
integer, intent(in) :: ph, en
real(pReal), dimension(6,6) :: C
end function phase_homogenizedC
module function phase_damage_C(C_homogenized,ph,en) result(C)
real(pReal), dimension(3,3,3,3), intent(in) :: C_homogenized
integer, intent(in) :: ph,en
real(pReal), dimension(3,3,3,3) :: C
end function phase_damage_C
module function phase_f_phi(phi,co,ce) result(f)
integer, intent(in) :: ce,co
@ -298,7 +295,6 @@ module phase
end interface
type(tDebugOptions) :: debugConstitutive
#if __INTEL_COMPILER >= 1900
public :: &
@ -377,19 +373,6 @@ subroutine phase_init
call damage_init
call thermal_init(phases)
phase_source_maxSizeDotState = 0
PhaseLoop2:do ph = 1,phases%length
!--------------------------------------------------------------------------------------------------
! partition and initialize state
plasticState(ph)%state = plasticState(ph)%state0
if(damageState(ph)%sizeState > 0) &
damageState(ph)%state = damageState(ph)%state0
enddo PhaseLoop2
phase_source_maxSizeDotState = maxval(damageState%sizeDotState)
phase_plasticity_maxSizeDotState = maxval(plasticState%sizeDotState)
end subroutine phase_init
@ -545,8 +528,7 @@ subroutine crystallite_init()
phases => config_material%get('phase')
do ph = 1, phases%length
if (damageState(ph)%sizeState > 0) &
allocate(damageState(ph)%subState0,source=damageState(ph)%state0) ! ToDo: hack
if (damageState(ph)%sizeState > 0) allocate(damageState(ph)%subState0,source=damageState(ph)%state0) ! ToDo: hack
enddo
print'(a42,1x,i10)', ' # of elements: ', eMax
@ -560,7 +542,7 @@ subroutine crystallite_init()
do ip = 1, size(material_phaseMemberAt,2)
do co = 1,homogenization_Nconstituents(material_homogenizationAt(el))
call crystallite_orientations(co,ip,el)
call plastic_dependentState(co,ip,el) ! update dependent state variables to be consistent with basic states
call plastic_dependentState(co,ip,el) ! update dependent state variables to be consistent with basic states
enddo
enddo
enddo
@ -576,9 +558,9 @@ end subroutine crystallite_init
subroutine crystallite_orientations(co,ip,el)
integer, intent(in) :: &
co, & !< counter in integration point component loop
ip, & !< counter in integration point loop
el !< counter in element loop
co, & !< counter in integration point component loop
ip, & !< counter in integration point loop
el !< counter in element loop
call crystallite_orientation(co,ip,el)%fromMatrix(transpose(math_rotationalPart(&

View File

@ -15,13 +15,15 @@ submodule(phase) damage
DAMAGE_ANISOBRITTLE_ID
end enum
integer :: phase_damage_maxSizeDotState
type :: tDataContainer
real(pReal), dimension(:), allocatable :: phi, d_phi_d_dot_phi
end type tDataContainer
integer(kind(DAMAGE_UNDEFINED_ID)), dimension(:), allocatable :: &
phase_source !< active sources mechanisms of each phase
phase_damage !< active sources mechanisms of each phase
type(tDataContainer), dimension(:), allocatable :: current
@ -126,17 +128,38 @@ module subroutine damage_init
enddo
allocate(phase_source(phases%length), source = DAMAGE_UNDEFINED_ID)
allocate(phase_damage(phases%length), source = DAMAGE_UNDEFINED_ID)
if (damage_active) then
where(isobrittle_init() ) phase_source = DAMAGE_ISOBRITTLE_ID
where(isoductile_init() ) phase_source = DAMAGE_ISODUCTILE_ID
where(anisobrittle_init()) phase_source = DAMAGE_ANISOBRITTLE_ID
where(isobrittle_init() ) phase_damage = DAMAGE_ISOBRITTLE_ID
where(isoductile_init() ) phase_damage = DAMAGE_ISODUCTILE_ID
where(anisobrittle_init()) phase_damage = DAMAGE_ANISOBRITTLE_ID
endif
phase_damage_maxSizeDotState = maxval(damageState%sizeDotState)
end subroutine damage_init
!--------------------------------------------------------------------------------------------------
!> @brief returns the degraded/modified elasticity matrix
!--------------------------------------------------------------------------------------------------
module function phase_damage_C(C_homogenized,ph,en) result(C)
real(pReal), dimension(3,3,3,3), intent(in) :: C_homogenized
integer, intent(in) :: ph,en
real(pReal), dimension(3,3,3,3) :: C
damageType: select case (phase_damage(ph))
case (DAMAGE_ISOBRITTLE_ID) damageType
C = C_homogenized * damage_phi(ph,en)**2
case default damageType
C = C_homogenized
end select damageType
end function phase_damage_C
!----------------------------------------------------------------------------------------------
!< @brief returns local part of nonlocal damage driving force
!----------------------------------------------------------------------------------------------
@ -155,7 +178,7 @@ module function phase_f_phi(phi,co,ce) result(f)
ph = material_phaseID(co,ce)
en = material_phaseEntry(co,ce)
select case(phase_source(ph))
select case(phase_damage(ph))
case(DAMAGE_ISOBRITTLE_ID,DAMAGE_ISODUCTILE_ID,DAMAGE_ANISOBRITTLE_ID)
f = 1.0_pReal &
- phi*damageState(ph)%state(1,en)
@ -187,9 +210,9 @@ module function integrateDamageState(dt,co,ip,el) result(broken)
size_so
real(pReal) :: &
zeta
real(pReal), dimension(phase_source_maxSizeDotState) :: &
real(pReal), dimension(phase_damage_maxSizeDotState) :: &
r ! state residuum
real(pReal), dimension(phase_source_maxSizeDotState,2) :: source_dotState
real(pReal), dimension(phase_damage_maxSizeDotState,2) :: source_dotState
logical :: &
converged_
@ -275,10 +298,10 @@ module subroutine damage_results(group,ph)
integer, intent(in) :: ph
if (phase_source(ph) /= DAMAGE_UNDEFINED_ID) &
if (phase_damage(ph) /= DAMAGE_UNDEFINED_ID) &
call results_closeGroup(results_addGroup(group//'damage'))
sourceType: select case (phase_source(ph))
sourceType: select case (phase_damage(ph))
case (DAMAGE_ISOBRITTLE_ID) sourceType
call isobrittle_results(ph,group//'damage/')
@ -309,7 +332,7 @@ function phase_damage_collectDotState(ph,me) result(broken)
if (damageState(ph)%sizeState > 0) then
sourceType: select case (phase_source(ph))
sourceType: select case (phase_damage(ph))
case (DAMAGE_ISODUCTILE_ID) sourceType
call isoductile_dotState(ph,me)
@ -376,7 +399,7 @@ function phase_damage_deltaState(Fe, ph, me) result(broken)
if (damageState(ph)%sizeState == 0) return
sourceType: select case (phase_source(ph))
sourceType: select case (phase_damage(ph))
case (DAMAGE_ISOBRITTLE_ID) sourceType
call isobrittle_deltaState(phase_homogenizedC(ph,me), Fe, ph,me)

View File

@ -66,14 +66,14 @@ module function isobrittle_init() result(mySources)
Nmembers = count(material_phaseID==ph)
call phase_allocateState(damageState(ph),Nmembers,1,1,1)
damageState(ph)%atol = src%get_asFloat('isoBrittle_atol',defaultVal=1.0e-3_pReal)
damageState(ph)%atol = src%get_asFloat('isobrittle_atol',defaultVal=1.0e-3_pReal)
if(any(damageState(ph)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isobrittle_atol'
end associate
!--------------------------------------------------------------------------------------------------
! exit if any parameter is out of range
if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(damage_isoBrittle)')
if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(damage_isobrittle)')
endif
enddo

View File

@ -70,14 +70,14 @@ module function isoductile_init() result(mySources)
Nmembers=count(material_phaseID==ph)
call phase_allocateState(damageState(ph),Nmembers,1,1,0)
damageState(ph)%atol = src%get_asFloat('isoDuctile_atol',defaultVal=1.0e-3_pReal)
damageState(ph)%atol = src%get_asFloat('isoductile_atol',defaultVal=1.0e-3_pReal)
if(any(damageState(ph)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isoductile_atol'
end associate
!--------------------------------------------------------------------------------------------------
! exit if any parameter is out of range
if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(damage_isoDuctile)')
if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(damage_isoductile)')
endif
enddo

View File

@ -5,10 +5,6 @@ submodule(phase) mechanical
enum, bind(c); enumerator :: &
ELASTICITY_UNDEFINED_ID, &
ELASTICITY_HOOKE_ID, &
STIFFNESS_DEGRADATION_UNDEFINED_ID, &
STIFFNESS_DEGRADATION_DAMAGE_ID, &
PLASTICITY_UNDEFINED_ID, &
PLASTICITY_NONE_ID, &
PLASTICITY_ISOTROPIC_ID, &
@ -23,11 +19,6 @@ submodule(phase) mechanical
KINEMATICS_THERMAL_EXPANSION_ID
end enum
integer(kind(ELASTICITY_UNDEFINED_ID)), dimension(:), allocatable :: &
phase_elasticity !< elasticity of each phase
integer(kind(STIFFNESS_DEGRADATION_UNDEFINED_ID)), dimension(:,:), allocatable :: &
phase_stiffnessDegradation !< active stiffness degradation mechanisms of each phase
type(tTensorContainer), dimension(:), allocatable :: &
! current value
phase_mechanical_Fe, &
@ -50,6 +41,7 @@ submodule(phase) mechanical
integer(kind(PLASTICITY_undefined_ID)), dimension(:), allocatable :: &
phase_plasticity !< plasticity of each phase
integer :: phase_plasticity_maxSizeDotState
interface
@ -57,9 +49,27 @@ submodule(phase) mechanical
class(tNode), pointer :: phases
end subroutine eigendeformation_init
module subroutine elastic_init(phases)
class(tNode), pointer :: phases
end subroutine elastic_init
module subroutine plastic_init
end subroutine plastic_init
module subroutine phase_hooke_SandItsTangents(S,dS_dFe,dS_dFi,Fe,Fi,ph,en)
integer, intent(in) :: &
ph, &
en
real(pReal), intent(in), dimension(3,3) :: &
Fe, & !< elastic deformation gradient
Fi !< intermediate deformation gradient
real(pReal), intent(out), dimension(3,3) :: &
S !< 2nd Piola-Kirchhoff stress tensor in lattice configuration
real(pReal), intent(out), dimension(3,3,3,3) :: &
dS_dFe, & !< derivative of 2nd P-K stress with respect to elastic deformation gradient
dS_dFi !< derivative of 2nd P-K stress with respect to intermediate deformation gradient
end subroutine phase_hooke_SandItsTangents
module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,ph,en)
real(pReal), dimension(3,3), intent(out) :: &
Li !< inleastic velocity gradient
@ -73,7 +83,6 @@ submodule(phase) mechanical
end subroutine plastic_isotropic_LiAndItsTangent
module function plastic_dotState(subdt,co,ip,el,ph,en) result(broken)
integer, intent(in) :: &
co, & !< component-ID of integration point
ip, & !< integration point
@ -198,17 +207,11 @@ module subroutine mechanical_init(materials,phases)
constituents, &
constituent, &
phase, &
mech, &
elastic, &
stiffDegradation
mech
print'(/,a)', ' <<<+- phase:mechanical init -+>>>'
!-------------------------------------------------------------------------------------------------
! initialize elasticity (hooke) !ToDO: Maybe move to elastic submodule along with function homogenizedC?
allocate(phase_elasticity(phases%length), source = ELASTICITY_undefined_ID)
allocate(phase_elasticityInstance(phases%length), source = 0)
allocate(phase_NstiffnessDegradations(phases%length),source=0)
allocate(output_constituent(phases%length))
allocate(phase_mechanical_Fe(phases%length))
@ -253,32 +256,8 @@ module subroutine mechanical_init(materials,phases)
#else
output_constituent(ph)%label = mech%get_as1dString('output',defaultVal=emptyStringArray)
#endif
elastic => mech%get('elastic')
if (IO_lc(elastic%get_asString('type')) == 'hooke') then ! accept small letter h for the moment
phase_elasticity(ph) = ELASTICITY_HOOKE_ID
else
call IO_error(200,ext_msg=elastic%get_asString('type'))
endif
stiffDegradation => mech%get('stiffness_degradation',defaultVal=emptyList) ! check for stiffness degradation mechanisms
phase_NstiffnessDegradations(ph) = stiffDegradation%length
enddo
allocate(phase_stiffnessDegradation(maxval(phase_NstiffnessDegradations),phases%length), &
source=STIFFNESS_DEGRADATION_undefined_ID)
if(maxVal(phase_NstiffnessDegradations)/=0) then
do ph = 1, phases%length
phase => phases%get(ph)
mech => phase%get('mechanical')
stiffDegradation => mech%get('stiffness_degradation',defaultVal=emptyList)
do stiffDegradationCtr = 1, stiffDegradation%length
if(stiffDegradation%get_asString(stiffDegradationCtr) == 'damage') &
phase_stiffnessDegradation(stiffDegradationCtr,ph) = STIFFNESS_DEGRADATION_damage_ID
enddo
enddo
endif
do el = 1, size(material_phaseMemberAt,3); do ip = 1, size(material_phaseMemberAt,2)
do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
material => materials%get(discretization_materialAt(el))
@ -306,6 +285,9 @@ module subroutine mechanical_init(materials,phases)
enddo; enddo
! initialize elasticity
call elastic_init(phases)
! initialize plasticity
allocate(plasticState(phases%length))
allocate(phase_plasticity(phases%length),source = PLASTICITY_undefined_ID)
@ -313,9 +295,11 @@ module subroutine mechanical_init(materials,phases)
call plastic_init()
do ph = 1, phases%length
phase_elasticityInstance(ph) = count(phase_elasticity(1:ph) == phase_elasticity(ph))
do ph = 1,phases%length
plasticState(ph)%state0 = plasticState(ph)%state
enddo
phase_plasticity_maxSizeDotState = maxval(plasticState%sizeDotState)
num_crystallite => config_numerics%get('crystallite',defaultVal=emptyDict)
@ -348,51 +332,6 @@ module subroutine mechanical_init(materials,phases)
end subroutine mechanical_init
!--------------------------------------------------------------------------------------------------
!> @brief returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to
!> the elastic and intermediate deformation gradients using Hooke's law
!--------------------------------------------------------------------------------------------------
subroutine phase_hooke_SandItsTangents(S, dS_dFe, dS_dFi, &
Fe, Fi, ph, en)
integer, intent(in) :: &
ph, &
en
real(pReal), intent(in), dimension(3,3) :: &
Fe, & !< elastic deformation gradient
Fi !< intermediate deformation gradient
real(pReal), intent(out), dimension(3,3) :: &
S !< 2nd Piola-Kirchhoff stress tensor in lattice configuration
real(pReal), intent(out), dimension(3,3,3,3) :: &
dS_dFe, & !< derivative of 2nd P-K stress with respect to elastic deformation gradient
dS_dFi !< derivative of 2nd P-K stress with respect to intermediate deformation gradient
real(pReal), dimension(3,3) :: E
real(pReal), dimension(3,3,3,3) :: C
integer :: &
d, & !< counter in degradation loop
i, j
C = math_66toSym3333(phase_homogenizedC(ph,en))
DegradationLoop: do d = 1, phase_NstiffnessDegradations(ph)
degradationType: select case(phase_stiffnessDegradation(d,ph))
case (STIFFNESS_DEGRADATION_damage_ID) degradationType
C = C * damage_phi(ph,en)**2
end select degradationType
enddo DegradationLoop
E = 0.5_pReal*(matmul(transpose(Fe),Fe)-math_I3) !< Green-Lagrange strain in unloaded configuration
S = math_mul3333xx33(C,matmul(matmul(transpose(Fi),E),Fi)) !< 2PK stress in lattice configuration in work conjugate with GL strain pulled back to lattice configuration
do i =1, 3;do j=1,3
dS_dFe(i,j,1:3,1:3) = matmul(Fe,matmul(matmul(Fi,C(i,j,1:3,1:3)),transpose(Fi))) !< dS_ij/dFe_kl = C_ijmn * Fi_lm * Fi_on * Fe_ko
dS_dFi(i,j,1:3,1:3) = 2.0_pReal*matmul(matmul(E,Fi),C(i,j,1:3,1:3)) !< dS_ij/dFi_kl = C_ijln * E_km * Fe_mn
enddo; enddo
end subroutine phase_hooke_SandItsTangents
module subroutine mechanical_results(group,ph)
character(len=*), intent(in) :: group
@ -1082,26 +1021,6 @@ module subroutine mechanical_forward()
end subroutine mechanical_forward
!--------------------------------------------------------------------------------------------------
!> @brief returns the homogenize elasticity matrix
!> ToDo: homogenizedC66 would be more consistent
!--------------------------------------------------------------------------------------------------
module function phase_homogenizedC(ph,en) result(C)
real(pReal), dimension(6,6) :: C
integer, intent(in) :: ph, en
plasticType: select case (phase_plasticity(ph))
case (PLASTICITY_DISLOTWIN_ID) plasticType
C = plastic_dislotwin_homogenizedC(ph,en)
case default plasticType
C = lattice_C66(1:6,1:6,ph)
end select plasticType
end function phase_homogenizedC
!--------------------------------------------------------------------------------------------------
!> @brief calculate stress (P)
!--------------------------------------------------------------------------------------------------

View File

@ -40,7 +40,6 @@ module function thermalexpansion_init(kinematics_length) result(myKinematics)
print'(/,a)', ' <<<+- phase:mechanical:eigen:thermalexpansion init -+>>>'
myKinematics = kinematics_active('thermalexpansion',kinematics_length)
print*, myKinematics
Ninstances = count(myKinematics)
print'(a,i2)', ' # phases: ',Ninstances; flush(IO_STDOUT)
if(Ninstances == 0) return
@ -101,16 +100,16 @@ module subroutine thermalexpansion_LiAndItsTangent(Li, dLi_dTstar, ph,me)
associate(prm => param(kinematics_thermal_expansion_instance(ph)))
Li = dot_T * ( &
prm%A(1:3,1:3,1)*(T - prm%T_ref)**0 & ! constant coefficient
+ prm%A(1:3,1:3,2)*(T - prm%T_ref)**1 & ! linear coefficient
+ prm%A(1:3,1:3,3)*(T - prm%T_ref)**2 & ! quadratic coefficient
prm%A(1:3,1:3,1)*(T - prm%T_ref)**0 & ! constant coefficient
+ prm%A(1:3,1:3,2)*(T - prm%T_ref)**1 & ! linear coefficient
+ prm%A(1:3,1:3,3)*(T - prm%T_ref)**2 & ! quadratic coefficient
) / &
(1.0_pReal &
+ prm%A(1:3,1:3,1)*(T - prm%T_ref)**1 / 1. &
+ prm%A(1:3,1:3,2)*(T - prm%T_ref)**2 / 2. &
+ prm%A(1:3,1:3,3)*(T - prm%T_ref)**3 / 3. &
)
end associate
end associate
dLi_dTstar = 0.0_pReal
end subroutine thermalexpansion_LiAndItsTangent

View File

@ -0,0 +1,102 @@
submodule(phase:mechanical) elastic
enum, bind(c); enumerator :: &
ELASTICITY_UNDEFINED_ID, &
ELASTICITY_HOOKE_ID
end enum
integer(kind(ELASTICITY_UNDEFINED_ID)), dimension(:), allocatable :: &
phase_elasticity !< elasticity of each phase
contains
module subroutine elastic_init(phases)
class(tNode), pointer :: &
phases
integer :: &
ph
class(tNode), pointer :: &
phase, &
mech, &
elastic
print'(/,a)', ' <<<+- phase:mechanical:elastic init -+>>>'
allocate(phase_elasticity(phases%length), source = ELASTICITY_undefined_ID)
do ph = 1, phases%length
phase => phases%get(ph)
mech => phase%get('mechanical')
elastic => mech%get('elastic')
if(IO_lc(elastic%get_asString('type')) == 'hooke') then ! accept small letter h for the moment
phase_elasticity(ph) = ELASTICITY_HOOKE_ID
else
call IO_error(200,ext_msg=elastic%get_asString('type'))
endif
enddo
end subroutine elastic_init
!--------------------------------------------------------------------------------------------------
!> @brief returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to
!> the elastic and intermediate deformation gradients using Hooke's law
!--------------------------------------------------------------------------------------------------
module subroutine phase_hooke_SandItsTangents(S, dS_dFe, dS_dFi, &
Fe, Fi, ph, en)
integer, intent(in) :: &
ph, &
en
real(pReal), intent(in), dimension(3,3) :: &
Fe, & !< elastic deformation gradient
Fi !< intermediate deformation gradient
real(pReal), intent(out), dimension(3,3) :: &
S !< 2nd Piola-Kirchhoff stress tensor in lattice configuration
real(pReal), intent(out), dimension(3,3,3,3) :: &
dS_dFe, & !< derivative of 2nd P-K stress with respect to elastic deformation gradient
dS_dFi !< derivative of 2nd P-K stress with respect to intermediate deformation gradient
real(pReal), dimension(3,3) :: E
real(pReal), dimension(3,3,3,3) :: C
integer :: &
d, & !< counter in degradation loop
i, j
C = math_66toSym3333(phase_homogenizedC(ph,en))
C = phase_damage_C(C,ph,en)
E = 0.5_pReal*(matmul(transpose(Fe),Fe)-math_I3) !< Green-Lagrange strain in unloaded configuration
S = math_mul3333xx33(C,matmul(matmul(transpose(Fi),E),Fi)) !< 2PK stress in lattice configuration in work conjugate with GL strain pulled back to lattice configuration
do i =1, 3;do j=1,3
dS_dFe(i,j,1:3,1:3) = matmul(Fe,matmul(matmul(Fi,C(i,j,1:3,1:3)),transpose(Fi))) !< dS_ij/dFe_kl = C_ijmn * Fi_lm * Fi_on * Fe_ko
dS_dFi(i,j,1:3,1:3) = 2.0_pReal*matmul(matmul(E,Fi),C(i,j,1:3,1:3)) !< dS_ij/dFi_kl = C_ijln * E_km * Fe_mn
enddo; enddo
end subroutine phase_hooke_SandItsTangents
!--------------------------------------------------------------------------------------------------
!> @brief returns the homogenized elasticity matrix
!> ToDo: homogenizedC66 would be more consistent
!--------------------------------------------------------------------------------------------------
module function phase_homogenizedC(ph,en) result(C)
real(pReal), dimension(6,6) :: C
integer, intent(in) :: ph, en
plasticType: select case (phase_plasticity(ph))
case (PLASTICITY_DISLOTWIN_ID) plasticType
C = plastic_dislotwin_homogenizedC(ph,en)
case default plasticType
C = lattice_C66(1:6,1:6,ph)
end select plasticType
end function phase_homogenizedC
end submodule elastic

View File

@ -254,8 +254,6 @@ module function plastic_dislotungsten_init() result(myPlasticity)
allocate(dst%Lambda_sl(prm%sum_N_sl,Nmembers), source=0.0_pReal)
allocate(dst%threshold_stress(prm%sum_N_sl,Nmembers), source=0.0_pReal)
plasticState(ph)%state0 = plasticState(ph)%state ! ToDo: this could be done centrally
end associate
!--------------------------------------------------------------------------------------------------
@ -414,7 +412,7 @@ module subroutine plastic_dislotungsten_results(ph,group)
'mobile dislocation density','1/m²')
case('rho_dip')
if(prm%sum_N_sl>0) call results_writeDataset(group,stt%rho_dip,trim(prm%output(o)), &
'dislocation dipole density''1/m²')
'dislocation dipole density','1/m²')
case('gamma_sl')
if(prm%sum_N_sl>0) call results_writeDataset(group,stt%gamma_sl,trim(prm%output(o)), &
'plastic shear','1')

View File

@ -467,8 +467,6 @@ module function plastic_dislotwin_init() result(myPlasticity)
allocate(dst%tau_r_tr (prm%sum_N_tr,Nmembers),source=0.0_pReal)
allocate(dst%V_tr (prm%sum_N_tr,Nmembers),source=0.0_pReal)
plasticState(ph)%state0 = plasticState(ph)%state ! ToDo: this could be done centrally
end associate
!--------------------------------------------------------------------------------------------------

View File

@ -140,8 +140,6 @@ module function plastic_isotropic_init() result(myPlasticity)
! global alias
plasticState(ph)%slipRate => plasticState(ph)%dotState(2:2,:)
plasticState(ph)%state0 = plasticState(ph)%state ! ToDo: this could be done centrally
end associate
!--------------------------------------------------------------------------------------------------

View File

@ -213,8 +213,6 @@ module function plastic_kinehardening_init() result(myPlasticity)
stt%gamma0 => plasticState(ph)%state (startIndex :endIndex ,:)
dlt%gamma0 => plasticState(ph)%deltaState(startIndex-o:endIndex-o,:)
plasticState(ph)%state0 = plasticState(ph)%state ! ToDo: this could be done centrally
end associate
!--------------------------------------------------------------------------------------------------

View File

@ -1767,22 +1767,22 @@ end subroutine kinetics
!> @brief returns copy of current dislocation densities from state
!> @details raw values is rectified
!--------------------------------------------------------------------------------------------------
pure function getRho(ph,en)
pure function getRho(ph,en) result(rho)
integer, intent(in) :: ph, en
real(pReal), dimension(param(ph)%sum_N_sl,10) :: getRho
real(pReal), dimension(param(ph)%sum_N_sl,10) :: rho
associate(prm => param(ph))
getRho = reshape(state(ph)%rho(:,en),[prm%sum_N_sl,10])
rho = reshape(state(ph)%rho(:,en),[prm%sum_N_sl,10])
! ensure positive densities (not for imm, they have a sign)
getRho(:,mob) = max(getRho(:,mob),0.0_pReal)
getRho(:,dip) = max(getRho(:,dip),0.0_pReal)
rho(:,mob) = max(rho(:,mob),0.0_pReal)
rho(:,dip) = max(rho(:,dip),0.0_pReal)
where(abs(getRho) < max(prm%rho_min/geom(ph)%V_0(en)**(2.0_pReal/3.0_pReal),prm%rho_significant)) &
getRho = 0.0_pReal
where(abs(rho) < max(prm%rho_min/geom(ph)%V_0(en)**(2.0_pReal/3.0_pReal),prm%rho_significant)) &
rho = 0.0_pReal
end associate
@ -1793,22 +1793,22 @@ end function getRho
!> @brief returns copy of current dislocation densities from state
!> @details raw values is rectified
!--------------------------------------------------------------------------------------------------
pure function getRho0(ph,en)
pure function getRho0(ph,en) result(rho0)
integer, intent(in) :: ph, en
real(pReal), dimension(param(ph)%sum_N_sl,10) :: getRho0
real(pReal), dimension(param(ph)%sum_N_sl,10) :: rho0
associate(prm => param(ph))
getRho0 = reshape(state0(ph)%rho(:,en),[prm%sum_N_sl,10])
rho0 = reshape(state0(ph)%rho(:,en),[prm%sum_N_sl,10])
! ensure positive densities (not for imm, they have a sign)
getRho0(:,mob) = max(getRho0(:,mob),0.0_pReal)
getRho0(:,dip) = max(getRho0(:,dip),0.0_pReal)
rho0(:,mob) = max(rho0(:,mob),0.0_pReal)
rho0(:,dip) = max(rho0(:,dip),0.0_pReal)
where (abs(getRho0) < max(prm%rho_min/geom(ph)%V_0(en)**(2.0_pReal/3.0_pReal),prm%rho_significant)) &
getRho0 = 0.0_pReal
where (abs(rho0) < max(prm%rho_min/geom(ph)%V_0(en)**(2.0_pReal/3.0_pReal),prm%rho_significant)) &
rho0 = 0.0_pReal
end associate

View File

@ -265,8 +265,6 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma'
plasticState(ph)%state0 = plasticState(ph)%state ! ToDo: this could be done centrally
end associate
!--------------------------------------------------------------------------------------------------

View File

@ -178,8 +178,7 @@ subroutine fromQuaternion(self,qu)
class(rotation), intent(out) :: self
real(pReal), dimension(4), intent(in) :: qu
if (dNeq(norm2(qu),1.0_pReal)) &
call IO_error(402,ext_msg='fromQuaternion')
if (dNeq(norm2(qu),1.0_pReal,1.0e-8_pReal)) call IO_error(402,ext_msg='fromQuaternion')
self%q = qu
@ -639,7 +638,7 @@ function om2ax(om) result(ax)
else
call dgeev('N','V',3,om_,3,Wr,Wi,devNull,3,VR,3,work,size(work,1),ierr)
if (ierr /= 0) error stop 'LAPACK error'
#if defined(__GFORTRAN__) && __GNUC__<9 || defined(__INTEL_COMPILER) && INTEL_COMPILER<1800 || defined(__PGI)
#if defined(__GFORTRAN__) && __GNUC__<9
i = maxloc(merge(1,0,cEq(cmplx(Wr,Wi,pReal),cmplx(1.0_pReal,0.0_pReal,pReal),tol=1.0e-14_pReal)),dim=1)
#else
i = findloc(cEq(cmplx(Wr,Wi,pReal),cmplx(1.0_pReal,0.0_pReal,pReal),tol=1.0e-14_pReal),.true.,dim=1) !find eigenvalue (1,0)
@ -1394,7 +1393,7 @@ end function conjugate_quaternion
!--------------------------------------------------------------------------------------------------
!> @brief Check correctness of some rotations functions.
!--------------------------------------------------------------------------------------------------
subroutine selfTest
subroutine selfTest()
type(rotation) :: R
real(pReal), dimension(4) :: qu, ax, ro
@ -1405,7 +1404,7 @@ subroutine selfTest
integer :: i
do i = 1, 10
do i = 1, 20
#if defined(__GFORTRAN__) && __GNUC__<9
if(i<7) cycle
@ -1433,6 +1432,7 @@ subroutine selfTest
sin(2.0_pReal*PI*x(1))*A]
if(qu(1)<0.0_pReal) qu = qu * (-1.0_pReal)
endif
if(.not. quaternion_equal(om2qu(qu2om(qu)),qu)) error stop 'om2qu/qu2om'
if(.not. quaternion_equal(eu2qu(qu2eu(qu)),qu)) error stop 'eu2qu/qu2eu'
if(.not. quaternion_equal(ax2qu(qu2ax(qu)),qu)) error stop 'ax2qu/qu2ax'
@ -1479,6 +1479,8 @@ subroutine selfTest
if(all(dNeq(R%rotTensor4(R%rotTensor4(t3333),active=.true.),t3333,1.0e-12_pReal))) &
error stop 'rotTensor4'
call R%fromQuaternion(qu * (1.0_pReal + merge(+5.e-9_pReal,-5.e-9_pReal, mod(i,2) == 0))) ! allow reasonable tolerance for ASCII/YAML
enddo
contains