using new rotation class

This commit is contained in:
Martin Diehl 2019-02-01 20:28:07 +01:00
parent 08009079ff
commit 9a4e9e62b6
4 changed files with 24 additions and 44 deletions

View File

@ -575,7 +575,6 @@ class Symmetry:
proper considers only vectors with z >= 0, hence uses two neighboring SSTs. proper considers only vectors with z >= 0, hence uses two neighboring SSTs.
Return inverse pole figure color if requested. Return inverse pole figure color if requested.
""" """
if self.lattice == 'cubic': if self.lattice == 'cubic':
basis = {'improper':np.array([ [-1. , 0. , 1. ], basis = {'improper':np.array([ [-1. , 0. , 1. ],
[ np.sqrt(2.) , -np.sqrt(2.) , 0. ], [ np.sqrt(2.) , -np.sqrt(2.) , 0. ],

View File

@ -348,7 +348,7 @@ end function constitutive_homogenizedC
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calls microstructure function of the different constitutive models !> @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: & use prec, only: &
pReal pReal
use material, only: & use material, only: &
@ -381,8 +381,6 @@ subroutine constitutive_microstructure(orientations, Fe, Fp, ipc, ip, el)
ho, & !< homogenization ho, & !< homogenization
tme, & !< thermal member position tme, & !< thermal member position
instance, of instance, of
real(pReal), intent(in), dimension(:,:,:,:) :: &
orientations !< crystal orientations as quaternions
ho = material_homogenizationAt(el) ho = material_homogenizationAt(el)
tme = thermalMapping(ho)%p(ip,el) tme = thermalMapping(ho)%p(ip,el)

View File

@ -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_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 crystallite_partionedTstar0_v !< 2nd Piola-Kirchhoff stress vector at start of homog inc ToDo: Should be called S, 3x3
type(rotation), dimension(:,:,:), allocatable, private :: & type(rotation), dimension(:,:,:), allocatable, private :: &
crystallite_ori, & !< orientation as quaternion crystallite_orientation, & !< orientation
crystallite_ori0 !< initial orientation as quaternion crystallite_orientation0 !< initial orientation
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
real(pReal), dimension(:,:,:,:,:), allocatable, public, protected :: & real(pReal), dimension(:,:,:,:,:), allocatable, public, protected :: &
crystallite_Fe, & !< current "elastic" def grad (end of converged time step) crystallite_Fe, & !< current "elastic" def grad (end of converged time step)
crystallite_P !< 1st Piola-Kirchhoff stress per grain 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_subdt(cMax,iMax,eMax), source=0.0_pReal)
allocate(crystallite_subFrac(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_subStep(cMax,iMax,eMax), source=0.0_pReal)
allocate(crystallite_orientation(4,cMax,iMax,eMax), source=0.0_pReal) allocate(crystallite_orientation(cMax,iMax,eMax))
allocate(crystallite_orientation0(4,cMax,iMax,eMax), source=0.0_pReal) allocate(crystallite_orientation0(cMax,iMax,eMax))
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_localPlasticity(cMax,iMax,eMax), source=.true.) allocate(crystallite_localPlasticity(cMax,iMax,eMax), source=.true.)
allocate(crystallite_requested(cMax,iMax,eMax), source=.false.) allocate(crystallite_requested(cMax,iMax,eMax), source=.false.)
allocate(crystallite_todo(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 crystallite_outputID(o,c) = orientation_ID
case ('grainrotation') outputName case ('grainrotation') outputName
crystallite_outputID(o,c) = grainrotation_ID crystallite_outputID(o,c) = grainrotation_ID
case ('eulerangles') outputName
crystallite_outputID(o,c) = eulerangles_ID
case ('defgrad','f') outputName case ('defgrad','f') outputName
crystallite_outputID(o,c) = defgrad_ID crystallite_outputID(o,c) = defgrad_ID
case ('fe') outputName case ('fe') outputName
@ -334,8 +325,6 @@ subroutine crystallite_init
mySize = 1_pInt mySize = 1_pInt
case(orientation_ID,grainrotation_ID) case(orientation_ID,grainrotation_ID)
mySize = 4_pInt 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) case(defgrad_ID,fe_ID,fp_ID,fi_ID,lp_ID,li_ID,p_ID,s_ID)
mySize = 9_pInt mySize = 9_pInt
case(elasmatrix_ID) case(elasmatrix_ID)
@ -401,13 +390,12 @@ subroutine crystallite_init
call crystallite_orientations() call crystallite_orientations()
crystallite_orientation0 = crystallite_orientation ! store initial orientations for calculation of grain rotations crystallite_orientation0 = crystallite_orientation ! store initial orientations for calculation of grain rotations
!$OMP PARALLEL DO !$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2) do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
do c = 1_pInt,homogenization_Ngrains(mesh_element(3,e)) do c = 1_pInt,homogenization_Ngrains(mesh_element(3,e))
call constitutive_microstructure(crystallite_orientation, & call constitutive_microstructure(crystallite_Fe(1:3,1:3,c,i,e), &
crystallite_Fe(1:3,1:3,c,i,e), &
crystallite_Fp(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 c,i,e) ! update dependent state variables to be consistent with basic states
enddo enddo
@ -908,10 +896,7 @@ subroutine crystallite_orientations
do e = FEsolving_execElem(1),FEsolving_execElem(2) do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
do c = 1_pInt,homogenization_Ngrains(mesh_element(3,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)))) call crystallite_orientation(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)
enddo; enddo; enddo enddo; enddo; enddo
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
@ -920,7 +905,7 @@ subroutine crystallite_orientations
!$OMP PARALLEL DO !$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2) do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) 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) call plastic_nonlocal_updateCompatibility(crystallite_orientation,i,e)
enddo; enddo enddo; enddo
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
@ -987,6 +972,8 @@ function crystallite_postResults(ipc, ip, el)
use constitutive, only: & use constitutive, only: &
constitutive_homogenizedC, & constitutive_homogenizedC, &
constitutive_postResults constitutive_postResults
use rotations, only: &
rotation
implicit none implicit none
integer(pInt), intent(in):: & integer(pInt), intent(in):: &
@ -1006,6 +993,7 @@ function crystallite_postResults(ipc, ip, el)
crystID, & crystID, &
mySize, & mySize, &
n n
type(rotation) :: rot
crystID = microstructure_crystallite(mesh_element(4,el)) 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) / real(homogenization_Ngrains(mesh_element(3,el)),pReal) ! grain volume (not fraction but absolute)
case (orientation_ID) case (orientation_ID)
mySize = 4_pInt mySize = 4_pInt
crystallite_postResults(c+1:c+mySize) = crystallite_ori(ipc,ip,el)%asQuaternion() crystallite_postResults(c+1:c+mySize) = crystallite_orientation(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
case (grainrotation_ID) case (grainrotation_ID)
rot = crystallite_orientation0(ipc,ip,el)%misorientation(crystallite_orientation(ipc,ip,el))
mySize = 4_pInt mySize = 4_pInt
crystallite_postResults(c+1:c+mySize) = & crystallite_postResults(c+1:c+mySize) = rot%asAxisAnglePair()
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+4) = inDeg * crystallite_postResults(c+4) ! angle in degree 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 ! 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 i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
do g = 1,homogenization_Ngrains(mesh_element(3,e)) do g = 1,homogenization_Ngrains(mesh_element(3,e))
if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) & if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) &
call constitutive_dependentState(crystallite_orientation, & call constitutive_dependentState(crystallite_Fe(1:3,1:3,g,i,e), &
crystallite_Fe(1:3,1:3,g,i,e), &
crystallite_Fp(1:3,1:3,g,i,e), & crystallite_Fp(1:3,1:3,g,i,e), &
g, i, e) g, i, e)
enddo; enddo; enddo enddo; enddo; enddo

