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 petsc-64bit-integer

This commit is contained in:
Martin Diehl 2022-01-12 16:37:22 +01:00
parent 56c22a91a5 731222d099
commit 4bfc814a53
35 changed files with 431 additions and 290 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
.gitlab-ci.yml
View File

@ -45,7 +45,7 @@ variables:
MPI_INTEL: "MPI/Intel/2022.0.1/IntelMPI/2021.5.0"
# ++++++++++++ PETSc ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
PETSC_GNU: "Libraries/PETSc/3.16.1/GNU-10-OpenMPI-4.1.1"
PETSC_INTELLLVM: "Libraries/PETSc/3.16.2/oneAPI-2022.0.1-IntelMPI-2021.5.0"
PETSC_INTELLLVM: "Libraries/PETSc/3.16.3/oneAPI-2022.0.1-IntelMPI-2021.5.0"
PETSC_INTEL: "Libraries/PETSc/3.16.2/Intel-2022.0.1-IntelMPI-2021.5.0"
# ++++++++++++ MSC Marc +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
MSC: "FEM/MSC/2021.3.1"

12
CMakeLists.txt
View File

@ -88,16 +88,12 @@ else()
message(FATAL_ERROR "Compiler type(CMAKE_Fortran_COMPILER_ID) not recognized")
endif()
file(STRINGS "$ENV{PETSC_DIR}/$ENV{PETSC_ARCH}/lib/petsc/conf/petscvariables" PETSC_EXTERNAL_LIB REGEX "PETSC_WITH_EXTERNAL_LIB = .*$?")
string(REGEX MATCHALL "-[lLW]([^\" ]+)" PETSC_EXTERNAL_LIB "${PETSC_EXTERNAL_LIB}")
list(REMOVE_DUPLICATES PETSC_EXTERNAL_LIB)
string(REPLACE ";" " " PETSC_EXTERNAL_LIB "${PETSC_EXTERNAL_LIB}")
file(STRINGS "$ENV{PETSC_DIR}/$ENV{PETSC_ARCH}/lib/petsc/conf/petscvariables" PETSC_EXTERNAL_LIB REGEX "PETSC_EXTERNAL_LIB_BASIC = .*$?")
string(REPLACE "PETSC_EXTERNAL_LIB_BASIC = " "" PETSC_EXTERNAL_LIB "${PETSC_EXTERNAL_LIB}")
message("PETSC_EXTERNAL_LIB:\n${PETSC_EXTERNAL_LIB}\n")
file(STRINGS "$ENV{PETSC_DIR}/$ENV{PETSC_ARCH}/lib/petsc/conf/petscvariables" PETSC_INCLUDES REGEX "PETSC_FC_INCLUDES = .*$?")
string(REGEX MATCHALL "-I([^\" ]+)" PETSC_INCLUDES "${PETSC_INCLUDES}")
list(REMOVE_DUPLICATES PETSC_INCLUDES)
string(REPLACE ";" " " PETSC_INCLUDES "${PETSC_INCLUDES}")
string(REPLACE "PETSC_FC_INCLUDES = " "" PETSC_INCLUDES "${PETSC_INCLUDES}")
message("PETSC_INCLUDES:\n${PETSC_INCLUDES}\n")
set(CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE} "${BUILDCMD_PRE} ${OPENMP_FLAGS} ${STANDARD_CHECK} ${OPTIMIZATION_FLAGS} ${COMPILE_FLAGS} ${PRECISION_FLAGS}")
@ -109,7 +105,7 @@ if(CMAKE_BUILD_TYPE STREQUAL "DEBUG")
endif()
set(CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE} "${CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE}} ${PETSC_INCLUDES} ${BUILDCMD_POST}")
set(CMAKE_Fortran_LINK_EXECUTABLE "${CMAKE_Fortran_LINK_EXECUTABLE} <OBJECTS> -o <TARGET> <LINK_LIBRARIES> ${PETSC_EXTERNAL_LIB} -lz ${BUILDCMD_POST}")
set(CMAKE_Fortran_LINK_EXECUTABLE "${CMAKE_Fortran_LINK_EXECUTABLE} <OBJECTS> -o <TARGET> <LINK_LIBRARIES> -L${PETSC_LIBRARY_DIRS} -lpetsc ${PETSC_EXTERNAL_LIB} -lz ${BUILDCMD_POST}")
message("Fortran Compiler Flags:\n${CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE}}\n")
message("C Compiler Flags:\n${CMAKE_C_FLAGS_${CMAKE_BUILD_TYPE}}\n")

6
Makefile
View File

@ -10,14 +10,12 @@ all: grid mesh
.PHONY: grid
grid:
@cmake -B build/grid -DDAMASK_SOLVER=grid -DCMAKE_INSTALL_PREFIX=${PWD} -DCMAKE_BUILD_TYPE=${BUILD_TYPE} -DBUILDCMD_POST=${BUILDCMD_POST} -DBUILDCMD_PRE=${BUILDCMD_PRE} -DOPTIMIZATION=${OPTIMIZATION} -DOPENMP=${OPENMP}
@cmake --build build/grid --parallel
@cmake --install build/grid
@cmake --build build/grid --parallel --target install
.PHONY: mesh
mesh:
@cmake -B build/mesh -DDAMASK_SOLVER=mesh -DCMAKE_INSTALL_PREFIX=${PWD} -DCMAKE_BUILD_TYPE=${BUILD_TYPE} -DBUILDCMD_POST=${BUILDCMD_POST} -DBUILDCMD_PRE=${BUILDCMD_PRE} -DOPTIMIZATION=${OPTIMIZATION} -DOPENMP=${OPENMP}
@cmake --build build/mesh --parallel
@cmake --install build/mesh
@cmake --build build/mesh --parallel --target install
.PHONY: clean
clean:

2
PRIVATE

@ -1 +1 @@
Subproject commit 96c32ba4237a51eaad92cd139e1a716ee5b32493
Subproject commit b898a8b5552bd9d1c555edc3d8134564dd32fe53

5
examples/config/Phase_Phenopowerlaw_BCC-Martensite.yaml
View File

@ -1,11 +1,6 @@
# Tasan et.al. 2015 Acta Materalia
# Tasan et.al. 2015 International Journal of Plasticity
# Diehl et.al. 2015 Meccanica
Martensite:
lattice: cI
mechanical:
elastic: {C_11: 417.4e+9, C_12: 242.4e+9, C_44: 211.1e+9, type: Hooke}
plastic:
N_sl: [12, 12]
a_sl: 2.0
dot_gamma_0_sl: 0.001

6
examples/config/phase/AISI304.yaml Normal file
View File

@ -0,0 +1,6 @@
references:
- H.M. Ledbetter
physica status solidi (a) 85(1):89-96, 1984
https://doi.org/10.1002/pssa.2210850111
lattice: cF
rho: 7937.0

7
examples/config/phase/mechanical/eigen/thermalexpansion_AISI304.yaml Normal file
View File

@ -0,0 +1,7 @@
type: thermalexpansion
references:
- R.H. Bogaard et al.
Thermochimica Acta 218:373-393, 1993
https://doi.org/10.1016/0040-6031(93)80437-F
A_11: 15.0e-6
T_ref: 300.0

8
examples/config/phase/mechanical/elastic/Hooke_AISI304.yaml Normal file
View File

@ -0,0 +1,8 @@
type: Hooke
references:
- H.M. Ledbetter
physica status solidi (a) 85(1):89-96, 1984
https://doi.org/10.1002/pssa.2210850111
C_11: 204.6e+9
C_12: 137.7e+9
C_44: 126.2e+9

8
examples/config/phase/mechanical/elastic/Hooke_SAE1050-martensite.yaml Normal file
View File

@ -0,0 +1,8 @@
type: Hooke
references:
- S.A. Kim and W.L. Johnson,
Materials Science & Engineering A 452-453:633-639, 2007,
https://doi.org/10.1016/j.msea.2006.11.147
C_11: 268.1e+9
C_12: 111.2e+9
C_44: 79.06e+9

3
examples/config/phase/mechanical/plastic/dislotwin_IF-steel.yaml
View File

@ -4,7 +4,8 @@ references:
International Journal of Plasticity 134:102779, 2020,
https://doi.org/10.1016/j.ijplas.2020.102779
- K. Sedighiani et al.,
Mechanics of Materials, submitted
Mechanics of Materials, 164:104117, 2022,
https://doi.org/10.1016/j.mechmat.2021.104117
output: [rho_dip, rho_mob]
N_sl: [12, 12]
b_sl: [2.49e-10, 2.49e-10]

