From 9a4e9e62b66a84abe610ee17f9aaa07f081e5ee7 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 1 Feb 2019 20:28:07 +0100 Subject: [PATCH] using new rotation class --- python/damask/orientation.py | 1 - src/constitutive.f90 | 4 +-- src/crystallite.f90 | 48 ++++++++++++------------------------ src/plastic_nonlocal.f90 | 15 ++++++----- 4 files changed, 24 insertions(+), 44 deletions(-) diff --git a/python/damask/orientation.py b/python/damask/orientation.py index 0386bc714..a1fe1f845 100644 --- a/python/damask/orientation.py +++ b/python/damask/orientation.py @@ -575,7 +575,6 @@ class Symmetry: proper considers only vectors with z >= 0, hence uses two neighboring SSTs. Return inverse pole figure color if requested. """ - if self.lattice == 'cubic': basis = {'improper':np.array([ [-1. , 0. , 1. ], [ np.sqrt(2.) , -np.sqrt(2.) , 0. ], diff --git a/src/constitutive.f90 b/src/constitutive.f90 index a0d7147a6..9483f2610 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -348,7 +348,7 @@ end function constitutive_homogenizedC !-------------------------------------------------------------------------------------------------- !> @brief calls microstructure function of the different constitutive models !-------------------------------------------------------------------------------------------------- -subroutine constitutive_microstructure(orientations, Fe, Fp, ipc, ip, el) +subroutine constitutive_microstructure(Fe, Fp, ipc, ip, el) use prec, only: & pReal use material, only: & @@ -381,8 +381,6 @@ subroutine constitutive_microstructure(orientations, Fe, Fp, ipc, ip, el) ho, & !< homogenization tme, & !< thermal member position instance, of - real(pReal), intent(in), dimension(:,:,:,:) :: & - orientations !< crystal orientations as quaternions ho = material_homogenizationAt(el) tme = thermalMapping(ho)%p(ip,el) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index e2f05284e..013c28a38 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -46,12 +46,8 @@ module crystallite crystallite_Tstar0_v, & !< 2nd Piola-Kirchhoff stress vector at start of FE inc ToDo: Should be called S, 3x3 crystallite_partionedTstar0_v !< 2nd Piola-Kirchhoff stress vector at start of homog inc ToDo: Should be called S, 3x3 type(rotation), dimension(:,:,:), allocatable, private :: & - crystallite_ori, & !< orientation as quaternion - crystallite_ori0 !< initial orientation as quaternion - real(pReal), dimension(:,:,:,:), allocatable, private :: & - crystallite_orientation, & !< orientation as quaternion - crystallite_orientation0, & !< initial orientation as quaternion - crystallite_rotation !< grain rotation away from initial orientation as axis-angle (in degrees) in crystal reference frame + crystallite_orientation, & !< orientation + crystallite_orientation0 !< initial orientation real(pReal), dimension(:,:,:,:,:), allocatable, public, protected :: & crystallite_Fe, & !< current "elastic" def grad (end of converged time step) crystallite_P !< 1st Piola-Kirchhoff stress per grain @@ -243,11 +239,8 @@ subroutine crystallite_init allocate(crystallite_subdt(cMax,iMax,eMax), source=0.0_pReal) allocate(crystallite_subFrac(cMax,iMax,eMax), source=0.0_pReal) allocate(crystallite_subStep(cMax,iMax,eMax), source=0.0_pReal) - allocate(crystallite_orientation(4,cMax,iMax,eMax), source=0.0_pReal) - allocate(crystallite_orientation0(4,cMax,iMax,eMax), source=0.0_pReal) - allocate(crystallite_ori(cMax,iMax,eMax)) - !allocate(crystallite_ori0(cMax,iMax,eMax)) - allocate(crystallite_rotation(4,cMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_orientation(cMax,iMax,eMax)) + allocate(crystallite_orientation0(cMax,iMax,eMax)) allocate(crystallite_localPlasticity(cMax,iMax,eMax), source=.true.) allocate(crystallite_requested(cMax,iMax,eMax), source=.false.) allocate(crystallite_todo(cMax,iMax,eMax), source=.false.) @@ -296,8 +289,6 @@ subroutine crystallite_init crystallite_outputID(o,c) = orientation_ID case ('grainrotation') outputName crystallite_outputID(o,c) = grainrotation_ID - case ('eulerangles') outputName - crystallite_outputID(o,c) = eulerangles_ID case ('defgrad','f') outputName crystallite_outputID(o,c) = defgrad_ID case ('fe') outputName @@ -334,8 +325,6 @@ subroutine crystallite_init mySize = 1_pInt case(orientation_ID,grainrotation_ID) mySize = 4_pInt - case(eulerangles_ID) - mySize = 3_pInt case(defgrad_ID,fe_ID,fp_ID,fi_ID,lp_ID,li_ID,p_ID,s_ID) mySize = 9_pInt case(elasmatrix_ID) @@ -401,13 +390,12 @@ subroutine crystallite_init call crystallite_orientations() crystallite_orientation0 = crystallite_orientation ! store initial orientations for calculation of grain rotations - + !$OMP PARALLEL DO do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) do c = 1_pInt,homogenization_Ngrains(mesh_element(3,e)) - call constitutive_microstructure(crystallite_orientation, & - crystallite_Fe(1:3,1:3,c,i,e), & + call constitutive_microstructure(crystallite_Fe(1:3,1:3,c,i,e), & crystallite_Fp(1:3,1:3,c,i,e), & c,i,e) ! update dependent state variables to be consistent with basic states enddo @@ -908,10 +896,7 @@ subroutine crystallite_orientations do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) do c = 1_pInt,homogenization_Ngrains(mesh_element(3,e)) - call crystallite_ori(c,i,e)%fromRotationMatrix(transpose(math_rotationalPart33(crystallite_Fe(1:3,1:3,c,i,e)))) - crystallite_orientation(1:4,c,i,e) = math_RtoQ(transpose(math_rotationalPart33(crystallite_Fe(1:3,1:3,c,i,e)))) - crystallite_rotation(1:4,c,i,e) = lattice_qDisorientation(crystallite_orientation0(1:4,c,i,e), &! active rotation from initial - crystallite_orientation(1:4,c,i,e)) ! to current orientation (with no symmetry) + call crystallite_orientation(c,i,e)%fromRotationMatrix(transpose(math_rotationalPart33(crystallite_Fe(1:3,1:3,c,i,e)))) enddo; enddo; enddo !$OMP END PARALLEL DO @@ -920,7 +905,7 @@ subroutine crystallite_orientations !$OMP PARALLEL DO do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) - if (plasticState(material_phase(1,i,e))%nonLocal) & ! if nonlocal model + if (plasticState(material_phase(1,i,e))%nonLocal) & ! if nonlocal model call plastic_nonlocal_updateCompatibility(crystallite_orientation,i,e) enddo; enddo !$OMP END PARALLEL DO @@ -987,6 +972,8 @@ function crystallite_postResults(ipc, ip, el) use constitutive, only: & constitutive_homogenizedC, & constitutive_postResults + use rotations, only: & + rotation implicit none integer(pInt), intent(in):: & @@ -1006,6 +993,7 @@ function crystallite_postResults(ipc, ip, el) crystID, & mySize, & n + type(rotation) :: rot crystID = microstructure_crystallite(mesh_element(4,el)) @@ -1029,15 +1017,12 @@ function crystallite_postResults(ipc, ip, el) / real(homogenization_Ngrains(mesh_element(3,el)),pReal) ! grain volume (not fraction but absolute) case (orientation_ID) mySize = 4_pInt - crystallite_postResults(c+1:c+mySize) = crystallite_ori(ipc,ip,el)%asQuaternion() - case (eulerangles_ID) - mySize = 3_pInt - crystallite_postResults(c+1:c+mySize) = inDeg & - * math_qToEuler(crystallite_orientation(1:4,ipc,ip,el)) ! grain orientation as Euler angles in degree + crystallite_postResults(c+1:c+mySize) = crystallite_orientation(ipc,ip,el)%asQuaternion() + case (grainrotation_ID) + rot = crystallite_orientation0(ipc,ip,el)%misorientation(crystallite_orientation(ipc,ip,el)) mySize = 4_pInt - crystallite_postResults(c+1:c+mySize) = & - math_qToEulerAxisAngle(crystallite_rotation(1:4,ipc,ip,el)) ! grain rotation away from initial orientation as axis-angle in sample reference coordinates + crystallite_postResults(c+1:c+mySize) = rot%asAxisAnglePair() crystallite_postResults(c+4) = inDeg * crystallite_postResults(c+4) ! angle in degree ! remark: tensor output is of the form 11,12,13, 21,22,23, 31,32,33 @@ -2222,8 +2207,7 @@ subroutine update_dependentState() do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) do g = 1,homogenization_Ngrains(mesh_element(3,e)) if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) & - call constitutive_dependentState(crystallite_orientation, & - crystallite_Fe(1:3,1:3,g,i,e), & + call constitutive_dependentState(crystallite_Fe(1:3,1:3,g,i,e), & crystallite_Fp(1:3,1:3,g,i,e), & g, i, e) enddo; enddo; enddo diff --git a/src/plastic_nonlocal.f90 b/src/plastic_nonlocal.f90 index e1355da8f..112592a8c 100644 --- a/src/plastic_nonlocal.f90 +++ b/src/plastic_nonlocal.f90 @@ -3005,10 +3005,9 @@ end subroutine plastic_nonlocal_dotState !* that sum up to a total of 1 are considered, all others are set to * !* zero. * !********************************************************************* -subroutine plastic_nonlocal_updateCompatibility(orientation,i,e) - -use math, only: math_mul3x3, & - math_qRot +subroutine plastic_nonlocal_updateCompatibility(orientation,i,e) +use math, only: math_mul3x3, math_qRot +use rotations, only: rotation use material, only: material_phase, & material_texture, & phase_localPlasticity, & @@ -3030,7 +3029,7 @@ implicit none !* input variables integer(pInt), intent(in) :: i, & ! ip index e ! element index -real(pReal), dimension(4,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & +type(rotation), dimension(1,mesh_maxNips,mesh_NcpElems), intent(in) :: & orientation ! crystal orientation in quaternions !* local variables @@ -3059,7 +3058,7 @@ real(pReal) my_compatibilitySum, & nThresholdValues logical, dimension(totalNslip(phase_plasticityInstance(material_phase(1,i,e)))) :: & belowThreshold - +type(rotation) :: rot Nneighbors = FE_NipNeighbors(FE_celltype(FE_geomtype(mesh_element(2,e)))) ph = material_phase(1,i,e) @@ -3129,8 +3128,8 @@ neighbors: do n = 1_pInt,Nneighbors !* Finally the smallest my_compatibility value is decreased until the sum is exactly equal to one. !* All values below the threshold are set to zero. else - absoluteMisorientation = lattice_qDisorientation(orientation(1:4,1,i,e), & - orientation(1:4,1,neighbor_i,neighbor_e)) ! no symmetry + rot = orientation(1,i,e)%misorientation(orientation(1,neighbor_i,neighbor_e)) + absoluteMisorientation = rot%asQuaternion() mySlipSystems: do s1 = 1_pInt,ns neighborSlipSystems: do s2 = 1_pInt,ns my_compatibility(1,s2,s1,n) = math_mul3x3(slipNormal(1:3,s1), math_qRot(absoluteMisorientation, slipNormal(1:3,s2))) &