Merge remote-tracking branch 'origin/development' into Fortran-simplifications

PRIVATE

@@ -1 +1 @@
-Subproject commit 65ec74c07052e77f35a4b5e80bf110aff1f5ae61
+Subproject commit 555f3e01f2b5cf43ade1bd48423b890adca21771

VERSION

@@ -1 +1 @@
-v3.0.0-alpha-216-gb51a1b7b
+v3.0.0-alpha-233-g190f8a82

python/damask/_rotation.py

@@ -438,6 +438,7 @@ class Rotation:
         if np.any(ax[...,3] < 0.0) or np.any(ax[...,3] > np.pi):
             raise ValueError('Axis angle rotation angle outside of [0..π].')
         if not np.all(np.isclose(np.linalg.norm(ax[...,0:3],axis=-1), 1.0)):
+            print(np.linalg.norm(ax[...,0:3],axis=-1))
             raise ValueError('Axis angle rotation axis is not of unit length.')
 
         return Rotation(Rotation._ax2qu(ax))
@@ -652,6 +653,84 @@ class Rotation:
     asAxisAngle    = as_axis_angle
     __mul__        = __matmul__
 
+
+    @staticmethod
+    def from_spherical_component(center,sigma,N=500,degrees=True,seed=None):
+        """
+        Calculate set of rotations with Gaussian distribution around center.
+
+        Parameters
+        ----------
+        center : Rotation
+            Central Rotation.
+        sigma : float
+            Standard deviation of (Gaussian) misorientation distribution.
+        N : int, optional
+            Number of samples, defaults to 500.
+        degrees : boolean, optional
+            sigma is given in degrees.
+        seed : {None, int, array_like[ints], SeedSequence, BitGenerator, Generator}, optional
+            A seed to initialize the BitGenerator. Defaults to None, i.e. unpredictable entropy
+            will be pulled from the OS.
+
+        """
+        rng = np.random.default_rng(seed)
+        sigma = np.radians(sigma) if degrees else sigma
+        u,Theta  = (rng.random((N,2)) * 2.0 * np.array([1,np.pi]) - np.array([1.0, 0])).T
+        omega = abs(rng.normal(scale=sigma,size=N))
+        p = np.column_stack([np.sqrt(1-u**2)*np.cos(Theta),
+                             np.sqrt(1-u**2)*np.sin(Theta),
+                             u, omega])
+
+        return  Rotation.from_axis_angle(p) @ center
+
+
+    @staticmethod
+    def from_fiber_component(alpha,beta,sigma=0.0,N=500,degrees=True,seed=None):
+        """
+        Calculate set of rotations with Gaussian distribution around direction.
+
+        Parameters
+        ----------
+        alpha : numpy.ndarray of size 2
+            Polar coordinates (phi from x,theta from z) of fiber direction in crystal frame.
+        beta : numpy.ndarray of size 2
+            Polar coordinates (phi from x,theta from z) of fiber direction in sample frame.
+        sigma : float, optional
+            Standard deviation of (Gaussian) misorientation distribution.
+            Defaults to 0.
+        N : int, optional
+            Number of samples, defaults to 500.
+        degrees : boolean, optional
+            sigma, alpha, and beta are given in degrees.
+        seed : {None, int, array_like[ints], SeedSequence, BitGenerator, Generator}, optional
+            A seed to initialize the BitGenerator. Defaults to None, i.e. unpredictable entropy
+            will be pulled from the OS.
+
+        """
+        rng = np.random.default_rng(seed)
+        sigma_,alpha_,beta_ = map(np.radians,(sigma,alpha,beta)) if degrees else (sigma,alpha,beta)
+
+        d_cr  = np.array([np.sin(alpha_[0])*np.cos(alpha_[1]), np.sin(alpha_[0])*np.sin(alpha_[1]), np.cos(alpha_[0])])
+        d_lab = np.array([np.sin( beta_[0])*np.cos( beta_[1]), np.sin( beta_[0])*np.sin( beta_[1]), np.cos( beta_[0])])
+        ax_align = np.append(np.cross(d_lab,d_cr), np.arccos(np.dot(d_lab,d_cr)))
+        if np.isclose(ax_align[3],0.0): ax_align[:3] = np.array([1,0,0])
+        R_align  = Rotation.from_axis_angle(ax_align if ax_align[3] > 0.0 else -ax_align,normalize=True) # rotate fiber axis from sample to crystal frame
+
+        u,Theta  = (rng.random((N,2)) * 2.0 * np.array([1,np.pi]) - np.array([1.0, 0])).T
+        omega  = abs(rng.normal(scale=sigma_,size=N))
+        p = np.column_stack([np.sqrt(1-u**2)*np.cos(Theta),
+                             np.sqrt(1-u**2)*np.sin(Theta),
+                             u, omega])
+        p[:,:3] = np.einsum('ij,...j',np.eye(3)-np.outer(d_lab,d_lab),p[:,:3])                      # remove component along fiber axis
+        f = np.column_stack((np.broadcast_to(d_lab,(N,3)),rng.random(N)*np.pi))
+        f[::2,:3] *= -1                                                                             # flip half the rotation axes to negative sense
+
+        return R_align.broadcast_to(N) \
+             @ Rotation.from_axis_angle(p,normalize=True) \
+             @ Rotation.from_axis_angle(f)
+
+
 ####################################################################################################
 # Code below available according to the following conditions on https://github.com/MarDiehl/3Drotations
 ####################################################################################################