9
examples/config/phase/thermal/AISI304.yaml Normal file
View File

@ -0,0 +1,9 @@
references:
- B.F. Blackwell et al.
Proceedings of 34th National Heat Transfer Conference 2000
https://www.osti.gov/servlets/purl/760791
- R.H. Bogaard et al.
Thermochimica Acta 218:373-393, 1993
https://doi.org/10.1016/0040-6031(93)80437-F
C_p: 470.0
K_11: 14.34

0
examples/config/phase/thermal/Steel-0.5C.yaml → examples/config/phase/thermal/steel-0.5C.yaml
View File

6
install/MarcMentat/apply_DAMASK_modifications.py
View File

@ -67,9 +67,7 @@ os.system(f'xvfb-run -a {executable} -compile {menu_file}')
print('setting file access rights...')
files = (glob.glob(str(marc_root/f'marc{marc_version}/tools/*_damask*')) +
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]')))
for file in files:
glob.glob(str(marc_root/f'mentat{marc_version}/bin/submit[4-6]'))):
os.chmod(file , 0o755)

71
processing/legacy/addDisplacement.py vendored
View File

@ -1,71 +0,0 @@
#!/usr/bin/env python3
import os
import sys
from io import StringIO
from optparse import OptionParser
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
# --------------------------------------------------------------------
# MAIN
# --------------------------------------------------------------------
parser = OptionParser(usage='%prog options [ASCIItable(s)]', description = """
Add displacments resulting from deformation gradient field.
Operates on periodic three-dimensional x,y,z-ordered data sets.
Outputs at cell centers or cell nodes (into separate file).
""", version = scriptID)
parser.add_option('-f',
'--defgrad',
dest = 'f',
metavar = 'string',
help = 'label of deformation gradient [%default]')
parser.add_option('-p',
'--pos', '--position',
dest = 'pos',
metavar = 'string',
help = 'label of coordinates [%default]')
parser.add_option('--nodal',
dest = 'nodal',
action = 'store_true',
help = 'output nodal (instead of cell-centered) displacements')
parser.set_defaults(f = 'f',
pos = 'pos',
)
(options,filenames) = parser.parse_args()
for name in filenames:
damask.util.report(scriptName,name)
table = damask.Table.load(StringIO(''.join(sys.stdin.read())) if name is None else name)
grid,size,origin = damask.grid_filters.cellsSizeOrigin_coordinates0_point(table.get(options.pos))
F = table.get(options.f).reshape(tuple(grid)+(-1,),order='F').reshape(tuple(grid)+(3,3))
if options.nodal:
damask.Table(damask.grid_filters.coordinates0_node(grid,size).reshape(-1,3,order='F'),
{'pos':(3,)})\
.add('avg({}).{}'.format(options.f,options.pos),
damask.grid_filters.displacement_avg_node(size,F).reshape(-1,3,order='F'),
scriptID+' '+' '.join(sys.argv[1:]))\
.add('fluct({}).{}'.format(options.f,options.pos),
damask.grid_filters.displacement_fluct_node(size,F).reshape(-1,3,order='F'),
scriptID+' '+' '.join(sys.argv[1:]))\
.save((sys.stdout if name is None else os.path.splitext(name)[0]+'_nodal.txt'))
else:
table.add('avg({}).{}'.format(options.f,options.pos),
damask.grid_filters.displacement_avg_point(size,F).reshape(-1,3,order='F'),
scriptID+' '+' '.join(sys.argv[1:]))\
.add('fluct({}).{}'.format(options.f,options.pos),
damask.grid_filters.displacement_fluct_point(size,F).reshape(-1,3,order='F'),
scriptID+' '+' '.join(sys.argv[1:]))\
.save((sys.stdout if name is None else name))

2
python/damask/VERSION
View File

@ -1 +1 @@
v3.0.0-alpha5-333-g01cd92755
v3.0.0-alpha5-375-g76fe2d2b3

1
python/damask/__init__.py
View File

@ -8,6 +8,7 @@ with open(_Path(__file__).parent/_Path('VERSION')) as _f:
version = _re.sub(r'^v','',_f.readline().strip())
__version__ = version
from . import _typehints # noqa
from . import util # noqa
from . import seeds # noqa
from . import tensor # noqa

38
python/damask/_colormap.py
View File

