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

Merge remote-tracking branch 'origin/development' into select-mu

This commit is contained in:
Sharan Roongta 2022-11-29 12:12:02 +01:00
parent 8b5fa37428 afa0a644ed
commit 599e4472e8
81 changed files with 3826 additions and 921 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
.gitlab-ci.yml
View File

@ -47,9 +47,9 @@ variables:
PETSC_INTELLLVM: "Libraries/PETSc/3.16.3/oneAPI-2022.0.1-IntelMPI-2021.5.0"
PETSC_INTEL: "Libraries/PETSc/3.16.5/Intel-2022.0.1-IntelMPI-2021.5.0"
# ++++++++++++ MSC Marc +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
MSC: "FEM/MSC/2022.1"
MSC: "FEM/MSC/2022.2"
IntelMarc: "Compiler/Intel/19.1.2 Libraries/IMKL/2020"
HDF5Marc: "HDF5/1.12.1/Intel-19.1.2"
HDF5Marc: "HDF5/1.12.2/Intel-19.1.2"
###################################################################################################

2
PRIVATE

@ -1 +1 @@
Subproject commit 7b1ad767256f796f56514d8888027499fe777132
Subproject commit 81f5f24d076a623e6052c234825c591267915285

2
VERSION
View File

@ -1 +1 @@
3.0.0-alpha7-86-g72408f389
3.0.0-alpha7-168-g741aeab1d

2
examples/config/numerics.yaml
View File

@ -80,4 +80,4 @@ commercialFEM:
generic:
random_seed: 0 # fixed seeding for pseudo-random number generator, Default 0: use random seed.
residualStiffness: 1.0e-6 # non-zero residual damage.
phi_min: 1.0e-6 # non-zero residual damage.

6
examples/config/phase/mechanical/eigen/thermalexpansion_AISI304.yaml
View File

@ -6,8 +6,8 @@ references:
https://doi.org/10.1016/0040-6031(93)80437-F,
fit to Fig. 6 (T_min=100K, T_max=1400K)
A_11: 2.068e-08
A_11,T: 1.579e-09
A_11,T^2: 3.449e-13
Alpha_11: 2.068e-08
Alpha_11,T: 1.579e-09
Alpha_11,T^2: 3.449e-13
T_ref: 293.15

2
examples/config/phase/mechanical/eigen/thermalexpansion_Al.yaml
View File

@ -4,4 +4,4 @@ references:
- https://en.wikipedia.org/wiki/Thermal_expansion,
293.15K
A_11: 23.1e-6
Alpha_11: 23.1e-6

2
examples/config/phase/mechanical/eigen/thermalexpansion_Au.yaml
View File

@ -4,4 +4,4 @@ references:
- https://en.wikipedia.org/wiki/Thermal_expansion,
293.15K
A_11: 14.e-6
Alpha_11: 14.e-6

6
examples/config/phase/mechanical/eigen/thermalexpansion_C35E.yaml
View File

@ -4,8 +4,8 @@ references:
- https://commons.wikimedia.org/wiki/File:Coefficient_dilatation_lineique_aciers.svg,
fit to image description (Scilab code)
A_11: 12.70371e-6
A_11,T: 7.54e-9
A_11,T^2: -1.0e-11
Alpha_11: 12.70371e-6
Alpha_11,T: 7.54e-9
Alpha_11,T^2: -1.0e-11
T_ref: 273.0

2
examples/config/phase/mechanical/eigen/thermalexpansion_Cu.yaml
View File

@ -4,4 +4,4 @@ references:
- https://en.wikipedia.org/wiki/Thermal_expansion,
293.15K
A_11: 17.e-6
Alpha_11: 17.e-6

2
examples/config/phase/mechanical/eigen/thermalexpansion_Fe.yaml
View File

@ -4,4 +4,4 @@ references:
- https://en.wikipedia.org/wiki/Thermal_expansion,
293.15K
A_11: 11.8e-6
Alpha_11: 11.8e-6

12
examples/config/phase/mechanical/eigen/thermalexpansion_Sn-beta.yaml
View File

@ -6,12 +6,12 @@ references:
https://doi.org/10.1107/S0365110X62000742,
fit to Tab. 2 (T_min=30ºC, T_max=210ºC)
A_11: 1.639e-05
A_11,T: 1.799e-08
A_11,T^2: 1.734e-10
Alpha_11: 1.639e-05
Alpha_11,T: 1.799e-08
Alpha_11,T^2: 1.734e-10
A_33: 3.263e-05
A_33,T: 1.387e-08
A_33,T^2: 5.794e-10
Alpha_33: 3.263e-05
Alpha_33,T: 1.387e-08
Alpha_33,T^2: 5.794e-10
T_ref: 293.15

2
examples/config/phase/mechanical/eigen/thermalexpansion_W.yaml
View File

@ -4,4 +4,4 @@ references:
- https://en.wikipedia.org/wiki/Thermal_expansion,
293.15K
A_11: 4.5e-6
Alpha_11: 4.5e-6

4
examples/config/phase/mechanical/eigen/thermalexpansion_X20Cr13.yaml
View File

@ -4,7 +4,7 @@ references:
- https://commons.wikimedia.org/wiki/File:Coefficient_dilatation_lineique_aciers.svg,
fit to image description (Scilab code)
A_11: 11.365e-6
A_11,T: 5.0e-9
Alpha_11: 11.365e-6
Alpha_11,T: 5.0e-9
T_ref: 273.0

2
install/MarcMentat/2020/Marc_tools/include_linux64.patch vendored
View File

@ -46,7 +46,7 @@
+
+# determine DAMASK version
+if test -n "$DAMASK_USER"; then
+ DAMASKROOT=`dirname $DAMASK_USER`/..
+ DAMASKROOT=`dirname $DAMASK_USER`/../..
+ read DAMASKVERSION < $DAMASKROOT/VERSION
+ DAMASKVERSION="'"$DAMASKVERSION"'"
+else

2
install/MarcMentat/2020/Marc_tools/run_damask_hmp.patch vendored
View File

@ -251,7 +251,7 @@
- usersub=$usersubname
- fi
-
+ userobj=$usermoext.o
+ userobj=$usernoext.o
fi
cat > $jid.runmarcscript << END4
if test "$user"

2
install/MarcMentat/2020/Marc_tools/run_damask_lmp.patch vendored
View File

@ -251,7 +251,7 @@
- usersub=$usersubname
- fi
-
+ userobj=$usermoext.o
+ userobj=$usernoext.o
fi
cat > $jid.runmarcscript << END4
if test "$user"

2
install/MarcMentat/2020/Marc_tools/run_damask_mp.patch vendored
View File

@ -251,7 +251,7 @@
- usersub=$usersubname
- fi
-
+ userobj=$usermoext.o
+ userobj=$usernoext.o
fi
cat > $jid.runmarcscript << END4
if test "$user"

2
install/MarcMentat/2020/Mentat_menus/job_run.ms.patch vendored
View File

@ -66,7 +66,7 @@
+ label {
+ position -32 +6
+ size 12 6
+ text "O2 / OpenMP"
+ text "O3 / OpenMP"
+ border_width 1
+ border_color black
+ }

2
install/MarcMentat/2021.2/Marc_tools/include_linux64.patch vendored
View File

@ -46,7 +46,7 @@
+
+# determine DAMASK version
+if test -n "$DAMASK_USER"; then
+ DAMASKROOT=`dirname $DAMASK_USER`/..
+ DAMASKROOT=`dirname $DAMASK_USER`/../..
+ read DAMASKVERSION < $DAMASKROOT/VERSION
+ DAMASKVERSION="'"$DAMASKVERSION"'"
+else

2
install/MarcMentat/2021.2/Marc_tools/run_damask_hmp.patch vendored
View File

@ -284,7 +284,7 @@
- usersub=$usersubname
- fi
-
+ userobj=$usermoext.o
+ userobj=$usernoext.o
fi
cat > $jid.runmarcscript << END4
if test "$user"

2
install/MarcMentat/2021.2/Marc_tools/run_damask_lmp.patch vendored
View File

@ -284,7 +284,7 @@
- usersub=$usersubname
- fi
-
+ userobj=$usermoext.o
+ userobj=$usernoext.o
fi
cat > $jid.runmarcscript << END4
if test "$user"

2
install/MarcMentat/2021.2/Marc_tools/run_damask_mp.patch vendored
View File

@ -284,7 +284,7 @@
- usersub=$usersubname
- fi
-
+ userobj=$usermoext.o
+ userobj=$usernoext.o
fi
cat > $jid.runmarcscript << END4
if test "$user"

2
install/MarcMentat/2021.2/Mentat_menus/job_run.ms.patch vendored
View File

@ -65,7 +65,7 @@
+ label {
+ position -32 +6
+ size 12 6
+ text "O2 / OpenMP"
+ text "O3 / OpenMP"
+ border_width 1
+ border_color black
+ }

2
install/MarcMentat/2021.3.1/Marc_tools/include_linux64.patch vendored
View File

@ -50,7 +50,7 @@
+# determine DAMASK version
+if test -n "$DAMASK_USER"; then
+ DAMASKROOT=`dirname $DAMASK_USER`/..
+ DAMASKROOT=`dirname $DAMASK_USER`/../..
+ read DAMASKVERSION < $DAMASKROOT/VERSION
+ DAMASKVERSION="'"$DAMASKVERSION"'"
+else

2
install/MarcMentat/2021.3.1/Marc_tools/run_damask_hmp.patch vendored
View File

@ -283,7 +283,7 @@
- usersub=$usersubname
- fi
-
+ userobj=$usermoext.o
+ userobj=$usernoext.o
fi
cat > $jid.runmarcscript << END4
if test "$user"

2
install/MarcMentat/2021.3.1/Marc_tools/run_damask_lmp.patch vendored
View File

@ -283,7 +283,7 @@
- usersub=$usersubname
- fi
-
+ userobj=$usermoext.o
+ userobj=$usernoext.o
fi
cat > $jid.runmarcscript << END4
if test "$user"

2
install/MarcMentat/2021.3.1/Marc_tools/run_damask_mp.patch vendored
View File

@ -301,7 +301,7 @@
- usersub=$usersubname
- fi
-
+ userobj=$usermoext.o
+ userobj=$usernoext.o
fi
cat > $jid.runmarcscript << END4
if test "$user"

2
install/MarcMentat/2021.3.1/Mentat_menus/job_run.ms.patch vendored
View File

@ -67,7 +67,7 @@
+ label {
+ position -32 +6
+ size 12 6
+ text "O2 / OpenMP"
+ text "O3 / OpenMP"
+ border_width 1
+ border_color black
+ }

2
install/MarcMentat/2022.1/Marc_tools/include_linux64.patch vendored
View File

@ -40,7 +40,7 @@
+# determine DAMASK version
+if test -n "$DAMASK_USER"; then
+ DAMASKROOT=`dirname $DAMASK_USER`/..
+ DAMASKROOT=`dirname $DAMASK_USER`/../..
+ read DAMASKVERSION < $DAMASKROOT/VERSION
+ DAMASKVERSION="'"$DAMASKVERSION"'"
+else

2
install/MarcMentat/2022.1/Marc_tools/run_damask_hmp.patch vendored
View File

@ -228,7 +228,7 @@
- usersub=$usersubname
- fi
-
+ userobj=$usermoext.o
+ userobj=$usernoext.o
fi
cat > $jid.runmarcscript << END4
if test "$user"

2
install/MarcMentat/2022.1/Marc_tools/run_damask_lmp.patch vendored
View File

@ -228,7 +228,7 @@
- usersub=$usersubname
- fi
-
+ userobj=$usermoext.o
+ userobj=$usernoext.o
fi
cat > $jid.runmarcscript << END4
if test "$user"

2
install/MarcMentat/2022.1/Marc_tools/run_damask_mp.patch vendored
View File

@ -228,7 +228,7 @@
- usersub=$usersubname
- fi
-
+ userobj=$usermoext.o
+ userobj=$usernoext.o
fi
cat > $jid.runmarcscript << END4
if test "$user"

2
install/MarcMentat/2022.1/Mentat_menus/job_run.ms.patch vendored
View File

@ -67,7 +67,7 @@
+ label {
+ position -32 +6
+ size 12 6
+ text "O2 / OpenMP"
+ text "O3 / OpenMP"
+ border_width 1
+ border_color black
+ }

49
install/MarcMentat/2022.2/Marc_tools/comp_damask_hmp.patch vendored Normal file
View File