python/tests/test_Rotation.py

@@ -2,6 +2,7 @@ import os
 
 import pytest
 import numpy as np
+from scipy import stats
 
 from damask import Rotation
 from damask import _rotation
@@ -920,3 +921,39 @@ class TestRotation:
         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.)
+
+
+    @pytest.mark.parametrize('sigma',[5,10,15,20])
+    @pytest.mark.parametrize('N',[1000,10000,100000])
+    def test_spherical_component(self,N,sigma):
+        c = Rotation.from_random()
+        o = Rotation.from_spherical_component(c,sigma,N)
+        _, angles = c.misorientation(o).as_axis_angle(pair=True,degrees=True)
+        angles[::2] *= -1                                                                           # flip angle for every second to symmetrize distribution
+
+        p = stats.normaltest(angles)[1]
+        sigma_out = np.std(angles)
+        print(f'\np: {p}, sigma ratio {sigma/sigma_out}')
+        assert (.9 < sigma/sigma_out < 1.1) and p > 0.001
+
+
+    @pytest.mark.parametrize('sigma',[5,10,15,20])
+    @pytest.mark.parametrize('N',[1000,10000,100000])
+    def test_from_fiber_component(self,N,sigma):
+        """https://en.wikipedia.org/wiki/Full_width_at_half_maximum."""
+        alpha = np.random.random(2)*np.pi
+        beta  = np.random.random(2)*np.pi
+
+        f_in_C = np.array([np.sin(alpha[0])*np.cos(alpha[1]), np.sin(alpha[0])*np.sin(alpha[1]), np.cos(alpha[0])])
+        f_in_S = np.array([np.sin(beta[0] )*np.cos(beta[1] ), np.sin(beta[0] )*np.sin(beta[1] ), np.cos(beta[0] )])
+        ax = np.append(np.cross(f_in_C,f_in_S), - np.arccos(np.dot(f_in_C,f_in_S)))
+        n = Rotation.from_axis_angle(ax if ax[3] > 0.0 else ax*-1.0 ,normalize=True)                # rotation to align fiber axis in crystal and sample system
+
+        o = Rotation.from_fiber_component(alpha,beta,np.radians(sigma),N,False)
+        angles = np.arccos(np.clip(np.dot(o@np.broadcast_to(f_in_S,(N,3)),n@f_in_S),-1,1))
+        dist   = np.array(angles) * (np.random.randint(0,2,N)*2-1)
+
+        p = stats.normaltest(dist)[1]
+        sigma_out = np.degrees(np.std(dist))
+        print(f'\np: {p}, sigma ratio {sigma/sigma_out}')
+        assert (.9 < sigma/sigma_out < 1.1) and p > 0.001

src/marc/discretization_marc.f90

@@ -723,7 +723,7 @@ subroutine inputRead_microstructure(microstructureAt,&
     endif
   enddo
 
-  if(any(microstructureAt < 1)) call IO_error(190)
+  if(any(microstructureAt < 1)) call IO_error(180)
 
 end subroutine inputRead_microstructure