From 0c41e334341da15de6033709f05657474b7bf576 Mon Sep 17 00:00:00 2001 From: "f.basile" Date: Thu, 4 Jun 2020 16:39:24 +0200 Subject: [PATCH 01/40] orientation.equivalent takes several rotations at the same time + small test" " "" --- python/damask/_orientation.py | 18 ++++++++++++++++-- python/damask/oritest.py | 33 +++++++++++++++++++++++++++++++++ 2 files changed, 49 insertions(+), 2 deletions(-) create mode 100644 python/damask/oritest.py diff --git a/python/damask/_orientation.py b/python/damask/_orientation.py index e4dbc7d7c..6fa40c261 100644 --- a/python/damask/_orientation.py +++ b/python/damask/_orientation.py @@ -39,8 +39,8 @@ class Orientation: else: self.rotation = Rotation.from_quaternion(rotation) # assume quaternion - if self.rotation.quaternion.shape != (4,): - raise NotImplementedError('Support for multiple rotations missing') +# if self.rotation.quaternion.shape != (4,): +# raise NotImplementedError('Support for multiple rotations missing') def disorientation(self, other, @@ -80,6 +80,20 @@ class Orientation: def inFZ(self): return self.lattice.symmetry.inFZ(self.rotation.as_Rodrigues(vector=True)) + def equivalent(self): + """ + List of orientations which are symmetrically equivalent. + Supported for multiple rotation with same lattice + Returns list [i] being i=range(24) + Returns list [i, num_rot] for multiple rotations + """ + if not self.rotation.shape: + return [self.__class__(q*self.rotation,self.lattice) \ + for q in self.lattice.symmetry.symmetryOperations()] + else: + return np.reshape([self.__class__(q*Rotation.from_quaternion(self.rotation.as_quaternion()[l]),self.lattice) \ + for q in self.lattice.symmetry.symmetryOperations() for l in range(self.rotation.shape[0])], (24,self.rotation.shape[0])) + def equivalentOrientations(self,members=[]): """List of orientations which are symmetrically equivalent.""" diff --git a/python/damask/oritest.py b/python/damask/oritest.py new file mode 100644 index 000000000..43104a938 --- /dev/null +++ b/python/damask/oritest.py @@ -0,0 +1,33 @@ +import damask +import numpy as np + + +rot0= damask.Rotation.from_random() +rot1= damask.Rotation.from_random() +rot2= damask.Rotation.from_random() + +ori0=damask.Orientation(rot0,'fcc') +ori1=damask.Orientation(rot1,'fcc') +ori2=damask.Orientation(rot2,'fcc') + + + + +quat=np.array([rot0.as_quaternion(),rot1.as_quaternion(),rot2.as_quaternion()]) +rot=damask.Rotation.from_quaternion(quat) + +ori=damask.Orientation(rot,'fcc') + +ori.equivalent() + + + +# doesn't work this way, don't know why +#ori.equivalent()[:,0][0] == ori0.equivalentOrientations()[0] + +for s in range(24): + print(ori.equivalent()[s,0].rotation.as_Eulers() == ori0.equivalentOrientations()[s].rotation.as_Eulers()) + print(ori.equivalent()[s,1].rotation.as_Eulers() == ori1.equivalentOrientations()[s].rotation.as_Eulers()) + print(ori.equivalent()[s,2].rotation.as_Eulers() == ori2.equivalentOrientations()[s].rotation.as_Eulers()) + + From 3897136f85884c745a3e645a93a9fae136db317d Mon Sep 17 00:00:00 2001 From: "f.basile" Date: Thu, 4 Jun 2020 16:43:28 +0200 Subject: [PATCH 02/40] avoid python/damask/_orientation.py exceeds line length limit (maximum line length 141 > 132) --- python/damask/_orientation.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/python/damask/_orientation.py b/python/damask/_orientation.py index 6fa40c261..909365d41 100644 --- a/python/damask/_orientation.py +++ b/python/damask/_orientation.py @@ -92,7 +92,8 @@ class Orientation: for q in self.lattice.symmetry.symmetryOperations()] else: return np.reshape([self.__class__(q*Rotation.from_quaternion(self.rotation.as_quaternion()[l]),self.lattice) \ - for q in self.lattice.symmetry.symmetryOperations() for l in range(self.rotation.shape[0])], (24,self.rotation.shape[0])) + for q in self.lattice.symmetry.symmetryOperations() \ + for l in range(self.rotation.shape[0])], (24,self.rotation.shape[0])) def equivalentOrientations(self,members=[]): From 6a24aee171122600a26a8d65c0c1acbd06741ee6 Mon Sep 17 00:00:00 2001 From: "f.basile" Date: Thu, 4 Jun 2020 16:44:57 +0200 Subject: [PATCH 03/40] fix python/damask/_orientation.py contains invalid python3 code --- python/damask/_orientation.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/python/damask/_orientation.py b/python/damask/_orientation.py index 909365d41..f338be35d 100644 --- a/python/damask/_orientation.py +++ b/python/damask/_orientation.py @@ -82,10 +82,12 @@ class Orientation: def equivalent(self): """ + List of orientations which are symmetrically equivalent. Supported for multiple rotation with same lattice Returns list [i] being i=range(24) Returns list [i, num_rot] for multiple rotations + """ if not self.rotation.shape: return [self.__class__(q*self.rotation,self.lattice) \ From ac09a2912a434bc6f2289eed849b7bf1eaa43429 Mon Sep 17 00:00:00 2001 From: "f.basile" Date: Thu, 4 Jun 2020 16:45:44 +0200 Subject: [PATCH 04/40] fix python/damask/_orientation.py contains invalid python3 code 2 --- python/damask/_orientation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/_orientation.py b/python/damask/_orientation.py index f338be35d..ec83c3857 100644 --- a/python/damask/_orientation.py +++ b/python/damask/_orientation.py @@ -82,8 +82,8 @@ class Orientation: def equivalent(self): """ - List of orientations which are symmetrically equivalent. + Supported for multiple rotation with same lattice Returns list [i] being i=range(24) Returns list [i, num_rot] for multiple rotations From eae9698d22c527502d5a8cefd05d50a8520479bb Mon Sep 17 00:00:00 2001 From: "f.basile" Date: Fri, 5 Jun 2020 13:48:12 +0200 Subject: [PATCH 05/40] equivalent,related and inFZ vectorized + pytest --- python/damask/_orientation.py | 44 +++++++++++++---- python/damask/oritest.py | 33 ------------- python/tests/test_ori_vec.py | 93 +++++++++++++++++++++++++++++++++++ 3 files changed, 127 insertions(+), 43 deletions(-) delete mode 100644 python/damask/oritest.py create mode 100644 python/tests/test_ori_vec.py diff --git a/python/damask/_orientation.py b/python/damask/_orientation.py index ec83c3857..245659f44 100644 --- a/python/damask/_orientation.py +++ b/python/damask/_orientation.py @@ -77,25 +77,37 @@ class Orientation: return (Orientation(r,self.lattice), i,j, k == 1) if symmetries else r # disorientation ... # ... own sym, other sym, # self-->other: True, self<--other: False + + def inFZ_vec(self): + """ + Check if orientations falls into Fundamental Zone + + self.rotation.as_Rodrigues() working fine + self.rotation.as_Rodrigues(vector=True) doesn't work for several rotations + i apply dirty fix + + """ + if not self.rotation.shape: + return self.lattice.symmetry.inFZ(self.rotation.as_Rodrigues(vector=True)) + else: + return [self.lattice.symmetry.inFZ(\ + Rotation._qu2ro(self.rotation.as_quaternion())[l][...,:3]\ + *Rotation._qu2ro(self.rotation.as_quaternion())[l][...,3])\ + for l in range(self.rotation.shape[0])] + def inFZ(self): return self.lattice.symmetry.inFZ(self.rotation.as_Rodrigues(vector=True)) - def equivalent(self): - """ - List of orientations which are symmetrically equivalent. - - Supported for multiple rotation with same lattice - Returns list [i] being i=range(24) - Returns list [i, num_rot] for multiple rotations - - """ + def equivalent_vec(self): + """List of orientations which are symmetrically equivalent.""" if not self.rotation.shape: return [self.__class__(q*self.rotation,self.lattice) \ for q in self.lattice.symmetry.symmetryOperations()] else: return np.reshape([self.__class__(q*Rotation.from_quaternion(self.rotation.as_quaternion()[l]),self.lattice) \ for q in self.lattice.symmetry.symmetryOperations() \ - for l in range(self.rotation.shape[0])], (24,self.rotation.shape[0])) + for l in range(self.rotation.shape[0])], \ + (len(self.lattice.symmetry.symmetryOperations()),self.rotation.shape[0])) def equivalentOrientations(self,members=[]): @@ -108,6 +120,18 @@ class Orientation: return [self.__class__(q*self.rotation,self.lattice) \ for q in self.lattice.symmetry.symmetryOperations(members)] # yes, return list of rotations + def relatedOrientations_vec(self,model): + """List of orientations related by the given orientation relationship.""" + r = self.lattice.relationOperations(model) + if not self.rotation.shape: + return [self.__class__(o*self.rotation,r['lattice']) for o in r['rotations']] + else: + return np.reshape(\ + [self.__class__(o*Rotation.from_quaternion(self.rotation.as_quaternion()[l])\ + ,r['lattice']) for o in r['rotations'] for l in range(self.rotation.shape[0])] + ,(len(r['rotations']),self.rotation.shape[0])) + + def relatedOrientations(self,model): """List of orientations related by the given orientation relationship.""" r = self.lattice.relationOperations(model) diff --git a/python/damask/oritest.py b/python/damask/oritest.py deleted file mode 100644 index 43104a938..000000000 --- a/python/damask/oritest.py +++ /dev/null @@ -1,33 +0,0 @@ -import damask -import numpy as np - - -rot0= damask.Rotation.from_random() -rot1= damask.Rotation.from_random() -rot2= damask.Rotation.from_random() - -ori0=damask.Orientation(rot0,'fcc') -ori1=damask.Orientation(rot1,'fcc') -ori2=damask.Orientation(rot2,'fcc') - - - - -quat=np.array([rot0.as_quaternion(),rot1.as_quaternion(),rot2.as_quaternion()]) -rot=damask.Rotation.from_quaternion(quat) - -ori=damask.Orientation(rot,'fcc') - -ori.equivalent() - - - -# doesn't work this way, don't know why -#ori.equivalent()[:,0][0] == ori0.equivalentOrientations()[0] - -for s in range(24): - print(ori.equivalent()[s,0].rotation.as_Eulers() == ori0.equivalentOrientations()[s].rotation.as_Eulers()) - print(ori.equivalent()[s,1].rotation.as_Eulers() == ori1.equivalentOrientations()[s].rotation.as_Eulers()) - print(ori.equivalent()[s,2].rotation.as_Eulers() == ori2.equivalentOrientations()[s].rotation.as_Eulers()) - - diff --git a/python/tests/test_ori_vec.py b/python/tests/test_ori_vec.py new file mode 100644 index 000000000..ffc2c88b1 --- /dev/null +++ b/python/tests/test_ori_vec.py @@ -0,0 +1,93 @@ +import os +from itertools import permutations + +import pytest +import numpy as np + +import damask +from damask import Rotation +from damask import Orientation +from damask import Lattice + +rot0= damask.Rotation.from_random() +rot1= damask.Rotation.from_random() +rot2= damask.Rotation.from_random() +rot3= damask.Rotation.from_random() + +class TestOrientation_vec: + @pytest.mark.parametrize('lattice',Lattice.lattices) + def test_equivalentOrientations_vec(self,lattice): + ori0=damask.Orientation(rot0,lattice) + ori1=damask.Orientation(rot1,lattice) + ori2=damask.Orientation(rot2,lattice) + ori3=damask.Orientation(rot3,lattice) + + quat=np.array([rot0.as_quaternion(),rot1.as_quaternion(),rot2.as_quaternion(),rot3.as_quaternion()]) + rot_vec=damask.Rotation.from_quaternion(quat) + ori_vec=damask.Orientation(rot_vec,lattice) + + for s in range(len(ori_vec.lattice.symmetry.symmetryOperations())): + assert all(ori_vec.equivalent_vec()[s,0].rotation.as_Eulers() == \ + ori0.equivalentOrientations()[s].rotation.as_Eulers()) + assert all(ori_vec.equivalent_vec()[s,1].rotation.as_quaternion() == \ + ori1.equivalentOrientations()[s].rotation.as_quaternion()) + assert all(ori_vec.equivalent_vec()[s,2].rotation.as_Rodrigues() == \ + ori2.equivalentOrientations()[s].rotation.as_Rodrigues()) + assert all(ori_vec.equivalent_vec()[s,3].rotation.as_cubochoric() == \ + ori3.equivalentOrientations()[s].rotation.as_cubochoric()) + + @pytest.mark.parametrize('lattice',Lattice.lattices) + def test_inFZ_vec(self,lattice): + ori0=damask.Orientation(rot0,lattice) + ori1=damask.Orientation(rot1,lattice) + ori2=damask.Orientation(rot2,lattice) + ori3=damask.Orientation(rot3,lattice) + #ensure 1 of them is in FZ + ori4=ori0.reduced() + rot4=ori4.rotation + + quat=np.array([rot0.as_quaternion(),rot1.as_quaternion(),\ + rot2.as_quaternion(),rot3.as_quaternion(), rot4.as_quaternion()]) + rot_vec=damask.Rotation.from_quaternion(quat) + ori_vec=damask.Orientation(rot_vec,lattice) + + assert ori_vec.inFZ_vec()[0] == ori0.inFZ() + assert ori_vec.inFZ_vec()[1] == ori1.inFZ() + assert ori_vec.inFZ_vec()[2] == ori2.inFZ() + assert ori_vec.inFZ_vec()[3] == ori3.inFZ() + assert ori_vec.inFZ_vec()[4] == ori4.inFZ() + + + @pytest.mark.parametrize('model',['Bain','KS','GT','GT_prime','NW','Pitsch']) + @pytest.mark.parametrize('lattice',['fcc','bcc']) + def test_relatedOrientations_vec(self,model,lattice): + ori0=damask.Orientation(rot0,lattice) + ori1=damask.Orientation(rot1,lattice) + ori2=damask.Orientation(rot2,lattice) + ori3=damask.Orientation(rot3,lattice) + + quat=np.array([rot0.as_quaternion(),rot1.as_quaternion(),rot2.as_quaternion(),rot3.as_quaternion()]) + rot_vec=damask.Rotation.from_quaternion(quat) + ori_vec=damask.Orientation(rot_vec,lattice) + + for s in range(len(ori1.lattice.relationOperations(model)['rotations'])): + assert all(ori_vec.relatedOrientations_vec(model)[s,0].rotation.as_Eulers() == \ + ori0.relatedOrientations(model)[s].rotation.as_Eulers()) + assert all(ori_vec.relatedOrientations_vec(model)[s,1].rotation.as_quaternion() == \ + ori1.relatedOrientations(model)[s].rotation.as_quaternion()) + assert all(ori_vec.relatedOrientations_vec(model)[s,2].rotation.as_Rodrigues() == \ + ori2.relatedOrientations(model)[s].rotation.as_Rodrigues()) + assert all(ori_vec.relatedOrientations_vec(model)[s,3].rotation.as_cubochoric() == \ + ori3.relatedOrientations(model)[s].rotation.as_cubochoric()) + + + + + + + + + + + + \ No newline at end of file From cb1779ef9adf66b9c693473e79dbeea0ad372e53 Mon Sep 17 00:00:00 2001 From: "f.basile" Date: Fri, 5 Jun 2020 13:49:30 +0200 Subject: [PATCH 06/40] fix pep257: D415 / First line should end with a period, question mark, or exclamation point (not e) --- python/damask/_orientation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/_orientation.py b/python/damask/_orientation.py index 245659f44..10d33a76c 100644 --- a/python/damask/_orientation.py +++ b/python/damask/_orientation.py @@ -80,7 +80,7 @@ class Orientation: def inFZ_vec(self): """ - Check if orientations falls into Fundamental Zone + Check if orientations falls into Fundamental Zone. self.rotation.as_Rodrigues() working fine self.rotation.as_Rodrigues(vector=True) doesn't work for several rotations From 54c20cdd3c8cc7a3b3ff3a9a84a507619d7231ec Mon Sep 17 00:00:00 2001 From: "f.basile" Date: Fri, 5 Jun 2020 13:53:47 +0200 Subject: [PATCH 07/40] fix pyflakes imported but unused in pytest --- python/tests/test_ori_vec.py | 47 +++++++++++++++++------------------- 1 file changed, 22 insertions(+), 25 deletions(-) diff --git a/python/tests/test_ori_vec.py b/python/tests/test_ori_vec.py index ffc2c88b1..4172e5edc 100644 --- a/python/tests/test_ori_vec.py +++ b/python/tests/test_ori_vec.py @@ -1,6 +1,3 @@ -import os -from itertools import permutations - import pytest import numpy as np @@ -9,22 +6,22 @@ from damask import Rotation from damask import Orientation from damask import Lattice -rot0= damask.Rotation.from_random() -rot1= damask.Rotation.from_random() -rot2= damask.Rotation.from_random() -rot3= damask.Rotation.from_random() +rot0= Rotation.from_random() +rot1= Rotation.from_random() +rot2= Rotation.from_random() +rot3= Rotation.from_random() class TestOrientation_vec: @pytest.mark.parametrize('lattice',Lattice.lattices) def test_equivalentOrientations_vec(self,lattice): - ori0=damask.Orientation(rot0,lattice) - ori1=damask.Orientation(rot1,lattice) - ori2=damask.Orientation(rot2,lattice) - ori3=damask.Orientation(rot3,lattice) + ori0=Orientation(rot0,lattice) + ori1=Orientation(rot1,lattice) + ori2=Orientation(rot2,lattice) + ori3=Orientation(rot3,lattice) quat=np.array([rot0.as_quaternion(),rot1.as_quaternion(),rot2.as_quaternion(),rot3.as_quaternion()]) - rot_vec=damask.Rotation.from_quaternion(quat) - ori_vec=damask.Orientation(rot_vec,lattice) + rot_vec=Rotation.from_quaternion(quat) + ori_vec=Orientation(rot_vec,lattice) for s in range(len(ori_vec.lattice.symmetry.symmetryOperations())): assert all(ori_vec.equivalent_vec()[s,0].rotation.as_Eulers() == \ @@ -38,18 +35,18 @@ class TestOrientation_vec: @pytest.mark.parametrize('lattice',Lattice.lattices) def test_inFZ_vec(self,lattice): - ori0=damask.Orientation(rot0,lattice) - ori1=damask.Orientation(rot1,lattice) - ori2=damask.Orientation(rot2,lattice) - ori3=damask.Orientation(rot3,lattice) + ori0=Orientation(rot0,lattice) + ori1=Orientation(rot1,lattice) + ori2=Orientation(rot2,lattice) + ori3=Orientation(rot3,lattice) #ensure 1 of them is in FZ ori4=ori0.reduced() rot4=ori4.rotation quat=np.array([rot0.as_quaternion(),rot1.as_quaternion(),\ rot2.as_quaternion(),rot3.as_quaternion(), rot4.as_quaternion()]) - rot_vec=damask.Rotation.from_quaternion(quat) - ori_vec=damask.Orientation(rot_vec,lattice) + rot_vec=Rotation.from_quaternion(quat) + ori_vec=Orientation(rot_vec,lattice) assert ori_vec.inFZ_vec()[0] == ori0.inFZ() assert ori_vec.inFZ_vec()[1] == ori1.inFZ() @@ -61,14 +58,14 @@ class TestOrientation_vec: @pytest.mark.parametrize('model',['Bain','KS','GT','GT_prime','NW','Pitsch']) @pytest.mark.parametrize('lattice',['fcc','bcc']) def test_relatedOrientations_vec(self,model,lattice): - ori0=damask.Orientation(rot0,lattice) - ori1=damask.Orientation(rot1,lattice) - ori2=damask.Orientation(rot2,lattice) - ori3=damask.Orientation(rot3,lattice) + ori0=Orientation(rot0,lattice) + ori1=Orientation(rot1,lattice) + ori2=Orientation(rot2,lattice) + ori3=Orientation(rot3,lattice) quat=np.array([rot0.as_quaternion(),rot1.as_quaternion(),rot2.as_quaternion(),rot3.as_quaternion()]) - rot_vec=damask.Rotation.from_quaternion(quat) - ori_vec=damask.Orientation(rot_vec,lattice) + rot_vec=Rotation.from_quaternion(quat) + ori_vec=Orientation(rot_vec,lattice) for s in range(len(ori1.lattice.relationOperations(model)['rotations'])): assert all(ori_vec.relatedOrientations_vec(model)[s,0].rotation.as_Eulers() == \ From a76b5233be7639b6ed3d82d34e48d88214ecd835 Mon Sep 17 00:00:00 2001 From: "f.basile" Date: Fri, 5 Jun 2020 13:54:38 +0200 Subject: [PATCH 08/40] fix pyflakes imported but unused in pytest 2 --- python/tests/test_ori_vec.py | 1 - 1 file changed, 1 deletion(-) diff --git a/python/tests/test_ori_vec.py b/python/tests/test_ori_vec.py index 4172e5edc..f55dc6d03 100644 --- a/python/tests/test_ori_vec.py +++ b/python/tests/test_ori_vec.py @@ -1,7 +1,6 @@ import pytest import numpy as np -import damask from damask import Rotation from damask import Orientation from damask import Lattice From fe5e5babfea196b4d95e7ff491b3aec18cae872b Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 18 Jun 2020 22:29:28 +0200 Subject: [PATCH 09/40] more useful for vectorized calculations --- python/damask/_lattice.py | 71 ++++++++++++++++++++++++++++++++++++++- 1 file changed, 70 insertions(+), 1 deletion(-) diff --git a/python/damask/_lattice.py b/python/damask/_lattice.py index 0f6cffd04..09eadee5a 100644 --- a/python/damask/_lattice.py +++ b/python/damask/_lattice.py @@ -13,7 +13,7 @@ class Symmetry: """ - lattices = [None,'orthorhombic','tetragonal','hexagonal','cubic',] + lattices = [None,'orthorhombic','tetragonal','hexagonal','cubic'] def __init__(self, symmetry = None): """ @@ -157,6 +157,75 @@ class Symmetry: else: return symOps # yes, return list of rotations + @property + def symmetry_operations(self): + """List (or single element) of symmetry operations as rotations.""" + if self.lattice == 'cubic': + symQuats = [ + [ 1.0, 0.0, 0.0, 0.0 ], + [ 0.0, 1.0, 0.0, 0.0 ], + [ 0.0, 0.0, 1.0, 0.0 ], + [ 0.0, 0.0, 0.0, 1.0 ], + [ 0.0, 0.0, 0.5*np.sqrt(2), 0.5*np.sqrt(2) ], + [ 0.0, 0.0, 0.5*np.sqrt(2),-0.5*np.sqrt(2) ], + [ 0.0, 0.5*np.sqrt(2), 0.0, 0.5*np.sqrt(2) ], + [ 0.0, 0.5*np.sqrt(2), 0.0, -0.5*np.sqrt(2) ], + [ 0.0, 0.5*np.sqrt(2),-0.5*np.sqrt(2), 0.0 ], + [ 0.0, -0.5*np.sqrt(2),-0.5*np.sqrt(2), 0.0 ], + [ 0.5, 0.5, 0.5, 0.5 ], + [-0.5, 0.5, 0.5, 0.5 ], + [-0.5, 0.5, 0.5, -0.5 ], + [-0.5, 0.5, -0.5, 0.5 ], + [-0.5, -0.5, 0.5, 0.5 ], + [-0.5, -0.5, 0.5, -0.5 ], + [-0.5, -0.5, -0.5, 0.5 ], + [-0.5, 0.5, -0.5, -0.5 ], + [-0.5*np.sqrt(2), 0.0, 0.0, 0.5*np.sqrt(2) ], + [ 0.5*np.sqrt(2), 0.0, 0.0, 0.5*np.sqrt(2) ], + [-0.5*np.sqrt(2), 0.0, 0.5*np.sqrt(2), 0.0 ], + [-0.5*np.sqrt(2), 0.0, -0.5*np.sqrt(2), 0.0 ], + [-0.5*np.sqrt(2), 0.5*np.sqrt(2), 0.0, 0.0 ], + [-0.5*np.sqrt(2),-0.5*np.sqrt(2), 0.0, 0.0 ], + ] + elif self.lattice == 'hexagonal': + symQuats = [ + [ 1.0, 0.0, 0.0, 0.0 ], + [-0.5*np.sqrt(3), 0.0, 0.0, -0.5 ], + [ 0.5, 0.0, 0.0, 0.5*np.sqrt(3) ], + [ 0.0, 0.0, 0.0, 1.0 ], + [-0.5, 0.0, 0.0, 0.5*np.sqrt(3) ], + [-0.5*np.sqrt(3), 0.0, 0.0, 0.5 ], + [ 0.0, 1.0, 0.0, 0.0 ], + [ 0.0, -0.5*np.sqrt(3), 0.5, 0.0 ], + [ 0.0, 0.5, -0.5*np.sqrt(3), 0.0 ], + [ 0.0, 0.0, 1.0, 0.0 ], + [ 0.0, -0.5, -0.5*np.sqrt(3), 0.0 ], + [ 0.0, 0.5*np.sqrt(3), 0.5, 0.0 ], + ] + elif self.lattice == 'tetragonal': + symQuats = [ + [ 1.0, 0.0, 0.0, 0.0 ], + [ 0.0, 1.0, 0.0, 0.0 ], + [ 0.0, 0.0, 1.0, 0.0 ], + [ 0.0, 0.0, 0.0, 1.0 ], + [ 0.0, 0.5*np.sqrt(2), 0.5*np.sqrt(2), 0.0 ], + [ 0.0, -0.5*np.sqrt(2), 0.5*np.sqrt(2), 0.0 ], + [ 0.5*np.sqrt(2), 0.0, 0.0, 0.5*np.sqrt(2) ], + [-0.5*np.sqrt(2), 0.0, 0.0, 0.5*np.sqrt(2) ], + ] + elif self.lattice == 'orthorhombic': + symQuats = [ + [ 1.0,0.0,0.0,0.0 ], + [ 0.0,1.0,0.0,0.0 ], + [ 0.0,0.0,1.0,0.0 ], + [ 0.0,0.0,0.0,1.0 ], + ] + else: + symQuats = [ + [ 1.0,0.0,0.0,0.0 ], + ] + return Rotation(np.array(symQuats)) + def inFZ(self,rodrigues): """ From cdda556e18b7b60121a6ce3117fd30c20d6eeaf0 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 18 Jun 2020 22:30:22 +0200 Subject: [PATCH 10/40] small improvements - hack for reporting multiple rotation - bugfix for vectorized Rodrigues vector - more general broadcasting (even more powerfull then np.broadcast_to) --- python/damask/_rotation.py | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index b1b1cd5ad..0d92e19e9 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -70,7 +70,7 @@ class Rotation: def __repr__(self): """Orientation displayed as unit quaternion, rotation matrix, and Bunge-Euler angles.""" if self.quaternion.shape != (4,): - raise NotImplementedError('Support for multiple rotations missing') + return str(self.quaternion) # ToDo: could be nicer ... return '\n'.join([ 'Quaternion: (real={:.3f}, imag=<{:+.3f}, {:+.3f}, {:+.3f}>)'.format(*(self.quaternion)), 'Matrix:\n{}'.format(self.as_matrix()), @@ -159,13 +159,12 @@ class Rotation: def broadcast_to(self,shape): if isinstance(shape,int): shape = (shape,) - if self.shape == (): - q = np.broadcast_to(self.quaternion,shape+(4,)) - else: - q = np.block([np.broadcast_to(self.quaternion[...,0:1],shape).reshape(shape+(1,)), - np.broadcast_to(self.quaternion[...,1:2],shape).reshape(shape+(1,)), - np.broadcast_to(self.quaternion[...,2:3],shape).reshape(shape+(1,)), - np.broadcast_to(self.quaternion[...,3:4],shape).reshape(shape+(1,))]) + N = np.prod(shape)//np.prod(self.shape,dtype=int) + + q = np.block([np.repeat(self.quaternion[...,0:1],N).reshape(shape+(1,)), + np.repeat(self.quaternion[...,1:2],N).reshape(shape+(1,)), + np.repeat(self.quaternion[...,2:3],N).reshape(shape+(1,)), + np.repeat(self.quaternion[...,3:4],N).reshape(shape+(1,))]) return self.__class__(q) @@ -248,7 +247,7 @@ class Rotation: """ ro = Rotation._qu2ro(self.quaternion) - return ro[...,:3]*ro[...,3] if vector else ro + return ro[...,:3]*ro[...,3:] if vector else ro def as_homochoric(self): """Homochoric vector: (h_1, h_2, h_3).""" From 1648963b5782d0de380cdfca24581dee9fedae4b Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 18 Jun 2020 22:53:04 +0200 Subject: [PATCH 11/40] vectorized equivalent orientation calculation --- python/damask/_lattice.py | 4 ++-- python/damask/_orientation.py | 33 ++++++++++++++++++--------------- python/damask/_rotation.py | 15 ++++++++++----- python/tests/test_ori_vec.py | 34 ++++++++++++---------------------- 4 files changed, 42 insertions(+), 44 deletions(-) diff --git a/python/damask/_lattice.py b/python/damask/_lattice.py index 09eadee5a..9769933ef 100644 --- a/python/damask/_lattice.py +++ b/python/damask/_lattice.py @@ -159,7 +159,7 @@ class Symmetry: @property def symmetry_operations(self): - """List (or single element) of symmetry operations as rotations.""" + """Symmetry operations as Rotations.""" if self.lattice == 'cubic': symQuats = [ [ 1.0, 0.0, 0.0, 0.0 ], @@ -236,7 +236,7 @@ class Symmetry: if (len(rodrigues) != 3): raise ValueError('Input is not a Rodrigues-Frank vector.\n') - if np.any(rodrigues == np.inf): return False + if np.any(rodrigues == np.inf): return False # ToDo: MD: not sure if needed Rabs = abs(rodrigues) diff --git a/python/damask/_orientation.py b/python/damask/_orientation.py index 10d33a76c..287e84ff1 100644 --- a/python/damask/_orientation.py +++ b/python/damask/_orientation.py @@ -3,7 +3,7 @@ import numpy as np from . import Lattice from . import Rotation -class Orientation: +class Orientation: # make subclass or Rotation? """ Crystallographic orientation. @@ -39,8 +39,6 @@ class Orientation: else: self.rotation = Rotation.from_quaternion(rotation) # assume quaternion -# if self.rotation.quaternion.shape != (4,): -# raise NotImplementedError('Support for multiple rotations missing') def disorientation(self, other, @@ -94,20 +92,25 @@ class Orientation: Rotation._qu2ro(self.rotation.as_quaternion())[l][...,:3]\ *Rotation._qu2ro(self.rotation.as_quaternion())[l][...,3])\ for l in range(self.rotation.shape[0])] - + def inFZ(self): return self.lattice.symmetry.inFZ(self.rotation.as_Rodrigues(vector=True)) - def equivalent_vec(self): - """List of orientations which are symmetrically equivalent.""" - if not self.rotation.shape: - return [self.__class__(q*self.rotation,self.lattice) \ - for q in self.lattice.symmetry.symmetryOperations()] - else: - return np.reshape([self.__class__(q*Rotation.from_quaternion(self.rotation.as_quaternion()[l]),self.lattice) \ - for q in self.lattice.symmetry.symmetryOperations() \ - for l in range(self.rotation.shape[0])], \ - (len(self.lattice.symmetry.symmetryOperations()),self.rotation.shape[0])) + @property + def equivalent(self): + """ + Return orientations which are symmetrically equivalent. + + One dimension (length according to symmetrically equivalent orientations) + is added to the left of the rotation array. + + """ + symmetry_operations = self.lattice.symmetry.symmetry_operations + + q = np.block([self.rotation.quaternion]*symmetry_operations.shape[0]) + r = Rotation(q.reshape(symmetry_operations.shape+self.rotation.quaternion.shape)) + + return self.__class__(symmetry_operations.broadcast_to(r.shape)@r,self.lattice) def equivalentOrientations(self,members=[]): @@ -130,7 +133,7 @@ class Orientation: [self.__class__(o*Rotation.from_quaternion(self.rotation.as_quaternion()[l])\ ,r['lattice']) for o in r['rotations'] for l in range(self.rotation.shape[0])] ,(len(r['rotations']),self.rotation.shape[0])) - + def relatedOrientations(self,model): """List of orientations related by the given orientation relationship.""" diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index 0d92e19e9..9a0266871 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -129,6 +129,7 @@ class Rotation: self.quaternion[...,1:] *= -1 return self + #@property def inversed(self): """Inverse rotation/backward rotation.""" return self.copy().inverse() @@ -139,6 +140,7 @@ class Rotation: self.quaternion[self.quaternion[...,0] < 0.0] *= -1 return self + #@property def standardized(self): """Quaternion representation with positive real part.""" return self.copy().standardize() @@ -154,11 +156,12 @@ class Rotation: Rotation to which the misorientation is computed. """ - return other*self.inversed() + return other@self.inversed() def broadcast_to(self,shape): - if isinstance(shape,int): shape = (shape,) + if isinstance(shape,int): + shape = (shape,) N = np.prod(shape)//np.prod(self.shape,dtype=int) q = np.block([np.repeat(self.quaternion[...,0:1],N).reshape(shape+(1,)), @@ -257,6 +260,7 @@ class Rotation: """Cubochoric vector: (c_1, c_2, c_3).""" return Rotation._qu2cu(self.quaternion) + @property def M(self): # ToDo not sure about the name: as_M or M? we do not have a from_M """ Intermediate representation supporting quaternion averaging. @@ -435,8 +439,8 @@ class Rotation: weights = np.ones(N,dtype='i') for i,(r,n) in enumerate(zip(rotations,weights)): - M = r.M() * n if i == 0 \ - else M + r.M() * n # noqa add (multiples) of this rotation to average noqa + M = r.M * n if i == 0 \ + else M + r.M * n # noqa add (multiples) of this rotation to average noqa eig, vec = np.linalg.eig(M/N) return Rotation.from_quaternion(np.real(vec.T[eig.argmax()]),accept_homomorph = True) @@ -461,7 +465,8 @@ class Rotation: # for compatibility (old names do not follow convention) - asM = M + def asM(self): + return self.M fromQuaternion = from_quaternion fromEulers = from_Eulers asAxisAngle = as_axis_angle diff --git a/python/tests/test_ori_vec.py b/python/tests/test_ori_vec.py index f55dc6d03..a13ad3a03 100644 --- a/python/tests/test_ori_vec.py +++ b/python/tests/test_ori_vec.py @@ -10,18 +10,19 @@ rot1= Rotation.from_random() rot2= Rotation.from_random() rot3= Rotation.from_random() -class TestOrientation_vec: +class TestOrientation_vec: + @pytest.mark.xfail @pytest.mark.parametrize('lattice',Lattice.lattices) def test_equivalentOrientations_vec(self,lattice): ori0=Orientation(rot0,lattice) ori1=Orientation(rot1,lattice) ori2=Orientation(rot2,lattice) ori3=Orientation(rot3,lattice) - + quat=np.array([rot0.as_quaternion(),rot1.as_quaternion(),rot2.as_quaternion(),rot3.as_quaternion()]) rot_vec=Rotation.from_quaternion(quat) ori_vec=Orientation(rot_vec,lattice) - + for s in range(len(ori_vec.lattice.symmetry.symmetryOperations())): assert all(ori_vec.equivalent_vec()[s,0].rotation.as_Eulers() == \ ori0.equivalentOrientations()[s].rotation.as_Eulers()) @@ -31,7 +32,7 @@ class TestOrientation_vec: ori2.equivalentOrientations()[s].rotation.as_Rodrigues()) assert all(ori_vec.equivalent_vec()[s,3].rotation.as_cubochoric() == \ ori3.equivalentOrientations()[s].rotation.as_cubochoric()) - + @pytest.mark.parametrize('lattice',Lattice.lattices) def test_inFZ_vec(self,lattice): ori0=Orientation(rot0,lattice) @@ -41,19 +42,19 @@ class TestOrientation_vec: #ensure 1 of them is in FZ ori4=ori0.reduced() rot4=ori4.rotation - + quat=np.array([rot0.as_quaternion(),rot1.as_quaternion(),\ rot2.as_quaternion(),rot3.as_quaternion(), rot4.as_quaternion()]) rot_vec=Rotation.from_quaternion(quat) ori_vec=Orientation(rot_vec,lattice) - + assert ori_vec.inFZ_vec()[0] == ori0.inFZ() assert ori_vec.inFZ_vec()[1] == ori1.inFZ() assert ori_vec.inFZ_vec()[2] == ori2.inFZ() assert ori_vec.inFZ_vec()[3] == ori3.inFZ() assert ori_vec.inFZ_vec()[4] == ori4.inFZ() - - + + @pytest.mark.parametrize('model',['Bain','KS','GT','GT_prime','NW','Pitsch']) @pytest.mark.parametrize('lattice',['fcc','bcc']) def test_relatedOrientations_vec(self,model,lattice): @@ -61,11 +62,11 @@ class TestOrientation_vec: ori1=Orientation(rot1,lattice) ori2=Orientation(rot2,lattice) ori3=Orientation(rot3,lattice) - + quat=np.array([rot0.as_quaternion(),rot1.as_quaternion(),rot2.as_quaternion(),rot3.as_quaternion()]) rot_vec=Rotation.from_quaternion(quat) ori_vec=Orientation(rot_vec,lattice) - + for s in range(len(ori1.lattice.relationOperations(model)['rotations'])): assert all(ori_vec.relatedOrientations_vec(model)[s,0].rotation.as_Eulers() == \ ori0.relatedOrientations(model)[s].rotation.as_Eulers()) @@ -75,15 +76,4 @@ class TestOrientation_vec: ori2.relatedOrientations(model)[s].rotation.as_Rodrigues()) assert all(ori_vec.relatedOrientations_vec(model)[s,3].rotation.as_cubochoric() == \ ori3.relatedOrientations(model)[s].rotation.as_cubochoric()) - - - - - - - - - - - - \ No newline at end of file + From 13bf7515cecd957d088b57ea878dbd02d9fbc639 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 19 Jun 2020 10:54:13 +0200 Subject: [PATCH 12/40] WIP (broken?): vectorized calculation of IPF color --- python/damask/_lattice.py | 89 ++++++++++++++++++++++++++++++++++- python/damask/_orientation.py | 10 ++++ 2 files changed, 97 insertions(+), 2 deletions(-) diff --git a/python/damask/_lattice.py b/python/damask/_lattice.py index 9769933ef..8582c8d34 100644 --- a/python/damask/_lattice.py +++ b/python/damask/_lattice.py @@ -374,8 +374,93 @@ class Symmetry: else: return inSST -# code derived from https://github.com/ezag/pyeuclid -# suggested reading: http://web.mit.edu/2.998/www/QuaternionReport1.pdf + + def in_SST(self, + vector, + proper = False, + color = False): + """ + Check whether given vector falls into standard stereographic triangle of own symmetry. + + proper considers only vectors with z >= 0, hence uses two neighboring SSTs. + Return inverse pole figure color if requested. + Bases are computed from + + >>> basis = {'cubic' : np.linalg.inv(np.array([[0.,0.,1.], # direction of red + ... [1.,0.,1.]/np.sqrt(2.), # direction of green + ... [1.,1.,1.]/np.sqrt(3.)]).T), # direction of blue + ... 'hexagonal' : np.linalg.inv(np.array([[0.,0.,1.], # direction of red + ... [1.,0.,0.], # direction of green + ... [np.sqrt(3.),1.,0.]/np.sqrt(4.)]).T), # direction of blue + ... 'tetragonal' : np.linalg.inv(np.array([[0.,0.,1.], # direction of red + ... [1.,0.,0.], # direction of green + ... [1.,1.,0.]/np.sqrt(2.)]).T), # direction of blue + ... 'orthorhombic': np.linalg.inv(np.array([[0.,0.,1.], # direction of red + ... [1.,0.,0.], # direction of green + ... [0.,1.,0.]]).T), # direction of blue + ... } + + """ + if self.lattice == 'cubic': + basis = {'improper':np.array([ [-1. , 0. , 1. ], + [ np.sqrt(2.) , -np.sqrt(2.) , 0. ], + [ 0. , np.sqrt(3.) , 0. ] ]), + 'proper':np.array([ [ 0. , -1. , 1. ], + [-np.sqrt(2.) , np.sqrt(2.) , 0. ], + [ np.sqrt(3.) , 0. , 0. ] ]), + } + elif self.lattice == 'hexagonal': + basis = {'improper':np.array([ [ 0. , 0. , 1. ], + [ 1. , -np.sqrt(3.) , 0. ], + [ 0. , 2. , 0. ] ]), + 'proper':np.array([ [ 0. , 0. , 1. ], + [-1. , np.sqrt(3.) , 0. ], + [ np.sqrt(3.) , -1. , 0. ] ]), + } + elif self.lattice == 'tetragonal': + basis = {'improper':np.array([ [ 0. , 0. , 1. ], + [ 1. , -1. , 0. ], + [ 0. , np.sqrt(2.) , 0. ] ]), + 'proper':np.array([ [ 0. , 0. , 1. ], + [-1. , 1. , 0. ], + [ np.sqrt(2.) , 0. , 0. ] ]), + } + elif self.lattice == 'orthorhombic': + basis = {'improper':np.array([ [ 0., 0., 1.], + [ 1., 0., 0.], + [ 0., 1., 0.] ]), + 'proper':np.array([ [ 0., 0., 1.], + [-1., 0., 0.], + [ 0., 1., 0.] ]), + } + else: # direct exit for unspecified symmetry + if color: + return (np.ones_like(vector[...,0],bool),np.zeros_like(vector)) + else: + return np.ones_like(vector[...,0],bool) + + b_p = np.broadcast_to(basis['proper'], vector.shape+(3,)) + if proper: + b_i = np.broadcast_to(basis['improper'],vector.shape+(3,)) + improper = np.all(np.around(np.einsum('...ji,...i',b_i,vector),12)>=0.0,axis=-1,keepdims=True) + theComponents = np.where(np.broadcast_to(improper,vector.shape), + np.around(np.einsum('...ji,...i',b_i,vector),12), + np.around(np.einsum('...ji,...i',b_p,vector),12)) + else: + vector_ = np.block([vector[...,0:2],np.abs(vector[...,2:3])]) # z component projects identical + theComponents = np.around(np.einsum('...ji,...i',b_p,vector_),12) + + in_SST = np.all(theComponents >= 0.0,axis=-1) + + if color: # have to return color array + with np.errstate(invalid='ignore',divide='ignore'): + rgb = (theComponents/np.linalg.norm(theComponents,axis=-1,keepdims=True))**0.5 # smoothen color ramps + rgb = np.minimum(1.,rgb) # limit to maximum intensity + rgb /= np.max(rgb,axis=-1,keepdims=True) # normalize to (HS)V = 1 + rgb[np.invert(np.broadcast_to(in_SST.reshape(vector[...,0].shape+(1,)),vector.shape))] = 0.0 + return (in_SST,rgb) + else: + return in_SST # ****************************************************************************************** diff --git a/python/damask/_orientation.py b/python/damask/_orientation.py index 287e84ff1..794bf9900 100644 --- a/python/damask/_orientation.py +++ b/python/damask/_orientation.py @@ -176,6 +176,16 @@ class Orientation: # make subclass or Rotation? return color + def IPF_color(self,axis): + """TSL color of inverse pole figure for given axis.""" + color = np.zeros(self.rotation.shape) + eq = self.equivalent + pole = eq.rotation @ np.broadcast_to(axis,eq.rotation.shape+(3,)) + in_SST, color = self.lattice.symmetry.in_SST(pole,color=True) + + return color[in_SST] + + @staticmethod def fromAverage(orientations, weights = []): From 262346ff5a7db8019bd83e37a1589d7c324afcbd Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 20 Jun 2020 16:34:19 +0200 Subject: [PATCH 13/40] polishing --- python/damask/_lattice.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/python/damask/_lattice.py b/python/damask/_lattice.py index 8582c8d34..0ebaa883d 100644 --- a/python/damask/_lattice.py +++ b/python/damask/_lattice.py @@ -224,7 +224,7 @@ class Symmetry: symQuats = [ [ 1.0,0.0,0.0,0.0 ], ] - return Rotation(np.array(symQuats)) + return np.array(symQuats) def inFZ(self,rodrigues): @@ -439,7 +439,7 @@ class Symmetry: else: return np.ones_like(vector[...,0],bool) - b_p = np.broadcast_to(basis['proper'], vector.shape+(3,)) + b_p = np.broadcast_to(basis['proper'], vector.shape+(3,)) if proper: b_i = np.broadcast_to(basis['improper'],vector.shape+(3,)) improper = np.all(np.around(np.einsum('...ji,...i',b_i,vector),12)>=0.0,axis=-1,keepdims=True) @@ -457,14 +457,14 @@ class Symmetry: rgb = (theComponents/np.linalg.norm(theComponents,axis=-1,keepdims=True))**0.5 # smoothen color ramps rgb = np.minimum(1.,rgb) # limit to maximum intensity rgb /= np.max(rgb,axis=-1,keepdims=True) # normalize to (HS)V = 1 - rgb[np.invert(np.broadcast_to(in_SST.reshape(vector[...,0].shape+(1,)),vector.shape))] = 0.0 + rgb[~np.broadcast_to(in_SST.reshape(vector[...,0].shape+(1,)),vector.shape)] = 0.0 return (in_SST,rgb) else: return in_SST # ****************************************************************************************** -class Lattice: +class Lattice: # ToDo: Make a subclass of Symmetry! """ Lattice system. From ebdb65d31ff6bbe9327556c7f3bbaf9581d20983 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 20 Jun 2020 16:35:22 +0200 Subject: [PATCH 14/40] standard broadcast_to behavior --- python/damask/_rotation.py | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index 9a0266871..694384de2 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -156,18 +156,18 @@ class Rotation: Rotation to which the misorientation is computed. """ - return other@self.inversed() + return other*self.inversed() def broadcast_to(self,shape): - if isinstance(shape,int): - shape = (shape,) - N = np.prod(shape)//np.prod(self.shape,dtype=int) - - q = np.block([np.repeat(self.quaternion[...,0:1],N).reshape(shape+(1,)), - np.repeat(self.quaternion[...,1:2],N).reshape(shape+(1,)), - np.repeat(self.quaternion[...,2:3],N).reshape(shape+(1,)), - np.repeat(self.quaternion[...,3:4],N).reshape(shape+(1,))]) + if isinstance(shape,int): shape = (shape,) + if self.shape == (): + q = np.broadcast_to(self.quaternion,shape+(4,)) + else: + q = np.block([np.broadcast_to(self.quaternion[...,0:1],shape).reshape(shape+(1,)), + np.broadcast_to(self.quaternion[...,1:2],shape).reshape(shape+(1,)), + np.broadcast_to(self.quaternion[...,2:3],shape).reshape(shape+(1,)), + np.broadcast_to(self.quaternion[...,3:4],shape).reshape(shape+(1,))]) return self.__class__(q) From 4dae3643c9c57127830b04126e279699625a41a1 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 20 Jun 2020 17:15:13 +0200 Subject: [PATCH 15/40] vectorized IPF color working results also uses the vectorized form. Still needs careful checking --- python/damask/_orientation.py | 30 +++++++++++++++++++----------- python/damask/_result.py | 21 ++++++++------------- 2 files changed, 27 insertions(+), 24 deletions(-) diff --git a/python/damask/_orientation.py b/python/damask/_orientation.py index 794bf9900..1704987b4 100644 --- a/python/damask/_orientation.py +++ b/python/damask/_orientation.py @@ -3,7 +3,7 @@ import numpy as np from . import Lattice from . import Rotation -class Orientation: # make subclass or Rotation? +class Orientation: # ToDo: make subclass of lattice and Rotation """ Crystallographic orientation. @@ -105,12 +105,14 @@ class Orientation: # make subclass or Rotation? is added to the left of the rotation array. """ - symmetry_operations = self.lattice.symmetry.symmetry_operations + s = self.lattice.symmetry.symmetry_operations + s = s.reshape(s.shape[:1]+(1,)*len(self.rotation.shape)+(4,)) + s = Rotation(np.broadcast_to(s,s.shape[:1]+self.rotation.quaternion.shape)) - q = np.block([self.rotation.quaternion]*symmetry_operations.shape[0]) - r = Rotation(q.reshape(symmetry_operations.shape+self.rotation.quaternion.shape)) + r = np.broadcast_to(self.rotation.quaternion,s.shape[:1]+self.rotation.quaternion.shape) + r = Rotation(r) - return self.__class__(symmetry_operations.broadcast_to(r.shape)@r,self.lattice) + return self.__class__(s@r,self.lattice) def equivalentOrientations(self,members=[]): @@ -156,10 +158,10 @@ class Orientation: # make subclass or Rotation? """Axis rotated according to orientation (using crystal symmetry to ensure location falls into SST).""" if SST: # pole requested to be within SST for i,o in enumerate(self.equivalentOrientations()): # test all symmetric equivalent quaternions - pole = o.rotation*axis # align crystal direction to axis + pole = o.rotation@axis # align crystal direction to axis if self.lattice.symmetry.inSST(pole,proper): break # found SST version else: - pole = self.rotation*axis # align crystal direction to axis + pole = self.rotation@axis # align crystal direction to axis return (pole,i if SST else 0) @@ -169,7 +171,7 @@ class Orientation: # make subclass or Rotation? color = np.zeros(3,'d') for o in self.equivalentOrientations(): - pole = o.rotation*axis # align crystal direction to axis + pole = o.rotation@axis # align crystal direction to axis inSST,color = self.lattice.symmetry.inSST(pole,color=True) if inSST: break @@ -178,12 +180,18 @@ class Orientation: # make subclass or Rotation? def IPF_color(self,axis): """TSL color of inverse pole figure for given axis.""" - color = np.zeros(self.rotation.shape) eq = self.equivalent - pole = eq.rotation @ np.broadcast_to(axis,eq.rotation.shape+(3,)) + pole = eq.rotation @ np.broadcast_to(axis/np.linalg.norm(axis),eq.rotation.shape+(3,)) in_SST, color = self.lattice.symmetry.in_SST(pole,color=True) - return color[in_SST] + # ignore duplicates (occur for highly symmetric orientations) + found = np.zeros_like(in_SST[1],dtype=bool) + c = np.empty(color.shape[1:]) + for s in range(in_SST.shape[0]): + c = np.where(np.expand_dims(np.logical_and(in_SST[s],~found),-1),color[s],c) + found = np.logical_or(in_SST[s],found) + + return c @staticmethod diff --git a/python/damask/_result.py b/python/damask/_result.py index 56c0eac17..e6dac9370 100644 --- a/python/damask/_result.py +++ b/python/damask/_result.py @@ -10,6 +10,7 @@ from functools import partial import h5py import numpy as np +from numpy.lib import recfunctions as rfn from . import VTK from . import Table @@ -267,7 +268,7 @@ class Result: def rename(self,name_old,name_new): """ Rename datasets. - + Parameters ---------- name_old : str @@ -733,23 +734,17 @@ class Result: @staticmethod def _add_IPFcolor(q,l): - d = np.array(l) - d_unit = d/np.linalg.norm(d) - m = util.scale_to_coprime(d) - colors = np.empty((len(q['data']),3),np.uint8) + m = util.scale_to_coprime(np.array(l)) - lattice = q['meta']['Lattice'] - - for i,qu in enumerate(q['data']): - o = Orientation(np.array([qu['w'],qu['x'],qu['y'],qu['z']]),lattice).reduced() - colors[i] = np.uint8(o.IPFcolor(d_unit)*255) + o = Orientation(Rotation(rfn.structured_to_unstructured(q['data'])), + lattice = q['meta']['Lattice']) return { - 'data': colors, + 'data': np.uint8(o.IPF_color(l)*255), 'label': 'IPFcolor_[{} {} {}]'.format(*m), 'meta' : { - 'Unit': 'RGB (8bit)', - 'Lattice': lattice, + 'Unit': '8-bit RGB', + 'Lattice': q['meta']['Lattice'], 'Description': 'Inverse Pole Figure (IPF) colors along sample direction [{} {} {}]'.format(*m), 'Creator': inspect.stack()[0][3][1:] } From e33895dd35d52e6b13cdf6360498311d7f492bb0 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 21 Jun 2020 10:37:09 +0200 Subject: [PATCH 16/40] [skip ci] better logic --- python/damask/_lattice.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/_lattice.py b/python/damask/_lattice.py index 0ebaa883d..251821868 100644 --- a/python/damask/_lattice.py +++ b/python/damask/_lattice.py @@ -457,7 +457,7 @@ class Symmetry: rgb = (theComponents/np.linalg.norm(theComponents,axis=-1,keepdims=True))**0.5 # smoothen color ramps rgb = np.minimum(1.,rgb) # limit to maximum intensity rgb /= np.max(rgb,axis=-1,keepdims=True) # normalize to (HS)V = 1 - rgb[~np.broadcast_to(in_SST.reshape(vector[...,0].shape+(1,)),vector.shape)] = 0.0 + rgb[np.broadcast_to(~in_SST.reshape(vector[...,0].shape+(1,)),vector.shape)] = 0.0 return (in_SST,rgb) else: return in_SST From 6fa5ae6ebf555e47923f6a5a62bfe1043915e313 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 22 Jun 2020 23:14:58 +0200 Subject: [PATCH 17/40] literature from Karo --- python/damask/_orientation.py | 1 + 1 file changed, 1 insertion(+) diff --git a/python/damask/_orientation.py b/python/damask/_orientation.py index 1704987b4..0e7b68954 100644 --- a/python/damask/_orientation.py +++ b/python/damask/_orientation.py @@ -198,6 +198,7 @@ class Orientation: # ToDo: make subclass of lattice and Rotation def fromAverage(orientations, weights = []): """Create orientation from average of list of orientations.""" + # further read: Orientation distribution analysis in deformed grains, https://doi.org/10.1107/S0021889801003077 if not all(isinstance(item, Orientation) for item in orientations): raise TypeError("Only instances of Orientation can be averaged.") From 352c4e95f1a8ce6cb6321aa10a324613e80544f9 Mon Sep 17 00:00:00 2001 From: "f.basile" Date: Sun, 28 Jun 2020 19:03:06 +0200 Subject: [PATCH 18/40] more vectorized --- python/damask/_orientation.py | 62 ++++++++++++++------------ python/tests/test_ori_vec.py | 83 +++++++++++++++++++++++++---------- 2 files changed, 95 insertions(+), 50 deletions(-) diff --git a/python/damask/_orientation.py b/python/damask/_orientation.py index 0e7b68954..e36713ee9 100644 --- a/python/damask/_orientation.py +++ b/python/damask/_orientation.py @@ -77,27 +77,19 @@ class Orientation: # ToDo: make subclass of lattice and Rotation # self-->other: True, self<--other: False def inFZ_vec(self): - """ - Check if orientations falls into Fundamental Zone. - - self.rotation.as_Rodrigues() working fine - self.rotation.as_Rodrigues(vector=True) doesn't work for several rotations - i apply dirty fix - - """ + """Check if orientations fall into Fundamental Zone.""" if not self.rotation.shape: return self.lattice.symmetry.inFZ(self.rotation.as_Rodrigues(vector=True)) else: return [self.lattice.symmetry.inFZ(\ - Rotation._qu2ro(self.rotation.as_quaternion())[l][...,:3]\ - *Rotation._qu2ro(self.rotation.as_quaternion())[l][...,3])\ - for l in range(self.rotation.shape[0])] + self.rotation.as_Rodrigues(vector=True)[l]) for l in range(self.rotation.shape[0])] + def inFZ(self): return self.lattice.symmetry.inFZ(self.rotation.as_Rodrigues(vector=True)) @property - def equivalent(self): + def equivalent_vec(self): """ Return orientations which are symmetrically equivalent. @@ -105,12 +97,12 @@ class Orientation: # ToDo: make subclass of lattice and Rotation is added to the left of the rotation array. """ - s = self.lattice.symmetry.symmetry_operations - s = s.reshape(s.shape[:1]+(1,)*len(self.rotation.shape)+(4,)) + s = self.lattice.symmetry.symmetry_operations #24 lines (sym) x 4 columns (quat) + s = s.reshape(s.shape[:1]+(1,)*len(self.rotation.shape)+(4,)) #reshape zo (24,1,4) s = Rotation(np.broadcast_to(s,s.shape[:1]+self.rotation.quaternion.shape)) - r = np.broadcast_to(self.rotation.quaternion,s.shape[:1]+self.rotation.quaternion.shape) - r = Rotation(r) + r = np.broadcast_to(self.rotation.quaternion,s.shape[:1]+self.rotation.quaternion.shape) #(24,NumRots,4) + r = Rotation(r) #(24, NumRot) return self.__class__(s@r,self.lattice) @@ -127,14 +119,17 @@ class Orientation: # ToDo: make subclass of lattice and Rotation def relatedOrientations_vec(self,model): """List of orientations related by the given orientation relationship.""" - r = self.lattice.relationOperations(model) - if not self.rotation.shape: - return [self.__class__(o*self.rotation,r['lattice']) for o in r['rotations']] - else: - return np.reshape(\ - [self.__class__(o*Rotation.from_quaternion(self.rotation.as_quaternion()[l])\ - ,r['lattice']) for o in r['rotations'] for l in range(self.rotation.shape[0])] - ,(len(r['rotations']),self.rotation.shape[0])) + h = self.lattice.relationOperations(model) + rot= h['rotations'] + op=np.array([o.as_quaternion() for o in rot]) + + s = op.reshape(op.shape[:1]+(1,)*len(self.rotation.shape)+(4,)) + s = Rotation(np.broadcast_to(s,s.shape[:1]+self.rotation.quaternion.shape)) + + r = np.broadcast_to(self.rotation.quaternion,s.shape[:1]+self.rotation.quaternion.shape) + r = Rotation(r) + + return self.__class__(s@r,h['lattice']) def relatedOrientations(self,model): @@ -142,6 +137,19 @@ class Orientation: # ToDo: make subclass of lattice and Rotation r = self.lattice.relationOperations(model) return [self.__class__(o*self.rotation,r['lattice']) for o in r['rotations']] + @property + def reduced_vec(self): + """Transform orientation to fall into fundamental zone according to symmetry.""" + equi=self.equivalent_vec.rotation #24 x rot x 3(rodrigues vector) + r= 1 if not self.rotation.shape else equi.shape[1] #number of rotations + quat=np.empty( [r , 4]) + for rot in range(r): + for sym in range(equi.shape[0]): + if self.lattice.symmetry.inFZ(equi.as_Rodrigues(vector=True)[sym,rot]) == True: + quat[rot]=equi.as_quaternion()[sym,rot] + return self.__class__(quat,self.lattice) + + def reduced(self): """Transform orientation to fall into fundamental zone according to symmetry.""" @@ -178,9 +186,9 @@ class Orientation: # ToDo: make subclass of lattice and Rotation return color - def IPF_color(self,axis): - """TSL color of inverse pole figure for given axis.""" - eq = self.equivalent + def IPFcolor_vec(self,axis): + """TSL color of inverse pole figure for given axis. Not for hex or triclinic lattices""" + eq = self.equivalent_vec pole = eq.rotation @ np.broadcast_to(axis/np.linalg.norm(axis),eq.rotation.shape+(3,)) in_SST, color = self.lattice.symmetry.in_SST(pole,color=True) diff --git a/python/tests/test_ori_vec.py b/python/tests/test_ori_vec.py index a13ad3a03..daae6c197 100644 --- a/python/tests/test_ori_vec.py +++ b/python/tests/test_ori_vec.py @@ -5,32 +5,36 @@ from damask import Rotation from damask import Orientation from damask import Lattice -rot0= Rotation.from_random() -rot1= Rotation.from_random() -rot2= Rotation.from_random() -rot3= Rotation.from_random() +rot0= Rotation.from_random() ; +rot1= Rotation.from_random() ; +rot2= Rotation.from_random() ; +rot3= Rotation.from_random() ; + +#disorientation + +#fromaverage +#average class TestOrientation_vec: - @pytest.mark.xfail + #@pytest.mark.xfail @pytest.mark.parametrize('lattice',Lattice.lattices) - def test_equivalentOrientations_vec(self,lattice): + def test_equivalent_vec(self,lattice): ori0=Orientation(rot0,lattice) ori1=Orientation(rot1,lattice) ori2=Orientation(rot2,lattice) ori3=Orientation(rot3,lattice) quat=np.array([rot0.as_quaternion(),rot1.as_quaternion(),rot2.as_quaternion(),rot3.as_quaternion()]) - rot_vec=Rotation.from_quaternion(quat) - ori_vec=Orientation(rot_vec,lattice) + ori_vec=Orientation(quat,lattice) for s in range(len(ori_vec.lattice.symmetry.symmetryOperations())): - assert all(ori_vec.equivalent_vec()[s,0].rotation.as_Eulers() == \ + assert all(ori_vec.equivalent_vec.rotation.as_Eulers()[s,0] == \ ori0.equivalentOrientations()[s].rotation.as_Eulers()) - assert all(ori_vec.equivalent_vec()[s,1].rotation.as_quaternion() == \ + assert all(ori_vec.equivalent_vec.rotation.as_quaternion()[s,1] == \ ori1.equivalentOrientations()[s].rotation.as_quaternion()) - assert all(ori_vec.equivalent_vec()[s,2].rotation.as_Rodrigues() == \ + assert all(ori_vec.equivalent_vec.rotation.as_Rodrigues()[s,2] == \ ori2.equivalentOrientations()[s].rotation.as_Rodrigues()) - assert all(ori_vec.equivalent_vec()[s,3].rotation.as_cubochoric() == \ + assert all(ori_vec.equivalent_vec.rotation.as_cubochoric()[s,3] == \ ori3.equivalentOrientations()[s].rotation.as_cubochoric()) @pytest.mark.parametrize('lattice',Lattice.lattices) @@ -39,14 +43,11 @@ class TestOrientation_vec: ori1=Orientation(rot1,lattice) ori2=Orientation(rot2,lattice) ori3=Orientation(rot3,lattice) - #ensure 1 of them is in FZ - ori4=ori0.reduced() - rot4=ori4.rotation + ori4=ori0.reduced() ; rot4=ori4.rotation #ensure 1 of them is in FZ quat=np.array([rot0.as_quaternion(),rot1.as_quaternion(),\ rot2.as_quaternion(),rot3.as_quaternion(), rot4.as_quaternion()]) - rot_vec=Rotation.from_quaternion(quat) - ori_vec=Orientation(rot_vec,lattice) + ori_vec=Orientation(quat,lattice) assert ori_vec.inFZ_vec()[0] == ori0.inFZ() assert ori_vec.inFZ_vec()[1] == ori1.inFZ() @@ -64,16 +65,52 @@ class TestOrientation_vec: ori3=Orientation(rot3,lattice) quat=np.array([rot0.as_quaternion(),rot1.as_quaternion(),rot2.as_quaternion(),rot3.as_quaternion()]) - rot_vec=Rotation.from_quaternion(quat) - ori_vec=Orientation(rot_vec,lattice) + ori_vec=Orientation(quat,lattice) + for s in range(len(ori1.lattice.relationOperations(model)['rotations'])): - assert all(ori_vec.relatedOrientations_vec(model)[s,0].rotation.as_Eulers() == \ + assert all(ori_vec.relatedOrientations_vec(model).rotation.as_Eulers()[s,0] == \ ori0.relatedOrientations(model)[s].rotation.as_Eulers()) - assert all(ori_vec.relatedOrientations_vec(model)[s,1].rotation.as_quaternion() == \ + assert all(ori_vec.relatedOrientations_vec(model).rotation.as_quaternion()[s,1] == \ ori1.relatedOrientations(model)[s].rotation.as_quaternion()) - assert all(ori_vec.relatedOrientations_vec(model)[s,2].rotation.as_Rodrigues() == \ + assert all(ori_vec.relatedOrientations_vec(model).rotation.as_Rodrigues()[s,2] == \ ori2.relatedOrientations(model)[s].rotation.as_Rodrigues()) - assert all(ori_vec.relatedOrientations_vec(model)[s,3].rotation.as_cubochoric() == \ + assert all(ori_vec.relatedOrientations_vec(model).rotation.as_cubochoric()[s,3] == \ ori3.relatedOrientations(model)[s].rotation.as_cubochoric()) + @pytest.mark.parametrize('lattice',Lattice.lattices) + def test_reduced_vec(self,lattice): + ori0=Orientation(rot0,lattice) + ori1=Orientation(rot1,lattice) + ori2=Orientation(rot2,lattice) + ori3=Orientation(rot3,lattice) + #ensure 1 of them is in FZ + ori4=ori0.reduced() + rot4=ori4.rotation + + quat=np.array([rot0.as_quaternion(),rot1.as_quaternion(),\ + rot2.as_quaternion(),rot3.as_quaternion(), rot4.as_quaternion()]) + ori_vec=Orientation(quat,lattice) + + assert all(ori_vec.reduced_vec.rotation.as_Eulers()[0] == ori0.reduced().rotation.as_Eulers() ) + assert all(ori_vec.reduced_vec.rotation.as_quaternion()[1] == ori1.reduced().rotation.as_quaternion() ) + assert all(ori_vec.reduced_vec.rotation.as_Rodrigues()[2] == ori2.reduced().rotation.as_Rodrigues() ) + assert all(ori_vec.reduced_vec.rotation.as_cubochoric()[3] == ori3.reduced().rotation.as_cubochoric() ) + assert all(ori_vec.reduced_vec.rotation.as_axis_angle()[4] == ori4.reduced().rotation.as_axis_angle() ) + + + @pytest.mark.parametrize('lattice',['bcc','fcc','bct']) + def test_IPFcolor_vec(self,lattice): + ori0=Orientation(rot0,lattice) + ori1=Orientation(rot1,lattice) + ori2=Orientation(rot2,lattice) + ori3=Orientation(rot3,lattice) + + quat=np.array([rot0.as_quaternion(),rot1.as_quaternion(),\ + rot2.as_quaternion(),rot3.as_quaternion()]) + ori_vec=Orientation(quat,lattice) + + assert np.allclose( ori_vec.IPFcolor_vec(np.array([0,0,1]))[0],ori0.IPFcolor(np.array([0,0,1]))) + assert np.allclose( ori_vec.IPFcolor_vec(np.array([0,2,1]))[1],ori1.IPFcolor(np.array([0,2,1]))) + assert np.allclose( ori_vec.IPFcolor_vec(np.array([0,3,1]))[2],ori2.IPFcolor(np.array([0,3,1]))) + assert np.allclose( ori_vec.IPFcolor_vec(np.array([4,0,1]))[3],ori3.IPFcolor(np.array([4,0,1]))) From 8484d2e6ccb0a6f85fbb1049010d66f8025f21d2 Mon Sep 17 00:00:00 2001 From: "f.basile" Date: Sun, 28 Jun 2020 19:05:10 +0200 Subject: [PATCH 19/40] fix github stuff --- python/damask/_orientation.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/python/damask/_orientation.py b/python/damask/_orientation.py index e36713ee9..c0952c5fd 100644 --- a/python/damask/_orientation.py +++ b/python/damask/_orientation.py @@ -145,7 +145,7 @@ class Orientation: # ToDo: make subclass of lattice and Rotation quat=np.empty( [r , 4]) for rot in range(r): for sym in range(equi.shape[0]): - if self.lattice.symmetry.inFZ(equi.as_Rodrigues(vector=True)[sym,rot]) == True: + if self.lattice.symmetry.inFZ(equi.as_Rodrigues(vector=True)[sym,rot]) is True: quat[rot]=equi.as_quaternion()[sym,rot] return self.__class__(quat,self.lattice) @@ -187,7 +187,7 @@ class Orientation: # ToDo: make subclass of lattice and Rotation def IPFcolor_vec(self,axis): - """TSL color of inverse pole figure for given axis. Not for hex or triclinic lattices""" + """TSL color of inverse pole figure for given axis. Not for hex or triclinic lattices.""" eq = self.equivalent_vec pole = eq.rotation @ np.broadcast_to(axis/np.linalg.norm(axis),eq.rotation.shape+(3,)) in_SST, color = self.lattice.symmetry.in_SST(pole,color=True) From a99f0164380a45699be7b8af6d9b3e1448560c5e Mon Sep 17 00:00:00 2001 From: "f.basile" Date: Sun, 28 Jun 2020 19:29:52 +0200 Subject: [PATCH 20/40] problem with if value is true / if value == True --- python/damask/_orientation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/_orientation.py b/python/damask/_orientation.py index c0952c5fd..adfa839d1 100644 --- a/python/damask/_orientation.py +++ b/python/damask/_orientation.py @@ -145,7 +145,7 @@ class Orientation: # ToDo: make subclass of lattice and Rotation quat=np.empty( [r , 4]) for rot in range(r): for sym in range(equi.shape[0]): - if self.lattice.symmetry.inFZ(equi.as_Rodrigues(vector=True)[sym,rot]) is True: + if self.lattice.symmetry.inFZ(equi.as_Rodrigues(vector=True)[sym,rot]) == True: quat[rot]=equi.as_quaternion()[sym,rot] return self.__class__(quat,self.lattice) From 4875191ffdda6bca74505d1aefb7e91bfb603ef8 Mon Sep 17 00:00:00 2001 From: "f.basile" Date: Sun, 28 Jun 2020 19:32:22 +0200 Subject: [PATCH 21/40] change if statement so github doesnt complain --- python/damask/_orientation.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/python/damask/_orientation.py b/python/damask/_orientation.py index adfa839d1..cc26acafc 100644 --- a/python/damask/_orientation.py +++ b/python/damask/_orientation.py @@ -145,8 +145,9 @@ class Orientation: # ToDo: make subclass of lattice and Rotation quat=np.empty( [r , 4]) for rot in range(r): for sym in range(equi.shape[0]): - if self.lattice.symmetry.inFZ(equi.as_Rodrigues(vector=True)[sym,rot]) == True: + if self.lattice.symmetry.inFZ(equi.as_Rodrigues(vector=True)[sym,rot]): quat[rot]=equi.as_quaternion()[sym,rot] + break return self.__class__(quat,self.lattice) From d06daec4cb5578f744d83fe91d7e31933947f6f3 Mon Sep 17 00:00:00 2001 From: "f.basile" Date: Mon, 29 Jun 2020 18:25:45 +0200 Subject: [PATCH 22/40] reducec vectorized is improved --- python/damask/_orientation.py | 28 ++++++++++++++++++---------- python/tests/test_ori_vec.py | 8 ++++---- 2 files changed, 22 insertions(+), 14 deletions(-) diff --git a/python/damask/_orientation.py b/python/damask/_orientation.py index cc26acafc..2b6f1f701 100644 --- a/python/damask/_orientation.py +++ b/python/damask/_orientation.py @@ -140,15 +140,23 @@ class Orientation: # ToDo: make subclass of lattice and Rotation @property def reduced_vec(self): """Transform orientation to fall into fundamental zone according to symmetry.""" - equi=self.equivalent_vec.rotation #24 x rot x 3(rodrigues vector) - r= 1 if not self.rotation.shape else equi.shape[1] #number of rotations - quat=np.empty( [r , 4]) - for rot in range(r): - for sym in range(equi.shape[0]): - if self.lattice.symmetry.inFZ(equi.as_Rodrigues(vector=True)[sym,rot]): - quat[rot]=equi.as_quaternion()[sym,rot] - break - return self.__class__(quat,self.lattice) + equi= self.equivalent_vec.rotation #equivalent orientations + r= 1 if not self.rotation.shape else equi.shape[1] #number of rotations + num_equi=equi.shape[0] #number of equivalente orientations + quat= np.reshape( equi.as_quaternion(), (r*num_equi,4) ,order='F') #equivalents are listed in intiuitive order + boolean=Orientation(quat, self.lattice).inFZ_vec() #check which ones are in FZ + if sum(boolean) == r: + return self.__class__(quat[boolean],self.lattice) + + else: + print('More than 1 equivalent orientation has been found for an orientation') + index=np.empty(r, dtype=int) + for l,h in enumerate(range(0,r*num_equi, num_equi)): + index[l]=np.where(boolean[h:h+num_equi])[0][0] + (l*num_equi) #get first index that is true then go check to next orientation + + return self.__class__(quat[index],self.lattice) + + @@ -187,7 +195,7 @@ class Orientation: # ToDo: make subclass of lattice and Rotation return color - def IPFcolor_vec(self,axis): + def IPF_color(self,axis): """TSL color of inverse pole figure for given axis. Not for hex or triclinic lattices.""" eq = self.equivalent_vec pole = eq.rotation @ np.broadcast_to(axis/np.linalg.norm(axis),eq.rotation.shape+(3,)) diff --git a/python/tests/test_ori_vec.py b/python/tests/test_ori_vec.py index daae6c197..6bbe37f63 100644 --- a/python/tests/test_ori_vec.py +++ b/python/tests/test_ori_vec.py @@ -110,7 +110,7 @@ class TestOrientation_vec: rot2.as_quaternion(),rot3.as_quaternion()]) ori_vec=Orientation(quat,lattice) - assert np.allclose( ori_vec.IPFcolor_vec(np.array([0,0,1]))[0],ori0.IPFcolor(np.array([0,0,1]))) - assert np.allclose( ori_vec.IPFcolor_vec(np.array([0,2,1]))[1],ori1.IPFcolor(np.array([0,2,1]))) - assert np.allclose( ori_vec.IPFcolor_vec(np.array([0,3,1]))[2],ori2.IPFcolor(np.array([0,3,1]))) - assert np.allclose( ori_vec.IPFcolor_vec(np.array([4,0,1]))[3],ori3.IPFcolor(np.array([4,0,1]))) + assert np.allclose( ori_vec.IPF_color(np.array([0,0,1]))[0],ori0.IPFcolor(np.array([0,0,1]))) + assert np.allclose( ori_vec.IPF_color(np.array([0,2,1]))[1],ori1.IPFcolor(np.array([0,2,1]))) + assert np.allclose( ori_vec.IPF_color(np.array([0,3,1]))[2],ori2.IPFcolor(np.array([0,3,1]))) + assert np.allclose( ori_vec.IPF_color(np.array([4,0,1]))[3],ori3.IPFcolor(np.array([4,0,1]))) From b8b34080fed48d5ba4e95dc783b624da214fa7f4 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 30 Jun 2020 12:16:47 +0200 Subject: [PATCH 23/40] enable array like slicing/iteration --- python/damask/_rotation.py | 19 ++++++++++++++++++- 1 file changed, 18 insertions(+), 1 deletion(-) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index 694384de2..a8d8a7928 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -52,6 +52,8 @@ class Rotation: Unit quaternion that follows the conventions. Use .from_quaternion to perform a sanity check. """ + if quaternion.shape[-1] != 4: + raise ValueError('Not a quaternion') self.quaternion = quaternion.copy() @@ -60,6 +62,7 @@ class Rotation: return self.quaternion.shape[:-1] + # ToDo: Check difference __copy__ vs __deepcopy__ def __copy__(self): """Copy.""" return self.__class__(self.quaternion) @@ -70,7 +73,7 @@ class Rotation: def __repr__(self): """Orientation displayed as unit quaternion, rotation matrix, and Bunge-Euler angles.""" if self.quaternion.shape != (4,): - return str(self.quaternion) # ToDo: could be nicer ... + return 'Quaternions:\n'+str(self.quaternion) # ToDo: could be nicer ... return '\n'.join([ 'Quaternion: (real={:.3f}, imag=<{:+.3f}, {:+.3f}, {:+.3f}>)'.format(*(self.quaternion)), 'Matrix:\n{}'.format(self.as_matrix()), @@ -78,6 +81,20 @@ class Rotation: ]) + def __len__(self): + return 0 if self.shape == () else len(self.shape) + + + def __getitem__(self,item): + if isinstance(item,tuple) and len(item) >= len(self): + raise IndexError('Too many indices') + return self.__class__(self.quaternion[item]) + + + def __len__(self): + return 0 if self.shape == () else self.shape[0] + + def __matmul__(self, other): """ Rotation of vector, second or fourth order tensor, or rotation object. From ce7018164fd615e09095df5c49e75b53fb29c5cf Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 30 Jun 2020 12:55:08 +0200 Subject: [PATCH 24/40] WIP: more reasonable naming --- python/damask/_lattice.py | 75 ++++++++++++++++++------------------ python/tests/test_Lattice.py | 22 +++++------ 2 files changed, 49 insertions(+), 48 deletions(-) diff --git a/python/damask/_lattice.py b/python/damask/_lattice.py index 251821868..22f1de320 100644 --- a/python/damask/_lattice.py +++ b/python/damask/_lattice.py @@ -5,7 +5,7 @@ from . import Rotation class Symmetry: """ - Symmetry operations for lattice systems. + Symmetry-related operations for crystal systems. References ---------- @@ -13,34 +13,35 @@ class Symmetry: """ - lattices = [None,'orthorhombic','tetragonal','hexagonal','cubic'] + crystal_systems = [None,'orthorhombic','tetragonal','hexagonal','cubic'] - def __init__(self, symmetry = None): + def __init__(self, system = None): """ Symmetry Definition. Parameters ---------- - symmetry : str, optional - label of the crystal system + system : {None,'orthorhombic','tetragonal','hexagonal','cubic'}, optional + Name of the crystal system. Defaults to 'None'. """ - if symmetry is not None and symmetry.lower() not in Symmetry.lattices: - raise KeyError('Symmetry/crystal system "{}" is unknown'.format(symmetry)) + if system is not None and system.lower() not in self.crystal_systems: + raise KeyError(f'Crystal system "{system}" is unknown') - self.lattice = symmetry.lower() if isinstance(symmetry,str) else symmetry + self.system = system.lower() if isinstance(system,str) else system + self.lattice = self.system # for compatibility def __copy__(self): """Copy.""" - return self.__class__(self.lattice) + return self.__class__(self.system) copy = __copy__ def __repr__(self): """Readable string.""" - return '{}'.format(self.lattice) + return '{}'.format(self.system) def __eq__(self, other): @@ -53,7 +54,7 @@ class Symmetry: Symmetry to check for equality. """ - return self.lattice == other.lattice + return self.system == other.system def __neq__(self, other): """ @@ -77,13 +78,13 @@ class Symmetry: Symmetry to check for for order. """ - myOrder = Symmetry.lattices.index(self.lattice) - otherOrder = Symmetry.lattices.index(other.lattice) + myOrder = self.crystal_systems.index(self.system) + otherOrder = self.crystal_systems.index(other.system) return (myOrder > otherOrder) - (myOrder < otherOrder) def symmetryOperations(self,members=[]): """List (or single element) of symmetry operations as rotations.""" - if self.lattice == 'cubic': + if self.system == 'cubic': symQuats = [ [ 1.0, 0.0, 0.0, 0.0 ], [ 0.0, 1.0, 0.0, 0.0 ], @@ -110,7 +111,7 @@ class Symmetry: [-0.5*np.sqrt(2), 0.5*np.sqrt(2), 0.0, 0.0 ], [-0.5*np.sqrt(2),-0.5*np.sqrt(2), 0.0, 0.0 ], ] - elif self.lattice == 'hexagonal': + elif self.system == 'hexagonal': symQuats = [ [ 1.0, 0.0, 0.0, 0.0 ], [-0.5*np.sqrt(3), 0.0, 0.0, -0.5 ], @@ -125,7 +126,7 @@ class Symmetry: [ 0.0, -0.5, -0.5*np.sqrt(3), 0.0 ], [ 0.0, 0.5*np.sqrt(3), 0.5, 0.0 ], ] - elif self.lattice == 'tetragonal': + elif self.system == 'tetragonal': symQuats = [ [ 1.0, 0.0, 0.0, 0.0 ], [ 0.0, 1.0, 0.0, 0.0 ], @@ -136,7 +137,7 @@ class Symmetry: [ 0.5*np.sqrt(2), 0.0, 0.0, 0.5*np.sqrt(2) ], [-0.5*np.sqrt(2), 0.0, 0.0, 0.5*np.sqrt(2) ], ] - elif self.lattice == 'orthorhombic': + elif self.system == 'orthorhombic': symQuats = [ [ 1.0,0.0,0.0,0.0 ], [ 0.0,1.0,0.0,0.0 ], @@ -160,7 +161,7 @@ class Symmetry: @property def symmetry_operations(self): """Symmetry operations as Rotations.""" - if self.lattice == 'cubic': + if self.system == 'cubic': symQuats = [ [ 1.0, 0.0, 0.0, 0.0 ], [ 0.0, 1.0, 0.0, 0.0 ], @@ -187,7 +188,7 @@ class Symmetry: [-0.5*np.sqrt(2), 0.5*np.sqrt(2), 0.0, 0.0 ], [-0.5*np.sqrt(2),-0.5*np.sqrt(2), 0.0, 0.0 ], ] - elif self.lattice == 'hexagonal': + elif self.system == 'hexagonal': symQuats = [ [ 1.0, 0.0, 0.0, 0.0 ], [-0.5*np.sqrt(3), 0.0, 0.0, -0.5 ], @@ -202,7 +203,7 @@ class Symmetry: [ 0.0, -0.5, -0.5*np.sqrt(3), 0.0 ], [ 0.0, 0.5*np.sqrt(3), 0.5, 0.0 ], ] - elif self.lattice == 'tetragonal': + elif self.system == 'tetragonal': symQuats = [ [ 1.0, 0.0, 0.0, 0.0 ], [ 0.0, 1.0, 0.0, 0.0 ], @@ -213,7 +214,7 @@ class Symmetry: [ 0.5*np.sqrt(2), 0.0, 0.0, 0.5*np.sqrt(2) ], [-0.5*np.sqrt(2), 0.0, 0.0, 0.5*np.sqrt(2) ], ] - elif self.lattice == 'orthorhombic': + elif self.system == 'orthorhombic': symQuats = [ [ 1.0,0.0,0.0,0.0 ], [ 0.0,1.0,0.0,0.0 ], @@ -240,21 +241,21 @@ class Symmetry: Rabs = abs(rodrigues) - if self.lattice == 'cubic': + if self.system == 'cubic': return np.sqrt(2.0)-1.0 >= Rabs[0] \ and np.sqrt(2.0)-1.0 >= Rabs[1] \ and np.sqrt(2.0)-1.0 >= Rabs[2] \ and 1.0 >= Rabs[0] + Rabs[1] + Rabs[2] - elif self.lattice == 'hexagonal': + elif self.system == 'hexagonal': return 1.0 >= Rabs[0] and 1.0 >= Rabs[1] and 1.0 >= Rabs[2] \ and 2.0 >= np.sqrt(3)*Rabs[0] + Rabs[1] \ and 2.0 >= np.sqrt(3)*Rabs[1] + Rabs[0] \ and 2.0 >= np.sqrt(3) + Rabs[2] - elif self.lattice == 'tetragonal': + elif self.system == 'tetragonal': return 1.0 >= Rabs[0] and 1.0 >= Rabs[1] \ and np.sqrt(2.0) >= Rabs[0] + Rabs[1] \ and np.sqrt(2.0) >= Rabs[2] + 1.0 - elif self.lattice == 'orthorhombic': + elif self.system == 'orthorhombic': return 1.0 >= Rabs[0] and 1.0 >= Rabs[1] and 1.0 >= Rabs[2] else: return True @@ -275,13 +276,13 @@ class Symmetry: R = rodrigues epsilon = 0.0 - if self.lattice == 'cubic': + if self.system == 'cubic': return R[0] >= R[1]+epsilon and R[1] >= R[2]+epsilon and R[2] >= epsilon - elif self.lattice == 'hexagonal': + elif self.system == 'hexagonal': return R[0] >= np.sqrt(3)*(R[1]-epsilon) and R[1] >= epsilon and R[2] >= epsilon - elif self.lattice == 'tetragonal': + elif self.system == 'tetragonal': return R[0] >= R[1]-epsilon and R[1] >= epsilon and R[2] >= epsilon - elif self.lattice == 'orthorhombic': + elif self.system == 'orthorhombic': return R[0] >= epsilon and R[1] >= epsilon and R[2] >= epsilon else: return True @@ -313,7 +314,7 @@ class Symmetry: ... } """ - if self.lattice == 'cubic': + if self.system == 'cubic': basis = {'improper':np.array([ [-1. , 0. , 1. ], [ np.sqrt(2.) , -np.sqrt(2.) , 0. ], [ 0. , np.sqrt(3.) , 0. ] ]), @@ -321,7 +322,7 @@ class Symmetry: [-np.sqrt(2.) , np.sqrt(2.) , 0. ], [ np.sqrt(3.) , 0. , 0. ] ]), } - elif self.lattice == 'hexagonal': + elif self.system == 'hexagonal': basis = {'improper':np.array([ [ 0. , 0. , 1. ], [ 1. , -np.sqrt(3.) , 0. ], [ 0. , 2. , 0. ] ]), @@ -329,7 +330,7 @@ class Symmetry: [-1. , np.sqrt(3.) , 0. ], [ np.sqrt(3.) , -1. , 0. ] ]), } - elif self.lattice == 'tetragonal': + elif self.system == 'tetragonal': basis = {'improper':np.array([ [ 0. , 0. , 1. ], [ 1. , -1. , 0. ], [ 0. , np.sqrt(2.) , 0. ] ]), @@ -337,7 +338,7 @@ class Symmetry: [-1. , 1. , 0. ], [ np.sqrt(2.) , 0. , 0. ] ]), } - elif self.lattice == 'orthorhombic': + elif self.system == 'orthorhombic': basis = {'improper':np.array([ [ 0., 0., 1.], [ 1., 0., 0.], [ 0., 1., 0.] ]), @@ -401,7 +402,7 @@ class Symmetry: ... } """ - if self.lattice == 'cubic': + if self.system == 'cubic': basis = {'improper':np.array([ [-1. , 0. , 1. ], [ np.sqrt(2.) , -np.sqrt(2.) , 0. ], [ 0. , np.sqrt(3.) , 0. ] ]), @@ -409,7 +410,7 @@ class Symmetry: [-np.sqrt(2.) , np.sqrt(2.) , 0. ], [ np.sqrt(3.) , 0. , 0. ] ]), } - elif self.lattice == 'hexagonal': + elif self.system == 'hexagonal': basis = {'improper':np.array([ [ 0. , 0. , 1. ], [ 1. , -np.sqrt(3.) , 0. ], [ 0. , 2. , 0. ] ]), @@ -417,7 +418,7 @@ class Symmetry: [-1. , np.sqrt(3.) , 0. ], [ np.sqrt(3.) , -1. , 0. ] ]), } - elif self.lattice == 'tetragonal': + elif self.system == 'tetragonal': basis = {'improper':np.array([ [ 0. , 0. , 1. ], [ 1. , -1. , 0. ], [ 0. , np.sqrt(2.) , 0. ] ]), @@ -425,7 +426,7 @@ class Symmetry: [-1. , 1. , 0. ], [ np.sqrt(2.) , 0. , 0. ] ]), } - elif self.lattice == 'orthorhombic': + elif self.system == 'orthorhombic': basis = {'improper':np.array([ [ 0., 0., 1.], [ 1., 0., 0.], [ 0., 1., 0.] ]), diff --git a/python/tests/test_Lattice.py b/python/tests/test_Lattice.py index fb8ce8ea5..3d9c35706 100644 --- a/python/tests/test_Lattice.py +++ b/python/tests/test_Lattice.py @@ -13,26 +13,26 @@ class TestSymmetry: s = Symmetry(invalid_symmetry) # noqa def test_equal(self): - symmetry = random.choice(Symmetry.lattices) + symmetry = random.choice(Symmetry.crystal_systems) print(symmetry) assert Symmetry(symmetry) == Symmetry(symmetry) def test_not_equal(self): - symmetries = random.sample(Symmetry.lattices,k=2) + symmetries = random.sample(Symmetry.crystal_systems,k=2) assert Symmetry(symmetries[0]) != Symmetry(symmetries[1]) - @pytest.mark.parametrize('lattice',Symmetry.lattices) - def test_inFZ(self,lattice): - assert Symmetry(lattice).inFZ(np.zeros(3)) + @pytest.mark.parametrize('system',Symmetry.crystal_systems) + def test_inFZ(self,system): + assert Symmetry(system).inFZ(np.zeros(3)) - @pytest.mark.parametrize('lattice',Symmetry.lattices) - def test_inDisorientationSST(self,lattice): - assert Symmetry(lattice).inDisorientationSST(np.zeros(3)) + @pytest.mark.parametrize('system',Symmetry.crystal_systems) + def test_inDisorientationSST(self,system): + assert Symmetry(system).inDisorientationSST(np.zeros(3)) - @pytest.mark.parametrize('lattice',Symmetry.lattices) + @pytest.mark.parametrize('system',Symmetry.crystal_systems) @pytest.mark.parametrize('proper',[True,False]) - def test_inSST(self,lattice,proper): - assert Symmetry(lattice).inSST(np.zeros(3),proper) + def test_inSST(self,system,proper): + assert Symmetry(system).inSST(np.zeros(3),proper) @pytest.mark.parametrize('function',['inFZ','inDisorientationSST']) def test_invalid_argument(self,function): From 9d94b521ad5feba47ec78fa2816a3ec5f5b27524 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 30 Jun 2020 13:31:58 +0200 Subject: [PATCH 25/40] polishing --- python/damask/_orientation.py | 12 ++++++------ python/damask/_rotation.py | 3 ++- python/tests/test_ori_vec.py | 8 ++++---- 3 files changed, 12 insertions(+), 11 deletions(-) diff --git a/python/damask/_orientation.py b/python/damask/_orientation.py index 2b6f1f701..6333bbfa0 100644 --- a/python/damask/_orientation.py +++ b/python/damask/_orientation.py @@ -125,10 +125,10 @@ class Orientation: # ToDo: make subclass of lattice and Rotation s = op.reshape(op.shape[:1]+(1,)*len(self.rotation.shape)+(4,)) s = Rotation(np.broadcast_to(s,s.shape[:1]+self.rotation.quaternion.shape)) - + r = np.broadcast_to(self.rotation.quaternion,s.shape[:1]+self.rotation.quaternion.shape) r = Rotation(r) - + return self.__class__(s@r,h['lattice']) @@ -142,7 +142,7 @@ class Orientation: # ToDo: make subclass of lattice and Rotation """Transform orientation to fall into fundamental zone according to symmetry.""" equi= self.equivalent_vec.rotation #equivalent orientations r= 1 if not self.rotation.shape else equi.shape[1] #number of rotations - num_equi=equi.shape[0] #number of equivalente orientations + num_equi=equi.shape[0] #number of equivalente orientations quat= np.reshape( equi.as_quaternion(), (r*num_equi,4) ,order='F') #equivalents are listed in intiuitive order boolean=Orientation(quat, self.lattice).inFZ_vec() #check which ones are in FZ if sum(boolean) == r: @@ -150,12 +150,12 @@ class Orientation: # ToDo: make subclass of lattice and Rotation else: print('More than 1 equivalent orientation has been found for an orientation') - index=np.empty(r, dtype=int) + index=np.empty(r, dtype=int) for l,h in enumerate(range(0,r*num_equi, num_equi)): index[l]=np.where(boolean[h:h+num_equi])[0][0] + (l*num_equi) #get first index that is true then go check to next orientation - + return self.__class__(quat[index],self.lattice) - + diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index a8d8a7928..686b144ae 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -267,7 +267,8 @@ class Rotation: """ ro = Rotation._qu2ro(self.quaternion) - return ro[...,:3]*ro[...,3:] if vector else ro + with np.errstate(invalid='ignore'): + return ro[...,:3]*ro[...,3:] if vector else ro def as_homochoric(self): """Homochoric vector: (h_1, h_2, h_3).""" diff --git a/python/tests/test_ori_vec.py b/python/tests/test_ori_vec.py index 6bbe37f63..ffe2923ff 100644 --- a/python/tests/test_ori_vec.py +++ b/python/tests/test_ori_vec.py @@ -5,10 +5,10 @@ from damask import Rotation from damask import Orientation from damask import Lattice -rot0= Rotation.from_random() ; -rot1= Rotation.from_random() ; -rot2= Rotation.from_random() ; -rot3= Rotation.from_random() ; +rot0= Rotation.from_random() +rot1= Rotation.from_random() +rot2= Rotation.from_random() +rot3= Rotation.from_random() #disorientation From c86e3e292c2b189eceb1cf49ce1dd8af7a9542ee Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 30 Jun 2020 13:55:09 +0200 Subject: [PATCH 26/40] WIP: cleaning namespace --- python/damask/_lattice.py | 17 +++++++++++------ python/damask/_orientation.py | 31 ++++++++++++++++++------------- python/tests/test_ori_vec.py | 8 ++++---- 3 files changed, 33 insertions(+), 23 deletions(-) diff --git a/python/damask/_lattice.py b/python/damask/_lattice.py index 22f1de320..ee718adcd 100644 --- a/python/damask/_lattice.py +++ b/python/damask/_lattice.py @@ -479,11 +479,11 @@ class Lattice: # ToDo: Make a subclass of Symmetry! """ lattices = { - 'triclinic':{'symmetry':None}, - 'bct':{'symmetry':'tetragonal'}, - 'hex':{'symmetry':'hexagonal'}, - 'fcc':{'symmetry':'cubic','c/a':1.0}, - 'bcc':{'symmetry':'cubic','c/a':1.0}, + 'triclinic':{'system':None}, + 'bct': {'system':'tetragonal'}, + 'hex': {'system':'hexagonal'}, + 'fcc': {'system':'cubic','c/a':1.0}, + 'bcc': {'system':'cubic','c/a':1.0}, } @@ -498,7 +498,12 @@ class Lattice: # ToDo: Make a subclass of Symmetry! """ self.lattice = lattice - self.symmetry = Symmetry(self.lattices[lattice]['symmetry']) + self.symmetry = Symmetry(self.lattices[lattice]['system']) + + # transition to subclass + self.system = self.symmetry.system + self.in_SST = self.symmetry.in_SST + self.inFZ = self.symmetry.inFZ def __repr__(self): diff --git a/python/damask/_orientation.py b/python/damask/_orientation.py index 6333bbfa0..753da4d4e 100644 --- a/python/damask/_orientation.py +++ b/python/damask/_orientation.py @@ -39,6 +39,11 @@ class Orientation: # ToDo: make subclass of lattice and Rotation else: self.rotation = Rotation.from_quaternion(rotation) # assume quaternion + def __getitem__(self,item): + if isinstance(item,tuple) and len(item) >= len(self): + raise IndexError('Too many indices') + return self.__class__(self.rotation[item],self.lattice) + def disorientation(self, other, @@ -66,7 +71,7 @@ class Orientation: # ToDo: make subclass of lattice and Rotation r = b*aInv for k in range(2): r.inverse() - breaker = self.lattice.symmetry.inFZ(r.as_Rodrigues(vector=True)) \ + breaker = self.lattice.inFZ(r.as_Rodrigues(vector=True)) \ and (not SST or other.lattice.symmetry.inDisorientationSST(r.as_Rodrigues(vector=True))) if breaker: break if breaker: break @@ -79,17 +84,17 @@ class Orientation: # ToDo: make subclass of lattice and Rotation def inFZ_vec(self): """Check if orientations fall into Fundamental Zone.""" if not self.rotation.shape: - return self.lattice.symmetry.inFZ(self.rotation.as_Rodrigues(vector=True)) + return self.lattice.inFZ(self.rotation.as_Rodrigues(vector=True)) else: - return [self.lattice.symmetry.inFZ(\ + return [self.lattice.inFZ(\ self.rotation.as_Rodrigues(vector=True)[l]) for l in range(self.rotation.shape[0])] def inFZ(self): - return self.lattice.symmetry.inFZ(self.rotation.as_Rodrigues(vector=True)) + return self.lattice.inFZ(self.rotation.as_Rodrigues(vector=True)) @property - def equivalent_vec(self): + def equivalent(self): """ Return orientations which are symmetrically equivalent. @@ -97,12 +102,12 @@ class Orientation: # ToDo: make subclass of lattice and Rotation is added to the left of the rotation array. """ - s = self.lattice.symmetry.symmetry_operations #24 lines (sym) x 4 columns (quat) - s = s.reshape(s.shape[:1]+(1,)*len(self.rotation.shape)+(4,)) #reshape zo (24,1,4) + s = self.lattice.symmetry.symmetry_operations + s = s.reshape(s.shape[:1]+(1,)*len(self.rotation.shape)+(4,)) s = Rotation(np.broadcast_to(s,s.shape[:1]+self.rotation.quaternion.shape)) - r = np.broadcast_to(self.rotation.quaternion,s.shape[:1]+self.rotation.quaternion.shape) #(24,NumRots,4) - r = Rotation(r) #(24, NumRot) + r = np.broadcast_to(self.rotation.quaternion,s.shape[:1]+self.rotation.quaternion.shape) + r = Rotation(r) return self.__class__(s@r,self.lattice) @@ -140,7 +145,7 @@ class Orientation: # ToDo: make subclass of lattice and Rotation @property def reduced_vec(self): """Transform orientation to fall into fundamental zone according to symmetry.""" - equi= self.equivalent_vec.rotation #equivalent orientations + equi= self.equivalent.rotation #equivalent orientations r= 1 if not self.rotation.shape else equi.shape[1] #number of rotations num_equi=equi.shape[0] #number of equivalente orientations quat= np.reshape( equi.as_quaternion(), (r*num_equi,4) ,order='F') #equivalents are listed in intiuitive order @@ -163,7 +168,7 @@ class Orientation: # ToDo: make subclass of lattice and Rotation def reduced(self): """Transform orientation to fall into fundamental zone according to symmetry.""" for me in self.equivalentOrientations(): - if self.lattice.symmetry.inFZ(me.rotation.as_Rodrigues(vector=True)): break + if self.lattice.inFZ(me.rotation.as_Rodrigues(vector=True)): break return self.__class__(me.rotation,self.lattice) @@ -197,9 +202,9 @@ class Orientation: # ToDo: make subclass of lattice and Rotation def IPF_color(self,axis): """TSL color of inverse pole figure for given axis. Not for hex or triclinic lattices.""" - eq = self.equivalent_vec + eq = self.equivalent pole = eq.rotation @ np.broadcast_to(axis/np.linalg.norm(axis),eq.rotation.shape+(3,)) - in_SST, color = self.lattice.symmetry.in_SST(pole,color=True) + in_SST, color = self.lattice.in_SST(pole,color=True) # ignore duplicates (occur for highly symmetric orientations) found = np.zeros_like(in_SST[1],dtype=bool) diff --git a/python/tests/test_ori_vec.py b/python/tests/test_ori_vec.py index ffe2923ff..ed64306c4 100644 --- a/python/tests/test_ori_vec.py +++ b/python/tests/test_ori_vec.py @@ -28,13 +28,13 @@ class TestOrientation_vec: ori_vec=Orientation(quat,lattice) for s in range(len(ori_vec.lattice.symmetry.symmetryOperations())): - assert all(ori_vec.equivalent_vec.rotation.as_Eulers()[s,0] == \ + assert all(ori_vec.equivalent.rotation.as_Eulers()[s,0] == \ ori0.equivalentOrientations()[s].rotation.as_Eulers()) - assert all(ori_vec.equivalent_vec.rotation.as_quaternion()[s,1] == \ + assert all(ori_vec.equivalent.rotation.as_quaternion()[s,1] == \ ori1.equivalentOrientations()[s].rotation.as_quaternion()) - assert all(ori_vec.equivalent_vec.rotation.as_Rodrigues()[s,2] == \ + assert all(ori_vec.equivalent.rotation.as_Rodrigues()[s,2] == \ ori2.equivalentOrientations()[s].rotation.as_Rodrigues()) - assert all(ori_vec.equivalent_vec.rotation.as_cubochoric()[s,3] == \ + assert all(ori_vec.equivalent.rotation.as_cubochoric()[s,3] == \ ori3.equivalentOrientations()[s].rotation.as_cubochoric()) @pytest.mark.parametrize('lattice',Lattice.lattices) From be21d1289d2d419c8132c9ee9a1320aaf590e94f Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 30 Jun 2020 16:33:04 +0200 Subject: [PATCH 27/40] using slicing method --- python/damask/_orientation.py | 2 +- python/tests/test_Orientation.py | 2 +- python/tests/test_ori_vec.py | 8 ++++---- 3 files changed, 6 insertions(+), 6 deletions(-) diff --git a/python/damask/_orientation.py b/python/damask/_orientation.py index 753da4d4e..7d0f11c49 100644 --- a/python/damask/_orientation.py +++ b/python/damask/_orientation.py @@ -192,7 +192,7 @@ class Orientation: # ToDo: make subclass of lattice and Rotation """TSL color of inverse pole figure for given axis.""" color = np.zeros(3,'d') - for o in self.equivalentOrientations(): + for o in self.equivalent: pole = o.rotation@axis # align crystal direction to axis inSST,color = self.lattice.symmetry.inSST(pole,color=True) if inSST: break diff --git a/python/tests/test_Orientation.py b/python/tests/test_Orientation.py index a8b8afdac..277fa2a4b 100644 --- a/python/tests/test_Orientation.py +++ b/python/tests/test_Orientation.py @@ -34,7 +34,7 @@ class TestOrientation: for rot in [Rotation.from_random() for r in range(n//100)]: R = damask.Orientation(rot,lattice) color = R.IPFcolor(direction) - for equivalent in R.equivalentOrientations(): + for equivalent in R.equivalent: assert np.allclose(color,R.IPFcolor(direction)) @pytest.mark.parametrize('model',['Bain','KS','GT','GT_prime','NW','Pitsch']) diff --git a/python/tests/test_ori_vec.py b/python/tests/test_ori_vec.py index ed64306c4..ff5fe80bc 100644 --- a/python/tests/test_ori_vec.py +++ b/python/tests/test_ori_vec.py @@ -29,13 +29,13 @@ class TestOrientation_vec: for s in range(len(ori_vec.lattice.symmetry.symmetryOperations())): assert all(ori_vec.equivalent.rotation.as_Eulers()[s,0] == \ - ori0.equivalentOrientations()[s].rotation.as_Eulers()) + ori0.equivalent[s].rotation.as_Eulers()) assert all(ori_vec.equivalent.rotation.as_quaternion()[s,1] == \ - ori1.equivalentOrientations()[s].rotation.as_quaternion()) + ori1.equivalent[s].rotation.as_quaternion()) assert all(ori_vec.equivalent.rotation.as_Rodrigues()[s,2] == \ - ori2.equivalentOrientations()[s].rotation.as_Rodrigues()) + ori2.equivalent[s].rotation.as_Rodrigues()) assert all(ori_vec.equivalent.rotation.as_cubochoric()[s,3] == \ - ori3.equivalentOrientations()[s].rotation.as_cubochoric()) + ori3.equivalent[s].rotation.as_cubochoric()) @pytest.mark.parametrize('lattice',Lattice.lattices) def test_inFZ_vec(self,lattice): From 3d6afff27a8a66cdd1c1a0b0a3921a8acb519c24 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 30 Jun 2020 17:30:29 +0200 Subject: [PATCH 28/40] clearer name --- python/tests/test_Rotation.py | 102 +++++++++++++++++----------------- 1 file changed, 51 insertions(+), 51 deletions(-) diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py index 5a1cd145d..b087cc774 100644 --- a/python/tests/test_Rotation.py +++ b/python/tests/test_Rotation.py @@ -6,15 +6,15 @@ import numpy as np from damask import Rotation from damask import _rotation - - -n = 1100 +n = 1000 atol=1.e-4 -scatter=1.e-2 + @pytest.fixture -def default(): +def set_of_rotations(): """A set of n random rotations.""" + n = 1100 + scatter=1.e-2 specials = np.array([ [1.0, 0.0, 0.0, 0.0], #---------------------- @@ -567,9 +567,9 @@ class TestRotation: (Rotation._qu2ro,Rotation._ro2qu), (Rotation._qu2ho,Rotation._ho2qu), (Rotation._qu2cu,Rotation._cu2qu)]) - def test_quaternion_internal(self,default,forward,backward): + def test_quaternion_internal(self,set_of_rotations,forward,backward): """Ensure invariance of conversion from quaternion and back.""" - for rot in default: + for rot in set_of_rotations: m = rot.as_quaternion() o = backward(forward(m)) ok = np.allclose(m,o,atol=atol) @@ -584,9 +584,9 @@ class TestRotation: (Rotation._om2ro,Rotation._ro2om), (Rotation._om2ho,Rotation._ho2om), (Rotation._om2cu,Rotation._cu2om)]) - def test_matrix_internal(self,default,forward,backward): + def test_matrix_internal(self,set_of_rotations,forward,backward): """Ensure invariance of conversion from rotation matrix and back.""" - for rot in default: + for rot in set_of_rotations: m = rot.as_matrix() o = backward(forward(m)) ok = np.allclose(m,o,atol=atol) @@ -599,9 +599,9 @@ class TestRotation: (Rotation._eu2ro,Rotation._ro2eu), (Rotation._eu2ho,Rotation._ho2eu), (Rotation._eu2cu,Rotation._cu2eu)]) - def test_Eulers_internal(self,default,forward,backward): + def test_Eulers_internal(self,set_of_rotations,forward,backward): """Ensure invariance of conversion from Euler angles and back.""" - for rot in default: + for rot in set_of_rotations: m = rot.as_Eulers() o = backward(forward(m)) u = np.array([np.pi*2,np.pi,np.pi*2]) @@ -619,9 +619,9 @@ class TestRotation: (Rotation._ax2ro,Rotation._ro2ax), (Rotation._ax2ho,Rotation._ho2ax), (Rotation._ax2cu,Rotation._cu2ax)]) - def test_axis_angle_internal(self,default,forward,backward): + def test_axis_angle_internal(self,set_of_rotations,forward,backward): """Ensure invariance of conversion from axis angle angles pair and back.""" - for rot in default: + for rot in set_of_rotations: m = rot.as_axis_angle() o = backward(forward(m)) ok = np.allclose(m,o,atol=atol) @@ -636,10 +636,10 @@ class TestRotation: (Rotation._ro2ax,Rotation._ax2ro), (Rotation._ro2ho,Rotation._ho2ro), (Rotation._ro2cu,Rotation._cu2ro)]) - def test_Rodrigues_internal(self,default,forward,backward): + def test_Rodrigues_internal(self,set_of_rotations,forward,backward): """Ensure invariance of conversion from Rodrigues-Frank vector and back.""" cutoff = np.tan(np.pi*.5*(1.-1e-4)) - for rot in default: + for rot in set_of_rotations: m = rot.as_Rodrigues() o = backward(forward(m)) ok = np.allclose(np.clip(m,None,cutoff),np.clip(o,None,cutoff),atol=atol) @@ -653,9 +653,9 @@ class TestRotation: (Rotation._ho2ax,Rotation._ax2ho), (Rotation._ho2ro,Rotation._ro2ho), (Rotation._ho2cu,Rotation._cu2ho)]) - def test_homochoric_internal(self,default,forward,backward): + def test_homochoric_internal(self,set_of_rotations,forward,backward): """Ensure invariance of conversion from homochoric vector and back.""" - for rot in default: + for rot in set_of_rotations: m = rot.as_homochoric() o = backward(forward(m)) ok = np.allclose(m,o,atol=atol) @@ -668,9 +668,9 @@ class TestRotation: (Rotation._cu2ax,Rotation._ax2cu), (Rotation._cu2ro,Rotation._ro2cu), (Rotation._cu2ho,Rotation._ho2cu)]) - def test_cubochoric_internal(self,default,forward,backward): + def test_cubochoric_internal(self,set_of_rotations,forward,backward): """Ensure invariance of conversion from cubochoric vector and back.""" - for rot in default: + for rot in set_of_rotations: m = rot.as_cubochoric() o = backward(forward(m)) ok = np.allclose(m,o,atol=atol) @@ -684,9 +684,9 @@ class TestRotation: (Rotation._qu2ax,qu2ax), (Rotation._qu2ro,qu2ro), (Rotation._qu2ho,qu2ho)]) - def test_quaternion_vectorization(self,default,vectorized,single): + def test_quaternion_vectorization(self,set_of_rotations,vectorized,single): """Check vectorized implementation for quaternion against single point calculation.""" - qu = np.array([rot.as_quaternion() for rot in default]) + qu = np.array([rot.as_quaternion() for rot in set_of_rotations]) vectorized(qu.reshape(qu.shape[0]//2,-1,4)) co = vectorized(qu) for q,c in zip(qu,co): @@ -697,9 +697,9 @@ class TestRotation: @pytest.mark.parametrize('vectorized, single',[(Rotation._om2qu,om2qu), (Rotation._om2eu,om2eu), (Rotation._om2ax,om2ax)]) - def test_matrix_vectorization(self,default,vectorized,single): + def test_matrix_vectorization(self,set_of_rotations,vectorized,single): """Check vectorized implementation for rotation matrix against single point calculation.""" - om = np.array([rot.as_matrix() for rot in default]) + om = np.array([rot.as_matrix() for rot in set_of_rotations]) vectorized(om.reshape(om.shape[0]//2,-1,3,3)) co = vectorized(om) for o,c in zip(om,co): @@ -710,9 +710,9 @@ class TestRotation: (Rotation._eu2om,eu2om), (Rotation._eu2ax,eu2ax), (Rotation._eu2ro,eu2ro)]) - def test_Eulers_vectorization(self,default,vectorized,single): + def test_Eulers_vectorization(self,set_of_rotations,vectorized,single): """Check vectorized implementation for Euler angles against single point calculation.""" - eu = np.array([rot.as_Eulers() for rot in default]) + eu = np.array([rot.as_Eulers() for rot in set_of_rotations]) vectorized(eu.reshape(eu.shape[0]//2,-1,3)) co = vectorized(eu) for e,c in zip(eu,co): @@ -723,9 +723,9 @@ class TestRotation: (Rotation._ax2om,ax2om), (Rotation._ax2ro,ax2ro), (Rotation._ax2ho,ax2ho)]) - def test_axis_angle_vectorization(self,default,vectorized,single): + def test_axis_angle_vectorization(self,set_of_rotations,vectorized,single): """Check vectorized implementation for axis angle pair against single point calculation.""" - ax = np.array([rot.as_axis_angle() for rot in default]) + ax = np.array([rot.as_axis_angle() for rot in set_of_rotations]) vectorized(ax.reshape(ax.shape[0]//2,-1,4)) co = vectorized(ax) for a,c in zip(ax,co): @@ -735,9 +735,9 @@ class TestRotation: @pytest.mark.parametrize('vectorized, single',[(Rotation._ro2ax,ro2ax), (Rotation._ro2ho,ro2ho)]) - def test_Rodrigues_vectorization(self,default,vectorized,single): + def test_Rodrigues_vectorization(self,set_of_rotations,vectorized,single): """Check vectorized implementation for Rodrigues-Frank vector against single point calculation.""" - ro = np.array([rot.as_Rodrigues() for rot in default]) + ro = np.array([rot.as_Rodrigues() for rot in set_of_rotations]) vectorized(ro.reshape(ro.shape[0]//2,-1,4)) co = vectorized(ro) for r,c in zip(ro,co): @@ -746,9 +746,9 @@ class TestRotation: @pytest.mark.parametrize('vectorized, single',[(Rotation._ho2ax,ho2ax), (Rotation._ho2cu,ho2cu)]) - def test_homochoric_vectorization(self,default,vectorized,single): + def test_homochoric_vectorization(self,set_of_rotations,vectorized,single): """Check vectorized implementation for homochoric vector against single point calculation.""" - ho = np.array([rot.as_homochoric() for rot in default]) + ho = np.array([rot.as_homochoric() for rot in set_of_rotations]) vectorized(ho.reshape(ho.shape[0]//2,-1,3)) co = vectorized(ho) for h,c in zip(ho,co): @@ -756,9 +756,9 @@ class TestRotation: assert np.allclose(single(h),c) and np.allclose(single(h),vectorized(h)) @pytest.mark.parametrize('vectorized, single',[(Rotation._cu2ho,cu2ho)]) - def test_cubochoric_vectorization(self,default,vectorized,single): + def test_cubochoric_vectorization(self,set_of_rotations,vectorized,single): """Check vectorized implementation for cubochoric vector against single point calculation.""" - cu = np.array([rot.as_cubochoric() for rot in default]) + cu = np.array([rot.as_cubochoric() for rot in set_of_rotations]) vectorized(cu.reshape(cu.shape[0]//2,-1,3)) co = vectorized(cu) for u,c in zip(cu,co): @@ -766,8 +766,8 @@ class TestRotation: assert np.allclose(single(u),c) and np.allclose(single(u),vectorized(u)) @pytest.mark.parametrize('degrees',[True,False]) - def test_Eulers(self,default,degrees): - for rot in default: + def test_Eulers(self,set_of_rotations,degrees): + for rot in set_of_rotations: m = rot.as_quaternion() o = Rotation.from_Eulers(rot.as_Eulers(degrees),degrees).as_quaternion() ok = np.allclose(m,o,atol=atol) @@ -779,9 +779,9 @@ class TestRotation: @pytest.mark.parametrize('P',[1,-1]) @pytest.mark.parametrize('normalise',[True,False]) @pytest.mark.parametrize('degrees',[True,False]) - def test_axis_angle(self,default,degrees,normalise,P): + def test_axis_angle(self,set_of_rotations,degrees,normalise,P): c = np.array([P*-1,P*-1,P*-1,1.]) - for rot in default: + for rot in set_of_rotations: m = rot.as_Eulers() o = Rotation.from_axis_angle(rot.as_axis_angle(degrees)*c,degrees,normalise,P).as_Eulers() u = np.array([np.pi*2,np.pi,np.pi*2]) @@ -793,8 +793,8 @@ class TestRotation: print(m,o,rot.as_quaternion()) assert ok and (np.zeros(3)-1.e-9 <= o).all() and (o <= np.array([np.pi*2.,np.pi,np.pi*2.])+1.e-9).all() - def test_matrix(self,default): - for rot in default: + def test_matrix(self,set_of_rotations): + for rot in set_of_rotations: m = rot.as_axis_angle() o = Rotation.from_axis_angle(rot.as_axis_angle()).as_axis_angle() ok = np.allclose(m,o,atol=atol) @@ -805,9 +805,9 @@ class TestRotation: @pytest.mark.parametrize('P',[1,-1]) @pytest.mark.parametrize('normalise',[True,False]) - def test_Rodrigues(self,default,normalise,P): + def test_Rodrigues(self,set_of_rotations,normalise,P): c = np.array([P*-1,P*-1,P*-1,1.]) - for rot in default: + for rot in set_of_rotations: m = rot.as_matrix() o = Rotation.from_Rodrigues(rot.as_Rodrigues()*c,normalise,P).as_matrix() ok = np.allclose(m,o,atol=atol) @@ -815,9 +815,9 @@ class TestRotation: assert ok and np.isclose(np.linalg.det(o),1.0) @pytest.mark.parametrize('P',[1,-1]) - def test_homochoric(self,default,P): + def test_homochoric(self,set_of_rotations,P): cutoff = np.tan(np.pi*.5*(1.-1e-4)) - for rot in default: + for rot in set_of_rotations: m = rot.as_Rodrigues() o = Rotation.from_homochoric(rot.as_homochoric()*P*-1,P).as_Rodrigues() ok = np.allclose(np.clip(m,None,cutoff),np.clip(o,None,cutoff),atol=atol) @@ -826,8 +826,8 @@ class TestRotation: assert ok and np.isclose(np.linalg.norm(o[:3]),1.0) @pytest.mark.parametrize('P',[1,-1]) - def test_cubochoric(self,default,P): - for rot in default: + def test_cubochoric(self,set_of_rotations,P): + for rot in set_of_rotations: m = rot.as_homochoric() o = Rotation.from_cubochoric(rot.as_cubochoric()*P*-1,P).as_homochoric() ok = np.allclose(m,o,atol=atol) @@ -836,9 +836,9 @@ class TestRotation: @pytest.mark.parametrize('P',[1,-1]) @pytest.mark.parametrize('accept_homomorph',[True,False]) - def test_quaternion(self,default,P,accept_homomorph): + 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) - for rot in default: + for rot in set_of_rotations: m = rot.as_cubochoric() o = Rotation.from_quaternion(rot.as_quaternion()*c,accept_homomorph,P).as_cubochoric() ok = np.allclose(m,o,atol=atol) @@ -848,8 +848,8 @@ class TestRotation: assert ok and o.max() < np.pi**(2./3.)*0.5+1.e-9 @pytest.mark.parametrize('reciprocal',[True,False]) - def test_basis(self,default,reciprocal): - for rot in default: + def test_basis(self,set_of_rotations,reciprocal): + for rot in set_of_rotations: om = rot.as_matrix() + 0.1*np.eye(3) rot = Rotation.from_basis(om,False,reciprocal=reciprocal) assert np.isclose(np.linalg.det(rot.as_matrix()),1.0) @@ -909,8 +909,8 @@ class TestRotation: @pytest.mark.parametrize('data',[np.random.rand(5,3), np.random.rand(5,3,3), np.random.rand(5,3,3,3,3)]) - def test_rotate_vectorization(self,default,data): - for rot in default: + def test_rotate_vectorization(self,set_of_rotations,data): + for rot in set_of_rotations: v = rot.broadcast_to((5,)) @ data for i in range(data.shape[0]): print(i-data[i]) From 6e27a140f66796e47072897c9d0659323127342f Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 30 Jun 2020 17:35:52 +0200 Subject: [PATCH 29/40] better split --- python/tests/conftest.py | 79 +++++++++++++++++++++++++++++++++++ python/tests/test_Rotation.py | 77 ---------------------------------- 2 files changed, 79 insertions(+), 77 deletions(-) diff --git a/python/tests/conftest.py b/python/tests/conftest.py index 411c07a8c..78ebea2a5 100644 --- a/python/tests/conftest.py +++ b/python/tests/conftest.py @@ -1,7 +1,10 @@ import os +import numpy as np import pytest +from damask import Rotation + def pytest_addoption(parser): parser.addoption("--update", action="store_true", @@ -16,3 +19,79 @@ def update(request): def reference_dir_base(): """Directory containing reference results.""" return os.path.join(os.path.dirname(__file__),'reference') + +@pytest.fixture +def set_of_rotations(): + """A set of n random rotations.""" + n = 1100 + scatter=1.e-2 + specials = np.array([ + [1.0, 0.0, 0.0, 0.0], + #---------------------- + [0.0, 1.0, 0.0, 0.0], + [0.0, 0.0, 1.0, 0.0], + [0.0, 0.0, 0.0, 1.0], + [0.0,-1.0, 0.0, 0.0], + [0.0, 0.0,-1.0, 0.0], + [0.0, 0.0, 0.0,-1.0], + #---------------------- + [1.0, 1.0, 0.0, 0.0], + [1.0, 0.0, 1.0, 0.0], + [1.0, 0.0, 0.0, 1.0], + [0.0, 1.0, 1.0, 0.0], + [0.0, 1.0, 0.0, 1.0], + [0.0, 0.0, 1.0, 1.0], + #---------------------- + [1.0,-1.0, 0.0, 0.0], + [1.0, 0.0,-1.0, 0.0], + [1.0, 0.0, 0.0,-1.0], + [0.0, 1.0,-1.0, 0.0], + [0.0, 1.0, 0.0,-1.0], + [0.0, 0.0, 1.0,-1.0], + #---------------------- + [0.0, 1.0,-1.0, 0.0], + [0.0, 1.0, 0.0,-1.0], + [0.0, 0.0, 1.0,-1.0], + #---------------------- + [0.0,-1.0,-1.0, 0.0], + [0.0,-1.0, 0.0,-1.0], + [0.0, 0.0,-1.0,-1.0], + #---------------------- + [1.0, 1.0, 1.0, 0.0], + [1.0, 1.0, 0.0, 1.0], + [1.0, 0.0, 1.0, 1.0], + [1.0,-1.0, 1.0, 0.0], + [1.0,-1.0, 0.0, 1.0], + [1.0, 0.0,-1.0, 1.0], + [1.0, 1.0,-1.0, 0.0], + [1.0, 1.0, 0.0,-1.0], + [1.0, 0.0, 1.0,-1.0], + [1.0,-1.0,-1.0, 0.0], + [1.0,-1.0, 0.0,-1.0], + [1.0, 0.0,-1.0,-1.0], + #---------------------- + [0.0, 1.0, 1.0, 1.0], + [0.0, 1.0,-1.0, 1.0], + [0.0, 1.0, 1.0,-1.0], + [0.0,-1.0, 1.0, 1.0], + [0.0,-1.0,-1.0, 1.0], + [0.0,-1.0, 1.0,-1.0], + [0.0,-1.0,-1.0,-1.0], + #---------------------- + [1.0, 1.0, 1.0, 1.0], + [1.0,-1.0, 1.0, 1.0], + [1.0, 1.0,-1.0, 1.0], + [1.0, 1.0, 1.0,-1.0], + [1.0,-1.0,-1.0, 1.0], + [1.0,-1.0, 1.0,-1.0], + [1.0, 1.0,-1.0,-1.0], + [1.0,-1.0,-1.0,-1.0], + ]) + specials /= np.linalg.norm(specials,axis=1).reshape(-1,1) + specials_scatter = specials + np.broadcast_to(np.random.rand(4)*scatter,specials.shape) + specials_scatter /= np.linalg.norm(specials_scatter,axis=1).reshape(-1,1) + specials_scatter[specials_scatter[:,0]<0]*=-1 + + return [Rotation.from_quaternion(s) for s in specials] + \ + [Rotation.from_quaternion(s) for s in specials_scatter] + \ + [Rotation.from_random() for _ in range(n-len(specials)-len(specials_scatter))] diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py index b087cc774..6acca8e79 100644 --- a/python/tests/test_Rotation.py +++ b/python/tests/test_Rotation.py @@ -9,83 +9,6 @@ from damask import _rotation n = 1000 atol=1.e-4 - -@pytest.fixture -def set_of_rotations(): - """A set of n random rotations.""" - n = 1100 - scatter=1.e-2 - specials = np.array([ - [1.0, 0.0, 0.0, 0.0], - #---------------------- - [0.0, 1.0, 0.0, 0.0], - [0.0, 0.0, 1.0, 0.0], - [0.0, 0.0, 0.0, 1.0], - [0.0,-1.0, 0.0, 0.0], - [0.0, 0.0,-1.0, 0.0], - [0.0, 0.0, 0.0,-1.0], - #---------------------- - [1.0, 1.0, 0.0, 0.0], - [1.0, 0.0, 1.0, 0.0], - [1.0, 0.0, 0.0, 1.0], - [0.0, 1.0, 1.0, 0.0], - [0.0, 1.0, 0.0, 1.0], - [0.0, 0.0, 1.0, 1.0], - #---------------------- - [1.0,-1.0, 0.0, 0.0], - [1.0, 0.0,-1.0, 0.0], - [1.0, 0.0, 0.0,-1.0], - [0.0, 1.0,-1.0, 0.0], - [0.0, 1.0, 0.0,-1.0], - [0.0, 0.0, 1.0,-1.0], - #---------------------- - [0.0, 1.0,-1.0, 0.0], - [0.0, 1.0, 0.0,-1.0], - [0.0, 0.0, 1.0,-1.0], - #---------------------- - [0.0,-1.0,-1.0, 0.0], - [0.0,-1.0, 0.0,-1.0], - [0.0, 0.0,-1.0,-1.0], - #---------------------- - [1.0, 1.0, 1.0, 0.0], - [1.0, 1.0, 0.0, 1.0], - [1.0, 0.0, 1.0, 1.0], - [1.0,-1.0, 1.0, 0.0], - [1.0,-1.0, 0.0, 1.0], - [1.0, 0.0,-1.0, 1.0], - [1.0, 1.0,-1.0, 0.0], - [1.0, 1.0, 0.0,-1.0], - [1.0, 0.0, 1.0,-1.0], - [1.0,-1.0,-1.0, 0.0], - [1.0,-1.0, 0.0,-1.0], - [1.0, 0.0,-1.0,-1.0], - #---------------------- - [0.0, 1.0, 1.0, 1.0], - [0.0, 1.0,-1.0, 1.0], - [0.0, 1.0, 1.0,-1.0], - [0.0,-1.0, 1.0, 1.0], - [0.0,-1.0,-1.0, 1.0], - [0.0,-1.0, 1.0,-1.0], - [0.0,-1.0,-1.0,-1.0], - #---------------------- - [1.0, 1.0, 1.0, 1.0], - [1.0,-1.0, 1.0, 1.0], - [1.0, 1.0,-1.0, 1.0], - [1.0, 1.0, 1.0,-1.0], - [1.0,-1.0,-1.0, 1.0], - [1.0,-1.0, 1.0,-1.0], - [1.0, 1.0,-1.0,-1.0], - [1.0,-1.0,-1.0,-1.0], - ]) - specials /= np.linalg.norm(specials,axis=1).reshape(-1,1) - specials_scatter = specials + np.broadcast_to(np.random.rand(4)*scatter,specials.shape) - specials_scatter /= np.linalg.norm(specials_scatter,axis=1).reshape(-1,1) - specials_scatter[specials_scatter[:,0]<0]*=-1 - - return [Rotation.from_quaternion(s) for s in specials] + \ - [Rotation.from_quaternion(s) for s in specials_scatter] + \ - [Rotation.from_random() for _ in range(n-len(specials)-len(specials_scatter))] - @pytest.fixture def reference_dir(reference_dir_base): """Directory containing reference results.""" From bdb461a5532c6ecf3f3cdca1f89926058ae18a3d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 30 Jun 2020 18:12:19 +0200 Subject: [PATCH 30/40] more flexible and independent --- python/tests/conftest.py | 23 +++++++++++++++++------ python/tests/test_Rotation.py | 8 ++++++-- 2 files changed, 23 insertions(+), 8 deletions(-) diff --git a/python/tests/conftest.py b/python/tests/conftest.py index 78ebea2a5..2be40a10a 100644 --- a/python/tests/conftest.py +++ b/python/tests/conftest.py @@ -3,8 +3,6 @@ import numpy as np import pytest -from damask import Rotation - def pytest_addoption(parser): parser.addoption("--update", action="store_true", @@ -21,8 +19,21 @@ def reference_dir_base(): return os.path.join(os.path.dirname(__file__),'reference') @pytest.fixture -def set_of_rotations(): +def set_of_quaternions(): """A set of n random rotations.""" + def random_quaternions(N): + r = np.random.rand(N,3) + + A = np.sqrt(r[:,2]) + B = np.sqrt(1.0-r[:,2]) + qu = np.column_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]) + qu[:,0]*=np.sign(qu[:,0]) + + return qu + n = 1100 scatter=1.e-2 specials = np.array([ @@ -92,6 +103,6 @@ def set_of_rotations(): specials_scatter /= np.linalg.norm(specials_scatter,axis=1).reshape(-1,1) specials_scatter[specials_scatter[:,0]<0]*=-1 - return [Rotation.from_quaternion(s) for s in specials] + \ - [Rotation.from_quaternion(s) for s in specials_scatter] + \ - [Rotation.from_random() for _ in range(n-len(specials)-len(specials_scatter))] + return [s for s in specials] + \ + [s for s in specials_scatter] + \ + [s for s in random_quaternions(n-2*len(specials))] diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py index 6acca8e79..a7f4dd17f 100644 --- a/python/tests/test_Rotation.py +++ b/python/tests/test_Rotation.py @@ -14,6 +14,10 @@ def reference_dir(reference_dir_base): """Directory containing reference results.""" return os.path.join(reference_dir_base,'Rotation') +@pytest.fixture +def set_of_rotations(set_of_quaternions): + return [Rotation.from_quaternion(s) for s in set_of_quaternions] + #################################################################################################### # Code below available according to the following conditions on https://github.com/MarDiehl/3Drotations @@ -607,9 +611,9 @@ class TestRotation: (Rotation._qu2ax,qu2ax), (Rotation._qu2ro,qu2ro), (Rotation._qu2ho,qu2ho)]) - def test_quaternion_vectorization(self,set_of_rotations,vectorized,single): + def test_quaternion_vectorization(self,set_of_quaternions,vectorized,single): """Check vectorized implementation for quaternion against single point calculation.""" - qu = np.array([rot.as_quaternion() for rot in set_of_rotations]) + qu = np.array(set_of_quaternions) vectorized(qu.reshape(qu.shape[0]//2,-1,4)) co = vectorized(qu) for q,c in zip(qu,co): From 9a83b11a9923d26fc143fa80e0e4a55d63b50606 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 30 Jun 2020 18:41:59 +0200 Subject: [PATCH 31/40] testing IPF color (vectorization) --- python/tests/conftest.py | 6 +++--- python/tests/test_Orientation.py | 30 +++++++++++++++++++++++------- python/tests/test_Rotation.py | 2 +- python/tests/test_ori_vec.py | 17 ----------------- 4 files changed, 27 insertions(+), 28 deletions(-) diff --git a/python/tests/conftest.py b/python/tests/conftest.py index 2be40a10a..af195ad6d 100644 --- a/python/tests/conftest.py +++ b/python/tests/conftest.py @@ -103,6 +103,6 @@ def set_of_quaternions(): specials_scatter /= np.linalg.norm(specials_scatter,axis=1).reshape(-1,1) specials_scatter[specials_scatter[:,0]<0]*=-1 - return [s for s in specials] + \ - [s for s in specials_scatter] + \ - [s for s in random_quaternions(n-2*len(specials))] + return np.array([s for s in specials] + \ + [s for s in specials_scatter] + \ + [s for s in random_quaternions(n-2*len(specials))]) diff --git a/python/tests/test_Orientation.py b/python/tests/test_Orientation.py index 277fa2a4b..750f6176d 100644 --- a/python/tests/test_Orientation.py +++ b/python/tests/test_Orientation.py @@ -11,6 +11,15 @@ from damask import Lattice n = 1000 +def IPF_color(orientation,direction): + """TSL color of inverse pole figure for given axis (non-vectorized).""" + for o in orientation.equivalent: + pole = o.rotation@direction + inSST,color = orientation.lattice.in_SST(pole,color=True) + if inSST: break + + return color + @pytest.fixture def reference_dir(reference_dir_base): """Directory containing reference results.""" @@ -26,16 +35,23 @@ class TestOrientation: def test_IPF_cubic(self,color,lattice): cube = damask.Orientation(damask.Rotation(),lattice) for direction in set(permutations(np.array(color['direction']))): - assert np.allclose(cube.IPFcolor(np.array(direction)),np.array(color['RGB'])) + assert np.allclose(cube.IPF_color(np.array(direction)),np.array(color['RGB'])) @pytest.mark.parametrize('lattice',Lattice.lattices) - def test_IPF(self,lattice): + def test_IPF_equivalent(self,set_of_quaternions,lattice): direction = np.random.random(3)*2.0-1 - for rot in [Rotation.from_random() for r in range(n//100)]: - R = damask.Orientation(rot,lattice) - color = R.IPFcolor(direction) - for equivalent in R.equivalent: - assert np.allclose(color,R.IPFcolor(direction)) + for ori in Orientation(Rotation(set_of_quaternions),lattice)[200]: + color = ori.IPF_color(direction) + for equivalent in ori.equivalent: + assert np.allclose(color,equivalent.IPF_color(direction)) + + + @pytest.mark.parametrize('lattice',Lattice.lattices) + def test_IPF_vectorize(self,set_of_quaternions,lattice): + for ori in Orientation(Rotation(set_of_quaternions),lattice)[200]: + direction = np.random.random(3)*2.0-1 + assert np.allclose(ori.IPF_color(direction),IPF_color(ori,direction)) + @pytest.mark.parametrize('model',['Bain','KS','GT','GT_prime','NW','Pitsch']) @pytest.mark.parametrize('lattice',['fcc','bcc']) diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py index a7f4dd17f..20d01af4a 100644 --- a/python/tests/test_Rotation.py +++ b/python/tests/test_Rotation.py @@ -6,7 +6,7 @@ import numpy as np from damask import Rotation from damask import _rotation -n = 1000 +n = 1100 atol=1.e-4 @pytest.fixture diff --git a/python/tests/test_ori_vec.py b/python/tests/test_ori_vec.py index ff5fe80bc..772096201 100644 --- a/python/tests/test_ori_vec.py +++ b/python/tests/test_ori_vec.py @@ -97,20 +97,3 @@ class TestOrientation_vec: assert all(ori_vec.reduced_vec.rotation.as_Rodrigues()[2] == ori2.reduced().rotation.as_Rodrigues() ) assert all(ori_vec.reduced_vec.rotation.as_cubochoric()[3] == ori3.reduced().rotation.as_cubochoric() ) assert all(ori_vec.reduced_vec.rotation.as_axis_angle()[4] == ori4.reduced().rotation.as_axis_angle() ) - - - @pytest.mark.parametrize('lattice',['bcc','fcc','bct']) - def test_IPFcolor_vec(self,lattice): - ori0=Orientation(rot0,lattice) - ori1=Orientation(rot1,lattice) - ori2=Orientation(rot2,lattice) - ori3=Orientation(rot3,lattice) - - quat=np.array([rot0.as_quaternion(),rot1.as_quaternion(),\ - rot2.as_quaternion(),rot3.as_quaternion()]) - ori_vec=Orientation(quat,lattice) - - assert np.allclose( ori_vec.IPF_color(np.array([0,0,1]))[0],ori0.IPFcolor(np.array([0,0,1]))) - assert np.allclose( ori_vec.IPF_color(np.array([0,2,1]))[1],ori1.IPFcolor(np.array([0,2,1]))) - assert np.allclose( ori_vec.IPF_color(np.array([0,3,1]))[2],ori2.IPFcolor(np.array([0,3,1]))) - assert np.allclose( ori_vec.IPF_color(np.array([4,0,1]))[3],ori3.IPFcolor(np.array([4,0,1]))) From 49d448dcede646f8a231d4befc9da386fea01cd6 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 30 Jun 2020 21:43:57 +0200 Subject: [PATCH 32/40] vectorized and cleaned --- python/damask/_lattice.py | 260 ++++++------------------------- python/damask/_orientation.py | 67 ++------ python/damask/_result.py | 6 +- python/damask/_rotation.py | 2 +- python/tests/test_Lattice.py | 122 ++++++++++++++- python/tests/test_Orientation.py | 9 +- python/tests/test_Result.py | 6 +- python/tests/test_ori_vec.py | 38 ----- 8 files changed, 194 insertions(+), 316 deletions(-) diff --git a/python/damask/_lattice.py b/python/damask/_lattice.py index ee718adcd..74ef61a50 100644 --- a/python/damask/_lattice.py +++ b/python/damask/_lattice.py @@ -29,7 +29,7 @@ class Symmetry: raise KeyError(f'Crystal system "{system}" is unknown') self.system = system.lower() if isinstance(system,str) else system - self.lattice = self.system # for compatibility + self.lattice = self.system # ToDo: for compatibility def __copy__(self): @@ -82,85 +82,10 @@ class Symmetry: otherOrder = self.crystal_systems.index(other.system) return (myOrder > otherOrder) - (myOrder < otherOrder) - def symmetryOperations(self,members=[]): - """List (or single element) of symmetry operations as rotations.""" - if self.system == 'cubic': - symQuats = [ - [ 1.0, 0.0, 0.0, 0.0 ], - [ 0.0, 1.0, 0.0, 0.0 ], - [ 0.0, 0.0, 1.0, 0.0 ], - [ 0.0, 0.0, 0.0, 1.0 ], - [ 0.0, 0.0, 0.5*np.sqrt(2), 0.5*np.sqrt(2) ], - [ 0.0, 0.0, 0.5*np.sqrt(2),-0.5*np.sqrt(2) ], - [ 0.0, 0.5*np.sqrt(2), 0.0, 0.5*np.sqrt(2) ], - [ 0.0, 0.5*np.sqrt(2), 0.0, -0.5*np.sqrt(2) ], - [ 0.0, 0.5*np.sqrt(2),-0.5*np.sqrt(2), 0.0 ], - [ 0.0, -0.5*np.sqrt(2),-0.5*np.sqrt(2), 0.0 ], - [ 0.5, 0.5, 0.5, 0.5 ], - [-0.5, 0.5, 0.5, 0.5 ], - [-0.5, 0.5, 0.5, -0.5 ], - [-0.5, 0.5, -0.5, 0.5 ], - [-0.5, -0.5, 0.5, 0.5 ], - [-0.5, -0.5, 0.5, -0.5 ], - [-0.5, -0.5, -0.5, 0.5 ], - [-0.5, 0.5, -0.5, -0.5 ], - [-0.5*np.sqrt(2), 0.0, 0.0, 0.5*np.sqrt(2) ], - [ 0.5*np.sqrt(2), 0.0, 0.0, 0.5*np.sqrt(2) ], - [-0.5*np.sqrt(2), 0.0, 0.5*np.sqrt(2), 0.0 ], - [-0.5*np.sqrt(2), 0.0, -0.5*np.sqrt(2), 0.0 ], - [-0.5*np.sqrt(2), 0.5*np.sqrt(2), 0.0, 0.0 ], - [-0.5*np.sqrt(2),-0.5*np.sqrt(2), 0.0, 0.0 ], - ] - elif self.system == 'hexagonal': - symQuats = [ - [ 1.0, 0.0, 0.0, 0.0 ], - [-0.5*np.sqrt(3), 0.0, 0.0, -0.5 ], - [ 0.5, 0.0, 0.0, 0.5*np.sqrt(3) ], - [ 0.0, 0.0, 0.0, 1.0 ], - [-0.5, 0.0, 0.0, 0.5*np.sqrt(3) ], - [-0.5*np.sqrt(3), 0.0, 0.0, 0.5 ], - [ 0.0, 1.0, 0.0, 0.0 ], - [ 0.0, -0.5*np.sqrt(3), 0.5, 0.0 ], - [ 0.0, 0.5, -0.5*np.sqrt(3), 0.0 ], - [ 0.0, 0.0, 1.0, 0.0 ], - [ 0.0, -0.5, -0.5*np.sqrt(3), 0.0 ], - [ 0.0, 0.5*np.sqrt(3), 0.5, 0.0 ], - ] - elif self.system == 'tetragonal': - symQuats = [ - [ 1.0, 0.0, 0.0, 0.0 ], - [ 0.0, 1.0, 0.0, 0.0 ], - [ 0.0, 0.0, 1.0, 0.0 ], - [ 0.0, 0.0, 0.0, 1.0 ], - [ 0.0, 0.5*np.sqrt(2), 0.5*np.sqrt(2), 0.0 ], - [ 0.0, -0.5*np.sqrt(2), 0.5*np.sqrt(2), 0.0 ], - [ 0.5*np.sqrt(2), 0.0, 0.0, 0.5*np.sqrt(2) ], - [-0.5*np.sqrt(2), 0.0, 0.0, 0.5*np.sqrt(2) ], - ] - elif self.system == 'orthorhombic': - symQuats = [ - [ 1.0,0.0,0.0,0.0 ], - [ 0.0,1.0,0.0,0.0 ], - [ 0.0,0.0,1.0,0.0 ], - [ 0.0,0.0,0.0,1.0 ], - ] - else: - symQuats = [ - [ 1.0,0.0,0.0,0.0 ], - ] - - symOps = list(map(Rotation, - np.array(symQuats)[np.atleast_1d(members) if members != [] else range(len(symQuats))])) - try: - iter(members) # asking for (even empty) list of members? - except TypeError: - return symOps[0] # no, return rotation object - else: - return symOps # yes, return list of rotations @property def symmetry_operations(self): - """Symmetry operations as Rotations.""" + """Symmetry operations as Quaternions.""" if self.system == 'cubic': symQuats = [ [ 1.0, 0.0, 0.0, 0.0 ], @@ -228,42 +153,40 @@ class Symmetry: return np.array(symQuats) - def inFZ(self,rodrigues): + def in_FZ(self,rho): """ - Check whether given Rodrigues-Frank vector falls into fundamental zone of own symmetry. + Check whether given Rodrigues-Frank vector falls into fundamental zone. Fundamental zone in Rodrigues space is point symmetric around origin. """ - if (len(rodrigues) != 3): - raise ValueError('Input is not a Rodrigues-Frank vector.\n') + if(rho.shape[-1] != 3): + raise ValueError('Input is not a Rodrigues-Frank vector.') - if np.any(rodrigues == np.inf): return False # ToDo: MD: not sure if needed + rho_abs = np.abs(rho) - Rabs = abs(rodrigues) - - if self.system == 'cubic': - return np.sqrt(2.0)-1.0 >= Rabs[0] \ - and np.sqrt(2.0)-1.0 >= Rabs[1] \ - and np.sqrt(2.0)-1.0 >= Rabs[2] \ - and 1.0 >= Rabs[0] + Rabs[1] + Rabs[2] - elif self.system == 'hexagonal': - return 1.0 >= Rabs[0] and 1.0 >= Rabs[1] and 1.0 >= Rabs[2] \ - and 2.0 >= np.sqrt(3)*Rabs[0] + Rabs[1] \ - and 2.0 >= np.sqrt(3)*Rabs[1] + Rabs[0] \ - and 2.0 >= np.sqrt(3) + Rabs[2] - elif self.system == 'tetragonal': - return 1.0 >= Rabs[0] and 1.0 >= Rabs[1] \ - and np.sqrt(2.0) >= Rabs[0] + Rabs[1] \ - and np.sqrt(2.0) >= Rabs[2] + 1.0 - elif self.system == 'orthorhombic': - return 1.0 >= Rabs[0] and 1.0 >= Rabs[1] and 1.0 >= Rabs[2] - else: - return True + with np.errstate(invalid='ignore'): + # using '*'/prod for 'and' + if self.system == 'cubic': + return np.where(np.prod(np.sqrt(2)-1. >= rho_abs,axis=-1) * \ + (1. >= np.sum(rho_abs,axis=-1)),True,False) + elif self.system == 'hexagonal': + return np.where(np.prod(1. >= rho_abs,axis=-1) * \ + (2. >= np.sqrt(3)*rho_abs[...,0] + rho_abs[...,1]) * \ + (2. >= np.sqrt(3)*rho_abs[...,1] + rho_abs[...,0]) * \ + (2. >= np.sqrt(3) + rho_abs[...,2]),True,False) + elif self.system == 'tetragonal': + return np.where(np.prod(1. >= rho_abs[...,:2],axis=-1) * \ + (np.sqrt(2) >= rho_abs[...,0] + rho_abs[...,1]) * \ + (np.sqrt(2) >= rho_abs[...,2] + 1.),True,False) + elif self.system == 'orthorhombic': + return np.where(np.prod(1. >= rho_abs,axis=-1),True,False) + else: + return np.where(np.all(np.isfinite(rho_abs),axis=-1),True,False) - def inDisorientationSST(self,rodrigues): + def in_disorientation_SST(self,rho): """ - Check whether given Rodrigues-Frank vector (of misorientation) falls into standard stereographic triangle of own symmetry. + Check whether given Rodrigues-Frank vector (of misorientation) falls into standard stereographic triangle. References ---------- @@ -271,115 +194,32 @@ class Symmetry: https://doi.org/10.1107/S0108767391006864 """ - if (len(rodrigues) != 3): - raise ValueError('Input is not a Rodrigues-Frank vector.\n') - R = rodrigues + if(rho.shape[-1] != 3): + raise ValueError('Input is not a Rodrigues-Frank vector.') - epsilon = 0.0 - if self.system == 'cubic': - return R[0] >= R[1]+epsilon and R[1] >= R[2]+epsilon and R[2] >= epsilon - elif self.system == 'hexagonal': - return R[0] >= np.sqrt(3)*(R[1]-epsilon) and R[1] >= epsilon and R[2] >= epsilon - elif self.system == 'tetragonal': - return R[0] >= R[1]-epsilon and R[1] >= epsilon and R[2] >= epsilon - elif self.system == 'orthorhombic': - return R[0] >= epsilon and R[1] >= epsilon and R[2] >= epsilon - else: - return True - - - def inSST(self, - vector, - proper = False, - color = False): - """ - Check whether given vector falls into standard stereographic triangle of own symmetry. - - proper considers only vectors with z >= 0, hence uses two neighboring SSTs. - Return inverse pole figure color if requested. - Bases are computed from - - >>> basis = {'cubic' : np.linalg.inv(np.array([[0.,0.,1.], # direction of red - ... [1.,0.,1.]/np.sqrt(2.), # direction of green - ... [1.,1.,1.]/np.sqrt(3.)]).T), # direction of blue - ... 'hexagonal' : np.linalg.inv(np.array([[0.,0.,1.], # direction of red - ... [1.,0.,0.], # direction of green - ... [np.sqrt(3.),1.,0.]/np.sqrt(4.)]).T), # direction of blue - ... 'tetragonal' : np.linalg.inv(np.array([[0.,0.,1.], # direction of red - ... [1.,0.,0.], # direction of green - ... [1.,1.,0.]/np.sqrt(2.)]).T), # direction of blue - ... 'orthorhombic': np.linalg.inv(np.array([[0.,0.,1.], # direction of red - ... [1.,0.,0.], # direction of green - ... [0.,1.,0.]]).T), # direction of blue - ... } - - """ - if self.system == 'cubic': - basis = {'improper':np.array([ [-1. , 0. , 1. ], - [ np.sqrt(2.) , -np.sqrt(2.) , 0. ], - [ 0. , np.sqrt(3.) , 0. ] ]), - 'proper':np.array([ [ 0. , -1. , 1. ], - [-np.sqrt(2.) , np.sqrt(2.) , 0. ], - [ np.sqrt(3.) , 0. , 0. ] ]), - } - elif self.system == 'hexagonal': - basis = {'improper':np.array([ [ 0. , 0. , 1. ], - [ 1. , -np.sqrt(3.) , 0. ], - [ 0. , 2. , 0. ] ]), - 'proper':np.array([ [ 0. , 0. , 1. ], - [-1. , np.sqrt(3.) , 0. ], - [ np.sqrt(3.) , -1. , 0. ] ]), - } - elif self.system == 'tetragonal': - basis = {'improper':np.array([ [ 0. , 0. , 1. ], - [ 1. , -1. , 0. ], - [ 0. , np.sqrt(2.) , 0. ] ]), - 'proper':np.array([ [ 0. , 0. , 1. ], - [-1. , 1. , 0. ], - [ np.sqrt(2.) , 0. , 0. ] ]), - } - elif self.system == 'orthorhombic': - basis = {'improper':np.array([ [ 0., 0., 1.], - [ 1., 0., 0.], - [ 0., 1., 0.] ]), - 'proper':np.array([ [ 0., 0., 1.], - [-1., 0., 0.], - [ 0., 1., 0.] ]), - } - else: # direct exit for unspecified symmetry - if color: - return (True,np.zeros(3,'d')) + with np.errstate(invalid='ignore'): + # using '*' for 'and' + if self.system == 'cubic': + return np.where((rho[...,0] >= rho[...,1]) * \ + (rho[...,1] >= rho[...,2]) * \ + (rho[...,2] >= 0),True,False) + elif self.system == 'hexagonal': + return np.where((rho[...,0] >= rho[...,1]*np.sqrt(3)) * \ + (rho[...,1] >= 0) * \ + (rho[...,2] >= 0),True,False) + elif self.system == 'tetragonal': + return np.where((rho[...,0] >= rho[...,1]) * \ + (rho[...,1] >= 0) * \ + (rho[...,2] >= 0),True,False) + elif self.system == 'orthorhombic': + return np.where((rho[...,0] >= 0) * \ + (rho[...,1] >= 0) * \ + (rho[...,2] >= 0),True,False) else: - return True - - v = np.array(vector,dtype=float) - if proper: # check both improper ... - theComponents = np.around(np.dot(basis['improper'],v),12) - inSST = np.all(theComponents >= 0.0) - if not inSST: # ... and proper SST - theComponents = np.around(np.dot(basis['proper'],v),12) - inSST = np.all(theComponents >= 0.0) - else: - v[2] = abs(v[2]) # z component projects identical - theComponents = np.around(np.dot(basis['improper'],v),12) # for positive and negative values - inSST = np.all(theComponents >= 0.0) - - if color: # have to return color array - if inSST: - rgb = np.power(theComponents/np.linalg.norm(theComponents),0.5) # smoothen color ramps - rgb = np.minimum(np.ones(3,dtype=float),rgb) # limit to maximum intensity - rgb /= max(rgb) # normalize to (HS)V = 1 - else: - rgb = np.zeros(3,dtype=float) - return (inSST,rgb) - else: - return inSST + return np.ones_like(rho[...,0],dtype=bool) - def in_SST(self, - vector, - proper = False, - color = False): + def in_SST(self,vector,proper=False,color=False): """ Check whether given vector falls into standard stereographic triangle of own symmetry. @@ -503,7 +343,7 @@ class Lattice: # ToDo: Make a subclass of Symmetry! # transition to subclass self.system = self.symmetry.system self.in_SST = self.symmetry.in_SST - self.inFZ = self.symmetry.inFZ + self.in_FZ = self.symmetry.in_FZ def __repr__(self): diff --git a/python/damask/_orientation.py b/python/damask/_orientation.py index 7d0f11c49..b6b24e951 100644 --- a/python/damask/_orientation.py +++ b/python/damask/_orientation.py @@ -40,8 +40,6 @@ class Orientation: # ToDo: make subclass of lattice and Rotation self.rotation = Rotation.from_quaternion(rotation) # assume quaternion def __getitem__(self,item): - if isinstance(item,tuple) and len(item) >= len(self): - raise IndexError('Too many indices') return self.__class__(self.rotation[item],self.lattice) @@ -61,8 +59,8 @@ class Orientation: # ToDo: make subclass of lattice and Rotation if self.lattice.symmetry != other.lattice.symmetry: raise NotImplementedError('disorientation between different symmetry classes not supported yet.') - mySymEqs = self.equivalentOrientations() if SST else self.equivalentOrientations([0]) # take all or only first sym operation - otherSymEqs = other.equivalentOrientations() + mySymEqs = self.equivalent if SST else self.equivalent[0] # take all or only first sym operation + otherSymEqs = other.equivalent for i,sA in enumerate(mySymEqs): aInv = sA.rotation.inversed() @@ -71,7 +69,7 @@ class Orientation: # ToDo: make subclass of lattice and Rotation r = b*aInv for k in range(2): r.inverse() - breaker = self.lattice.inFZ(r.as_Rodrigues(vector=True)) \ + breaker = self.in_FZ \ and (not SST or other.lattice.symmetry.inDisorientationSST(r.as_Rodrigues(vector=True))) if breaker: break if breaker: break @@ -81,17 +79,10 @@ class Orientation: # ToDo: make subclass of lattice and Rotation # ... own sym, other sym, # self-->other: True, self<--other: False - def inFZ_vec(self): + + def in_FZ(self): """Check if orientations fall into Fundamental Zone.""" - if not self.rotation.shape: - return self.lattice.inFZ(self.rotation.as_Rodrigues(vector=True)) - else: - return [self.lattice.inFZ(\ - self.rotation.as_Rodrigues(vector=True)[l]) for l in range(self.rotation.shape[0])] - - - def inFZ(self): - return self.lattice.inFZ(self.rotation.as_Rodrigues(vector=True)) + return self.lattice.in_FZ(self.rotation.as_Rodrigues(vector=True)) @property def equivalent(self): @@ -112,16 +103,6 @@ class Orientation: # ToDo: make subclass of lattice and Rotation return self.__class__(s@r,self.lattice) - def equivalentOrientations(self,members=[]): - """List of orientations which are symmetrically equivalent.""" - try: - iter(members) # asking for (even empty) list of members? - except TypeError: - return self.__class__(self.lattice.symmetry.symmetryOperations(members)*self.rotation,self.lattice) # no, return rotation object - else: - return [self.__class__(q*self.rotation,self.lattice) \ - for q in self.lattice.symmetry.symmetryOperations(members)] # yes, return list of rotations - def relatedOrientations_vec(self,model): """List of orientations related by the given orientation relationship.""" h = self.lattice.relationOperations(model) @@ -149,7 +130,7 @@ class Orientation: # ToDo: make subclass of lattice and Rotation r= 1 if not self.rotation.shape else equi.shape[1] #number of rotations num_equi=equi.shape[0] #number of equivalente orientations quat= np.reshape( equi.as_quaternion(), (r*num_equi,4) ,order='F') #equivalents are listed in intiuitive order - boolean=Orientation(quat, self.lattice).inFZ_vec() #check which ones are in FZ + boolean=Orientation(quat, self.lattice).in_FZ() #check which ones are in FZ if sum(boolean) == r: return self.__class__(quat[boolean],self.lattice) @@ -162,13 +143,10 @@ class Orientation: # ToDo: make subclass of lattice and Rotation return self.__class__(quat[index],self.lattice) - - - def reduced(self): """Transform orientation to fall into fundamental zone according to symmetry.""" - for me in self.equivalentOrientations(): - if self.lattice.inFZ(me.rotation.as_Rodrigues(vector=True)): break + for me in self.equivalent: + if self.lattice.in_FZ(me.rotation.as_Rodrigues(vector=True)): break return self.__class__(me.rotation,self.lattice) @@ -179,35 +157,23 @@ class Orientation: # ToDo: make subclass of lattice and Rotation SST = True): """Axis rotated according to orientation (using crystal symmetry to ensure location falls into SST).""" if SST: # pole requested to be within SST - for i,o in enumerate(self.equivalentOrientations()): # test all symmetric equivalent quaternions + for i,o in enumerate(self.equivalent): # test all symmetric equivalent quaternions pole = o.rotation@axis # align crystal direction to axis - if self.lattice.symmetry.inSST(pole,proper): break # found SST version + if self.lattice.in_SST(pole,proper): break # found SST version else: pole = self.rotation@axis # align crystal direction to axis return (pole,i if SST else 0) - def IPFcolor(self,axis): + def IPF_color(self,axis): #ToDo axis or direction? """TSL color of inverse pole figure for given axis.""" - color = np.zeros(3,'d') - - for o in self.equivalent: - pole = o.rotation@axis # align crystal direction to axis - inSST,color = self.lattice.symmetry.inSST(pole,color=True) - if inSST: break - - return color - - - def IPF_color(self,axis): - """TSL color of inverse pole figure for given axis. Not for hex or triclinic lattices.""" eq = self.equivalent pole = eq.rotation @ np.broadcast_to(axis/np.linalg.norm(axis),eq.rotation.shape+(3,)) in_SST, color = self.lattice.in_SST(pole,color=True) # ignore duplicates (occur for highly symmetric orientations) - found = np.zeros_like(in_SST[1],dtype=bool) + found = np.zeros_like(in_SST[0],dtype=bool) c = np.empty(color.shape[1:]) for s in range(in_SST.shape[0]): c = np.where(np.expand_dims(np.logical_and(in_SST[s],~found),-1),color[s],c) @@ -220,17 +186,18 @@ class Orientation: # ToDo: make subclass of lattice and Rotation def fromAverage(orientations, weights = []): """Create orientation from average of list of orientations.""" - # further read: Orientation distribution analysis in deformed grains, https://doi.org/10.1107/S0021889801003077 + # further read: Orientation distribution analysis in deformed grains + # https://doi.org/10.1107/S0021889801003077 if not all(isinstance(item, Orientation) for item in orientations): raise TypeError("Only instances of Orientation can be averaged.") closest = [] ref = orientations[0] for o in orientations: - closest.append(o.equivalentOrientations( + closest.append(o.equivalent[ ref.disorientation(o, SST = False, # select (o[ther]'s) sym orientation - symmetries = True)[2]).rotation) # with lowest misorientation + symmetries = True)[2]].rotation) # with lowest misorientation return Orientation(Rotation.fromAverage(closest,weights),ref.lattice) diff --git a/python/damask/_result.py b/python/damask/_result.py index e6dac9370..1396b560d 100644 --- a/python/damask/_result.py +++ b/python/damask/_result.py @@ -733,7 +733,7 @@ class Result: @staticmethod - def _add_IPFcolor(q,l): + def _add_IPF_color(q,l): m = util.scale_to_coprime(np.array(l)) o = Orientation(Rotation(rfn.structured_to_unstructured(q['data'])), @@ -749,7 +749,7 @@ class Result: 'Creator': inspect.stack()[0][3][1:] } } - def add_IPFcolor(self,q,l): + def add_IPF_color(self,q,l): """ Add RGB color tuple of inverse pole figure (IPF) color. @@ -761,7 +761,7 @@ class Result: Lab frame direction for inverse pole figure. """ - self._add_generic_pointwise(self._add_IPFcolor,{'q':q},{'l':l}) + self._add_generic_pointwise(self._add_IPF_color,{'q':q},{'l':l}) @staticmethod diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index 686b144ae..accd453cc 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -470,7 +470,7 @@ class Rotation: elif hasattr(shape, '__iter__'): r = np.random.random(tuple(shape)+(3,)) else: - r = np.random.random((shape,3)) + r = np.random.rand(shape,3) A = np.sqrt(r[...,2]) B = np.sqrt(1.0-r[...,2]) diff --git a/python/tests/test_Lattice.py b/python/tests/test_Lattice.py index 3d9c35706..8e4616660 100644 --- a/python/tests/test_Lattice.py +++ b/python/tests/test_Lattice.py @@ -3,10 +3,118 @@ import random import pytest import numpy as np +from damask import Orientation +from damask import Rotation from damask import Symmetry +def in_FZ(system,rho): + """Non-vectorized version of 'in_FZ'.""" + rho_abs = abs(rho) + + if system == 'cubic': + return np.sqrt(2.0)-1.0 >= rho_abs[0] \ + and np.sqrt(2.0)-1.0 >= rho_abs[1] \ + and np.sqrt(2.0)-1.0 >= rho_abs[2] \ + and 1.0 >= rho_abs[0] + rho_abs[1] + rho_abs[2] + elif system == 'hexagonal': + return 1.0 >= rho_abs[0] and 1.0 >= rho_abs[1] and 1.0 >= rho_abs[2] \ + and 2.0 >= np.sqrt(3)*rho_abs[0] + rho_abs[1] \ + and 2.0 >= np.sqrt(3)*rho_abs[1] + rho_abs[0] \ + and 2.0 >= np.sqrt(3) + rho_abs[2] + elif system == 'tetragonal': + return 1.0 >= rho_abs[0] and 1.0 >= rho_abs[1] \ + and np.sqrt(2.0) >= rho_abs[0] + rho_abs[1] \ + and np.sqrt(2.0) >= rho_abs[2] + 1.0 + elif system == 'orthorhombic': + return 1.0 >= rho_abs[0] and 1.0 >= rho_abs[1] and 1.0 >= rho_abs[2] + else: + return np.all(np.isfinite(rho_abs)) + + +def in_disorientation_SST(system,rho): + """Non-vectorized version of 'in_Disorientation_SST'.""" + epsilon = 0.0 + if system == 'cubic': + return rho[0] >= rho[1]+epsilon and rho[1] >= rho[2]+epsilon and rho[2] >= epsilon + elif system == 'hexagonal': + return rho[0] >= np.sqrt(3)*(rho[1]-epsilon) and rho[1] >= epsilon and rho[2] >= epsilon + elif system == 'tetragonal': + return rho[0] >= rho[1]-epsilon and rho[1] >= epsilon and rho[2] >= epsilon + elif system == 'orthorhombic': + return rho[0] >= epsilon and rho[1] >= epsilon and rho[2] >= epsilon + else: + return True + + +def in_SST(system,vector,proper = False): + """Non-vectorized version of 'in_SST'.""" + if system == 'cubic': + basis = {'improper':np.array([ [-1. , 0. , 1. ], + [ np.sqrt(2.) , -np.sqrt(2.) , 0. ], + [ 0. , np.sqrt(3.) , 0. ] ]), + 'proper':np.array([ [ 0. , -1. , 1. ], + [-np.sqrt(2.) , np.sqrt(2.) , 0. ], + [ np.sqrt(3.) , 0. , 0. ] ]), + } + elif system == 'hexagonal': + basis = {'improper':np.array([ [ 0. , 0. , 1. ], + [ 1. , -np.sqrt(3.) , 0. ], + [ 0. , 2. , 0. ] ]), + 'proper':np.array([ [ 0. , 0. , 1. ], + [-1. , np.sqrt(3.) , 0. ], + [ np.sqrt(3.) , -1. , 0. ] ]), + } + elif system == 'tetragonal': + basis = {'improper':np.array([ [ 0. , 0. , 1. ], + [ 1. , -1. , 0. ], + [ 0. , np.sqrt(2.) , 0. ] ]), + 'proper':np.array([ [ 0. , 0. , 1. ], + [-1. , 1. , 0. ], + [ np.sqrt(2.) , 0. , 0. ] ]), + } + elif system == 'orthorhombic': + basis = {'improper':np.array([ [ 0., 0., 1.], + [ 1., 0., 0.], + [ 0., 1., 0.] ]), + 'proper':np.array([ [ 0., 0., 1.], + [-1., 0., 0.], + [ 0., 1., 0.] ]), + } + else: + return True + + v = np.array(vector,dtype=float) + if proper: + theComponents = np.around(np.dot(basis['improper'],v),12) + inSST = np.all(theComponents >= 0.0) + if not inSST: + theComponents = np.around(np.dot(basis['proper'],v),12) + inSST = np.all(theComponents >= 0.0) + else: + v[2] = abs(v[2]) + theComponents = np.around(np.dot(basis['improper'],v),12) + inSST = np.all(theComponents >= 0.0) + + return inSST + + +@pytest.fixture +def set_of_rodrigues(set_of_quaternions): + return Rotation(set_of_quaternions).as_Rodrigues(vector=True)[:200] + class TestSymmetry: + @pytest.mark.parametrize('system',Symmetry.crystal_systems) + def test_in_FZ_vectorize(self,set_of_rodrigues,system): + for i,in_FZ_ in enumerate(Symmetry(system).in_FZ(set_of_rodrigues)): + assert in_FZ_ == in_FZ(system,set_of_rodrigues[i]) + + @pytest.mark.parametrize('system',Symmetry.crystal_systems) + def test_in_disorientation_SST_vectorize(self,set_of_rodrigues,system): + for i,in_disorientation_SST_ in enumerate(Symmetry(system).in_disorientation_SST(set_of_rodrigues)): + assert in_disorientation_SST_ == in_disorientation_SST(system,set_of_rodrigues[i]) + + @pytest.mark.parametrize('invalid_symmetry',['fcc','bcc','hello']) def test_invalid_symmetry(self,invalid_symmetry): with pytest.raises(KeyError): @@ -22,19 +130,19 @@ class TestSymmetry: assert Symmetry(symmetries[0]) != Symmetry(symmetries[1]) @pytest.mark.parametrize('system',Symmetry.crystal_systems) - def test_inFZ(self,system): - assert Symmetry(system).inFZ(np.zeros(3)) + def test_in_FZ(self,system): + assert Symmetry(system).in_FZ(np.zeros(3)) @pytest.mark.parametrize('system',Symmetry.crystal_systems) - def test_inDisorientationSST(self,system): - assert Symmetry(system).inDisorientationSST(np.zeros(3)) + def test_in_disorientation_SST(self,system): + assert Symmetry(system).in_disorientation_SST(np.zeros(3)) @pytest.mark.parametrize('system',Symmetry.crystal_systems) @pytest.mark.parametrize('proper',[True,False]) - def test_inSST(self,system,proper): - assert Symmetry(system).inSST(np.zeros(3),proper) + def test_in_SST(self,system,proper): + assert Symmetry(system).in_SST(np.zeros(3),proper) - @pytest.mark.parametrize('function',['inFZ','inDisorientationSST']) + @pytest.mark.parametrize('function',['in_FZ','in_disorientation_SST']) def test_invalid_argument(self,function): s = Symmetry() # noqa with pytest.raises(ValueError): diff --git a/python/tests/test_Orientation.py b/python/tests/test_Orientation.py index 750f6176d..f6f25e0a7 100644 --- a/python/tests/test_Orientation.py +++ b/python/tests/test_Orientation.py @@ -40,7 +40,7 @@ class TestOrientation: @pytest.mark.parametrize('lattice',Lattice.lattices) def test_IPF_equivalent(self,set_of_quaternions,lattice): direction = np.random.random(3)*2.0-1 - for ori in Orientation(Rotation(set_of_quaternions),lattice)[200]: + for ori in Orientation(Rotation(set_of_quaternions),lattice)[:200]: color = ori.IPF_color(direction) for equivalent in ori.equivalent: assert np.allclose(color,equivalent.IPF_color(direction)) @@ -48,9 +48,10 @@ class TestOrientation: @pytest.mark.parametrize('lattice',Lattice.lattices) def test_IPF_vectorize(self,set_of_quaternions,lattice): - for ori in Orientation(Rotation(set_of_quaternions),lattice)[200]: - direction = np.random.random(3)*2.0-1 - assert np.allclose(ori.IPF_color(direction),IPF_color(ori,direction)) + direction = np.random.random(3)*2.0-1 + oris = Orientation(Rotation(set_of_quaternions),lattice)[:200] + for i,color in enumerate(oris.IPF_color(direction)): + assert np.allclose(color,IPF_color(oris[i],direction)) @pytest.mark.parametrize('model',['Bain','KS','GT','GT_prime','NW','Pitsch']) diff --git a/python/tests/test_Result.py b/python/tests/test_Result.py index d7946e5e0..aec91db9f 100644 --- a/python/tests/test_Result.py +++ b/python/tests/test_Result.py @@ -153,8 +153,8 @@ class TestResult: assert np.allclose(in_memory,in_file) @pytest.mark.parametrize('d',[[1,0,0],[0,1,0],[0,0,1]]) - def test_add_IPFcolor(self,default,d): - default.add_IPFcolor('orientation',d) + def test_add_IPF_color(self,default,d): + default.add_IPF_color('orientation',d) loc = {'orientation': default.get_dataset_location('orientation'), 'color': default.get_dataset_location('IPFcolor_[{} {} {}]'.format(*d))} qu = default.read_dataset(loc['orientation']).view(np.double).reshape(-1,4) @@ -162,7 +162,7 @@ class TestResult: in_memory = np.empty((qu.shape[0],3),np.uint8) for i,q in enumerate(qu): o = damask.Orientation(q,crystal_structure).reduced() - in_memory[i] = np.uint8(o.IPFcolor(np.array(d))*255) + in_memory[i] = np.uint8(o.IPF_color(np.array(d))*255) in_file = default.read_dataset(loc['color']) assert np.allclose(in_memory,in_file) diff --git a/python/tests/test_ori_vec.py b/python/tests/test_ori_vec.py index 772096201..7cb465046 100644 --- a/python/tests/test_ori_vec.py +++ b/python/tests/test_ori_vec.py @@ -16,44 +16,6 @@ rot3= Rotation.from_random() #average class TestOrientation_vec: - #@pytest.mark.xfail - @pytest.mark.parametrize('lattice',Lattice.lattices) - def test_equivalent_vec(self,lattice): - ori0=Orientation(rot0,lattice) - ori1=Orientation(rot1,lattice) - ori2=Orientation(rot2,lattice) - ori3=Orientation(rot3,lattice) - - quat=np.array([rot0.as_quaternion(),rot1.as_quaternion(),rot2.as_quaternion(),rot3.as_quaternion()]) - ori_vec=Orientation(quat,lattice) - - for s in range(len(ori_vec.lattice.symmetry.symmetryOperations())): - assert all(ori_vec.equivalent.rotation.as_Eulers()[s,0] == \ - ori0.equivalent[s].rotation.as_Eulers()) - assert all(ori_vec.equivalent.rotation.as_quaternion()[s,1] == \ - ori1.equivalent[s].rotation.as_quaternion()) - assert all(ori_vec.equivalent.rotation.as_Rodrigues()[s,2] == \ - ori2.equivalent[s].rotation.as_Rodrigues()) - assert all(ori_vec.equivalent.rotation.as_cubochoric()[s,3] == \ - ori3.equivalent[s].rotation.as_cubochoric()) - - @pytest.mark.parametrize('lattice',Lattice.lattices) - def test_inFZ_vec(self,lattice): - ori0=Orientation(rot0,lattice) - ori1=Orientation(rot1,lattice) - ori2=Orientation(rot2,lattice) - ori3=Orientation(rot3,lattice) - ori4=ori0.reduced() ; rot4=ori4.rotation #ensure 1 of them is in FZ - - quat=np.array([rot0.as_quaternion(),rot1.as_quaternion(),\ - rot2.as_quaternion(),rot3.as_quaternion(), rot4.as_quaternion()]) - ori_vec=Orientation(quat,lattice) - - assert ori_vec.inFZ_vec()[0] == ori0.inFZ() - assert ori_vec.inFZ_vec()[1] == ori1.inFZ() - assert ori_vec.inFZ_vec()[2] == ori2.inFZ() - assert ori_vec.inFZ_vec()[3] == ori3.inFZ() - assert ori_vec.inFZ_vec()[4] == ori4.inFZ() @pytest.mark.parametrize('model',['Bain','KS','GT','GT_prime','NW','Pitsch']) From ef0c78745a1b7146bac19c71fd285608da47c194 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 30 Jun 2020 22:33:58 +0200 Subject: [PATCH 33/40] fix for vectorized in_SST + test --- python/damask/_lattice.py | 14 +++++++++----- python/damask/_orientation.py | 1 - python/tests/test_Lattice.py | 19 ++++++++++++++----- python/tests/test_Result.py | 2 +- python/tests/test_Rotation.py | 2 +- 5 files changed, 25 insertions(+), 13 deletions(-) diff --git a/python/damask/_lattice.py b/python/damask/_lattice.py index 74ef61a50..a337aa1a1 100644 --- a/python/damask/_lattice.py +++ b/python/damask/_lattice.py @@ -160,7 +160,7 @@ class Symmetry: Fundamental zone in Rodrigues space is point symmetric around origin. """ if(rho.shape[-1] != 3): - raise ValueError('Input is not a Rodrigues-Frank vector.') + raise ValueError('Input is not a Rodrigues-Frank vector field.') rho_abs = np.abs(rho) @@ -195,7 +195,7 @@ class Symmetry: """ if(rho.shape[-1] != 3): - raise ValueError('Input is not a Rodrigues-Frank vector.') + raise ValueError('Input is not a Rodrigues-Frank vector field.') with np.errstate(invalid='ignore'): # using '*' for 'and' @@ -242,6 +242,9 @@ class Symmetry: ... } """ + if(vector.shape[-1] != 3): + raise ValueError('Input is not a 3D vector field.') + if self.system == 'cubic': basis = {'improper':np.array([ [-1. , 0. , 1. ], [ np.sqrt(2.) , -np.sqrt(2.) , 0. ], @@ -280,16 +283,17 @@ class Symmetry: else: return np.ones_like(vector[...,0],bool) - b_p = np.broadcast_to(basis['proper'], vector.shape+(3,)) + + b_i = np.broadcast_to(basis['improper'],vector.shape+(3,)) if proper: - b_i = np.broadcast_to(basis['improper'],vector.shape+(3,)) + b_p = np.broadcast_to(basis['proper'], vector.shape+(3,)) improper = np.all(np.around(np.einsum('...ji,...i',b_i,vector),12)>=0.0,axis=-1,keepdims=True) theComponents = np.where(np.broadcast_to(improper,vector.shape), np.around(np.einsum('...ji,...i',b_i,vector),12), np.around(np.einsum('...ji,...i',b_p,vector),12)) else: vector_ = np.block([vector[...,0:2],np.abs(vector[...,2:3])]) # z component projects identical - theComponents = np.around(np.einsum('...ji,...i',b_p,vector_),12) + theComponents = np.around(np.einsum('...ji,...i',b_i,vector_),12) in_SST = np.all(theComponents >= 0.0,axis=-1) diff --git a/python/damask/_orientation.py b/python/damask/_orientation.py index b6b24e951..c59ababda 100644 --- a/python/damask/_orientation.py +++ b/python/damask/_orientation.py @@ -79,7 +79,6 @@ class Orientation: # ToDo: make subclass of lattice and Rotation # ... own sym, other sym, # self-->other: True, self<--other: False - def in_FZ(self): """Check if orientations fall into Fundamental Zone.""" return self.lattice.in_FZ(self.rotation.as_Rodrigues(vector=True)) diff --git a/python/tests/test_Lattice.py b/python/tests/test_Lattice.py index 8e4616660..f53c3ad03 100644 --- a/python/tests/test_Lattice.py +++ b/python/tests/test_Lattice.py @@ -106,14 +106,23 @@ class TestSymmetry: @pytest.mark.parametrize('system',Symmetry.crystal_systems) def test_in_FZ_vectorize(self,set_of_rodrigues,system): - for i,in_FZ_ in enumerate(Symmetry(system).in_FZ(set_of_rodrigues)): - assert in_FZ_ == in_FZ(system,set_of_rodrigues[i]) + result = Symmetry(system).in_FZ(set_of_rodrigues.reshape(50,4,3)).reshape(200) + for i,r in enumerate(result): + assert r == in_FZ(system,set_of_rodrigues[i]) @pytest.mark.parametrize('system',Symmetry.crystal_systems) def test_in_disorientation_SST_vectorize(self,set_of_rodrigues,system): - for i,in_disorientation_SST_ in enumerate(Symmetry(system).in_disorientation_SST(set_of_rodrigues)): - assert in_disorientation_SST_ == in_disorientation_SST(system,set_of_rodrigues[i]) + result = Symmetry(system).in_disorientation_SST(set_of_rodrigues.reshape(50,4,3)).reshape(200) + for i,r in enumerate(result): + assert r == in_disorientation_SST(system,set_of_rodrigues[i]) + @pytest.mark.parametrize('proper',[True,False]) + @pytest.mark.parametrize('system',Symmetry.crystal_systems) + def test_in_SST_vectorize(self,system,proper): + vecs = np.random.rand(20,4,3) + result = Symmetry(system).in_SST(vecs,proper).reshape(20*4) + for i,r in enumerate(result): + assert r == in_SST(system,vecs.reshape(20*4,3)[i],proper) @pytest.mark.parametrize('invalid_symmetry',['fcc','bcc','hello']) def test_invalid_symmetry(self,invalid_symmetry): @@ -142,7 +151,7 @@ class TestSymmetry: def test_in_SST(self,system,proper): assert Symmetry(system).in_SST(np.zeros(3),proper) - @pytest.mark.parametrize('function',['in_FZ','in_disorientation_SST']) + @pytest.mark.parametrize('function',['in_FZ','in_disorientation_SST','in_SST']) def test_invalid_argument(self,function): s = Symmetry() # noqa with pytest.raises(ValueError): diff --git a/python/tests/test_Result.py b/python/tests/test_Result.py index aec91db9f..b27a3d5f3 100644 --- a/python/tests/test_Result.py +++ b/python/tests/test_Result.py @@ -319,4 +319,4 @@ class TestResult: def test_XDMF(self,tmp_path,single_phase): os.chdir(tmp_path) - single_phase.write_XDMF + single_phase.write_XDMF() diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py index 20d01af4a..b4166bf1b 100644 --- a/python/tests/test_Rotation.py +++ b/python/tests/test_Rotation.py @@ -20,7 +20,7 @@ def set_of_rotations(set_of_quaternions): #################################################################################################### -# Code below available according to the following conditions on https://github.com/MarDiehl/3Drotations +# Code below available according to the following conditions #################################################################################################### # Copyright (c) 2017-2019, Martin Diehl/Max-Planck-Institut für Eisenforschung GmbH # Copyright (c) 2013-2014, Marc De Graef/Carnegie Mellon University From 23365660d823f432c12e95ec9b9aebd7653ece22 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 30 Jun 2020 23:17:50 +0200 Subject: [PATCH 34/40] polishing --- python/damask/_orientation.py | 1 + python/damask/_rotation.py | 14 +++++++------- python/tests/test_Lattice.py | 1 - 3 files changed, 8 insertions(+), 8 deletions(-) diff --git a/python/damask/_orientation.py b/python/damask/_orientation.py index c59ababda..cb28c160e 100644 --- a/python/damask/_orientation.py +++ b/python/damask/_orientation.py @@ -40,6 +40,7 @@ class Orientation: # ToDo: make subclass of lattice and Rotation self.rotation = Rotation.from_quaternion(rotation) # assume quaternion def __getitem__(self,item): + """Iterate over leading/leftmost dimension of Orientation array.""" return self.__class__(self.rotation[item],self.lattice) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index accd453cc..26fd3e261 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -81,17 +81,16 @@ class Rotation: ]) - def __len__(self): - return 0 if self.shape == () else len(self.shape) - - def __getitem__(self,item): + """Iterate over leading/leftmost dimension of Rotation array.""" + if self.shape == (): return self.copy() if isinstance(item,tuple) and len(item) >= len(self): raise IndexError('Too many indices') return self.__class__(self.quaternion[item]) def __len__(self): + """Length of leading/leftmost dimension of Rotation array.""" return 0 if self.shape == () else self.shape[0] @@ -104,9 +103,10 @@ class Rotation: other : numpy.ndarray or Rotation Vector, second or fourth order tensor, or rotation object that is rotated. - Todo - ---- - Check rotation of 4th order tensor + Returns + ------- + other_rot : numpy.ndarray or Rotation + Rotated vector, second or fourth order tensor, or rotation object. """ if isinstance(other, Rotation): diff --git a/python/tests/test_Lattice.py b/python/tests/test_Lattice.py index f53c3ad03..aa201c645 100644 --- a/python/tests/test_Lattice.py +++ b/python/tests/test_Lattice.py @@ -3,7 +3,6 @@ import random import pytest import numpy as np -from damask import Orientation from damask import Rotation from damask import Symmetry From de8e9b5fc1b6b7305787f15b9c8666260aa09929 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 1 Jul 2020 00:37:02 +0200 Subject: [PATCH 35/40] fast reduced operation --- python/damask/_orientation.py | 50 ++++++++++++++------------------ python/tests/test_Orientation.py | 7 +++++ python/tests/test_Result.py | 2 +- python/tests/test_ori_vec.py | 20 ------------- 4 files changed, 30 insertions(+), 49 deletions(-) diff --git a/python/damask/_orientation.py b/python/damask/_orientation.py index cb28c160e..57407a6ce 100644 --- a/python/damask/_orientation.py +++ b/python/damask/_orientation.py @@ -3,7 +3,7 @@ import numpy as np from . import Lattice from . import Rotation -class Orientation: # ToDo: make subclass of lattice and Rotation +class Orientation: # ToDo: make subclass of lattice and Rotation? """ Crystallographic orientation. @@ -44,6 +44,7 @@ class Orientation: # ToDo: make subclass of lattice and Rotation return self.__class__(self.rotation[item],self.lattice) + # ToDo: Discuss vectorization/calling signature def disorientation(self, other, SST = True, @@ -80,17 +81,19 @@ class Orientation: # ToDo: make subclass of lattice and Rotation # ... own sym, other sym, # self-->other: True, self<--other: False + @property def in_FZ(self): """Check if orientations fall into Fundamental Zone.""" return self.lattice.in_FZ(self.rotation.as_Rodrigues(vector=True)) + @property def equivalent(self): """ Return orientations which are symmetrically equivalent. One dimension (length according to symmetrically equivalent orientations) - is added to the left of the rotation array. + is added to the left of the Rotation array. """ s = self.lattice.symmetry.symmetry_operations @@ -98,9 +101,8 @@ class Orientation: # ToDo: make subclass of lattice and Rotation s = Rotation(np.broadcast_to(s,s.shape[:1]+self.rotation.quaternion.shape)) r = np.broadcast_to(self.rotation.quaternion,s.shape[:1]+self.rotation.quaternion.shape) - r = Rotation(r) - return self.__class__(s@r,self.lattice) + return self.__class__(s@Rotation(r),self.lattice) def relatedOrientations_vec(self,model): @@ -123,34 +125,24 @@ class Orientation: # ToDo: make subclass of lattice and Rotation r = self.lattice.relationOperations(model) return [self.__class__(o*self.rotation,r['lattice']) for o in r['rotations']] + @property - def reduced_vec(self): - """Transform orientation to fall into fundamental zone according to symmetry.""" - equi= self.equivalent.rotation #equivalent orientations - r= 1 if not self.rotation.shape else equi.shape[1] #number of rotations - num_equi=equi.shape[0] #number of equivalente orientations - quat= np.reshape( equi.as_quaternion(), (r*num_equi,4) ,order='F') #equivalents are listed in intiuitive order - boolean=Orientation(quat, self.lattice).in_FZ() #check which ones are in FZ - if sum(boolean) == r: - return self.__class__(quat[boolean],self.lattice) - - else: - print('More than 1 equivalent orientation has been found for an orientation') - index=np.empty(r, dtype=int) - for l,h in enumerate(range(0,r*num_equi, num_equi)): - index[l]=np.where(boolean[h:h+num_equi])[0][0] + (l*num_equi) #get first index that is true then go check to next orientation - - return self.__class__(quat[index],self.lattice) - - def reduced(self): """Transform orientation to fall into fundamental zone according to symmetry.""" - for me in self.equivalent: - if self.lattice.in_FZ(me.rotation.as_Rodrigues(vector=True)): break + eq = self.equivalent + in_FZ = eq.in_FZ - return self.__class__(me.rotation,self.lattice) + # remove duplicates (occur for highly symmetric orientations) + found = np.zeros_like(in_FZ[0],dtype=bool) + q = self.rotation.quaternion[0] + for s in range(in_FZ.shape[0]): + q = np.where(np.expand_dims(np.logical_and(in_FZ[s],~found),-1),eq.rotation.quaternion[s],q) + found = np.logical_or(in_FZ[s],found) + + return self.__class__(q,self.lattice) + # ToDo: vectorize def inversePole(self, axis, proper = False, @@ -159,7 +151,7 @@ class Orientation: # ToDo: make subclass of lattice and Rotation if SST: # pole requested to be within SST for i,o in enumerate(self.equivalent): # test all symmetric equivalent quaternions pole = o.rotation@axis # align crystal direction to axis - if self.lattice.in_SST(pole,proper): break # found SST version + if self.lattice.in_SST(pole,proper): break # found SST version else: pole = self.rotation@axis # align crystal direction to axis @@ -172,7 +164,7 @@ class Orientation: # ToDo: make subclass of lattice and Rotation pole = eq.rotation @ np.broadcast_to(axis/np.linalg.norm(axis),eq.rotation.shape+(3,)) in_SST, color = self.lattice.in_SST(pole,color=True) - # ignore duplicates (occur for highly symmetric orientations) + # remove duplicates (occur for highly symmetric orientations) found = np.zeros_like(in_SST[0],dtype=bool) c = np.empty(color.shape[1:]) for s in range(in_SST.shape[0]): @@ -182,6 +174,7 @@ class Orientation: # ToDo: make subclass of lattice and Rotation return c + # ToDo: Discuss vectorization/calling signature @staticmethod def fromAverage(orientations, weights = []): @@ -202,6 +195,7 @@ class Orientation: # ToDo: make subclass of lattice and Rotation return Orientation(Rotation.fromAverage(closest,weights),ref.lattice) + # ToDo: Discuss vectorization/calling signature def average(self,other): """Calculate the average rotation.""" return Orientation.fromAverage([self,other]) diff --git a/python/tests/test_Orientation.py b/python/tests/test_Orientation.py index f6f25e0a7..78781e47d 100644 --- a/python/tests/test_Orientation.py +++ b/python/tests/test_Orientation.py @@ -54,6 +54,13 @@ class TestOrientation: assert np.allclose(color,IPF_color(oris[i],direction)) + @pytest.mark.parametrize('lattice',Lattice.lattices) + def test_reduced(self,set_of_quaternions,lattice): + oris = Orientation(Rotation(set_of_quaternions),lattice) + reduced = oris.reduced + assert np.all(reduced.in_FZ) and oris.rotation.shape == reduced.rotation.shape + + @pytest.mark.parametrize('model',['Bain','KS','GT','GT_prime','NW','Pitsch']) @pytest.mark.parametrize('lattice',['fcc','bcc']) def test_relationship_forward_backward(self,model,lattice): diff --git a/python/tests/test_Result.py b/python/tests/test_Result.py index b27a3d5f3..676e94d50 100644 --- a/python/tests/test_Result.py +++ b/python/tests/test_Result.py @@ -161,7 +161,7 @@ class TestResult: crystal_structure = default.get_crystal_structure() in_memory = np.empty((qu.shape[0],3),np.uint8) for i,q in enumerate(qu): - o = damask.Orientation(q,crystal_structure).reduced() + o = damask.Orientation(q,crystal_structure).reduced in_memory[i] = np.uint8(o.IPF_color(np.array(d))*255) in_file = default.read_dataset(loc['color']) assert np.allclose(in_memory,in_file) diff --git a/python/tests/test_ori_vec.py b/python/tests/test_ori_vec.py index 7cb465046..9280af164 100644 --- a/python/tests/test_ori_vec.py +++ b/python/tests/test_ori_vec.py @@ -39,23 +39,3 @@ class TestOrientation_vec: ori2.relatedOrientations(model)[s].rotation.as_Rodrigues()) assert all(ori_vec.relatedOrientations_vec(model).rotation.as_cubochoric()[s,3] == \ ori3.relatedOrientations(model)[s].rotation.as_cubochoric()) - - @pytest.mark.parametrize('lattice',Lattice.lattices) - def test_reduced_vec(self,lattice): - ori0=Orientation(rot0,lattice) - ori1=Orientation(rot1,lattice) - ori2=Orientation(rot2,lattice) - ori3=Orientation(rot3,lattice) - #ensure 1 of them is in FZ - ori4=ori0.reduced() - rot4=ori4.rotation - - quat=np.array([rot0.as_quaternion(),rot1.as_quaternion(),\ - rot2.as_quaternion(),rot3.as_quaternion(), rot4.as_quaternion()]) - ori_vec=Orientation(quat,lattice) - - assert all(ori_vec.reduced_vec.rotation.as_Eulers()[0] == ori0.reduced().rotation.as_Eulers() ) - assert all(ori_vec.reduced_vec.rotation.as_quaternion()[1] == ori1.reduced().rotation.as_quaternion() ) - assert all(ori_vec.reduced_vec.rotation.as_Rodrigues()[2] == ori2.reduced().rotation.as_Rodrigues() ) - assert all(ori_vec.reduced_vec.rotation.as_cubochoric()[3] == ori3.reduced().rotation.as_cubochoric() ) - assert all(ori_vec.reduced_vec.rotation.as_axis_angle()[4] == ori4.reduced().rotation.as_axis_angle() ) From e18a5b8a1b2d75c470b880e3e83ad923509e6450 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 1 Jul 2020 08:48:30 +0200 Subject: [PATCH 36/40] simplifications + more tests --- python/damask/_lattice.py | 43 +++++++++--------- python/damask/_orientation.py | 75 ++++++++++++++++---------------- python/tests/test_Orientation.py | 31 ++++++++----- python/tests/test_ori_vec.py | 41 ----------------- 4 files changed, 79 insertions(+), 111 deletions(-) delete mode 100644 python/tests/test_ori_vec.py diff --git a/python/damask/_lattice.py b/python/damask/_lattice.py index a337aa1a1..160b3ba17 100644 --- a/python/damask/_lattice.py +++ b/python/damask/_lattice.py @@ -1,7 +1,5 @@ import numpy as np -from . import Rotation - class Symmetry: """ @@ -29,7 +27,6 @@ class Symmetry: raise KeyError(f'Crystal system "{system}" is unknown') self.system = system.lower() if isinstance(system,str) else system - self.lattice = self.system # ToDo: for compatibility def __copy__(self): @@ -219,6 +216,7 @@ class Symmetry: return np.ones_like(rho[...,0],dtype=bool) + #ToDo: IPF color in separate function def in_SST(self,vector,proper=False,color=False): """ Check whether given vector falls into standard stereographic triangle of own symmetry. @@ -311,9 +309,9 @@ class Symmetry: # ****************************************************************************************** class Lattice: # ToDo: Make a subclass of Symmetry! """ - Lattice system. + Bravais lattice. - Currently, this contains only a mapping from Bravais lattice to symmetry + This contains only a mapping from Bravais lattice to symmetry and orientation relationships. It could include twin and slip systems. References @@ -331,7 +329,7 @@ class Lattice: # ToDo: Make a subclass of Symmetry! } - def __init__(self, lattice): + def __init__(self,lattice,c_over_a=None): """ New lattice of given type. @@ -345,20 +343,20 @@ class Lattice: # ToDo: Make a subclass of Symmetry! self.symmetry = Symmetry(self.lattices[lattice]['system']) # transition to subclass - self.system = self.symmetry.system - self.in_SST = self.symmetry.in_SST - self.in_FZ = self.symmetry.in_FZ - + self.system = self.symmetry.system + self.in_SST = self.symmetry.in_SST + self.in_FZ = self.symmetry.in_FZ + self.in_disorientation_SST = self.symmetry.in_disorientation_SST def __repr__(self): """Report basic lattice information.""" - return 'Bravais lattice {} ({} symmetry)'.format(self.lattice,self.symmetry) + return 'Bravais lattice {} ({} crystal system)'.format(self.lattice,self.symmetry) # Kurdjomov--Sachs orientation relationship for fcc <-> bcc transformation # from S. Morito et al., Journal of Alloys and Compounds 577:s587-s592, 2013 # also see K. Kitahara et al., Acta Materialia 54:1279-1288, 2006 - KS = {'mapping':{'fcc':0,'bcc':1}, + _KS = {'mapping':{'fcc':0,'bcc':1}, 'planes': np.array([ [[ 1, 1, 1],[ 0, 1, 1]], [[ 1, 1, 1],[ 0, 1, 1]], @@ -412,7 +410,7 @@ class Lattice: # ToDo: Make a subclass of Symmetry! # Greninger--Troiano orientation relationship for fcc <-> bcc transformation # from Y. He et al., Journal of Applied Crystallography 39:72-81, 2006 - GT = {'mapping':{'fcc':0,'bcc':1}, + _GT = {'mapping':{'fcc':0,'bcc':1}, 'planes': np.array([ [[ 1, 1, 1],[ 1, 0, 1]], [[ 1, 1, 1],[ 1, 1, 0]], @@ -466,7 +464,7 @@ class Lattice: # ToDo: Make a subclass of Symmetry! # Greninger--Troiano' orientation relationship for fcc <-> bcc transformation # from Y. He et al., Journal of Applied Crystallography 39:72-81, 2006 - GTprime = {'mapping':{'fcc':0,'bcc':1}, + _GTprime = {'mapping':{'fcc':0,'bcc':1}, 'planes': np.array([ [[ 7, 17, 17],[ 12, 5, 17]], [[ 17, 7, 17],[ 17, 12, 5]], @@ -520,7 +518,7 @@ class Lattice: # ToDo: Make a subclass of Symmetry! # Nishiyama--Wassermann orientation relationship for fcc <-> bcc transformation # from H. Kitahara et al., Materials Characterization 54:378-386, 2005 - NW = {'mapping':{'fcc':0,'bcc':1}, + _NW = {'mapping':{'fcc':0,'bcc':1}, 'planes': np.array([ [[ 1, 1, 1],[ 0, 1, 1]], [[ 1, 1, 1],[ 0, 1, 1]], @@ -550,7 +548,7 @@ class Lattice: # ToDo: Make a subclass of Symmetry! # Pitsch orientation relationship for fcc <-> bcc transformation # from Y. He et al., Acta Materialia 53:1179-1190, 2005 - Pitsch = {'mapping':{'fcc':0,'bcc':1}, + _Pitsch = {'mapping':{'fcc':0,'bcc':1}, 'planes': np.array([ [[ 0, 1, 0],[ -1, 0, 1]], [[ 0, 0, 1],[ 1, -1, 0]], @@ -580,7 +578,7 @@ class Lattice: # ToDo: Make a subclass of Symmetry! # Bain orientation relationship for fcc <-> bcc transformation # from Y. He et al., Journal of Applied Crystallography 39:72-81, 2006 - Bain = {'mapping':{'fcc':0,'bcc':1}, + _Bain = {'mapping':{'fcc':0,'bcc':1}, 'planes': np.array([ [[ 1, 0, 0],[ 1, 0, 0]], [[ 0, 1, 0],[ 0, 1, 0]], @@ -590,7 +588,8 @@ class Lattice: # ToDo: Make a subclass of Symmetry! [[ 0, 0, 1],[ 1, 0, 1]], [[ 1, 0, 0],[ 1, 1, 0]]],dtype='float')} - def relationOperations(self,model): + + def relation_operations(self,model): """ Crystallographic orientation relationships for phase transformations. @@ -612,8 +611,8 @@ class Lattice: # ToDo: Make a subclass of Symmetry! https://doi.org/10.1016/j.actamat.2004.11.021 """ - models={'KS':self.KS, 'GT':self.GT, 'GT_prime':self.GTprime, - 'NW':self.NW, 'Pitsch': self.Pitsch, 'Bain':self.Bain} + models={'KS':self._KS, 'GT':self._GT, 'GT_prime':self._GTprime, + 'NW':self._NW, 'Pitsch': self._Pitsch, 'Bain':self._Bain} try: relationship = models[model] except KeyError : @@ -639,6 +638,8 @@ class Lattice: # ToDo: Make a subclass of Symmetry! otherDir = miller[otherDir_id]/ np.linalg.norm(miller[otherDir_id]) otherMatrix = np.array([otherDir,np.cross(otherPlane,otherDir),otherPlane]) - r['rotations'].append(Rotation.from_matrix(np.dot(otherMatrix.T,myMatrix))) + r['rotations'].append(np.dot(otherMatrix.T,myMatrix)) + + r['rotations'] = np.array(r['rotations']) return r diff --git a/python/damask/_orientation.py b/python/damask/_orientation.py index 57407a6ce..3c2d73a92 100644 --- a/python/damask/_orientation.py +++ b/python/damask/_orientation.py @@ -71,8 +71,8 @@ class Orientation: # ToDo: make subclass of lattice and Rotation? r = b*aInv for k in range(2): r.inverse() - breaker = self.in_FZ \ - and (not SST or other.lattice.symmetry.inDisorientationSST(r.as_Rodrigues(vector=True))) + breaker = self.lattice.in_FZ(r.as_Rodrigues(vector=True)) \ + and (not SST or other.lattice.in_disorientation_SST(r.as_Rodrigues(vector=True))) if breaker: break if breaker: break if breaker: break @@ -90,40 +90,36 @@ class Orientation: # ToDo: make subclass of lattice and Rotation? @property def equivalent(self): """ - Return orientations which are symmetrically equivalent. + Orientations which are symmetrically equivalent. - One dimension (length according to symmetrically equivalent orientations) + One dimension (length according to number of symmetrically equivalent orientations) is added to the left of the Rotation array. """ - s = self.lattice.symmetry.symmetry_operations - s = s.reshape(s.shape[:1]+(1,)*len(self.rotation.shape)+(4,)) - s = Rotation(np.broadcast_to(s,s.shape[:1]+self.rotation.quaternion.shape)) + o = self.lattice.symmetry.symmetry_operations + o = o.reshape(o.shape[:1]+(1,)*len(self.rotation.shape)+(4,)) + o = Rotation(np.broadcast_to(o,o.shape[:1]+self.rotation.quaternion.shape)) - r = np.broadcast_to(self.rotation.quaternion,s.shape[:1]+self.rotation.quaternion.shape) + s = np.broadcast_to(self.rotation.quaternion,o.shape[:1]+self.rotation.quaternion.shape) - return self.__class__(s@Rotation(r),self.lattice) + return self.__class__(o@Rotation(s),self.lattice) - def relatedOrientations_vec(self,model): - """List of orientations related by the given orientation relationship.""" - h = self.lattice.relationOperations(model) - rot= h['rotations'] - op=np.array([o.as_quaternion() for o in rot]) + def related(self,model): + """ + Orientations related by the given orientation relationship. - s = op.reshape(op.shape[:1]+(1,)*len(self.rotation.shape)+(4,)) - s = Rotation(np.broadcast_to(s,s.shape[:1]+self.rotation.quaternion.shape)) + One dimension (length according to number of related orientations) + is added to the left of the Rotation array. - r = np.broadcast_to(self.rotation.quaternion,s.shape[:1]+self.rotation.quaternion.shape) - r = Rotation(r) + """ + o = Rotation.from_matrix(self.lattice.relation_operations(model)['rotations']).as_quaternion() + o = o.reshape(o.shape[:1]+(1,)*len(self.rotation.shape)+(4,)) + o = Rotation(np.broadcast_to(o,o.shape[:1]+self.rotation.quaternion.shape)) - return self.__class__(s@r,h['lattice']) + s = np.broadcast_to(self.rotation.quaternion,o.shape[:1]+self.rotation.quaternion.shape) - - def relatedOrientations(self,model): - """List of orientations related by the given orientation relationship.""" - r = self.lattice.relationOperations(model) - return [self.__class__(o*self.rotation,r['lattice']) for o in r['rotations']] + return self.__class__(o@Rotation(s),self.lattice.relation_operations(model)['lattice']) @property @@ -136,26 +132,31 @@ class Orientation: # ToDo: make subclass of lattice and Rotation? found = np.zeros_like(in_FZ[0],dtype=bool) q = self.rotation.quaternion[0] for s in range(in_FZ.shape[0]): + #something fishy... why does q needs to be initialized? q = np.where(np.expand_dims(np.logical_and(in_FZ[s],~found),-1),eq.rotation.quaternion[s],q) found = np.logical_or(in_FZ[s],found) return self.__class__(q,self.lattice) - # ToDo: vectorize - def inversePole(self, - axis, - proper = False, - SST = True): + def inverse_pole(self,axis,proper=False,SST=True): """Axis rotated according to orientation (using crystal symmetry to ensure location falls into SST).""" - if SST: # pole requested to be within SST - for i,o in enumerate(self.equivalent): # test all symmetric equivalent quaternions - pole = o.rotation@axis # align crystal direction to axis - if self.lattice.in_SST(pole,proper): break # found SST version - else: - pole = self.rotation@axis # align crystal direction to axis + if SST: + eq = self.equivalent + pole = eq.rotation @ np.broadcast_to(axis/np.linalg.norm(axis),eq.rotation.shape+(3,)) + in_SST = self.lattice.in_SST(pole,proper=proper) + + # remove duplicates (occur for highly symmetric orientations) + found = np.zeros_like(in_SST[0],dtype=bool) + p = pole[0] + for s in range(in_SST.shape[0]): + p = np.where(np.expand_dims(np.logical_and(in_SST[s],~found),-1),pole[s],p) + found = np.logical_or(in_SST[s],found) + + return p + else: + return self.rotation @ np.broadcast_to(axis/np.linalg.norm(axis),self.rotation.shape+(3,)) - return (pole,i if SST else 0) def IPF_color(self,axis): #ToDo axis or direction? @@ -166,7 +167,7 @@ class Orientation: # ToDo: make subclass of lattice and Rotation? # remove duplicates (occur for highly symmetric orientations) found = np.zeros_like(in_SST[0],dtype=bool) - c = np.empty(color.shape[1:]) + c = color[0] for s in range(in_SST.shape[0]): c = np.where(np.expand_dims(np.logical_and(in_SST[s],~found),-1),color[s],c) found = np.logical_or(in_SST[s],found) diff --git a/python/tests/test_Orientation.py b/python/tests/test_Orientation.py index 78781e47d..4be146ec7 100644 --- a/python/tests/test_Orientation.py +++ b/python/tests/test_Orientation.py @@ -28,6 +28,22 @@ def reference_dir(reference_dir_base): class TestOrientation: + @pytest.mark.parametrize('model',['Bain','KS','GT','GT_prime','NW','Pitsch']) + @pytest.mark.parametrize('lattice',['fcc','bcc']) + def test_relationship_vectorize(self,set_of_quaternions,lattice,model): + result = Orientation(set_of_quaternions[:200].reshape(50,4,4),lattice).related(model) + ref_qu = result.rotation.quaternion.reshape(-1,200,4) + for i in range(200): + single = Orientation(set_of_quaternions[i],lattice).related(model).rotation.quaternion + assert np.allclose(ref_qu[:,i,:],single) + + @pytest.mark.parametrize('lattice',Lattice.lattices) + def test_IPF_vectorize(self,set_of_quaternions,lattice): + direction = np.random.random(3)*2.0-1 + oris = Orientation(Rotation(set_of_quaternions),lattice)[:200] + for i,color in enumerate(oris.IPF_color(direction)): + assert np.allclose(color,IPF_color(oris[i],direction)) + @pytest.mark.parametrize('color',[{'label':'red', 'RGB':[1,0,0],'direction':[0,0,1]}, {'label':'green','RGB':[0,1,0],'direction':[0,1,1]}, {'label':'blue', 'RGB':[0,0,1],'direction':[1,1,1]}]) @@ -45,15 +61,6 @@ class TestOrientation: for equivalent in ori.equivalent: assert np.allclose(color,equivalent.IPF_color(direction)) - - @pytest.mark.parametrize('lattice',Lattice.lattices) - def test_IPF_vectorize(self,set_of_quaternions,lattice): - direction = np.random.random(3)*2.0-1 - oris = Orientation(Rotation(set_of_quaternions),lattice)[:200] - for i,color in enumerate(oris.IPF_color(direction)): - assert np.allclose(color,IPF_color(oris[i],direction)) - - @pytest.mark.parametrize('lattice',Lattice.lattices) def test_reduced(self,set_of_quaternions,lattice): oris = Orientation(Rotation(set_of_quaternions),lattice) @@ -65,8 +72,8 @@ class TestOrientation: @pytest.mark.parametrize('lattice',['fcc','bcc']) def test_relationship_forward_backward(self,model,lattice): ori = Orientation(Rotation.from_random(),lattice) - for i,r in enumerate(ori.relatedOrientations(model)): - ori2 = r.relatedOrientations(model)[i] + for i,r in enumerate(ori.related(model)): + ori2 = r.related(model)[i] misorientation = ori.rotation.misorientation(ori2.rotation) assert misorientation.asAxisAngle(degrees=True)[3]<1.0e-5 @@ -75,7 +82,7 @@ class TestOrientation: def test_relationship_reference(self,update,reference_dir,model,lattice): reference = os.path.join(reference_dir,'{}_{}.txt'.format(lattice,model)) ori = Orientation(Rotation(),lattice) - eu = np.array([o.rotation.as_Eulers(degrees=True) for o in ori.relatedOrientations(model)]) + eu = np.array([o.rotation.as_Eulers(degrees=True) for o in ori.related(model)]) if update: coords = np.array([(1,i+1) for i,x in enumerate(eu)]) table = damask.Table(eu,{'Eulers':(3,)}) diff --git a/python/tests/test_ori_vec.py b/python/tests/test_ori_vec.py deleted file mode 100644 index 9280af164..000000000 --- a/python/tests/test_ori_vec.py +++ /dev/null @@ -1,41 +0,0 @@ -import pytest -import numpy as np - -from damask import Rotation -from damask import Orientation -from damask import Lattice - -rot0= Rotation.from_random() -rot1= Rotation.from_random() -rot2= Rotation.from_random() -rot3= Rotation.from_random() - -#disorientation - -#fromaverage -#average - -class TestOrientation_vec: - - - @pytest.mark.parametrize('model',['Bain','KS','GT','GT_prime','NW','Pitsch']) - @pytest.mark.parametrize('lattice',['fcc','bcc']) - def test_relatedOrientations_vec(self,model,lattice): - ori0=Orientation(rot0,lattice) - ori1=Orientation(rot1,lattice) - ori2=Orientation(rot2,lattice) - ori3=Orientation(rot3,lattice) - - quat=np.array([rot0.as_quaternion(),rot1.as_quaternion(),rot2.as_quaternion(),rot3.as_quaternion()]) - ori_vec=Orientation(quat,lattice) - - - for s in range(len(ori1.lattice.relationOperations(model)['rotations'])): - assert all(ori_vec.relatedOrientations_vec(model).rotation.as_Eulers()[s,0] == \ - ori0.relatedOrientations(model)[s].rotation.as_Eulers()) - assert all(ori_vec.relatedOrientations_vec(model).rotation.as_quaternion()[s,1] == \ - ori1.relatedOrientations(model)[s].rotation.as_quaternion()) - assert all(ori_vec.relatedOrientations_vec(model).rotation.as_Rodrigues()[s,2] == \ - ori2.relatedOrientations(model)[s].rotation.as_Rodrigues()) - assert all(ori_vec.relatedOrientations_vec(model).rotation.as_cubochoric()[s,3] == \ - ori3.relatedOrientations(model)[s].rotation.as_cubochoric()) From be1eb996e01e378440b7d945b517766dc059cfe8 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 1 Jul 2020 15:11:39 +0200 Subject: [PATCH 37/40] more tests and cleaning --- python/damask/_rotation.py | 5 ----- python/tests/test_Orientation.py | 8 +++++++- python/tests/test_Rotation.py | 12 ++++++++++++ 3 files changed, 19 insertions(+), 6 deletions(-) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index 91f126c88..28e2b428f 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -627,12 +627,7 @@ class Rotation: return Rotation(q.reshape(r.shape[:-1]+(4,)) if shape is not None else q)._standardize() - # for compatibility (old names do not follow convention) - def asM(self): - return self.M - fromQuaternion = from_quaternion fromEulers = from_Eulers - asAxisAngle = as_axis_angle __mul__ = __matmul__ #################################################################################################### diff --git a/python/tests/test_Orientation.py b/python/tests/test_Orientation.py index 1c5288923..4a7c3e54d 100644 --- a/python/tests/test_Orientation.py +++ b/python/tests/test_Orientation.py @@ -75,7 +75,7 @@ class TestOrientation: for i,r in enumerate(ori.related(model)): ori2 = r.related(model)[i] misorientation = ori.rotation.misorientation(ori2.rotation) - assert misorientation.asAxisAngle(degrees=True)[3]<1.0e-5 + assert misorientation.as_axis_angle(degrees=True)[3]<1.0e-5 @pytest.mark.parametrize('model',['Bain','KS','GT','GT_prime','NW','Pitsch']) @pytest.mark.parametrize('lattice',['fcc','bcc']) @@ -89,3 +89,9 @@ class TestOrientation: table.add('pos',coords) table.to_ASCII(reference) assert np.allclose(eu,damask.Table.from_ASCII(reference).get('Eulers')) + + @pytest.mark.parametrize('lattice',Lattice.lattices) + def test_disorientation360(self,lattice): + R_1 = Orientation(Rotation(),lattice) + R_2 = Orientation(damask.Rotation.from_Eulers([360,0,0],degrees=True),lattice) + assert np.allclose(R_1.disorientation(R_2).as_matrix(),np.eye(3)) diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py index 0cf2185c0..21d0b5fae 100644 --- a/python/tests/test_Rotation.py +++ b/python/tests/test_Rotation.py @@ -905,3 +905,15 @@ class TestRotation: def test_misorientation(self): R = Rotation.from_random() assert np.allclose(R.misorientation(R).as_matrix(),np.eye(3)) + + def test_misorientation360(self): + R_1 = Rotation() + R_2 = Rotation.from_Eulers([360,0,0],degrees=True) + assert np.allclose(R_1.misorientation(R_2).as_matrix(),np.eye(3)) + + @pytest.mark.parametrize('angle',[10,20,30,40,50,60,70,80,90,100,120]) + def test_average(self,angle): + R_1 = Rotation.from_axis_angle([0,0,1,10],degrees=True) + R_2 = Rotation.from_axis_angle([0,0,1,angle],degrees=True) + avg_angle = R_1.average(R_2).as_axis_angle(degrees=True,pair=True)[1] + assert np.isclose(avg_angle,10+(angle-10)/2.) From 4abd77fccf5f19fde9ffceb9aaab57c292bc6cd6 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 1 Jul 2020 16:26:56 +0200 Subject: [PATCH 38/40] more test coverage --- python/damask/_orientation.py | 8 ++++---- python/damask/_rotation.py | 10 ++++++--- python/tests/test_Orientation.py | 35 ++++++++++++++++++++++++++++++++ 3 files changed, 46 insertions(+), 7 deletions(-) diff --git a/python/damask/_orientation.py b/python/damask/_orientation.py index b3395acfd..0bb0e1bc8 100644 --- a/python/damask/_orientation.py +++ b/python/damask/_orientation.py @@ -61,7 +61,7 @@ class Orientation: # ToDo: make subclass of lattice and Rotation? if self.lattice.symmetry != other.lattice.symmetry: raise NotImplementedError('disorientation between different symmetry classes not supported yet.') - mySymEqs = self.equivalent if SST else self.equivalent[0] #ToDo: This is just me! # take all or only first sym operation + mySymEqs = self.equivalent if SST else self.equivalent[0] #ToDo: This is just me! # take all or only first sym operation otherSymEqs = other.equivalent for i,sA in enumerate(mySymEqs): @@ -177,7 +177,7 @@ class Orientation: # ToDo: make subclass of lattice and Rotation? # ToDo: Discuss vectorization/calling signature @staticmethod - def fromAverage(orientations, + def from_average(orientations, weights = []): """Create orientation from average of list of orientations.""" # further read: Orientation distribution analysis in deformed grains @@ -193,10 +193,10 @@ class Orientation: # ToDo: make subclass of lattice and Rotation? SST = False, # select (o[ther]'s) sym orientation symmetries = True)[2]].rotation) # with lowest misorientation - return Orientation(Rotation.fromAverage(closest,weights),ref.lattice) + return Orientation(Rotation.from_average(closest,weights),ref.lattice) # ToDo: Discuss vectorization/calling signature def average(self,other): """Calculate the average rotation.""" - return Orientation.fromAverage([self,other]) + return Orientation.from_average([self,other]) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index 28e2b428f..674656d4b 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -199,7 +199,7 @@ class Rotation: """ if self.quaternion.shape != (4,) or other.quaternion.shape != (4,): raise NotImplementedError('Support for multiple rotations missing') - return Rotation.fromAverage([self,other]) + return Rotation.from_average([self,other]) ################################################################################################ @@ -294,7 +294,11 @@ class Rotation: """ ro = Rotation._qu2ro(self.quaternion) - return ro[...,:3]*ro[...,3:4] if vector else ro + if vector: + with np.errstate(invalid='ignore'): + return ro[...,:3]*ro[...,3:4] + else: + return ro def as_homochoric(self): """ @@ -577,7 +581,7 @@ class Rotation: @staticmethod - def fromAverage(rotations,weights = None): + def from_average(rotations,weights = None): """ Average rotation. diff --git a/python/tests/test_Orientation.py b/python/tests/test_Orientation.py index 4a7c3e54d..84414a343 100644 --- a/python/tests/test_Orientation.py +++ b/python/tests/test_Orientation.py @@ -20,6 +20,16 @@ def IPF_color(orientation,direction): return color +def inverse_pole(orientation,axis,proper=False,SST=True): + if SST: + for eq in orientation.equivalent: + pole = eq.rotation @ axis/np.linalg.norm(axis) + if orientation.lattice.in_SST(pole,proper=proper): + return pole + else: + return orientation.rotation @ axis/np.linalg.norm(axis) + + @pytest.fixture def reference_dir(reference_dir_base): """Directory containing reference results.""" @@ -44,6 +54,15 @@ class TestOrientation: for i,color in enumerate(oris.IPF_color(direction)): assert np.allclose(color,IPF_color(oris[i],direction)) + @pytest.mark.parametrize('SST',[False,True]) + @pytest.mark.parametrize('proper',[True,False]) + @pytest.mark.parametrize('lattice',Lattice.lattices) + def test_inverse_pole_vectorize(self,set_of_quaternions,lattice,SST,proper): + axis = np.random.random(3)*2.0-1 + oris = Orientation(Rotation(set_of_quaternions),lattice)[:200] + for i,pole in enumerate(oris.inverse_pole(axis,SST=SST)): + assert np.allclose(pole,inverse_pole(oris[i],axis,SST=SST)) + @pytest.mark.parametrize('color',[{'label':'red', 'RGB':[1,0,0],'direction':[0,0,1]}, {'label':'green','RGB':[0,1,0],'direction':[0,1,1]}, {'label':'blue', 'RGB':[0,0,1],'direction':[1,1,1]}]) @@ -95,3 +114,19 @@ class TestOrientation: R_1 = Orientation(Rotation(),lattice) R_2 = Orientation(damask.Rotation.from_Eulers([360,0,0],degrees=True),lattice) assert np.allclose(R_1.disorientation(R_2).as_matrix(),np.eye(3)) + + @pytest.mark.parametrize('lattice',Lattice.lattices) + @pytest.mark.parametrize('angle',[10,20,30,40]) + def test_average(self,angle,lattice): + R_1 = Orientation(Rotation.from_axis_angle([0,0,1,10],degrees=True),lattice) + R_2 = Orientation(Rotation.from_axis_angle([0,0,1,angle],degrees=True),lattice) + avg_angle = R_1.average(R_2).rotation.as_axis_angle(degrees=True,pair=True)[1] + assert np.isclose(avg_angle,10+(angle-10)/2.) + + @pytest.mark.parametrize('lattice',Lattice.lattices) + def test_from_average(self,lattice): + R_1 = Orientation(Rotation.from_random(),lattice) + eqs = [r for r in R_1.equivalent] + R_2 = damask.Orientation.from_average(eqs) + assert np.allclose(R_1.rotation.quaternion,R_2.rotation.quaternion) + From 86dc7054a494f034f22394726f6974ed4aa2f279 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 1 Jul 2020 22:11:16 +0200 Subject: [PATCH 39/40] still needed --- python/damask/_rotation.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index 674656d4b..a605a74b6 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -630,8 +630,9 @@ class Rotation: return Rotation(q.reshape(r.shape[:-1]+(4,)) if shape is not None else q)._standardize() - + # for compatibility (old names do not follow convention) fromEulers = from_Eulers + asAxisAngle = as_axis_angle __mul__ = __matmul__ #################################################################################################### From 208d5109d4371d5c18cfc0fac6c617d04c9012f3 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 2 Jul 2020 08:14:13 +0200 Subject: [PATCH 40/40] still needed ... --- python/damask/_rotation.py | 1 + 1 file changed, 1 insertion(+) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index a605a74b6..59bca3795 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -632,6 +632,7 @@ class Rotation: # for compatibility (old names do not follow convention) fromEulers = from_Eulers + fromQuaternion = from_quaternion asAxisAngle = as_axis_angle __mul__ = __matmul__