View File

@ -3005,10 +3005,9 @@ end subroutine plastic_nonlocal_dotState
!* that sum up to a total of 1 are considered, all others are set to * !* that sum up to a total of 1 are considered, all others are set to *
!* zero. * !* zero. *
!********************************************************************* !*********************************************************************
subroutine plastic_nonlocal_updateCompatibility(orientation,i,e) subroutine plastic_nonlocal_updateCompatibility(orientation,i,e)
use math, only: math_mul3x3, math_qRot
use math, only: math_mul3x3, & use rotations, only: rotation
math_qRot
use material, only: material_phase, & use material, only: material_phase, &
material_texture, & material_texture, &
phase_localPlasticity, & phase_localPlasticity, &
@ -3030,7 +3029,7 @@ implicit none
!* input variables !* input variables
integer(pInt), intent(in) :: i, & ! ip index integer(pInt), intent(in) :: i, & ! ip index
e ! element 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 orientation ! crystal orientation in quaternions
!* local variables !* local variables
@ -3059,7 +3058,7 @@ real(pReal) my_compatibilitySum, &
nThresholdValues nThresholdValues
logical, dimension(totalNslip(phase_plasticityInstance(material_phase(1,i,e)))) :: & logical, dimension(totalNslip(phase_plasticityInstance(material_phase(1,i,e)))) :: &
belowThreshold belowThreshold
type(rotation) :: rot
Nneighbors = FE_NipNeighbors(FE_celltype(FE_geomtype(mesh_element(2,e)))) Nneighbors = FE_NipNeighbors(FE_celltype(FE_geomtype(mesh_element(2,e))))
ph = material_phase(1,i,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. !* 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. !* All values below the threshold are set to zero.
else else
absoluteMisorientation = lattice_qDisorientation(orientation(1:4,1,i,e), & rot = orientation(1,i,e)%misorientation(orientation(1,neighbor_i,neighbor_e))
orientation(1:4,1,neighbor_i,neighbor_e)) ! no symmetry absoluteMisorientation = rot%asQuaternion()
mySlipSystems: do s1 = 1_pInt,ns mySlipSystems: do s1 = 1_pInt,ns
neighborSlipSystems: do s2 = 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))) & my_compatibility(1,s2,s1,n) = math_mul3x3(slipNormal(1:3,s1), math_qRot(absoluteMisorientation, slipNormal(1:3,s2))) &