From 753fbb70fd7401745f9a3259b5c124affd5e0df6 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 16 Jun 2020 13:27:21 +0200 Subject: [PATCH 01/13] cleaning --- src/discretization.f90 | 2 +- src/mesh/discretization_mesh.f90 | 146 ++++++++++++++----------------- 2 files changed, 67 insertions(+), 81 deletions(-) diff --git a/src/discretization.f90 b/src/discretization.f90 index e9a987912..16f76b158 100644 --- a/src/discretization.f90 +++ b/src/discretization.f90 @@ -49,7 +49,7 @@ subroutine discretization_init(homogenizationAt,microstructureAt,& IPcoords0, & NodeCoords0 integer, optional, intent(in) :: & - sharedNodesBegin + sharedNodesBegin ! index of first node shared among different processes (MPI) write(6,'(/,a)') ' <<<+- discretization init -+>>>'; flush(6) diff --git a/src/mesh/discretization_mesh.f90 b/src/mesh/discretization_mesh.f90 index 0880de115..329205f77 100644 --- a/src/mesh/discretization_mesh.f90 +++ b/src/mesh/discretization_mesh.f90 @@ -8,53 +8,50 @@ module discretization_mesh #include #include #include - use PETScdmplex - use PETScdmda - use PETScis + use PETScdmplex + use PETScdmda + use PETScis - use DAMASK_interface - use IO - use debug - use discretization - use numerics - use FEsolving - use FEM_quadrature - use prec - - implicit none - private - - integer, public, protected :: & - mesh_Nboundaries, & - mesh_NcpElemsGlobal + use DAMASK_interface + use IO + use debug + use discretization + use numerics + use FEsolving + use FEM_quadrature + use prec - integer :: & - mesh_NcpElems !< total number of CP elements in mesh + implicit none + private + + integer, public, protected :: & + mesh_Nboundaries, & + mesh_NcpElemsGlobal + + integer :: & + mesh_NcpElems !< total number of CP elements in mesh !!!! BEGIN DEPRECATED !!!!! - integer, public, protected :: & - mesh_maxNips !< max number of IPs in any CP element + integer, public, protected :: & + mesh_maxNips !< max number of IPs in any CP element !!!! BEGIN DEPRECATED !!!!! - integer, dimension(:,:), allocatable :: & - mesh_element !DEPRECATED + real(pReal), dimension(:,:), allocatable :: & + mesh_ipVolume, & !< volume associated with IP (initially!) + mesh_node0 !< node x,y,z coordinates (initially!) - real(pReal), dimension(:,:), allocatable :: & - mesh_ipVolume, & !< volume associated with IP (initially!) - mesh_node0 !< node x,y,z coordinates (initially!) - - real(pReal), dimension(:,:,:), allocatable :: & - mesh_ipCoordinates !< IP x,y,z coordinates (after deformation!) + real(pReal), dimension(:,:,:), allocatable :: & + mesh_ipCoordinates !< IP x,y,z coordinates (after deformation!) - DM, public :: geomMesh - - PetscInt, dimension(:), allocatable, public, protected :: & - mesh_boundaries + DM, public :: geomMesh - public :: & - discretization_mesh_init, & - mesh_FEM_build_ipVolumes, & - mesh_FEM_build_ipCoordinates + PetscInt, dimension(:), allocatable, public, protected :: & + mesh_boundaries + + public :: & + discretization_mesh_init, & + mesh_FEM_build_ipVolumes, & + mesh_FEM_build_ipCoordinates contains @@ -67,39 +64,32 @@ subroutine discretization_mesh_init(restart) logical, intent(in) :: restart - integer, dimension(1), parameter:: FE_geomtype = [1] !< geometry type of particular element type - integer, dimension(1) :: FE_Nips !< number of IPs in a specific type of element - integer, allocatable, dimension(:) :: chunkPos integer :: dimPlex, & mesh_Nnodes, & !< total number of nodes in mesh j, l - integer, parameter :: & - mesh_ElemType=1 !< Element type of the mesh (only support homogeneous meshes) PetscSF :: sf DM :: globalMesh PetscInt :: nFaceSets PetscInt, pointer, dimension(:) :: pFaceSets character(len=pStringLen), dimension(:), allocatable :: fileContent - IS :: faceSetIS + IS :: faceSetIS PetscErrorCode :: ierr + integer, dimension(:,:), allocatable :: & + mesh_element + - write(6,'(/,a)') ' <<<+- mesh init -+>>>' - ! read in file call DMPlexCreateFromFile(PETSC_COMM_WORLD,geometryFile,PETSC_TRUE,globalMesh,ierr) CHKERRQ(ierr) - ! get spatial dimension (2 or 3?) call DMGetDimension(globalMesh,dimPlex,ierr) CHKERRQ(ierr) - write(6,*) 'dimension',dimPlex;flush(6) call DMGetStratumSize(globalMesh,'depth',dimPlex,mesh_NcpElemsGlobal,ierr) CHKERRQ(ierr) ! get number of IDs in face sets (for boundary conditions?) call DMGetLabelSize(globalMesh,'Face Sets',mesh_Nboundaries,ierr) CHKERRQ(ierr) - write(6,*) 'number of "Face Sets"',mesh_Nboundaries;flush(6) call MPI_Bcast(mesh_Nboundaries,1,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr) call MPI_Bcast(mesh_NcpElemsGlobal,1,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr) call MPI_Bcast(dimPlex,1,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr) @@ -142,48 +132,44 @@ subroutine discretization_mesh_init(restart) enddo call DMClone(globalMesh,geomMesh,ierr) CHKERRQ(ierr) - else + else call DMPlexDistribute(globalMesh,0,sf,geomMesh,ierr) CHKERRQ(ierr) - endif + endif call DMDestroy(globalMesh,ierr); CHKERRQ(ierr) - + call DMGetStratumSize(geomMesh,'depth',dimPlex,mesh_NcpElems,ierr) CHKERRQ(ierr) call DMGetStratumSize(geomMesh,'depth',0,mesh_Nnodes,ierr) CHKERRQ(ierr) - FE_Nips(FE_geomtype(1)) = FEM_nQuadrature(dimPlex,integrationOrder) - mesh_maxNips = FE_Nips(1) - - write(6,*) 'mesh_maxNips',mesh_maxNips + mesh_maxNips = FEM_nQuadrature(dimPlex,integrationOrder) + call mesh_FEM_build_ipCoordinates(dimPlex,FEM_quadrature_points(dimPlex,integrationOrder)%p) call mesh_FEM_build_ipVolumes(dimPlex) - + allocate (mesh_element (4,mesh_NcpElems)); mesh_element = 0 do j = 1, mesh_NcpElems - mesh_element( 1,j) = -1 ! DEPRECATED - mesh_element( 2,j) = mesh_elemType ! elem type + mesh_element( 1,j) = -1 + mesh_element( 2,j) = -1 mesh_element( 3,j) = 1 ! homogenization call DMGetLabelValue(geomMesh,'material',j-1,mesh_element(4,j),ierr) CHKERRQ(ierr) - end do + end do + + if (debug_e < 1 .or. debug_e > mesh_NcpElems) call IO_error(602,ext_msg='element') + if (debug_i < 1 .or. debug_i > mesh_maxNips) call IO_error(602,ext_msg='IP') + + FEsolving_execElem = [1,mesh_NcpElems] ! parallel loop bounds set to comprise all DAMASK elements + FEsolving_execIP = [1,mesh_maxNips] - if (debug_e < 1 .or. debug_e > mesh_NcpElems) & - call IO_error(602,ext_msg='element') ! selected element does not exist - if (debug_i < 1 .or. debug_i > FE_Nips(FE_geomtype(mesh_element(2,debug_e)))) & - call IO_error(602,ext_msg='IP') ! selected element does not have requested IP - - FEsolving_execElem = [ 1,mesh_NcpElems ] ! parallel loop bounds set to comprise all DAMASK elements - FEsolving_execIP = [1,FE_Nips(FE_geomtype(mesh_element(2,1)))] - allocate(mesh_node0(3,mesh_Nnodes),source=0.0_pReal) call discretization_init(mesh_element(3,:),mesh_element(4,:),& reshape(mesh_ipCoordinates,[3,mesh_maxNips*mesh_NcpElems]), & mesh_node0) - + end subroutine discretization_mesh_init @@ -191,23 +177,23 @@ end subroutine discretization_mesh_init !> @brief Calculates IP volume. Allocates global array 'mesh_ipVolume' !-------------------------------------------------------------------------------------------------- subroutine mesh_FEM_build_ipVolumes(dimPlex) - + PetscInt,intent(in):: dimPlex PetscReal :: vol PetscReal, pointer,dimension(:) :: pCent, pNorm PetscInt :: cellStart, cellEnd, cell PetscErrorCode :: ierr - + allocate(mesh_ipVolume(mesh_maxNips,mesh_NcpElems),source=0.0_pReal) - + call DMPlexGetHeightStratum(geomMesh,0,cellStart,cellEnd,ierr); CHKERRQ(ierr) allocate(pCent(dimPlex)) allocate(pNorm(dimPlex)) do cell = cellStart, cellEnd-1 - call DMPlexComputeCellGeometryFVM(geomMesh,cell,vol,pCent,pNorm,ierr) + call DMPlexComputeCellGeometryFVM(geomMesh,cell,vol,pCent,pNorm,ierr) CHKERRQ(ierr) mesh_ipVolume(:,cell+1) = vol/real(mesh_maxNips,pReal) - enddo + enddo end subroutine mesh_FEM_build_ipVolumes @@ -219,20 +205,20 @@ subroutine mesh_FEM_build_ipCoordinates(dimPlex,qPoints) PetscInt, intent(in) :: dimPlex PetscReal, intent(in) :: qPoints(mesh_maxNips*dimPlex) - + PetscReal, pointer,dimension(:) :: pV0, pCellJ, pInvcellJ PetscReal :: detJ PetscInt :: cellStart, cellEnd, cell, qPt, dirI, dirJ, qOffset PetscErrorCode :: ierr - - + + allocate(mesh_ipCoordinates(3,mesh_maxNips,mesh_NcpElems),source=0.0_pReal) - + allocate(pV0(dimPlex)) allocatE(pCellJ(dimPlex**2)) allocatE(pinvCellJ(dimPlex**2)) call DMPlexGetHeightStratum(geomMesh,0,cellStart,cellEnd,ierr); CHKERRQ(ierr) - do cell = cellStart, cellEnd-1 !< loop over all elements + do cell = cellStart, cellEnd-1 !< loop over all elements call DMPlexComputeCellGeometryAffineFEM(geomMesh,cell,pV0,pCellJ,pInvcellJ,detJ,ierr) CHKERRQ(ierr) qOffset = 0 @@ -246,7 +232,7 @@ subroutine mesh_FEM_build_ipCoordinates(dimPlex,qPoints) enddo qOffset = qOffset + dimPlex enddo - enddo + enddo end subroutine mesh_FEM_build_ipCoordinates From e0d4ee44a340a7adf327732c3fb9e921ee917b2a Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 16 Jun 2020 13:59:59 +0200 Subject: [PATCH 02/13] better name --- src/mesh/discretization_mesh.f90 | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/src/mesh/discretization_mesh.f90 b/src/mesh/discretization_mesh.f90 index 329205f77..f462f7196 100644 --- a/src/mesh/discretization_mesh.f90 +++ b/src/mesh/discretization_mesh.f90 @@ -75,8 +75,9 @@ subroutine discretization_mesh_init(restart) character(len=pStringLen), dimension(:), allocatable :: fileContent IS :: faceSetIS PetscErrorCode :: ierr - integer, dimension(:,:), allocatable :: & - mesh_element + integer, dimension(:), allocatable :: & + homogenizationAt, & + microstructureAt write(6,'(/,a)') ' <<<+- mesh init -+>>>' @@ -149,12 +150,10 @@ subroutine discretization_mesh_init(restart) call mesh_FEM_build_ipCoordinates(dimPlex,FEM_quadrature_points(dimPlex,integrationOrder)%p) call mesh_FEM_build_ipVolumes(dimPlex) - allocate (mesh_element (4,mesh_NcpElems)); mesh_element = 0 + allocate(microstructureAt(mesh_NcpElems)) + allocate(homogenizationAt(mesh_NcpElems),source=1) do j = 1, mesh_NcpElems - mesh_element( 1,j) = -1 - mesh_element( 2,j) = -1 - mesh_element( 3,j) = 1 ! homogenization - call DMGetLabelValue(geomMesh,'material',j-1,mesh_element(4,j),ierr) + call DMGetLabelValue(geomMesh,'material',j-1,microstructureAt(j),ierr) CHKERRQ(ierr) end do @@ -166,7 +165,7 @@ subroutine discretization_mesh_init(restart) allocate(mesh_node0(3,mesh_Nnodes),source=0.0_pReal) - call discretization_init(mesh_element(3,:),mesh_element(4,:),& + call discretization_init(microstructureAt,homogenizationAt,& reshape(mesh_ipCoordinates,[3,mesh_maxNips*mesh_NcpElems]), & mesh_node0) From d0d9245707a665952042db5a2d2090af8a574939 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 3 Jun 2020 15:02:39 +0200 Subject: [PATCH 03/13] clearer intention --- python/tests/test_Result.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/tests/test_Result.py b/python/tests/test_Result.py index d7946e5e0..38c0d84cc 100644 --- a/python/tests/test_Result.py +++ b/python/tests/test_Result.py @@ -302,8 +302,8 @@ class TestResult: @pytest.mark.parametrize('allowed',['off','on']) def test_rename(self,default,allowed): - F = default.read_dataset(default.get_dataset_location('F')) if allowed == 'on': + F = default.read_dataset(default.get_dataset_location('F')) default.allow_modification() default.rename('F','new_name') assert np.all(F == default.read_dataset(default.get_dataset_location('new_name'))) From 829896390ca3089d2f9e2e1d7d9d3b2659f2878e Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 16 Jun 2020 23:37:24 +0200 Subject: [PATCH 04/13] hopefully not needed any more --- src/crystallite.f90 | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 929fec862..e85f9e517 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -287,11 +287,9 @@ end subroutine crystallite_init !-------------------------------------------------------------------------------------------------- !> @brief calculate stress (P) !-------------------------------------------------------------------------------------------------- -function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC) +function crystallite_stress() logical, dimension(discretization_nIP,discretization_nElem) :: crystallite_stress - real(pReal), intent(in), optional :: & - dummyArgumentToPreventInternalCompilerErrorWithGCC real(pReal) :: & formerSubStep integer :: & From 62384b583677347de327158d431511f5ee6347da Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 16 Jun 2020 23:43:15 +0200 Subject: [PATCH 05/13] bugfix: invalid description/unit --- src/constitutive_plastic_dislotwin.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/constitutive_plastic_dislotwin.f90 b/src/constitutive_plastic_dislotwin.f90 index 7c7d24ab8..dc40d269e 100644 --- a/src/constitutive_plastic_dislotwin.f90 +++ b/src/constitutive_plastic_dislotwin.f90 @@ -830,7 +830,7 @@ module subroutine plastic_dislotwin_results(instance,group) 'mobile dislocation density','1/m²') case('rho_dip') if(prm%sum_N_sl>0) call results_writeDataset(group,stt%rho_dip,'rho_dip',& - 'dislocation dipole density''1/m²') + 'dislocation dipole density','1/m²') case('gamma_sl') if(prm%sum_N_sl>0) call results_writeDataset(group,stt%gamma_sl,'gamma_sl',& 'plastic shear','1') From 8fa68101b8e8c13774c4b28be67b3f49dcefd65f Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 17 Jun 2020 22:28:29 +0200 Subject: [PATCH 06/13] not needed --- examples/ConfigFiles/Phase_DisloUCLA_Tungsten.config | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/ConfigFiles/Phase_DisloUCLA_Tungsten.config b/examples/ConfigFiles/Phase_DisloUCLA_Tungsten.config index 6920aee2b..2cdd5c58f 100644 --- a/examples/ConfigFiles/Phase_DisloUCLA_Tungsten.config +++ b/examples/ConfigFiles/Phase_DisloUCLA_Tungsten.config @@ -27,7 +27,7 @@ SolidSolutionStrength 0.0 # Strength due to elements in solid solution ### Dislocation glide parameters ### #per family -Nslip 12 0 +Nslip 12 slipburgers 2.72e-10 # Burgers vector of slip system [m] rhoedge0 1.0e12 # Initial edge dislocation density [m/m**3] rhoedgedip0 1.0 # Initial edged dipole dislocation density [m/m**3] From 055fa64f5f0270feb4f41e4384586b6a70928f57 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 19 Jun 2020 12:25:46 +0200 Subject: [PATCH 07/13] better readable --- python/damask/_rotation.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index b1b1cd5ad..7e9a03c3a 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -530,7 +530,7 @@ class Rotation: np.zeros(qu.shape[:-1]+(2,))]), np.where(np.abs(q03_s) < 1.0e-8, np.block([np.arctan2( 2.0*qu[...,1:2]*qu[...,2:3],qu[...,1:2]**2-qu[...,2:3]**2), - np.broadcast_to(np.pi,qu.shape[:-1]+(1,)), + np.broadcast_to(np.pi,qu[...,0:1].shape), np.zeros(qu.shape[:-1]+(1,))]), np.block([np.arctan2((-_P*q02+q13)*chi, (-_P*q01-q23)*chi), np.arctan2( 2.0*chi, q03_s-q12_s ), @@ -553,7 +553,7 @@ class Rotation: s = np.sign(qu[...,0:1])/np.sqrt(qu[...,1:2]**2+qu[...,2:3]**2+qu[...,3:4]**2) omega = 2.0 * np.arccos(np.clip(qu[...,0:1],-1.0,1.0)) ax = np.where(np.broadcast_to(qu[...,0:1] < 1.0e-8,qu.shape), - np.block([qu[...,1:4],np.broadcast_to(np.pi,qu.shape[:-1]+(1,))]), + np.block([qu[...,1:4],np.broadcast_to(np.pi,qu[...,0:1].shape)]), np.block([qu[...,1:4]*s,omega])) ax[np.isclose(qu[...,0],1.,rtol=0.0)] = [0.0, 0.0, 1.0, 0.0] return ax @@ -564,7 +564,7 @@ class Rotation: with np.errstate(invalid='ignore',divide='ignore'): s = np.linalg.norm(qu[...,1:4],axis=-1,keepdims=True) ro = np.where(np.broadcast_to(np.abs(qu[...,0:1]) < 1.0e-12,qu.shape), - np.block([qu[...,1:2], qu[...,2:3], qu[...,3:4], np.broadcast_to(np.inf,qu.shape[:-1]+(1,))]), + np.block([qu[...,1:2], qu[...,2:3], qu[...,3:4], np.broadcast_to(np.inf,qu[...,0:1].shape)]), np.block([qu[...,1:2]/s,qu[...,2:3]/s,qu[...,3:4]/s, np.tan(np.arccos(np.clip(qu[...,0:1],-1.0,1.0))) ]) From d93ed2bc5c9dfa4a48c4539c74d4dee48087c6db Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 20 Jun 2020 12:20:43 +0200 Subject: [PATCH 08/13] several improvements: - vectorized from_directions - more tests (96% coverage, only random functionality is untested) - updated documentation, folloing numpy standard - inverse operator '~' introduced --- python/damask/_rotation.py | 293 ++++++++++++++++++++++++++-------- python/tests/test_Rotation.py | 41 ++++- 2 files changed, 263 insertions(+), 71 deletions(-) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index 7e9a03c3a..a33b4bc51 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -20,8 +20,8 @@ class Rotation: when viewing from the end point of the rotation axis towards the origin. - rotations will be interpreted in the passive sense. - Euler angle triplets are implemented using the Bunge convention, - with the angular ranges as [0, 2π],[0, π],[0, 2π]. - - the rotation angle ω is limited to the interval [0, π]. + with the angular ranges as [0,2π], [0,π], [0,2π]. + - the rotation angle ω is limited to the interval [0,π]. - the real part of a quaternion is positive, Re(q) > 0 - P = -1 (as default). @@ -49,7 +49,8 @@ class Rotation: Parameters ---------- quaternion : numpy.ndarray, optional - Unit quaternion that follows the conventions. Use .from_quaternion to perform a sanity check. + Unit quaternion in positive real hemisphere. + Use .from_quaternion to perform a sanity check. """ self.quaternion = quaternion.copy() @@ -73,7 +74,7 @@ class Rotation: raise NotImplementedError('Support for multiple rotations missing') return '\n'.join([ 'Quaternion: (real={:.3f}, imag=<{:+.3f}, {:+.3f}, {:+.3f}>)'.format(*(self.quaternion)), - 'Matrix:\n{}'.format(self.as_matrix()), + 'Matrix:\n{}'.format(np.round(self.as_matrix(),8)), 'Bunge Eulers / deg: ({:3.2f}, {:3.2f}, {:3.2f})'.format(*self.as_Eulers(degrees=True)), ]) @@ -87,10 +88,6 @@ 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 - """ if isinstance(other, Rotation): q_m = self.quaternion[...,0:1] @@ -99,7 +96,7 @@ class Rotation: p_o = other.quaternion[...,1:] q = (q_m*q_o - np.einsum('...i,...i',p_m,p_o).reshape(self.shape+(1,))) p = q_m*p_o + q_o*p_m + _P * np.cross(p_m,p_o) - return self.__class__(np.block([q,p])).standardize() + return self.__class__(np.block([q,p]))._standardize() elif isinstance(other,np.ndarray): if self.shape + (3,) == other.shape: @@ -124,24 +121,24 @@ class Rotation: else: raise TypeError('Cannot rotate {}'.format(type(other))) - def inverse(self): - """In-place inverse rotation/backward rotation.""" - self.quaternion[...,1:] *= -1 - return self - def inversed(self): - """Inverse rotation/backward rotation.""" - return self.copy().inverse() - - - def standardize(self): - """In-place quaternion representation with positive real part.""" + def _standardize(self): + """Standardize (ensure positive real hemisphere).""" self.quaternion[self.quaternion[...,0] < 0.0] *= -1 return self - def standardized(self): - """Quaternion representation with positive real part.""" - return self.copy().standardize() + def inverse(self): + """In-place inverse rotation (backward rotation).""" + self.quaternion[...,1:] *= -1 + return self + + def __invert__(self): + """Inverse rotation (backward rotation).""" + return self.copy().inverse() + + def inversed(self): + """Inverse rotation (backward rotation).""" + return ~ self def misorientation(self,other): @@ -154,7 +151,7 @@ class Rotation: Rotation to which the misorientation is computed. """ - return other*self.inversed() + return other@~self def broadcast_to(self,shape): @@ -169,7 +166,7 @@ class Rotation: return self.__class__(q) - def average(self,other): + def average(self,other): #ToDo: discuss calling for vectors """ Calculate the average rotation. @@ -189,25 +186,31 @@ class Rotation: def as_quaternion(self): """ - Unit quaternion [q, p_1, p_2, p_3]. + Represent as unit quaternion. - Parameters - ---------- - quaternion : bool, optional - return quaternion as DAMASK object. + Returns + ------- + q : numpy.ndarray of shape (...,4) + Unit quaternion in positive real hemisphere: (q_0, q_1, q_2, q_3), |q|=1, q_0 ≥ 0. """ - return self.quaternion + return self.quaternion.copy() def as_Eulers(self, degrees = False): """ - Bunge-Euler angles: (φ_1, ϕ, φ_2). + Represent as Bunge-Euler angles. Parameters ---------- degrees : bool, optional - return angles in degrees. + Return angles in degrees. + + Returns + ------- + phi : numpy.ndarray of shape (...,3) + Bunge-Euler angles: (φ_1, ϕ, φ_2), φ_1 ∈ [0,2π], ϕ ∈ [0,π], φ_2 ∈ [0,2π] + unless degrees == True: φ_1 ∈ [0,360], ϕ ∈ [0,180], φ_2 ∈ [0,360] """ eu = Rotation._qu2eu(self.quaternion) @@ -218,14 +221,21 @@ class Rotation: degrees = False, pair = False): """ - Axis angle representation [n_1, n_2, n_3, ω] unless pair == True: ([n_1, n_2, n_3], ω). + Represent as axis angle pair. Parameters ---------- degrees : bool, optional - return rotation angle in degrees. + Return rotation angle in degrees. Defaults to False. pair : bool, optional - return tuple of axis and angle. + Return tuple of axis and angle. Defaults to False. + + Returns + ------- + axis_angle : numpy.ndarray of shape (...,4) unless pair == True: + tuple containing numpy.ndarray of shapes (...,3) and (...) + Axis angle pair: (n_1, n_2, n_3, ω), |n| = 1 and ω ∈ [0,π] + unless degrees = True: ω ∈ [0,180]. """ ax = Rotation._qu2ax(self.quaternion) @@ -233,29 +243,60 @@ class Rotation: return (ax[...,:3],ax[...,3]) if pair else ax def as_matrix(self): - """Rotation matrix.""" + """ + Represent as rotation matrix. + + Returns + ------- + R : numpy.ndarray of shape (...,3,3) + Rotation matrix R, det(R) = 1, R.T∙R=I. + + """ return Rotation._qu2om(self.quaternion) def as_Rodrigues(self, vector = False): """ - Rodrigues-Frank vector representation [n_1, n_2, n_3, tan(ω/2)] unless vector == True: [n_1, n_2, n_3] * tan(ω/2). + Represent as Rodrigues-Frank vector with separated axis and angle argument. Parameters ---------- vector : bool, optional - return as actual Rodrigues--Frank vector, i.e. rotation axis scaled by tan(ω/2). + Return as actual Rodrigues-Frank vector, i.e. axis + and angle argument are not separated. + + Returns + ------- + rho : numpy.ndarray of shape (...,4) unless vector == True: + numpy.ndarray of shape (...,3) + Rodrigues-Frank vector: [n_1, n_2, n_3, tan(ω/2)], |n| = 1 and ω ∈ [0,π]. """ ro = Rotation._qu2ro(self.quaternion) return ro[...,:3]*ro[...,3] if vector else ro def as_homochoric(self): - """Homochoric vector: (h_1, h_2, h_3).""" + """ + Represent as homochoric vector. + + Returns + ------- + h : numpy.ndarray of shape (...,3) + Homochoric vector: (h_1, h_2, h_3), |h| < 1/2*π^(2/3). + + """ return Rotation._qu2ho(self.quaternion) def as_cubochoric(self): - """Cubochoric vector: (c_1, c_2, c_3).""" + """ + Represent as cubocoric vector. + + Returns + ------- + c : numpy.ndarray of shape (...,3) + Cubochoric vector: (c_1, c_2, c_3), max(c_i) < 1/2*π^(2/3). + + """ return Rotation._qu2cu(self.quaternion) def M(self): # ToDo not sure about the name: as_M or M? we do not have a from_M @@ -275,18 +316,34 @@ class Rotation: # Static constructors. The input data needs to follow the conventions, options allow to # relax the conventions. @staticmethod - def from_quaternion(quaternion, + def from_quaternion(q, accept_homomorph = False, P = -1, - acceptHomomorph = None): + acceptHomomorph = None): # old name (for compatibility) + """ + Initialize from quaternion. + Parameters + ---------- + q: numpy.ndarray of shape (...,4) + Unit quaternion in positive real hemisphere: (q_0, q_1, q_2, q_3), + |q|=1, q_0 ≥ 0. + accept_homomorph: boolean, optional + Allow homomorphic variants, i.e. q_0 < 0 (negative real hemisphere). + Defaults to False. + P: integer ∈ {-1,1}, optional + Convention used. Defaults to -1. + + """ if acceptHomomorph is not None: - accept_homomorph = acceptHomomorph - qu = np.array(quaternion,dtype=float) + accept_homomorph = acceptHomomorph # for compatibility + qu = np.array(q,dtype=float) if qu.shape[:-2:-1] != (4,): raise ValueError('Invalid shape.') + if abs(P) != 1: + raise ValueError('P ∉ {-1,1}') - if P > 0: qu[...,1:4] *= -1 # convert from P=1 to P=-1 + if P == 1: qu[...,1:4] *= -1 if accept_homomorph: qu[qu[...,0] < 0.0] *= -1 else: @@ -298,10 +355,21 @@ class Rotation: return Rotation(qu) @staticmethod - def from_Eulers(eulers, + def from_Eulers(phi, degrees = False): + """ + Initialize from Bunge-Euler angles. - eu = np.array(eulers,dtype=float) + Parameters + ---------- + phi: numpy.ndarray of shape (...,3) + Bunge-Euler angles: (φ_1, ϕ, φ_2), φ_1 ∈ [0,2π], ϕ ∈ [0,π], φ_2 ∈ [0,2π] + unless degrees == True: φ_1 ∈ [0,360], ϕ ∈ [0,180], φ_2 ∈ [0,360]. + degrees: boolean, optional + Bunge-Euler angles are given in degrees. Defaults to False. + + """ + eu = np.array(phi,dtype=float) if eu.shape[:-2:-1] != (3,): raise ValueError('Invalid shape.') @@ -316,12 +384,29 @@ class Rotation: degrees = False, normalise = False, P = -1): + """ + Initialize from Axis angle pair. + Parameters + ---------- + axis_angle: numpy.ndarray of shape (...,4) + Axis angle pair: [n_1, n_2, n_3, ω], |n| = 1 and ω ∈ [0,π] + unless degrees = True: ω ∈ [0,180]. + degrees: boolean, optional + Angle ω is given in degrees. Defaults to False. + normalize: boolean, optional + Allow |n| ≠ 1. Defaults to False. + P: integer ∈ {-1,1}, optional + Convention used. Defaults to -1. + + """ ax = np.array(axis_angle,dtype=float) if ax.shape[:-2:-1] != (4,): raise ValueError('Invalid shape.') + if abs(P) != 1: + raise ValueError('P ∉ {-1,1}') - if P > 0: ax[...,0:3] *= -1 # convert from P=1 to P=-1 + if P == 1: ax[...,0:3] *= -1 if degrees: ax[..., 3] = np.radians(ax[...,3]) if normalise: ax[...,0:3] /= np.linalg.norm(ax[...,0:3],axis=-1) if np.any(ax[...,3] < 0.0) or np.any(ax[...,3] > np.pi): @@ -335,7 +420,19 @@ class Rotation: def from_basis(basis, orthonormal = True, reciprocal = False): + """ + Initialize from tbd. + Parameters + ---------- + basis: numpy.ndarray of shape (...,3,3) + tbd + orthonormal: boolean, optional + tbd. Defaults to True. + reciprocal: boolean, optional + tbd. Defaults to False. + + """ om = np.array(basis,dtype=float) if om.shape[:-3:-1] != (3,3): raise ValueError('Invalid shape.') @@ -356,20 +453,64 @@ class Rotation: return Rotation(Rotation._om2qu(om)) @staticmethod - def from_matrix(om): + def from_directions(hkl,uvw): + """ + Initialize from pair of directions/planes. - return Rotation.from_basis(om) + Parameters + ---------- + hkl: numpy.ndarray of shape (...,3) + Direction parallel to z direction, i.e. (h k l) || (0,0,1). + uvw: numpy.ndarray of shape (...,3) + Direction parallel to x direction, i.e. || (1,0,0). + + """ + hkl_ = hkl/np.linalg.norm(hkl,axis=-1,keepdims=True) + uvw_ = uvw/np.linalg.norm(uvw,axis=-1,keepdims=True) + v_1 = np.block([uvw_,np.cross(hkl_,uvw_),hkl_]).reshape(hkl_.shape+(3,)) + v_2 = np.block([uvw_,np.cross(uvw_,hkl_),hkl_]).reshape(hkl_.shape+(3,)) + R = np.where(np.broadcast_to(np.expand_dims(np.expand_dims(np.linalg.det(v_1)>0,-1),-1),v_1.shape), + v_1,v_2) + return Rotation.from_basis(np.swapaxes(R,axis2=-2,axis1=-1)) @staticmethod - def from_Rodrigues(rodrigues, + def from_matrix(R): + """ + Initialize from rotation matrix. + + Parameters + ---------- + R: numpy.ndarray of shape (...,3,3) + Rotation matrix: det(R) = 1, R.T∙R=I. + + """ + return Rotation.from_basis(R) + + @staticmethod + def from_Rodrigues(rho, normalise = False, P = -1): + """ + Initialize from Rodrigues-Frank vector. - ro = np.array(rodrigues,dtype=float) + Parameters + ---------- + rho: numpy.ndarray of shape (...,4) + Rodrigues-Frank vector (angle separated from axis). + (n_1, n_2, n_3, tan(ω/2)), |n| = 1 and ω ∈ [0,π]. + normalize: boolean, optional + Allow |n| ≠ 1. Defaults to False. + P: integer ∈ {-1,1}, optional + Convention used. Defaults to -1. + + """ + ro = np.array(rho,dtype=float) if ro.shape[:-2:-1] != (4,): raise ValueError('Invalid shape.') + if abs(P) != 1: + raise ValueError('P ∉ {-1,1}') - if P > 0: ro[...,0:3] *= -1 # convert from P=1 to P=-1 + if P == 1: ro[...,0:3] *= -1 if normalise: ro[...,0:3] /= np.linalg.norm(ro[...,0:3],axis=-1) if np.any(ro[...,3] < 0.0): raise ValueError('Rodrigues vector rotation angle not positive.') @@ -379,14 +520,26 @@ class Rotation: return Rotation(Rotation._ro2qu(ro)) @staticmethod - def from_homochoric(homochoric, + def from_homochoric(h, P = -1): + """ + Initialize from homochoric vector. - ho = np.array(homochoric,dtype=float) + Parameters + ---------- + h: numpy.ndarray of shape (...,3) + Homochoric vector: (h_1, h_2, h_3), |h| < (3/4*π)^(1/3). + P: integer ∈ {-1,1}, optional + Convention used. Defaults to -1. + + """ + ho = np.array(h,dtype=float) if ho.shape[:-2:-1] != (3,): raise ValueError('Invalid shape.') + if abs(P) != 1: + raise ValueError('P ∉ {-1,1}') - if P > 0: ho *= -1 # convert from P=1 to P=-1 + if P == 1: ho *= -1 if np.any(np.linalg.norm(ho,axis=-1) >_R1+1e-9): raise ValueError('Homochoric coordinate outside of the sphere.') @@ -394,18 +547,30 @@ class Rotation: return Rotation(Rotation._ho2qu(ho)) @staticmethod - def from_cubochoric(cubochoric, - P = -1): + def from_cubochoric(c, + P = -1): + """ + Initialize from cubochoric vector. - cu = np.array(cubochoric,dtype=float) + Parameters + ---------- + c: numpy.ndarray of shape (...,3) + Cubochoric vector: (c_1, c_2, c_3), max(c_i) < 1/2*π^(2/3). + P: integer ∈ {-1,1}, optional + Convention used. Defaults to -1. + + """ + cu = np.array(c,dtype=float) if cu.shape[:-2:-1] != (3,): raise ValueError('Invalid shape.') + if abs(P) != 1: + raise ValueError('P ∉ {-1,1}') - if np.abs(np.max(cu))>np.pi**(2./3.) * 0.5+1e-9: - raise ValueError('Cubochoric coordinate outside of the cube: {} {} {}.'.format(*cu)) + if np.abs(np.max(cu)) > np.pi**(2./3.) * 0.5+1e-9: + raise ValueError('Cubochoric coordinate outside of the cube.') ho = Rotation._cu2ho(cu) - if P > 0: ho *= -1 # convert from P=1 to P=-1 + if P == 1: ho *= -1 return Rotation(Rotation._ho2qu(ho)) @@ -458,7 +623,7 @@ class Rotation: np.cos(2.0*np.pi*r[...,1])*B, np.sin(2.0*np.pi*r[...,0])*A],axis=-1) - return Rotation(q.reshape(r.shape[:-1]+(4,)) if shape is not None else q).standardize() + return Rotation(q.reshape(r.shape[:-1]+(4,)) if shape is not None else q)._standardize() # for compatibility (old names do not follow convention) diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py index 5a1cd145d..a6b4fea69 100644 --- a/python/tests/test_Rotation.py +++ b/python/tests/test_Rotation.py @@ -14,7 +14,7 @@ scatter=1.e-2 @pytest.fixture def default(): - """A set of n random rotations.""" + """A set of n rotations (corner cases and random).""" specials = np.array([ [1.0, 0.0, 0.0, 0.0], #---------------------- @@ -170,7 +170,7 @@ def qu2ax(qu): Modified version of the original formulation, should be numerically more stable """ - if np.isclose(qu[0],1.,rtol=0.0): # set axis to [001] if the angle is 0/360 + if np.isclose(qu[0],1.,rtol=0.0): # set axis to [001] if the angle is 0/360 ax = np.array([ 0.0, 0.0, 1.0, 0.0 ]) elif qu[0] > 1.e-8: s = np.sign(qu[0])/np.sqrt(qu[1]**2+qu[2]**2+qu[3]**2) @@ -534,7 +534,7 @@ def mul(me, other): other_p = other.quaternion[1:] R = me.__class__(np.append(me_q*other_q - np.dot(me_p,other_p), me_q*other_p + other_q*me_p + _P * np.cross(me_p,other_p))) - return R.standardize() + return R._standardize() elif isinstance(other, np.ndarray): if other.shape == (3,): A = me.quaternion[0]**2.0 - np.dot(me.quaternion[1:],me.quaternion[1:]) @@ -554,7 +554,7 @@ def mul(me, other): axes = ((0, 2, 4, 6), (0, 1, 2, 3)) return np.tensordot(RRRR, other, axes) else: - raise ValueError('Can only rotate vectors, 2nd order ternsors, and 4th order tensors') + raise ValueError('Can only rotate vectors, 2nd order tensors, and 4th order tensors') else: raise TypeError('Cannot rotate {}'.format(type(other))) @@ -854,6 +854,15 @@ class TestRotation: rot = Rotation.from_basis(om,False,reciprocal=reciprocal) assert np.isclose(np.linalg.det(rot.as_matrix()),1.0) + def test_directions(self): + hkl = np.array([0.,0.,1.]) + uvw = np.array([1.,0.,0.]) + assert np.allclose(Rotation.from_directions(hkl,uvw).as_matrix(),np.eye(3)) + + @pytest.mark.parametrize('shape',[None,1,(4,4)]) + def test_random(self,shape): + Rotation.from_random(shape) + @pytest.mark.parametrize('function',[Rotation.from_quaternion, Rotation.from_Eulers, Rotation.from_axis_angle, @@ -866,6 +875,16 @@ class TestRotation: with pytest.raises(ValueError): function(invalid_shape) + @pytest.mark.parametrize('fr,to',[(Rotation.from_quaternion,'as_quaternion'), + (Rotation.from_axis_angle,'as_axis_angle'), + (Rotation.from_Rodrigues, 'as_Rodrigues'), + (Rotation.from_homochoric,'as_homochoric'), + (Rotation.from_cubochoric,'as_cubochoric')]) + def test_invalid_P(self,fr,to): + R = Rotation.from_random(np.random.randint(8,32,(3))) # noqa + with pytest.raises(ValueError): + fr(eval('R.{}()'.format(to)),P=-30) + @pytest.mark.parametrize('shape',[None,(3,),(4,2)]) def test_broadcast(self,shape): rot = Rotation.from_random(shape) @@ -932,14 +951,18 @@ class TestRotation: phi_2 = 2*np.pi - phi_1 R_1 = Rotation.from_Eulers(np.array([phi_1,0.,0.])) R_2 = Rotation.from_Eulers(np.array([0.,0.,phi_2])) - assert np.allclose(data,R_2*(R_1*data)) + assert np.allclose(data,R_2@(R_1@data)) + + def test_rotate_inverse(self): + R = Rotation.from_random() + assert np.allclose(np.eye(3),(R.inversed()@R).as_matrix()) @pytest.mark.parametrize('data',[np.random.rand(3), np.random.rand(3,3), np.random.rand(3,3,3,3)]) - def test_rotate_inverse(self,data): + def test_rotate_inverse_array(self,data): R = Rotation.from_random() - assert np.allclose(data,R.inversed()*(R*data)) + assert np.allclose(data,R.inversed()@(R@data)) @pytest.mark.parametrize('data',[np.random.rand(4), np.random.rand(3,2), @@ -956,3 +979,7 @@ class TestRotation: R = Rotation.from_random() with pytest.raises(TypeError): R*data + + def test_misorientation(self): + R = Rotation.from_random() + assert np.allclose(R.misorientation(R).as_matrix(),np.eye(3)) From 3f63a4fdbc56f96dae3f87952b91c9c7d8a667db Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 20 Jun 2020 18:13:34 +0200 Subject: [PATCH 09/13] [skip ci] typo --- python/damask/_rotation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index a33b4bc51..09aeab256 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -289,7 +289,7 @@ class Rotation: def as_cubochoric(self): """ - Represent as cubocoric vector. + Represent as cubochoric vector. Returns ------- From 15b43bcebf9ae932d54e5a623b82177351214569 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 20 Jun 2020 19:57:49 +0200 Subject: [PATCH 10/13] from_directions is not general, removed polishing --- python/damask/_rotation.py | 67 ++++++++++++----------------------- python/tests/test_Rotation.py | 5 --- 2 files changed, 23 insertions(+), 49 deletions(-) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index 09aeab256..7c469a2b5 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -325,13 +325,13 @@ class Rotation: Parameters ---------- - q: numpy.ndarray of shape (...,4) + q : numpy.ndarray of shape (...,4) Unit quaternion in positive real hemisphere: (q_0, q_1, q_2, q_3), |q|=1, q_0 ≥ 0. - accept_homomorph: boolean, optional + accept_homomorph : boolean, optional Allow homomorphic variants, i.e. q_0 < 0 (negative real hemisphere). Defaults to False. - P: integer ∈ {-1,1}, optional + P : integer ∈ {-1,1}, optional Convention used. Defaults to -1. """ @@ -362,10 +362,10 @@ class Rotation: Parameters ---------- - phi: numpy.ndarray of shape (...,3) + phi : numpy.ndarray of shape (...,3) Bunge-Euler angles: (φ_1, ϕ, φ_2), φ_1 ∈ [0,2π], ϕ ∈ [0,π], φ_2 ∈ [0,2π] unless degrees == True: φ_1 ∈ [0,360], ϕ ∈ [0,180], φ_2 ∈ [0,360]. - degrees: boolean, optional + degrees : boolean, optional Bunge-Euler angles are given in degrees. Defaults to False. """ @@ -389,14 +389,14 @@ class Rotation: Parameters ---------- - axis_angle: numpy.ndarray of shape (...,4) + axis_angle : numpy.ndarray of shape (...,4) Axis angle pair: [n_1, n_2, n_3, ω], |n| = 1 and ω ∈ [0,π] unless degrees = True: ω ∈ [0,180]. - degrees: boolean, optional + degrees : boolean, optional Angle ω is given in degrees. Defaults to False. normalize: boolean, optional Allow |n| ≠ 1. Defaults to False. - P: integer ∈ {-1,1}, optional + P : integer ∈ {-1,1}, optional Convention used. Defaults to -1. """ @@ -421,16 +421,16 @@ class Rotation: orthonormal = True, reciprocal = False): """ - Initialize from tbd. + Initialize from lattice basis vectors. Parameters ---------- - basis: numpy.ndarray of shape (...,3,3) - tbd - orthonormal: boolean, optional - tbd. Defaults to True. - reciprocal: boolean, optional - tbd. Defaults to False. + basis : numpy.ndarray of shape (...,3,3) + Three lattice basis vectors in three dimensions. + orthonormal : boolean, optional + Basis is strictly orthonormal, i.e. is free of stretch components. Defaults to True. + reciprocal : boolean, optional + Basis vectors are given in reciprocal (instead of real) space. Defaults to False. """ om = np.array(basis,dtype=float) @@ -452,27 +452,6 @@ class Rotation: return Rotation(Rotation._om2qu(om)) - @staticmethod - def from_directions(hkl,uvw): - """ - Initialize from pair of directions/planes. - - Parameters - ---------- - hkl: numpy.ndarray of shape (...,3) - Direction parallel to z direction, i.e. (h k l) || (0,0,1). - uvw: numpy.ndarray of shape (...,3) - Direction parallel to x direction, i.e. || (1,0,0). - - """ - hkl_ = hkl/np.linalg.norm(hkl,axis=-1,keepdims=True) - uvw_ = uvw/np.linalg.norm(uvw,axis=-1,keepdims=True) - v_1 = np.block([uvw_,np.cross(hkl_,uvw_),hkl_]).reshape(hkl_.shape+(3,)) - v_2 = np.block([uvw_,np.cross(uvw_,hkl_),hkl_]).reshape(hkl_.shape+(3,)) - R = np.where(np.broadcast_to(np.expand_dims(np.expand_dims(np.linalg.det(v_1)>0,-1),-1),v_1.shape), - v_1,v_2) - return Rotation.from_basis(np.swapaxes(R,axis2=-2,axis1=-1)) - @staticmethod def from_matrix(R): """ @@ -480,7 +459,7 @@ class Rotation: Parameters ---------- - R: numpy.ndarray of shape (...,3,3) + R : numpy.ndarray of shape (...,3,3) Rotation matrix: det(R) = 1, R.T∙R=I. """ @@ -495,12 +474,12 @@ class Rotation: Parameters ---------- - rho: numpy.ndarray of shape (...,4) + rho : numpy.ndarray of shape (...,4) Rodrigues-Frank vector (angle separated from axis). (n_1, n_2, n_3, tan(ω/2)), |n| = 1 and ω ∈ [0,π]. - normalize: boolean, optional + normalize : boolean, optional Allow |n| ≠ 1. Defaults to False. - P: integer ∈ {-1,1}, optional + P : integer ∈ {-1,1}, optional Convention used. Defaults to -1. """ @@ -527,9 +506,9 @@ class Rotation: Parameters ---------- - h: numpy.ndarray of shape (...,3) + h : numpy.ndarray of shape (...,3) Homochoric vector: (h_1, h_2, h_3), |h| < (3/4*π)^(1/3). - P: integer ∈ {-1,1}, optional + P : integer ∈ {-1,1}, optional Convention used. Defaults to -1. """ @@ -554,9 +533,9 @@ class Rotation: Parameters ---------- - c: numpy.ndarray of shape (...,3) + c : numpy.ndarray of shape (...,3) Cubochoric vector: (c_1, c_2, c_3), max(c_i) < 1/2*π^(2/3). - P: integer ∈ {-1,1}, optional + P : integer ∈ {-1,1}, optional Convention used. Defaults to -1. """ diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py index a6b4fea69..933719d6f 100644 --- a/python/tests/test_Rotation.py +++ b/python/tests/test_Rotation.py @@ -854,11 +854,6 @@ class TestRotation: rot = Rotation.from_basis(om,False,reciprocal=reciprocal) assert np.isclose(np.linalg.det(rot.as_matrix()),1.0) - def test_directions(self): - hkl = np.array([0.,0.,1.]) - uvw = np.array([1.,0.,0.]) - assert np.allclose(Rotation.from_directions(hkl,uvw).as_matrix(),np.eye(3)) - @pytest.mark.parametrize('shape',[None,1,(4,4)]) def test_random(self,shape): Rotation.from_random(shape) From 4c5939ef233ff4021d0f836f836df3d21ae25c91 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 20 Jun 2020 22:51:00 +0200 Subject: [PATCH 11/13] small polishing --- python/damask/util.py | 14 +++++++------- python/tests/test_Rotation.py | 4 ++-- python/tests/test_util.py | 11 +++++++++++ src/discretization.f90 | 2 +- 4 files changed, 21 insertions(+), 10 deletions(-) create mode 100644 python/tests/test_util.py diff --git a/python/damask/util.py b/python/damask/util.py index cb1d8757d..a3fc71d3a 100644 --- a/python/damask/util.py +++ b/python/damask/util.py @@ -93,7 +93,7 @@ def strikeout(what): def execute(cmd, - streamIn = None, + stream_in = None, wd = './', env = None): """ @@ -103,7 +103,7 @@ def execute(cmd, ---------- cmd : str Command to be executed. - streanIn :, optional + stream_in : file object, optional Input (via pipe) for executed process. wd : str, optional Working directory of process. Defaults to ./ . @@ -119,14 +119,14 @@ def execute(cmd, stderr = subprocess.PIPE, stdin = subprocess.PIPE, env = myEnv) - out,error = [i for i in (process.communicate() if streamIn is None - else process.communicate(streamIn.read().encode('utf-8')))] + stdout, stderr = [i for i in (process.communicate() if stream_in is None + else process.communicate(stream_in.read().encode('utf-8')))] os.chdir(initialPath) - out = out.decode('utf-8').replace('\x08','') - error = error.decode('utf-8').replace('\x08','') + stdout = stdout.decode('utf-8').replace('\x08','') + stderr = stderr.decode('utf-8').replace('\x08','') if process.returncode != 0: raise RuntimeError('{} failed with returncode {}'.format(cmd,process.returncode)) - return out,error + return stdout, stderr def show_progress(iterable,N_iter=None,prefix='',bar_length=50): diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py index 933719d6f..86e5fa2a7 100644 --- a/python/tests/test_Rotation.py +++ b/python/tests/test_Rotation.py @@ -950,14 +950,14 @@ class TestRotation: def test_rotate_inverse(self): R = Rotation.from_random() - assert np.allclose(np.eye(3),(R.inversed()@R).as_matrix()) + assert np.allclose(np.eye(3),(~R@R).as_matrix()) @pytest.mark.parametrize('data',[np.random.rand(3), np.random.rand(3,3), np.random.rand(3,3,3,3)]) def test_rotate_inverse_array(self,data): R = Rotation.from_random() - assert np.allclose(data,R.inversed()@(R@data)) + assert np.allclose(data,~R@(R@data)) @pytest.mark.parametrize('data',[np.random.rand(4), np.random.rand(3,2), diff --git a/python/tests/test_util.py b/python/tests/test_util.py new file mode 100644 index 000000000..e67478d82 --- /dev/null +++ b/python/tests/test_util.py @@ -0,0 +1,11 @@ +from damask import util + +class TestUtil: + + def test_execute_direct(self): + out,err = util.execute('echo test') + assert out=='test\n' and err=='' + + def test_execute_env(self): + out,err = util.execute('sh -c "echo $test_for_execute"',env={'test_for_execute':'test'}) + assert out=='test\n' and err=='' diff --git a/src/discretization.f90 b/src/discretization.f90 index 16f76b158..e52784c6d 100644 --- a/src/discretization.f90 +++ b/src/discretization.f90 @@ -49,7 +49,7 @@ subroutine discretization_init(homogenizationAt,microstructureAt,& IPcoords0, & NodeCoords0 integer, optional, intent(in) :: & - sharedNodesBegin ! index of first node shared among different processes (MPI) + sharedNodesBegin !< index of first node shared among different processes (MPI) write(6,'(/,a)') ' <<<+- discretization init -+>>>'; flush(6) From d4efadb333f864b4ebef645373cdbeccb2356e83 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 21 Jun 2020 10:03:52 +0200 Subject: [PATCH 12/13] should be availabe outside of this module --- src/quaternions.f90 | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/quaternions.f90 b/src/quaternions.f90 index f5f5276e8..c337bd681 100644 --- a/src/quaternions.f90 +++ b/src/quaternions.f90 @@ -101,6 +101,8 @@ module quaternions assignment(=), & conjg, aimag, & log, exp, & + abs, dot_product, & + inverse, & real contains From c6a5bb8a3b873a66ba968339d76c754bc63c0d4d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 21 Jun 2020 10:04:45 +0200 Subject: [PATCH 13/13] is 2020 --- python/damask/_rotation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index 7c469a2b5..f433b2b37 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -615,7 +615,7 @@ class Rotation: #################################################################################################### # Code below available according to the following conditions on https://github.com/MarDiehl/3Drotations #################################################################################################### -# Copyright (c) 2017-2019, Martin Diehl/Max-Planck-Institut für Eisenforschung GmbH +# Copyright (c) 2017-2020, Martin Diehl/Max-Planck-Institut für Eisenforschung GmbH # Copyright (c) 2013-2014, Marc De Graef/Carnegie Mellon University # All rights reserved. #