@ -3,13 +3,9 @@ import json
import functools
import colorsys
from pathlib import Path
from typing import Sequence, Union, TextIO
from typing import Union, TextIO
import numpy as np
try:
from numpy.typing import ArrayLike
except ImportError:
ArrayLike = Union[np.ndarray,Sequence[float]] # type: ignore
import scipy.interpolate as interp
import matplotlib as mpl
if os.name == 'posix' and 'DISPLAY' not in os.environ:
@ -18,6 +14,7 @@ import matplotlib.pyplot as plt
from matplotlib import cm
from PIL import Image
from ._typehints import FloatSequence, FileHandle
from . import util
from . import Table
@ -82,8 +79,8 @@ class Colormap(mpl.colors.ListedColormap):
@staticmethod
def from_range(low: ArrayLike,
high: ArrayLike,
def from_range(low: FloatSequence,
high: FloatSequence,
name: str = 'DAMASK colormap',
N: int = 256,
model: str = 'rgb') -> 'Colormap':
@ -197,7 +194,7 @@ class Colormap(mpl.colors.ListedColormap):
def at(self,
fraction : Union[float,Sequence[float]]) -> np.ndarray:
fraction : Union[float,FloatSequence]) -> np.ndarray:
"""
Interpolate color at fraction.
@ -229,14 +226,14 @@ class Colormap(mpl.colors.ListedColormap):
def shade(self,
field: np.ndarray,
bounds: ArrayLike = None,
bounds: FloatSequence = None,
gap: float = None) -> Image:
"""
Generate PIL image of 2D field using colormap.
Parameters
----------
field : numpy.array, shape (:,:)
field : numpy.ndarray, shape (:,:)
Data to be shaded.
bounds : sequence of float, len (2), optional
Value range (left,right) spanned by colormap.
@ -296,7 +293,7 @@ class Colormap(mpl.colors.ListedColormap):
def _get_file_handle(self,
fname: Union[TextIO, str, Path, None],
fname: Union[FileHandle, None],
suffix: str = '') -> TextIO:
"""
Provide file handle.
@ -323,7 +320,7 @@ class Colormap(mpl.colors.ListedColormap):
return fname
def save_paraview(self, fname: Union[TextIO, str, Path] = None):
def save_paraview(self, fname: FileHandle = None):
"""
Save as JSON file for use in Paraview.
@ -350,7 +347,7 @@ class Colormap(mpl.colors.ListedColormap):
fhandle.write('\n')
def save_ASCII(self, fname: Union[TextIO, str, Path] = None):
def save_ASCII(self, fname: FileHandle = None):
"""
Save as ASCII file.
@ -365,7 +362,7 @@ class Colormap(mpl.colors.ListedColormap):
t.save(self._get_file_handle(fname,'.txt'))
def save_GOM(self, fname: Union[TextIO, str, Path] = None):
def save_GOM(self, fname: FileHandle = None):
"""
Save as ASCII file for use in GOM Aramis.
@ -385,7 +382,7 @@ class Colormap(mpl.colors.ListedColormap):
self._get_file_handle(fname,'.legend').write(GOM_str)
def save_gmsh(self, fname: Union[TextIO, str, Path] = None):
def save_gmsh(self, fname: FileHandle = None):
"""
Save as ASCII file for use in gmsh.
@ -616,7 +613,7 @@ class Colormap(mpl.colors.ListedColormap):
@staticmethod
def _lab2xyz(lab: np.ndarray, ref_white: np.ndarray = None) -> np.ndarray:
def _lab2xyz(lab: np.ndarray, ref_white: np.ndarray = _REF_WHITE) -> np.ndarray:
"""
CIE Lab to CIE Xyz.
@ -624,6 +621,8 @@ class Colormap(mpl.colors.ListedColormap):
----------
lab : numpy.ndarray, shape (3)
CIE lab values.
ref_white : numpy.ndarray, shape (3)
Reference white, default value is the standard 2° observer for D65.
Returns
-------
@ -642,10 +641,10 @@ class Colormap(mpl.colors.ListedColormap):
f_x**3. if f_x**3. > _EPS else (116.*f_x-16.)/_KAPPA,
((lab[0]+16.)/116.)**3 if lab[0]>_KAPPA*_EPS else lab[0]/_KAPPA,
f_z**3. if f_z**3. > _EPS else (116.*f_z-16.)/_KAPPA
])*(ref_white if ref_white is not None else _REF_WHITE)
])*ref_white
@staticmethod
def _xyz2lab(xyz: np.ndarray, ref_white: np.ndarray = None) -> np.ndarray:
def _xyz2lab(xyz: np.ndarray, ref_white: np.ndarray = _REF_WHITE) -> np.ndarray:
"""
CIE Xyz to CIE Lab.
@ -653,6 +652,8 @@ class Colormap(mpl.colors.ListedColormap):
----------
xyz : numpy.ndarray, shape (3)
CIE Xyz values.
ref_white : numpy.ndarray, shape (3)
Reference white, default value is the standard 2° observer for D65.
Returns
-------
@ -664,7 +665,6 @@ class Colormap(mpl.colors.ListedColormap):
http://www.brucelindbloom.com/index.html?Eqn_Lab_to_XYZ.html
"""
ref_white = ref_white if ref_white is not None else _REF_WHITE
f = np.where(xyz/ref_white > _EPS,(xyz/ref_white)**(1./3.),(_KAPPA*xyz/ref_white+16.)/116.)
return np.array([

80
python/damask/_crystal.py
View File

@ -114,12 +114,13 @@ class Crystal():
def __repr__(self):
"""Represent."""
return '\n'.join([f'Crystal family {self.family}']
+ ([] if self.lattice is None else [f'Bravais lattice {self.lattice}']+
list(map(lambda x:f'{x[0]}: {x[1]:.5g}',
zip(['a','b','c','α','β','γ',],
self.parameters))))
)
family = f'Crystal family: {self.family}'
return family if self.lattice is None else \
'\n'.join([family,
f'Bravais lattice: {self.lattice}',
'a={:.5g}m, b={:.5g}m, c={:.5g}m'.format(*self.parameters[:3]),
'α={:.5g}°, β={:.5g}°, γ={:.5g}°'.format(*np.degrees(self.parameters[3:]))])
def __eq__(self,other):
"""
@ -543,6 +544,73 @@ class Crystal():
[ 1,-2, 1,-3, 1,-2, 1, 2],
[ 2,-1,-1,-3, 2,-1,-1, 2]])]
},
'tI': {
'slip': [np.array([
[+0,+0,+1, +1,+0,+0],
[+0,+0,+1, +0,+1,+0]]),
np.array([
[+0,+0,+1, +1,+1,+0],
[+0,+0,+1, -1,+1,+0]]),
np.array([
[+0,+1,+0, +1,+0,+0],
[+1,+0,+0, +0,+1,+0]]),
np.array([
[+1,-1,+1, +1,+1,+0],
[+1,-1,-1, +1,+1,+0],
[-1,-1,-1, -1,+1,+0],
[-1,-1,+1, -1,+1,+0]]),
np.array([
[+1,-1,+0, +1,+1,+0],
[+1,+1,+0, +1,-1,+0]]),
np.array([
[+0,+1,+1, +1,+0,+0],
[+0,-1,+1, +1,+0,+0],
[-1,+0,+1, +0,+1,+0],
[+1,+0,+1, +0,+1,+0]]),
np.array([
[+0,+1,+0, +0,+0,+1],
[+1,+0,+0, +0,+0,+1]]),
np.array([
[+1,+1,+0, +0,+0,+1],
[-1,+1,+0, +0,+0,+1]]),
np.array([
[+0,+1,-1, +0,+1,+1],
[+0,-1,-1, +0,-1,+1],
[-1,+0,-1, -1,+0,+1],
[+1,+0,-1, +1,+0,+1]]),
np.array([
[+1,-1,+1, +0,+1,+1],
[+1,+1,-1, +0,+1,+1],
[+1,+1,+1, +0,+1,-1],
[-1,+1,+1, +0,+1,-1],
[+1,-1,-1, +1,+0,+1],
[-1,-1,+1, +1,+0,+1],
[+1,+1,+1, +1,+0,-1],
[+1,-1,+1, +1,+0,-1]]),
np.array([
[+1,+0,+0, +0,+1,+1],
[+1,+0,+0, +0,+1,-1],
[+0,+1,+0, +1,+0,+1],
[+0,+1,+0, +1,+0,-1]]),
np.array([
[+0,+1,-1, +2,+1,+1],
[+0,-1,-1, +2,-1,+1],
[+1,+0,-1, +1,+2,+1],
[-1,+0,-1, -1,+2,+1],
[+0,+1,-1, -2,+1,+1],
[+0,-1,-1, -2,-1,+1],
[-1,+0,-1, -1,-2,+1],
[+1,+0,-1, +1,-2,+1]]),
np.array([
[-1,+1,+1, +2,+1,+1],
[-1,-1,+1, +2,-1,+1],
[+1,-1,+1, +1,+2,+1],
[-1,-1,+1, -1,+2,+1],
[+1,+1,+1, -2,+1,+1],
[+1,-1,+1, -2,-1,+1],
[-1,+1,+1, -1,-2,+1],
[+1,+1,+1, +1,-2,+1]])]
}
}
master = _kinematics[self.lattice][mode]
if self.lattice == 'hP':

11
python/damask/_orientation.py
View File

@ -514,6 +514,17 @@ class Orientation(Rotation,Crystal):
[ 0.07359167 -0.36505797 0.92807163]]
Bunge Eulers / deg: (11.40, 21.86, 0.60)
Plot a sample from the Mackenzie distribution.
>>> import matplotlib.pyplot as plt
>>> import damask
>>> N = 10000
>>> a = damask.Orientation.from_random(shape=N,family='cubic')
>>> b = damask.Orientation.from_random(shape=N,family='cubic')
>>> d = a.disorientation(b).as_axis_angle(degrees=True,pair=True)[1]
>>> plt.hist(d,25)
>>> plt.show()
"""
if self.family != other.family:
raise NotImplementedError('disorientation between different crystal families')

11
python/damask/_typehints.py Normal file
View File

@ -0,0 +1,11 @@
"""Functionality for typehints."""
from typing import Sequence, Union, TextIO
from pathlib import Path
import numpy as np
FloatSequence = Union[np.ndarray,Sequence[float]]
IntSequence = Union[np.ndarray,Sequence[int]]
FileHandle = Union[TextIO, str, Path]

169
python/damask/grid_filters.py
View File

@ -12,21 +12,23 @@ the following operations are required for tensorial data:
"""
from typing import Sequence, Tuple, Union
from typing import Tuple as _Tuple
from scipy import spatial as _spatial
import numpy as _np
from ._typehints import FloatSequence as _FloatSequence, IntSequence as _IntSequence
def _ks(size: _np.ndarray, cells: Union[_np.ndarray,Sequence[int]], first_order: bool = False) -> _np.ndarray:
def _ks(size: _FloatSequence, cells: _IntSequence, first_order: bool = False) -> _np.ndarray:
"""
Get wave numbers operator.
Parameters
----------
size : numpy.ndarray of shape (3)
size : sequence of float, len (3)
Physical size of the periodic field.
cells : numpy.ndarray of shape (3)
cells : sequence of int, len (3)
Number of cells.
first_order : bool, optional
Correction for first order derivatives, defaults to False.
@ -45,20 +47,20 @@ def _ks(size: _np.ndarray, cells: Union[_np.ndarray,Sequence[int]], first_order:
return _np.stack(_np.meshgrid(k_sk,k_sj,k_si,indexing = 'ij'), axis=-1)
def curl(size: _np.ndarray, f: _np.ndarray) -> _np.ndarray:
def curl(size: _FloatSequence, f: _np.ndarray) -> _np.ndarray:
u"""
Calculate curl of a vector or tensor field in Fourier space.
Parameters
----------
size : numpy.ndarray of shape (3)
size : sequence of float, len (3)
Physical size of the periodic field.
f : numpy.ndarray of shape (:,:,:,3) or (:,:,:,3,3)
f : numpy.ndarray, shape (:,:,:,3) or (:,:,:,3,3)
Periodic field of which the curl is calculated.
Returns
-------
∇ × f : numpy.ndarray
∇ × f : numpy.ndarray, shape (:,:,:,3) or (:,:,:,3,3)
Curl of f.
"""
@ -76,20 +78,20 @@ def curl(size: _np.ndarray, f: _np.ndarray) -> _np.ndarray:
return _np.fft.irfftn(curl_,axes=(0,1,2),s=f.shape[:3])
def divergence(size: _np.ndarray, f: _np.ndarray) -> _np.ndarray:
def divergence(size: _FloatSequence, f: _np.ndarray) -> _np.ndarray:
u"""
Calculate divergence of a vector or tensor field in Fourier space.
Parameters
----------
size : numpy.ndarray of shape (3)
size : sequence of float, len (3)
Physical size of the periodic field.
f : numpy.ndarray of shape (:,:,:,3) or (:,:,:,3,3)
f : numpy.ndarray, shape (:,:,:,3) or (:,:,:,3,3)
Periodic field of which the divergence is calculated.
Returns
-------
∇ · f : numpy.ndarray
∇ · f : numpy.ndarray, shape (:,:,:,1) or (:,:,:,3)
Divergence of f.
"""
@ -103,20 +105,20 @@ def divergence(size: _np.ndarray, f: _np.ndarray) -> _np.ndarray:
return _np.fft.irfftn(div_,axes=(0,1,2),s=f.shape[:3])
def gradient(size: _np.ndarray, f: _np.ndarray) -> _np.ndarray:
def gradient(size: _FloatSequence, f: _np.ndarray) -> _np.ndarray:
u"""
Calculate gradient of a scalar or vector field in Fourier space.
Parameters
----------
size : numpy.ndarray of shape (3)
size : sequence of float, len (3)
Physical size of the periodic field.
f : numpy.ndarray of shape (:,:,:,1) or (:,:,:,3)
f : numpy.ndarray, shape (:,:,:,1) or (:,:,:,3)
Periodic field of which the gradient is calculated.
Returns
-------
∇ f : numpy.ndarray
∇ f : numpy.ndarray, shape (:,:,:,3) or (:,:,:,3,3)
Divergence of f.
"""
@ -130,29 +132,30 @@ def gradient(size: _np.ndarray, f: _np.ndarray) -> _np.ndarray:
return _np.fft.irfftn(grad_,axes=(0,1,2),s=f.shape[:3])
def coordinates0_point(cells: Union[ _np.ndarray,Sequence[int]],
size: _np.ndarray,
origin: _np.ndarray = _np.zeros(3)) -> _np.ndarray:
def coordinates0_point(cells: _IntSequence,
size: _FloatSequence,
origin: _FloatSequence = _np.zeros(3)) -> _np.ndarray:
"""
Cell center positions (undeformed).
Parameters
----------
cells : numpy.ndarray of shape (3)
cells : sequence of int, len (3)
Number of cells.
size : numpy.ndarray of shape (3)
size : sequence of float, len (3)
Physical size of the periodic field.
origin : numpy.ndarray, optional
origin : sequence of float, len(3), optional
Physical origin of the periodic field. Defaults to [0.0,0.0,0.0].
Returns
-------
x_p_0 : numpy.ndarray
x_p_0 : numpy.ndarray, shape (:,:,:,3)
Undeformed cell center coordinates.
"""
start = origin + size/_np.array(cells)*.5
end = origin + size - size/_np.array(cells)*.5
size_ = _np.array(size,float)
start = origin + size_/_np.array(cells,int)*.5
end = origin + size_ - size_/_np.array(cells,int)*.5
return _np.stack(_np.meshgrid(_np.linspace(start[0],end[0],cells[0]),
_np.linspace(start[1],end[1],cells[1]),
@ -160,24 +163,24 @@ def coordinates0_point(cells: Union[ _np.ndarray,Sequence[int]],
axis = -1)
def displacement_fluct_point(size: _np.ndarray, F: _np.ndarray) -> _np.ndarray:
def displacement_fluct_point(size: _FloatSequence, F: _np.ndarray) -> _np.ndarray:
"""
Cell center displacement field from fluctuation part of the deformation gradient field.
Parameters
----------
size : numpy.ndarray of shape (3)
size : sequence of float, len (3)
Physical size of the periodic field.
F : numpy.ndarray
F : numpy.ndarray, shape (:,:,:,3,3)
Deformation gradient field.
Returns
-------
u_p_fluct : numpy.ndarray
u_p_fluct : numpy.ndarray, shape (:,:,:,3)
Fluctuating part of the cell center displacements.
"""
integrator = 0.5j*size/_np.pi
integrator = 0.5j*_np.array(size,float)/_np.pi
k_s = _ks(size,F.shape[:3],False)
k_s_squared = _np.einsum('...l,...l',k_s,k_s)
@ -192,20 +195,20 @@ def displacement_fluct_point(size: _np.ndarray, F: _np.ndarray) -> _np.ndarray:
return _np.fft.irfftn(displacement,axes=(0,1,2),s=F.shape[:3])
def displacement_avg_point(size: _np.ndarray, F: _np.ndarray) -> _np.ndarray:
def displacement_avg_point(size: _FloatSequence, F: _np.ndarray) -> _np.ndarray:
"""
Cell center displacement field from average part of the deformation gradient field.
Parameters
----------
size : numpy.ndarray of shape (3)
size : sequence of float, len (3)
Physical size of the periodic field.
F : numpy.ndarray
F : numpy.ndarray, shape (:,:,:,3,3)
Deformation gradient field.
Returns
-------
u_p_avg : numpy.ndarray
u_p_avg : numpy.ndarray, shape (:,:,:,3)
Average part of the cell center displacements.
"""
@ -213,42 +216,42 @@ def displacement_avg_point(size: _np.ndarray, F: _np.ndarray) -> _np.ndarray:
return _np.einsum('ml,ijkl->ijkm',F_avg - _np.eye(3),coordinates0_point(F.shape[:3],size))
def displacement_point(size: _np.ndarray, F: _np.ndarray) -> _np.ndarray:
def displacement_point(size: _FloatSequence, F: _np.ndarray) -> _np.ndarray:
"""
Cell center displacement field from deformation gradient field.
Parameters
----------
size : numpy.ndarray of shape (3)
size : sequence of float, len (3)
Physical size of the periodic field.
F : numpy.ndarray
F : numpy.ndarray, shape (:,:,:,3,3)
Deformation gradient field.
Returns
-------
u_p : numpy.ndarray
u_p : numpy.ndarray, shape (:,:,:,3)
Cell center displacements.
"""
return displacement_avg_point(size,F) + displacement_fluct_point(size,F)
def coordinates_point(size: _np.ndarray, F: _np.ndarray, origin: _np.ndarray = _np.zeros(3)) -> _np.ndarray:
def coordinates_point(size: _FloatSequence, F: _np.ndarray, origin: _FloatSequence = _np.zeros(3)) -> _np.ndarray:
"""
Cell center positions.
Parameters
----------
size : numpy.ndarray of shape (3)
size : sequence of float, len (3)
Physical size of the periodic field.
F : numpy.ndarray
F : numpy.ndarray, shape (:,:,:,3,3)
Deformation gradient field.
origin : numpy.ndarray of shape (3), optional
origin : sequence of float, len(3), optional
Physical origin of the periodic field. Defaults to [0.0,0.0,0.0].
Returns
-------
x_p : numpy.ndarray
x_p : numpy.ndarray, shape (:,:,:,3)
Cell center coordinates.
"""
@ -256,14 +259,14 @@ def coordinates_point(size: _np.ndarray, F: _np.ndarray, origin: _np.ndarray = _
def cellsSizeOrigin_coordinates0_point(coordinates0: _np.ndarray,
ordered: bool = True) -> Tuple[_np.ndarray,_np.ndarray,_np.ndarray]:
ordered: bool = True) -> _Tuple[_np.ndarray,_np.ndarray,_np.ndarray]:
"""
Return grid 'DNA', i.e. cells, size, and origin from 1D array of point positions.
Parameters
----------
coordinates0 : numpy.ndarray of shape (:,3)
Undeformed cell coordinates.
coordinates0 : numpy.ndarray, shape (:,3)
Undeformed cell center coordinates.
ordered : bool, optional
Expect coordinates0 data to be ordered (x fast, z slow).
Defaults to True.
@ -277,7 +280,7 @@ def cellsSizeOrigin_coordinates0_point(coordinates0: _np.ndarray,
coords = [_np.unique(coordinates0[:,i]) for i in range(3)]
mincorner = _np.array(list(map(min,coords)))
maxcorner = _np.array(list(map(max,coords)))
cells = _np.array(list(map(len,coords)),'i')
cells = _np.array(list(map(len,coords)),int)
size = cells/_np.maximum(cells-1,1) * (maxcorner-mincorner)
delta = size/cells
origin = mincorner - delta*.5
@ -305,24 +308,24 @@ def cellsSizeOrigin_coordinates0_point(coordinates0: _np.ndarray,
return (cells,size,origin)
def coordinates0_node(cells: Union[_np.ndarray,Sequence[int]],
size: _np.ndarray,
origin: _np.ndarray = _np.zeros(3)) -> _np.ndarray:
def coordinates0_node(cells: _IntSequence,
size: _FloatSequence,
origin: _FloatSequence = _np.zeros(3)) -> _np.ndarray:
"""
Nodal positions (undeformed).
Parameters
----------
cells : numpy.ndarray of shape (3)
cells : sequence of int, len (3)
Number of cells.
size : numpy.ndarray of shape (3)
size : sequence of float, len (3)
Physical size of the periodic field.
origin : numpy.ndarray of shape (3), optional
origin : sequence of float, len(3), optional
Physical origin of the periodic field. Defaults to [0.0,0.0,0.0].
Returns
-------
x_n_0 : numpy.ndarray
x_n_0 : numpy.ndarray, shape (:,:,:,3)
Undeformed nodal coordinates.
"""
@ -332,40 +335,40 @@ def coordinates0_node(cells: Union[_np.ndarray,Sequence[int]],
axis = -1)
def displacement_fluct_node(size: _np.ndarray, F: _np.ndarray) -> _np.ndarray:
def displacement_fluct_node(size: _FloatSequence, F: _np.ndarray) -> _np.ndarray:
"""
Nodal displacement field from fluctuation part of the deformation gradient field.
Parameters
----------
size : numpy.ndarray of shape (3)
size : sequence of float, len (3)
Physical size of the periodic field.
F : numpy.ndarray
F : numpy.ndarray, shape (:,:,:,3,3)
Deformation gradient field.
Returns
-------
u_n_fluct : numpy.ndarray
u_n_fluct : numpy.ndarray, shape (:,:,:,3)
Fluctuating part of the nodal displacements.
"""
return point_to_node(displacement_fluct_point(size,F))
def displacement_avg_node(size: _np.ndarray, F: _np.ndarray) -> _np.ndarray:
def displacement_avg_node(size: _FloatSequence, F: _np.ndarray) -> _np.ndarray:
"""
Nodal displacement field from average part of the deformation gradient field.
Parameters
----------
size : numpy.ndarray of shape (3)
size : sequence of float, len (3)
Physical size of the periodic field.
F : numpy.ndarray
F : numpy.ndarray, shape (:,:,:,3,3)
Deformation gradient field.
Returns
-------
u_n_avg : numpy.ndarray
u_n_avg : numpy.ndarray, shape (:,:,:,3)
Average part of the nodal displacements.
"""
@ -373,42 +376,42 @@ def displacement_avg_node(size: _np.ndarray, F: _np.ndarray) -> _np.ndarray:
return _np.einsum('ml,ijkl->ijkm',F_avg - _np.eye(3),coordinates0_node(F.shape[:3],size))
def displacement_node(size: _np.ndarray, F: _np.ndarray) -> _np.ndarray:
def displacement_node(size: _FloatSequence, F: _np.ndarray) -> _np.ndarray:
"""
Nodal displacement field from deformation gradient field.
Parameters
----------
size : numpy.ndarray of shape (3)
size : sequence of float, len (3)
Physical size of the periodic field.
F : numpy.ndarray
F : numpy.ndarray, shape (:,:,:,3,3)
Deformation gradient field.
Returns
-------
u_p : numpy.ndarray
u_p : numpy.ndarray, shape (:,:,:,3)
Nodal displacements.
"""
return displacement_avg_node(size,F) + displacement_fluct_node(size,F)
def coordinates_node(size: _np.ndarray, F: _np.ndarray, origin: _np.ndarray = _np.zeros(3)) -> _np.ndarray:
def coordinates_node(size: _FloatSequence, F: _np.ndarray, origin: _FloatSequence = _np.zeros(3)) -> _np.ndarray:
"""
Nodal positions.
Parameters
----------
size : numpy.ndarray of shape (3)
size : sequence of float, len (3)
Physical size of the periodic field.
F : numpy.ndarray
F : numpy.ndarray, shape (:,:,:,3,3)
Deformation gradient field.
origin : numpy.ndarray of shape (3), optional
origin : sequence of float, len(3), optional
Physical origin of the periodic field. Defaults to [0.0,0.0,0.0].
Returns
-------
x_n : numpy.ndarray
x_n : numpy.ndarray, shape (:,:,:,3)
Nodal coordinates.
"""
@ -416,13 +419,13 @@ def coordinates_node(size: _np.ndarray, F: _np.ndarray, origin: _np.ndarray = _n
def cellsSizeOrigin_coordinates0_node(coordinates0: _np.ndarray,
ordered: bool = True) -> Tuple[_np.ndarray,_np.ndarray,_np.ndarray]:
ordered: bool = True) -> _Tuple[_np.ndarray,_np.ndarray,_np.ndarray]:
"""
Return grid 'DNA', i.e. cells, size, and origin from 1D array of nodal positions.
Parameters
----------
coordinates0 : numpy.ndarray of shape (:,3)
coordinates0 : numpy.ndarray, shape (:,3)
Undeformed nodal coordinates.
ordered : bool, optional
Expect coordinates0 data to be ordered (x fast, z slow).
@ -437,7 +440,7 @@ def cellsSizeOrigin_coordinates0_node(coordinates0: _np.ndarray,
coords = [_np.unique(coordinates0[:,i]) for i in range(3)]
mincorner = _np.array(list(map(min,coords)))
maxcorner = _np.array(list(map(max,coords)))
cells = _np.array(list(map(len,coords)),'i') - 1
cells = _np.array(list(map(len,coords)),int) - 1
size = maxcorner-mincorner
origin = mincorner
@ -463,12 +466,12 @@ def point_to_node(cell_data: _np.ndarray) -> _np.ndarray:
Parameters
----------
cell_data : numpy.ndarray of shape (:,:,:,...)
cell_data : numpy.ndarray, shape (:,:,:,...)
Data defined on the cell centers of a periodic grid.
Returns
-------
node_data : numpy.ndarray of shape (:,:,:,...)
node_data : numpy.ndarray, shape (:,:,:,...)
Data defined on the nodes of a periodic grid.
"""
@ -485,12 +488,12 @@ def node_to_point(node_data: _np.ndarray) -> _np.ndarray:
Parameters
----------
node_data : numpy.ndarray of shape (:,:,:,...)
node_data : numpy.ndarray, shape (:,:,:,...)
Data defined on the nodes of a periodic grid.
Returns
-------
cell_data : numpy.ndarray of shape (:,:,:,...)
cell_data : numpy.ndarray, shape (:,:,:,...)
Data defined on the cell centers of a periodic grid.
"""
@ -507,7 +510,7 @@ def coordinates0_valid(coordinates0: _np.ndarray) -> bool:
Parameters
----------
coordinates0 : numpy.ndarray
coordinates0 : numpy.ndarray, shape (:,3)
Array of undeformed cell coordinates.
Returns
@ -523,17 +526,17 @@ def coordinates0_valid(coordinates0: _np.ndarray) -> bool:
return False
def regrid(size: _np.ndarray, F: _np.ndarray, cells: Union[_np.ndarray,Sequence[int]]) -> _np.ndarray:
def regrid(size: _FloatSequence, F: _np.ndarray, cells: _IntSequence) -> _np.ndarray:
"""
Return mapping from coordinates in deformed configuration to a regular grid.
Parameters
----------
size : numpy.ndarray of shape (3)
size : sequence of float, len (3)
Physical size.
F : numpy.ndarray of shape (:,:,:,3,3)
F : numpy.ndarray, shape (:,:,:,3,3), shape (:,:,:,3,3)
Deformation gradient field.
cells : numpy.ndarray of shape (3)
cells : sequence of int, len (3)
Cell count along x,y,z of remapping grid.
"""

4
python/damask/mechanics.py
View File

@ -5,7 +5,7 @@ All routines operate on numpy.ndarrays of shape (...,3,3).
"""
from typing import Sequence
from typing import Sequence as _Sequence
import numpy as _np
@ -243,7 +243,7 @@ def stretch_right(T: _np.ndarray) -> _np.ndarray:
return _polar_decomposition(T,'U')[0]
def _polar_decomposition(T: _np.ndarray, requested: Sequence[str]) -> tuple:
def _polar_decomposition(T: _np.ndarray, requested: _Sequence[str]) -> tuple:
"""
Perform singular value decomposition.

33
python/damask/seeds.py
View File

@ -1,25 +1,27 @@
"""Functionality for generation of seed points for Voronoi or Laguerre tessellation."""
from typing import Sequence,Tuple
from typing import Tuple as _Tuple
from scipy import spatial as _spatial
import numpy as _np
from ._typehints import FloatSequence as _FloatSequence, IntSequence as _IntSequence
from . import util as _util
from . import grid_filters as _grid_filters
def from_random(size: _np.ndarray, N_seeds: int, cells: _np.ndarray = None, rng_seed=None) -> _np.ndarray:
def from_random(size: _FloatSequence, N_seeds: int, cells: _IntSequence = None,
rng_seed=None) -> _np.ndarray:
"""
Place seeds randomly in space.
Parameters
----------
size : numpy.ndarray of shape (3)
size : sequence of float, len (3)
Physical size of the seeding domain.
N_seeds : int
Number of seeds.
cells : numpy.ndarray of shape (3), optional.
cells : sequence of int, len (3), optional.
If given, ensures that each seed results in a grain when a standard Voronoi
tessellation is performed using the given grid resolution (i.e. size/cells).
rng_seed : {None, int, array_like[ints], SeedSequence, BitGenerator, Generator}, optional
@ -28,29 +30,30 @@ def from_random(size: _np.ndarray, N_seeds: int, cells: _np.ndarray = None, rng_
Returns
-------
coords : numpy.ndarray of shape (N_seeds,3)
coords : numpy.ndarray, shape (N_seeds,3)
Seed coordinates in 3D space.
"""
size_ = _np.array(size,float)
rng = _np.random.default_rng(rng_seed)
if cells is None:
coords = rng.random((N_seeds,3)) * size
coords = rng.random((N_seeds,3)) * size_
else:
grid_coords = _grid_filters.coordinates0_point(cells,size).reshape(-1,3,order='F')
coords = grid_coords[rng.choice(_np.prod(cells),N_seeds, replace=False)] \
+ _np.broadcast_to(size/cells,(N_seeds,3))*(rng.random((N_seeds,3))*.5-.25) # wobble without leaving cells
+ _np.broadcast_to(size_/_np.array(cells,int),(N_seeds,3))*(rng.random((N_seeds,3))*.5-.25) # wobble w/o leaving grid
return coords
def from_Poisson_disc(size: _np.ndarray, N_seeds: int, N_candidates: int, distance: float,
def from_Poisson_disc(size: _FloatSequence, N_seeds: int, N_candidates: int, distance: float,
periodic: bool = True, rng_seed=None) -> _np.ndarray:
"""
Place seeds according to a Poisson disc distribution.
Parameters
----------
size : numpy.ndarray of shape (3)
size : sequence of float, len (3)
Physical size of the seeding domain.
N_seeds : int
Number of seeds.
@ -66,13 +69,13 @@ def from_Poisson_disc(size: _np.ndarray, N_seeds: int, N_candidates: int, distan
Returns
-------
coords : numpy.ndarray of shape (N_seeds,3)
coords : numpy.ndarray, shape (N_seeds,3)
Seed coordinates in 3D space.
"""
rng = _np.random.default_rng(rng_seed)
coords = _np.empty((N_seeds,3))
coords[0] = rng.random(3) * size
coords[0] = rng.random(3) * _np.array(size,float)
s = 1
i = 0
@ -96,8 +99,8 @@ def from_Poisson_disc(size: _np.ndarray, N_seeds: int, N_candidates: int, distan
return coords
def from_grid(grid, selection: Sequence[int] = None,
invert: bool = False, average: bool = False, periodic: bool = True) -> Tuple[_np.ndarray, _np.ndarray]:
def from_grid(grid, selection: _IntSequence = None,
invert: bool = False, average: bool = False, periodic: bool = True) -> _Tuple[_np.ndarray, _np.ndarray]:
"""
Create seeds from grid description.
@ -105,7 +108,7 @@ def from_grid(grid, selection: Sequence[int] = None,
----------
grid : damask.Grid
Grid from which the material IDs are used as seeds.
selection : iterable of integers, optional
selection : sequence of int, optional
Material IDs to consider.
invert : boolean, false
Consider all material IDs except those in selection. Defaults to False.
@ -116,7 +119,7 @@ def from_grid(grid, selection: Sequence[int] = None,
Returns
-------
coords, materials : numpy.ndarray of shape (:,3), numpy.ndarray of shape (:)
coords, materials : numpy.ndarray, shape (:,3); numpy.ndarray, shape (:)
Seed coordinates in 3D space, material IDs.
"""

66
python/damask/util.py
View File

@ -22,7 +22,7 @@ __all__=[
'natural_sort',
'show_progress',
'scale_to_coprime',
'project_stereographic',
'project_equal_angle', 'project_equal_area',
'hybrid_IA',
'execution_stamp',
'shapeshifter', 'shapeblender',
@ -267,13 +267,13 @@ def scale_to_coprime(v):
return m
def project_stereographic(vector,direction='z',normalize=True,keepdims=False):
def project_equal_angle(vector,direction='z',normalize=True,keepdims=False):
"""
Apply stereographic projection to vector.
Apply equal-angle projection to vector.
Parameters
----------
vector : numpy.ndarray of shape (...,3)
vector : numpy.ndarray, shape (...,3)
Vector coordinates to be projected.
direction : str
Projection direction 'x', 'y', or 'z'.
@ -281,32 +281,74 @@ def project_stereographic(vector,direction='z',normalize=True,keepdims=False):
normalize : bool
Ensure unit length of input vector. Defaults to True.
keepdims : bool
Maintain three-dimensional output coordinates.
Default two-dimensional output uses right-handed frame spanned by
Maintain three-dimensional output coordinates. Defaults to False.
Two-dimensional output uses right-handed frame spanned by
the next and next-next axis relative to the projection direction,
e.g. x-y when projecting along z and z-x when projecting along y.
Returns
-------
coordinates : numpy.ndarray of shape (...,2 | 3)
coordinates : numpy.ndarray, shape (...,2 | 3)
Projected coordinates.
Examples
--------
>>> import damask
>>> import numpy as np
>>> project_stereographic(np.ones(3))
>>> project_equal_angle(np.ones(3))
[0.3660254, 0.3660254]
>>> project_stereographic(np.ones(3),direction='x',normalize=False,keepdims=True)
>>> project_equal_angle(np.ones(3),direction='x',normalize=False,keepdims=True)
[0, 0.5, 0.5]
>>> project_stereographic([0,1,1],direction='y',normalize=True,keepdims=False)
>>> project_equal_angle([0,1,1],direction='y',normalize=True,keepdims=False)
[0.41421356, 0]
"""
shift = 'zyx'.index(direction)
v_ = np.roll(vector/np.linalg.norm(vector,axis=-1,keepdims=True) if normalize else vector,
v = np.roll(vector/np.linalg.norm(vector,axis=-1,keepdims=True) if normalize else vector,
shift,axis=-1)
return np.roll(np.block([v_[...,:2]/(1+np.abs(v_[...,2:3])),np.zeros_like(v_[...,2:3])]),
return np.roll(np.block([v[...,:2]/(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 project_equal_area(vector,direction='z',normalize=True,keepdims=False):
"""
Apply equal-area projection to vector.
Parameters
----------
vector : numpy.ndarray, shape (...,3)
Vector coordinates to be projected.
direction : str
Projection direction 'x', 'y', or 'z'.
Defaults to 'z'.
normalize : bool
Ensure unit length of input vector. Defaults to True.
keepdims : bool
Maintain three-dimensional output coordinates. Defaults to False.
Two-dimensional output uses right-handed frame spanned by
the next and next-next axis relative to the projection direction,
e.g. x-y when projecting along z and z-x when projecting along y.
Returns
-------
coordinates : numpy.ndarray, shape (...,2 | 3)
Projected coordinates.
Examples
--------
>>> import damask
>>> import numpy as np
>>> project_equal_area(np.ones(3))
[0.45970084, 0.45970084]
>>> project_equal_area(np.ones(3),direction='x',normalize=False,keepdims=True)
[0.0, 0.70710678, 0.70710678]
>>> project_equal_area([0,1,1],direction='y',normalize=True,keepdims=False)
[0.5411961, 0.0]
"""
shift = 'zyx'.index(direction)
v = np.roll(vector/np.linalg.norm(vector,axis=-1,keepdims=True) if normalize else vector,
shift,axis=-1)
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]

20
python/tests/test_Crystal.py
View File

@ -79,3 +79,23 @@ class TestCrystal:
a=a,b=b,c=c,
alpha=alpha,beta=beta,gamma=gamma)
assert np.allclose(points,c.lattice_points)
@pytest.mark.parametrize('crystal,length',
[(Crystal(lattice='cF'),[12,6]),
(Crystal(lattice='cI'),[12,12,24]),
(Crystal(lattice='hP'),[3,3,6,12,6]),
(Crystal(lattice='tI',c=1.2),[2,2,2,4,2,4,2,2,4,8,4,8,8])
])
def test_N_slip(self,crystal,length):
assert [len(s) for s in crystal.kinematics('slip')['direction']] == length
assert [len(s) for s in crystal.kinematics('slip')['plane']] == length
@pytest.mark.parametrize('crystal,length',
[(Crystal(lattice='cF'),[12]),
(Crystal(lattice='cI'),[12]),
(Crystal(lattice='hP'),[6,6,6,6]),
])
def test_N_twin(self,crystal,length):
assert [len(s) for s in crystal.kinematics('twin')['direction']] == length
assert [len(s) for s in crystal.kinematics('twin')['plane']] == length

11
python/tests/test_grid_filters.py
View File

@ -2,6 +2,8 @@ import pytest
import numpy as np
from damask import grid_filters
from damask import Grid
from damask import seeds
class TestGridFilters:
@ -139,12 +141,19 @@ class TestGridFilters:
else:
function(unordered,mode)
def test_regrid(self):
def test_regrid_identity(self):
size = np.random.random(3)
cells = np.random.randint(8,32,(3))
F = np.broadcast_to(np.eye(3), tuple(cells)+(3,3))
assert all(grid_filters.regrid(size,F,cells) == np.arange(cells.prod()))
def test_regrid_double_cells(self):
size = np.random.random(3)
cells = np.random.randint(8,32,(3))
g = Grid.from_Voronoi_tessellation(cells,size,seeds.from_random(size,10))
F = np.broadcast_to(np.eye(3), tuple(cells)+(3,3))
assert all(g.scale(cells*2).material.flatten() ==
g.material.flatten()[grid_filters.regrid(size,F,cells*2)])
@pytest.mark.parametrize('differential_operator',[grid_filters.curl,
grid_filters.divergence,

17
python/tests/test_util.py
View File

@ -59,8 +59,21 @@ class TestUtil:
([1,1,0],'x',False,False,[0.5,0]),
([1,1,1],'y',True, True, [0.3660254, 0,0.3660254]),
])
def test_project_stereographic(self,point,direction,normalize,keepdims,answer):
assert np.allclose(util.project_stereographic(np.array(point),direction=direction,
def test_project_equal_angle(self,point,direction,normalize,keepdims,answer):
assert np.allclose(util.project_equal_angle(np.array(point),direction=direction,
normalize=normalize,keepdims=keepdims),answer)
@pytest.mark.parametrize('point,direction,normalize,keepdims,answer',
[
([1,0,0],'z',False,True, [1,0,0]),
([1,0,0],'z',True, False,[1,0]),
([0,1,1],'z',False,True, [0,0.70710678,0]),
([0,1,1],'y',True, False,[0.5411961,0]),
([1,1,0],'x',False,False,[0.70710678,0]),
([1,1,1],'y',True, True, [0.45970084,0,0.45970084]),
])
def test_project_equal_area(self,point,direction,normalize,keepdims,answer):
assert np.allclose(util.project_equal_area(np.array(point),direction=direction,
normalize=normalize,keepdims=keepdims),answer)
@pytest.mark.parametrize('fro,to,mode,answer',

14
src/LAPACK_interface.f90
View File

@ -6,7 +6,7 @@
module LAPACK_interface
interface
subroutine dgeev(jobvl,jobvr,n,a,lda,wr,wi,vl,ldvl,vr,ldvr,work,lwork,info)
pure subroutine dgeev(jobvl,jobvr,n,a,lda,wr,wi,vl,ldvl,vr,ldvr,work,lwork,info)
use prec
character, intent(in) :: jobvl,jobvr
integer, intent(in) :: n,lda,ldvl,ldvr,lwork
@ -18,16 +18,16 @@ module LAPACK_interface
integer, intent(out) :: info
end subroutine dgeev
subroutine dgesv(n,nrhs,a,lda,ipiv,b,ldb,info)
pure subroutine dgesv(n,nrhs,a,lda,ipiv,b,ldb,info)
use prec
integer, intent(in) :: n,nrhs,lda,ldb
real(pReal), intent(inout), dimension(lda,n) :: a
integer, intent(out), dimension(n) :: ipiv
real(pReal), intent(out), dimension(ldb,nrhs) :: b
real(pReal), intent(inout), dimension(ldb,nrhs) :: b
integer, intent(out) :: info
end subroutine dgesv
subroutine dgetrf(m,n,a,lda,ipiv,info)
pure subroutine dgetrf(m,n,a,lda,ipiv,info)
use prec
integer, intent(in) :: m,n,lda
real(pReal), intent(inout), dimension(lda,n) :: a
@ -35,16 +35,16 @@ module LAPACK_interface
integer, intent(out) :: info
end subroutine dgetrf
subroutine dgetri(n,a,lda,ipiv,work,lwork,info)
pure subroutine dgetri(n,a,lda,ipiv,work,lwork,info)
use prec
integer, intent(in) :: n,lda,lwork
real(pReal), intent(inout), dimension(lda,n) :: a
integer, intent(out), dimension(n) :: ipiv
integer, intent(in), dimension(n) :: ipiv
real(pReal), intent(out), dimension(max(1,lwork)) :: work
integer, intent(out) :: info
end subroutine dgetri
subroutine dsyev(jobz,uplo,n,a,lda,w,work,lwork,info)
pure subroutine dsyev(jobz,uplo,n,a,lda,w,work,lwork,info)
use prec
character, intent(in) :: jobz,uplo
integer, intent(in) :: n,lda,lwork

5
src/constants.f90
View File

@ -9,7 +9,8 @@ module constants
public
real(pReal), parameter :: &
T_ROOM = 300.0_pReal, & !< Room temperature in K
K_B = 1.38e-23_pReal !< Boltzmann constant in J/Kelvin
T_ROOM = 300.0_pReal, & !< Room temperature in K. ToDo: IUPAC: 298.15
K_B = 1.38e-23_pReal, & !< Boltzmann constant in J/Kelvin
N_A = 6.02214076e23_pReal !< Avogadro constant in 1/mol
end module constants

4
src/lattice.f90
View File

@ -2070,7 +2070,7 @@ end function getlabels
!> @brief Equivalent Poisson's ratio (ν)
!> @details https://doi.org/10.1143/JPSJ.20.635
!--------------------------------------------------------------------------------------------------
function lattice_equivalent_nu(C,assumption) result(nu)
pure function lattice_equivalent_nu(C,assumption) result(nu)
real(pReal), dimension(6,6), intent(in) :: C !< Stiffness tensor (Voigt notation)
character(len=5), intent(in) :: assumption !< Assumption ('Voigt' = isostrain, 'Reuss' = isostress)
@ -2103,7 +2103,7 @@ end function lattice_equivalent_nu
!> @brief Equivalent shear modulus (μ)
!> @details https://doi.org/10.1143/JPSJ.20.635
!--------------------------------------------------------------------------------------------------
function lattice_equivalent_mu(C,assumption) result(mu)
pure function lattice_equivalent_mu(C,assumption) result(mu)
real(pReal), dimension(6,6), intent(in) :: C !< Stiffness tensor (Voigt notation)
character(len=5), intent(in) :: assumption !< Assumption ('Voigt' = isostrain, 'Reuss' = isostress)

22
src/math.f90
View File

@ -512,7 +512,7 @@ end subroutine math_invert33
!--------------------------------------------------------------------------------------------------
!> @brief Inversion of symmetriced 3x3x3x3 matrix
!--------------------------------------------------------------------------------------------------
function math_invSym3333(A)
pure function math_invSym3333(A)
real(pReal),dimension(3,3,3,3) :: math_invSym3333
@ -538,7 +538,7 @@ end function math_invSym3333
!--------------------------------------------------------------------------------------------------
!> @brief invert quadratic matrix of arbitrary dimension
!--------------------------------------------------------------------------------------------------
subroutine math_invert(InvA, error, A)
pure subroutine math_invert(InvA, error, A)
real(pReal), dimension(:,:), intent(in) :: A
real(pReal), dimension(size(A,1),size(A,1)), intent(out) :: invA
@ -996,7 +996,7 @@ end subroutine math_normal
!--------------------------------------------------------------------------------------------------
!> @brief eigenvalues and eigenvectors of symmetric matrix
!--------------------------------------------------------------------------------------------------
subroutine math_eigh(w,v,error,m)
pure subroutine math_eigh(w,v,error,m)
real(pReal), dimension(:,:), intent(in) :: m !< quadratic matrix to compute eigenvectors and values of
real(pReal), dimension(size(m,1)), intent(out) :: w !< eigenvalues
@ -1021,7 +1021,7 @@ end subroutine math_eigh
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @details See http://arxiv.org/abs/physics/0610206 (DSYEVH3)
!--------------------------------------------------------------------------------------------------
subroutine math_eigh33(w,v,m)
pure subroutine math_eigh33(w,v,m)
real(pReal), dimension(3,3),intent(in) :: m !< 3x3 matrix to compute eigenvectors and values of
real(pReal), dimension(3), intent(out) :: w !< eigenvalues
@ -1114,7 +1114,7 @@ end function math_rotationalPart
!> @brief Eigenvalues of symmetric matrix
! will return NaN on error
!--------------------------------------------------------------------------------------------------
function math_eigvalsh(m)
pure function math_eigvalsh(m)
real(pReal), dimension(:,:), intent(in) :: m !< symmetric matrix to compute eigenvalues of
real(pReal), dimension(size(m,1)) :: math_eigvalsh
@ -1137,7 +1137,7 @@ end function math_eigvalsh
!> but apparently more stable solution and has general LAPACK powered version for arbritrary sized
!> matrices as fallback
!--------------------------------------------------------------------------------------------------
function math_eigvalsh33(m)
pure function math_eigvalsh33(m)
real(pReal), intent(in), dimension(3,3) :: m !< 3x3 symmetric matrix to compute eigenvalues of
real(pReal), dimension(3) :: math_eigvalsh33,I
@ -1432,9 +1432,11 @@ subroutine selfTest
error stop 'math_LeviCivita'
normal_distribution: block
real(pReal), dimension(500000) :: r
integer, parameter :: N = 1000000
real(pReal), dimension(:), allocatable :: r
real(pReal) :: mu, sigma
allocate(r(N))
call random_number(mu)
call random_number(sigma)
@ -1443,11 +1445,11 @@ subroutine selfTest
call math_normal(r,mu,sigma)
if (abs(mu -sum(r)/real(size(r),pReal))>5.0e-2_pReal) &
if (abs(mu -sum(r)/real(N,pReal))>5.0e-2_pReal) &
error stop 'math_normal(mu)'
mu = sum(r)/real(size(r),pReal)
if (abs(sigma**2 -1.0_pReal/real(size(r)-1,pReal) * sum((r-mu)**2))/sigma > 5.0e-2_pReal) &
mu = sum(r)/real(N,pReal)
if (abs(sigma**2 -1.0_pReal/real(N-1,pReal) * sum((r-mu)**2))/sigma > 5.0e-2_pReal) &
error stop 'math_normal(sigma)'
end block normal_distribution

6
src/phase_mechanical.f90
View File

@ -168,17 +168,17 @@ submodule(phase) mechanical
integer, intent(in) :: ph,en
end function plastic_dislotwin_homogenizedC
module function elastic_C66(ph,en) result(C66)
pure module function elastic_C66(ph,en) result(C66)
real(pReal), dimension(6,6) :: C66
integer, intent(in) :: ph, en
end function elastic_C66
module function elastic_mu(ph,en) result(mu)
pure module function elastic_mu(ph,en) result(mu)
real(pReal) :: mu
integer, intent(in) :: ph, en
end function elastic_mu
module function elastic_nu(ph,en) result(nu)
pure module function elastic_nu(ph,en) result(nu)
real(pReal) :: nu
integer, intent(in) :: ph, en
end function elastic_nu

8
src/phase_mechanical_elastic.f90
View File

@ -30,7 +30,7 @@ module subroutine elastic_init(phases)
phase, &
mech, &
elastic
logical :: thermal_active
print'(/,1x,a)', '<<<+- phase:mechanical:elastic init -+>>>'
print'(/,1x,a)', '<<<+- phase:mechanical:elastic:Hooke init -+>>>'
@ -86,7 +86,7 @@ end subroutine elastic_init
!--------------------------------------------------------------------------------------------------
!> @brief return 6x6 elasticity tensor
!--------------------------------------------------------------------------------------------------
module function elastic_C66(ph,en) result(C66)
pure module function elastic_C66(ph,en) result(C66)
integer, intent(in) :: &
ph, &
@ -140,7 +140,7 @@ end function elastic_C66
!--------------------------------------------------------------------------------------------------
!> @brief return shear modulus
!--------------------------------------------------------------------------------------------------
module function elastic_mu(ph,en) result(mu)
pure module function elastic_mu(ph,en) result(mu)
integer, intent(in) :: &
ph, &
@ -157,7 +157,7 @@ end function elastic_mu
!--------------------------------------------------------------------------------------------------
!> @brief return Poisson ratio
!--------------------------------------------------------------------------------------------------
module function elastic_nu(ph,en) result(nu)
pure module function elastic_nu(ph,en) result(nu)
integer, intent(in) :: &
ph, &

2
src/phase_mechanical_plastic_kinehardening.f90
View File

@ -86,6 +86,8 @@ module function plastic_kinehardening_init() result(myPlasticity)
print'(/,1x,a)', '<<<+- phase:mechanical:plastic:kinehardening init -+>>>'
print'(/,a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT)
print'(/,1x,a)', 'J.A. Wollmershauser et al., International Journal of Fatigue 36:181–193, 2012'
print'( 1x,a)', 'https://doi.org/10.1016/j.ijfatigue.2011.07.008'
phases => config_material%get('phase')
allocate(param(phases%length))

8
src/rotations.f90
View File

@ -372,7 +372,7 @@ end function rotTensor4
!---------------------------------------------------------------------------------------------------
!> @brief Rotate a rank-4 tensor in Voigt 6x6 notation passively (default) or actively.
!> @brief Rotate a rank-4 stiffness tensor in Voigt 6x6 notation passively (default) or actively.
!> @details: https://scicomp.stackexchange.com/questions/35600
!! ToDo: Need to check active/passive !!!
!---------------------------------------------------------------------------------------------------
@ -393,11 +393,11 @@ pure function rotStiffness(self,C,active) result(cRot)
R = self%asMatrix()
endif
M = reshape([R(1,1)**2.0_pReal, R(2,1)**2.0_pReal, R(3,1)**2.0_pReal, &
M = reshape([R(1,1)**2, R(2,1)**2, R(3,1)**2, &
R(2,1)*R(3,1), R(1,1)*R(3,1), R(1,1)*R(2,1), &
R(1,2)**2.0_pReal, R(2,2)**2.0_pReal, R(3,2)**2.0_pReal, &
R(1,2)**2, R(2,2)**2, R(3,2)**2, &
R(2,2)*R(3,2), R(1,2)*R(3,2), R(1,2)*R(2,2), &
R(1,3)**2.0_pReal, R(2,3)**2.0_pReal, R(3,3)**2.0_pReal, &
R(1,3)**2, R(2,3)**2, R(3,3)**2, &
R(2,3)*R(3,3), R(1,3)*R(3,3), R(1,3)*R(2,3), &
2.0_pReal*R(1,2)*R(1,3), 2.0_pReal*R(2,2)*R(2,3), 2.0_pReal*R(3,2)*R(3,3), &
R(2,2)*R(3,3)+R(2,3)*R(3,2), R(1,2)*R(3,3)+R(1,3)*R(3,2), R(1,2)*R(2,3)+R(1,3)*R(2,2), &