From baf04ad5ebc3b9ccb7e9155e753312b6abe5d0c4 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 28 Sep 2023 11:52:34 +0200 Subject: [PATCH 01/13] TMPDIR controls the behavior of mktemp --- .gitlab-ci.yml | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 0a2439020..5ddfcc27e 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -89,10 +89,10 @@ unittest_GNU_DEBUG: stage: compile script: - module load ${COMPILER_GNU} ${MPI_GNU} ${PETSC_GNU} - - TMPDIR=$(mktemp -d) - - cmake -B ${TMPDIR} -DDAMASK_SOLVER=test -DCMAKE_INSTALL_PREFIX=${TMPDIR} -DCMAKE_BUILD_TYPE=RELEASE -DBUILDCMD_POST=-coverage - - cmake --build ${TMPDIR} --target install - - cd ${TMPDIR} + - TEMPDIR=$(mktemp -d) + - cmake -B ${TEMPDIR} -DDAMASK_SOLVER=test -DCMAKE_INSTALL_PREFIX=${TEMPDIR} -DCMAKE_BUILD_TYPE=RELEASE -DBUILDCMD_POST=-coverage + - cmake --build ${TEMPDIR} --target install + - cd ${TEMPDIR} - ./bin/DAMASK_test - find -name \*.gcda -not -path "**/test/*" | xargs gcov @@ -100,10 +100,10 @@ unittest_GNU_RELEASE: stage: compile script: - module load ${COMPILER_GNU} ${MPI_GNU} ${PETSC_GNU} - - TMPDIR=$(mktemp -d) - - cmake -B ${TMPDIR} -DDAMASK_SOLVER=test -DCMAKE_INSTALL_PREFIX=${TMPDIR} -DCMAKE_BUILD_TYPE=RELEASE -DBUILDCMD_POST=-coverage - - cmake --build ${TMPDIR} --target install - - cd ${TMPDIR} + - TEMPDIR=$(mktemp -d) + - cmake -B ${TEMPDIR} -DDAMASK_SOLVER=test -DCMAKE_INSTALL_PREFIX=${TEMPDIR} -DCMAKE_BUILD_TYPE=RELEASE -DBUILDCMD_POST=-coverage + - cmake --build ${TEMPDIR} --target install + - cd ${TEMPDIR} - ./bin/DAMASK_test - find -name \*.gcda -not -path "**/test/*" | xargs gcov @@ -111,10 +111,10 @@ unittest_GNU_PERFORMANCE: stage: compile script: - module load ${COMPILER_GNU} ${MPI_GNU} ${PETSC_GNU} - - TMPDIR=$(mktemp -d) - - cmake -B ${TMPDIR} -DDAMASK_SOLVER=test -DCMAKE_INSTALL_PREFIX=${TMPDIR} -DCMAKE_BUILD_TYPE=PERFORMANCE -DBUILDCMD_POST=-coverage - - cmake --build ${TMPDIR} --target install - - cd ${TMPDIR} + - TEMPDIR=$(mktemp -d) + - cmake -B ${TEMPDIR} -DDAMASK_SOLVER=test -DCMAKE_INSTALL_PREFIX=${TEMPDIR} -DCMAKE_BUILD_TYPE=PERFORMANCE -DBUILDCMD_POST=-coverage + - cmake --build ${TEMPDIR} --target install + - cd ${TEMPDIR} - ./bin/DAMASK_test - find -name \*.gcda -not -path "**/test/*" | xargs gcov From 05e675bbdb38de09f3337cf24279d0d6b53ec680 Mon Sep 17 00:00:00 2001 From: Test User Date: Fri, 29 Sep 2023 18:53:37 +0200 Subject: [PATCH 02/13] [skip ci] updated version information after successful test of v3.0.0-alpha7-867-g026dd3c14 --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index cd37d7af8..baa2c1a48 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -3.0.0-alpha7-864-g9cf37c493 +3.0.0-alpha7-867-g026dd3c14 From a0dc25c16ea89878981093134504a1e10353cac1 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 25 Sep 2023 15:12:00 +0200 Subject: [PATCH 03/13] more specific typehints --- python/damask/_result.py | 70 ++++++++++++++++++------------------- python/damask/_typehints.py | 14 +++++++- 2 files changed, 48 insertions(+), 36 deletions(-) diff --git a/python/damask/_result.py b/python/damask/_result.py index f36e1e83a..5f55fd13a 100644 --- a/python/damask/_result.py +++ b/python/damask/_result.py @@ -22,7 +22,7 @@ from . import grid_filters from . import mechanics from . import tensor from . import util -from ._typehints import FloatSequence, IntSequence +from ._typehints import FloatSequence, IntSequence, DADF5Dataset h5py3 = h5py.__version__[0] == '3' @@ -36,7 +36,7 @@ def _read(dataset: h5py._hl.dataset.Dataset) -> np.ndarray: return np.array(dataset,dtype=dtype) def _match(requested, - existing: h5py._hl.base.KeysViewHDF5) -> List[Any]: + existing: h5py._hl.base.KeysViewHDF5) -> List[str]: """Find matches among two sets of labels.""" def flatten_list(list_of_lists): return [e for e_ in list_of_lists for e in e_] @@ -609,7 +609,7 @@ class Result: Name of scalar, vector, or tensor dataset to take absolute value of. """ - def absolute(x: Dict[str, Any]) -> Dict[str, Any]: + def absolute(x: DADF5Dataset) -> DADF5Dataset: return { 'data': np.abs(x['data']), 'label': f'|{x["label"]}|', @@ -671,7 +671,7 @@ class Result: ... 'Mises equivalent of the Cauchy stress') """ - def calculation(**kwargs) -> Dict[str, Any]: + def calculation(**kwargs) -> DADF5Dataset: formula = kwargs['formula'] for d in re.findall(r'#(.*?)#',formula): formula = formula.replace(f'#{d}#',f"kwargs['{d}']['data']") @@ -712,7 +712,7 @@ class Result: """ - def stress_Cauchy(P: Dict[str, Any], F: Dict[str, Any]) -> Dict[str, Any]: + def stress_Cauchy(P: DADF5Dataset, F: DADF5Dataset) -> DADF5Dataset: return { 'data': mechanics.stress_Cauchy(P['data'],F['data']), 'label': 'sigma', @@ -747,7 +747,7 @@ class Result: """ - def determinant(T: Dict[str, Any]) -> Dict[str, Any]: + def determinant(T: DADF5Dataset) -> DADF5Dataset: return { 'data': np.linalg.det(T['data']), 'label': f"det({T['label']})", @@ -780,7 +780,7 @@ class Result: """ - def deviator(T: Dict[str, Any]) -> Dict[str, Any]: + def deviator(T: DADF5Dataset) -> DADF5Dataset: return { 'data': tensor.deviatoric(T['data']), 'label': f"s_{T['label']}", @@ -817,7 +817,7 @@ class Result: """ - def eigenval(T_sym: Dict[str, Any], eigenvalue: Literal['max, mid, min']) -> Dict[str, Any]: + def eigenval(T_sym: DADF5Dataset, eigenvalue: Literal['max, mid, min']) -> DADF5Dataset: if eigenvalue == 'max': label,p = 'maximum',2 elif eigenvalue == 'mid': @@ -856,7 +856,7 @@ class Result: """ - def eigenvector(T_sym: Dict[str, Any], eigenvalue: Literal['max', 'mid', 'min']) -> Dict[str, Any]: + def eigenvector(T_sym: DADF5Dataset, eigenvalue: Literal['max', 'mid', 'min']) -> DADF5Dataset: if eigenvalue == 'max': label,p = 'maximum',2 elif eigenvalue == 'mid': @@ -904,13 +904,13 @@ class Result: """ - def IPF_color(l: FloatSequence, q: Dict[str, Any]) -> Dict[str, Any]: + def IPF_color(l: FloatSequence, q: DADF5Dataset) -> DADF5Dataset: m = util.scale_to_coprime(np.array(l)) lattice = q['meta']['lattice'] o = Orientation(rotation = q['data'],lattice=lattice) return { - 'data': np.uint8(o.IPF_color(l)*255), + 'data': (o.IPF_color(l)*255).astype(np.uint8), 'label': 'IPFcolor_({} {} {})'.format(*m), 'meta' : { 'unit': '8-bit RGB', @@ -933,7 +933,7 @@ class Result: Name of symmetric tensor dataset. """ - def maximum_shear(T_sym: Dict[str, Any]) -> Dict[str, Any]: + def maximum_shear(T_sym: DADF5Dataset) -> DADF5Dataset: return { 'data': mechanics.maximum_shear(T_sym['data']), 'label': f"max_shear({T_sym['label']})", @@ -976,7 +976,7 @@ class Result: >>> r.add_equivalent_Mises('epsilon_V^0.0(F)') """ - def equivalent_Mises(T_sym: Dict[str, Any], kind: str) -> Dict[str, Any]: + def equivalent_Mises(T_sym: DADF5Dataset, kind: str) -> DADF5Dataset: k = kind if k is None: if T_sym['meta']['unit'] == '1': @@ -1014,7 +1014,7 @@ class Result: Order of the norm. inf means NumPy's inf object. For details refer to numpy.linalg.norm. """ - def norm(x: Dict[str, Any], ord: Union[int, float, Literal['fro', 'nuc']]) -> Dict[str, Any]: + def norm(x: DADF5Dataset, ord: Union[int, float, Literal['fro', 'nuc']]) -> DADF5Dataset: o = ord if len(x['data'].shape) == 2: axis: Union[int, Tuple[int, int]] = 1 @@ -1062,7 +1062,7 @@ class Result: is taken into account. """ - def stress_second_Piola_Kirchhoff(P: Dict[str, Any], F: Dict[str, Any]) -> Dict[str, Any]: + def stress_second_Piola_Kirchhoff(P: DADF5Dataset, F: DADF5Dataset) -> DADF5Dataset: return { 'data': mechanics.stress_second_Piola_Kirchhoff(P['data'],F['data']), 'label': 'S', @@ -1104,12 +1104,11 @@ class Result: Defaults to True. """ - def pole(q: Dict[str, Any], - uvw: FloatSequence, - hkl: FloatSequence, - with_symmetry: bool, - normalize: bool) -> Dict[str, Any]: - c = q['meta']['c/a'] if 'c/a' in q['meta'] else 1 + def pole(q: DADF5Dataset, + uvw: FloatSequence, hkl: FloatSequence, + with_symmetry: bool, + normalize: bool) -> DADF5Dataset: + c = q['meta']['c/a'] if 'c/a' in q['meta'] else 1.0 brackets = ['[]','()','⟨⟩','{}'][(uvw is None)*1+with_symmetry*2] label = 'p^' + '{}{} {} {}{}'.format(brackets[0], *(uvw if uvw else hkl), @@ -1149,7 +1148,7 @@ class Result: >>> r.add_rotation('F') """ - def rotation(F: Dict[str, Any]) -> Dict[str, Any]: + def rotation(F: DADF5Dataset) -> DADF5Dataset: return { 'data': mechanics.rotation(F['data']).as_matrix(), 'label': f"R({F['label']})", @@ -1181,7 +1180,7 @@ class Result: >>> r.add_spherical('sigma') """ - def spherical(T: Dict[str, Any]) -> Dict[str, Any]: + def spherical(T: DADF5Dataset) -> DADF5Dataset: return { 'data': tensor.spherical(T['data'],False), 'label': f"p_{T['label']}", @@ -1255,7 +1254,7 @@ class Result: | https://de.wikipedia.org/wiki/Verzerrungstensor """ - def strain(F: Dict[str, Any], t: Literal['V', 'U'], m: float) -> Dict[str, Any]: + def strain(F: DADF5Dataset, t: Literal['V', 'U'], m: float) -> DADF5Dataset: side = 'left' if t == 'V' else 'right' return { 'data': mechanics.strain(F['data'],t,m), @@ -1286,7 +1285,7 @@ class Result: Defaults to 'V'. """ - def stretch_tensor(F: Dict[str, Any], t: str) -> Dict[str, Any]: + def stretch_tensor(F: DADF5Dataset, t: str) -> DADF5Dataset: return { 'data': (mechanics.stretch_left if t.upper() == 'V' else mechanics.stretch_right)(F['data']), 'label': f"{t}({F['label']})", @@ -1316,7 +1315,7 @@ class Result: i.e. fields resulting from the grid solver. """ - def curl(f: Dict[str, Any], size: np.ndarray) -> Dict[str, Any]: + def curl(f: DADF5Dataset, size: np.ndarray) -> DADF5Dataset: return { 'data': grid_filters.curl(size,f['data']), 'label': f"curl({f['label']})", @@ -1345,7 +1344,7 @@ class Result: i.e. fields resulting from the grid solver. """ - def divergence(f: Dict[str, Any], size: np.ndarray) -> Dict[str, Any]: + def divergence(f: DADF5Dataset, size: np.ndarray) -> DADF5Dataset: return { 'data': grid_filters.divergence(size,f['data']), 'label': f"divergence({f['label']})", @@ -1374,7 +1373,7 @@ class Result: i.e. fields resulting from the grid solver. """ - def gradient(f: Dict[str, Any], size: np.ndarray) -> Dict[str, Any]: + def gradient(f: DADF5Dataset, size: np.ndarray) -> DADF5Dataset: return { 'data': grid_filters.gradient(size,f['data'] if len(f['data'].shape) == 4 else \ f['data'].reshape(f['data'].shape+(1,))), @@ -1390,7 +1389,7 @@ class Result: def _add_generic_grid(self, - func: Callable, + func: Callable[..., DADF5Dataset], datasets: Dict[str, str], args: Dict[str, str] = {}, constituents = None): @@ -1441,7 +1440,7 @@ class Result: now.strftime('%Y-%m-%d %H:%M:%S%z').encode() for l,v in r['meta'].items(): - h5_dataset.attrs[l.lower()]=v if h5py3 else v.encode() + h5_dataset.attrs[l.lower()]=v.encode() if not h5py3 and type(v) is str else v creator = h5_dataset.attrs['creator'] if h5py3 else \ h5_dataset.attrs['creator'].decode() h5_dataset.attrs['creator'] = f'damask.Result.{creator} v{damask.version}' if h5py3 else \ @@ -1451,8 +1450,8 @@ class Result: def _add_generic_pointwise(self, - func: Callable, - datasets: Dict[str, Any], + func: Callable[..., DADF5Dataset], + datasets: Dict[str, str], args: Dict[str, Any] = {}): """ General function to add pointwise data. @@ -1471,9 +1470,9 @@ class Result: """ def job_pointwise(group: str, - callback: Callable, + callback: Callable[..., DADF5Dataset], datasets: Dict[str, str], - args: Dict[str, str]) -> Union[None, Any]: + args: Dict[str, str]) -> Union[None, DADF5Dataset]: try: datasets_in = {} with h5py.File(self.fname,'r') as f: @@ -1561,7 +1560,7 @@ class Result: def get(self, output: Union[str, List[str]] = '*', flatten: bool = True, - prune: bool = True) -> Optional[Dict[str,Any]]: + prune: bool = True) -> Union[None,Dict[str,Any]]: """ Collect data per phase/homogenization reflecting the group/folder structure in the DADF5 file. @@ -1930,6 +1929,7 @@ class Result: v.save(vtk_dir/f'{self.fname.stem}_inc{inc.split(prefix_inc)[-1].zfill(N_digits)}', parallel=parallel) + def export_DADF5(self, fname, output: Union[str, List[str]] = '*', diff --git a/python/damask/_typehints.py b/python/damask/_typehints.py index 29ab9f3ad..5dbd80cae 100644 --- a/python/damask/_typehints.py +++ b/python/damask/_typehints.py @@ -1,6 +1,6 @@ """Functionality for typehints.""" -from typing import Sequence, Union, Literal, TextIO +from typing import Sequence, Union, TypedDict, Literal, TextIO from pathlib import Path import numpy as np @@ -16,3 +16,15 @@ CrystalKinematics = Literal['slip', 'twin'] NumpyRngSeed = Union[int, IntSequence, np.random.SeedSequence, np.random.Generator] # BitGenerator does not exists in older numpy versions #NumpyRngSeed = Union[int, IntSequence, np.random.SeedSequence, np.random.BitGenerator, np.random.Generator] + +# https://peps.python.org/pep-0655/ +# Metadata = TypedDict('Metadata', {'unit': str, 'description': str, 'creator': str, 'lattice': NotRequired[str]}) +_Metadata = TypedDict('_Metadata', {'lattice': str, 'c/a': float}, total=False) + +class Metadata(_Metadata): + unit: str + description: str + creator: str + + +DADF5Dataset = TypedDict('DADF5Dataset', {'data': np.ndarray, 'label': str, 'meta': Metadata}) From 983a376f45d8d37e5486ca398c234757eb599a57 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Wed, 4 Oct 2023 09:10:11 -0400 Subject: [PATCH 04/13] bugfix; serial string use --- python/damask/_result.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/python/damask/_result.py b/python/damask/_result.py index 5f55fd13a..6893af1bc 100644 --- a/python/damask/_result.py +++ b/python/damask/_result.py @@ -817,7 +817,7 @@ class Result: """ - def eigenval(T_sym: DADF5Dataset, eigenvalue: Literal['max, mid, min']) -> DADF5Dataset: + def eigenval(T_sym: DADF5Dataset, eigenvalue: Literal['max', 'mid', 'min']) -> DADF5Dataset: if eigenvalue == 'max': label,p = 'maximum',2 elif eigenvalue == 'mid': @@ -1261,7 +1261,7 @@ class Result: 'label': f"epsilon_{t}^{m}({F['label']})", 'meta': { 'unit': F['meta']['unit'], - 'description': f'Seth-Hill strain tensor of order {m} based on {side} stretch tensor '+\ + 'description': f'Seth-Hill strain tensor of order {m} based on {side} stretch tensor ' f"of {F['label']} ({F['meta']['description']})", 'creator': 'add_strain' } @@ -1291,8 +1291,8 @@ class Result: 'label': f"{t}({F['label']})", 'meta': { 'unit': F['meta']['unit'], - 'description': f"{'left' if t.upper() == 'V' else 'right'} stretch tensor "\ - +f"of {F['label']} ({F['meta']['description']})", # noqa + 'description': f"{'left' if t.upper() == 'V' else 'right'} stretch tensor " + f"of {F['label']} ({F['meta']['description']})", # noqa 'creator': 'add_stretch_tensor' } } From 550c757cdc87d913f7be9747a14d0efab1321e85 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Wed, 4 Oct 2023 09:34:23 -0400 Subject: [PATCH 05/13] overall smart blending; more robust disFZ --- python/damask/_orientation.py | 84 +++++++++++++++++--------------- python/damask/_rotation.py | 19 ++++++-- python/damask/util.py | 24 +++++---- python/tests/test_Orientation.py | 4 +- python/tests/test_util.py | 4 +- 5 files changed, 80 insertions(+), 55 deletions(-) diff --git a/python/damask/_orientation.py b/python/damask/_orientation.py index 15f8464fd..ef71366c8 100644 --- a/python/damask/_orientation.py +++ b/python/damask/_orientation.py @@ -434,28 +434,29 @@ class Orientation(Rotation,Crystal): https://doi.org/10.1107/S0108767391006864 """ - rho = self.as_Rodrigues_vector(compact=True)*(1.0-1.0e-9) - with np.errstate(invalid='ignore'): + rho = self.as_Rodrigues_vector(compact=True) + rho_ = rho + np.linalg.norm(rho,axis=-1,keepdims=True)*1e-9 + if self.family == 'cubic': - return ((rho[...,0] >= rho[...,1]) & - (rho[...,1] >= rho[...,2]) & - (rho[...,2] >= 0)).astype(bool) + return ((rho_[...,0] >= rho[...,1]) & + (rho_[...,1] >= rho[...,2]) & + (rho_[...,2] >= 0)).astype(bool) if self.family == 'hexagonal': - return ((rho[...,0] >= rho[...,1]*np.sqrt(3)) & - (rho[...,1] >= 0) & - (rho[...,2] >= 0)).astype(bool) + return ((rho_[...,0] >= rho[...,1]*np.sqrt(3)) & + (rho_[...,1] >= 0) & + (rho_[...,2] >= 0)).astype(bool) if self.family == 'tetragonal': - return ((rho[...,0] >= rho[...,1]) & - (rho[...,1] >= 0) & - (rho[...,2] >= 0)).astype(bool) + return ((rho_[...,0] >= rho[...,1]) & + (rho_[...,1] >= 0) & + (rho_[...,2] >= 0)).astype(bool) if self.family == 'orthorhombic': - return ((rho[...,0] >= 0) & - (rho[...,1] >= 0) & - (rho[...,2] >= 0)).astype(bool) + return ((rho_[...,0] >= 0) & + (rho_[...,1] >= 0) & + (rho_[...,2] >= 0)).astype(bool) if self.family == 'monoclinic': - return ((rho[...,1] >= 0) & - (rho[...,2] >= 0)).astype(bool) + return ((rho_[...,1] >= 0) & + (rho_[...,2] >= 0)).astype(bool) return np.ones_like(rho[...,0],dtype=bool) @@ -521,20 +522,21 @@ class Orientation(Rotation,Crystal): if self.family != other.family: raise NotImplementedError('disorientation between different crystal families') - blend = util.shapeblender(self.shape,other.shape) - s = self.equivalent - o = other.equivalent + blend = util.shapeblender( self.shape,other.shape) + s_m = util.shapeshifter( self.shape,blend,mode='right') + s_o = util.shapeshifter(other.shape,blend,mode='left') - s_ = s.reshape((s.shape[0],1)+ self.shape).broadcast_to((s.shape[0],o.shape[0])+blend,mode='right') - o_ = o.reshape((1,o.shape[0])+other.shape).broadcast_to((s.shape[0],o.shape[0])+blend,mode='right') - r_ = s_.misorientation(o_) + s = self.broadcast_to(s_m).equivalent + o = other.broadcast_to(s_o).equivalent + + r_ = s[:,np.newaxis,...].misorientation(o[np.newaxis,:,...]) _r = ~r_ - forward = r_.in_FZ & r_.in_disorientation_FZ - reverse = _r.in_FZ & _r.in_disorientation_FZ - ok = forward | reverse - ok &= (np.cumsum(ok.reshape((-1,)+ok.shape[2:]),axis=0) == 1).reshape(ok.shape) - r = np.where(np.any(forward[...,np.newaxis],axis=(0,1),keepdims=True), + forward = r_.in_disorientation_FZ + reverse = _r.in_disorientation_FZ + ok = (forward | reverse) & r_.in_FZ + ok &= (np.cumsum(ok.reshape((-1,)+blend),axis=0) == 1).reshape(ok.shape) + r = np.where(np.any((ok&forward)[...,np.newaxis],axis=(0,1),keepdims=True), r_.quaternion, _r.quaternion) loc = np.where(ok) @@ -584,6 +586,7 @@ class Orientation(Rotation,Crystal): np.argmin(m,axis=0)[np.newaxis,...,np.newaxis], axis=0), axis=0)) + return ((self.copy(Rotation(r).average(weights)),self.copy(Rotation(r))) if return_cloud else self.copy(Rotation(r).average(weights)) ) @@ -620,17 +623,19 @@ class Orientation(Rotation,Crystal): vector_ = np.array(vector,float) if vector_.shape[-1] != 3: raise ValueError('input is not a field of three-dimensional vectors') - eq = self.equivalent - blend = util.shapeblender(eq.shape,vector_.shape[:-1]) - poles = eq.broadcast_to(blend,mode='right') @ np.broadcast_to(vector_,blend+(3,)) + + blend = util.shapeblender( self.shape,vector_.shape[:-1]) + eq = self.broadcast_to(util.shapeshifter( self.shape,blend,mode='right')).equivalent + poles = np.atleast_2d(eq @ np.broadcast_to(vector_,(1,)+blend+(3,))) ok = self.in_SST(poles,proper=proper) ok &= np.cumsum(ok,axis=0) == 1 loc = np.where(ok) sort = 0 if len(loc) == 1 else np.lexsort(loc[:0:-1]) + return ( - (poles[ok][sort].reshape(blend[1:]+(3,)), (np.vstack(loc[:1]).T)[sort].reshape(blend[1:])) + (poles[ok][sort].reshape(blend+(3,)), (np.vstack(loc[:1]).T)[sort].reshape(blend)) if return_operators else - poles[ok][sort].reshape(blend[1:]+(3,)) + poles[ok][sort].reshape(blend+(3,)) ) @@ -795,16 +800,17 @@ class Orientation(Rotation,Crystal): """ v = self.to_frame(uvw=uvw,hkl=hkl) - blend = util.shapeblender(self.shape,v.shape[:-1]) + s_v = v.shape[:-1] + blend = util.shapeblender(self.shape,s_v) if normalize: - v /= np.linalg.norm(v,axis=-1,keepdims=len(v.shape)>1) + v /= np.linalg.norm(v,axis=-1,keepdims=len(s_v)>0) if with_symmetry: sym_ops = self.symmetry_operations - shape = v.shape[:-1]+sym_ops.shape + s_v += sym_ops.shape blend += sym_ops.shape - v = sym_ops.broadcast_to(shape) \ - @ np.broadcast_to(v.reshape(util.shapeshifter(v.shape,shape+(3,))),shape+(3,)) - return ~(self.broadcast_to(blend))@np.broadcast_to(v,blend+(3,)) + v = sym_ops.broadcast_to(s_v) @ v[...,np.newaxis,:] + + return ~(self.broadcast_to(blend)) @ np.broadcast_to(v,blend+(3,)) def Schmid(self, *, @@ -854,6 +860,7 @@ class Orientation(Rotation,Crystal): p/np.linalg.norm(p,axis=1,keepdims=True)) shape = P.shape[0:1]+self.shape+(3,3) + return ~self.broadcast_to(shape[:-2]) \ @ np.broadcast_to(P.reshape(util.shapeshifter(P.shape,shape)),shape) @@ -897,6 +904,7 @@ class Orientation(Rotation,Crystal): lattice,o = self.relation_operations(model) target = Crystal(lattice=lattice) o = o.broadcast_to(o.shape+self.shape,mode='right') + return Orientation(rotation=o*Rotation(self.quaternion).broadcast_to(o.shape,mode='left'), lattice=lattice, b = self.b if target.ratio['b'] is None else self.a*target.ratio['b'], diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index c03bff4b2..21197f105 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -295,6 +295,7 @@ class Rotation: ---------- other : Rotation, shape (self.shape) Rotation for composition. + Compatible innermost dimensions will blend. Returns ------- @@ -303,10 +304,15 @@ class Rotation: """ if isinstance(other,Rotation): - q_m = self.quaternion[...,0:1] - p_m = self.quaternion[...,1:] - q_o = other.quaternion[...,0:1] - p_o = other.quaternion[...,1:] + blend = util.shapeblender( self.shape,other.shape) + s_m = util.shapeshifter( self.shape,blend,mode='right') + s_o = util.shapeshifter(other.shape,blend,mode='left') + + q_m = self.broadcast_to(s_m).quaternion[...,0:1] + p_m = self.broadcast_to(s_m).quaternion[...,1:] + q_o = other.broadcast_to(s_o).quaternion[...,0:1] + p_o = other.broadcast_to(s_o).quaternion[...,1:] + qmo = q_m*q_o q = (qmo - np.einsum('...i,...i',p_m,p_o).reshape(qmo.shape)) p = q_m*p_o + q_o*p_m + _P * np.cross(p_m,p_o) @@ -325,6 +331,7 @@ class Rotation: ---------- other : Rotation, shape (self.shape) Rotation for composition. + Compatible innermost dimensions will blend. """ return self*other @@ -341,6 +348,7 @@ class Rotation: ---------- other : damask.Rotation, shape (self.shape) Rotation to invert for composition. + Compatible innermost dimensions will blend. Returns ------- @@ -434,9 +442,10 @@ class Rotation: """ if isinstance(other, np.ndarray): - obs = util.shapeblender(self.shape,other.shape,keep_ones=False)[len(self.shape):] + obs = util.shapeblender(self.shape,other.shape)[len(self.shape):] for l in [4,2,1]: if obs[-l:] == l*(3,): + print(f'rotate {l}') bs = util.shapeblender(self.shape,other.shape[:-l],False) self_ = self.broadcast_to(bs) if self.shape != bs else self if l==1: diff --git a/python/damask/util.py b/python/damask/util.py index 242c259b3..513cab87f 100644 --- a/python/damask/util.py +++ b/python/damask/util.py @@ -513,7 +513,7 @@ def shapeshifter(fro: _Tuple[int, ...], def shapeblender(a: _Tuple[int, ...], b: _Tuple[int, ...], - keep_ones: bool = True) -> _Tuple[int, ...]: + keep_ones: bool = False) -> _Tuple[int, ...]: """ Return a shape that overlaps the rightmost entries of 'a' with the leftmost of 'b'. @@ -525,20 +525,24 @@ def shapeblender(a: _Tuple[int, ...], Shape of second array. keep_ones : bool, optional Treat innermost '1's as literal value instead of dimensional placeholder. - Defaults to True. + Defaults to False. Examples -------- - >>> shapeblender((4,4,3),(3,2,1)) - (4,4,3,2,1) - >>> shapeblender((1,2),(1,2,3)) - (1,2,3) - >>> shapeblender((1,),(2,2,1)) - (1,2,2,1) - >>> shapeblender((1,),(2,2,1),False) - (2,2,1) >>> shapeblender((3,2),(3,2)) (3,2) + >>> shapeblender((4,3),(3,2)) + (4,3,2) + >>> shapeblender((4,4),(3,2)) + (4,4,3,2) + >>> shapeblender((1,2),(1,2,3)) + (1,2,3) + >>> shapeblender((),(2,2,1)) + (2,2,1) + >>> shapeblender((1,),(2,2,1)) + (2,2,1) + >>> shapeblender((1,),(2,2,1),True) + (1,2,2,1) """ def is_broadcastable(a,b): diff --git a/python/tests/test_Orientation.py b/python/tests/test_Orientation.py index 9836ee16e..6de2c99b1 100644 --- a/python/tests/test_Orientation.py +++ b/python/tests/test_Orientation.py @@ -358,7 +358,9 @@ class TestOrientation: a=a,b=b,c=c, alpha=alpha,beta=beta,gamma=gamma) assert o.to_pole(**{kw:vector,'with_symmetry':with_symmetry}).shape \ - == o.shape + vector.shape[:-1] + (o.symmetry_operations.shape if with_symmetry else ()) + vector.shape[-1:] + == util.shapeblender(o.shape,vector.shape[:-1]) \ + + (o.symmetry_operations.shape if with_symmetry else ()) \ + + vector.shape[-1:] @pytest.mark.parametrize('lattice',['hP','cI','cF']) #tI not included yet def test_Schmid(self,update,res_path,lattice): diff --git a/python/tests/test_util.py b/python/tests/test_util.py index d5ec36d75..ce5d1b3a9 100644 --- a/python/tests/test_util.py +++ b/python/tests/test_util.py @@ -136,11 +136,13 @@ class TestUtil: ((1,),(7,),False,(7,)), ((1,),(7,),True,(1,7)), ((2,),(2,2),False,(2,2)), - ((1,2),(2,2),False,(2,2)), + ((1,3),(2,3),False,(2,3)), ((1,1,2),(2,2),False,(1,2,2)), ((1,1,2),(2,2),True,(1,1,2,2)), ((1,2,3),(2,3,4),False,(1,2,3,4)), ((1,2,3),(1,2,3),False,(1,2,3)), + ((2,3,1,1),(2,3),False,(2,3,2,3)), + ((2,3,1,1),(2,3),True,(2,3,1,1,2,3)), ]) def test_shapeblender(self,a,b,ones,answer): assert util.shapeblender(a,b,ones) == answer From afe2094c65aa31ffaae83cfe6d975f725a7ce3dc Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Wed, 4 Oct 2023 10:29:55 -0400 Subject: [PATCH 06/13] remove debug print --- python/damask/_rotation.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index 21197f105..7286b3ecf 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -112,7 +112,7 @@ class Rotation: def __getitem__(self, - item: Union[Tuple[int], int, bool, np.bool_, np.ndarray]): + item: Union[Tuple[Union[None, int, slice, ellipsis]], int, bool, np.bool_, np.ndarray]): """ Return self[item]. @@ -445,7 +445,6 @@ class Rotation: obs = util.shapeblender(self.shape,other.shape)[len(self.shape):] for l in [4,2,1]: if obs[-l:] == l*(3,): - print(f'rotate {l}') bs = util.shapeblender(self.shape,other.shape[:-l],False) self_ = self.broadcast_to(bs) if self.shape != bs else self if l==1: From 634f9fd1f54d21ccc9c93fa6528eb9947365fd43 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Wed, 4 Oct 2023 10:30:25 -0400 Subject: [PATCH 07/13] extended slicing types, but have to ignore..! --- python/damask/_orientation.py | 2 +- python/damask/_rotation.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/python/damask/_orientation.py b/python/damask/_orientation.py index ef71366c8..03ef16213 100644 --- a/python/damask/_orientation.py +++ b/python/damask/_orientation.py @@ -529,7 +529,7 @@ class Orientation(Rotation,Crystal): s = self.broadcast_to(s_m).equivalent o = other.broadcast_to(s_o).equivalent - r_ = s[:,np.newaxis,...].misorientation(o[np.newaxis,:,...]) + r_ = s[:,np.newaxis,...].misorientation(o[np.newaxis,:,...]) # type: ignore[index] _r = ~r_ forward = r_.in_disorientation_FZ diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index 7286b3ecf..f4f6a65f3 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -112,7 +112,7 @@ class Rotation: def __getitem__(self, - item: Union[Tuple[Union[None, int, slice, ellipsis]], int, bool, np.bool_, np.ndarray]): + item: Union[Tuple[Union[None, int, slice]], int, bool, np.bool_, np.ndarray]): """ Return self[item]. From b26398319b65b69b590553e53df1316778b4a59d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 5 Oct 2023 16:59:43 +0200 Subject: [PATCH 08/13] testing YAML_types --- src/YAML_parse.f90 | 6 +++--- src/YAML_types.f90 | 7 ++++--- src/crystal.f90 | 2 +- src/test/DAMASK_test.f90 | 5 +++++ src/test/test_YAML_types.f90 | 17 +++++++++++++++++ 5 files changed, 30 insertions(+), 7 deletions(-) create mode 100644 src/test/test_YAML_types.f90 diff --git a/src/YAML_parse.f90 b/src/YAML_parse.f90 index 527a9b3fc..74d9ab083 100644 --- a/src/YAML_parse.f90 +++ b/src/YAML_parse.f90 @@ -48,7 +48,7 @@ subroutine YAML_parse_init() #ifdef FYAML print'(/,1x,a)', 'libfyaml powered' #else - call selfTest() + call YAML_parse_selfTest() #endif end subroutine YAML_parse_init @@ -870,7 +870,7 @@ end function to_flow !-------------------------------------------------------------------------------------------------- !> @brief Check correctness of some YAML functions. !-------------------------------------------------------------------------------------------------- -subroutine selfTest() +subroutine YAML_parse_selfTest() if (indentDepth(' a') /= 1) error stop 'indentDepth' if (indentDepth('a') /= 0) error stop 'indentDepth' @@ -1031,7 +1031,7 @@ subroutine selfTest() end block parse -end subroutine selfTest +end subroutine YAML_parse_selfTest #endif end module YAML_parse diff --git a/src/YAML_types.f90 b/src/YAML_types.f90 index 5d7411dd3..e9904e3c5 100644 --- a/src/YAML_types.f90 +++ b/src/YAML_types.f90 @@ -150,6 +150,7 @@ module YAML_types public :: & YAML_types_init, & + YAML_types_selfTest, & #ifdef __GFORTRAN__ output_as1dStr, & !ToDo: Hack for GNU. Remove later #endif @@ -164,7 +165,7 @@ subroutine YAML_types_init print'(/,1x,a)', '<<<+- YAML_types init -+>>>' - call selfTest() + call YAML_types_selfTest() end subroutine YAML_types_init @@ -172,7 +173,7 @@ end subroutine YAML_types_init !-------------------------------------------------------------------------------------------------- !> @brief Check correctness of some type bound procedures. !-------------------------------------------------------------------------------------------------- -subroutine selfTest() +subroutine YAML_types_selfTest() scalar: block type(tScalar), target :: s @@ -266,7 +267,7 @@ subroutine selfTest() end block dict -end subroutine selfTest +end subroutine YAML_types_selfTest !--------------------------------------------------------------------------------------------------- diff --git a/src/crystal.f90 b/src/crystal.f90 index b405f5afe..32040bce4 100644 --- a/src/crystal.f90 +++ b/src/crystal.f90 @@ -2223,7 +2223,7 @@ end function crystal_isotropic_mu !-------------------------------------------------------------------------------------------------- !> @brief Check correctness of some crystal functions. !-------------------------------------------------------------------------------------------------- -subroutine crystal_selfTest +subroutine crystal_selfTest() real(pREAL), dimension(:,:,:), allocatable :: CoSy real(pREAL), dimension(:,:), allocatable :: system diff --git a/src/test/DAMASK_test.f90 b/src/test/DAMASK_test.f90 index 003519ba3..e2566be85 100644 --- a/src/test/DAMASK_test.f90 +++ b/src/test/DAMASK_test.f90 @@ -12,6 +12,7 @@ program DAMASK_test use test_crystal use test_rotations use test_IO + use test_YAML_types use test_HDF5_utilities external :: quit @@ -57,6 +58,10 @@ program DAMASK_test call test_IO_run() write(IO_STDOUT,fmt='(a)') ok + write(IO_STDOUT,fmt=fmt, advance='no') 'YAML_types','...' + call test_YAML_types_run() + write(IO_STDOUT,fmt='(a)') ok + write(IO_STDOUT,fmt=fmt, advance='no') 'HDF5_utilities','...' call test_HDF5_utilities_run() write(IO_STDOUT,fmt='(a)') ok diff --git a/src/test/test_YAML_types.f90 b/src/test/test_YAML_types.f90 new file mode 100644 index 000000000..4c337afb1 --- /dev/null +++ b/src/test/test_YAML_types.f90 @@ -0,0 +1,17 @@ +module test_YAML_types + use YAML_types + + implicit none(type,external) + + private + public :: test_YAML_types_run + + contains + +subroutine test_YAML_types_run() + + call YAML_types_selfTest() + +end subroutine test_YAML_types_run + +end module test_YAML_types From 23ec58fb4a40408e560fcba88d9bdd5e7f0d9f65 Mon Sep 17 00:00:00 2001 From: Test User Date: Fri, 6 Oct 2023 14:10:11 +0200 Subject: [PATCH 09/13] [skip ci] updated version information after successful test of v3.0.0-alpha7-871-g4f2c726b9 --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index baa2c1a48..de776ef6a 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -3.0.0-alpha7-867-g026dd3c14 +3.0.0-alpha7-871-g4f2c726b9 From 51a35dd17ca2176f671a99b3cef17590b3cfff54 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Fri, 6 Oct 2023 09:41:36 -0400 Subject: [PATCH 10/13] correct in_disorientation_FZ --- python/damask/_orientation.py | 37 ++++++++++++----------------------- 1 file changed, 13 insertions(+), 24 deletions(-) diff --git a/python/damask/_orientation.py b/python/damask/_orientation.py index 03ef16213..7f5ea1cb5 100644 --- a/python/damask/_orientation.py +++ b/python/damask/_orientation.py @@ -434,31 +434,20 @@ class Orientation(Rotation,Crystal): https://doi.org/10.1107/S0108767391006864 """ - with np.errstate(invalid='ignore'): - rho = self.as_Rodrigues_vector(compact=True) - rho_ = rho + np.linalg.norm(rho,axis=-1,keepdims=True)*1e-9 + def larger_or_equal(v,c): + return ((np.isclose(c[0],v[...,0]) | (v[...,0] > c[0])) & + (np.isclose(c[1],v[...,1]) | (v[...,1] > c[1])) & + (np.isclose(c[2],v[...,2]) | (v[...,2] > c[2]))).astype(bool) - if self.family == 'cubic': - return ((rho_[...,0] >= rho[...,1]) & - (rho_[...,1] >= rho[...,2]) & - (rho_[...,2] >= 0)).astype(bool) - if self.family == 'hexagonal': - return ((rho_[...,0] >= rho[...,1]*np.sqrt(3)) & - (rho_[...,1] >= 0) & - (rho_[...,2] >= 0)).astype(bool) - if self.family == 'tetragonal': - return ((rho_[...,0] >= rho[...,1]) & - (rho_[...,1] >= 0) & - (rho_[...,2] >= 0)).astype(bool) - if self.family == 'orthorhombic': - return ((rho_[...,0] >= 0) & - (rho_[...,1] >= 0) & - (rho_[...,2] >= 0)).astype(bool) - if self.family == 'monoclinic': - return ((rho_[...,1] >= 0) & - (rho_[...,2] >= 0)).astype(bool) + rho = self.as_Rodrigues_vector(compact=True) + return larger_or_equal(rho, + [rho[...,1], rho[...,2],0] if self.family == 'cubic' + else [rho[...,1]*np.sqrt(3),0, 0] if self.family == 'hexagonal' + else [rho[...,1], 0, 0] if self.family == 'tetragonal' + else [0, 0, 0] if self.family == 'orthorhombic' + else [-np.inf, 0, 0] if self.family == 'monoclinic' + else [-np.inf, -np.inf, -np.inf]) & self.in_FZ - return np.ones_like(rho[...,0],dtype=bool) def disorientation(self, other: 'Orientation', @@ -534,7 +523,7 @@ class Orientation(Rotation,Crystal): forward = r_.in_disorientation_FZ reverse = _r.in_disorientation_FZ - ok = (forward | reverse) & r_.in_FZ + ok = forward | reverse ok &= (np.cumsum(ok.reshape((-1,)+blend),axis=0) == 1).reshape(ok.shape) r = np.where(np.any((ok&forward)[...,np.newaxis],axis=(0,1),keepdims=True), r_.quaternion, From 0c7384b95dcb8ea9092ecee8897c0fb798fed3f5 Mon Sep 17 00:00:00 2001 From: Test User Date: Sat, 7 Oct 2023 02:22:05 +0200 Subject: [PATCH 11/13] [skip ci] updated version information after successful test of v3.0.0-alpha7-877-g5fc3b23cc --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index de776ef6a..fdf6ecbb0 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -3.0.0-alpha7-871-g4f2c726b9 +3.0.0-alpha7-877-g5fc3b23cc From c48cee1af2f9fa8ef623826d107f7c4a160c974c Mon Sep 17 00:00:00 2001 From: Test User Date: Mon, 9 Oct 2023 17:09:09 +0200 Subject: [PATCH 12/13] [skip ci] updated version information after successful test of v3.0.0-alpha7-880-g3644fe586 --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index fdf6ecbb0..2dfbe91c2 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -3.0.0-alpha7-877-g5fc3b23cc +3.0.0-alpha7-880-g3644fe586 From 57e02bb991033281e2280ae9b1a0f5115d7a99ab Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 8 Oct 2023 09:48:53 +0200 Subject: [PATCH 13/13] PETSc 3.20 --- CMakeLists.txt | 2 +- src/CLI.f90 | 2 +- src/mesh/mesh_mech_FEM.f90 | 3 +++ 3 files changed, 5 insertions(+), 2 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index de04a95e3..60a05cfc7 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -11,7 +11,7 @@ endif() project(Prerequisites LANGUAGES) set(ENV{PKG_CONFIG_PATH} "$ENV{PETSC_DIR}/$ENV{PETSC_ARCH}/lib/pkgconfig:$ENV{PKG_CONFIG_PATH}") pkg_check_modules(PETSC_MIN REQUIRED PETSc>=3.12.0 QUIET) #CMake does not support version range -pkg_check_modules(PETSC REQUIRED PETSc<3.20.0) +pkg_check_modules(PETSC REQUIRED PETSc<3.21.0) pkg_get_variable(CMAKE_Fortran_COMPILER PETSc fcompiler) pkg_get_variable(CMAKE_C_COMPILER PETSc ccompiler) diff --git a/src/CLI.f90 b/src/CLI.f90 index 675f3c3bc..968186f3d 100644 --- a/src/CLI.f90 +++ b/src/CLI.f90 @@ -6,7 +6,7 @@ !> @brief Parse command line interface for PETSc-based solvers !-------------------------------------------------------------------------------------------------- #define PETSC_MINOR_MIN 12 -#define PETSC_MINOR_MAX 19 +#define PETSC_MINOR_MAX 20 module CLI use, intrinsic :: ISO_fortran_env diff --git a/src/mesh/mesh_mech_FEM.f90 b/src/mesh/mesh_mech_FEM.f90 index c566990e2..67b7859ba 100644 --- a/src/mesh/mesh_mech_FEM.f90 +++ b/src/mesh/mesh_mech_FEM.f90 @@ -74,6 +74,9 @@ module mesh_mechanical_FEM external :: & ! ToDo: write interfaces #if defined(PETSC_USE_64BIT_INDICES) || PETSC_VERSION_MINOR < 17 ISDestroy, & +#endif +#if PETSC_VERSION_MINOR > 18 + DMAddField, & #endif PetscSectionGetNumFields, & PetscFESetQuadrature, &