diff --git a/PRIVATE b/PRIVATE index 65ec74c07..555f3e01f 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 65ec74c07052e77f35a4b5e80bf110aff1f5ae61 +Subproject commit 555f3e01f2b5cf43ade1bd48423b890adca21771 diff --git a/VERSION b/VERSION index fa691fc66..bebbc79d0 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v3.0.0-alpha-216-gb51a1b7b +v3.0.0-alpha-233-g190f8a82 diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index 461f879af..d9237c7e2 100644 --- a/python/damask/_rotation.py +++ b/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 #################################################################################################### diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py index 1696ef247..5435895a2 100644 --- a/python/tests/test_Rotation.py +++ b/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 diff --git a/src/marc/discretization_marc.f90 b/src/marc/discretization_marc.f90 index 30d1bd9fc..cbfc40995 100644 --- a/src/marc/discretization_marc.f90 +++ b/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