From d6378ec9bcf5319a6561443f0483e52b8c957d1e Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Tue, 15 Sep 2020 16:34:19 -0400 Subject: [PATCH] shortened from_fiber_component algorithm --- python/damask/_rotation.py | 66 +++++++++++++------------------------- 1 file changed, 23 insertions(+), 43 deletions(-) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index 425e7aeef..8950ee38c 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -434,10 +434,11 @@ class Rotation: if P == 1: ax[...,0:3] *= -1 if degrees: ax[..., 3] = np.radians(ax[...,3]) - if normalize: ax[...,0:3] /= np.linalg.norm(ax[...,0:3],axis=-1) + if normalize: ax[...,0:3] /= np.linalg.norm(ax[...,0:3],axis=-1,keepdims=True) 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)) @@ -516,7 +517,7 @@ class Rotation: raise ValueError('P ∉ {-1,1}') if P == 1: ro[...,0:3] *= -1 - if normalize: ro[...,0:3] /= np.linalg.norm(ro[...,0:3],axis=-1) + if normalize: ro[...,0:3] /= np.linalg.norm(ro[...,0:3],axis=-1,keepdims=True) if np.any(ro[...,3] < 0.0): raise ValueError('Rodrigues vector rotation angle not positive.') if not np.all(np.isclose(np.linalg.norm(ro[...,0:3],axis=-1), 1.0)): @@ -686,7 +687,7 @@ class Rotation: @staticmethod - def from_fiber_component(alpha,beta,FWHM,N_samples=500,degrees=True,seed=None): + def from_fiber_component(alpha,beta,FWHM=0.0,N=500,degrees=True,seed=None): """ Calculate set of rotations with Gaussian distribution around direction. @@ -701,8 +702,9 @@ class Rotation: tbd. beta : numpy.ndarray of size 2 tbd. - FWHM : float + FWHM : float, optional Full width at half maximum of the Gaussian distribution. + Defaults to 0. N_samples : int, optional Number of samples, defaults to 500. degrees : boolean, optional @@ -713,48 +715,26 @@ class Rotation: """ rng = np.random.default_rng(seed) - FWHM_,alpha_,beta_ = np.radians((FWHM,alpha,beta)) if degrees else (FWHM,alpha,beta) + FWHM_,alpha_,beta_ = map(np.radians,(FWHM,alpha,beta)) if degrees else (FWHM,alpha,beta) - 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))) - R_align = 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 + 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) # rotation to align fiber axis in crystal and sample system - rotations = [] - for i in range(N_samples): - rnd = rng.random(3) - ax = np.append(f_in_S,rnd[0]*np.pi*2.0) - if ax[3] > np.pi: ax = np.append(ax[:3]*-1,np.pi*2-ax[3]) - R = R_align @ Rotation.from_axis_angle(ax) # rotation (0..360deg) perpendicular to fiber axis + u,v,b = (np.random.random((N,3)) * 2 * np.array([1,np.pi,np.pi]) - np.array([1,0,np.pi])).T + a = abs(np.random.normal(scale=FWHM_,size=N)) + p = np.vstack((np.sqrt(1-u**2)*np.cos(v), + np.sqrt(1-u**2)*np.sin(v), + u, + a)).T - if FWHM_ > np.radians(0.1): - - i_smallest = np.argmin(np.abs(f_in_S)) - i_non_smallest = list(filter(lambda x: x!=i_smallest, [0,1,2])) - - s = f_in_S[i_smallest] - a = f_in_S[i_non_smallest] - x = np.sum(a**2) - - u = np.empty(3) - while True: # rejection sampling - angle = (rnd[1] - 0.5)*4 *FWHM_ - # solve cos(angle) = dot_product(fInS,u) for u. This is underdetermined, hence assume that - # they share the smallest component. - c = np.cos(angle) - s**2 - u[i_non_smallest[1]] = -(2.0*c*a[1] + np.sqrt(4*((c*a[1])**2.0-x*(c**2.0-a[0]**2*(1.0-s**2)))))/(2*x) - u[i_non_smallest[0]] = np.sqrt(1.0-u[i_non_smallest[1]]**2.0-s**2.0) - u[i_smallest] = s - - if (rnd[2] <= np.exp(-4.0*np.log(2.0)*(angle/FWHM_)**2)): - ax = np.append(np.cross(u,f_in_S),angle) - R = R @ Rotation.from_axis_angle(ax if ax[3] > 0.0 else ax*-1.0,normalize=True) # tilt around direction of smallest component - break - else: - rnd = rng.random(3) - rotations.append(R.as_quaternion()) - - return Rotation.from_quaternion(rotations) + f = np.hstack((np.broadcast_to(d_cr,(N,3)),b.reshape(N,1))) + f[f[:,3]<0] *= -1. + return Rotation.from_axis_angle(p) \ + * Rotation.from_axis_angle(f) \ + * R_align