@ -0,0 +1,49 @@
---
+++
@@ -6,18 +6,27 @@
DIR=$1
user=$3
program=$4
+usernoext=$user
+usernoext=`dirname $usernoext`/`$BASENAME $usernoext .f`
+usernoext=`dirname $usernoext`/`$BASENAME $usernoext .F`
+usernoext=`dirname $usernoext`/`$BASENAME $usernoext .for`
+usernoext=`dirname $usernoext`/`$BASENAME $usernoext .f90`
+
+# add BLAS options for linking
+ BLAS="%BLAS%"
+
. $DIR/tools/include
DIRJOB=$2
cd $DIRJOB
-echo "Compiling and linking user subroutine $user.f on host `hostname`"
+echo "Compiling and linking user subroutine $user on host `hostname`"
echo "program: $program"
- $FORTRAN $user.f || \
+ $DFORTHIGHMP $user || \
{
- echo "$0: compile failed for $user.f"
+ echo "$0: compile failed for $user"
exit 1
}
/bin/rm $program 2>/dev/null
- userobj=$user.o
+ userobj=$usernoext.o
$LOAD ${program} $DIR/lib/main.o\
@@ -33,9 +42,13 @@
$TKLIBS \
$MRCLIBS \
$METISLIBS \
+ $BLAS \
$SYSLIBS || \
{
- echo "$0: link failed for $user.o on host `hostname`"
+ echo "$0: link failed for $usernoext.o on host `hostname`"
exit 1
}
/bin/rm $userobj
+ /bin/rm $DIRJOB/*.mod
+ /bin/rm $DIRJOB/*.smod
+ /bin/rm $DIRJOB/*_genmod.f90

49
install/MarcMentat/2022.2/Marc_tools/comp_damask_lmp.patch vendored Normal file
View File

@ -0,0 +1,49 @@
---
+++
@@ -6,18 +6,27 @@
DIR=$1
user=$3
program=$4
+usernoext=$user
+usernoext=`dirname $usernoext`/`$BASENAME $usernoext .f`
+usernoext=`dirname $usernoext`/`$BASENAME $usernoext .F`
+usernoext=`dirname $usernoext`/`$BASENAME $usernoext .for`
+usernoext=`dirname $usernoext`/`$BASENAME $usernoext .f90`
+
+# add BLAS options for linking
+ BLAS="%BLAS%"
+
. $DIR/tools/include
DIRJOB=$2
cd $DIRJOB
-echo "Compiling and linking user subroutine $user.f on host `hostname`"
+echo "Compiling and linking user subroutine $user on host `hostname`"
echo "program: $program"
- $FORTRAN $user.f || \
+ $DFORTRANLOWMP $user || \
{
- echo "$0: compile failed for $user.f"
+ echo "$0: compile failed for $user"
exit 1
}
/bin/rm $program 2>/dev/null
- userobj=$user.o
+ userobj=$usernoext.o
$LOAD ${program} $DIR/lib/main.o\
@@ -33,9 +42,13 @@
$TKLIBS \
$MRCLIBS \
$METISLIBS \
+ $BLAS \
$SYSLIBS || \
{
- echo "$0: link failed for $user.o on host `hostname`"
+ echo "$0: link failed for $usernoext.o on host `hostname`"
exit 1
}
/bin/rm $userobj
+ /bin/rm $DIRJOB/*.mod
+ /bin/rm $DIRJOB/*.smod
+ /bin/rm $DIRJOB/*_genmod.f90

49
install/MarcMentat/2022.2/Marc_tools/comp_damask_mp.patch vendored Normal file
View File

@ -0,0 +1,49 @@
---
+++
@@ -6,18 +6,27 @@
DIR=$1
user=$3
program=$4
+usernoext=$user
+usernoext=`dirname $usernoext`/`$BASENAME $usernoext .f`
+usernoext=`dirname $usernoext`/`$BASENAME $usernoext .F`
+usernoext=`dirname $usernoext`/`$BASENAME $usernoext .for`
+usernoext=`dirname $usernoext`/`$BASENAME $usernoext .f90`
+
+# add BLAS options for linking
+ BLAS="%BLAS%"
+
. $DIR/tools/include
DIRJOB=$2
cd $DIRJOB
-echo "Compiling and linking user subroutine $user.f on host `hostname`"
+echo "Compiling and linking user subroutine $user on host `hostname`"
echo "program: $program"
- $FORTRAN $user.f || \
+ $DFORTRANMP $user || \
{
- echo "$0: compile failed for $user.f"
+ echo "$0: compile failed for $user"
exit 1
}
/bin/rm $program 2>/dev/null
- userobj=$user.o
+ userobj=$usernoext.o
$LOAD ${program} $DIR/lib/main.o\
@@ -33,9 +42,13 @@
$TKLIBS \
$MRCLIBS \
$METISLIBS \
+ $BLAS \
$SYSLIBS || \
{
- echo "$0: link failed for $user.o on host `hostname`"
+ echo "$0: link failed for $usernoext.o on host `hostname`"
exit 1
}
/bin/rm $userobj
+ /bin/rm $DIRJOB/*.mod
+ /bin/rm $DIRJOB/*.smod
+ /bin/rm $DIRJOB/*_genmod.f90

75
install/MarcMentat/2022.2/Marc_tools/include_linux64.patch vendored Normal file
View File

@ -0,0 +1,75 @@
---
+++
@@ -166,6 +166,15 @@
MARC_COSIM_LIB="$MSCCOSIM_HOME/CoSim$MSCCOSIM_VERSION/Dcosim$MSCCOSIM_VERSION/lib"
fi
+# DAMASK uses the HDF5 compiler wrapper around the Intel compiler
+H5FC=$(h5fc -shlib -show)
+if [[ "$H5FC" == *"$dir is"* ]]; then
+ H5FC=$(echo $(echo "$H5FC" | tail -n1) | sed -e "s/\-shlib/-fPIC -qopenmp/g")
+ H5FC=${H5FC%-lmpifort*}
+fi
+HDF5_LIB=${H5FC//*ifort/}
+FCOMP="$H5FC"
+
# AEM
if test "$MARCDLLOUTDIR" = ""; then
DLLOUTDIR="$MARC_LIB"
@@ -594,7 +603,7 @@
PROFILE=" $PROFILE -pg"
fi
-FORT_OPT="-c -assume byterecl -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"
@@ -607,7 +616,7 @@
FORT_OPT=" $FORT_OPT -save -zero"
fi
if test "$MARCHDF_HDF" = "HDF"; then
- FORT_OPT="$FORT_OPT -DMARCHDF_HDF=$MARCHDF_HDF $HDF_INCLUDE"
+ FORT_OPT="$FORT_OPT -DMARCHDF=$MARCHDF_HDF"
fi
FORTLOW="$FCOMP $FORT_OPT $PROFILE -O0 $I8FFLAGS -I$MARC_SOURCE/common \
@@ -621,6 +630,29 @@
# for compiling free form f90 files. high opt, integer(4)
FORTF90="$FCOMP -c -O3"
+# determine DAMASK version
+if test -n "$DAMASK_USER"; then
+ DAMASKROOT=`dirname $DAMASK_USER`/../..
+ read DAMASKVERSION < $DAMASKROOT/VERSION
+ DAMASKVERSION="'"$DAMASKVERSION"'"
+else
+ DAMASKVERSION="'N/A'"
+fi
+
+# DAMASK compiler calls
+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=2022.1 -DDAMASKVERSION=$DAMASKVERSION \
+ -qopenmp -qopenmp-threadprivate=compat\
+ $MUMPS_INCLUDE $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS -I$KDTREE2_MOD -I$MARC_MOD"
+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=2022.1 -DDAMASKVERSION=$DAMASKVERSION \
+ -qopenmp -qopenmp-threadprivate=compat\
+ $MUMPS_INCLUDE $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS -I$KDTREE2_MOD -I$MARC_MOD"
+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=2022.1 -DDAMASKVERSION=$DAMASKVERSION \
+ -qopenmp -qopenmp-threadprivate=compat\
+ $MUMPS_INCLUDE $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS -I$KDTREE2_MOD -I$MARC_MOD"
+
if test "$MARCDEBUG" = "ON"
then
FORTLOW="$FCOMP $FORT_OPT $PROFILE $I8FFLAGS -I$MARC_SOURCE/common \
@@ -778,7 +810,7 @@
SOLVERLIBS="${BCSSOLVERLIBS} ${VKISOLVERLIBS} ${CASISOLVERLIBS} ${MF2SOLVERLIBS} \
-L$MARC_MKL \
- $MARC_LIB/blas_src.a ${ACSI_LIB}/ACSI_MarcLib.a $KDTREE2_LIB/libkdtree2.a $MARC_LIB/libtetmeshinterface.a $MARC_LIB/libcaefatigueinterface.a -L$MARC_LIB -lmkl_blacs_intelmpi_ilp64 -lmkl_scalapack_ilp64 -lmkl_intel_ilp64 -lmkl_intel_thread -lmkl_core -liomp5 -ltetmesh -lmeshgems -lmg-tetra -lmeshgems_stubs $HDF_LIBS $SOLVER2LIBS"
+ $MARC_LIB/blas_src.a ${ACSI_LIB}/ACSI_MarcLib.a $KDTREE2_LIB/libkdtree2.a $MARC_LIB/libtetmeshinterface.a $MARC_LIB/libcaefatigueinterface.a -L$MARC_LIB -lmkl_blacs_intelmpi_ilp64 -lmkl_scalapack_ilp64 -lmkl_intel_ilp64 -lmkl_intel_thread -lmkl_core -liomp5 -ltetmesh -lmeshgems -lmg-tetra -lmeshgems_stubs $HDF5_LIB $SOLVER2LIBS"
SOLVERLIBS_DLL=${SOLVERLIBS}
if test "$AEM_DLL" -eq 1

517
install/MarcMentat/2022.2/Marc_tools/run_damask_hmp.patch vendored Normal file
View File

@ -0,0 +1,517 @@
---
+++
@@ -136,6 +136,11 @@
# is created. For job running in the background, the log #
# file is always created. Default is "yes" #
##############################################################################
+# remove all Mentat paths from LD_LIBRARY_PATH
+LD_LIBRARY_PATH=:$LD_LIBRARY_PATH:
+LD_LIBRARY_PATH=${LD_LIBRARY_PATH//+([!(:)])mentat2022.2+([!(:)])/:}
+LD_LIBRARY_PATH=${LD_LIBRARY_PATH//+([(:)])/:}
+LD_LIBRARY_PATH=${LD_LIBRARY_PATH#:}; LD_LIBRARY_PATH=${LD_LIBRARY_PATH%:}
# set DIR to the directory in which this script is
REALCOM="`/bin/ls -l $0 |awk '{ print $NF; }'`"
DIR=`dirname $REALCOM`
@@ -302,7 +307,23 @@
. "$DIR/getarch"
+
+# getting user subroutine file name
+found=0
+for i in "$@"; do
+ if test $found = 1; then
+ DAMASK_USER=$i
+ found=0
+ fi
+ case $i in
+ -u* | -U*)
+ found=1
+ ;;
+ esac
+done
+# sourcing include_linux64 (needs DAMASK_USER to be set)
. $MARC_INCLUDE
+
#
#
@@ -405,7 +426,7 @@
did=
vid=
user=
-usersubname=
+usernoext=
objs=
qid=background
cpu=
@@ -676,50 +697,19 @@
esac
;;
-u* | -U*)
- user=`dirname $value`/`$BASENAME $value .f`
- usersubname=$user
- basefile=`$BASENAME $value`
- if test ${basefile##*.} = f
- then
- user=`dirname $value`/`$BASENAME $value .f`
- usersubname=$user.f
- elif test ${basefile##*.} = F
- then
- user=`dirname $value`/`$BASENAME $value .F`
- usersubname=$user.F
- elif test ${basefile##*.} = f90
- then
- user=`dirname $value`/`$BASENAME $value .f90`
- usersubname=$user.f90
- elif test ${basefile##*.} = F90
- then
- user=`dirname $value`/`$BASENAME $value .F90`
- usersubname=$user.F90
- fi
+ user=$value
case $user in
\/*)
;;
*)
user=`pwd`/$user
- usersubname=`pwd`/$usersubname
;;
esac
- if test ! -f $usersubname
- then
- if test -f $usersubname.f
- then
- usersubname=$usersubname.f
- elif test -f $usersubname.F
- then
- usersubname=$usersubname.F
- elif test -f $usersubname.f90
- then
- usersubname=$usersubname.f90
- elif test -f $usersubname.F90
- then
- usersubname=$usersubname.F90
- fi
- fi
+ usernoext=$user
+ usernoext=`dirname $usernoext`/`$BASENAME $usernoext .f`
+ usernoext=`dirname $usernoext`/`$BASENAME $usernoext .F`
+ usernoext=`dirname $usernoext`/`$BASENAME $usernoext .for`
+ usernoext=`dirname $usernoext`/`$BASENAME $usernoext .f90`
;;
-obj | -OBJ)
objs="$value"
@@ -1207,12 +1197,12 @@
fi
fi
fi
- if test "$usersubname"
+ if test "$user"
then
- if test ! -f $usersubname
+ if test ! -f $user
then
error="$error
-user subroutine file $usersubname not accessible"
+user subroutine file $user not accessible"
fi
fi
if test "$objs"
@@ -1531,7 +1521,7 @@
Marc shared lib : $progdll
Version type : $mode
Job ID : $DIRJID/$jid$extra_job_info
-User subroutine name : $usersubname
+User subroutine name : $user
User objects/libs : $objs
Restart file job ID : $rid
Substructure file ID : $sid
@@ -1564,7 +1554,7 @@
Marc shared lib : $progdll
Version type : $mode
Job ID : $DIRJID/$jid$extra_job_info
-User subroutine name : $usersubname
+User subroutine name : $user
User objects/libs : $objs
Restart file job ID : $rid
Substructure file ID : $sid
@@ -1687,7 +1677,7 @@
;;
esac
fi
- $ECHO "User subroutine name ($usersubname)? $ECHOTXT"
+ $ECHO "User subroutine name ($user)? $ECHOTXT"
read value
if test "$value"
then
@@ -1696,50 +1686,19 @@
user=
;;
*)
- user=`dirname $value`/`$BASENAME $value .f`
- usersubname=$user
- basefile=`$BASENAME $value`
- if test ${basefile##*.} = f
- then
- user=`dirname $value`/`$BASENAME $value .f`
- usersubname=$user.f
- elif test ${basefile##*.} = F
- then
- user=`dirname $value`/`$BASENAME $value .F`
- usersubname=$user.F
- elif test ${basefile##*.} = f90
- then
- user=`dirname $value`/`$BASENAME $value .f90`
- usersubname=$user.f90
- elif test ${basefile##*.} = F90
- then
- user=`dirname $value`/`$BASENAME $value .F90`
- usersubname=$user.F90
- fi
+ user=$value
case $user in
\/*)
;;
*)
user=`pwd`/$user
- usersubname=`pwd`/$usersubname
;;
esac
- if test ! -f $usersubname
- then
- if test -f $usersubname.f
- then
- usersubname=$usersubname.f
- elif test -f $usersubname.F
- then
- usersubname=$usersubname.F
- elif test -f $usersubname.f90
- then
- usersubname=$usersubname.f90
- elif test -f $usersubname.F90
- then
- usersubname=$usersubname.F90
- fi
- fi
+ usernoext=$user
+ usernoext=`dirname $usernoext`/`$BASENAME $usernoext .f`
+ usernoext=`dirname $usernoext`/`$BASENAME $usernoext .F`
+ usernoext=`dirname $usernoext`/`$BASENAME $usernoext .for`
+ usernoext=`dirname $usernoext`/`$BASENAME $usernoext .f90`
;;
esac
fi
@@ -2274,11 +2233,12 @@
#
# user subroutine used
#
+# add DAMASK options for linking
+ DAMASK="-lstdc++"
if test "$user"
then
-# program=$user.marc
- program=$DIRJOB/`$BASENAME $user .f`.marc
+ program=$usernoext.marc
case $program in
\/* | \.\/*)
bd=
@@ -2391,7 +2351,7 @@
fi
if test "$user"
then
- execpath=$DIRJOB/`$BASENAME $user .f`.marc
+ execpath=$usernoext.marc
usersub=1
fi
export execpath
@@ -3274,44 +3234,27 @@
echo
if test "$user"
then
- userobj=$DIRJOB/`$BASENAME $user .f`.o
- basefile=`$BASENAME $usersubname`
- if test ${basefile##*.} = f
- then
- usersub=$DIRJOB/`$BASENAME $user .f`.F
- ln -sf "$user.f" "$usersub"
- else
- usersub=$usersubname
- fi
-
+ userobj=$usernoext.o
fi
cat > $jid.runmarcscript << END4
if test "$user"
then
- if test ${basefile##*.} = f
- then
- ln -sf "$user.f" "$usersub"
- fi
if test $MACHINENAME = "CRAY"
then
- $FORTRAN $usersub || \
+ $DFORTHIGHMP $user || \
{
- echo "$0: compile failed for $user.f"
+ echo "$0: compile failed for $user"
exit 1
}
/bin/rm $program 2>/dev/null
else
- $FORTRAN $usersub -o $userobj || \
+ $DFORTHIGHMP $user -o $userobj || \
{
- echo "$0: compile failed for $user.f"
+ echo "$0: compile failed for $user"
exit 1
}
/bin/rm $program 2>/dev/null
fi
- if test ${basefile##*.} = f
- then
- /bin/rm -f "$usersub"
- fi
fi
@@ -3331,6 +3274,7 @@
$TKLIBS \
$MRCLIBS \
$METISLIBS \
+ $DAMASK \
$SFLIB \
$OPENSSL_LIB \
$SYSLIBS \
@@ -3344,6 +3288,9 @@
prgsav=yes
fi
/bin/rm $userobj 2>/dev/null
+/bin/rm $DIRJOB/*.mod 2>/dev/null
+/bin/rm $DIRJOB/*.smod 2>/dev/null
+/bin/rm $DIRJOB/*_genmod.f90 2>/dev/null
#
# run marc
@@ -3390,7 +3337,7 @@
fi
else
if test $cpdll = yes; then
- filename=`basename $usersubname .f`
+ filename=$usernoext
/bin/cp $DIRJOB/$marcdll $DIRJOB/${filename}_$marcdll 2>/dev/null
fi
if test $rmdll = yes
@@ -3556,7 +3503,7 @@
# first copy over the user sub if local directories
if test ${dirstatus[$counter]} = "local"
then
- $RCP $user.f $i:$DIR1/
+ $RCP $user $i:$DIR1/
fi
# do the compilation on the other machine
if test ${dirstatus[$counter]} = "shared"
@@ -3569,21 +3516,21 @@
remoteuser=$DIR1/`$BASENAME $user`
$RSH $i /bin/rm $remoteprog 2> /dev/null
echo
- $RSH $i $DIR2/tools/comp_user $DIR2 $DIR1 $remoteuser $remoteprog
+ $RSH $i $DIR2/tools/comp_damask_hmp $DIR2 $DIR1 $remoteuser $remoteprog
# check if successful, the new executable should be there
line=`$RSH $i /bin/ls $remoteprog 2> /dev/null`
if test "$line"
then
echo compilation and linking successful on host $i
else
- echo "$0: compile failed for $user.f on host $i"
+ echo "$0: compile failed for $user on host $i"
echo " $PRODUCT Exit number 3"
exit 1
fi
# remove the user subroutine on remote machine
if test ${dirstatus[$counter]} = "local"
then
- $RSH $i /bin/rm $remoteuser.f 2> /dev/null
+ $RSH $i /bin/rm $remoteuser 2> /dev/null
fi
fi
fi
@@ -3593,39 +3540,27 @@
if test "$userhost"
then
echo
- echo "Compiling and linking user subroutine $user.f on host `hostname`"
- fi
- userobj=$DIRJOB/`$BASENAME $user .f`.o
- basefile=`$BASENAME $usersubname`
- if test ${basefile##*.} = f
- then
- usersub=$DIRJOB/`$BASENAME $user .f`.F
- ln -sf "$user.f" "$usersub"
- else
- usersub=$usersubname
+ echo "Compiling and linking user subroutine $user on host `hostname`"
fi
+ userobj=$usernoext.o
if test $MACHINENAME = "CRAY"
then
- $FORTRAN $usersub || \
+ $DFORTHIGHMP $user || \
{
- echo "$0: compile failed for $user.f"
+ echo "$0: compile failed for $user"
echo " $PRODUCT Exit number 3"
exit 1
}
/bin/rm $program 2>/dev/null
else
- $FORTRAN $usersub -o $userobj || \
+ $DFORTHIGHMP $user -o $userobj || \
{
- echo "$0: compile failed for $user.f"
+ echo "$0: compile failed for $user"
echo " $PRODUCT Exit number 3"
exit 1
}
/bin/rm $program 2>/dev/null
fi
- if test ${basefile##*.} = f
- then
- /bin/rm -f "$usersub"
- fi
fi # if test $user
@@ -3645,6 +3580,7 @@
$TKLIBS \
$MRCLIBS \
$METISLIBS \
+ $DAMASK \
$SFLIB \
$OPENSSL_LIB \
$SYSLIBS \
@@ -3686,6 +3622,9 @@
prgsav=yes
fi # if test $link
/bin/rm $userobj 2>/dev/null
+/bin/rm $DIRJOB/*.mod 2>/dev/null
+/bin/rm $DIRJOB/*.smod 2>/dev/null
+/bin/rm $DIRJOB/*_genmod.f90 2>/dev/null
#
# run marc
@@ -3779,7 +3718,7 @@
else
#dllrun >0
if test $cpdll = yes; then
- filename=`basename $usersubname .f`
+ filename=$usernoext
/bin/cp $DIRJOB/$marcdll $DIRJOB/${filename}_$marcdll 2>/dev/null
fi
if test $rmdll = yes;then
@@ -3904,7 +3843,7 @@
# first copy over the user sub if local directories
if test ${dirstatus[$counter]} = "local"
then
- $RCP $user.f $i:$DIR1/
+ $RCP $user $i:$DIR1/
fi
# do the compilation on the other machine
if test ${dirstatus[$counter]} = "shared"
@@ -3917,20 +3856,20 @@
remoteuser=$DIR1/`$BASENAME $user`
$RSH $i /bin/rm $remoteprog 2> /dev/null
echo
- $RSH $i $DIR2/tools/comp_user $DIR2 $DIR1 $remoteuser $remoteprog
+ $RSH $i $DIR2/tools/comp_damask_hmp $DIR2 $DIR1 $remoteuser $remoteprog
# check if successful, the new executable should be there
line=`$RSH $i /bin/ls $remoteprog 2> /dev/null`
if test "$line"
then
echo compilation and linking successful on host $i
else
- echo "$0: compile failed for $user.f on host $i"
+ echo "$0: compile failed for $user on host $i"
exit 1
fi
# remove the user subroutine on remote machine
if test ${dirstatus[$counter]} = "local"
then
- $RSH $i /bin/rm $remoteuser.f 2> /dev/null
+ $RSH $i /bin/rm $remoteuser 2> /dev/null
fi
fi
fi
@@ -3940,37 +3879,25 @@
if test "$userhost"
then
echo
- echo "Compiling and linking user subroutine $user.f on host `hostname`"
- fi
- userobj=$DIRJOB/`$BASENAME $user .f`.o
- basefile=`$BASENAME $usersubname`
- if test ${basefile##*.} = f
- then
- usersub=$DIRJOB/`$BASENAME $user .f`.F
- ln -sf "$user.f" "$usersub"
- else
- usersub=$usersubname
+ echo "Compiling and linking user subroutine $user on host `hostname`"
fi
+ userobj=$usernoext.o
if test $MACHINENAME = "CRAY"
then
- $FORTRAN $usersub || \
+ $DFORTHIGHMP $user || \
{
- echo "$0: compile failed for $user.f"
+ echo "$0: compile failed for $user"
exit 1
}
/bin/rm $program 2>/dev/null
else
- $FORTRAN $usersub -o $userobj || \
+ $DFORTHIGHMP $user -o $userobj || \
{
- echo "$0: compile failed for $user.f"
+ echo "$0: compile failed for $user"
exit 1
}
/bin/rm $program 2>/dev/null
fi
- if test ${basefile##*.} = f
- then
- /bin/rm -f "$usersub"
- fi
fi # if test $user
@@ -3990,6 +3917,7 @@
$TKLIBS \
$MRCLIBS \
$METISLIBS \
+ $DAMASK \
$SFLIB \
$OPENSSL_LIB \
$SYSLIBS \
@@ -4030,7 +3958,9 @@
prgsav=yes
fi # if test $link
/bin/rm $userobj 2>/dev/null
-
+/bin/rm $DIRJOB/*.mod 2>/dev/null
+/bin/rm $DIRJOB/*.smod 2>/dev/null
+/bin/rm $DIRJOB/*_genmod.f90 2>/dev/null
# done if no job id given
if test -z "$jid"
then
@@ -4149,7 +4079,7 @@
else
#dllrun >0
if test $cpdll = yes; then
- filename=`basename $usersubname .f`
+ filename=$usernoext
/bin/cp $DIRJOB/$marcdll $DIRJOB/${filename}_$marcdll 2>/dev/null
fi
if test $rmdll = yes;then

517
install/MarcMentat/2022.2/Marc_tools/run_damask_lmp.patch vendored Normal file
View File

@ -0,0 +1,517 @@
---
+++
@@ -136,6 +136,11 @@
# is created. For job running in the background, the log #
# file is always created. Default is "yes" #
##############################################################################
+# remove all Mentat paths from LD_LIBRARY_PATH
+LD_LIBRARY_PATH=:$LD_LIBRARY_PATH:
+LD_LIBRARY_PATH=${LD_LIBRARY_PATH//+([!(:)])mentat2022.2+([!(:)])/:}
+LD_LIBRARY_PATH=${LD_LIBRARY_PATH//+([(:)])/:}
+LD_LIBRARY_PATH=${LD_LIBRARY_PATH#:}; LD_LIBRARY_PATH=${LD_LIBRARY_PATH%:}
# set DIR to the directory in which this script is
REALCOM="`/bin/ls -l $0 |awk '{ print $NF; }'`"
DIR=`dirname $REALCOM`
@@ -302,7 +307,23 @@
. "$DIR/getarch"
+
+# getting user subroutine file name
+found=0
+for i in "$@"; do
+ if test $found = 1; then
+ DAMASK_USER=$i
+ found=0
+ fi
+ case $i in
+ -u* | -U*)
+ found=1
+ ;;
+ esac
+done
+# sourcing include_linux64 (needs DAMASK_USER to be set)
. $MARC_INCLUDE
+
#
#
@@ -405,7 +426,7 @@
did=
vid=
user=
-usersubname=
+usernoext=
objs=
qid=background
cpu=
@@ -676,50 +697,19 @@
esac
;;
-u* | -U*)
- user=`dirname $value`/`$BASENAME $value .f`
- usersubname=$user
- basefile=`$BASENAME $value`
- if test ${basefile##*.} = f
- then
- user=`dirname $value`/`$BASENAME $value .f`
- usersubname=$user.f
- elif test ${basefile##*.} = F
- then
- user=`dirname $value`/`$BASENAME $value .F`
- usersubname=$user.F
- elif test ${basefile##*.} = f90
- then
- user=`dirname $value`/`$BASENAME $value .f90`
- usersubname=$user.f90
- elif test ${basefile##*.} = F90
- then
- user=`dirname $value`/`$BASENAME $value .F90`
- usersubname=$user.F90
- fi
+ user=$value
case $user in
\/*)
;;
*)
user=`pwd`/$user
- usersubname=`pwd`/$usersubname
;;
esac
- if test ! -f $usersubname
- then
- if test -f $usersubname.f
- then
- usersubname=$usersubname.f
- elif test -f $usersubname.F
- then
- usersubname=$usersubname.F
- elif test -f $usersubname.f90
- then
- usersubname=$usersubname.f90
- elif test -f $usersubname.F90
- then
- usersubname=$usersubname.F90
- fi
- fi
+ usernoext=$user
+ usernoext=`dirname $usernoext`/`$BASENAME $usernoext .f`
+ usernoext=`dirname $usernoext`/`$BASENAME $usernoext .F`
+ usernoext=`dirname $usernoext`/`$BASENAME $usernoext .for`
+ usernoext=`dirname $usernoext`/`$BASENAME $usernoext .f90`
;;
-obj | -OBJ)
objs="$value"
@@ -1207,12 +1197,12 @@
fi
fi
fi
- if test "$usersubname"
+ if test "$user"
then
- if test ! -f $usersubname
+ if test ! -f $user
then
error="$error
-user subroutine file $usersubname not accessible"
+user subroutine file $user not accessible"
fi
fi
if test "$objs"
@@ -1531,7 +1521,7 @@
Marc shared lib : $progdll
Version type : $mode
Job ID : $DIRJID/$jid$extra_job_info
-User subroutine name : $usersubname
+User subroutine name : $user
User objects/libs : $objs
Restart file job ID : $rid
Substructure file ID : $sid
@@ -1564,7 +1554,7 @@
Marc shared lib : $progdll
Version type : $mode
Job ID : $DIRJID/$jid$extra_job_info
-User subroutine name : $usersubname
+User subroutine name : $user
User objects/libs : $objs
Restart file job ID : $rid
Substructure file ID : $sid
@@ -1687,7 +1677,7 @@
;;
esac
fi
- $ECHO "User subroutine name ($usersubname)? $ECHOTXT"
+ $ECHO "User subroutine name ($user)? $ECHOTXT"
read value
if test "$value"
then
@@ -1696,50 +1686,19 @@
user=
;;
*)
- user=`dirname $value`/`$BASENAME $value .f`
- usersubname=$user
- basefile=`$BASENAME $value`
- if test ${basefile##*.} = f
- then
- user=`dirname $value`/`$BASENAME $value .f`
- usersubname=$user.f
- elif test ${basefile##*.} = F
- then
- user=`dirname $value`/`$BASENAME $value .F`
- usersubname=$user.F
- elif test ${basefile##*.} = f90
- then
- user=`dirname $value`/`$BASENAME $value .f90`
- usersubname=$user.f90
- elif test ${basefile##*.} = F90
- then
- user=`dirname $value`/`$BASENAME $value .F90`
- usersubname=$user.F90
- fi
+ user=$value
case $user in
\/*)
;;
*)
user=`pwd`/$user
- usersubname=`pwd`/$usersubname
;;
esac
- if test ! -f $usersubname
- then
- if test -f $usersubname.f
- then
- usersubname=$usersubname.f
- elif test -f $usersubname.F
- then
- usersubname=$usersubname.F
- elif test -f $usersubname.f90
- then
- usersubname=$usersubname.f90
- elif test -f $usersubname.F90
- then
- usersubname=$usersubname.F90
- fi
- fi
+ usernoext=$user
+ usernoext=`dirname $usernoext`/`$BASENAME $usernoext .f`
+ usernoext=`dirname $usernoext`/`$BASENAME $usernoext .F`
+ usernoext=`dirname $usernoext`/`$BASENAME $usernoext .for`
+ usernoext=`dirname $usernoext`/`$BASENAME $usernoext .f90`
;;
esac
fi
@@ -2274,11 +2233,12 @@
#
# user subroutine used
#
+# add DAMASK options for linking
+ DAMASK="-lstdc++"
if test "$user"
then
-# program=$user.marc
- program=$DIRJOB/`$BASENAME $user .f`.marc
+ program=$usernoext.marc
case $program in
\/* | \.\/*)
bd=
@@ -2391,7 +2351,7 @@
fi
if test "$user"
then
- execpath=$DIRJOB/`$BASENAME $user .f`.marc
+ execpath=$usernoext.marc
usersub=1
fi
export execpath
@@ -3274,44 +3234,27 @@
echo
if test "$user"
then
- userobj=$DIRJOB/`$BASENAME $user .f`.o
- basefile=`$BASENAME $usersubname`
- if test ${basefile##*.} = f
- then
- usersub=$DIRJOB/`$BASENAME $user .f`.F
- ln -sf "$user.f" "$usersub"
- else
- usersub=$usersubname
- fi
-
+ userobj=$usernoext.o
fi
cat > $jid.runmarcscript << END4
if test "$user"
then
- if test ${basefile##*.} = f
- then
- ln -sf "$user.f" "$usersub"
- fi
if test $MACHINENAME = "CRAY"
then
- $FORTRAN $usersub || \
+ $DFORTLOWMP $user || \
{
- echo "$0: compile failed for $user.f"
+ echo "$0: compile failed for $user"
exit 1
}
/bin/rm $program 2>/dev/null
else
- $FORTRAN $usersub -o $userobj || \
+ $DFORTLOWMP $user -o $userobj || \
{
- echo "$0: compile failed for $user.f"
+ echo "$0: compile failed for $user"
exit 1
}
/bin/rm $program 2>/dev/null
fi
- if test ${basefile##*.} = f
- then
- /bin/rm -f "$usersub"
- fi
fi
@@ -3331,6 +3274,7 @@
$TKLIBS \
$MRCLIBS \
$METISLIBS \
+ $DAMASK \
$SFLIB \
$OPENSSL_LIB \
$SYSLIBS \
@@ -3344,6 +3288,9 @@
prgsav=yes
fi
/bin/rm $userobj 2>/dev/null
+/bin/rm $DIRJOB/*.mod 2>/dev/null
+/bin/rm $DIRJOB/*.smod 2>/dev/null
+/bin/rm $DIRJOB/*_genmod.f90 2>/dev/null
#
# run marc
@@ -3390,7 +3337,7 @@
fi
else
if test $cpdll = yes; then
- filename=`basename $usersubname .f`
+ filename=$usernoext
/bin/cp $DIRJOB/$marcdll $DIRJOB/${filename}_$marcdll 2>/dev/null
fi
if test $rmdll = yes
@@ -3556,7 +3503,7 @@
# first copy over the user sub if local directories
if test ${dirstatus[$counter]} = "local"
then
- $RCP $user.f $i:$DIR1/
+ $RCP $user $i:$DIR1/
fi
# do the compilation on the other machine
if test ${dirstatus[$counter]} = "shared"
@@ -3569,21 +3516,21 @@
remoteuser=$DIR1/`$BASENAME $user`
$RSH $i /bin/rm $remoteprog 2> /dev/null
echo
- $RSH $i $DIR2/tools/comp_user $DIR2 $DIR1 $remoteuser $remoteprog
+ $RSH $i $DIR2/tools/comp_damask_lmp $DIR2 $DIR1 $remoteuser $remoteprog
# check if successful, the new executable should be there
line=`$RSH $i /bin/ls $remoteprog 2> /dev/null`
if test "$line"
then
echo compilation and linking successful on host $i
else
- echo "$0: compile failed for $user.f on host $i"
+ echo "$0: compile failed for $user on host $i"
echo " $PRODUCT Exit number 3"
exit 1
fi
# remove the user subroutine on remote machine
if test ${dirstatus[$counter]} = "local"
then
- $RSH $i /bin/rm $remoteuser.f 2> /dev/null
+ $RSH $i /bin/rm $remoteuser 2> /dev/null
fi
fi
fi
@@ -3593,39 +3540,27 @@
if test "$userhost"
then
echo
- echo "Compiling and linking user subroutine $user.f on host `hostname`"
- fi
- userobj=$DIRJOB/`$BASENAME $user .f`.o
- basefile=`$BASENAME $usersubname`
- if test ${basefile##*.} = f
- then
- usersub=$DIRJOB/`$BASENAME $user .f`.F
- ln -sf "$user.f" "$usersub"
- else
- usersub=$usersubname
+ echo "Compiling and linking user subroutine $user on host `hostname`"
fi
+ userobj=$usernoext.o
if test $MACHINENAME = "CRAY"
then
- $FORTRAN $usersub || \
+ $DFORTLOWMP $user || \
{
- echo "$0: compile failed for $user.f"
+ echo "$0: compile failed for $user"
echo " $PRODUCT Exit number 3"
exit 1
}
/bin/rm $program 2>/dev/null
else
- $FORTRAN $usersub -o $userobj || \
+ $DFORTLOWMP $user -o $userobj || \
{
- echo "$0: compile failed for $user.f"
+ echo "$0: compile failed for $user"
echo " $PRODUCT Exit number 3"
exit 1
}
/bin/rm $program 2>/dev/null
fi
- if test ${basefile##*.} = f
- then
- /bin/rm -f "$usersub"
- fi
fi # if test $user
@@ -3645,6 +3580,7 @@
$TKLIBS \
$MRCLIBS \
$METISLIBS \
+ $DAMASK \
$SFLIB \
$OPENSSL_LIB \
$SYSLIBS \
@@ -3686,6 +3622,9 @@
prgsav=yes
fi # if test $link
/bin/rm $userobj 2>/dev/null
+/bin/rm $DIRJOB/*.mod 2>/dev/null
+/bin/rm $DIRJOB/*.smod 2>/dev/null
+/bin/rm $DIRJOB/*_genmod.f90 2>/dev/null
#
# run marc
@@ -3779,7 +3718,7 @@
else
#dllrun >0
if test $cpdll = yes; then
- filename=`basename $usersubname .f`
+ filename=$usernoext
/bin/cp $DIRJOB/$marcdll $DIRJOB/${filename}_$marcdll 2>/dev/null
fi
if test $rmdll = yes;then
@@ -3904,7 +3843,7 @@
# first copy over the user sub if local directories
if test ${dirstatus[$counter]} = "local"
then
- $RCP $user.f $i:$DIR1/
+ $RCP $user $i:$DIR1/
fi
# do the compilation on the other machine
if test ${dirstatus[$counter]} = "shared"
@@ -3917,20 +3856,20 @@
remoteuser=$DIR1/`$BASENAME $user`
$RSH $i /bin/rm $remoteprog 2> /dev/null
echo
- $RSH $i $DIR2/tools/comp_user $DIR2 $DIR1 $remoteuser $remoteprog
+ $RSH $i $DIR2/tools/comp_damask_lmp $DIR2 $DIR1 $remoteuser $remoteprog
# check if successful, the new executable should be there
line=`$RSH $i /bin/ls $remoteprog 2> /dev/null`
if test "$line"
then
echo compilation and linking successful on host $i
else
- echo "$0: compile failed for $user.f on host $i"
+ echo "$0: compile failed for $user on host $i"
exit 1
fi
# remove the user subroutine on remote machine
if test ${dirstatus[$counter]} = "local"
then
- $RSH $i /bin/rm $remoteuser.f 2> /dev/null
+ $RSH $i /bin/rm $remoteuser 2> /dev/null
fi
fi
fi
@@ -3940,37 +3879,25 @@
if test "$userhost"
then
echo
- echo "Compiling and linking user subroutine $user.f on host `hostname`"
- fi
- userobj=$DIRJOB/`$BASENAME $user .f`.o
- basefile=`$BASENAME $usersubname`
- if test ${basefile##*.} = f
- then
- usersub=$DIRJOB/`$BASENAME $user .f`.F
- ln -sf "$user.f" "$usersub"
- else
- usersub=$usersubname
+ echo "Compiling and linking user subroutine $user on host `hostname`"
fi
+ userobj=$usernoext.o
if test $MACHINENAME = "CRAY"
then
- $FORTRAN $usersub || \
+ $DFORTLOWMP $user || \
{
- echo "$0: compile failed for $user.f"
+ echo "$0: compile failed for $user"
exit 1
}
/bin/rm $program 2>/dev/null
else
- $FORTRAN $usersub -o $userobj || \
+ $DFORTLOWMP $user -o $userobj || \
{
- echo "$0: compile failed for $user.f"
+ echo "$0: compile failed for $user"
exit 1
}
/bin/rm $program 2>/dev/null
fi
- if test ${basefile##*.} = f
- then
- /bin/rm -f "$usersub"
- fi
fi # if test $user
@@ -3990,6 +3917,7 @@
$TKLIBS \
$MRCLIBS \
$METISLIBS \
+ $DAMASK \
$SFLIB \
$OPENSSL_LIB \
$SYSLIBS \
@@ -4030,7 +3958,9 @@
prgsav=yes
fi # if test $link
/bin/rm $userobj 2>/dev/null
-
+/bin/rm $DIRJOB/*.mod 2>/dev/null
+/bin/rm $DIRJOB/*.smod 2>/dev/null
+/bin/rm $DIRJOB/*_genmod.f90 2>/dev/null
# done if no job id given
if test -z "$jid"
then
@@ -4149,7 +4079,7 @@
else
#dllrun >0
if test $cpdll = yes; then
- filename=`basename $usersubname .f`
+ filename=$usernoext
/bin/cp $DIRJOB/$marcdll $DIRJOB/${filename}_$marcdll 2>/dev/null
fi
if test $rmdll = yes;then

517
install/MarcMentat/2022.2/Marc_tools/run_damask_mp.patch vendored Normal file
View File

@ -0,0 +1,517 @@
---
+++
@@ -136,6 +136,11 @@
# is created. For job running in the background, the log #
# file is always created. Default is "yes" #
##############################################################################
+# remove all Mentat paths from LD_LIBRARY_PATH
+LD_LIBRARY_PATH=:$LD_LIBRARY_PATH:
+LD_LIBRARY_PATH=${LD_LIBRARY_PATH//+([!(:)])mentat2022.2+([!(:)])/:}
+LD_LIBRARY_PATH=${LD_LIBRARY_PATH//+([(:)])/:}
+LD_LIBRARY_PATH=${LD_LIBRARY_PATH#:}; LD_LIBRARY_PATH=${LD_LIBRARY_PATH%:}
# set DIR to the directory in which this script is
REALCOM="`/bin/ls -l $0 |awk '{ print $NF; }'`"
DIR=`dirname $REALCOM`
@@ -302,7 +307,23 @@
. "$DIR/getarch"
+
+# getting user subroutine file name
+found=0
+for i in "$@"; do
+ if test $found = 1; then
+ DAMASK_USER=$i
+ found=0
+ fi
+ case $i in
+ -u* | -U*)
+ found=1
+ ;;
+ esac
+done
+# sourcing include_linux64 (needs DAMASK_USER to be set)
. $MARC_INCLUDE
+
#
#
@@ -405,7 +426,7 @@
did=
vid=
user=
-usersubname=
+usernoext=
objs=
qid=background
cpu=
@@ -676,50 +697,19 @@
esac
;;
-u* | -U*)
- user=`dirname $value`/`$BASENAME $value .f`
- usersubname=$user
- basefile=`$BASENAME $value`
- if test ${basefile##*.} = f
- then
- user=`dirname $value`/`$BASENAME $value .f`
- usersubname=$user.f
- elif test ${basefile##*.} = F
- then
- user=`dirname $value`/`$BASENAME $value .F`
- usersubname=$user.F
- elif test ${basefile##*.} = f90
- then
- user=`dirname $value`/`$BASENAME $value .f90`
- usersubname=$user.f90
- elif test ${basefile##*.} = F90
- then
- user=`dirname $value`/`$BASENAME $value .F90`
- usersubname=$user.F90
- fi
+ user=$value
case $user in
\/*)
;;
*)
user=`pwd`/$user
- usersubname=`pwd`/$usersubname
;;
esac
- if test ! -f $usersubname
- then
- if test -f $usersubname.f
- then
- usersubname=$usersubname.f
- elif test -f $usersubname.F
- then
- usersubname=$usersubname.F
- elif test -f $usersubname.f90
- then
- usersubname=$usersubname.f90
- elif test -f $usersubname.F90
- then
- usersubname=$usersubname.F90
- fi
- fi
+ usernoext=$user
+ usernoext=`dirname $usernoext`/`$BASENAME $usernoext .f`
+ usernoext=`dirname $usernoext`/`$BASENAME $usernoext .F`
+ usernoext=`dirname $usernoext`/`$BASENAME $usernoext .for`
+ usernoext=`dirname $usernoext`/`$BASENAME $usernoext .f90`
;;
-obj | -OBJ)
objs="$value"
@@ -1207,12 +1197,12 @@
fi
fi
fi
- if test "$usersubname"
+ if test "$user"
then
- if test ! -f $usersubname
+ if test ! -f $user
then
error="$error
-user subroutine file $usersubname not accessible"
+user subroutine file $user not accessible"
fi
fi
if test "$objs"
@@ -1531,7 +1521,7 @@
Marc shared lib : $progdll
Version type : $mode
Job ID : $DIRJID/$jid$extra_job_info
-User subroutine name : $usersubname
+User subroutine name : $user
User objects/libs : $objs
Restart file job ID : $rid
Substructure file ID : $sid
@@ -1564,7 +1554,7 @@
Marc shared lib : $progdll
Version type : $mode
Job ID : $DIRJID/$jid$extra_job_info
-User subroutine name : $usersubname
+User subroutine name : $user
User objects/libs : $objs
Restart file job ID : $rid
Substructure file ID : $sid
@@ -1687,7 +1677,7 @@
;;
esac
fi
- $ECHO "User subroutine name ($usersubname)? $ECHOTXT"
+ $ECHO "User subroutine name ($user)? $ECHOTXT"
read value
if test "$value"
then
@@ -1696,50 +1686,19 @@
user=
;;
*)
- user=`dirname $value`/`$BASENAME $value .f`
- usersubname=$user
- basefile=`$BASENAME $value`
- if test ${basefile##*.} = f
- then
- user=`dirname $value`/`$BASENAME $value .f`
- usersubname=$user.f
- elif test ${basefile##*.} = F
- then
- user=`dirname $value`/`$BASENAME $value .F`
- usersubname=$user.F
- elif test ${basefile##*.} = f90
- then
- user=`dirname $value`/`$BASENAME $value .f90`
- usersubname=$user.f90
- elif test ${basefile##*.} = F90
- then
- user=`dirname $value`/`$BASENAME $value .F90`
- usersubname=$user.F90
- fi
+ user=$value
case $user in
\/*)
;;
*)
user=`pwd`/$user
- usersubname=`pwd`/$usersubname
;;
esac
- if test ! -f $usersubname
- then
- if test -f $usersubname.f
- then
- usersubname=$usersubname.f
- elif test -f $usersubname.F
- then
- usersubname=$usersubname.F
- elif test -f $usersubname.f90
- then
- usersubname=$usersubname.f90
- elif test -f $usersubname.F90
- then
- usersubname=$usersubname.F90
- fi
- fi
+ usernoext=$user
+ usernoext=`dirname $usernoext`/`$BASENAME $usernoext .f`
+ usernoext=`dirname $usernoext`/`$BASENAME $usernoext .F`
+ usernoext=`dirname $usernoext`/`$BASENAME $usernoext .for`
+ usernoext=`dirname $usernoext`/`$BASENAME $usernoext .f90`
;;
esac
fi
@@ -2274,11 +2233,12 @@
#
# user subroutine used
#
+# add DAMASK options for linking
+ DAMASK="-lstdc++"
if test "$user"
then
-# program=$user.marc
- program=$DIRJOB/`$BASENAME $user .f`.marc
+ program=$usernoext.marc
case $program in
\/* | \.\/*)
bd=
@@ -2391,7 +2351,7 @@
fi
if test "$user"
then
- execpath=$DIRJOB/`$BASENAME $user .f`.marc
+ execpath=$usernoext.marc
usersub=1
fi
export execpath
@@ -3274,44 +3234,27 @@
echo
if test "$user"
then
- userobj=$DIRJOB/`$BASENAME $user .f`.o
- basefile=`$BASENAME $usersubname`
- if test ${basefile##*.} = f
- then
- usersub=$DIRJOB/`$BASENAME $user .f`.F
- ln -sf "$user.f" "$usersub"
- else
- usersub=$usersubname
- fi
-
+ userobj=$usernoext.o
fi
cat > $jid.runmarcscript << END4
if test "$user"
then
- if test ${basefile##*.} = f
- then
- ln -sf "$user.f" "$usersub"
- fi
if test $MACHINENAME = "CRAY"
then
- $FORTRAN $usersub || \
+ $DFORTRANMP $user || \
{
- echo "$0: compile failed for $user.f"
+ echo "$0: compile failed for $user"
exit 1
}
/bin/rm $program 2>/dev/null
else
- $FORTRAN $usersub -o $userobj || \
+ $DFORTRANMP $user -o $userobj || \
{
- echo "$0: compile failed for $user.f"
+ echo "$0: compile failed for $user"
exit 1
}
/bin/rm $program 2>/dev/null
fi
- if test ${basefile##*.} = f
- then
- /bin/rm -f "$usersub"
- fi
fi
@@ -3331,6 +3274,7 @@
$TKLIBS \
$MRCLIBS \
$METISLIBS \
+ $DAMASK \
$SFLIB \
$OPENSSL_LIB \
$SYSLIBS \
@@ -3344,6 +3288,9 @@
prgsav=yes
fi
/bin/rm $userobj 2>/dev/null
+/bin/rm $DIRJOB/*.mod 2>/dev/null
+/bin/rm $DIRJOB/*.smod 2>/dev/null
+/bin/rm $DIRJOB/*_genmod.f90 2>/dev/null
#
# run marc
@@ -3390,7 +3337,7 @@
fi
else
if test $cpdll = yes; then
- filename=`basename $usersubname .f`
+ filename=$usernoext
/bin/cp $DIRJOB/$marcdll $DIRJOB/${filename}_$marcdll 2>/dev/null
fi
if test $rmdll = yes
@@ -3556,7 +3503,7 @@
# first copy over the user sub if local directories
if test ${dirstatus[$counter]} = "local"
then
- $RCP $user.f $i:$DIR1/
+ $RCP $user $i:$DIR1/
fi
# do the compilation on the other machine
if test ${dirstatus[$counter]} = "shared"
@@ -3569,21 +3516,21 @@
remoteuser=$DIR1/`$BASENAME $user`
$RSH $i /bin/rm $remoteprog 2> /dev/null
echo
- $RSH $i $DIR2/tools/comp_user $DIR2 $DIR1 $remoteuser $remoteprog
+ $RSH $i $DIR2/tools/comp_damask_mp $DIR2 $DIR1 $remoteuser $remoteprog
# check if successful, the new executable should be there
line=`$RSH $i /bin/ls $remoteprog 2> /dev/null`
if test "$line"
then
echo compilation and linking successful on host $i
else
- echo "$0: compile failed for $user.f on host $i"
+ echo "$0: compile failed for $user on host $i"
echo " $PRODUCT Exit number 3"
exit 1
fi
# remove the user subroutine on remote machine
if test ${dirstatus[$counter]} = "local"
then
- $RSH $i /bin/rm $remoteuser.f 2> /dev/null
+ $RSH $i /bin/rm $remoteuser 2> /dev/null
fi
fi
fi
@@ -3593,39 +3540,27 @@
if test "$userhost"
then
echo
- echo "Compiling and linking user subroutine $user.f on host `hostname`"
- fi
- userobj=$DIRJOB/`$BASENAME $user .f`.o
- basefile=`$BASENAME $usersubname`
- if test ${basefile##*.} = f
- then
- usersub=$DIRJOB/`$BASENAME $user .f`.F
- ln -sf "$user.f" "$usersub"
- else
- usersub=$usersubname
+ echo "Compiling and linking user subroutine $user on host `hostname`"
fi
+ userobj=$usernoext.o
if test $MACHINENAME = "CRAY"
then
- $FORTRAN $usersub || \
+ $DFORTRANMP $user || \
{
- echo "$0: compile failed for $user.f"
+ echo "$0: compile failed for $user"
echo " $PRODUCT Exit number 3"
exit 1
}
/bin/rm $program 2>/dev/null
else
- $FORTRAN $usersub -o $userobj || \
+ $DFORTRANMP $user -o $userobj || \
{
- echo "$0: compile failed for $user.f"
+ echo "$0: compile failed for $user"
echo " $PRODUCT Exit number 3"
exit 1
}
/bin/rm $program 2>/dev/null
fi
- if test ${basefile##*.} = f
- then
- /bin/rm -f "$usersub"
- fi
fi # if test $user
@@ -3645,6 +3580,7 @@
$TKLIBS \
$MRCLIBS \
$METISLIBS \
+ $DAMASK \
$SFLIB \
$OPENSSL_LIB \
$SYSLIBS \
@@ -3686,6 +3622,9 @@
prgsav=yes
fi # if test $link
/bin/rm $userobj 2>/dev/null
+/bin/rm $DIRJOB/*.mod 2>/dev/null
+/bin/rm $DIRJOB/*.smod 2>/dev/null
+/bin/rm $DIRJOB/*_genmod.f90 2>/dev/null
#
# run marc
@@ -3779,7 +3718,7 @@
else
#dllrun >0
if test $cpdll = yes; then
- filename=`basename $usersubname .f`
+ filename=$usernoext
/bin/cp $DIRJOB/$marcdll $DIRJOB/${filename}_$marcdll 2>/dev/null
fi
if test $rmdll = yes;then
@@ -3904,7 +3843,7 @@
# first copy over the user sub if local directories
if test ${dirstatus[$counter]} = "local"
then
- $RCP $user.f $i:$DIR1/
+ $RCP $user $i:$DIR1/
fi
# do the compilation on the other machine
if test ${dirstatus[$counter]} = "shared"
@@ -3917,20 +3856,20 @@
remoteuser=$DIR1/`$BASENAME $user`
$RSH $i /bin/rm $remoteprog 2> /dev/null
echo
- $RSH $i $DIR2/tools/comp_user $DIR2 $DIR1 $remoteuser $remoteprog
+ $RSH $i $DIR2/tools/comp_damask_mp $DIR2 $DIR1 $remoteuser $remoteprog
# check if successful, the new executable should be there
line=`$RSH $i /bin/ls $remoteprog 2> /dev/null`
if test "$line"
then
echo compilation and linking successful on host $i
else
- echo "$0: compile failed for $user.f on host $i"
+ echo "$0: compile failed for $user on host $i"
exit 1
fi
# remove the user subroutine on remote machine
if test ${dirstatus[$counter]} = "local"
then
- $RSH $i /bin/rm $remoteuser.f 2> /dev/null
+ $RSH $i /bin/rm $remoteuser 2> /dev/null
fi
fi
fi
@@ -3940,37 +3879,25 @@
if test "$userhost"
then
echo
- echo "Compiling and linking user subroutine $user.f on host `hostname`"
- fi
- userobj=$DIRJOB/`$BASENAME $user .f`.o
- basefile=`$BASENAME $usersubname`
- if test ${basefile##*.} = f
- then
- usersub=$DIRJOB/`$BASENAME $user .f`.F
- ln -sf "$user.f" "$usersub"
- else
- usersub=$usersubname
+ echo "Compiling and linking user subroutine $user on host `hostname`"
fi
+ userobj=$usernoext.o
if test $MACHINENAME = "CRAY"
then
- $FORTRAN $usersub || \
+ $DFORTRANMP $user || \
{
- echo "$0: compile failed for $user.f"
+ echo "$0: compile failed for $user"
exit 1
}
/bin/rm $program 2>/dev/null
else
- $FORTRAN $usersub -o $userobj || \
+ $DFORTRANMP $user -o $userobj || \
{
- echo "$0: compile failed for $user.f"
+ echo "$0: compile failed for $user"
exit 1
}
/bin/rm $program 2>/dev/null
fi
- if test ${basefile##*.} = f
- then
- /bin/rm -f "$usersub"
- fi
fi # if test $user
@@ -3990,6 +3917,7 @@
$TKLIBS \
$MRCLIBS \
$METISLIBS \
+ $DAMASK \
$SFLIB \
$OPENSSL_LIB \
$SYSLIBS \
@@ -4030,7 +3958,9 @@
prgsav=yes
fi # if test $link
/bin/rm $userobj 2>/dev/null
-
+/bin/rm $DIRJOB/*.mod 2>/dev/null
+/bin/rm $DIRJOB/*.smod 2>/dev/null
+/bin/rm $DIRJOB/*_genmod.f90 2>/dev/null
# done if no job id given
if test -z "$jid"
then
@@ -4149,7 +4079,7 @@
else
#dllrun >0
if test $cpdll = yes; then
- filename=`basename $usersubname .f`
+ filename=$usernoext
/bin/cp $DIRJOB/$marcdll $DIRJOB/${filename}_$marcdll 2>/dev/null
fi
if test $rmdll = yes;then

24
install/MarcMentat/2022.2/Mentat_bin/edit_window.patch vendored Normal file
View File

@ -0,0 +1,24 @@
---
+++
@@ -1,18 +1,5 @@
#!/bin/sh
-# This script opens a window running an editor. The default window is an
-# xterm, and the default editor is vi. These may be customized.
+# This script opens a window running an editor.
+# The command to invoke the editor is specified during DAMASK installation
-dir=
-for d in /usr/bin /usr/bin/X11; do
- if test -x "$d/xterm"; then
- dir="$d"
- break
- fi
-done
-
-if test -z "$dir"; then
- echo "$0: Could not find xterm"
- exit 1
-fi
-
-"$dir/xterm" -T "vi $*" -n "vi $*" -e vi $*
+%EDITOR% $*

0
install/MarcMentat/2022.2/Mentat_bin/kill4.patch vendored Normal file
View File

0
install/MarcMentat/2022.2/Mentat_bin/kill5.patch vendored Normal file
View File

0
install/MarcMentat/2022.2/Mentat_bin/kill6.patch vendored Normal file
View File

38
install/MarcMentat/2022.2/Mentat_bin/submit4.patch vendored Normal file
View File

@ -0,0 +1,38 @@
---
+++
@@ -63,10 +63,10 @@
if [ "$slv" != "" -a "$slv" != "marc" -a "$slv" != "datfit" ]; then
slv="-iam sfm"
fi
-if [ "$slv" == "marc" ]; then
+if [ "$slv" = "marc" ]; then
slv=""
fi
-if [ "$slv" == "datfit" ]; then
+if [ "$slv" = "datfit" ]; then
slv="-iam datfit"
fi
@@ -91,6 +91,7 @@
srcfile="-u $srcfile -save y"
;;
runsaved)
+ srcfile=${srcfile%.*}".marc"
srcfile="-prog $srcfile"
;;
esac
@@ -189,12 +190,12 @@
unset PYTHONPATH
if [ "$doe_first" = "-" ]; then # submit of regular Marc job
- "${DIR}/tools/run_marc" $slv -j $job -v n -b y $nprocds $nprocd \
+ "${DIR}/tools/run_damask_hmp" $slv -j $job -v n -b y $nprocds $nprocd \
$srcfile $restart $postfile $viewfactorsfile $hostfile \
$compat $copy_datfile $copy_postfile $scr_dir $dcoup \
$assem_recov_nthread $nthread $nsolver $mode $gpu > /dev/null 2>&1
else # submit of a DoE Marc job
- "${DIR}/tools/run_marc" $slv -j $job -v n -b n $nprocds $nprocd \
+ "${DIR}/tools/run_damask_hmp" $slv -j $job -v n -b n $nprocds $nprocd \
$srcfile $restart $postfile $viewfactorsfile $hostfile \
$compat $copy_datfile $copy_postfile $scr_dir $dcoup \
$assem_recov_nthread $nthread $nsolver $mode $gpu

38
install/MarcMentat/2022.2/Mentat_bin/submit5.patch vendored Normal file
View File

@ -0,0 +1,38 @@
---
+++
@@ -63,10 +63,10 @@
if [ "$slv" != "" -a "$slv" != "marc" -a "$slv" != "datfit" ]; then
slv="-iam sfm"
fi
-if [ "$slv" == "marc" ]; then
+if [ "$slv" = "marc" ]; then
slv=""
fi
-if [ "$slv" == "datfit" ]; then
+if [ "$slv" = "datfit" ]; then
slv="-iam datfit"
fi
@@ -91,6 +91,7 @@
srcfile="-u $srcfile -save y"
;;
runsaved)
+ srcfile=${srcfile%.*}".marc"
srcfile="-prog $srcfile"
;;
esac
@@ -189,12 +190,12 @@
unset PYTHONPATH
if [ "$doe_first" = "-" ]; then # submit of regular Marc job
- "${DIR}/tools/run_marc" $slv -j $job -v n -b y $nprocds $nprocd \
+ "${DIR}/tools/run_damask_mp" $slv -j $job -v n -b y $nprocds $nprocd \
$srcfile $restart $postfile $viewfactorsfile $hostfile \
$compat $copy_datfile $copy_postfile $scr_dir $dcoup \
$assem_recov_nthread $nthread $nsolver $mode $gpu > /dev/null 2>&1
else # submit of a DoE Marc job
- "${DIR}/tools/run_marc" $slv -j $job -v n -b n $nprocds $nprocd \
+ "${DIR}/tools/run_damask_mp" $slv -j $job -v n -b n $nprocds $nprocd \
$srcfile $restart $postfile $viewfactorsfile $hostfile \
$compat $copy_datfile $copy_postfile $scr_dir $dcoup \
$assem_recov_nthread $nthread $nsolver $mode $gpu

38
install/MarcMentat/2022.2/Mentat_bin/submit6.patch vendored Normal file
View File

@ -0,0 +1,38 @@
---
+++
@@ -63,10 +63,10 @@
if [ "$slv" != "" -a "$slv" != "marc" -a "$slv" != "datfit" ]; then
slv="-iam sfm"
fi
-if [ "$slv" == "marc" ]; then
+if [ "$slv" = "marc" ]; then
slv=""
fi
-if [ "$slv" == "datfit" ]; then
+if [ "$slv" = "datfit" ]; then
slv="-iam datfit"
fi
@@ -91,6 +91,7 @@
srcfile="-u $srcfile -save y"
;;
runsaved)
+ srcfile=${srcfile%.*}".marc"
srcfile="-prog $srcfile"
;;
esac
@@ -189,12 +190,12 @@
unset PYTHONPATH
if [ "$doe_first" = "-" ]; then # submit of regular Marc job
- "${DIR}/tools/run_marc" $slv -j $job -v n -b y $nprocds $nprocd \
+ "${DIR}/tools/run_damask_lmp" $slv -j $job -v n -b y $nprocds $nprocd \
$srcfile $restart $postfile $viewfactorsfile $hostfile \
$compat $copy_datfile $copy_postfile $scr_dir $dcoup \
$assem_recov_nthread $nthread $nsolver $mode $gpu > /dev/null 2>&1
else # submit of a DoE Marc job
- "${DIR}/tools/run_marc" $slv -j $job -v n -b n $nprocds $nprocd \
+ "${DIR}/tools/run_damask_lmp" $slv -j $job -v n -b n $nprocds $nprocd \
$srcfile $restart $postfile $viewfactorsfile $hostfile \
$compat $copy_datfile $copy_postfile $scr_dir $dcoup \
$assem_recov_nthread $nthread $nsolver $mode $gpu

158
install/MarcMentat/2022.2/Mentat_menus/job_run.ms.patch vendored Normal file
View File

@ -0,0 +1,158 @@
---
+++
@@ -261,11 +261,18 @@
}
button {
position +25 =
- size 25 4
+ size 18 4
text "ADVANCED JOB SUBMISSION"
help "job_run#Job Submission And Control"
popmenu job_submit_adv_pm
}
+ button {
+ position +18 =
+ size 7 4
+ text "DAMASK"
+ help "damask_run#Job Submission And Control"
+ popmenu damask
+ }
button {
position 0 +4
size 16 4
@@ -1202,6 +1209,135 @@
}
+#--------------------------------------------------------------------------------------------------
+popmenu damask {
+
+#ifdef QT_MENTAT
+ text "DAMASK.MPIE.DE"
+#endif
+
+ group {
+#ifndef QT_MENTAT
+ label {
+ position 0 0
+ size 50 4
+ text "DAMASK.MPIE.DE"
+ }
+#endif
+
+ label {
+ position 1 6
+ size 13 6
+ text "Optimzation"
+ border_width 1
+ border_color black
+ }
+
+ label {
+ position +13 =
+ size 20 6
+ text "write Input"
+ border_width 1
+ border_color black
+ }
+
+ label {
+ position +18 =
+ size 30 6
+ text "do not write Inp."
+ border_width 1
+ border_color black
+ }
+
+ label {
+ position -32 +6
+ size 12 6
+ text "O3 / OpenMP"
+ border_width 1
+ border_color black
+ }
+
+ popdown {
+ position +12 =
+ size 20 6
+ text "Submit"
+ command "*submit_job 4 *monitor_job"
+ }
+
+ popdown {
+ position +20 =
+ size 20 6
+ text "Execute"
+ command "*execute_job 4 *monitor_job"
+ }
+
+ label {
+ position -32 +6
+ size 12 6
+ text "O1 / OpenMP"
+ border_width 1
+ border_color black
+ }
+
+ popdown {
+ position +12 =
+ size 20 6
+ text "Submit"
+ command "*submit_job 5 *monitor_job"
+ }
+
+ popdown {
+ position +20 =
+ size 20 6
+ text "Execute"
+ command "*execute_job 5 *monitor_job"
+ }
+
+ label {
+ position -32 +6
+ size 12 6
+ text "O0 / OpenMP"
+ border_width 1
+ border_color black
+ }
+
+ popdown {
+ position +12 =
+ size 20 6
+ text "Submit"
+ command "*submit_job 6 *monitor_job"
+ }
+
+ popdown {
+ position +20 =
+ size 20 6
+ text "Execute"
+ command "*execute_job 6 *monitor_job"
+ }
+
+ popdown {
+ position 19 +8
+ size 12 8
+ text "CANCEL"
+ }
+}
+
+ window {
+ parent mentat
+ origin 38 8
+#ifdef DCOM
+ size 50 100
+#else
+ size 50 94
+#endif
+ background_color body
+ border_width 1
+ border_color border
+ buffering single
+ }
+ mode permanent
+}
+
#--------------------------------------------------------------------------------------------------
popmenu job_exit_msg_pm {

47
install/MarcMentat/apply_DAMASK_modifications.py → install/MarcMentat/MSC_modifications.py vendored
View File

@ -9,9 +9,12 @@ from pathlib import Path
import subprocess
import shlex
sys.path.append(str(Path(__file__).parents[2]/'python/damask'))
sys.path.append(str(Path(__file__).resolve().parents[2]/'python/damask'))
import solver
APPLY = 'install'
RESTORE = 'uninstall'
def copy_and_patch(patch,orig,editor):
try:
shutil.copyfile(orig,orig.parent/patch.stem)
@ -25,9 +28,11 @@ def copy_and_patch(patch,orig,editor):
parser = argparse.ArgumentParser(
description='Apply DAMASK modification to MSC Marc/Mentat',
description=f'{APPLY.capitalize()} or {RESTORE} DAMASK modifications to MSC Marc/Mentat',
formatter_class=argparse.ArgumentDefaultsHelpFormatter)
parser.add_argument('command', metavar='command', nargs='?', default=APPLY, choices=[APPLY,RESTORE],
help=f'Mode of operation {[APPLY,RESTORE]}')
parser.add_argument('--editor', dest='editor', metavar='string', default='vi',
help='Name of the editor (executable) used by Marc Mentat')
parser.add_argument('--marc-root', dest='marc_root', metavar='string',
@ -40,10 +45,10 @@ parser.add_argument('--damask-root', dest='damask_root', metavar = 'string',
default=solver._marc._damask_root,
help='DAMASK root directory')
args = parser.parse_args()
marc_root = Path(args.marc_root).expanduser()
damask_root = Path(args.damask_root).expanduser()
marc_root = Path(args.marc_root).expanduser()
marc_version = args.marc_version
matches = {'Marc_tools': [['comp_user','comp_damask_*mp'],
@ -54,23 +59,41 @@ matches = {'Marc_tools': [['comp_user','comp_damask_*mp'],
['kill1','kill?']],
'Mentat_menus':[['job_run.ms','job_run.ms']]}
for cmd in ['patch','xvfb-run']:
for cmd in ['xvfb-run','patch'] if args.command == APPLY else ['xvfb-run'] if args.command == RESTORE else []:
try:
subprocess.run(shlex.split(f'{cmd} --help'))
subprocess.run([cmd,'--help'], capture_output=True)
except FileNotFoundError:
print(f'"{cmd}" not found, please install')
sys.exit()
if args.command == APPLY:
print('patching files...')
for directory in glob.glob(str(damask_root/'install/MarcMentat'/marc_version/'*')):
for orig, mods in matches[Path(directory).name]:
product,subfolder = (marc_root/Path(directory)).name.split('_')
product,subfolder = (marc_root/directory).name.split('_')
orig = marc_root/f'{product.lower()}{marc_version}/{subfolder}/{orig}'
for patch in glob.glob(f'{directory}/{mods}.patch'):
copy_and_patch(Path(patch),orig,args.editor)
elif args.command == RESTORE:
print('deleting modified files...')
for file in (glob.glob(str(marc_root/f'marc{marc_version}/tools/*_damask*')) +
glob.glob(str(marc_root/f'mentat{marc_version}/bin/kill[4-6]')) +
glob.glob(str(marc_root/f'mentat{marc_version}/bin/submit[4-6]*'))):
os.remove(file)
print('restoring original files...')
for file in (glob.glob(str(marc_root/f'marc{marc_version}/tools/include_linux64.orig')) +
glob.glob(str(marc_root/f'mentat{marc_version}/bin/edit_window.orig')) +
glob.glob(str(marc_root/f'mentat{marc_version}/menus/job_run.ms.orig'))):
shutil.copyfile(file,Path(file).with_suffix(''))
os.remove(file)
else:
print('skipping file modifications...')
print('compiling Mentat menu binaries...')
executable = marc_root/f'mentat{marc_version}/bin/mentat'
@ -79,8 +102,10 @@ subprocess.run(shlex.split(f'xvfb-run -a {executable} -compile {menu_file}'))
print('setting file access rights...')
for file in (glob.glob(str(marc_root/f'marc{marc_version}/tools/*_damask*')) +
glob.glob(str(marc_root/f'mentat{marc_version}/bin/edit_window')) +
for file in (
(glob.glob(str(marc_root/f'marc{marc_version}/tools/*_damask*')) +
glob.glob(str(marc_root/f'mentat{marc_version}/bin/kill[4-6]')) +
glob.glob(str(marc_root/f'mentat{marc_version}/bin/submit[4-6]'))):
glob.glob(str(marc_root/f'mentat{marc_version}/bin/submit[4-6]')) if args.command == APPLY else []) +
glob.glob(str(marc_root/f'mentat{marc_version}/bin/edit_window'))
):
os.chmod(file , 0o755)

15
install/MarcMentat/installation.txt vendored
View File

@ -1,7 +1,12 @@
Install DAMASK modifications to use DAMASK_marc
This is for the Linux64 version of Marc/Mentat
Install or uninstall DAMASK modifications to MSC Marc/Mentat.
This is for the Linux64 version of Marc/Mentat.
Refer to http://damask.mpie.de for complete installation instructions.
usage: MSC_modifications.py [-h] [--editor string] [--marc-root string] [--marc-version string] [--damask-root string]
[command]
where command is either 'install' or 'uninstall'
Refer to https://damask.mpie.de/installation/source_code.html#msc-marc for complete installation instructions.
Usually you will need to be root for this to work!
@ -9,6 +14,6 @@ See Marc and Mentat Release Guide for List of Build and Supported Platforms.
The Intel Fortran compiler needs to be installed.
1) Install Marc, Mentat and Documentation as usual
1) Install Marc, Mentat, and Documentation as usual.
Run the test example including subroutines to confirm that the installation of both Marc/Mentat and the Intel Fortran Compiler is ok.
2) Run the apply_DAMASK_modifications.py script from this directory.
2) Run the MSC_modifications.py script from this directory.

25
python/damask/_colormap.py
View File

@ -2,7 +2,8 @@ import os
import json
import functools
import colorsys
from typing import Union, TextIO
from typing import Optional, Union, TextIO
from itertools import chain
import numpy as np
import scipy.interpolate as interp
@ -248,7 +249,7 @@ class Colormap(mpl.colors.ListedColormap):
Parameters
----------
fraction : float or sequence of float
fraction : (sequence of) float
Fractional coordinate(s) to evaluate Colormap at.
Returns
@ -274,8 +275,8 @@ class Colormap(mpl.colors.ListedColormap):
def shade(self,
field: np.ndarray,
bounds: FloatSequence = None,
gap: float = None) -> Image:
bounds: Optional[FloatSequence] = None,
gap: Optional[float] = None) -> Image:
"""
Generate PIL image of 2D field using colormap.
@ -314,7 +315,7 @@ class Colormap(mpl.colors.ListedColormap):
def reversed(self,
name: str = None) -> 'Colormap':
name: Optional[str] = None) -> 'Colormap':
"""
Reverse.
@ -363,7 +364,7 @@ class Colormap(mpl.colors.ListedColormap):
def save_paraview(self,
fname: FileHandle = None):
fname: Optional[FileHandle] = None):
"""
Save as JSON file for use in Paraview.
@ -373,16 +374,12 @@ class Colormap(mpl.colors.ListedColormap):
File to store results. Defaults to colormap name + '.json'.
"""
colors = []
for i,c in enumerate(np.round(self.colors,6).tolist()):
colors+=[i]+c
out = [{
'Creator':util.execution_stamp('Colormap'),
'ColorSpace':'RGB',
'Name':self.name,
'DefaultMap':True,
'RGBPoints':colors
'RGBPoints':list(chain.from_iterable([(i,*c) for i,c in enumerate(self.colors.round(6))]))
}]
fhandle = self._get_file_handle(fname,'.json')
@ -391,7 +388,7 @@ class Colormap(mpl.colors.ListedColormap):
def save_ASCII(self,
fname: FileHandle = None):
fname: Optional[FileHandle] = None):
"""
Save as ASCII file.
@ -406,7 +403,7 @@ class Colormap(mpl.colors.ListedColormap):
t.save(self._get_file_handle(fname,'.txt'))
def save_GOM(self, fname: FileHandle = None):
def save_GOM(self, fname: Optional[FileHandle] = None):
"""
Save as ASCII file for use in GOM Aramis.
@ -427,7 +424,7 @@ class Colormap(mpl.colors.ListedColormap):
def save_gmsh(self,
fname: FileHandle = None):
fname: Optional[FileHandle] = None):
"""
Save as ASCII file for use in gmsh.

6
python/damask/_config.py
View File

@ -2,7 +2,7 @@ import copy
from io import StringIO
from collections.abc import Iterable
import abc
from typing import Union, Dict, Any, Type, TypeVar
from typing import Optional, Union, Dict, Any, Type, TypeVar
import numpy as np
import yaml
@ -21,7 +21,7 @@ class NiceDumper(yaml.SafeDumper):
"""Make YAML readable for humans."""
def write_line_break(self,
data: str = None):
data: Optional[str] = None):
super().write_line_break(data)
if len(self.indents) == 1:
@ -53,7 +53,7 @@ class Config(dict):
"""YAML-based configuration."""
def __init__(self,
yml: Union[str, Dict[str, Any]] = None,
yml: Union[None, str, Dict[str, Any]] = None,
**kwargs):
"""Initialize from YAML, dict, or key=value pairs."""
if isinstance(yml,str):

190
python/damask/_configmaterial.py
View File

@ -1,6 +1,6 @@
import numpy as np
import h5py
from typing import Sequence, Dict, Any, Collection
from typing import Optional, Union, Sequence, Dict, Any, Collection
from ._typehints import FileHandle
from . import Config
@ -22,7 +22,7 @@ class ConfigMaterial(Config):
"""
def __init__(self,
d: Dict[str, Any] = None,
d: Optional[Dict[str, Any]] = None,
**kwargs):
"""
New material configuration.
@ -83,13 +83,13 @@ class ConfigMaterial(Config):
@staticmethod
def load_DREAM3D(fname: str,
grain_data: str = None,
cell_data: str = None,
grain_data: Optional[str] = None,
cell_data: Optional[str] = None,
cell_ensemble_data: str = 'CellEnsembleData',
phases: str = 'Phases',
Euler_angles: str = 'EulerAngles',
phase_names: str = 'PhaseName',
base_group: str = None) -> 'ConfigMaterial':
base_group: Optional[str] = None) -> 'ConfigMaterial':
"""
Load DREAM.3D (HDF5) file.
@ -160,7 +160,7 @@ class ConfigMaterial(Config):
pass
base_config = ConfigMaterial({'phase':{k if isinstance(k,int) else str(k):'t.b.d.' for k in np.unique(phase)},
base_config = ConfigMaterial({'phase':{k if isinstance(k,int) else str(k): None for k in np.unique(phase)},
'homogenization':{'direct':{'N_constituents':1}}})
constituent = {k:np.atleast_1d(v[idx].squeeze()) for k,v in zip(['O','phase'],[O,phase])}
@ -193,10 +193,10 @@ class ConfigMaterial(Config):
>>> import damask.ConfigMaterial as cm
>>> t = damask.Table.load('small.txt')
>>> t
pos pos pos qu qu qu qu phase homog
3:pos pos pos 4:qu qu qu qu phase homog
0 0 0 0 0.19 0.8 0.24 -0.51 Aluminum SX
1 1 0 0 0.8 0.19 0.24 -0.51 Steel SX
1 1 1 0 0.8 0.19 0.24 -0.51 Steel SX
2 1 1 0 0.8 0.19 0.24 -0.51 Steel SX
>>> cm.from_table(t,O='qu',phase='phase',homogenization='homog')
material:
- constituents:
@ -209,8 +209,8 @@ class ConfigMaterial(Config):
v: 1.0
phase: Steel
homogenization: SX
homogenization: {}
phase: {}
homogenization: {SX: null}
phase: {Aluminum: null, Steel: null}
>>> cm.from_table(t,O='qu',phase='phase',homogenization='single_crystal')
material:
@ -224,8 +224,8 @@ class ConfigMaterial(Config):
v: 1.0
phase: Steel
homogenization: single_crystal
homogenization: {}
phase: {}
homogenization: {single_crystal: null}
phase: {Aluminum: null, Steel: null}
"""
kwargs_ = {k:table.get(v) if v in table.labels else np.atleast_2d([v]*len(table)).T for k,v in kwargs.items()}
@ -233,6 +233,8 @@ class ConfigMaterial(Config):
_,idx = np.unique(np.hstack(list(kwargs_.values())),return_index=True,axis=0)
idx = np.sort(idx)
kwargs_ = {k:np.atleast_1d(v[idx].squeeze()) for k,v in kwargs_.items()}
for what in ['phase','homogenization']:
if what not in kwargs_: kwargs_[what] = what+'_label'
return ConfigMaterial().material_add(**kwargs_)
@ -243,7 +245,7 @@ class ConfigMaterial(Config):
Check for completeness.
Only the general file layout is considered.
This check does not consider whether parameters for
This check does not consider whether specific parameters for
a particular phase/homogenization model are missing.
Returns
@ -252,52 +254,56 @@ class ConfigMaterial(Config):
Whether the material.yaml definition is complete.
"""
def LabeledList(label,items):
return f'{label.capitalize()}{"s" if len(items)>1 else ""} {util.srepr(items,",",quote=True)}'
ok = True
for top_level in ['homogenization','phase','material']:
ok &= top_level in self
if top_level not in self: print(f'{top_level} entry missing')
msg = []
all = set(['homogenization','phase','material'])
miss = set([item for item in all if item not in self])
empty = set([item for item in all-miss if self[item] is None])
if miss:
msg.append(f'{LabeledList("top-level",miss)} missing')
ok = False
if empty:
msg.append(f'{LabeledList("top-level",empty)} empty')
if ok:
ok &= len(self['material']) > 0
if len(self['material']) < 1: print('Incomplete material definition')
if len(self['material']) < 1: msg.append('No materials defined')
if ok:
homogenization = set()
phase = set()
for i,v in enumerate(self['material']):
if 'homogenization' in v:
homogenization.add(v['homogenization'])
else:
print(f'No homogenization specified in material {i}')
msg.append(f'No homogenization specified for material {i}')
ok = False
if 'constituents' in v:
for ii,vv in enumerate(v['constituents']):
if 'O' not in vv:
print('No orientation specified in constituent {ii} of material {i}')
msg.append(f'No orientation specified for constituent {ii} of material {i}')
ok = False
if 'phase' in vv:
phase.add(vv['phase'])
else:
print(f'No phase specified in constituent {ii} of material {i}')
msg.append(f'No phase specified for constituent {ii} of material {i}')
ok = False
for k,v in self['phase'].items():
if 'lattice' not in v:
print(f'No lattice specified in phase {k}')
for v,other in {'phase':phase,
'homogenization':homogenization}.items():
me = set([] if v in empty else self[v])
if _miss := other - me:
msg.append(f'{LabeledList(v,_miss)} missing')
ok = False
if len(_empty := [item for item in me if self[v][item] is None]) > 0:
msg.append(f'{LabeledList(v,_empty)} undefined')
ok = False
for k,v in self['homogenization'].items():
if 'N_constituents' not in v:
print(f'No. of constituents not specified in homogenization {k}')
ok = False
if phase - set(self['phase']):
print(f'Phase(s) {phase-set(self["phase"])} missing')
ok = False
if homogenization - set(self['homogenization']):
print(f'Homogenization(s) {homogenization-set(self["homogenization"])} missing')
ok = False
print(util.srepr(msg))
return ok
@ -320,7 +326,7 @@ class ConfigMaterial(Config):
if 'phase' in self:
for k,v in self['phase'].items():
if 'lattice' in v:
if v is not None and 'lattice' in v:
try:
Orientation(lattice=v['lattice'])
except KeyError:
@ -348,8 +354,8 @@ class ConfigMaterial(Config):
def material_rename_phase(self,
mapping: Dict[str, str],
ID: Sequence[int] = None,
constituent: Sequence[int] = None) -> 'ConfigMaterial':
ID: Optional[Sequence[int]] = None,
constituent: Optional[Sequence[int]] = None) -> 'ConfigMaterial':
"""
Change phase name in material.
@ -382,7 +388,7 @@ class ConfigMaterial(Config):
def material_rename_homogenization(self,
mapping: Dict[str, str],
ID: Sequence[int] = None) -> 'ConfigMaterial':
ID: Optional[Sequence[int]] = None) -> 'ConfigMaterial':
"""
Change homogenization name in material.
@ -418,6 +424,8 @@ class ConfigMaterial(Config):
----------
**kwargs
Key-value pairs.
First index of array-like values runs over materials,
whereas second index runs over constituents.
Returns
-------
@ -426,13 +434,12 @@ class ConfigMaterial(Config):
Examples
--------
Create a dual-phase steel microstructure for micromechanical simulations:
Create two grains of ferrite and one grain of martensite, each with random orientation:
>>> import numpy as np
>>> import damask
>>> m = damask.ConfigMaterial()
>>> m = m.material_add(phase = ['Ferrite','Martensite'],
... O = damask.Rotation.from_random(2),
>>> m = m.material_add(phase = ['Ferrite','Martensite','Ferrite'],
... O = damask.Rotation.from_random(3),
... homogenization = 'SX')
>>> m
material:
@ -446,60 +453,91 @@ class ConfigMaterial(Config):
v: 1.0
phase: Martensite
homogenization: SX
homogenization: {}
phase: {}
- constituents:
- O: [0.47925185, -0.04294454, 0.78760173, -0.3849116 ]
v: 1.0
phase: Ferrite
homogenization: SX
homogenization: {SX: null}
phase: {Ferrite: null, Martensite: null}
Create a duplex stainless steel microstructure for forming simulations:
Create hundred materials that each approximate a duplex stainless steel microstructure
with three austenite and one relatively bigger ferrite grain of random orientation each:
>>> import numpy as np
>>> import damask
>>> m = damask.ConfigMaterial()
>>> m = m.material_add(phase = np.array(['Austenite','Ferrite']).reshape(1,2),
... O = damask.Rotation.from_random((2,2)),
... v = np.array([0.2,0.8]).reshape(1,2),
>>> m = m.material_add(phase = np.array(['Austenite']*3+['Ferrite']),
... O = damask.Rotation.from_random((100,4)),
... v = np.array([0.2]*3+[0.4]),
... homogenization = 'Taylor')
>>> m
material:
- constituents:
- phase: Austenite
O: [0.659802978293224, 0.6953785848195171, 0.22426295326327111, -0.17554139512785227]
v: 0.2
- phase: Ferrite
O: [0.49356745891301596, 0.2841806579193434, -0.7487679215072818, -0.339085707289975]
v: 0.8
- v: 0.2
phase: Austenite
O: [0.46183665006602664, 0.2215160420973196, -0.5594313187331139, 0.6516702781083836]
- v: 0.2
phase: Austenite
O: [0.11321658382410027, 0.6354079414360444, 0.00562701344273936, 0.7638108992590535]
- v: 0.2
phase: Austenite
O: [0.050991978809077604, 0.8069522034362003, -0.11352928955610851, -0.5773552285027659]
- v: 0.4
phase: Ferrite
O: [0.9460076150721788, 0.15880754622367604, -0.0069841062241482385, -0.28249066842661014]
homogenization: Taylor
.
.
.
- constituents:
- phase: Austenite
O: [0.26542221365204055, 0.7268854930702071, 0.4474726435701472, -0.44828201137283735]
v: 0.2
- phase: Ferrite
O: [0.6545817158479885, -0.08004812803625233, -0.6226561293931374, 0.4212059104577611]
v: 0.8
- v: 0.2
phase: Austenite
O: [0.12531400788494199, -0.18637769037997565, 0.31737548053338394, -0.9213210951197429]
- v: 0.2
phase: Austenite
O: [0.37453930577161404, -0.33529507696450805, -0.3266564259130028, -0.800370601162502]
- v: 0.2
phase: Austenite
O: [0.035776891752713764, -0.720706371010592, -0.4540438656728926, -0.5226342017569017]
- v: 0.4
phase: Ferrite
O: [0.6782596727966124, -0.20800082041703685, -0.138636083554039, 0.6909989227925536]
homogenization: Taylor
homogenization: {}
phase: {}
homogenization: {Taylor: null}
phase: {Austenite: null, Ferrite: null}
"""
N,n,shaped = 1,1,{}
_constituent_properties = ['phase','O','v','V_e']
_dim = {'O':(4,),'V_e':(3,3,)}
_ex = dict((k, -len(v)) for k, v in _dim.items())
N,n = 1,1
shaped : Dict[str, Union[None,np.ndarray]] = \
{'v': None,
'phase': None,
'homogenization': None,
}
map_dim = {'O':-1,'V_e':-2}
for k,v in kwargs.items():
shaped[k] = np.array(v)
s = shaped[k].shape[:map_dim.get(k,None)]
s = shaped[k].shape[:_ex.get(k,None)] # type: ignore
N = max(N,s[0]) if len(s)>0 else N
n = max(n,s[1]) if len(s)>1 else n
shaped['v'] = np.array(1./n) if shaped['v'] is None else shaped['v']
mat: Sequence[dict] = [{'constituents':[{} for _ in range(n)]} for _ in range(N)]
if 'v' not in kwargs:
shaped['v'] = np.broadcast_to(1/n,(N,n))
map_shape = {'O':(N,n,4),'V_e':(N,n,3,3)}
for k,v in shaped.items():
target = map_shape.get(k,(N,n))
obj = np.broadcast_to(v.reshape(util.shapeshifter(v.shape, target, mode = 'right')), target)
target = (N,n) + _dim.get(k,())
obj = np.broadcast_to(np.array(v).reshape(util.shapeshifter(() if v is None else v.shape,
target,
mode = 'right')),
target)
for i in range(N):
if k in ['phase','O','v','V_e']:
if k in _constituent_properties:
for j in range(n):
mat[i]['constituents'][j][k] = obj[i,j].item() if isinstance(obj[i,j],np.generic) else obj[i,j]
else:
@ -508,4 +546,8 @@ class ConfigMaterial(Config):
dup = self.copy()
dup['material'] = dup['material'] + mat if 'material' in dup else mat
for what in [item for item in ['phase','homogenization'] if shaped[item] is not None]:
for k in np.unique(shaped[what]): # type: ignore
if k not in dup[what]: dup[what][str(k)] = None
return dup

18
python/damask/_crystal.py
View File

@ -1,4 +1,4 @@
from typing import Union, Dict, List, Tuple
from typing import Optional, Union, Dict, List, Tuple
import numpy as np
@ -324,10 +324,10 @@ class Crystal():
"""
def __init__(self, *,
family: CrystalFamily = None,
lattice: CrystalLattice = None,
a: float = None, b: float = None, c: float = None,
alpha: float = None, beta: float = None, gamma: float = None,
family: Optional[CrystalFamily] = None,
lattice: Optional[CrystalLattice] = None,
a: Optional[float] = None, b: Optional[float] = None, c: Optional[float] = None,
alpha: Optional[float] = None, beta: Optional[float] = None, gamma: Optional[float] = None,
degrees: bool = False):
"""
New representation of a crystal.
@ -690,8 +690,8 @@ class Crystal():
self.lattice[-1],None),dtype=float)
def to_lattice(self, *,
direction: FloatSequence = None,
plane: FloatSequence = None) -> np.ndarray:
direction: Optional[FloatSequence] = None,
plane: Optional[FloatSequence] = None) -> np.ndarray:
"""
Calculate lattice vector corresponding to crystal frame direction or plane normal.
@ -717,8 +717,8 @@ class Crystal():
def to_frame(self, *,
uvw: FloatSequence = None,
hkl: FloatSequence = None) -> np.ndarray:
uvw: Optional[FloatSequence] = None,
hkl: Optional[FloatSequence] = None) -> np.ndarray:
"""
Calculate crystal frame vector corresponding to lattice direction [uvw] or plane normal (hkl).

57
python/damask/_grid.py
View File

@ -4,7 +4,7 @@ import warnings
import multiprocessing as mp
from functools import partial
import typing
from typing import Union, Optional, TextIO, List, Sequence, Dict
from typing import Optional, Union, TextIO, List, Sequence, Dict
from pathlib import Path
import numpy as np
@ -34,8 +34,8 @@ class Grid:
material: np.ndarray,
size: FloatSequence,
origin: FloatSequence = np.zeros(3),
initial_conditions: Dict[str,np.ndarray] = None,
comments: Union[str, Sequence[str]] = None):
initial_conditions: Optional[Dict[str,np.ndarray]] = None,
comments: Union[None, str, Sequence[str]] = None):
"""
New geometry definition for grid solvers.
@ -50,7 +50,7 @@ class Grid:
Coordinates of grid origin in meter. Defaults to [0.0,0.0,0.0].
initial_conditions : dictionary, optional
Labels and values of the inital conditions at each material point.
comments : str or sequence of str, optional
comments : (sequence of) str, optional
Additional, human-readable information, e.g. history of operations.
"""
@ -348,9 +348,11 @@ class Grid:
@staticmethod
def load_DREAM3D(fname: Union[str, Path],
feature_IDs: str = None, cell_data: str = None,
phases: str = 'Phases', Euler_angles: str = 'EulerAngles',
base_group: str = None) -> 'Grid':
feature_IDs: Optional[str] = None,
cell_data: Optional[str] = None,
phases: str = 'Phases',
Euler_angles: str = 'EulerAngles',
base_group: Optional[str] = None) -> 'Grid':
"""
Load DREAM.3D (HDF5) file.
@ -427,7 +429,7 @@ class Grid:
coordinates : str
Label of the vector column containing the spatial coordinates.
Need to be ordered (1./x fast, 3./z slow).
labels : str or sequence of str
labels : (sequence of) str
Label(s) of the columns containing the material definition.
Each unique combination of values results in one material ID.
@ -463,7 +465,7 @@ class Grid:
size: FloatSequence,
seeds: np.ndarray,
weights: FloatSequence,
material: IntSequence = None,
material: Optional[IntSequence] = None,
periodic: bool = True):
"""
Create grid from Laguerre tessellation.
@ -520,7 +522,7 @@ class Grid:
def from_Voronoi_tessellation(cells: IntSequence,
size: FloatSequence,
seeds: np.ndarray,
material: IntSequence = None,
material: Optional[IntSequence] = None,
periodic: bool = True) -> 'Grid':
"""
Create grid from Voronoi tessellation.
@ -763,9 +765,9 @@ class Grid:
def canvas(self,
cells: IntSequence = None,
offset: IntSequence = None,
fill: int = None) -> 'Grid':
cells: Optional[IntSequence] = None,
offset: Optional[IntSequence] = None,
fill: Optional[int] = None) -> 'Grid':
"""
Crop or enlarge/pad grid.
@ -901,7 +903,7 @@ class Grid:
def rotate(self,
R: Rotation,
fill: int = None) -> 'Grid':
fill: Optional[int] = None) -> 'Grid':
"""
Rotate grid (and pad if required).
@ -972,15 +974,16 @@ class Grid:
# materials: 1
"""
options = ('nearest',False,None)
orig = tuple(map(np.linspace,self.origin + self.size/self.cells*.5,
self.origin + self.size - self.size/self.cells*.5,self.cells))
interpolator = partial(interpolate.RegularGridInterpolator,
points=orig,method='nearest',bounds_error=False,fill_value=None)
new = grid_filters.coordinates0_point(cells,self.size,self.origin)
return Grid(material = interpolate.RegularGridInterpolator(orig,self.material,*options)(new).astype(int),
return Grid(material = interpolator(values=self.material)(new).astype(int),
size = self.size,
origin = self.origin,
initial_conditions = {k: interpolate.RegularGridInterpolator(orig,v,*options)(new)
initial_conditions = {k: interpolator(values=v)(new)
for k,v in self.initial_conditions.items()},
comments = self.comments+[util.execution_stamp('Grid','scale')],
)
@ -1043,9 +1046,9 @@ class Grid:
Parameters
----------
from_material : int or sequence of int
from_material : (sequence of) int
Material indices to be substituted.
to_material : int or sequence of int
to_material : (sequence of) int
New material indices.
Returns
@ -1092,10 +1095,10 @@ class Grid:
def clean(self,
distance: float = np.sqrt(3),
selection: IntCollection = None,
selection: Optional[IntCollection] = None,
invert_selection: bool = False,
periodic: bool = True,
rng_seed: NumpyRngSeed = None) -> 'Grid':
rng_seed: Optional[NumpyRngSeed] = None) -> 'Grid':
"""
Smooth grid by selecting most frequent material ID within given stencil at each location.
@ -1104,7 +1107,7 @@ class Grid:
distance : float, optional
Voxel distance checked for presence of other materials.
Defaults to sqrt(3).
selection : int or collection of int, optional
selection : (collection of) int, optional
Material IDs to consider. Defaults to all.
invert_selection : bool, optional
Consider all material IDs except those in selection. Defaults to False.
@ -1162,7 +1165,7 @@ class Grid:
dimension: Union[FloatSequence, IntSequence],
center: Union[FloatSequence, IntSequence],
exponent: Union[FloatSequence, float],
fill: int = None,
fill: Optional[int] = None,
R: Rotation = Rotation(),
inverse: bool = False,
periodic: bool = True) -> 'Grid':
@ -1179,7 +1182,7 @@ class Grid:
Center of the primitive.
If given as integers, cell centers are addressed.
If given as floats, physical coordinates are addressed.
exponent : float or sequence of float, len (3)
exponent : (sequence of) float, len (3)
Exponents for the three axes.
0 gives octahedron (ǀxǀ^(2^0) + ǀyǀ^(2^0) + ǀzǀ^(2^0) < 1)
1 gives sphere (ǀxǀ^(2^1) + ǀyǀ^(2^1) + ǀzǀ^(2^1) < 1)
@ -1253,8 +1256,8 @@ class Grid:
def vicinity_offset(self,
distance: float = np.sqrt(3),
offset: int = None,
selection: IntCollection = None,
offset: Optional[int] = None,
selection: Optional[IntCollection] = None,
invert_selection: bool = False,
periodic: bool = True) -> 'Grid':
"""
@ -1271,7 +1274,7 @@ class Grid:
offset : int, optional
Offset (positive or negative) to tag material IDs.
Defaults to material.max()+1.
selection : int or collection of int, optional
selection : (collection of) int, optional
Material IDs that trigger an offset.
Defaults to any other than own material ID.
invert_selection : bool, optional

68
python/damask/_orientation.py
View File

@ -1,6 +1,6 @@
import inspect
import copy
from typing import Union, Callable, Dict, Any, Tuple, TypeVar
from typing import Optional, Union, Callable, Dict, Any, Tuple, TypeVar
import numpy as np
@ -94,14 +94,14 @@ class Orientation(Rotation,Crystal):
"""
@util.extend_docstring(_parameter_doc)
@util.extend_docstring(extra_parameters=_parameter_doc)
def __init__(self,
rotation: Union[FloatSequence, Rotation] = np.array([1.,0.,0.,0.]),
*,
family: CrystalFamily = None,
lattice: CrystalLattice = None,
a: float = None, b: float = None, c: float = None,
alpha: float = None, beta: float = None, gamma: float = None,
family: Optional[CrystalFamily] = None,
lattice: Optional[CrystalLattice] = None,
a: Optional[float] = None, b: Optional[float] = None, c: Optional[float] = None,
alpha: Optional[float] = None, beta: Optional[float] = None, gamma: Optional[float] = None,
degrees: bool = False):
"""
New orientation.
@ -131,7 +131,7 @@ class Orientation(Rotation,Crystal):
def __copy__(self: MyType,
rotation: Union[FloatSequence, Rotation] = None) -> MyType:
rotation: Union[None, FloatSequence, Rotation] = None) -> MyType:
"""
Return deepcopy(self).
@ -300,84 +300,95 @@ class Orientation(Rotation,Crystal):
@classmethod
@util.extended_docstring(Rotation.from_random, _parameter_doc)
@util.extend_docstring(Rotation.from_random,
extra_parameters=_parameter_doc)
def from_random(cls, **kwargs) -> 'Orientation':
kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_random)
return cls(rotation=Rotation.from_random(**kwargs_rot),**kwargs_ori)
@classmethod
@util.extended_docstring(Rotation.from_quaternion,_parameter_doc)
@util.extend_docstring(Rotation.from_quaternion,
extra_parameters=_parameter_doc)
def from_quaternion(cls, **kwargs) -> 'Orientation':
kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_quaternion)
return cls(rotation=Rotation.from_quaternion(**kwargs_rot),**kwargs_ori)
@classmethod
@util.extended_docstring(Rotation.from_Euler_angles,_parameter_doc)
@util.extend_docstring(Rotation.from_Euler_angles,
extra_parameters=_parameter_doc)
def from_Euler_angles(cls, **kwargs) -> 'Orientation':
kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_Euler_angles)
return cls(rotation=Rotation.from_Euler_angles(**kwargs_rot),**kwargs_ori)
@classmethod
@util.extended_docstring(Rotation.from_axis_angle,_parameter_doc)
@util.extend_docstring(Rotation.from_axis_angle,
extra_parameters=_parameter_doc)
def from_axis_angle(cls, **kwargs) -> 'Orientation':
kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_axis_angle)
return cls(rotation=Rotation.from_axis_angle(**kwargs_rot),**kwargs_ori)
@classmethod
@util.extended_docstring(Rotation.from_basis,_parameter_doc)
@util.extend_docstring(Rotation.from_basis,
extra_parameters=_parameter_doc)
def from_basis(cls, **kwargs) -> 'Orientation':
kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_basis)
return cls(rotation=Rotation.from_basis(**kwargs_rot),**kwargs_ori)
@classmethod
@util.extended_docstring(Rotation.from_matrix,_parameter_doc)
@util.extend_docstring(Rotation.from_matrix,
extra_parameters=_parameter_doc)
def from_matrix(cls, **kwargs) -> 'Orientation':
kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_matrix)
return cls(rotation=Rotation.from_matrix(**kwargs_rot),**kwargs_ori)
@classmethod
@util.extended_docstring(Rotation.from_Rodrigues_vector,_parameter_doc)
@util.extend_docstring(Rotation.from_Rodrigues_vector,
extra_parameters=_parameter_doc)
def from_Rodrigues_vector(cls, **kwargs) -> 'Orientation':
kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_Rodrigues_vector)
return cls(rotation=Rotation.from_Rodrigues_vector(**kwargs_rot),**kwargs_ori)
@classmethod
@util.extended_docstring(Rotation.from_homochoric,_parameter_doc)
@util.extend_docstring(Rotation.from_homochoric,
extra_parameters=_parameter_doc)
def from_homochoric(cls, **kwargs) -> 'Orientation':
kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_homochoric)
return cls(rotation=Rotation.from_homochoric(**kwargs_rot),**kwargs_ori)
@classmethod
@util.extended_docstring(Rotation.from_cubochoric,_parameter_doc)
@util.extend_docstring(Rotation.from_cubochoric,
extra_parameters=_parameter_doc)
def from_cubochoric(cls, **kwargs) -> 'Orientation':
kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_cubochoric)
return cls(rotation=Rotation.from_cubochoric(**kwargs_rot),**kwargs_ori)
@classmethod
@util.extended_docstring(Rotation.from_spherical_component,_parameter_doc)
@util.extend_docstring(Rotation.from_spherical_component,
extra_parameters=_parameter_doc)
def from_spherical_component(cls, **kwargs) -> 'Orientation':
kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_spherical_component)
return cls(rotation=Rotation.from_spherical_component(**kwargs_rot),**kwargs_ori)
@classmethod
@util.extended_docstring(Rotation.from_fiber_component,_parameter_doc)
@util.extend_docstring(Rotation.from_fiber_component,
extra_parameters=_parameter_doc)
def from_fiber_component(cls, **kwargs) -> 'Orientation':
kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_fiber_component)
return cls(rotation=Rotation.from_fiber_component(**kwargs_rot),**kwargs_ori)
@classmethod
@util.extend_docstring(_parameter_doc)
@util.extend_docstring(extra_parameters=_parameter_doc)
def from_directions(cls,
uvw: FloatSequence,
hkl: FloatSequence,
@ -392,6 +403,10 @@ class Orientation(Rotation,Crystal):
hkl : numpy.ndarray, shape (...,3)
Lattice plane normal aligned with lab frame z-direction.
Returns
-------
new : damask.Orientation
"""
o = cls(**kwargs)
x = o.to_frame(uvw=uvw)
@ -538,8 +553,7 @@ class Orientation(Rotation,Crystal):
Notes
-----
Currently requires same crystal family for both orientations.
For extension to cases with differing symmetry see A. Heinz and P. Neumann 1991 and 10.1107/S0021889808016373.
Requires same crystal family for both orientations.
Examples
--------
@ -569,6 +583,8 @@ class Orientation(Rotation,Crystal):
>>> plt.show()
"""
# For extension to cases with differing symmetry see
# https://doi.org/10.1107/S0021889808016373 and https://doi.org/10.1107/S0108767391006864
if self.family != other.family:
raise NotImplementedError('disorientation between different crystal families')
@ -601,7 +617,7 @@ class Orientation(Rotation,Crystal):
def average(self,
weights: FloatSequence = None,
weights: Optional[FloatSequence] = None,
return_cloud: bool = False):
"""
Return orientation average over last dimension.
@ -803,8 +819,8 @@ class Orientation(Rotation,Crystal):
# functions that require lattice, not just family
def to_pole(self, *,
uvw: FloatSequence = None,
hkl: FloatSequence = None,
uvw: Optional[FloatSequence] = None,
hkl: Optional[FloatSequence] = None,
with_symmetry: bool = False,
normalize: bool = True) -> np.ndarray:
"""
@ -845,8 +861,8 @@ class Orientation(Rotation,Crystal):
def Schmid(self, *,
N_slip: IntSequence = None,
N_twin: IntSequence = None) -> np.ndarray:
N_slip: Optional[IntSequence] = None,
N_twin: Optional[IntSequence] = None) -> np.ndarray:
u"""
Calculate Schmid matrix P = d ⨂ n in the lab frame for selected deformation systems.

72
python/damask/_result.py
View File

@ -11,7 +11,7 @@ from pathlib import Path
from functools import partial
from collections import defaultdict
from collections.abc import Iterable
from typing import Union, Callable, Any, Sequence, Literal, Dict, List, Tuple, Optional
from typing import Optional, Union, Callable, Any, Sequence, Literal, Dict, List, Tuple
import h5py
import numpy as np
@ -189,11 +189,11 @@ class Result:
def _manage_view(self,
action: Literal['set', 'add', 'del'],
increments: Union[int, Sequence[int], str, Sequence[str], bool] = None,
times: Union[float, Sequence[float], str, Sequence[str], bool] = None,
phases: Union[str, Sequence[str], bool] = None,
homogenizations: Union[str, Sequence[str], bool] = None,
fields: Union[str, Sequence[str], bool] = None) -> "Result":
increments: Union[None, int, Sequence[int], str, Sequence[str], bool] = None,
times: Union[None, float, Sequence[float], str, Sequence[str], bool] = None,
phases: Union[None, str, Sequence[str], bool] = None,
homogenizations: Union[None, str, Sequence[str], bool] = None,
fields: Union[None, str, Sequence[str], bool] = None) -> "Result":
"""
Manages the visibility of the groups.
@ -256,8 +256,8 @@ class Result:
def increments_in_range(self,
start: Union[str, int] = None,
end: Union[str, int] = None) -> Sequence[int]:
start: Union[None, str, int] = None,
end: Union[None, str, int] = None) -> Sequence[int]:
"""
Get all increments within a given range.
@ -280,8 +280,8 @@ class Result:
return [i for i in self.incs if s <= i <= e]
def times_in_range(self,
start: float = None,
end: float = None) -> Sequence[float]:
start: Optional[float] = None,
end: Optional[float] = None) -> Sequence[float]:
"""
Get times of all increments within a given time range.
@ -304,12 +304,12 @@ class Result:
def view(self,*,
increments: Union[int, Sequence[int], str, Sequence[str], bool] = None,
times: Union[float, Sequence[float], str, Sequence[str], bool] = None,
phases: Union[str, Sequence[str], bool] = None,
homogenizations: Union[str, Sequence[str], bool] = None,
fields: Union[str, Sequence[str], bool] = None,
protected: bool = None) -> "Result":
increments: Union[None, int, Sequence[int], str, Sequence[str], bool] = None,
times: Union[None, float, Sequence[float], str, Sequence[str], bool] = None,
phases: Union[None, str, Sequence[str], bool] = None,
homogenizations: Union[None, str, Sequence[str], bool] = None,
fields: Union[None, str, Sequence[str], bool] = None,
protected: Optional[bool] = None) -> "Result":
"""
Set view.
@ -361,11 +361,11 @@ class Result:
def view_more(self,*,
increments: Union[int, Sequence[int], str, Sequence[str], bool] = None,
times: Union[float, Sequence[float], str, Sequence[str], bool] = None,
phases: Union[str, Sequence[str], bool] = None,
homogenizations: Union[str, Sequence[str], bool] = None,
fields: Union[str, Sequence[str], bool] = None) -> "Result":
increments: Union[None, int, Sequence[int], str, Sequence[str], bool] = None,
times: Union[None, float, Sequence[float], str, Sequence[str], bool] = None,
phases: Union[None, str, Sequence[str], bool] = None,
homogenizations: Union[None, str, Sequence[str], bool] = None,
fields: Union[None, str, Sequence[str], bool] = None) -> "Result":
"""
Add to view.
@ -404,11 +404,11 @@ class Result:
def view_less(self,*,
increments: Union[int, Sequence[int], str, Sequence[str], bool] = None,
times: Union[float, Sequence[float], str, Sequence[str], bool] = None,
phases: Union[str, Sequence[str], bool] = None,
homogenizations: Union[str, Sequence[str], bool] = None,
fields: Union[str, Sequence[str], bool] = None) -> "Result":
increments: Union[None, int, Sequence[int], str, Sequence[str], bool] = None,
times: Union[None, float, Sequence[float], str, Sequence[str], bool] = None,
phases: Union[None, str, Sequence[str], bool] = None,
homogenizations: Union[None, str, Sequence[str], bool] = None,
fields: Union[None, str, Sequence[str], bool] = None) -> "Result":
"""
Remove from view.
@ -650,7 +650,7 @@ class Result:
formula: str,
name: str,
unit: str = 'n/a',
description: str = None):
description: Optional[str] = None):
"""
Add result of a general formula.
@ -966,7 +966,7 @@ class Result:
}
def add_equivalent_Mises(self,
T_sym: str,
kind: str = None):
kind: Optional[str] = None):
"""
Add the equivalent Mises stress or strain of a symmetric tensor.
@ -1021,7 +1021,7 @@ class Result:
}
def add_norm(self,
x: str,
ord: Union[int, float, Literal['fro', 'nuc']] = None):
ord: Union[None, int, float, Literal['fro', 'nuc']] = None):
"""
Add the norm of vector or tensor.
@ -1101,8 +1101,8 @@ class Result:
def add_pole(self,
q: str = 'O',
*,
uvw: FloatSequence = None,
hkl: FloatSequence = None,
uvw: Optional[FloatSequence] = None,
hkl: Optional[FloatSequence] = None,
with_symmetry: bool = False,
normalize: bool = True):
"""
@ -1593,7 +1593,7 @@ class Result:
output: Union[str, List[str]] = '*',
flatten: bool = True,
prune: bool = True,
constituents: IntSequence = None,
constituents: Optional[IntSequence] = None,
fill_float: float = np.nan,
fill_int: int = 0) -> Optional[Dict[str,Any]]:
"""
@ -1683,7 +1683,7 @@ class Result:
def export_XDMF(self,
output: Union[str, List[str]] = '*',
target_dir: Union[str, Path] = None,
target_dir: Union[None, str, Path] = None,
absolute_path: bool = False):
"""
Write XDMF file to directly visualize data from DADF5 file.
@ -1811,8 +1811,8 @@ class Result:
def export_VTK(self,
output: Union[str,List[str]] = '*',
mode: str = 'cell',
constituents: IntSequence = None,
target_dir: Union[str, Path] = None,
constituents: Optional[IntSequence] = None,
target_dir: Union[None, str, Path] = None,
fill_float: float = np.nan,
fill_int: int = 0,
parallel: bool = True):
@ -1958,7 +1958,7 @@ class Result:
def export_simulation_setup(self,
output: Union[str, List[str]] = '*',
target_dir: Union[str, Path] = None,
target_dir: Union[None, str, Path] = None,
overwrite: bool = False,
):
"""

441
python/damask/_rotation.py
View File

@ -1,5 +1,5 @@
import copy
from typing import Union, Sequence, Tuple, Literal, List, TypeVar
from typing import Optional, Union, Sequence, Tuple, Literal, List, TypeVar
import numpy as np
@ -26,11 +26,8 @@ class Rotation:
- Coordinate frames are right-handed.
- A rotation angle ω is taken to be positive for a counterclockwise rotation
when viewing from the end point of the rotation axis towards the origin.
- Rotations will be interpreted in the passive sense.
- Euler angle triplets are implemented using the Bunge convention,
with angular ranges of [0,2π], [0,π], [0,2π].
- The rotation angle ω is limited to the interval [0,π].
- The real part of a quaternion is positive, Re(q) ≥ 0
- Rotations will be interpreted in the passive sense, i.e. as rotation of
the coordinate frame.
- P = -1 (as default).
Examples
@ -66,7 +63,7 @@ class Rotation:
__slots__ = ['quaternion']
def __init__(self,
rotation: Union[FloatSequence, 'Rotation'] = np.array([1.0,0.0,0.0,0.0])):
rotation: Union[FloatSequence, 'Rotation'] = np.array([1.,0.,0.,0.])):
"""
New rotation.
@ -82,7 +79,7 @@ class Rotation:
if isinstance(rotation,Rotation):
self.quaternion = rotation.quaternion.copy()
elif np.array(rotation).shape[-1] == 4:
self.quaternion = np.array(rotation)
self.quaternion = np.array(rotation,dtype=float)
else:
raise TypeError('"rotation" is neither a Rotation nor a quaternion')
@ -99,7 +96,7 @@ class Rotation:
def __copy__(self: MyType,
rotation: Union[FloatSequence, 'Rotation'] = None) -> MyType:
rotation: Union[None, FloatSequence, 'Rotation'] = None) -> MyType:
"""
Return deepcopy(self).
@ -141,7 +138,7 @@ class Rotation:
"""
return NotImplemented if not isinstance(other, Rotation) else \
np.logical_or(np.all(self.quaternion == other.quaternion,axis=-1),
np.all(self.quaternion == -1.0*other.quaternion,axis=-1))
np.all(self.quaternion == -1.*other.quaternion,axis=-1))
def __ne__(self,
@ -161,8 +158,8 @@ class Rotation:
def isclose(self: MyType,
other: MyType,
rtol: float = 1e-5,
atol: float = 1e-8,
rtol: float = 1.e-5,
atol: float = 1.e-8,
equal_nan: bool = True) -> bool:
"""
Report where values are approximately equal to corresponding ones of other Rotation.
@ -187,13 +184,13 @@ class Rotation:
s = self.quaternion
o = other.quaternion
return np.logical_or(np.all(np.isclose(s, o,rtol,atol,equal_nan),axis=-1),
np.all(np.isclose(s,-1.0*o,rtol,atol,equal_nan),axis=-1))
np.all(np.isclose(s,-1.*o,rtol,atol,equal_nan),axis=-1))
def allclose(self: MyType,
other: MyType,
rtol: float = 1e-5,
atol: float = 1e-8,
rtol: float = 1.e-5,
atol: float = 1.e-8,
equal_nan: bool = True) -> Union[np.bool_, bool]:
"""
Test whether all values are approximately equal to corresponding ones of other Rotation.
@ -250,7 +247,7 @@ class Rotation:
"""
dup = self.copy()
dup.quaternion[...,1:] *= -1
dup.quaternion[...,1:] *= -1.
return dup
@ -393,9 +390,9 @@ class Rotation:
if self.shape + (3,) == other.shape:
q_m = self.quaternion[...,0]
p_m = self.quaternion[...,1:]
A = q_m**2.0 - np.einsum('...i,...i',p_m,p_m)
B = 2.0 * np.einsum('...i,...i',p_m,other)
C = 2.0 * _P * q_m
A = q_m**2 - np.einsum('...i,...i',p_m,p_m)
B = 2. * np.einsum('...i,...i',p_m,other)
C = 2. * _P * q_m
return np.block([(A * other[...,i]).reshape(self.shape+(1,)) +
(B * p_m[...,i]).reshape(self.shape+(1,)) +
(C * ( p_m[...,(i+1)%3]*other[...,(i+2)%3]\
@ -419,7 +416,7 @@ class Rotation:
def _standardize(self: MyType) -> MyType:
"""Standardize quaternion (ensure positive real hemisphere)."""
self.quaternion[self.quaternion[...,0] < 0.0] *= -1
self.quaternion[self.quaternion[...,0] < 0.] *= -1.
return self
@ -468,7 +465,7 @@ class Rotation:
Parameters
----------
shape : int or sequence of ints
shape : (sequence of) int
New shape, number of elements needs to match the original shape.
If an integer is supplied, then the result will be a 1-D array of that length.
order : {'C', 'F', 'A'}, optional
@ -496,7 +493,7 @@ class Rotation:
Parameters
----------
shape : int or sequence of ints
shape : (sequence of) int
Shape of broadcasted array, needs to be compatible with the original shape.
mode : str, optional
Where to preferentially locate missing dimensions.
@ -514,7 +511,7 @@ class Rotation:
def average(self: MyType,
weights: FloatSequence = None) -> MyType:
weights: Optional[FloatSequence] = None) -> MyType:
"""
Average along last array dimension.
@ -540,7 +537,7 @@ class Rotation:
weights_ = np.ones(self.shape,dtype=float) if weights is None else np.array(weights,float)
eig, vec = np.linalg.eig(np.sum(_M(self.quaternion) * weights_[...,np.newaxis,np.newaxis],axis=-3) \
eig, vec = np.linalg.eig(np.sum(_M(self.quaternion) * weights_[...,np.newaxis,np.newaxis],axis=-3)
/np.sum( weights_[...,np.newaxis,np.newaxis],axis=-3))
return self.copy(Rotation.from_quaternion(np.real(
@ -751,6 +748,7 @@ class Rotation:
@staticmethod
def from_quaternion(q: Union[Sequence[FloatSequence], np.ndarray],
accept_homomorph: bool = False,
normalize: bool = False,
P: Literal[1, -1] = -1) -> 'Rotation':
"""
Initialize from quaternion.
@ -762,23 +760,29 @@ class Rotation:
accept_homomorph : bool, optional
Allow homomorphic variants, i.e. q_0 < 0 (negative real hemisphere).
Defaults to False.
normalize: bool, optional
Allow ǀqǀ ≠ 1. Defaults to False.
P : int ∈ {-1,1}, optional
Sign convention. Defaults to -1.
Returns
-------
new : damask.Rotation
"""
qu = np.array(q,dtype=float)
if qu.shape[:-2:-1] != (4,):
raise ValueError('invalid shape')
if abs(P) != 1:
raise ValueError('P ∉ {-1,1}')
if qu.shape[:-2:-1] != (4,): raise ValueError('invalid shape')
if abs(P) != 1: raise ValueError('P ∉ {-1,1}')
qu[...,1:4] *= -P
if accept_homomorph:
qu[qu[...,0] < 0.0] *= -1
else:
if np.any(qu[...,0] < 0.0):
qu[qu[...,0]<0.] *= -1.
elif np.any(qu[...,0] < 0.):
raise ValueError('quaternion with negative first (real) component')
if not np.all(np.isclose(np.linalg.norm(qu,axis=-1), 1.0,rtol=0.0)):
if normalize:
qu /= np.linalg.norm(qu,axis=-1,keepdims=True)
elif not np.allclose(np.linalg.norm(qu,axis=-1),1.,rtol=1.e-8):
raise ValueError('quaternion is not of unit length')
return Rotation(qu)
@ -797,17 +801,20 @@ class Rotation:
degrees : bool, optional
Euler angles are given in degrees. Defaults to False.
Returns
-------
new : damask.Rotation
Notes
-----
Bunge Euler angles correspond to a rotation axis sequence of z–x'–z''.
"""
eu = np.array(phi,dtype=float)
if eu.shape[:-2:-1] != (3,):
raise ValueError('invalid shape')
if eu.shape[:-2:-1] != (3,): raise ValueError('invalid shape')
eu = np.radians(eu) if degrees else eu
if np.any(eu < 0.0) or np.any(eu > 2.0*np.pi) or np.any(eu[...,1] > np.pi): # ToDo: No separate check for PHI
if np.any(eu < 0.) or np.any(eu > np.pi*np.array([2.,1.,2.])):
raise ValueError('Euler angles outside of [0..2π],[0..π],[0..2π]')
return Rotation(Rotation._eu2qu(eu))
@ -832,20 +839,23 @@ class Rotation:
P : int ∈ {-1,1}, optional
Sign convention. Defaults to -1.
Returns
-------
new : damask.Rotation
"""
ax = np.array(axis_angle,dtype=float)
if ax.shape[:-2:-1] != (4,):
raise ValueError('invalid shape')
if abs(P) != 1:
raise ValueError('P ∉ {-1,1}')
if ax.shape[:-2:-1] != (4,): raise ValueError('invalid shape')
if abs(P) != 1: raise ValueError('P ∉ {-1,1}')
ax[...,0:3] *= -P
if degrees: ax[..., 3] = np.radians(ax[...,3])
if normalize: ax[...,0:3] /= np.linalg.norm(ax[...,0:3],axis=-1,keepdims=True)
if np.any(ax[...,3] < 0.0) or np.any(ax[...,3] > np.pi):
if np.any(ax[...,3] < 0.) or np.any(ax[...,3] > np.pi):
raise ValueError('axis–angle rotation angle outside of [0..π]')
if not np.all(np.isclose(np.linalg.norm(ax[...,0:3],axis=-1), 1.0)):
print(np.linalg.norm(ax[...,0:3],axis=-1))
if normalize:
ax[...,0:3] /= np.linalg.norm(ax[...,0:3],axis=-1,keepdims=True)
elif not np.allclose(np.linalg.norm(ax[...,0:3],axis=-1),1.):
raise ValueError('axis–angle rotation axis is not of unit length')
return Rotation(Rotation._ax2qu(ax))
@ -866,24 +876,29 @@ class Rotation:
reciprocal : bool, optional
Basis vectors are given in reciprocal (instead of real) space. Defaults to False.
Returns
-------
new : damask.Rotation
"""
om = np.array(basis,dtype=float)
if om.shape[-2:] != (3,3):
raise ValueError('invalid shape')
if om.shape[-2:] != (3,3): raise ValueError('invalid shape')
if reciprocal:
om = np.linalg.inv(tensor.transpose(om)/np.pi) # transform reciprocal basis set
orthonormal = False # contains stretch
if not orthonormal:
(U,S,Vh) = np.linalg.svd(om) # singular value decomposition
om = np.einsum('...ij,...jl',U,Vh)
if not np.all(np.isclose(np.linalg.det(om),1.0)):
raise ValueError('orientation matrix has determinant ≠ 1')
if not np.all(np.isclose(np.einsum('...i,...i',om[...,0],om[...,1]), 0.0)) \
or not np.all(np.isclose(np.einsum('...i,...i',om[...,1],om[...,2]), 0.0)) \
or not np.all(np.isclose(np.einsum('...i,...i',om[...,2],om[...,0]), 0.0)):
elif not np.allclose(np.einsum('...i,...i',om[...,0],om[...,1]),0.) \
or not np.allclose(np.einsum('...i,...i',om[...,1],om[...,2]),0.) \
or not np.allclose(np.einsum('...i,...i',om[...,2],om[...,0]),0.):
raise ValueError('orientation matrix is not orthogonal')
if not np.allclose(np.linalg.det(om),1.):
raise ValueError('orientation matrix has determinant ≠ 1')
return Rotation(Rotation._om2qu(om))
@staticmethod
@ -896,6 +911,10 @@ class Rotation:
R : numpy.ndarray, shape (...,3,3)
Rotation matrix with det(R) = 1 and R.T ∙ R = I.
Returns
-------
new : damask.Rotation
"""
return Rotation.from_basis(R)
@ -912,9 +931,13 @@ class Rotation:
b : numpy.ndarray, shape (...,2,3)
Corresponding three-dimensional vectors of second basis.
Returns
-------
new : damask.Rotation
"""
a_ = np.array(a)
b_ = np.array(b)
a_ = np.array(a,dtype=float)
b_ = np.array(b,dtype=float)
if a_.shape[-2:] != (2,3) or b_.shape[-2:] != (2,3) or a_.shape != b_.shape:
raise ValueError('invalid shape')
am = np.stack([ a_[...,0,:],
@ -944,18 +967,21 @@ class Rotation:
P : int ∈ {-1,1}, optional
Sign convention. Defaults to -1.
Returns
-------
new : damask.Rotation
"""
ro = np.array(rho,dtype=float)
if ro.shape[:-2:-1] != (4,):
raise ValueError('invalid shape')
if abs(P) != 1:
raise ValueError('P ∉ {-1,1}')
if ro.shape[:-2:-1] != (4,): raise ValueError('invalid shape')
if abs(P) != 1: raise ValueError('P ∉ {-1,1}')
ro[...,0:3] *= -P
if normalize: ro[...,0:3] /= np.linalg.norm(ro[...,0:3],axis=-1,keepdims=True)
if np.any(ro[...,3] < 0.0):
raise ValueError('Rodrigues vector rotation angle is negative')
if not np.all(np.isclose(np.linalg.norm(ro[...,0:3],axis=-1), 1.0)):
if np.any(ro[...,3] < 0.): raise ValueError('Rodrigues vector rotation angle is negative')
if normalize:
ro[...,0:3] /= np.linalg.norm(ro[...,0:3],axis=-1,keepdims=True)
elif not np.allclose(np.linalg.norm(ro[...,0:3],axis=-1),1.):
raise ValueError('Rodrigues vector rotation axis is not of unit length')
return Rotation(Rotation._ro2qu(ro))
@ -973,16 +999,18 @@ class Rotation:
P : int ∈ {-1,1}, optional
Sign convention. Defaults to -1.
Returns
-------
new : damask.Rotation
"""
ho = np.array(h,dtype=float)
if ho.shape[:-2:-1] != (3,):
raise ValueError('invalid shape')
if abs(P) != 1:
raise ValueError('P ∉ {-1,1}')
if ho.shape[:-2:-1] != (3,): raise ValueError('invalid shape')
if abs(P) != 1: raise ValueError('P ∉ {-1,1}')
ho *= -P
if np.any(np.linalg.norm(ho,axis=-1) >_R1+1e-9):
if np.any(np.linalg.norm(ho,axis=-1) > _R1+1.e-9):
raise ValueError('homochoric coordinate outside of the sphere')
return Rotation(Rotation._ho2qu(ho))
@ -1000,14 +1028,15 @@ class Rotation:
P : int ∈ {-1,1}, optional
Sign convention. Defaults to -1.
Returns
-------
new : damask.Rotation
"""
cu = np.array(x,dtype=float)
if cu.shape[:-2:-1] != (3,):
raise ValueError('invalid shape')
if abs(P) != 1:
raise ValueError('P ∉ {-1,1}')
if np.abs(np.max(cu)) > np.pi**(2./3.) * 0.5+1e-9:
if cu.shape[:-2:-1] != (3,): raise ValueError('invalid shape')
if abs(P) != 1: raise ValueError('P ∉ {-1,1}')
if np.abs(np.max(cu)) > np.pi**(2./3.) * 0.5+1.e-9:
raise ValueError('cubochoric coordinate outside of the cube')
ho = -P * Rotation._cu2ho(cu)
@ -1016,29 +1045,33 @@ class Rotation:
@staticmethod
def from_random(shape: Union[int, IntSequence] = None,
rng_seed: NumpyRngSeed = None) -> 'Rotation':
def from_random(shape: Union[None, int, IntSequence] = None,
rng_seed: Optional[NumpyRngSeed] = None) -> 'Rotation':
"""
Initialize with samples from a uniform distribution.
Parameters
----------
shape : int or sequence of ints, optional
shape : (sequence of) int, optional
Shape of the returned array. Defaults to None, which gives a scalar.
rng_seed : {None, int, array_like[ints], SeedSequence, BitGenerator, Generator}, optional
A seed to initialize the BitGenerator.
Defaults to None, i.e. unpredictable entropy will be pulled from the OS.
Returns
-------
new : damask.Rotation
"""
rng = np.random.default_rng(rng_seed)
r = rng.random(3 if shape is None else tuple(shape)+(3,) if hasattr(shape, '__iter__') else (shape,3)) # type: ignore
A = np.sqrt(r[...,2])
B = np.sqrt(1.0-r[...,2])
q = np.stack([np.cos(2.0*np.pi*r[...,0])*A,
np.sin(2.0*np.pi*r[...,1])*B,
np.cos(2.0*np.pi*r[...,1])*B,
np.sin(2.0*np.pi*r[...,0])*A],axis=-1)
B = np.sqrt(1.-r[...,2])
q = np.stack([np.cos(2.*np.pi*r[...,0])*A,
np.sin(2.*np.pi*r[...,1])*B,
np.cos(2.*np.pi*r[...,1])*B,
np.sin(2.*np.pi*r[...,0])*A],axis=-1)
return Rotation(q if shape is None else q.reshape(r.shape[:-1]+(4,)))._standardize()
@ -1046,10 +1079,10 @@ class Rotation:
@staticmethod
def from_ODF(weights: np.ndarray,
phi: np.ndarray,
shape: Union[int, IntSequence] = None,
shape: Union[None, int, IntSequence] = None,
degrees: bool = False,
fractions: bool = True,
rng_seed: NumpyRngSeed = None) -> 'Rotation':
rng_seed: Optional[NumpyRngSeed] = None) -> 'Rotation':
"""
Initialize with samples from a binned orientation distribution function (ODF).
@ -1059,7 +1092,7 @@ class Rotation:
Texture intensity values (probability density or volume fraction) at Euler space grid points.
phi : numpy.ndarray, shape (n,3)
Grid coordinates in Euler space at which weights are defined.
shape : int or sequence of ints, optional
shape : (sequence of) int, optional
Shape of the returned array. Defaults to None, which gives a scalar.
degrees : bool, optional
Euler space grid coordinates are in degrees. Defaults to True.
@ -1070,6 +1103,10 @@ class Rotation:
A seed to initialize the BitGenerator.
Defaults to None, i.e. unpredictable entropy will be pulled from the OS.
Returns
-------
new : damask.Rotation
Notes
-----
Due to the distortion of Euler space in the vicinity of ϕ = 0, probability densities, p, defined on
@ -1087,21 +1124,21 @@ class Rotation:
phi_sorted = eu[np.lexsort((eu[:,0],eu[:,1],eu[:,2]))]
steps,size,_ = grid_filters.cellsSizeOrigin_coordinates0_point(phi_sorted)
delta = np.radians(size/steps) if deg else size/steps
return delta[0]*2.0*np.sin(delta[1]/2.0)*delta[2] / 8.0 / np.pi**2 * np.sin(np.radians(eu[:,1]) if deg else eu[:,1])
return delta[0]*2.*np.sin(delta[1]/2.)*delta[2] / 8. / np.pi**2 * np.sin(np.radians(eu[:,1]) if deg else eu[:,1])
dg = 1.0 if fractions else _dg(phi,degrees)
dV_V = dg * np.maximum(0.0,weights.squeeze())
dg = 1. if fractions else _dg(phi,degrees)
dV_V = dg * np.maximum(0.,weights.squeeze())
N = 1 if shape is None else np.prod(shape)
N = 1 if shape is None else np.prod(shape).astype(int)
return Rotation.from_Euler_angles(phi[util.hybrid_IA(dV_V,N,rng_seed)],degrees).reshape(() if shape is None else shape)
@staticmethod
def from_spherical_component(center: 'Rotation',
sigma: float,
shape: Union[int, IntSequence] = None,
shape: Union[None, int, IntSequence] = None,
degrees: bool = False,
rng_seed: NumpyRngSeed = None) -> 'Rotation':
rng_seed: Optional[NumpyRngSeed] = None) -> 'Rotation':
"""
Initialize with samples from a Gaussian distribution around a given center.
@ -1111,7 +1148,7 @@ class Rotation:
Central rotation.
sigma : float
Standard deviation of (Gaussian) misorientation distribution.
shape : int or sequence of ints, optional
shape : (sequence of) int, optional
Shape of the returned array. Defaults to None, which gives a scalar.
degrees : bool, optional
sigma is given in degrees. Defaults to False.
@ -1125,24 +1162,24 @@ class Rotation:
200 orientations:
>>> import damask
>>> center = damask.Rotation.from_Euler_angles([35.0,45.0,0.0],degrees=True)
>>> brass = damask.Rotation.from_spherical_component(center=center,sigma=3.0,shape=200,degrees=True)
>>> center = damask.Rotation.from_Euler_angles([35.,45.,0.],degrees=True)
>>> brass = damask.Rotation.from_spherical_component(center=center,sigma=3.,shape=200,degrees=True)
Create a Goss texture consisting of
100 orientations:
>>> import damask
>>> center = damask.Rotation.from_Euler_angles([0.0,45.0,0.0],degrees=True)
>>> goss = damask.Rotation.from_spherical_component(center=center,sigma=3.0,shape=100,degrees=True)
>>> center = damask.Rotation.from_Euler_angles([0.,45.,0.],degrees=True)
>>> goss = damask.Rotation.from_spherical_component(center=center,sigma=3.,shape=100,degrees=True)
"""
rng = np.random.default_rng(rng_seed)
sigma = np.radians(sigma) if degrees else sigma
N = 1 if shape is None else np.prod(shape)
u,Theta = (rng.random((N,2)) * 2.0 * np.array([1,np.pi]) - np.array([1.0, 0])).T
u,Theta = (rng.random((N,2)) * 2. * np.array([1.,np.pi]) - np.array([1.,0.])).T
omega = abs(rng.normal(scale=sigma,size=N))
p = np.column_stack([np.sqrt(1-u**2)*np.cos(Theta),
np.sqrt(1-u**2)*np.sin(Theta),
p = np.column_stack([np.sqrt(1.-u**2)*np.cos(Theta),
np.sqrt(1.-u**2)*np.sin(Theta),
u, omega])
return Rotation.from_axis_angle(p).reshape(() if shape is None else shape) * center
@ -1151,10 +1188,10 @@ class Rotation:
@staticmethod
def from_fiber_component(crystal: IntSequence,
sample: IntSequence,
sigma: float = 0.0,
shape: Union[int, IntSequence] = None,
sigma: float = 0.,
shape: Union[None, int, IntSequence] = None,
degrees: bool = False,
rng_seed: NumpyRngSeed = None):
rng_seed: Optional[NumpyRngSeed] = None) -> 'Rotation':
"""
Initialize with samples from a Gaussian distribution around a given direction.
@ -1169,7 +1206,7 @@ class Rotation:
sigma : float, optional
Standard deviation of (Gaussian) misorientation distribution.
Defaults to 0.
shape : int or sequence of ints, optional
shape : (sequence of) int, optional
Shape of the returned array. Defaults to None, which gives a scalar.
degrees : bool, optional
sigma and polar coordinates are given in degrees. Defaults to False.
@ -1177,6 +1214,10 @@ class Rotation:
A seed to initialize the BitGenerator.
Defaults to None, i.e. unpredictable entropy will be pulled from the OS.
Returns
-------
new : damask.Rotation
Notes
-----
The crystal direction for (θ=0,φ=0) is [0 0 1],
@ -1201,7 +1242,7 @@ class Rotation:
100 orientations:
>>> import damask
>>> gamma = damask.Rotation.from_fiber_component([54.7,45.0],[0.,0.],shape=100,degrees=True)
>>> gamma = damask.Rotation.from_fiber_component([54.7,45.],[0.,0.],shape=100,degrees=True)
"""
rng = np.random.default_rng(rng_seed)
@ -1211,18 +1252,18 @@ class Rotation:
d_cr = np.array([np.sin(alpha[0])*np.cos(alpha[1]), np.sin(alpha[0])*np.sin(alpha[1]), np.cos(alpha[0])])
d_lab = np.array([np.sin( beta[0])*np.cos( beta[1]), np.sin( beta[0])*np.sin( beta[1]), np.cos( beta[0])])
ax_align = np.append(np.cross(d_lab,d_cr), np.arccos(np.dot(d_lab,d_cr)))
if np.isclose(ax_align[3],0.0): ax_align[:3] = np.array([1,0,0])
R_align = Rotation.from_axis_angle(ax_align if ax_align[3] > 0.0 else -ax_align,normalize=True) # rotate fiber axis from sample to crystal frame
if np.isclose(ax_align[3],0.): ax_align[:3] = np.array([1.,0.,0.])
R_align = Rotation.from_axis_angle(ax_align if ax_align[3] > 0. else -ax_align,normalize=True) # rotate fiber axis from sample to crystal frame
N = 1 if shape is None else np.prod(shape)
u,Theta = (rng.random((N,2)) * 2.0 * np.array([1,np.pi]) - np.array([1.0, 0])).T
N = 1 if shape is None else np.prod(shape).astype(int)
u,Theta = (rng.random((N,2)) * 2. * np.array([1.,np.pi]) - np.array([1.,0.])).T
omega = abs(rng.normal(scale=sigma_,size=N))
p = np.column_stack([np.sqrt(1-u**2)*np.cos(Theta),
np.sqrt(1-u**2)*np.sin(Theta),
p = np.column_stack([np.sqrt(1.-u**2)*np.cos(Theta),
np.sqrt(1.-u**2)*np.sin(Theta),
u, omega])
p[:,:3] = np.einsum('ij,...j',np.eye(3)-np.outer(d_lab,d_lab),p[:,:3]) # remove component along fiber axis
f = np.column_stack((np.broadcast_to(d_lab,(N,3)),rng.random(N)*np.pi))
f[::2,:3] *= -1 # flip half the rotation axes to negative sense
f[::2,:3] *= -1. # flip half the rotation axes to negative sense
return (R_align.broadcast_to(N)
* Rotation.from_axis_angle(p,normalize=True)
@ -1263,15 +1304,15 @@ class Rotation:
@staticmethod
def _qu2om(qu: np.ndarray) -> np.ndarray:
qq = qu[...,0:1]**2-(qu[...,1:2]**2 + qu[...,2:3]**2 + qu[...,3:4]**2)
om = np.block([qq + 2.0*qu[...,1:2]**2,
2.0*(qu[...,2:3]*qu[...,1:2]-_P*qu[...,0:1]*qu[...,3:4]),
2.0*(qu[...,3:4]*qu[...,1:2]+_P*qu[...,0:1]*qu[...,2:3]),
2.0*(qu[...,1:2]*qu[...,2:3]+_P*qu[...,0:1]*qu[...,3:4]),
qq + 2.0*qu[...,2:3]**2,
2.0*(qu[...,3:4]*qu[...,2:3]-_P*qu[...,0:1]*qu[...,1:2]),
2.0*(qu[...,1:2]*qu[...,3:4]-_P*qu[...,0:1]*qu[...,2:3]),
2.0*(qu[...,2:3]*qu[...,3:4]+_P*qu[...,0:1]*qu[...,1:2]),
qq + 2.0*qu[...,3:4]**2,
om = np.block([qq + 2.*qu[...,1:2]**2,
2.*(qu[...,2:3]*qu[...,1:2]-_P*qu[...,0:1]*qu[...,3:4]),
2.*(qu[...,3:4]*qu[...,1:2]+_P*qu[...,0:1]*qu[...,2:3]),
2.*(qu[...,1:2]*qu[...,2:3]+_P*qu[...,0:1]*qu[...,3:4]),
qq + 2.*qu[...,2:3]**2,
2.*(qu[...,3:4]*qu[...,2:3]-_P*qu[...,0:1]*qu[...,1:2]),
2.*(qu[...,1:2]*qu[...,3:4]-_P*qu[...,0:1]*qu[...,2:3]),
2.*(qu[...,2:3]*qu[...,3:4]+_P*qu[...,0:1]*qu[...,1:2]),
qq + 2.*qu[...,3:4]**2,
]).reshape(qu.shape[:-1]+(3,3))
return om
@ -1286,22 +1327,20 @@ class Rotation:
q12_s = qu[...,1:2]**2+qu[...,2:3]**2
chi = np.sqrt(q03_s*q12_s)
eu = np.where(np.abs(q12_s) < 1.0e-8,
np.block([np.arctan2(-_P*2.0*qu[...,0:1]*qu[...,3:4],qu[...,0:1]**2-qu[...,3:4]**2),
eu = np.where(np.abs(q12_s) < 1.e-8,
np.block([np.arctan2(-_P*2.*qu[...,0:1]*qu[...,3:4],qu[...,0:1]**2-qu[...,3:4]**2),
np.zeros(qu.shape[:-1]+(2,))]),
np.where(np.abs(q03_s) < 1.0e-8,
np.block([np.arctan2( 2.0*qu[...,1:2]*qu[...,2:3],qu[...,1:2]**2-qu[...,2:3]**2),
np.where(np.abs(q03_s) < 1.e-8,
np.block([np.arctan2( 2.*qu[...,1:2]*qu[...,2:3],qu[...,1:2]**2-qu[...,2:3]**2),
np.broadcast_to(np.pi,qu[...,0:1].shape),
np.zeros(qu.shape[:-1]+(1,))]),
np.block([np.arctan2((-_P*q02+q13)*chi, (-_P*q01-q23)*chi),
np.arctan2( 2.0*chi, q03_s-q12_s ),
np.arctan2( 2.*chi, q03_s-q12_s ),
np.arctan2(( _P*q02+q13)*chi, (-_P*q01+q23)*chi)])
)
)
# reduce Euler angles to definition range
eu[np.abs(eu)<1.e-6] = 0.0
eu = np.where(eu<0, (eu+2.0*np.pi)%np.array([2.0*np.pi,np.pi,2.0*np.pi]),eu) # needed?
return eu
eu[np.abs(eu) < 1.e-6] = 0.
return np.where(eu < 0., eu%(np.pi*np.array([2.,1.,2.])),eu)
@staticmethod
def _qu2ax(qu: np.ndarray) -> np.ndarray:
@ -1312,11 +1351,11 @@ class Rotation:
"""
with np.errstate(invalid='ignore',divide='ignore'):
s = np.sign(qu[...,0:1])/np.sqrt(qu[...,1:2]**2+qu[...,2:3]**2+qu[...,3:4]**2)
omega = 2.0 * np.arccos(np.clip(qu[...,0:1],-1.0,1.0))
ax = np.where(np.broadcast_to(qu[...,0:1] < 1.0e-8,qu.shape),
omega = 2. * np.arccos(np.clip(qu[...,0:1],-1.,1.))
ax = np.where(np.broadcast_to(qu[...,0:1] < 1.e-8,qu.shape),
np.block([qu[...,1:4],np.broadcast_to(np.pi,qu[...,0:1].shape)]),
np.block([qu[...,1:4]*s,omega]))
ax[np.isclose(qu[...,0],1.,rtol=0.0)] = [0.0, 0.0, 1.0, 0.0]
ax[np.isclose(qu[...,0],1.,rtol=0.)] = np.array([0.,0.,1.,0.])
return ax
@staticmethod
@ -1324,23 +1363,23 @@ class Rotation:
"""Quaternion to Rodrigues–Frank vector."""
with np.errstate(invalid='ignore',divide='ignore'):
s = np.linalg.norm(qu[...,1:4],axis=-1,keepdims=True)
ro = np.where(np.broadcast_to(np.abs(qu[...,0:1]) < 1.0e-12,qu.shape),
ro = np.where(np.broadcast_to(np.abs(qu[...,0:1]) < 1.e-12,qu.shape),
np.block([qu[...,1:2], qu[...,2:3], qu[...,3:4], np.broadcast_to(np.inf,qu[...,0:1].shape)]),
np.block([qu[...,1:2]/s,qu[...,2:3]/s,qu[...,3:4]/s,
np.tan(np.arccos(np.clip(qu[...,0:1],-1.0,1.0)))
np.tan(np.arccos(np.clip(qu[...,0:1],-1.,1.)))
])
)
ro[np.abs(s).squeeze(-1) < 1.0e-12] = [0.0,0.0,_P,0.0]
ro[np.abs(s).squeeze(-1) < 1.e-12] = np.array([0.,0.,_P,0.])
return ro
@staticmethod
def _qu2ho(qu: np.ndarray) -> np.ndarray:
"""Quaternion to homochoric vector."""
with np.errstate(invalid='ignore'):
omega = 2.0 * np.arccos(np.clip(qu[...,0:1],-1.0,1.0))
ho = np.where(np.abs(omega) < 1.0e-12,
omega = 2. * np.arccos(np.clip(qu[...,0:1],-1.,1.))
ho = np.where(np.abs(omega) < 1.e-12,
np.zeros(3),
qu[...,1:4]/np.linalg.norm(qu[...,1:4],axis=-1,keepdims=True) \
qu[...,1:4]/np.linalg.norm(qu[...,1:4],axis=-1,keepdims=True)
* (0.75*(omega - np.sin(omega)))**(1./3.))
return ho
@ -1362,12 +1401,12 @@ class Rotation:
trace = om[...,0,0:1] + om[...,1,1:2] + om[...,2,2:3]
with np.errstate(invalid='ignore',divide='ignore'):
s = [
0.5 / np.sqrt( 1.0 + trace),
2.0 * np.sqrt( 1.0 + om[...,0,0:1] - om[...,1,1:2] - om[...,2,2:3]),
2.0 * np.sqrt( 1.0 + om[...,1,1:2] - om[...,2,2:3] - om[...,0,0:1]),
2.0 * np.sqrt( 1.0 + om[...,2,2:3] - om[...,0,0:1] - om[...,1,1:2] )
]
s = np.array([
0.5 / np.sqrt( 1. + trace),
2. * np.sqrt( 1. + om[...,0,0:1] - om[...,1,1:2] - om[...,2,2:3]),
2. * np.sqrt( 1. + om[...,1,1:2] - om[...,2,2:3] - om[...,0,0:1]),
2. * np.sqrt( 1. + om[...,2,2:3] - om[...,0,0:1] - om[...,1,1:2] )
])
qu = np.where(trace>0,
np.block([0.25 / s[0],
(om[...,2,1:2] - om[...,1,2:3] ) * s[0],
@ -1389,18 +1428,18 @@ class Rotation:
0.25 * s[3]]),
)
)
)*np.array([1,_P,_P,_P])
qu[qu[...,0]<0] *=-1
)*np.array([1.,_P,_P,_P])
qu[qu[...,0] < 0.] *= -1.
return qu
@staticmethod
def _om2eu(om: np.ndarray) -> np.ndarray:
"""Rotation matrix to Bunge Euler angles."""
with np.errstate(invalid='ignore',divide='ignore'):
zeta = 1.0/np.sqrt(1.0-om[...,2,2:3]**2)
eu = np.where(np.isclose(np.abs(om[...,2,2:3]),1.0,0.0),
zeta = 1./np.sqrt(1.-om[...,2,2:3]**2)
eu = np.where(np.isclose(np.abs(om[...,2,2:3]),1.,0.),
np.block([np.arctan2(om[...,0,1:2],om[...,0,0:1]),
np.pi*0.5*(1-om[...,2,2:3]),
np.pi*0.5*(1.-om[...,2,2:3]),
np.zeros(om.shape[:-2]+(1,)),
]),
np.block([np.arctan2(om[...,2,0:1]*zeta,-om[...,2,1:2]*zeta),
@ -1409,8 +1448,7 @@ class Rotation:
])
)
eu[np.abs(eu) < 1.e-8] = 0.0
eu = np.where(eu<0, (eu+2.0*np.pi)%np.array([2.0*np.pi,np.pi,2.0*np.pi]),eu)
return eu
return np.where(eu < 0., eu%(np.pi*np.array([2.,1.,2.])),eu)
@staticmethod
def _om2ax(om: np.ndarray) -> np.ndarray:
@ -1419,18 +1457,18 @@ class Rotation:
om[...,2,0:1]-om[...,0,2:3],
om[...,0,1:2]-om[...,1,0:1]
])
t = 0.5*(om.trace(axis2=-2,axis1=-1) -1.0).reshape(om.shape[:-2]+(1,))
t = 0.5*(om.trace(axis2=-2,axis1=-1) -1.).reshape(om.shape[:-2]+(1,))
w,vr = np.linalg.eig(om)
# mask duplicated real eigenvalues
w[np.isclose(w[...,0],1.0+0.0j),1:] = 0.
w[np.isclose(w[...,1],1.0+0.0j),2:] = 0.
w[np.isclose(w[...,0],1.+0.j),1:] = 0.
w[np.isclose(w[...,1],1.+0.j),2:] = 0.
vr = np.swapaxes(vr,-1,-2)
ax = np.where(np.abs(diag_delta)<1e-13,
np.real(vr[np.isclose(w,1.0+0.0j)]).reshape(om.shape[:-2]+(3,)),
np.abs(np.real(vr[np.isclose(w,1.0+0.0j)]).reshape(om.shape[:-2]+(3,))) \
ax = np.where(np.abs(diag_delta)<1.e-13,
np.real(vr[np.isclose(w,1.+0.j)]).reshape(om.shape[:-2]+(3,)),
np.abs(np.real(vr[np.isclose(w,1.+0.j)]).reshape(om.shape[:-2]+(3,)))
*np.sign(diag_delta))
ax = np.block([ax,np.arccos(np.clip(t,-1.0,1.0))])
ax[np.abs(ax[...,3])<1.e-8] = [ 0.0, 0.0, 1.0, 0.0]
ax = np.block([ax,np.arccos(np.clip(t,-1.,1.))])
ax[np.abs(ax[...,3]) < 1.e-8] = np.array([0.,0.,1.,0.])
return ax
@staticmethod
@ -1460,7 +1498,7 @@ class Rotation:
-_P*sPhi*np.cos(ee[...,0:1]-ee[...,2:3]),
-_P*sPhi*np.sin(ee[...,0:1]-ee[...,2:3]),
-_P*cPhi*np.sin(ee[...,0:1]+ee[...,2:3])])
qu[qu[...,0]<0.0]*=-1
qu[qu[...,0] < 0.] *= -1.
return qu
@staticmethod
@ -1478,7 +1516,7 @@ class Rotation:
-c[...,0:1]*s[...,1:2],
+c[...,1:2]
]).reshape(eu.shape[:-1]+(3,3))
om[np.abs(om)<1.e-12] = 0.0
om[np.abs(om) < 1.e-12] = 0.
return om
@staticmethod
@ -1488,16 +1526,16 @@ class Rotation:
sigma = 0.5*(eu[...,0:1]+eu[...,2:3])
delta = 0.5*(eu[...,0:1]-eu[...,2:3])
tau = np.linalg.norm(np.block([t,np.sin(sigma)]),axis=-1,keepdims=True)
alpha = np.where(np.abs(np.cos(sigma))<1.e-12,np.pi,2.0*np.arctan(tau/np.cos(sigma)))
alpha = np.where(np.abs(np.cos(sigma))<1.e-12,np.pi,2.*np.arctan(tau/np.cos(sigma)))
with np.errstate(invalid='ignore',divide='ignore'):
ax = np.where(np.broadcast_to(np.abs(alpha)<1.0e-12,eu.shape[:-1]+(4,)),
[0.0,0.0,1.0,0.0],
ax = np.where(np.broadcast_to(np.abs(alpha)<1.e-12,eu.shape[:-1]+(4,)),
[0.,0.,1.,0.],
np.block([-_P/tau*t*np.cos(delta),
-_P/tau*t*np.sin(delta),
-_P/tau* np.sin(sigma),
alpha
]))
ax[(alpha<0.0).squeeze()] *=-1
ax[(alpha<0.).squeeze()] *= -1.
return ax
@staticmethod
@ -1506,7 +1544,7 @@ class Rotation:
ax = Rotation._eu2ax(eu)
ro = np.block([ax[...,:3],np.tan(ax[...,3:4]*.5)])
ro[ax[...,3] >= np.pi,3] = np.inf
ro[np.abs(ax[...,3])<1.e-16] = [ 0.0, 0.0, _P, 0.0 ]
ro[np.abs(ax[...,3])<1.e-16] = np.array([0.,0.,_P,0.])
return ro
@staticmethod
@ -1526,7 +1564,7 @@ class Rotation:
"""Axis–angle pair to quaternion."""
c = np.cos(ax[...,3:4]*.5)
s = np.sin(ax[...,3:4]*.5)
qu = np.where(np.abs(ax[...,3:4])<1.e-6,[1.0, 0.0, 0.0, 0.0],np.block([c, ax[...,:3]*s]))
qu = np.where(np.abs(ax[...,3:4]) < 1.e-6,[1.,0.,0.,0.],np.block([c,ax[...,:3]*s]))
return qu
@staticmethod
@ -1544,7 +1582,7 @@ class Rotation:
omc*ax[...,0:1]*ax[...,2:3] + s*ax[...,1:2],
omc*ax[...,1:2]*ax[...,2:3] - s*ax[...,0:1],
c+omc*ax[...,2:3]**2]).reshape(ax.shape[:-1]+(3,3))
return om if _P < 0.0 else np.swapaxes(om,-1,-2)
return om if _P < 0. else np.swapaxes(om,-1,-2)
@staticmethod
def _ax2eu(ax: np.ndarray) -> np.ndarray:
@ -1559,15 +1597,14 @@ class Rotation:
np.inf,
np.tan(ax[...,3:4]*0.5))
])
ro[np.abs(ax[...,3])<1.e-6] = [.0,.0,_P,.0]
ro[np.abs(ax[...,3]) < 1.e-6] = np.array([.0,.0,_P,.0])
return ro
@staticmethod
def _ax2ho(ax: np.ndarray) -> np.ndarray:
"""Axis–angle pair to homochoric vector."""
f = (0.75 * ( ax[...,3:4] - np.sin(ax[...,3:4]) ))**(1.0/3.0)
ho = ax[...,:3] * f
return ho
f = (0.75 * ( ax[...,3:4] - np.sin(ax[...,3:4]) ))**(1./3.)
return ax[...,:3] * f
@staticmethod
def _ax2cu(ax: np.ndarray) -> np.ndarray:
@ -1598,16 +1635,15 @@ class Rotation:
ax = np.where(np.isfinite(ro[...,3:4]),
np.block([ro[...,0:3]*np.linalg.norm(ro[...,0:3],axis=-1,keepdims=True),2.*np.arctan(ro[...,3:4])]),
np.block([ro[...,0:3],np.broadcast_to(np.pi,ro[...,3:4].shape)]))
ax[np.abs(ro[...,3]) < 1.e-8] = np.array([ 0.0, 0.0, 1.0, 0.0 ])
ax[np.abs(ro[...,3]) < 1.e-8] = np.array([0.,0.,1.,0.])
return ax
@staticmethod
def _ro2ho(ro: np.ndarray) -> np.ndarray:
"""Rodrigues–Frank vector to homochoric vector."""
f = np.where(np.isfinite(ro[...,3:4]),2.0*np.arctan(ro[...,3:4]) -np.sin(2.0*np.arctan(ro[...,3:4])),np.pi)
ho = np.where(np.broadcast_to(np.sum(ro[...,0:3]**2.0,axis=-1,keepdims=True) < 1.e-8,ro[...,0:3].shape),
np.zeros(3), ro[...,0:3]* (0.75*f)**(1.0/3.0))
return ho
f = np.where(np.isfinite(ro[...,3:4]),2.*np.arctan(ro[...,3:4]) -np.sin(2.*np.arctan(ro[...,3:4])),np.pi)
return np.where(np.broadcast_to(np.sum(ro[...,0:3]**2,axis=-1,keepdims=True) < 1.e-8,ro[...,0:3].shape),
np.zeros(3), ro[...,0:3]* (0.75*f)**(1./3.))
@staticmethod
def _ro2cu(ro: np.ndarray) -> np.ndarray:
@ -1641,13 +1677,12 @@ class Rotation:
+9.528014229335313e-6, -5.660288876265125e-6, +1.2844901692764126e-6,
+1.1255185726258763e-6, -1.3834391419956455e-6, +7.513691751164847e-7,
-2.401996891720091e-7, +4.386887017466388e-8, -3.5917775353564864e-9])
hmag_squared = np.sum(ho**2.,axis=-1,keepdims=True)
hmag_squared = np.sum(ho**2,axis=-1,keepdims=True)
s = np.sum(tfit*hmag_squared**np.arange(len(tfit)),axis=-1,keepdims=True)
with np.errstate(invalid='ignore'):
ax = np.where(np.broadcast_to(np.abs(hmag_squared)<1.e-8,ho.shape[:-1]+(4,)),
[ 0.0, 0.0, 1.0, 0.0 ],
np.block([ho/np.sqrt(hmag_squared),2.0*np.arccos(np.clip(s,-1.0,1.0))]))
return ax
return np.where(np.broadcast_to(np.abs(hmag_squared)<1.e-8,ho.shape[:-1]+(4,)),
[0.,0.,1.,0.],
np.block([ho/np.sqrt(hmag_squared),2.*np.arccos(np.clip(s,-1.,1.))]))
@staticmethod
def _ho2ro(ho: np.ndarray) -> np.ndarray:
@ -1671,27 +1706,25 @@ class Rotation:
with np.errstate(invalid='ignore',divide='ignore'):
# inverse M_3
xyz2 = xyz3[...,0:2] * np.sqrt( 2.0*rs/(rs+np.abs(xyz3[...,2:3])) )
xyz2 = xyz3[...,0:2] * np.sqrt( 2.*rs/(rs+np.abs(xyz3[...,2:3])) )
qxy = np.sum(xyz2**2,axis=-1,keepdims=True)
q2 = qxy + np.max(np.abs(xyz2),axis=-1,keepdims=True)**2
sq2 = np.sqrt(q2)
q = (_beta/np.sqrt(2.0)/_R1) * np.sqrt(q2*qxy/(q2-np.max(np.abs(xyz2),axis=-1,keepdims=True)*sq2))
q = (_beta/np.sqrt(2.)/_R1) * np.sqrt(q2*qxy/(q2-np.max(np.abs(xyz2),axis=-1,keepdims=True)*sq2))
tt = np.clip((np.min(np.abs(xyz2),axis=-1,keepdims=True)**2\
+np.max(np.abs(xyz2),axis=-1,keepdims=True)*sq2)/np.sqrt(2.0)/qxy,-1.0,1.0)
+np.max(np.abs(xyz2),axis=-1,keepdims=True)*sq2)/np.sqrt(2.)/qxy,-1.,1.)
T_inv = np.where(np.abs(xyz2[...,1:2]) <= np.abs(xyz2[...,0:1]),
np.block([np.ones_like(tt),np.arccos(tt)/np.pi*12.0]),
np.block([np.arccos(tt)/np.pi*12.0,np.ones_like(tt)]))*q
T_inv[xyz2<0.0] *= -1.0
T_inv[np.broadcast_to(np.isclose(qxy,0.0,rtol=0.0,atol=1.0e-12),T_inv.shape)] = 0.0
cu = np.block([T_inv, np.where(xyz3[...,2:3]<0.0,-np.ones_like(xyz3[...,2:3]),np.ones_like(xyz3[...,2:3])) \
* rs/np.sqrt(6.0/np.pi),
np.block([np.ones_like(tt),np.arccos(tt)/np.pi*12.]),
np.block([np.arccos(tt)/np.pi*12.,np.ones_like(tt)]))*q
T_inv[xyz2<0.] *= -1.
T_inv[np.broadcast_to(np.isclose(qxy,0.,rtol=0.,atol=1.e-12),T_inv.shape)] = 0.
cu = np.block([T_inv, np.where(xyz3[...,2:3]<0.,-np.ones_like(xyz3[...,2:3]),np.ones_like(xyz3[...,2:3])) \
* rs/np.sqrt(6./np.pi),
])/ _sc
cu[np.isclose(np.sum(np.abs(ho),axis=-1),0.0,rtol=0.0,atol=1.0e-16)] = 0.0
cu = np.take_along_axis(cu,Rotation._get_pyramid_order(ho,'backward'),-1)
return cu
cu[np.isclose(np.sum(np.abs(ho),axis=-1),0.,rtol=0.,atol=1.e-16)] = 0.
return np.take_along_axis(cu,Rotation._get_pyramid_order(ho,'backward'),-1)
#---------- Cubochoric ----------
@staticmethod
@ -1734,32 +1767,30 @@ class Rotation:
# get pyramid and scale by grid parameter ratio
XYZ = np.take_along_axis(cu,Rotation._get_pyramid_order(cu,'forward'),-1) * _sc
order = np.abs(XYZ[...,1:2]) <= np.abs(XYZ[...,0:1])
q = np.pi/12.0 * np.where(order,XYZ[...,1:2],XYZ[...,0:1]) \
q = np.pi/12. * np.where(order,XYZ[...,1:2],XYZ[...,0:1]) \
/ np.where(order,XYZ[...,0:1],XYZ[...,1:2])
c = np.cos(q)
s = np.sin(q)
q = _R1*2.0**0.25/_beta/ np.sqrt(np.sqrt(2.0)-c) \
q = _R1*2.**0.25/_beta/ np.sqrt(np.sqrt(2.)-c) \
* np.where(order,XYZ[...,0:1],XYZ[...,1:2])
T = np.block([ (np.sqrt(2.0)*c - 1.0), np.sqrt(2.0) * s]) * q
T = np.block([(np.sqrt(2.)*c - 1.), np.sqrt(2.) * s]) * q
# transform to sphere grid (inverse Lambert)
c = np.sum(T**2,axis=-1,keepdims=True)
s = c * np.pi/24.0 /XYZ[...,2:3]**2
c = c * np.sqrt(np.pi/24.0)/XYZ[...,2:3]
q = np.sqrt( 1.0 - s)
s = c * np.pi/24. /XYZ[...,2:3]**2
c = c * np.sqrt(np.pi/24.)/XYZ[...,2:3]
q = np.sqrt( 1. - s)
ho = np.where(np.isclose(np.sum(np.abs(XYZ[...,0:2]),axis=-1,keepdims=True),0.0,rtol=0.0,atol=1.0e-16),
np.block([np.zeros_like(XYZ[...,0:2]),np.sqrt(6.0/np.pi) *XYZ[...,2:3]]),
ho = np.where(np.isclose(np.sum(np.abs(XYZ[...,0:2]),axis=-1,keepdims=True),0.,rtol=0.,atol=1.e-16),
np.block([np.zeros_like(XYZ[...,0:2]),np.sqrt(6./np.pi)*XYZ[...,2:3]]),
np.block([np.where(order,T[...,0:1],T[...,1:2])*q,
np.where(order,T[...,1:2],T[...,0:1])*q,
np.sqrt(6.0/np.pi) * XYZ[...,2:3] - c])
np.sqrt(6./np.pi) * XYZ[...,2:3] - c])
)
ho[np.isclose(np.sum(np.abs(cu),axis=-1),0.0,rtol=0.0,atol=1.0e-16)] = 0.0
ho = np.take_along_axis(ho,Rotation._get_pyramid_order(cu,'backward'),-1)
return ho
ho[np.isclose(np.sum(np.abs(cu),axis=-1),0.,rtol=0.,atol=1.e-16)] = 0.
return np.take_along_axis(ho,Rotation._get_pyramid_order(cu,'backward'),-1)
@staticmethod

18
python/damask/_table.py
View File

@ -1,6 +1,6 @@
import re
import copy
from typing import Union, Tuple, List, Iterable
from typing import Optional, Union, Tuple, List, Iterable
import pandas as pd
import numpy as np
@ -13,8 +13,8 @@ class Table:
def __init__(self,
shapes: dict = {},
data: np.ndarray = None,
comments: Union[str, Iterable[str]] = None):
data: Optional[np.ndarray] = None,
comments: Union[None, str, Iterable[str]] = None):
"""
New spreadsheet.
@ -25,7 +25,7 @@ class Table:
For instance, 'F':(3,3) for a deformation gradient, or 'r':(1,) for a scalar.
data : numpy.ndarray or pandas.DataFrame, optional
Data. Existing column labels of a pandas.DataFrame will be replaced.
comments : str or iterable of str, optional
comments : (iterable of) str, optional
Additional, human-readable information.
"""
@ -188,7 +188,7 @@ class Table:
def _add_comment(self,
label: str,
shape: Tuple[int, ...],
info: str = None):
info: Optional[str] = None):
if info is not None:
specific = f'{label}{" "+str(shape) if np.prod(shape,dtype=np.int64) > 1 else ""}: {info}'
general = util.execution_stamp('Table')
@ -383,7 +383,7 @@ class Table:
def set(self,
label: str,
data: np.ndarray,
info: str = None) -> 'Table':
info: Optional[str] = None) -> 'Table':
"""
Add new or replace existing column data.
@ -458,15 +458,15 @@ class Table:
def rename(self,
old: Union[str, Iterable[str]],
new: Union[str, Iterable[str]],
info: str = None) -> 'Table':
info: Optional[str] = None) -> 'Table':
"""
Rename column data.
Parameters
----------
label_old : str or iterable of str
label_old : (iterable of) str
Old column label(s).
label_new : str or iterable of str
label_new : (iterable of) str
New column label(s).
Returns

25
python/damask/_vtk.py
View File

@ -1,7 +1,7 @@
import os
import multiprocessing as mp
from pathlib import Path
from typing import Union, Literal, List, Sequence
from typing import Optional, Union, Literal, List, Sequence
import numpy as np
import vtk
@ -110,13 +110,16 @@ class VTK:
Parameters
----------
comments : str or sequence of str
comments : (sequence of) str
Comments.
"""
s = vtk.vtkStringArray()
s.SetName('comments')
for c in util.tail_repack(comments,self.comments):
comments_ = util.tail_repack(comments,self.comments) if comments[:len(self.comments)] == self.comments else \
[comments] if isinstance(comments,str) else \
comments
for c in comments_:
s.InsertNextValue(c)
self.vtk_data.GetFieldData().AddArray(s)
@ -286,7 +289,7 @@ class VTK:
@staticmethod
def load(fname: Union[str, Path],
dataset_type: Literal['ImageData', 'UnstructuredGrid', 'PolyData', 'RectilinearGrid'] = None) -> 'VTK':
dataset_type: Literal[None, 'ImageData', 'UnstructuredGrid', 'PolyData', 'RectilinearGrid'] = None) -> 'VTK':
"""
Load from VTK file.
@ -409,11 +412,11 @@ class VTK:
# Check https://blog.kitware.com/ghost-and-blanking-visibility-changes/ for missing data
def set(self,
label: str = None,
data: Union[np.ndarray, np.ma.MaskedArray] = None,
info: str = None,
label: Optional[str] = None,
data: Union[None, np.ndarray, np.ma.MaskedArray] = None,
info: Optional[str] = None,
*,
table: 'Table' = None):
table: Optional['Table'] = None):
"""
Add new or replace existing point or cell data.
@ -533,7 +536,7 @@ class VTK:
def show(self,
label: str = None,
label: Optional[str] = None,
colormap: Union[Colormap, str] = 'cividis'):
"""
Render.
@ -547,9 +550,11 @@ class VTK:
Notes
-----
See http://compilatrix.com/article/vtk-1 for further ideas.
The first component is shown when visualizing vector datasets
(this includes tensor datasets because they are flattened).
"""
# See http://compilatrix.com/article/vtk-1 for possible improvements.
try:
import wx
_ = wx.App(False) # noqa

13
python/damask/seeds.py
View File

@ -1,7 +1,6 @@
"""Functionality for generation of seed points for Voronoi or Laguerre tessellation."""
from typing import Tuple as _Tuple
from typing import Optional as _Optional, Tuple as _Tuple
from scipy import spatial as _spatial
import numpy as _np
@ -14,8 +13,8 @@ from . import grid_filters as _grid_filters
def from_random(size: _FloatSequence,
N_seeds: int,
cells: _IntSequence = None,
rng_seed: _NumpyRngSeed = None) -> _np.ndarray:
cells: _Optional[_IntSequence] = None,
rng_seed: _Optional[_NumpyRngSeed] = None) -> _np.ndarray:
"""
Place seeds randomly in space.
@ -55,7 +54,7 @@ def from_Poisson_disc(size: _FloatSequence,
N_candidates: int,
distance: float,
periodic: bool = True,
rng_seed: _NumpyRngSeed = None) -> _np.ndarray:
rng_seed: _Optional[_NumpyRngSeed] = None) -> _np.ndarray:
"""
Place seeds according to a Poisson disc distribution.
@ -107,7 +106,7 @@ def from_Poisson_disc(size: _FloatSequence,
def from_grid(grid,
selection: _IntCollection = None,
selection: _Optional[_IntCollection] = None,
invert_selection: bool = False,
average: bool = False,
periodic: bool = True) -> _Tuple[_np.ndarray, _np.ndarray]:
@ -118,7 +117,7 @@ def from_grid(grid,
----------
grid : damask.Grid
Grid from which the material IDs are used as seeds.
selection : int or collection of int, optional
selection : (collection of) int, optional
Material IDs to consider.
invert_selection : bool, optional
Consider all material IDs except those in selection. Defaults to False.

2
python/damask/solver/_marc.py
View File

@ -3,7 +3,7 @@ import shlex
import re
from pathlib import Path
_marc_version = '2022.1'
_marc_version = '2022.2'
_marc_root = '/opt/msc'
_damask_root = str(Path(__file__).parents[3])

163
python/damask/util.py
View File

@ -10,8 +10,9 @@ import signal as _signal
import fractions as _fractions
from collections import abc as _abc
from functools import reduce as _reduce, partial as _partial
from typing import Callable as _Callable, Union as _Union, Iterable as _Iterable, Sequence as _Sequence, Dict as _Dict, \
List as _List, Tuple as _Tuple, Literal as _Literal, Any as _Any, Collection as _Collection, TextIO as _TextIO
from typing import Optional as _Optional, Callable as _Callable, Union as _Union, Iterable as _Iterable, \
Sequence as _Sequence, Dict as _Dict, List as _List, Tuple as _Tuple, Literal as _Literal, \
Any as _Any, Collection as _Collection, TextIO as _TextIO
from pathlib import Path as _Path
import numpy as _np
@ -40,29 +41,33 @@ _colors = {
# Functions
####################################################################################################
def srepr(msg,
glue: str = '\n') -> str:
glue: str = '\n',
quote: bool = False) -> str:
r"""
Join items with glue string.
Join (quoted) items with glue string.
Parameters
----------
msg : object with __repr__ or sequence of objects with __repr__
msg : (sequence of) object with __repr__
Items to join.
glue : str, optional
Glue used for joining operation. Defaults to '\n'.
quote : bool, optional
Quote items. Defaults to False.
Returns
-------
joined : str
String representation of the joined items.
String representation of the joined and quoted items.
"""
q = '"' if quote else ''
if (not hasattr(msg, 'strip') and
(hasattr(msg, '__getitem__') or
hasattr(msg, '__iter__'))):
return glue.join(str(x) for x in msg)
return glue.join(q+str(x)+q for x in msg)
else:
return msg if isinstance(msg,str) else repr(msg)
return q+(msg if isinstance(msg,str) else repr(msg))+q
def emph(msg) -> str:
@ -71,7 +76,7 @@ def emph(msg) -> str:
Parameters
----------
msg : object with __repr__ or sequence of objects with __repr__
msg : (sequence of) object with __repr__
Message to format.
Returns
@ -88,7 +93,7 @@ def deemph(msg) -> str:
Parameters
----------
msg : object with __repr__ or sequence of objects with __repr__
msg : (sequence of) object with __repr__
Message to format.
Returns
@ -105,7 +110,7 @@ def warn(msg) -> str:
Parameters
----------
msg : object with __repr__ or sequence of objects with __repr__
msg : (sequence of) object with __repr__
Message to format.
Returns
@ -122,7 +127,7 @@ def strikeout(msg) -> str:
Parameters
----------
msg : object with __repr__ or iterable of objects with __repr__
msg : (iterable of) object with __repr__
Message to format.
Returns
@ -136,8 +141,8 @@ def strikeout(msg) -> str:
def run(cmd: str,
wd: str = './',
env: _Dict[str, str] = None,
timeout: int = None) -> _Tuple[str, str]:
env: _Optional[_Dict[str, str]] = None,
timeout: _Optional[int] = None) -> _Tuple[str, str]:
"""
Run a command.
@ -210,6 +215,14 @@ def open_text(fname: _FileHandle,
open(_Path(fname).expanduser(),mode,newline=('\n' if mode == 'w' else None))
def execution_stamp(class_name: str,
function_name: _Optional[str] = None) -> str:
"""Timestamp the execution of a (function within a) class."""
now = _datetime.datetime.now().astimezone().strftime('%Y-%m-%d %H:%M:%S%z')
_function_name = '' if function_name is None else f'.{function_name}'
return f'damask.{class_name}{_function_name} v{_version} ({now})'
def natural_sort(key: str) -> _List[_Union[int, str]]:
"""
Natural sort.
@ -226,7 +239,7 @@ def natural_sort(key: str) -> _List[_Union[int, str]]:
def show_progress(iterable: _Iterable,
N_iter: int = None,
N_iter: _Optional[int] = None,
prefix: str = '',
bar_length: int = 50) -> _Any:
"""
@ -403,36 +416,35 @@ def project_equal_area(vector: _np.ndarray,
return _np.roll(_np.block([v[...,:2]/_np.sqrt(1.0+_np.abs(v[...,2:3])),_np.zeros_like(v[...,2:3])]),
-shift if keepdims else 0,axis=-1)[...,:3 if keepdims else 2]
def execution_stamp(class_name: str,
function_name: str = None) -> str:
"""Timestamp the execution of a (function within a) class."""
now = _datetime.datetime.now().astimezone().strftime('%Y-%m-%d %H:%M:%S%z')
_function_name = '' if function_name is None else f'.{function_name}'
return f'damask.{class_name}{_function_name} v{_version} ({now})'
def hybrid_IA(dist: _np.ndarray,
def hybrid_IA(dist: _FloatSequence,
N: int,
rng_seed: _NumpyRngSeed = None) -> _np.ndarray:
rng_seed: _Optional[_NumpyRngSeed] = None) -> _np.ndarray:
"""
Hybrid integer approximation.
Parameters
----------
dist : numpy.ndarray
Distribution to be approximated
Distribution to be approximated.
N : int
Number of samples to draw.
rng_seed : {None, int, array_like[ints], SeedSequence, BitGenerator, Generator}, optional
A seed to initialize the BitGenerator. Defaults to None.
If None, then fresh, unpredictable entropy will be pulled from the OS.
Returns
-------
hist : numpy.ndarray, shape (N)
Integer approximation of the distribution.
"""
N_opt_samples,N_inv_samples = (max(_np.count_nonzero(dist),N),0) # random subsampling if too little samples requested
N_opt_samples = max(_np.count_nonzero(dist),N) # random subsampling if too little samples requested
N_inv_samples = 0
scale_,scale,inc_factor = (0.0,float(N_opt_samples),1.0)
while (not _np.isclose(scale, scale_)) and (N_inv_samples != N_opt_samples):
repeats = _np.rint(scale*dist).astype(_np.int64)
repeats = _np.rint(scale*_np.array(dist)).astype(_np.int64)
N_inv_samples = _np.sum(repeats)
scale_,scale,inc_factor = (scale,scale+inc_factor*0.5*(scale - scale_), inc_factor*2.0) \
if N_inv_samples < N_opt_samples else \
@ -528,37 +540,82 @@ def shapeblender(a: _Tuple[int, ...],
return a + b[i:]
def extend_docstring(extra_docstring: str) -> _Callable:
def _docstringer(docstring: _Union[str, _Callable],
extra_parameters: _Optional[str] = None,
# extra_examples: _Optional[str] = None,
# extra_notes: _Optional[str] = None,
return_type: _Union[None, str, _Callable] = None) -> str:
"""
Decorator: Append to function's docstring.
Extend a docstring.
Parameters
----------
extra_docstring : str
Docstring to append.
docstring : str or callable, optional
Docstring (of callable) to extend.
extra_parameters : str, optional
Additional information to append to Parameters section.
return_type : str or callable, optional
Type of return variable.
"""
def _decorator(func):
func.__doc__ += extra_docstring
return func
return _decorator
docstring_ = str( docstring if isinstance(docstring,str)
else docstring.__doc__ if hasattr(docstring,'__doc__')
else '')
d = dict(Parameters=extra_parameters,
# Examples=extra_examples,
# Notes=extra_notes,
)
for key,extra in [(k,v) for (k,v) in d.items() if v is not None]:
if not (heading := _re.search(fr'^([ ]*){key}\s*\n\1{"-"*len(key)}',
docstring_,flags=_re.MULTILINE)):
raise RuntimeError(f"Docstring {docstring_} lacks a correctly formatted {key} section to insert values into")
content = [line for line in extra.split('\n') if line.strip()]
indent = len(heading.group(1))
shift = min([len(line)-len(line.lstrip(' '))-indent for line in content])
extra = '\n'.join([(line[shift:] if shift > 0 else
f'{" "*-shift}{line}') for line in content])
docstring_ = _re.sub(fr'(^([ ]*){key}\s*\n\2{"-"*len(key)}[\n ]*[A-Za-z0-9_ ]*: ([^\n]+\n)*)',
fr'\1{extra}\n',
docstring_,flags=_re.MULTILINE)
if return_type is None:
return docstring_
else:
if isinstance(return_type,str):
return_type_ = return_type
else:
return_class = return_type.__annotations__.get('return','')
return_type_ = (_sys.modules[return_type.__module__].__name__.split('.')[0]
+'.'
+(return_class.__name__ if not isinstance(return_class,str) else return_class)
)
def extended_docstring(f: _Callable,
extra_docstring: str) -> _Callable:
return _re.sub(r'(^([ ]*)Returns\s*\n\2-------\s*\n[ ]*[A-Za-z0-9_ ]*: )(.*)\n',
fr'\1{return_type_}\n',
docstring_,flags=_re.MULTILINE)
def extend_docstring(docstring: _Union[None, str, _Callable] = None,
extra_parameters: _Optional[str] = None) -> _Callable:
"""
Decorator: Combine another function's docstring with a given docstring.
Decorator: Extend the function's docstring.
Parameters
----------
f : function
Function of which the docstring is taken.
extra_docstring : str
Docstring to append.
docstring : str or callable, optional
Docstring to extend. Defaults to that of decorated function.
extra_parameters : str, optional
Additional information to append to Parameters section.
Notes
-----
Return type will become own type if docstring is callable.
"""
def _decorator(func):
func.__doc__ = f.__doc__ + extra_docstring
func.__doc__ = _docstringer(func.__doc__ if docstring is None else docstring,
extra_parameters,
func if isinstance(docstring,_Callable) else None,
)
return func
return _decorator
@ -622,8 +679,8 @@ def DREAM3D_cell_data_group(fname: _Union[str, _Path]) -> str:
def Bravais_to_Miller(*,
uvtw: _np.ndarray = None,
hkil: _np.ndarray = None) -> _np.ndarray:
uvtw: _Optional[_np.ndarray] = None,
hkil: _Optional[_np.ndarray] = None) -> _np.ndarray:
"""
Transform 4 Miller–Bravais indices to 3 Miller indices of crystal direction [uvw] or plane normal (hkl).
@ -649,10 +706,9 @@ def Bravais_to_Miller(*,
[0,0,0,1]]))
return _np.einsum('il,...l',basis,axis)
def Miller_to_Bravais(*,
uvw: _np.ndarray = None,
hkl: _np.ndarray = None) -> _np.ndarray:
uvw: _Optional[_np.ndarray] = None,
hkl: _Optional[_np.ndarray] = None) -> _np.ndarray:
"""
Transform 3 Miller indices to 4 Miller–Bravais indices of crystal direction [uvtw] or plane normal (hkil).
@ -706,7 +762,6 @@ def dict_prune(d: _Dict) -> _Dict:
return new
def dict_flatten(d: _Dict) -> _Dict:
"""
Recursively remove keys of single-entry dictionaries.
@ -738,7 +793,7 @@ def tail_repack(extended: _Union[str, _Sequence[str]],
Parameters
----------
extended : str or list of str
extended : (sequence of) str
Extended string list with potentially autosplitted tailing string relative to `existing`.
existing : list of str
Base string list.
@ -756,9 +811,9 @@ def tail_repack(extended: _Union[str, _Sequence[str]],
['a','new','shiny','e','n','t','r','y']
"""
return [extended] if isinstance(extended,str) else existing + \
([''.join(extended[len(existing):])] if _np.prod([len(i) for i in extended[len(existing):]]) == 1 else
list(extended[len(existing):]))
new = extended[len(existing):]
return [extended] if isinstance(extended,str) else \
existing + list([''.join(new)] if _np.prod([len(i) for i in new]) == 1 else new)
def aslist(arg: _Union[_IntCollection, int, None]) -> _List:
@ -767,7 +822,7 @@ def aslist(arg: _Union[_IntCollection, int, None]) -> _List:
Parameters
----------
arg : int or collection of int or None
arg : (collection of) int or None
Entity to transform into list.
Returns

2
python/tests/reference/ConfigMaterial/measured.material.yaml vendored
View File

@ -1,4 +1,4 @@
phase: {'1': t.b.d., '2': t.b.d.}
phase: {'1': null, '2': null}
homogenization:
direct: {N_constituents: 1}

30
python/tests/test_ConfigMaterial.py
View File

@ -65,17 +65,6 @@ class TestConfigMaterial:
del material_config['material'][0]['homogenization']
assert not material_config.is_complete
def test_incomplete_homogenization_N_constituents(self,ref_path):
material_config = ConfigMaterial.load(ref_path/'material.yaml')
for h in material_config['homogenization'].keys():
del material_config['homogenization'][h]['N_constituents']
assert not material_config.is_complete
def test_incomplete_phase_lattice(self,ref_path):
material_config = ConfigMaterial.load(ref_path/'material.yaml')
del material_config['phase']['Aluminum']['lattice']
assert not material_config.is_complete
def test_incomplete_wrong_phase(self,ref_path):
material_config = ConfigMaterial.load(ref_path/'material.yaml')
new = material_config.material_rename_phase({'Steel':'FeNbC'})
@ -86,6 +75,16 @@ class TestConfigMaterial:
new = material_config.material_rename_homogenization({'Taylor':'isostrain'})
assert not new.is_complete
def test_empty_phase(self,ref_path):
material_config = ConfigMaterial.load(ref_path/'material.yaml')
material_config['phase'] = None
assert not material_config.is_complete
def test_empty_homogenization(self,ref_path):
material_config = ConfigMaterial.load(ref_path/'material.yaml')
material_config['homogenization'] = None
assert not material_config.is_complete
def test_from_table(self):
N = np.random.randint(3,10)
a = np.vstack((np.hstack((np.arange(N),np.arange(N)[::-1])),
@ -98,6 +97,15 @@ class TestConfigMaterial:
for i,m in enumerate(c['material']):
assert m['homogenization'] == 1 and (m['constituents'][0]['O'] == [1,0,1,1]).all()
def test_updated_dicts(self,ref_path):
m1 = ConfigMaterial().material_add(phase=['Aluminum'],O=[1.0,0.0,0.0,0.0],homogenization='SX')
m2 = ConfigMaterial.load(ref_path/'material.yaml')
for k in m2['phase']:
m2 = m2.material_add(phase=[k],O=[1.0,0.0,0.0,0.0],homogenization='SX')
assert not m2['phase'].get(k) is None
assert m1['phase'].get('Aluminum') is None
assert m1['homogenization'].get('SX') is None
def test_from_table_with_constant(self):
N = np.random.randint(3,10)
a = np.vstack((np.hstack((np.arange(N),np.arange(N)[::-1])),

35
python/tests/test_Rotation.py
View File

@ -708,6 +708,7 @@ class TestRotation:
@pytest.mark.parametrize('degrees',[True,False])
def test_axis_angle(self,set_of_rotations,degrees,normalize,P):
c = np.array([P*-1,P*-1,P*-1,1.])
c[:3] *= 0.9 if normalize else 1.0
for rot in set_of_rotations:
m = rot.as_Euler_angles()
o = Rotation.from_axis_angle(rot.as_axis_angle(degrees)*c,degrees,normalize,P).as_Euler_angles()
@ -730,16 +731,30 @@ class TestRotation:
assert ok and np.isclose(np.linalg.norm(o[:3]),1.0) \
and o[3]<=np.pi+1.e-9, f'{m},{o},{rot.as_quaternion()}'
def test_parallel(self,set_of_rotations):
a = np.array([[1.0,0.0,0.0],
[0.0,1.0,0.0]])
for rot in set_of_rotations:
assert rot.allclose(Rotation.from_parallel(a,rot.broadcast_to((2,))@a))
@pytest.mark.parametrize('P',[1,-1])
@pytest.mark.parametrize('normalize',[True,False])
def test_Rodrigues(self,set_of_rotations,normalize,P):
c = np.array([P*-1,P*-1,P*-1,1.])
c[:3] *= 0.9 if normalize else 1.0
for rot in set_of_rotations:
m = rot.as_matrix()
o = Rotation.from_Rodrigues_vector(rot.as_Rodrigues_vector()*c,normalize,P).as_matrix()
ok = np.allclose(m,o,atol=atol)
assert ok and np.isclose(np.linalg.det(o),1.0), f'{m},{o}'
def test_Rodrigues_compact(self,set_of_rotations):
for rot in set_of_rotations:
c = rot.as_Rodrigues_vector(compact=True)
r = rot.as_Rodrigues_vector(compact=False)
assert np.allclose(r[:3]*r[3], c, equal_nan=True)
@pytest.mark.parametrize('P',[1,-1])
def test_homochoric(self,set_of_rotations,P):
cutoff = np.tan(np.pi*.5*(1.-1e-4))
@ -760,11 +775,12 @@ class TestRotation:
@pytest.mark.parametrize('P',[1,-1])
@pytest.mark.parametrize('accept_homomorph',[True,False])
def test_quaternion(self,set_of_rotations,P,accept_homomorph):
c = np.array([1,P*-1,P*-1,P*-1]) * (-1 if accept_homomorph else 1)
@pytest.mark.parametrize('normalize',[True,False])
def test_quaternion(self,set_of_rotations,P,accept_homomorph,normalize):
c = np.array([1,P*-1,P*-1,P*-1]) * (-1 if accept_homomorph else 1) * (0.9 if normalize else 1.0)
for rot in set_of_rotations:
m = rot.as_cubochoric()
o = Rotation.from_quaternion(rot.as_quaternion()*c,accept_homomorph,P).as_cubochoric()
o = Rotation.from_quaternion(rot.as_quaternion()*c,accept_homomorph,normalize,P).as_cubochoric()
ok = np.allclose(m,o,atol=atol)
if np.count_nonzero(np.isclose(np.abs(o),np.pi**(2./3.)*.5)):
ok |= np.allclose(m*-1.,o,atol=atol)
@ -889,6 +905,15 @@ class TestRotation:
with pytest.raises(ValueError):
fr(eval(f'R.{to}()'),P=-30)
def test_invalid_multiplication(self):
rot = Rotation.from_random()
with pytest.raises(TypeError):
rot@Rotation.from_random()
with pytest.raises(TypeError):
rot@[1,2,3,4]
@pytest.mark.parametrize('shape',[None,(3,),(4,2)])
def test_broadcast(self,shape):
rot = Rotation.from_random(shape)
@ -902,6 +927,7 @@ class TestRotation:
@pytest.mark.parametrize('function,invalid',[(Rotation.from_quaternion, np.array([-1,0,0,0])),
(Rotation.from_quaternion, np.array([1,1,1,0])),
(Rotation.from_Euler_angles, np.array([1,4,0])),
(Rotation.from_Euler_angles, np.array([-1,0,0])),
(Rotation.from_axis_angle, np.array([1,0,0,4])),
(Rotation.from_axis_angle, np.array([1,1,0,1])),
(Rotation.from_matrix, np.random.rand(3,3)),
@ -909,7 +935,8 @@ class TestRotation:
(Rotation.from_Rodrigues_vector, np.array([1,0,0,-1])),
(Rotation.from_Rodrigues_vector, np.array([1,1,0,1])),
(Rotation.from_homochoric, np.array([2,2,2])),
(Rotation.from_cubochoric, np.array([1.1,0,0])) ])
(Rotation.from_cubochoric, np.array([1.1,0,0])),
])
def test_invalid_value(self,function,invalid):
with pytest.raises(ValueError):
function(invalid)

143
python/tests/test_util.py
View File

@ -8,7 +8,6 @@ import h5py
from damask import util
class TestUtil:
@pytest.mark.xfail(sys.platform == 'win32', reason='echo is not a Windows command')
@ -44,7 +43,7 @@ class TestUtil:
@pytest.mark.parametrize('rv',[stats.rayleigh(),stats.weibull_min(1.2),stats.halfnorm(),stats.pareto(2.62)])
def test_hybridIA(self,rv):
def test_hybridIA_distribution(self,rv):
bins = np.linspace(0,10,100000)
centers = (bins[1:]+bins[:-1])/2
N_samples = bins.shape[0]-1000
@ -53,6 +52,21 @@ class TestUtil:
dist_sampled = np.histogram(centers[selected],bins)[0]/N_samples*np.sum(dist)
assert np.sqrt(((dist - dist_sampled) ** 2).mean()) < .025 and selected.shape[0]==N_samples
def test_hybridIA_constant(self):
N_bins = np.random.randint(20,400)
m = np.random.randint(1,20)
N_samples = m * N_bins
dist = np.ones(N_bins)*np.random.rand()
assert np.all(np.sort(util.hybrid_IA(dist,N_samples))==np.arange(N_samples).astype(int)//m)
def test_hybridIA_linear(self):
N_points = np.random.randint(10,200)
m = np.random.randint(1,20)
dist = np.arange(N_points)
N_samples = m * np.sum(dist)
assert np.all(np.bincount(util.hybrid_IA(dist*np.random.rand(),N_samples)) == dist*m)
@pytest.mark.parametrize('point,direction,normalize,keepdims,answer',
[
([1,0,0],'z',False,True, [1,0,0]),
@ -208,3 +222,128 @@ class TestUtil:
@pytest.mark.parametrize('kw_Miller,kw_Bravais',[('uvw','uvtw'),('hkl','hkil')])
def test_Bravais_Miller_Bravais(self,vector,kw_Miller,kw_Bravais):
assert np.all(vector == util.Miller_to_Bravais(**{kw_Miller:util.Bravais_to_Miller(**{kw_Bravais:vector})}))
@pytest.mark.parametrize('extra_parameters',["""
p2 : str, optional
p2 description 1
p2 description 2
""",
"""
p2 : str, optional
p2 description 1
p2 description 2
""",
"""
p2 : str, optional
p2 description 1
p2 description 2
"""])
@pytest.mark.parametrize('invalid_docstring',["""
Function description
Parameters ----------
p0 : numpy.ndarray, shape (...,4)
p0 description 1
p0 description 2
p1 : int, optional
p1 description
Remaining description
""",
"""
Function description
Parameters
----------
p0 : numpy.ndarray, shape (...,4)
p0 description 1
p0 description 2
p1 : int, optional
p1 description
Remaining description
""",])
def test_extend_docstring_parameters(self,extra_parameters,invalid_docstring):
test_docstring = """
Function description
Parameters
----------
p0 : numpy.ndarray, shape (...,4)
p0 description 1
p0 description 2
p1 : int, optional
p1 description
Remaining description
"""
invalid_docstring = """
Function description
Parameters ----------
p0 : numpy.ndarray, shape (...,4)
p0 description 1
p0 description 2
p1 : int, optional
p1 description
Remaining description
"""
expected = """
Function description
Parameters
----------
p0 : numpy.ndarray, shape (...,4)
p0 description 1
p0 description 2
p1 : int, optional
p1 description
p2 : str, optional
p2 description 1
p2 description 2
Remaining description
""".split("\n")
assert expected == util._docstringer(test_docstring,extra_parameters).split('\n')
with pytest.raises(RuntimeError):
util._docstringer(invalid_docstring,extra_parameters)
def test_replace_docstring_return_type(self):
class TestClassOriginal:
pass
def original_func() -> TestClassOriginal:
pass
class TestClassDecorated:
def decorated_func_bound(self) -> 'TestClassDecorated':
pass
def decorated_func() -> TestClassDecorated:
pass
original_func.__doc__ = """
Function description/Parameters
Returns
-------
Return value : test_util.TestClassOriginal
Remaining description
"""
expected = """
Function description/Parameters
Returns
-------
Return value : test_util.TestClassDecorated
Remaining description
"""
assert expected == util._docstringer(original_func,return_type=decorated_func)
assert expected == util._docstringer(original_func,return_type=TestClassDecorated.decorated_func_bound)

2
src/Marc/DAMASK_Marc.f90
View File

@ -413,6 +413,8 @@ subroutine uedinc(inc,incsub)
use discretization_Marc
implicit none(type,external)
external :: nodvar
integer(pI64), intent(in) :: inc, incsub
integer :: n, nqncomp, nqdatatype

466
src/Marc/include/concom2022.2 vendored Normal file
View File

@ -0,0 +1,466 @@
! common block definition file taken from respective MSC.Marc release and reformated to free format
!***********************************************************************
!
! File: concom.cmn
!
! MSC.Marc include file
!
integer &
iacous, iasmbl, iautth, ibear, icompl, iconj, icreep, ideva, idyn, idynt,&
ielas, ielcma, ielect, iform, ifour, iharm, ihcps, iheat, iheatt, ihresp,&
ijoule, ilem, ilnmom, iloren, inc, incext, incsub, ipass, iplres, ipois,&
ipoist, irpflo, ismall, ismalt, isoil, ispect, ispnow, istore, iswep, ithcrp,&
itherm, iupblg, iupdat, jacflg, jel, jparks, largst, lfond, loadup, loaduq,&
lodcor, lovl, lsub, magnet, ncycle, newtnt, newton, noshr, linear, ivscpl,&
icrpim, iradrt, ipshft, itshr, iangin, iupmdr, iconjf, jincfl, jpermg, jhour,&
isolvr, jritz, jtable, jshell, jdoubl, jform, jcentr, imini, kautth, iautof,&
ibukty, iassum, icnstd, icnstt, kmakmas, imethvp, iradrte, iradrtp, iupdate, iupdatp,&
ncycnt, marmen , idynme, ihavca, ispf, kmini, imixex, largtt, kdoela, iautofg,&
ipshftp, idntrc, ipore, jtablm, jtablc, isnecma, itrnspo, imsdif, jtrnspo, mcnear,&
imech, imecht, ielcmat, ielectt, magnett, imsdift, noplas, jtabls, jactch, jtablth,&
kgmsto , jpzo, ifricsh, iremkin, iremfor, ishearp, jspf, machining, jlshell, icompsol,&
iupblgfo, jcondir, nstcrp, nactive, ipassref, nstspnt, ibeart, icheckmpc, noline, icuring,&
ishrink, ioffsflg, isetoff, ioffsetm,iharmt, inc_incdat, iautspc, ibrake, icbush, istream_input,&
iprsinp, ivlsinp, ifirst_time,ipin_m, jgnstr_glb, imarc_return,iqvcinp, nqvceid, istpnx, imicro1,&
iaxisymm, jbreakglue,iglstif, jfastasm,iwear, iwearcf, imixmeth, ielcmadyn, idinout, igena_meth,&
magf_meth, non_assumed, iredoboudry, ioffsz0,icomplt, mesh_dual, iactrp, mgnewton, iusedens,igsigd0,&
iaem, icosim, inodels, nlharm, iampini, iphasetr, inonlcl, inonlct, iforminp,ispecerror,&
icsprg, imol, imolt, idatafit,iharmpar, inclcase, imultifreq,init_elas, ifatig, iftgmat,&
nchybrid, ibuckle, iexpande
dimension :: ideva(60)
integer num_concom
parameter(num_concom=263)
common/marc_concom/&
iacous, iasmbl, iautth, ibear, icompl, iconj, icreep, ideva, idyn, idynt,&
ielas, ielcma, ielect, iform, ifour, iharm, ihcps, iheat, iheatt, ihresp,&
ijoule, ilem, ilnmom, iloren, inc, incext, incsub, ipass, iplres, ipois,&
ipoist, irpflo, ismall, ismalt, isoil, ispect, ispnow, istore, iswep, ithcrp,&
itherm, iupblg, iupdat, jacflg, jel, jparks, largst, lfond, loadup, loaduq,&
lodcor, lovl, lsub, magnet, ncycle, newtnt, newton, noshr, linear, ivscpl,&
icrpim, iradrt, ipshft, itshr, iangin, iupmdr, iconjf, jincfl, jpermg, jhour,&
isolvr, jritz, jtable, jshell, jdoubl, jform, jcentr, imini, kautth, iautof,&
ibukty, iassum, icnstd, icnstt, kmakmas, imethvp, iradrte, iradrtp, iupdate, iupdatp,&
ncycnt, marmen, idynme, ihavca, ispf, kmini, imixex, largtt, kdoela, iautofg,&
ipshftp, idntrc, ipore, jtablm, jtablc, isnecma, itrnspo, imsdif, jtrnspo, mcnear,&
imech, imecht, ielcmat, ielectt, magnett, imsdift, noplas, jtabls, jactch, jtablth,&
kgmsto , jpzo, ifricsh, iremkin, iremfor, ishearp, jspf, machining, jlshell, icompsol,&
iupblgfo, jcondir, nstcrp, nactive, ipassref, nstspnt, ibeart, icheckmpc, noline, icuring,&
ishrink, ioffsflg, isetoff, ioffsetm,iharmt, inc_incdat, iautspc, ibrake, icbush, istream_input,&
iprsinp, ivlsinp, ifirst_time,ipin_m, jgnstr_glb, imarc_return,iqvcinp, nqvceid, istpnx, imicro1,&
iaxisymm, jbreakglue,iglstif, jfastasm,iwear, iwearcf, imixmeth, ielcmadyn, idinout, igena_meth,&
magf_meth, non_assumed, iredoboudry, ioffsz0,icomplt, mesh_dual, iactrp, mgnewton, iusedens,igsigd0,&
iaem, icosim, inodels, nlharm, iampini, iphasetr, inonlcl, inonlct, iforminp,ispecerror,&
icsprg, imol, imolt, idatafit,iharmpar, inclcase, imultifreq,init_elas, ifatig, iftgmat,&
nchybrid, ibuckle, iexpande
!
! comments of variables:
!
! iacous Control flag for acoustic analysis. Input data.
! iacous=1 modal acoustic analysis.
! iacous=2 harmonic acoustic-structural analysis.
! iasmbl Control flag to indicate that operator matrix should be
! recalculated.
! iautth Control flag for AUTO THERM option.
! ibear Control flag for bearing analysis. Input data.
! icompl Control variable to indicate that a complex analysis is
! being performed. Either a Harmonic analysis with damping,
! or a harmonic electro-magnetic analysis. Input data.
! iconj Flag for EBE conjugate gradient solver (=solver 1, retired)
! Also used for VKI iterative solver.
! icreep Control flag for creep analysis. Input data.
! ideva(60) - debug print out flag
! 1 print element stiffness matrices, mass matrix
! 2 output matrices used in tying
! 3 force the solution of a nonpositive definite matrix
! 4 print info of connections to each node
! 5 info of gap convergence, internal heat generated, contact
! touching and separation
! 6 nodal value array during rezoning
! 7 tying info in CONRAD GAP option, fluid element numbers in
! CHANNEL option
! 8 output incremental displacements in local coord. system
! 9 latent heat output
! 10 stress-strain in local coord. system
! 11 additional info on interlaminar stress
! 12 output right hand side and solution vector
! 13 info of CPU resources used and memory available on NT
! 14 info of mesh adaption process, 2D outline information
! info of penetration checking for remeshing
! save .fem files after afmesh3d meshing
! print local adaptivity info
! 15 surface energy balance flag
! 16 print info regarding pyrolysis
! 17 print info of "streamline topology"
! 18 print mesh data changes after remeshing
! 19 print material flow stress data read in from *.mat file
! if unit flag is on, print out flow stress after conversion
! 20 print information on table input
! 21 print out information regarding kinematic boundary conditions
! 22 print out information regarding dist loads, point loads, film
! and foundations
! 23 print out information about automatic domain decomposition
! 24 print out iteration information in SuperForm status report file
! 25 print out information for ablation
! 26 print out information for films - Table input
! 27 print out the tying forces
! 28 print out for CASI solver, convection,
! 29 DDM single file debug printout
! 30 print out cavity debug info
! 31 print out welding related info
! 32 prints categorized DDM memory usage
! 33 print out the cutting info regarding machining feature
! 34 print out the list of quantities which can be defined via a table
! and for each quantity the supported independent variables
! 35 print out detailed coupling region info
! 36 print out solver debug info level 1 (Least Detailed)
! 37 print out solver debug info level 1 (Medium Detailed)
! 38 print out solver debug info level 1 (Very Detailed)
! 39 print detailed memory allocation info
! 40 print out marc-adams debug info
! 41 output rezone mapping post file for debugging
! 42 output post file after calling oprofos() for debugging
! 43 debug printout for vcct
! 44 debug printout for progressive failure
! 45 print out automatically generated midside node coordinates (arecrd)
! 46 print out message about routine and location, where the ibort is raised (ibort_inc)
! 47 print out summary message of element variables on a
! group-basis after all the automatic changes have been
! made (em_ellibp)
! 48 Automatically generate check results based on max and min vals.
! These vals are stored in the checkr file, which is inserted
! into the *dat file by the generate_check_results script from /marc/tools
! 49 Automatically generate check results based on the real calculated values
! at the sppecified check result locations.
! These vals are stored in the checkr file, which is inserted
! into the *dat file by the update_check_results script from /marc/tools
! 50 generate a file containing the resistance or capacity matrix;
! this file can be used to compare results with a reference file
! 51 print out detailed information for segment-to-segment contact
! 52 print out detailed relative displacement information
! for uniaxial sliding contact
! 53 print out detailed sliding direction information for
! uniaxial sliding contact
! 54 print out detailed information for edges attached to a curve
! 55 print information related to viscoelasticity calculations
! 56 print out detailed information for element coloring for multithreading
! 57 print out extra overheads due to multi-threading.
! These overhead includes (i) time and (ii) memory.
! The memory report will be summed over all the children.
! 58 debug output for ELSTO usage
! 59 print out contact body forces and nodes in contact
!
! idyn Control flag for dynamics. Input data.
! 1 = eigenvalue extraction and / or modal superposition
! 2 = Newmark Beta and Single Step Houbolt (ssh with idynme=1)
! 3 = Houbolt
! 4 = Central difference
! 5 = Newer central difference
! idynt Copy of idyn at begining of increment
! ielas Control flag for ELASTIC analysis. Input data.
! Set by user or automatically turned on by Fourier option.
! Implies that each load case is treated separately.
! In Adaptive meshing analysis , forces re-analysis until
! convergence obtained.
! Also seriously misused to indicate no convergence.
! = 1 elastic option with fourier analysis
! = 2 elastic option without fourier analysis
! =-1 no convergence in recycles or max # increments reached
! Set to 1 if ELASTIC or SUBSTRUC parameter cards are used,
! or if fourier option is used.
! Then set to 2 if not fourier analysis.
! ielcma Control flag for electromagnetic analysis. Input data.
! ielcma = 1 Harmonic formulation
! ielcma = 2 Transient formulation
! ielect Control flag for electrostatic option. Input data.
! iform Control flag indicating that contact will be performed.
! ifour Control flag for Fourier analysis.
! 0 = Odd and even terms.
! 1 = symmetric (cosine) terms
! 2 = antisymmetric (sine) terms.
! iharm Control flag to indicate that a harmonic analysis will
! be performed. May change between passes.
! ihcps Control flag for coupled thermal - stress analysis.
! iheat Control flag for heat transfer analysis. Input data.
! iheatt Permanent control flag for heat transfer analysis.
! Note in coupled analysis iheatt will remain as one,
! but iheat will be zero in stress pass.
! ihresp Control flag to indicate to perform a harmonic subincrement.
! ijoule Control flag for Joule heating.
! ilem Control flag to determin which vector is to be transformed.
! Control flag to see where one is:
! ilem = 1 - elem.f
! ilem = 2 - initst.f
! ilem = 3 - pressr.f
! ilem = 3 - fstif.f
! ilem = 4 - jflux.f
! ilem = 4 - strass.f
! ilem = 5 - mass.f
! ilem = 5 - osolty.f
! ilnmom Control flag for soil - pore pressure calculation. Input data.
! ilnmom = 0 - perform only pore pressure calculation.
! = 1 - couples pore pressure - displacement analysis
! iloren Control flag for DeLorenzi J-Integral evaluation. Input data.
! inc Increment number.
! incext Control flag indicating that currently working on a
! subincrement.
! Could be due to harmonics , damping component (bearing),
! stiffness component (bearing), auto therm creep or
! old viscoplaticity
! incsub Sub-increment number.
! inonlcl control flag for nonlocal pass
! inonlct permanent control flag for nonlocal pass
! ipass Control flag for which part of coupled analysis.
! ipass = -1 - reset to base values
! ipass = 0 - do nothing
! ipass = 1 - stress part
! ipass = 2 - heat transfer part
! 3 - fluid pass
! 4 - joule heating pass
! 5 - pore pressure pass
! 6 - electrostatic pass
! 7 - magnetostatic pass
! 8 - electromagnetic pass
! 9 - diffusion pass
! ipass = 10 - nonlocal part
! iplres Flag indicating that either second matrix is stored.
! dynamic analysis - mass matrix
! heat transfer - specific heat matrix
! buckle - initial stress stiffness
! ipois Control flag indicating Poisson type analysis
! ipois = 1 for heat transfer
! = 1 for heat transfer part of coupled
! = 1 for bearing
! = 1 for electrostatic
! = 1 for magnetostatic
! = 1 for nonlocal part
! ipoist Permanent copy of ipois. In coupled analysis , ipois = 0
! in stress portion, yet ipoist will still =1.
! irpflo global flag for rigid plastic flow analysis
! = 1 eularian formulation
! = 2 regular formulation; rigid material present in the analysis
! ismall control flag to indicate small displacement analysis. input data.
! ismall = 0 - large disp included.
! ismall = 1 - small displacement.
! the flag is changing between passes.
! ismalt permanent copy of ismall . in heat transfer portion of
! coupled analysis ismall =0 , but ismalt remains the same.
! isoil control flag indicating that soil / pore pressure
! calculation . input data.
! ispect control flag for response spectrum calculation. input data.
! ispnow control flag to indicate to perform a spectrum response
! calculation now.
! istore store stresses flag.
! istore = 0 in elem.f and if first pass of creep
! convergence checking in ogetst.f
! or harmonic analysis or thruc.f if not
! converged.
! iswep control flag for eigenvalue analysis.
! iswep=1 - go do extraction process
! ithcrp control flag for auto therm creep option. input data.
! itherm control flag for either temperature dependent material
! properties and/or thermal loads.
! iupblg control flag for follower force option. input data.
! iupdat control flag for update lagrange option for current element.
! jacflg control flag for lanczos iteration method. input data.
! jel control flag indicating that total load applied in
! increment, ignore previous solution.
! jel = 1 in increment 0
! = 1 if elastic or fourier
! = 1 in subincrements with elastic and adaptive
! jparks control flag for j integral by parks method. input data.
! largst control flag for finite strain plasticity. input data.
! lfond control variable that indicates if doing elastic
! foundation or film calculation. influences whether
! this is volumetric or surface integration.
! loadup control flag that indicates that nonlinearity occurred
! during previous increment.
! loaduq control flag that indicates that nonlinearity occurred.
! lodcor control flag for switching on the residual load correction.
! notice in input stage lodcor=0 means no loadcor,
! after omarc lodcor=1 means no loadcor
! lovl control flag for determining which "overlay" is to
! be called from ellib.
! lovl = 1 omarc
! = 2 oaread
! = 3 opress
! = 4 oasemb
! = 5 osolty
! = 6 ogetst
! = 7 oscinc
! = 8 odynam
! = 9 opmesh
! = 10 omesh2
! = 11 osetz
! = 12 oass
! = 13 oincdt
! = 14 oasmas
! = 15 ofluas
! = 16 ofluso
! = 17 oshtra
! = 18 ocass
! = 19 osoltc
! = 20 orezon
! = 21 otest
! = 22 oeigen
! lsub control variable to determine which part of element
! assembly function is being done.
! lsub = 1 - no longer used
! = 2 - beta*
! = 3 - cons*
! = 4 - ldef*
! = 5 - posw*
! = 6 - theta*
! = 7 - tmarx*
! = 8 - geom*
! magnet control flag for magnetostatic analysis. input data.
! ncycle cycle number. accumulated in osolty.f
! note first time through oasemb.f , ncycle = 0.
! newtnt control flag for permanent copy of newton.
! newton iteration type. input data.
! newton : = 1 full newton raphson
! 2 modified newton raphson
! 3 newton raphson with strain correct.
! 4 direct substitution
! 5 direct substitution followed by n.r.
! 6 direct substitution with line search
! 7 full newton raphson with secant initial stress
! 8 secant method
! 9 full newton raphson with line search
! noshr control flag for calculation interlaminar shears for
! elements 22,45, and 75. input data.
!ees
!
! jactch = 1 or 2 if elements are activated or deactivated
! = 3 if elements are adaptively remeshed or rezoned
! = 0 normally / reset to 0 when assembly is done
! ifricsh = 0 call to fricsh in otest not needed
! = 1 call to fricsh (nodal friction) in otest needed
! iremkin = 0 remove deactivated kinematic boundary conditions
! immediately - only in new input format (this is default)
! = 1 remove deactivated kinematic boundary conditions
! gradually - only in new input format
! iremfor = 0 remove force boundary conditions immediately -
! only in new input format (this is default)
! = 1 remove force boundary conditions gradually -
! only in new input format (this is default)
! ishearp set to 1 if shear panel elements are present in the model
!
! jspf = 0 not in spf loadcase
! > 0 in spf loadcase (jspf=1 during first increment)
! machining = 1 if the metal cutting feature is used, for memory allocation purpose
! = 0 (default) if no metal cutting feature required
!
! jlshell = 1 if there is a shell element in the mesh
! icompsol = 1 if there is a composite solid element in the mesh
! iupblgfo = 1 if follower force for point loads
! jcondir = 1 if contact priority option is used
! nstcrp = 0 (default) steady state creep flag (undocumented feature.
! if not 0, turns off special ncycle = 0 code in radial.f)
! nactive = number of active passes, if =1 then it's not a coupled analysis
! ipassref = reference ipass, if not in a multiphysics pass ipass=ipassref
! icheckmpc = value of mpc-check parameter option
! noline = set to 1 in osolty if no line seacrh should be done in ogetst
! icuring = set to 1 if the curing is included for the heat transfer analysis.
! ishrink = set to 1 if shrinkage strain is included for mechancial analysis.
! ioffsflg = 1 for small displacement beam/shell offsets
! = 2 for large displacement beam/shell offsets
! isetoff = 0 - do not apply beam/shell offsets
! = 1 - apply beam/shell offsets
! ioffsetm = min. value of offset flag
! iharmt = 1 global flag if a coupled analysis contains an harmonic pass
! inc_incdat = flag to record increment number of a new loadcase in incdat.f
! iautspc = flag for AutoSPC option
! ibrake = brake squeal in this increment
! icbush = set to 1 if cbush elements present in model
! istream_input = set to 1 for streaming input calling Marc as library
! iprsinp = set to 1 if pressure input, introduced so other variables
! such as h could be a function of pressure
! ivlsinp = set to 1 if velocity input, introduced so other variables
! such as h could be a function of velocity
! ipin_m = # of beam element with PIN flag
! jgnstr_glb = global control over pre or fast integrated composite shells
! imarc_return = Marc return flag for streaming input control
! iqvcimp = if non-zero, then the number of QVECT boundary conditions
! nqvceid = number of QVECT boundary conditions, where emisivity/absorbtion id entered
! istpnx = 1 if to stop at end of increment
! imicro1 = 1 if micro1 interface is used
! iaxisymm = set to 1 if axisymmetric analysis
! jbreakglue = set to 1 if breaking glued option is used
! iglstif = 1 if ddm and global stiffness matrix formed (sgi solver 6 or solver9)
! jfastasm = 1 do fast assembly using SuperForm code
! iwear = set to 1 if wear model, set to 2 if wear model and coordinates updated
! iwearcf = set to 1 to store nodal coefficient of friction for wear calculation
! imixmeth = set=1 then use nonlinear mixture material - allocate memory
! ielcmadyn = flag for magnetodynamics
! 0 - electromagnetics using newmark beta
! 1 - transient magnetics using backward euler
! idinout = flag to control if inside out elements should be deactivated
! igena_meth = 0 - generalized alpha parameters depend on whether or not contact
! is flagged (dynamic,7)
! 10 - generalized alpha parameters are optimized for a contact
! analysis (dynamic,8)
! 11 - generalized alpha parameters are optimized for an analysis
! without contact (dynamic,8)
! magf_meth = - Method to compute force in magnetostatic - structural
! = 1 - Virtual work method based on finite difference for the force computation
! = 2 - Maxwell stress tensor
! = 3 - Virtual work method based on local derivative for the force computation
! non_assumed = 1 no assumed strain formulation (forced)
! iredoboudry set to 1 if contact boundary needs to be recalculated
! ioffsz0 = 1 if composite are used with reference position.ne.0
! icomplt = 1 global flag if a coupled analysis contains an complex pass
! mesh_dual = 1 two independent meshes are used in magnetodynamic/thermal/structural
! one for magnetodynamic and the other for the remaining passes
! iactrp = 1 in an analysis with global remeshing, include inactive
! rigid bodies on post file
! mgnewton = 1 Use full Newton Raphson iteration for magnetostatic pass
!
! iusedens > 0 if mass density is used in the analysis (dynamics, mass dependent loading)
! igsigd0 = 1 set varselem(igsigd) to zero in next oasemb
! iaem = 1 if marc is called from aem (0 - off - default)
! icosim = 1 if marc is used in co-simulation analysis with ADAMS using the CosimEngine
! = 2 if marc is used in co-simulation analysis with ADAMS using the ACSI interface
! = 3 if marc is used in co-simulation analysis with scFLOW using the CosimEngine
! = 4 if marc is used in co-simulation analysis with scFLOW and ADAMS using the CosimEngine
! inodels = 1 nodal integration elements 239/240/241 present
! nlharm = 0 harmonic subincrements are linear
! = 1 harmonic subincrements are nonlinear
! iampini = 0 amplitude of previous harmonic subinc is initial estimate (default)
! = 1 zero amplitude is initial estimate
! iphasetr = 1 phase transformation material model is used
! iforminp flag indicating that contact is switched on via the CONTACT
! option in the input file (as opposed to the case that contact
! is switched on internally due to cyclic symmetry or model
! section creation)
! ispecerror = a+10*b (only for spectrum response analysis with missing mass option)
! a=0 or a=1 (modal shape with non-zero shift)
! b=0 or b=1 (recover with new assembly of stiffness matrix)
! icsprg = set to 1 if spring elements present in model
! imol Control flag for molecualr diffusion pass
! imolt Permanent control flag for molecualr diffusion pass
! Note in coupled analysis imolt will remain as one,
! but imol will be zero in stress pass or thermal pass.
! idatafit = run Marc to fit parameters
! iharmpar = 1 if harmonic parameter option is used
! inclcase load case increment use for cyclic plasticity data fitting
! imultifreq flag to indicate how many harmonic magnetodynamic passes are computed in coupled
! magnetodynamic/thermal(/structural) analyses.
! 0 or 1 one pass 2 two passes 3 or more is not supported
! init_elas use elastic stress-strain law as the material tangent for
! the first cycle of an increment
! ifatig packed integer telling which fatigue mode is active
! 1 = elastomer
! 10 = stress-life
! 100 = strain-life
! = 2 strain-life fatigue
! iftgmat = 0 no fatigue material properties in the dat file
! = 1 fatigue material properties in the dat file
! nchybrid cycle count used for hybrid contact; meant to force an extra iteration
! if the overlap for a node in hybrid contact is too large
! ibuckle buckle parameter option is active
! iexpande set to 1 if expanded elements (248, 249, 250 or 251) are
! present, 0 otherwise
!
!***********************************************************************
!$omp threadprivate(/marc_concom/)
!!

73
src/Marc/include/creeps2022.2 vendored Normal file
View File

@ -0,0 +1,73 @@
! common block definition file taken from respective MSC.Marc release and reformated to free format
!***********************************************************************
!
! File: creeps.cmn
!
! MSC.Marc include file
!
real(pReal) cptim,timinc,timinc_p,timinc_s,timincm,timinc_a,timinc_b
integer icfte,icfst,icfeq,icftm,icetem,mcreep,jcreep,icpa,icftmp,icfstr,&
icfqcp,icfcpm,icrppr,icrcha,icpb,iicpmt,iicpa
real(pReal) time_beg_lcase,time_beg_inc,fractol,time_beg_pst
real(pReal) fraction_donn,timinc_ol2
!
integer num_creepsr,num_creepsi,num_creeps2r,ncrp_arry
parameter(num_creepsr=7)
parameter(num_creepsi=17)
parameter(num_creeps2r=6)
parameter(ncrp_arry=7)
common/marc_creeps/cptim,timinc,timinc_p,timinc_s,timincm,timinc_a,timinc_b,icfte,icfst,&
icfeq,icftm,icetem,mcreep,jcreep,icpa,icftmp,icfstr,icfqcp,icfcpm,icrppr,icrcha,icpb,iicpmt,iicpa
common/marc_creeps2/time_beg_lcase,time_beg_inc,fractol,time_beg_pst,fraction_donn,timinc_ol2
!
! cptim Total time at begining of increment.
! timinc Incremental time for this step.
! icfte Local copy number of slopes of creep strain rate function
! versus temperature. Is -1 if exponent law used.
! icfst Local copy number of slopes of creep strain rate function
! versus equivalent stress. Is -1 if exponent law used.
! icfeq Local copy number of slopes of creep strain rate function
! versus equivalent strain. Is -1 if exponent law used.
! icftm Local copy number of slopes of creep strain rate function
! versus time. Is -1 if exponent law used.
! icetem Element number that needs to be checked for creep convergence
! or, if negative, the number of elements that need to
! be checked. In the latter case the elements to check
! are stored in ielcp.
! mcreep Maximum nuber of iterations for explicit creep.
! jcreep Counter of number of iterations for explicit creep
! procedure. jcreep must be .le. mcreep
! icpa(1-6) Pointer to constants in creep strain rate expression.
! icftmp Pointer to temperature dependent creep strain rate data.
! icfstr Pointer to equivalent stress dependent creep strain rate data.
! icfqcp Pointer to equivalent creep strain dependent creep strain
! rate data.
! icfcpm Pointer to equivalent creep strain rate dependent
! creep strain rate data.
! icrppr Permanent copy of icreep
! icrcha Control flag for creep convergence checking , if set to
! 1 then testing on absolute change in stress and creep
! strain, not relative testing. Input data.
! icpb(1-4) Pointer to storage of material id cross reference numbers.
! iicpmt creep law type ID
! =1 - power law
! =2 - solder
! =3 - steady-creep
! =4 - hyperbolic steady-creep
! iicpa Pointer to table IDs for constants in creep strain rate
! expression
!
!
! time_beg_lcase time at the beginning of the current load case
! time_beg_inc time at the beginning of the current increment
! fractol fraction of loadcase or increment time when we
! consider it to be finished
! time_beg_pst time corresponding to first increment to be
! read in from thermal post file for auto step
!
! timinc_old Time step of the previous increment
!
!***********************************************************************
!!$omp threadprivate(/marc_creeps/)
!!$omp threadprivate(/marc_creeps2/)
!!

69
src/YAML_types.f90
View File

@ -196,6 +196,7 @@ subroutine selfTest
s1 = '1'
s2 = '2'
allocate(l)
if (l%contains('1')) error stop 'empty tList_contains'
call l%append(s1)
call l%append(s2)
if (any(l%as1dInt() /= [1,2])) error stop 'tList_as1dInt'
@ -204,7 +205,9 @@ subroutine selfTest
s2 = 'false'
if (any(l%as1dBool() .neqv. [.true.,.false.])) error stop 'tList_as1dBool'
if (any(l%as1dString() /= ['true ','false'])) error stop 'tList_as1dString'
if (l%asFormattedString() /= '[true, false]') error stop 'tScalar_asFormattedString'
if (l%asFormattedString() /= '[true, false]') error stop 'tList_asFormattedString'
if ( .not. l%contains('true') &
.or. .not. l%contains('false')) error stop 'tList_contains'
end block list
@ -226,6 +229,8 @@ subroutine selfTest
s3 = '3'
s4 = '4'
allocate(d)
if (d%contains('one-two')) error stop 'empty tDict_contains'
if (d%get_asInt('one-two',defaultVal=-1) /= -1) error stop 'empty tDict_get'
call d%set('one-two',l)
call d%set('three',s3)
call d%set('four',s4)
@ -233,6 +238,13 @@ subroutine selfTest
error stop 'tDict_asFormattedString'
if (d%get_asInt('three') /= 3) error stop 'tDict_get_asInt'
if (any(d%get_as1dInt('one-two') /= [1,2])) error stop 'tDict_get_as1dInt'
call d%set('one-two',s4)
if (d%asFormattedString() /= '{one-two: 4, three: 3, four: 4}') &
error stop 'tDict_set overwrite'
if ( .not. d%contains('one-two') &
.or. .not. d%contains('three') &
.or. .not. d%contains('four') &
) error stop 'tDict_contains'
end block dict
@ -407,7 +419,7 @@ recursive function tList_asFormattedString(self) result(str)
str = '['
item => self%first
do i = 1, self%length -1
do i = 2, self%length
str = str//item%node%asFormattedString()//', '
item => item%next
end do
@ -593,15 +605,14 @@ function tList_contains(self,k) result(exists)
type(tScalar), pointer :: scalar
exists = .false.
item => self%first
do j = 1, self%length
exists = .false.
j = 1
do while (j <= self%length .and. .not. exists)
scalar => item%node%asScalar()
if (scalar%value == k) then
exists = .true.
exit
endif
exists = scalar%value == k
item => item%next
j = j + 1
end do
end function tList_contains
@ -620,7 +631,8 @@ function tList_get(self,i) result(node)
integer :: j
if (i < 1 .or. i > self%length) call IO_error(150,ext_msg='tList_get @ '//IO_intAsString(i))
if (i < 1 .or. i > self%length) call IO_error(150,ext_msg='tList_get @ '//IO_intAsString(i) &
//' of '//IO_intAsString(self%length) )
item => self%first
do j = 2, i
item => item%next
@ -854,7 +866,7 @@ recursive function tDict_asFormattedString(self) result(str)
str = '{'
item => self%first
do i = 1, self%length -1
do i = 2, self%length
str = str//trim(item%key)//': '//item%node%asFormattedString()//', '
item => item%next
end do
@ -881,8 +893,7 @@ subroutine tDict_set(self,key,node)
self%length = 1
else
item => self%first
searchExisting: do while (associated(item%next))
if (item%key == key) exit
searchExisting: do while (associated(item%next) .and. item%key /= key)
item => item%next
end do searchExisting
if (item%key /= key) then
@ -935,9 +946,10 @@ function tDict_key(self,i) result(key)
type(tItem), pointer :: item
if (i < 1 .or. i > self%length) call IO_error(150,ext_msg='tDict_key @ '//IO_intAsString(i))
if (i < 1 .or. i > self%length) call IO_error(150,ext_msg='tDict_key @ '//IO_intAsString(i) &
//' of '//IO_intAsString(self%length) )
item => self%first
do j = 1, i-1
do j = 2, i
item => item%next
end do
@ -986,11 +998,10 @@ function tDict_contains(self,k) result(exists)
exists = .false.
do j=1, self%length
if (self%key(j) == k) then
exists = .true.
return
end if
j = 1
do while(j <= self%length .and. .not. exists)
exists = self%key(j) == k
j = j + 1
end do
end function tDict_contains
@ -1008,27 +1019,21 @@ function tDict_get(self,k,defaultVal) result(node)
type(tItem), pointer :: item
integer :: j
logical :: found
found = present(defaultVal)
if (found) node => defaultVal
j = 1
item => self%first
do while(j <= self%length)
do j=1, self%length
if (item%key == k) then
found = .true.
exit
node => item%node
return
end if
item => item%next
j = j + 1
end do
if (.not. found) then
call IO_error(143,ext_msg=k)
if (present(defaultVal)) then
node => defaultVal
else
if (associated(item)) node => item%node
call IO_error(143,ext_msg=k)
end if
end function tDict_get

2
src/grid/DAMASK_grid.f90
View File

@ -106,8 +106,6 @@ program DAMASK_grid
external :: &
quit
class(tNode), pointer :: &
tmp
type(tDict), pointer :: &
config_load, &
num_grid, &

42
src/grid/grid_damage_spectral.f90
View File

@ -33,7 +33,7 @@ module grid_damage_spectral
integer :: &
itmax !< maximum number of iterations
real(pReal) :: &
residualStiffness, & !< non-zero residual damage
phi_min, & !< non-zero residual damage
eps_damage_atol, & !< absolute tolerance for damage evolution
eps_damage_rtol !< relative tolerance for damage evolution
end type tNumerics
@ -95,9 +95,9 @@ subroutine grid_damage_spectral_init()
num%eps_damage_rtol = num_grid%get_asFloat ('eps_damage_rtol',defaultVal=1.0e-6_pReal)
num_generic => config_numerics%get_dict('generic',defaultVal=emptyDict)
num%residualStiffness = num_generic%get_asFloat('residualStiffness', defaultVal=1.0e-6_pReal)
num%phi_min = num_generic%get_asFloat('phi_min', defaultVal=1.0e-6_pReal)
if (num%residualStiffness < 0.0_pReal) call IO_error(301,ext_msg='residualStiffness')
if (num%phi_min < 0.0_pReal) call IO_error(301,ext_msg='phi_min')
if (num%itmax <= 1) call IO_error(301,ext_msg='itmax')
if (num%eps_damage_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_damage_atol')
if (num%eps_damage_rtol <= 0.0_pReal) call IO_error(301,ext_msg='eps_damage_rtol')
@ -253,11 +253,13 @@ end function grid_damage_spectral_solution
subroutine grid_damage_spectral_forward(cutBack)
logical, intent(in) :: cutBack
integer :: i, j, k, ce
DM :: dm_local
PetscScalar, dimension(:,:,:), pointer :: phi_PETSc
PetscErrorCode :: err_PETSc
if (cutBack) then
phi_current = phi_lastInc
phi_stagInc = phi_lastInc
@ -284,7 +286,7 @@ end subroutine grid_damage_spectral_forward
!--------------------------------------------------------------------------------------------------
!> @brief forms the spectral damage residual vector
!> @brief Construct the residual vector.
!--------------------------------------------------------------------------------------------------
subroutine formResidual(in,x_scal,r,dummy,err_PETSc)
@ -297,48 +299,34 @@ subroutine formResidual(in,x_scal,r,dummy,err_PETSc)
X_RANGE,Y_RANGE,Z_RANGE), intent(out) :: &
r
PetscObject :: dummy
PetscErrorCode :: err_PETSc
PetscErrorCode, intent(out) :: err_PETSc
integer :: i, j, k, ce
real(pReal), dimension(3,cells(1),cells(2),cells3) :: vectorField
phi_current = x_scal
!--------------------------------------------------------------------------------------------------
! evaluate polarization field
scalarField_real = 0.0_pReal
scalarField_real(1:cells(1),1:cells(2),1:cells3) = phi_current
call utilities_FFTscalarForward
call utilities_fourierScalarGradient !< calculate gradient of damage field
call utilities_FFTvectorBackward
vectorField = utilities_ScalarGradient(phi_current)
ce = 0
do k = 1, cells3; do j = 1, cells(2); do i = 1,cells(1)
ce = ce + 1
vectorField_real(1:3,i,j,k) = matmul(homogenization_K_phi(ce) - K_ref, vectorField_real(1:3,i,j,k))
vectorField(1:3,i,j,k) = matmul(homogenization_K_phi(ce) - K_ref, vectorField(1:3,i,j,k))
end do; end do; end do
call utilities_FFTvectorForward
call utilities_fourierVectorDivergence !< calculate damage divergence in fourier field
call utilities_FFTscalarBackward
r = utilities_VectorDivergence(vectorField)
ce = 0
do k = 1, cells3; do j = 1, cells(2); do i = 1,cells(1)
ce = ce + 1
scalarField_real(i,j,k) = params%Delta_t*(scalarField_real(i,j,k) + homogenization_f_phi(phi_current(i,j,k),ce)) &
r(i,j,k) = params%Delta_t*(r(i,j,k) + homogenization_f_phi(phi_current(i,j,k),ce)) &
+ homogenization_mu_phi(ce)*(phi_lastInc(i,j,k) - phi_current(i,j,k)) &
+ mu_ref*phi_current(i,j,k)
end do; end do; end do
!--------------------------------------------------------------------------------------------------
! convolution of damage field with green operator
call utilities_FFTscalarForward
call utilities_fourierGreenConvolution(K_ref, mu_ref, params%Delta_t)
call utilities_FFTscalarBackward
where(scalarField_real(1:cells(1),1:cells(2),1:cells3) > phi_lastInc) &
scalarField_real(1:cells(1),1:cells(2),1:cells3) = phi_lastInc
where(scalarField_real(1:cells(1),1:cells(2),1:cells3) < num%residualStiffness) &
scalarField_real(1:cells(1),1:cells(2),1:cells3) = num%residualStiffness
!--------------------------------------------------------------------------------------------------
! constructing residual
r = scalarField_real(1:cells(1),1:cells(2),1:cells3) - phi_current
r = max(min(utilities_GreenConvolution(r, K_ref, mu_ref, params%Delta_t),phi_lastInc),num%phi_min) &
- phi_current
err_PETSc = 0
end subroutine formResidual

32
src/grid/grid_mech_spectral_basic.f90
View File

@ -491,7 +491,7 @@ end subroutine converged
!--------------------------------------------------------------------------------------------------
!> @brief forms the residual vector
!> @brief Construct the residual vector.
!--------------------------------------------------------------------------------------------------
subroutine formResidual(in, F, &
r, dummy, err_PETSc)
@ -501,15 +501,17 @@ subroutine formResidual(in, F, &
intent(in) :: F !< deformation gradient field
PetscScalar, dimension(3,3,X_RANGE,Y_RANGE,Z_RANGE), &
intent(out) :: r !< residuum field
PetscObject :: dummy
PetscErrorCode :: err_PETSc
real(pReal), dimension(3,3) :: &
deltaF_aim
PetscInt :: &
PETScIter, &
nfuncs
PetscObject :: dummy
PetscErrorCode :: err_PETSc
integer(MPI_INTEGER_KIND) :: err_MPI
call SNESGetNumberFunctionEvals(SNES_mechanical,nfuncs,err_PETSc)
CHKERRQ(err_PETSc)
call SNESGetIterationNumber(SNES_mechanical,PETScIter,err_PETSc)
@ -517,8 +519,6 @@ subroutine formResidual(in, F, &
if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1 ! new increment
!--------------------------------------------------------------------------------------------------
! begin of new iteration
newIteration: if (totalIter <= PETScIter) then
totalIter = totalIter + 1
print'(1x,a,3(a,i0))', trim(incInfo), ' @ Iteration ', num%itmin, '≤',totalIter, '≤', num%itmax
@ -529,32 +529,20 @@ subroutine formResidual(in, F, &
flush(IO_STDOUT)
end if newIteration
!--------------------------------------------------------------------------------------------------
! evaluate constitutive response
call utilities_constitutiveResponse(r, & ! residuum gets field of first PK stress (to save memory)
associate (P => r)
call utilities_constitutiveResponse(P, &
P_av,C_volAvg,C_minMaxAvg, &
F,params%Delta_t,params%rotation_BC)
call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
err_div = utilities_divergenceRMS(P)
end associate
!--------------------------------------------------------------------------------------------------
! stress BC handling
deltaF_aim = math_mul3333xx33(S, P_av - P_aim) ! S = 0.0 for no bc
F_aim = F_aim - deltaF_aim
err_BC = maxval(abs(merge(.0_pReal,P_av - P_aim,params%stress_mask)))
!--------------------------------------------------------------------------------------------------
! updated deformation gradient using fix point algorithm of basic scheme
tensorField_real = 0.0_pReal
tensorField_real(1:3,1:3,1:cells(1),1:cells(2),1:cells3) = r ! store fPK field for subsequent FFT forward transform
call utilities_FFTtensorForward ! FFT forward of global "tensorField_real"
err_div = utilities_divergenceRMS() ! divRMS of tensorField_fourier for later use
call utilities_fourierGammaConvolution(params%rotation_BC%rotate(deltaF_aim,active=.true.)) ! convolution of Gamma and tensorField_fourier
call utilities_FFTtensorBackward ! FFT backward of global tensorField_fourier
!--------------------------------------------------------------------------------------------------
! constructing residual
r = tensorField_real(1:3,1:3,1:cells(1),1:cells(2),1:cells3) ! Gamma*P gives correction towards div(P) = 0, so needs to be zero, too
r = utilities_GammaConvolution(r,params%rotation_BC%rotate(deltaF_aim,active=.true.))
end subroutine formResidual

76
src/grid/grid_mech_spectral_polarisation.f90
View File

@ -541,7 +541,7 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dumm
err_div/divTol, ' (',err_div, ' / m, tol = ',divTol,')'
print '(1x,a,f12.2,a,es8.2,a,es9.2,a)', 'error curl = ', &
err_curl/curlTol,' (',err_curl,' -, tol = ',curlTol,')'
print '(1x,a,f12.2,a,es8.2,a,es9.2,a)', 'error stress BC = ', &
print '(1x,a,f12.2,a,es8.2,a,es9.2,a)', 'error mech BC = ', &
err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')'
print'(/,1x,a)', '==========================================================================='
flush(IO_STDOUT)
@ -551,7 +551,7 @@ end subroutine converged
!--------------------------------------------------------------------------------------------------
!> @brief forms the residual vector
!> @brief Construct the residual vector.
!--------------------------------------------------------------------------------------------------
subroutine formResidual(in, FandF_tau, &
r, dummy,err_PETSc)
@ -561,6 +561,9 @@ subroutine formResidual(in, FandF_tau, &
target, intent(in) :: FandF_tau
PetscScalar, dimension(3,3,2,X_RANGE,Y_RANGE,Z_RANGE),&
target, intent(out) :: r !< residuum field
PetscObject :: dummy
PetscErrorCode :: err_PETSc
PetscScalar, pointer, dimension(:,:,:,:,:) :: &
F, &
F_tau, &
@ -569,13 +572,10 @@ subroutine formResidual(in, FandF_tau, &
PetscInt :: &
PETScIter, &
nfuncs
PetscObject :: dummy
PetscErrorCode :: err_PETSc
integer(MPI_INTEGER_KIND) :: err_MPI
integer :: &
i, j, k, e
!---------------------------------------------------------------------------------------------------
F => FandF_tau(1:3,1:3,1,&
XG_RANGE,YG_RANGE,ZG_RANGE)
@ -597,8 +597,6 @@ subroutine formResidual(in, FandF_tau, &
if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1 ! new increment
!--------------------------------------------------------------------------------------------------
! begin of new iteration
newIteration: if (totalIter <= PETScIter) then
totalIter = totalIter + 1
print'(1x,a,3(a,i0))', trim(incInfo), ' @ Iteration ', num%itmin, '≤',totalIter, '≤', num%itmax
@ -609,63 +607,53 @@ subroutine formResidual(in, FandF_tau, &
flush(IO_STDOUT)
end if newIteration
!--------------------------------------------------------------------------------------------------
!
tensorField_real = 0.0_pReal
do k = 1, cells3; do j = 1, cells(2); do i = 1, cells(1)
tensorField_real(1:3,1:3,i,j,k) = &
r_F_tau(1:3,1:3,i,j,k) = &
num%beta*math_mul3333xx33(C_scale,F(1:3,1:3,i,j,k) - math_I3) -&
num%alpha*matmul(F(1:3,1:3,i,j,k), &
math_mul3333xx33(C_scale,F_tau(1:3,1:3,i,j,k) - F(1:3,1:3,i,j,k) - math_I3))
end do; end do; end do
r_F_tau = num%beta*F &
- utilities_GammaConvolution(r_F_tau,params%rotation_BC%rotate(num%beta*F_aim,active=.true.))
!--------------------------------------------------------------------------------------------------
! doing convolution in Fourier space
call utilities_FFTtensorForward
call utilities_fourierGammaConvolution(params%rotation_BC%rotate(num%beta*F_aim,active=.true.))
call utilities_FFTtensorBackward
err_curl = utilities_curlRMS(F)
!--------------------------------------------------------------------------------------------------
! constructing residual
r_F_tau = num%beta*F - tensorField_real(1:3,1:3,1:cells(1),1:cells(2),1:cells3)
!--------------------------------------------------------------------------------------------------
! evaluate constitutive response
call utilities_constitutiveResponse(r_F, & ! "residuum" gets field of first PK stress (to save memory)
#ifdef __GFORTRAN__
call utilities_constitutiveResponse(r_F, &
#else
associate (P => r_F)
call utilities_constitutiveResponse(P, &
#endif
P_av,C_volAvg,C_minMaxAvg, &
F - r_F_tau/num%beta,params%Delta_t,params%rotation_BC)
call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI)
!--------------------------------------------------------------------------------------------------
! stress BC handling
F_aim = F_aim - math_mul3333xx33(S, P_av - P_aim) ! S = 0.0 for no bc
err_BC = maxval(abs(merge(math_mul3333xx33(C_scale,F_aim-params%rotation_BC%rotate(F_av)), &
P_av-P_aim, &
params%stress_mask)))
! calculate divergence
tensorField_real = 0.0_pReal
tensorField_real(1:3,1:3,1:cells(1),1:cells(2),1:cells3) = r_F !< stress field in disguise
call utilities_FFTtensorForward
err_div = utilities_divergenceRMS() !< root mean squared error in divergence of stress
!--------------------------------------------------------------------------------------------------
! constructing residual
#ifdef __GFORTRAN__
err_div = utilities_divergenceRMS(r_F)
#else
err_div = utilities_divergenceRMS(P)
#endif
e = 0
do k = 1, cells3; do j = 1, cells(2); do i = 1, cells(1)
e = e + 1
r_F(1:3,1:3,i,j,k) = &
math_mul3333xx33(math_invSym3333(homogenization_dPdF(1:3,1:3,1:3,1:3,e) + C_scale), &
#ifdef __GFORTRAN__
r_F(1:3,1:3,i,j,k) - matmul(F(1:3,1:3,i,j,k), &
#else
P(1:3,1:3,i,j,k) - matmul(F(1:3,1:3,i,j,k), &
#endif
math_mul3333xx33(C_scale,F_tau(1:3,1:3,i,j,k) - F(1:3,1:3,i,j,k) - math_I3))) &
+ r_F_tau(1:3,1:3,i,j,k)
end do; end do; end do
#ifndef __GFORTRAN__
end associate
#endif
F_aim = F_aim - math_mul3333xx33(S, P_av - P_aim) ! S = 0.0 for no bc
err_BC = maxval(abs(merge(math_mul3333xx33(C_scale,F_aim-params%rotation_BC%rotate(F_av)), &
P_av-P_aim, &
params%stress_mask)))
!--------------------------------------------------------------------------------------------------
! calculating curl
tensorField_real = 0.0_pReal
tensorField_real(1:3,1:3,1:cells(1),1:cells(2),1:cells3) = F
call utilities_FFTtensorForward
err_curl = utilities_curlRMS()
end subroutine formResidual

32
src/grid/grid_thermal_spectral.f90
View File

@ -242,11 +242,13 @@ end function grid_thermal_spectral_solution
subroutine grid_thermal_spectral_forward(cutBack)
logical, intent(in) :: cutBack
integer :: i, j, k, ce
DM :: dm_local
PetscScalar, dimension(:,:,:), pointer :: T_PETSc
PetscErrorCode :: err_PETSc
if (cutBack) then
T_current = T_lastInc
T_stagInc = T_lastInc
@ -307,7 +309,7 @@ end subroutine grid_thermal_spectral_restartWrite
!--------------------------------------------------------------------------------------------------
!> @brief forms the spectral thermal residual vector
!> @brief Construct the residual vector.
!--------------------------------------------------------------------------------------------------
subroutine formResidual(in,x_scal,r,dummy,err_PETSc)
@ -320,42 +322,34 @@ subroutine formResidual(in,x_scal,r,dummy,err_PETSc)
X_RANGE,Y_RANGE,Z_RANGE), intent(out) :: &
r
PetscObject :: dummy
PetscErrorCode :: err_PETSc
PetscErrorCode, intent(out) :: err_PETSc
integer :: i, j, k, ce
real(pReal), dimension(3,cells(1),cells(2),cells3) :: vectorField
T_current = x_scal
!--------------------------------------------------------------------------------------------------
! evaluate polarization field
scalarField_real = 0.0_pReal
scalarField_real(1:cells(1),1:cells(2),1:cells3) = T_current
call utilities_FFTscalarForward
call utilities_fourierScalarGradient !< calculate gradient of temperature field
call utilities_FFTvectorBackward
vectorField = utilities_ScalarGradient(T_current)
ce = 0
do k = 1, cells3; do j = 1, cells(2); do i = 1,cells(1)
ce = ce + 1
vectorField_real(1:3,i,j,k) = matmul(homogenization_K_T(ce) - K_ref, vectorField_real(1:3,i,j,k))
vectorField(1:3,i,j,k) = matmul(homogenization_K_T(ce) - K_ref, vectorField(1:3,i,j,k))
end do; end do; end do
call utilities_FFTvectorForward
call utilities_fourierVectorDivergence !< calculate temperature divergence in fourier field
call utilities_FFTscalarBackward
r = utilities_VectorDivergence(vectorField)
ce = 0
do k = 1, cells3; do j = 1, cells(2); do i = 1,cells(1)
ce = ce + 1
scalarField_real(i,j,k) = params%Delta_t*(scalarField_real(i,j,k) + homogenization_f_T(ce)) &
r(i,j,k) = params%Delta_t*(r(i,j,k) + homogenization_f_T(ce)) &
+ homogenization_mu_T(ce) * (T_lastInc(i,j,k) - T_current(i,j,k)) &
+ mu_ref*T_current(i,j,k)
end do; end do; end do
!--------------------------------------------------------------------------------------------------
! convolution of temperature field with green operator
call utilities_FFTscalarForward
call utilities_fourierGreenConvolution(K_ref, mu_ref, params%Delta_t)
call utilities_FFTscalarBackward
!--------------------------------------------------------------------------------------------------
! constructing residual
r = T_current - scalarField_real(1:cells(1),1:cells(2),1:cells3)
r = T_current &
- utilities_GreenConvolution(r, K_ref, mu_ref, params%Delta_t)
err_PETSc = 0
end subroutine formResidual

220
src/grid/spectral_utilities.f90
View File

@ -42,9 +42,9 @@ module spectral_utilities
!--------------------------------------------------------------------------------------------------
! variables storing information for spectral method and FFTW
real(C_DOUBLE), public, dimension(:,:,:,:,:), pointer :: tensorField_real !< tensor field in real space
real(C_DOUBLE), public, dimension(:,:,:,:), pointer :: vectorField_real !< vector field in real space
real(C_DOUBLE), public, dimension(:,:,:), pointer :: scalarField_real !< scalar field in real space
real(C_DOUBLE), dimension(:,:,:,:,:), pointer :: tensorField_real !< tensor field in real space
real(C_DOUBLE), dimension(:,:,:,:), pointer :: vectorField_real !< vector field in real space
real(C_DOUBLE), dimension(:,:,:), pointer :: scalarField_real !< scalar field in real space
complex(C_DOUBLE_COMPLEX), dimension(:,:,:,:,:), pointer :: tensorField_fourier !< tensor field in Fourier space
complex(C_DOUBLE_COMPLEX), dimension(:,:,:,:), pointer :: vectorField_fourier !< vector field in Fourier space
complex(C_DOUBLE_COMPLEX), dimension(:,:,:), pointer :: scalarField_fourier !< scalar field in Fourier space
@ -116,18 +116,12 @@ module spectral_utilities
public :: &
spectral_utilities_init, &
utilities_updateGamma, &
utilities_FFTtensorForward, &
utilities_FFTtensorBackward, &
utilities_FFTvectorForward, &
utilities_FFTvectorBackward, &
utilities_FFTscalarForward, &
utilities_FFTscalarBackward, &
utilities_fourierGammaConvolution, &
utilities_fourierGreenConvolution, &
utilities_GammaConvolution, &
utilities_GreenConvolution, &
utilities_divergenceRMS, &
utilities_curlRMS, &
utilities_fourierScalarGradient, &
utilities_fourierVectorDivergence, &
utilities_ScalarGradient, &
utilities_VectorDivergence, &
utilities_maskedCompliance, &
utilities_constitutiveResponse, &
utilities_calculateRate, &
@ -385,6 +379,7 @@ end subroutine spectral_utilities_init
subroutine utilities_updateGamma(C)
real(pReal), intent(in), dimension(3,3,3,3) :: C !< input stiffness to store as reference stiffness
complex(pReal), dimension(3,3) :: temp33_cmplx, xiDyad_cmplx
real(pReal), dimension(6,6) :: A, A_inv
integer :: &
@ -392,7 +387,8 @@ subroutine utilities_updateGamma(C)
l, m, n, o
logical :: err
C_ref = C
C_ref = C/wgt
if (.not. num%memory_efficient) then
gamma_hat = cmplx(0.0_pReal,0.0_pReal,pReal) ! for the singular point and any non invertible A
@ -434,68 +430,6 @@ subroutine utilities_updateGamma(C)
end subroutine utilities_updateGamma
!--------------------------------------------------------------------------------------------------
!> @brief forward FFT of data in field_real to field_fourier
!> @details Does an unweighted FFT transform from real to complex. Extra padding entries are set
! to 0.0
!--------------------------------------------------------------------------------------------------
subroutine utilities_FFTtensorForward()
tensorField_real(1:3,1:3,cells(1)+1:cells1Red*2,:,:) = 0.0_pReal
call fftw_mpi_execute_dft_r2c(planTensorForth,tensorField_real,tensorField_fourier)
end subroutine utilities_FFTtensorForward
!--------------------------------------------------------------------------------------------------
!> @brief backward FFT of data in field_fourier to field_real
!> @details Does an weighted inverse FFT transform from complex to real
!--------------------------------------------------------------------------------------------------
subroutine utilities_FFTtensorBackward()
call fftw_mpi_execute_dft_c2r(planTensorBack,tensorField_fourier,tensorField_real)
tensorField_real = tensorField_real * wgt ! normalize the result by number of elements
end subroutine utilities_FFTtensorBackward
!--------------------------------------------------------------------------------------------------
!> @brief forward FFT of data in scalarField_real to scalarField_fourier
!> @details Does an unweighted FFT transform from real to complex. Extra padding entries are set
! to 0.0
!--------------------------------------------------------------------------------------------------
subroutine utilities_FFTscalarForward()
scalarField_real(cells(1)+1:cells1Red*2,:,:) = 0.0_pReal
call fftw_mpi_execute_dft_r2c(planScalarForth,scalarField_real,scalarField_fourier)
end subroutine utilities_FFTscalarForward
!--------------------------------------------------------------------------------------------------
!> @brief backward FFT of data in scalarField_fourier to scalarField_real
!> @details Does an weighted inverse FFT transform from complex to real
!--------------------------------------------------------------------------------------------------
subroutine utilities_FFTscalarBackward()
call fftw_mpi_execute_dft_c2r(planScalarBack,scalarField_fourier,scalarField_real)
scalarField_real = scalarField_real * wgt ! normalize the result by number of elements
end subroutine utilities_FFTscalarBackward
!--------------------------------------------------------------------------------------------------
!> @brief forward FFT of data in field_real to field_fourier with highest freqs. removed
!> @details Does an unweighted FFT transform from real to complex. Extra padding entries are set
! to 0.0
!--------------------------------------------------------------------------------------------------
subroutine utilities_FFTvectorForward()
vectorField_real(1:3,cells(1)+1:cells1Red*2,:,:) = 0.0_pReal
call fftw_mpi_execute_dft_r2c(planVectorForth,vectorField_real,vectorField_fourier)
end subroutine utilities_FFTvectorForward
!--------------------------------------------------------------------------------------------------
!> @brief backward FFT of data in field_fourier to field_real
!> @details Does an weighted inverse FFT transform from complex to real
@ -511,12 +445,14 @@ end subroutine utilities_FFTvectorBackward
!--------------------------------------------------------------------------------------------------
!> @brief doing convolution gamma_hat * field_real, ensuring that average value = fieldAim
!--------------------------------------------------------------------------------------------------
subroutine utilities_fourierGammaConvolution(fieldAim)
function utilities_GammaConvolution(field, fieldAim) result(gammaField)
real(pReal), intent(in), dimension(3,3,cells(1),cells(2),cells3) :: field
real(pReal), intent(in), dimension(3,3) :: fieldAim !< desired average value of the field after convolution
real(pReal), dimension(3,3,cells(1),cells(2),cells3) :: gammaField
complex(pReal), dimension(3,3) :: temp33_cmplx, xiDyad_cmplx
real(pReal), dimension(6,6) :: A, A_inv
integer :: &
i, j, k, &
l, m, n, o
@ -526,8 +462,10 @@ subroutine utilities_fourierGammaConvolution(fieldAim)
print'(/,1x,a)', '... doing gamma convolution ...............................................'
flush(IO_STDOUT)
!--------------------------------------------------------------------------------------------------
! do the actual spectral method calculation (mechanical equilibrium)
tensorField_real(1:3,1:3,cells(1)+1:cells1Red*2,1:cells(2),1:cells3) = 0.0_pReal
tensorField_real(1:3,1:3,1:cells(1), 1:cells(2),1:cells3) = field
call fftw_mpi_execute_dft_r2c(planTensorForth,tensorField_real,tensorField_fourier)
memoryEfficient: if (num%memory_efficient) then
!$OMP PARALLEL DO PRIVATE(l,m,n,o,temp33_cmplx,xiDyad_cmplx,A,A_inv,err,gamma_hat)
do j = 1, cells2; do k = 1, cells(3); do i = 1, cells1Red
@ -586,39 +524,53 @@ subroutine utilities_fourierGammaConvolution(fieldAim)
!$OMP END PARALLEL DO
end if memoryEfficient
if (cells3Offset == 0) tensorField_fourier(1:3,1:3,1,1,1) = cmplx(fieldAim/wgt,0.0_pReal,pReal)
if (cells3Offset == 0) tensorField_fourier(1:3,1:3,1,1,1) = cmplx(fieldAim,0.0_pReal,pReal)
end subroutine utilities_fourierGammaConvolution
call fftw_mpi_execute_dft_c2r(planTensorBack,tensorField_fourier,tensorField_real)
gammaField = tensorField_real(1:3,1:3,1:cells(1),1:cells(2),1:cells3)
end function utilities_GammaConvolution
!--------------------------------------------------------------------------------------------------
!> @brief doing convolution DamageGreenOp_hat * field_real
!> @brief Convolution of Greens' operator for damage/thermal.
!--------------------------------------------------------------------------------------------------
subroutine utilities_fourierGreenConvolution(D_ref, mu_ref, Delta_t)
function utilities_GreenConvolution(field, D_ref, mu_ref, Delta_t) result(greenField)
real(pReal), intent(in), dimension(cells(1),cells(2),cells3) :: field
real(pReal), dimension(3,3), intent(in) :: D_ref
real(pReal), intent(in) :: mu_ref, Delta_t
real(pReal), dimension(cells(1),cells(2),cells3) :: greenField
complex(pReal) :: GreenOp_hat
integer :: i, j, k
!--------------------------------------------------------------------------------------------------
! do the actual spectral method calculation
scalarField_real(cells(1)+1:cells1Red*2,1:cells(2),1:cells3) = 0.0_pReal
scalarField_real(1:cells(1), 1:cells(2),1:cells3) = field
call fftw_mpi_execute_dft_r2c(planScalarForth,scalarField_real,scalarField_fourier)
!$OMP PARALLEL DO PRIVATE(GreenOp_hat)
do j = 1, cells2; do k = 1, cells(3); do i = 1, cells1Red
GreenOp_hat = cmplx(1.0_pReal,0.0_pReal,pReal) &
GreenOp_hat = cmplx(wgt,0.0_pReal,pReal) &
/ (cmplx(mu_ref,0.0_pReal,pReal) + cmplx(Delta_t,0.0_pReal,pReal) &
* sum(conjg(xi1st(1:3,i,k,j))* matmul(cmplx(D_ref,0.0_pReal,pReal),xi1st(1:3,i,k,j))))
scalarField_fourier(i,k,j) = scalarField_fourier(i,k,j)*GreenOp_hat
end do; end do; end do
!$OMP END PARALLEL DO
end subroutine utilities_fourierGreenConvolution
call fftw_mpi_execute_dft_c2r(planScalarBack,scalarField_fourier,scalarField_real)
greenField = scalarField_real(1:cells(1),1:cells(2),1:cells3)
end function utilities_GreenConvolution
!--------------------------------------------------------------------------------------------------
!> @brief calculate root mean square of divergence of field_fourier
!> @brief Calculate root mean square of divergence.
!--------------------------------------------------------------------------------------------------
real(pReal) function utilities_divergenceRMS()
real(pReal) function utilities_divergenceRMS(tensorField)
real(pReal), dimension(3,3,cells(1),cells(2),cells3), intent(in) :: tensorField
integer :: i, j, k
integer(MPI_INTEGER_KIND) :: err_MPI
@ -628,6 +580,10 @@ real(pReal) function utilities_divergenceRMS()
print'(/,1x,a)', '... calculating divergence ................................................'
flush(IO_STDOUT)
tensorField_real(1:3,1:3,cells(1)+1:cells1Red*2,1:cells(2),1:cells3) = 0.0_pReal
tensorField_real(1:3,1:3,1:cells(1), 1:cells(2),1:cells3) = tensorField
call fftw_mpi_execute_dft_r2c(planTensorForth,tensorField_real,tensorField_fourier)
rescaledGeom = cmplx(geomSize/scaledGeomSize,0.0_pReal,pReal)
!--------------------------------------------------------------------------------------------------
@ -660,9 +616,11 @@ end function utilities_divergenceRMS
!--------------------------------------------------------------------------------------------------
!> @brief calculate max of curl of field_fourier
!> @brief Calculate root mean square of curl.
!--------------------------------------------------------------------------------------------------
real(pReal) function utilities_curlRMS()
real(pReal) function utilities_curlRMS(tensorField)
real(pReal), dimension(3,3,cells(1),cells(2),cells3), intent(in) :: tensorField
integer :: i, j, k, l
integer(MPI_INTEGER_KIND) :: err_MPI
@ -673,6 +631,10 @@ real(pReal) function utilities_curlRMS()
print'(/,1x,a)', '... calculating curl ......................................................'
flush(IO_STDOUT)
tensorField_real(1:3,1:3,cells(1)+1:cells1Red*2,1:cells(2),1:cells3) = 0.0_pReal
tensorField_real(1:3,1:3,1:cells(1), 1:cells(2),1:cells3) = tensorField
call fftw_mpi_execute_dft_r2c(planTensorForth,tensorField_real,tensorField_fourier)
rescaledGeom = cmplx(geomSize/scaledGeomSize,0.0_pReal,pReal)
!--------------------------------------------------------------------------------------------------
@ -723,7 +685,7 @@ end function utilities_curlRMS
!--------------------------------------------------------------------------------------------------
!> @brief calculates mask compliance tensor used to adjust F to fullfill stress BC
!> @brief Calculate masked compliance tensor used to adjust F to fullfill stress BC.
!--------------------------------------------------------------------------------------------------
function utilities_maskedCompliance(rot_BC,mask_stress,C)
@ -793,29 +755,46 @@ end function utilities_maskedCompliance
!--------------------------------------------------------------------------------------------------
!> @brief calculate scalar gradient in fourier field
!> @brief Calculate gradient of scalar field.
!--------------------------------------------------------------------------------------------------
subroutine utilities_fourierScalarGradient()
function utilities_ScalarGradient(field) result(grad)
real(pReal), intent(in), dimension( cells(1),cells(2),cells3) :: field
real(pReal), dimension(3,cells(1),cells(2),cells3) :: grad
integer :: i, j, k
scalarField_real(cells(1)+1:cells1Red*2,1:cells(2),1:cells3) = 0.0_pReal
scalarField_real(1:cells(1), 1:cells(2),1:cells3) = field
call fftw_mpi_execute_dft_r2c(planScalarForth,scalarField_real,scalarField_fourier)
do j = 1, cells2; do k = 1, cells(3); do i = 1,cells1Red
vectorField_fourier(1:3,i,k,j) = scalarField_fourier(i,k,j)*xi1st(1:3,i,k,j) ! ToDo: no -conjg?
end do; end do; end do
call fftw_mpi_execute_dft_c2r(planVectorBack,vectorField_fourier,vectorField_real)
grad = vectorField_real(1:3,1:cells(1),1:cells(2),1:cells3)*wgt
end subroutine utilities_fourierScalarGradient
end function utilities_ScalarGradient
!--------------------------------------------------------------------------------------------------
!> @brief calculate vector divergence in fourier field
!> @brief Calculate divergence of vector field.
!--------------------------------------------------------------------------------------------------
subroutine utilities_fourierVectorDivergence()
function utilities_VectorDivergence(field) result(div)
real(pReal), intent(in), dimension(3,cells(1),cells(2),cells3) :: field
real(pReal), dimension( cells(1),cells(2),cells3) :: div
vectorField_real(1:3,cells(1)+1:cells1Red*2,1:cells(2),1:cells3) = 0.0_pReal
vectorField_real(1:3,1:cells(1), 1:cells(2),1:cells3) = field
call fftw_mpi_execute_dft_r2c(planVectorForth,vectorField_real,vectorField_fourier)
scalarField_fourier(1:cells1Red,1:cells(3),1:cells2) = sum(vectorField_fourier(1:3,1:cells1Red,1:cells(3),1:cells2) &
*conjg(-xi1st),1)
call fftw_mpi_execute_dft_c2r(planScalarBack,scalarField_fourier,scalarField_real)
div = scalarField_real(1:cells(1),1:cells(2),1:cells3)*wgt
end subroutine utilities_fourierVectorDivergence
end function utilities_VectorDivergence
!--------------------------------------------------------------------------------------------------
@ -1046,11 +1025,19 @@ subroutine utilities_updateCoords(F)
step = geomSize/real(cells, pReal)
tensorField_real(1:3,1:3,1:cells(1), 1:cells(2),1:cells3) = F
tensorField_real(1:3,1:3,cells(1)+1:cells1Red*2,1:cells(2),1:cells3) = 0.0_pReal
call fftw_mpi_execute_dft_r2c(planTensorForth,tensorField_real,tensorField_fourier)
!--------------------------------------------------------------------------------------------------
! average F
if (cells3Offset == 0) Favg = tensorField_fourier(1:3,1:3,1,1,1)%re*wgt
call MPI_Bcast(Favg,9_MPI_INTEGER_KIND,MPI_DOUBLE,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
!--------------------------------------------------------------------------------------------------
! integration in Fourier space to get fluctuations of cell center discplacements
tensorField_real(1:3,1:3,1:cells(1),1:cells(2),1:cells3) = F
call utilities_FFTtensorForward()
!$OMP PARALLEL DO
do j = 1, cells2; do k = 1, cells(3); do i = 1, cells1Red
if (any([i,j+cells2Offset,k] /= 1)) then
@ -1062,13 +1049,8 @@ subroutine utilities_updateCoords(F)
end do; end do; end do
!$OMP END PARALLEL DO
call utilities_FFTvectorBackward()
!--------------------------------------------------------------------------------------------------
! average F
if (cells3Offset == 0) Favg = tensorField_fourier(1:3,1:3,1,1,1)%re*wgt
call MPI_Bcast(Favg,9_MPI_INTEGER_KIND,MPI_DOUBLE,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
call fftw_mpi_execute_dft_c2r(planVectorBack,vectorField_fourier,vectorField_real)
vectorField_real = vectorField_real * wgt ! normalize the result by number of elements
!--------------------------------------------------------------------------------------------------
! pad cell center fluctuations along z-direction (needed when running MPI simulation)
@ -1136,7 +1118,7 @@ subroutine utilities_saveReferenceStiffness()
open(newunit=fileUnit, file=getSolverJobName()//'.C_ref',&
status='replace',access='stream',action='write',iostat=ierr)
if (ierr /=0) call IO_error(100,ext_msg='could not open file '//getSolverJobName()//'.C_ref')
write(fileUnit) C_ref
write(fileUnit) C_ref*wgt
close(fileUnit)
end if
@ -1156,38 +1138,38 @@ subroutine selfTest()
call random_number(tensorField_real)
tensorField_real(1:3,1:3,cells(1)+1:cells1Red*2,:,:) = 0.0_pReal
tensorField_real_ = tensorField_real
call utilities_FFTtensorForward()
call fftw_mpi_execute_dft_r2c(planTensorForth,tensorField_real,tensorField_fourier)
if (worldsize==1) then
if (any(dNeq(sum(sum(sum(tensorField_real_,dim=5),dim=4),dim=3)/tensorField_fourier(:,:,1,1,1)%re,1.0_pReal,1.0e-12_pReal))) &
error stop 'tensorField avg'
endif
call utilities_FFTtensorBackward()
call fftw_mpi_execute_dft_c2r(planTensorBack,tensorField_fourier,tensorField_real)
tensorField_real(1:3,1:3,cells(1)+1:cells1Red*2,:,:) = 0.0_pReal
if (maxval(abs(tensorField_real_ - tensorField_real))>5.0e-15_pReal) error stop 'tensorField'
if (maxval(abs(tensorField_real_ - tensorField_real*wgt))>5.0e-15_pReal) error stop 'tensorField'
call random_number(vectorField_real)
vectorField_real(1:3,cells(1)+1:cells1Red*2,:,:) = 0.0_pReal
vectorField_real_ = vectorField_real
call utilities_FFTvectorForward()
call fftw_mpi_execute_dft_r2c(planVectorForth,vectorField_real,vectorField_fourier)
if (worldsize==1) then
if (any(dNeq(sum(sum(sum(vectorField_real_,dim=4),dim=3),dim=2)/vectorField_fourier(:,1,1,1)%re,1.0_pReal,1.0e-12_pReal))) &
error stop 'vector avg'
endif
call utilities_FFTvectorBackward()
call fftw_mpi_execute_dft_c2r(planVectorBack,vectorField_fourier,vectorField_real)
vectorField_real(1:3,cells(1)+1:cells1Red*2,:,:) = 0.0_pReal
if (maxval(abs(vectorField_real_ - vectorField_real))>5.0e-15_pReal) error stop 'vectorField'
if (maxval(abs(vectorField_real_ - vectorField_real*wgt))>5.0e-15_pReal) error stop 'vectorField'
call random_number(scalarField_real)
scalarField_real(cells(1)+1:cells1Red*2,:,:) = 0.0_pReal
scalarField_real_ = scalarField_real
call utilities_FFTscalarForward()
call fftw_mpi_execute_dft_r2c(planScalarForth,scalarField_real,scalarField_fourier)
if (worldsize==1) then
if (dNeq(sum(sum(sum(scalarField_real_,dim=3),dim=2),dim=1)/scalarField_fourier(1,1,1)%re,1.0_pReal,1.0e-12_pReal)) &
error stop 'scalar avg'
endif
call utilities_FFTscalarBackward()
call fftw_mpi_execute_dft_c2r(planScalarBack,scalarField_fourier,scalarField_real)
scalarField_real(cells(1)+1:cells1Red*2,:,:) = 0.0_pReal
if (maxval(abs(scalarField_real_ - scalarField_real))>5.0e-15_pReal) error stop 'scalarField'
if (maxval(abs(scalarField_real_ - scalarField_real*wgt))>5.0e-15_pReal) error stop 'scalarField'
end subroutine selfTest

20
src/phase_mechanical_eigen_thermalexpansion.f90
View File

@ -9,8 +9,8 @@ submodule(phase:eigen) thermalexpansion
type :: tParameters
type(tPolynomial) :: &
A_11, &
A_33
Alpha_11, &
Alpha_33
end type tParameters
type(tParameters), dimension(:), allocatable :: param
@ -57,9 +57,9 @@ module function thermalexpansion_init(kinematics_length) result(myKinematics)
if (myKinematics(k,p)) then
associate(prm => param(kinematics_thermal_expansion_instance(p)))
prm%A_11 = polynomial(kinematics%get_dict(k),'A_11','T')
prm%Alpha_11 = polynomial(kinematics%get_dict(k),'Alpha_11','T')
if (any(phase_lattice(p) == ['hP','tI'])) &
prm%A_33 = polynomial(kinematics%get_dict(k),'A_33','T')
prm%Alpha_33 = polynomial(kinematics%get_dict(k),'Alpha_33','T')
end associate
end if
end do
@ -80,7 +80,7 @@ module subroutine thermalexpansion_LiAndItsTangent(Li, dLi_dTstar, ph,me)
dLi_dTstar !< derivative of Li with respect to Tstar (4th-order tensor defined to be zero)
real(pReal) :: T, dot_T
real(pReal), dimension(3,3) :: A
real(pReal), dimension(3,3) :: Alpha
T = thermal_T(ph,me)
@ -88,11 +88,11 @@ module subroutine thermalexpansion_LiAndItsTangent(Li, dLi_dTstar, ph,me)
associate(prm => param(kinematics_thermal_expansion_instance(ph)))
A = 0.0_pReal
A(1,1) = prm%A_11%at(T)
if (any(phase_lattice(ph) == ['hP','tI'])) A(3,3) = prm%A_33%at(T)
A = lattice_symmetrize_33(A,phase_lattice(ph))
Li = dot_T * A
Alpha = 0.0_pReal
Alpha(1,1) = prm%Alpha_11%at(T)
if (any(phase_lattice(ph) == ['hP','tI'])) Alpha(3,3) = prm%Alpha_33%at(T)
Alpha = lattice_symmetrize_33(Alpha,phase_lattice(ph))
Li = dot_T * Alpha
end associate
dLi_dTstar = 0.0_pReal

12
src/phase_mechanical_elastic.f90
View File

@ -45,17 +45,17 @@ module subroutine elastic_init(phases)
associate(prm => param(ph))
prm%C_11 = polynomial(elastic%asDict(),'C_11','T')
prm%C_12 = polynomial(elastic%asDict(),'C_12','T')
prm%C_44 = polynomial(elastic%asDict(),'C_44','T')
prm%C_11 = polynomial(elastic,'C_11','T')
prm%C_12 = polynomial(elastic,'C_12','T')
prm%C_44 = polynomial(elastic,'C_44','T')
if (any(phase_lattice(ph) == ['hP','tI'])) then
prm%C_13 = polynomial(elastic%asDict(),'C_13','T')
prm%C_33 = polynomial(elastic%asDict(),'C_33','T')
prm%C_13 = polynomial(elastic,'C_13','T')
prm%C_33 = polynomial(elastic,'C_33','T')
end if
if (phase_lattice(ph) == 'tI') &
prm%C_66 = polynomial(elastic%asDict(),'C_66','T')
prm%C_66 = polynomial(elastic,'C_66','T')
end associate
end do

4
src/phase_mechanical_plastic_dislotwin.f90
View File

@ -303,7 +303,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
prm%b_tr = math_expand(prm%b_tr,prm%N_tr)
prm%i_tr = pl%get_asFloat('i_tr')
prm%Delta_G = polynomial(pl%asDict(),'Delta_G','T')
prm%Delta_G = polynomial(pl,'Delta_G','T')
prm%L_tr = pl%get_asFloat('L_tr')
a_cF = prm%b_tr(1)*sqrt(6.0_pReal) ! b_tr is Shockley partial
prm%h = 5.0_pReal * a_cF/sqrt(3.0_pReal)
@ -361,7 +361,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
end if
if (prm%sum_N_tw + prm%sum_N_tr > 0 .or. prm%extendedDislocations) &
prm%Gamma_sf = polynomial(pl%asDict(),'Gamma_sf','T')
prm%Gamma_sf = polynomial(pl,'Gamma_sf','T')
slipAndTwinActive: if (prm%sum_N_sl * prm%sum_N_tw > 0) then
prm%h_sl_tw = lattice_interaction_SlipByTwin(N_sl,prm%N_tw,pl%get_as1dFloat('h_sl-tw'), &

10
src/polynomials.f90
View File

@ -157,7 +157,7 @@ subroutine selfTest()
'C,T^4: '//trim(adjustl(coef_s(5)))//IO_EOL//&
'T_ref: '//trim(adjustl(x_ref_s))//IO_EOL
dict => YAML_parse_str_asDict(trim(YAML_s))
p2 = polynomial(dict%asDict(),'C','T')
p2 = polynomial(dict,'C','T')
if (dNeq(p1%at(x),p2%at(x),1.0e-6_pReal)) error stop 'polynomials: init'
y = coef(1)+coef(2)*(x-x_ref)+coef(3)*(x-x_ref)**2+coef(4)*(x-x_ref)**3+coef(5)*(x-x_ref)**4
if (dNeq(p1%at(x),y,1.0e-6_pReal)) error stop 'polynomials: eval(full)'
@ -166,28 +166,28 @@ subroutine selfTest()
'C,T: '//trim(adjustl(coef_s(2)))//IO_EOL//&
'T_ref: '//trim(adjustl(x_ref_s))//IO_EOL
dict => YAML_parse_str_asDict(trim(YAML_s))
p1 = polynomial(dict%asDict(),'C','T')
p1 = polynomial(dict,'C','T')
if (dNeq(p1%at(x_ref+x),-p1%at(x_ref-x),1.0e-10_pReal)) error stop 'polynomials: eval(linear)'
YAML_s = 'C: 0.0'//IO_EOL//&
'C,T^2: '//trim(adjustl(coef_s(3)))//IO_EOL//&
'T_ref: '//trim(adjustl(x_ref_s))//IO_EOL
dict => YAML_parse_str_asDict(trim(YAML_s))
p1 = polynomial(dict%asDict(),'C','T')
p1 = polynomial(dict,'C','T')
if (dNeq(p1%at(x_ref+x),p1%at(x_ref-x),1e-10_pReal)) error stop 'polynomials: eval(quadratic)'
YAML_s = 'Y: '//trim(adjustl(coef_s(1)))//IO_EOL//&
'Y,X^3: '//trim(adjustl(coef_s(2)))//IO_EOL//&
'X_ref: '//trim(adjustl(x_ref_s))//IO_EOL
dict => YAML_parse_str_asDict(trim(YAML_s))
p1 = polynomial(dict%asDict(),'Y','X')
p1 = polynomial(dict,'Y','X')
if (dNeq(p1%at(x_ref+x)-coef(1),-(p1%at(x_ref-x)-coef(1)),1.0e-8_pReal)) error stop 'polynomials: eval(cubic)'
YAML_s = 'Y: '//trim(adjustl(coef_s(1)))//IO_EOL//&
'Y,X^4: '//trim(adjustl(coef_s(2)))//IO_EOL//&
'X_ref: '//trim(adjustl(x_ref_s))//IO_EOL
dict => YAML_parse_str_asDict(trim(YAML_s))
p1 = polynomial(dict%asDict(),'Y','X')
p1 = polynomial(dict,'Y','X')
if (dNeq(p1%at(x_ref+x),p1%at(x_ref-x),1.0e-6_pReal)) error stop 'polynomials: eval(quartic)'