diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90 index e83cf3283..51c97456b 100644 --- a/src/grid/DAMASK_grid.f90 +++ b/src/grid/DAMASK_grid.f90 @@ -28,7 +28,6 @@ program DAMASK_spectral use grid_damage_spectral use grid_thermal_spectral use results - use rotations implicit none @@ -78,7 +77,6 @@ program DAMASK_spectral character(len=6) :: loadcase_string character(len=1024) :: & incInfo - type(rotation) :: R type(tLoadCase), allocatable, dimension(:) :: loadCases !< array of all load cases type(tLoadCase) :: newLoadCase type(tSolutionState), allocatable, dimension(:) :: solres @@ -189,6 +187,7 @@ program DAMASK_spectral newLoadCase%ID(field) = FIELD_DAMAGE_ID endif damageActive + call newLoadCase%rot%fromEulers(real([0.0,0.0,0.0],pReal)) readIn: do i = 1, chunkPos(1) select case (IO_lc(IO_stringValue(line,chunkPos,i))) case('fdot','dotf','l','f') ! assign values for the deformation BC matrix @@ -244,14 +243,13 @@ program DAMASK_spectral do j = 1, 3 temp_valueVector(j) = IO_floatValue(line,chunkPos,i+k+j) enddo - call R%fromEulers(temp_valueVector(1:3),degrees=(l==1)) - newLoadCase%rotation = R%asMatrix() + call newLoadCase%rot%fromEulers(temp_valueVector(1:3),degrees=(l==1)) case('rotation','rot') ! assign values for the rotation matrix temp_valueVector = 0.0_pReal do j = 1, 9 temp_valueVector(j) = IO_floatValue(line,chunkPos,i+j) enddo - newLoadCase%rotation = math_9to33(temp_valueVector) + call newLoadCase%rot%fromMatrix(math_9to33(temp_valueVector)) end select enddo readIn @@ -295,14 +293,12 @@ program DAMASK_spectral endif enddo; write(6,'(/)',advance='no') enddo - if (any(abs(matmul(newLoadCase%rotation, & - transpose(newLoadCase%rotation))-math_I3) > & - reshape(spread(tol_math_check,1,9),[ 3,3]))& - .or. abs(math_det33(newLoadCase%rotation)) > & - 1.0_pReal + tol_math_check) errorID = 846 ! given rotation matrix contains strain - if (any(dNeq(newLoadCase%rotation, math_I3))) & + if (any(abs(matmul(newLoadCase%rot%asMatrix(), & + transpose(newLoadCase%rot%asMatrix()))-math_I3) > & + reshape(spread(tol_math_check,1,9),[ 3,3]))) errorID = 846 ! given rotation matrix contains strain + if (any(dNeq(newLoadCase%rot%asMatrix(), math_I3))) & write(6,'(2x,a,/,3(3(3x,f12.7,1x)/))',advance='no') 'rotation of loadframe:',& - transpose(newLoadCase%rotation) + transpose(newLoadCase%rot%asMatrix()) if (newLoadCase%time < 0.0_pReal) errorID = 834 ! negative time increment write(6,'(2x,a,f12.6)') 'time: ', newLoadCase%time if (newLoadCase%incs < 1) errorID = 835 ! non-positive incs count @@ -469,7 +465,7 @@ program DAMASK_spectral cutBack,guess,timeinc,timeIncOld,remainingLoadCaseTime, & deformation_BC = loadCases(currentLoadCase)%deformation, & stress_BC = loadCases(currentLoadCase)%stress, & - rotation_BC = loadCases(currentLoadCase)%rotation) + rotation_BC = loadCases(currentLoadCase)%rot) case(FIELD_THERMAL_ID); call grid_thermal_spectral_forward(cutBack) case(FIELD_DAMAGE_ID); call grid_damage_spectral_forward(cutBack) @@ -488,7 +484,7 @@ program DAMASK_spectral solres(field) = mech_solution (& incInfo,timeinc,timeIncOld, & stress_BC = loadCases(currentLoadCase)%stress, & - rotation_BC = loadCases(currentLoadCase)%rotation) + rotation_BC = loadCases(currentLoadCase)%rot) case(FIELD_THERMAL_ID) solres(field) = grid_thermal_spectral_solution(timeinc,timeIncOld) diff --git a/src/grid/grid_mech_FEM.f90 b/src/grid/grid_mech_FEM.f90 index 0c3844fcf..e87e02d09 100644 --- a/src/grid/grid_mech_FEM.f90 +++ b/src/grid/grid_mech_FEM.f90 @@ -207,8 +207,7 @@ subroutine grid_mech_FEM_init call utilities_updateCoords(F) call utilities_constitutiveResponse(P_current,temp33_Real,C_volAvg,devNull, & ! stress field, stress avg, global average of stiffness and (min+max)/2 F, & ! target F - 0.0_pReal, & ! time increment - math_I3) ! no rotation of boundary condition + 0.0_pReal) ! time increment call DMDAVecRestoreArrayF90(mech_grid,solution_current,u_current,ierr) CHKERRQ(ierr) call DMDAVecRestoreArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr) @@ -242,7 +241,8 @@ function grid_mech_FEM_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation timeinc_old !< time increment of last successful increment type(tBoundaryCondition), intent(in) :: & stress_BC - real(pReal), dimension(3,3), intent(in) :: rotation_BC + type(rotation), intent(in) :: & + rotation_BC type(tSolutionState) :: & solution !-------------------------------------------------------------------------------------------------- @@ -297,7 +297,7 @@ subroutine grid_mech_FEM_forward(cutBack,guess,timeinc,timeinc_old,loadCaseTime, type(tBoundaryCondition), intent(in) :: & stress_BC, & deformation_BC - real(pReal), dimension(3,3), intent(in) :: & + type(rotation), intent(in) :: & rotation_BC PetscErrorCode :: ierr PetscScalar, pointer, dimension(:,:,:,:) :: & @@ -482,7 +482,7 @@ subroutine formResidual(da_local,x_local, & trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter+1, '≤', itmax if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) & write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & - ' deformation gradient aim (lab) =', transpose(math_rotate_backward33(F_aim,params%rotation_BC)) + ' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotTensor2(F_aim,active=.true.)) write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & ' deformation gradient aim =', transpose(F_aim) flush(6) @@ -498,7 +498,7 @@ subroutine formResidual(da_local,x_local, & x_elem(ctr,1:3) = x_scal(0:2,i+ii,j+jj,k+kk) enddo; enddo; enddo ii = i-xstart+1; jj = j-ystart+1; kk = k-zstart+1 - F(1:3,1:3,ii,jj,kk) = math_rotate_backward33(F_aim,params%rotation_BC) + transpose(matmul(BMat,x_elem)) + F(1:3,1:3,ii,jj,kk) = params%rotation_BC%rotTensor2(F_aim,active=.true.) + transpose(matmul(BMat,x_elem)) enddo; enddo; enddo call DMDAVecRestoreArrayF90(da_local,x_local,x_scal,ierr);CHKERRQ(ierr) diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index 3dc978dc6..36caabf90 100644 --- a/src/grid/grid_mech_spectral_basic.f90 +++ b/src/grid/grid_mech_spectral_basic.f90 @@ -173,8 +173,7 @@ subroutine grid_mech_spectral_basic_init call Utilities_updateCoords(reshape(F,shape(F_lastInc))) call Utilities_constitutiveResponse(P,temp33_Real,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2 reshape(F,shape(F_lastInc)), & ! target F - 0.0_pReal, & ! time increment - math_I3) ! no rotation of boundary condition + 0.0_pReal) ! time increment call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! deassociate pointer restartRead2: if (interface_restartInc > 0) then @@ -212,7 +211,8 @@ function grid_mech_spectral_basic_solution(incInfoIn,timeinc,timeinc_old,stress_ timeinc_old !< time increment of last successful increment type(tBoundaryCondition), intent(in) :: & stress_BC - real(pReal), dimension(3,3), intent(in) :: rotation_BC + type(rotation), intent(in) :: & + rotation_BC type(tSolutionState) :: & solution !-------------------------------------------------------------------------------------------------- @@ -269,7 +269,7 @@ subroutine grid_mech_spectral_basic_forward(cutBack,guess,timeinc,timeinc_old,lo type(tBoundaryCondition), intent(in) :: & stress_BC, & deformation_BC - real(pReal), dimension(3,3), intent(in) :: & + type(rotation), intent(in) :: & rotation_BC PetscErrorCode :: ierr PetscScalar, dimension(:,:,:,:), pointer :: F @@ -299,9 +299,9 @@ subroutine grid_mech_spectral_basic_forward(cutBack,guess,timeinc,timeinc_old,lo F_aimDot + deformation_BC%maskFloat * (deformation_BC%values - F_aim_lastInc)/loadCaseTime endif - Fdot = utilities_calculateRate(guess, & - F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid3]),timeinc_old, & - math_rotate_backward33(F_aimDot,rotation_BC)) + Fdot = utilities_calculateRate(guess, & + F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid3]),timeinc_old, & + rotation_BC%rotTensor2(F_aimDot,active=.true.)) F_lastInc = reshape(F,[3,3,grid(1),grid(2),grid3]) materialpoint_F0 = reshape(F, [3,3,1,product(grid(1:2))*grid3]) @@ -311,7 +311,7 @@ subroutine grid_mech_spectral_basic_forward(cutBack,guess,timeinc,timeinc_old,lo ! update average and local deformation gradients F_aim = F_aim_lastInc + F_aimDot * timeinc F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot, & ! estimate of F at end of time+timeinc that matches rotated F_aim on average - math_rotate_backward33(F_aim,rotation_BC)),[9,grid(1),grid(2),grid3]) + rotation_BC%rotTensor2(F_aim,active=.true.)),[9,grid(1),grid(2),grid3]) call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) end subroutine grid_mech_spectral_basic_forward @@ -446,7 +446,7 @@ subroutine formResidual(in, F, & trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) & write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & - ' deformation gradient aim (lab) =', transpose(math_rotate_backward33(F_aim,params%rotation_BC)) + ' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotTensor2(F_aim,active=.true.)) write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & ' deformation gradient aim =', transpose(F_aim) flush(6) @@ -471,7 +471,7 @@ subroutine formResidual(in, F, & tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = residuum ! store fPK field for subsequent FFT forward transform call utilities_FFTtensorForward ! FFT forward of global "tensorField_real" err_div = Utilities_divergenceRMS() ! divRMS of tensorField_fourier for later use - call utilities_fourierGammaConvolution(math_rotate_backward33(deltaF_aim,params%rotation_BC)) ! convolution of Gamma and tensorField_fourier, with arg + call utilities_fourierGammaConvolution(params%rotation_BC%rotTensor2(deltaF_aim,active=.true.)) ! convolution of Gamma and tensorField_fourier call utilities_FFTtensorBackward ! FFT backward of global tensorField_fourier !-------------------------------------------------------------------------------------------------- diff --git a/src/grid/grid_mech_spectral_polarisation.f90 b/src/grid/grid_mech_spectral_polarisation.f90 index c1b5d79c9..07d9a5afc 100644 --- a/src/grid/grid_mech_spectral_polarisation.f90 +++ b/src/grid/grid_mech_spectral_polarisation.f90 @@ -14,6 +14,7 @@ module grid_mech_spectral_polarisation use DAMASK_interface use HDF5_utilities use math + use rotations use spectral_utilities use IO use FEsolving @@ -186,8 +187,7 @@ subroutine grid_mech_spectral_polarisation_init call Utilities_updateCoords(reshape(F,shape(F_lastInc))) call Utilities_constitutiveResponse(P,temp33_Real,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2 reshape(F,shape(F_lastInc)), & ! target F - 0.0_pReal, & ! time increment - math_I3) ! no rotation of boundary condition + 0.0_pReal) ! time increment call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) ! deassociate pointer restartRead2: if (interface_restartInc > 0) then @@ -222,12 +222,13 @@ function grid_mech_spectral_polarisation_solution(incInfoIn,timeinc,timeinc_old, ! input data for solution character(len=*), intent(in) :: & incInfoIn - real(pReal), intent(in) :: & + real(pReal), intent(in) :: & timeinc, & !< time increment of current solution timeinc_old !< time increment of last successful increment type(tBoundaryCondition), intent(in) :: & stress_BC - real(pReal), dimension(3,3), intent(in) :: rotation_BC + type(rotation), intent(in) :: & + rotation_BC type(tSolutionState) :: & solution !-------------------------------------------------------------------------------------------------- @@ -288,7 +289,7 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc type(tBoundaryCondition), intent(in) :: & stress_BC, & deformation_BC - real(pReal), dimension(3,3), intent(in) ::& + type(rotation), intent(in) :: & rotation_BC PetscErrorCode :: ierr PetscScalar, dimension(:,:,:,:), pointer :: FandF_tau, F, F_tau @@ -324,10 +325,10 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc Fdot = utilities_calculateRate(guess, & F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid3]),timeinc_old, & - math_rotate_backward33(F_aimDot,rotation_BC)) + rotation_BC%rotTensor2(F_aimDot,active=.true.)) F_tauDot = utilities_calculateRate(guess, & F_tau_lastInc,reshape(F_tau,[3,3,grid(1),grid(2),grid3]), timeinc_old, & - math_rotate_backward33(F_aimDot,rotation_BC)) + rotation_BC%rotTensor2(F_aimDot,active=.true.)) F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid3]) F_tau_lastInc = reshape(F_tau,[3,3,grid(1),grid(2),grid3]) @@ -338,7 +339,7 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc ! update average and local deformation gradients F_aim = F_aim_lastInc + F_aimDot * timeinc F = reshape(utilities_forwardField(timeinc,F_lastInc,Fdot, & ! estimate of F at end of time+timeinc that matches rotated F_aim on average - math_rotate_backward33(F_aim,rotation_BC)),& + rotation_BC%rotTensor2(F_aim,active=.true.)),& [9,grid(1),grid(2),grid3]) if (guess) then F_tau = reshape(Utilities_forwardField(timeinc,F_tau_lastInc,F_taudot), & @@ -349,8 +350,8 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc F_lambda33 = math_mul3333xx33(S_scale,matmul(F_lambda33, & math_mul3333xx33(C_scale,& matmul(transpose(F_lambda33),& - F_lambda33)-math_I3))*0.5_pReal)& - + math_I3 + F_lambda33)-math_I3))*0.5_pReal) & + + math_I3 F_tau(1:9,i,j,k) = reshape(F_lambda33,[9])+F(1:9,i,j,k) enddo; enddo; enddo endif @@ -514,7 +515,7 @@ subroutine formResidual(in, FandF_tau, & trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) & write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & - ' deformation gradient aim (lab) =', transpose(math_rotate_backward33(F_aim,params%rotation_BC)) + ' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotTensor2(F_aim,active=.true.)) write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & ' deformation gradient aim =', transpose(F_aim) flush(6) @@ -533,7 +534,7 @@ subroutine formResidual(in, FandF_tau, & !-------------------------------------------------------------------------------------------------- ! doing convolution in Fourier space call utilities_FFTtensorForward - call utilities_fourierGammaConvolution(math_rotate_backward33(polarBeta*F_aim,params%rotation_BC)) + call utilities_fourierGammaConvolution(params%rotation_BC%rotTensor2(polarBeta*F_aim,active=.true.)) call utilities_FFTtensorBackward !-------------------------------------------------------------------------------------------------- @@ -551,7 +552,7 @@ subroutine formResidual(in, FandF_tau, & ! stress BC handling F_aim = F_aim - math_mul3333xx33(S, ((P_av - params%stress_BC))) ! S = 0.0 for no bc err_BC = maxval(abs((1.0_pReal-params%stress_mask) * math_mul3333xx33(C_scale,F_aim & - -math_rotate_forward33(F_av,params%rotation_BC)) + & + -params%rotation_BC%rotTensor2(F_av)) + & params%stress_mask * (P_av-params%stress_BC))) ! mask = 0.0 for no bc ! calculate divergence tensorField_real = 0.0_pReal diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index 8b0e430f3..a1265556a 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -10,6 +10,7 @@ module spectral_utilities use prec use math + use rotations use IO use mesh_grid use numerics @@ -90,7 +91,7 @@ module spectral_utilities end type tBoundaryCondition type, public :: tLoadCase - real(pReal), dimension (3,3) :: rotation = math_I3 !< rotation of BC + type(rotation) :: rot !< rotation of BC type(tBoundaryCondition) :: stress, & !< stress BC deformation !< deformation BC (Fdot or L) real(pReal) :: time = 0.0_pReal !< length of increment @@ -103,7 +104,8 @@ module spectral_utilities end type tLoadCase type, public :: tSolutionParams !< @todo use here the type definition for a full loadcase - real(pReal), dimension(3,3) :: stress_mask, stress_BC, rotation_BC + real(pReal), dimension(3,3) :: stress_mask, stress_BC + type(rotation) :: rotation_BC real(pReal) :: timeinc real(pReal) :: timeincOld end type tSolutionParams @@ -684,10 +686,11 @@ end function utilities_curlRMS !-------------------------------------------------------------------------------------------------- function utilities_maskedCompliance(rot_BC,mask_stress,C) - real(pReal), dimension(3,3,3,3) :: utilities_maskedCompliance !< masked compliance - real(pReal), intent(in) , dimension(3,3,3,3) :: C !< current average stiffness - real(pReal), intent(in) , dimension(3,3) :: rot_BC !< rotation of load frame - logical, intent(in), dimension(3,3) :: mask_stress !< mask of stress BC + real(pReal), dimension(3,3,3,3) :: utilities_maskedCompliance !< masked compliance + real(pReal), intent(in), dimension(3,3,3,3) :: C !< current average stiffness + type(rotation), intent(in) :: rot_BC !< rotation of load frame + logical, intent(in), dimension(3,3) :: mask_stress !< mask of stress BC + integer :: j, k, m, n logical, dimension(9) :: mask_stressVector real(pReal), dimension(9,9) :: temp99_Real @@ -705,7 +708,7 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C) allocate (c_reduced(size_reduced,size_reduced), source =0.0_pReal) allocate (s_reduced(size_reduced,size_reduced), source =0.0_pReal) allocate (sTimesC(size_reduced,size_reduced), source =0.0_pReal) - temp99_Real = math_3333to99(math_rotate_forward3333(C,rot_BC)) + temp99_Real = math_3333to99(rot_BC%rotTensor4(C)) if(debugGeneral) then write(6,'(/,a)') ' ... updating masked compliance ............................................' @@ -834,12 +837,12 @@ end subroutine utilities_fourierTensorDivergence subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,& F,timeinc,rotation_BC) - real(pReal),intent(out), dimension(3,3,3,3) :: C_volAvg, C_minmaxAvg !< average stiffness - real(pReal),intent(out), dimension(3,3) :: P_av !< average PK stress - real(pReal),intent(out), dimension(3,3,grid(1),grid(2),grid3) :: P !< PK stress - real(pReal), intent(in), dimension(3,3,grid(1),grid(2),grid3) :: F !< deformation gradient target - real(pReal), intent(in) :: timeinc !< loading time - real(pReal), intent(in), dimension(3,3) :: rotation_BC !< rotation of load frame + real(pReal), intent(out), dimension(3,3,3,3) :: C_volAvg, C_minmaxAvg !< average stiffness + real(pReal), intent(out), dimension(3,3) :: P_av !< average PK stress + real(pReal), intent(out), dimension(3,3,grid(1),grid(2),grid3) :: P !< PK stress + real(pReal), intent(in), dimension(3,3,grid(1),grid(2),grid3) :: F !< deformation gradient target + real(pReal), intent(in) :: timeinc !< loading time + type(rotation), intent(in), optional :: rotation_BC !< rotation of load frame integer :: & @@ -861,7 +864,8 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,& if (debugRotation) & write(6,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress (lab) / MPa =',& transpose(P_av)*1.e-6_pReal - P_av = math_rotate_forward33(P_av,rotation_BC) + if(present(rotation_BC)) & + P_av = rotation_BC%rotTensor2(P_av) write(6,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress / MPa =',& transpose(P_av)*1.e-6_pReal flush(6) diff --git a/src/math.f90 b/src/math.f90 index 7eba379d2..0b06c9186 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -838,33 +838,6 @@ pure function math_Voigt66to3333(m66) end function math_Voigt66to3333 -!-------------------------------------------------------------------------------------------------- -!> @brief action of a quaternion on a vector (rotate vector v with Q) -!> @details deprecated -!-------------------------------------------------------------------------------------------------- -pure function math_qRot(Q,v) - - real(pReal), dimension(4), intent(in) :: Q - real(pReal), dimension(3), intent(in) :: v - real(pReal), dimension(3) :: math_qRot - real(pReal), dimension(4,4) :: T - integer :: i, j - - do i = 1,4 - do j = 1,i - T(i,j) = Q(i) * Q(j) - enddo - enddo - - math_qRot = [-v(1)*(T(3,3)+T(4,4)) + v(2)*(T(3,2)-T(4,1)) + v(3)*(T(4,2)+T(3,1)), & - v(1)*(T(3,2)+T(4,1)) - v(2)*(T(2,2)+T(4,4)) + v(3)*(T(4,3)-T(2,1)), & - v(1)*(T(4,2)-T(3,1)) + v(2)*(T(4,3)+T(2,1)) - v(3)*(T(2,2)+T(3,3))] - - math_qRot = 2.0_pReal * math_qRot + v - -end function math_qRot - - !-------------------------------------------------------------------------------------------------- !> @brief rotation matrix from Bunge-Euler (3-1-3) angles (in radians) !> @details deprecated @@ -1328,55 +1301,6 @@ real(pReal) pure function math_areaTriangle(v1,v2,v3) end function math_areaTriangle -!-------------------------------------------------------------------------------------------------- -!> @brief rotate 33 tensor forward -!> @details deprecated -!-------------------------------------------------------------------------------------------------- -pure function math_rotate_forward33(tensor,R) - - real(pReal), dimension(3,3) :: math_rotate_forward33 - real(pReal), dimension(3,3), intent(in) :: tensor, R - - math_rotate_forward33 = matmul(R,matmul(tensor,transpose(R))) - -end function math_rotate_forward33 - - -!-------------------------------------------------------------------------------------------------- -!> @brief rotate 33 tensor backward -!> @details deprecated -!-------------------------------------------------------------------------------------------------- -pure function math_rotate_backward33(tensor,R) - - real(pReal), dimension(3,3) :: math_rotate_backward33 - real(pReal), dimension(3,3), intent(in) :: tensor, R - - math_rotate_backward33 = matmul(transpose(R),matmul(tensor,R)) - -end function math_rotate_backward33 - - -!-------------------------------------------------------------------------------------------------- -!> @brief rotate 3333 tensor C'_ijkl=g_im*g_jn*g_ko*g_lp*C_mnop -!> @details deprecated -!-------------------------------------------------------------------------------------------------- -pure function math_rotate_forward3333(tensor,R) - - real(pReal), dimension(3,3,3,3) :: math_rotate_forward3333 - real(pReal), dimension(3,3), intent(in) :: R - real(pReal), dimension(3,3,3,3), intent(in) :: tensor - integer :: i,j,k,l,m,n,o,p - - math_rotate_forward3333 = 0.0_pReal - do i = 1,3;do j = 1,3;do k = 1,3;do l = 1,3 - do m = 1,3;do n = 1,3;do o = 1,3;do p = 1,3 - math_rotate_forward3333(i,j,k,l) = math_rotate_forward3333(i,j,k,l) & - + R(i,m) * R(j,n) * R(k,o) * R(l,p) * tensor(m,n,o,p) - enddo; enddo; enddo; enddo; enddo; enddo; enddo; enddo - -end function math_rotate_forward3333 - - !-------------------------------------------------------------------------------------------------- !> @brief limits a scalar value to a certain range (either one or two sided) ! Will return NaN if left > right diff --git a/src/plastic_nonlocal.f90 b/src/plastic_nonlocal.f90 index 2e6734c9d..5375aba49 100644 --- a/src/plastic_nonlocal.f90 +++ b/src/plastic_nonlocal.f90 @@ -24,10 +24,10 @@ module plastic_nonlocal implicit none private - real(pReal), parameter, private :: & + real(pReal), parameter :: & KB = 1.38e-23_pReal !< Physical parameter, Boltzmann constant in J/Kelvin - character(len=64), dimension(:,:), allocatable, target, public :: & + character(len=64), dimension(:,:), allocatable :: & plastic_nonlocal_output !< name of each post result output ! storage order of dislocation types @@ -50,18 +50,18 @@ module plastic_nonlocal mob_scr_neg = 4 !< mobile screw positive ! BEGIN DEPRECATES - integer, dimension(:,:,:), allocatable, private :: & + integer, dimension(:,:,:), allocatable :: & iRhoU, & !< state indices for unblocked density iRhoB, & !< state indices for blocked density iRhoD, & !< state indices for dipole density iV, & !< state indices for dislcation velocities iD !< state indices for stable dipole height - integer, dimension(:), allocatable, private, protected :: & + integer, dimension(:), allocatable :: & totalNslip !< total number of active slip systems for each instance !END DEPRECATED - real(pReal), dimension(:,:,:,:,:,:), allocatable, private :: & - compatibility !< slip system compatibility between me and my neighbors + real(pReal), dimension(:,:,:,:,:,:), allocatable :: & + compatibility !< slip system compatibility between me and my neighbors enum, bind(c) enumerator :: & @@ -89,7 +89,7 @@ module plastic_nonlocal gamma_ID end enum - type, private :: tParameters !< container type for internal constitutive parameters + type :: tParameters !< container type for internal constitutive parameters real(pReal) :: & atomicVolume, & !< atomic volume Dsd0, & !< prefactor for self-diffusion coefficient @@ -139,19 +139,19 @@ module plastic_nonlocal interactionSlipSlip ,& !< coefficients for slip-slip interaction forestProjection_Edge, & !< matrix of forest projections of edge dislocations forestProjection_Screw !< matrix of forest projections of screw dislocations - real(pReal), dimension(:), allocatable, private :: & + real(pReal), dimension(:), allocatable :: & nonSchmidCoeff - real(pReal), dimension(:,:,:), allocatable, private :: & + real(pReal), dimension(:,:,:), allocatable :: & Schmid, & !< Schmid contribution nonSchmid_pos, & nonSchmid_neg !< combined projection of Schmid and non-Schmid contributions to the resolved shear stress (only for screws) integer :: & totalNslip - integer, dimension(:) ,allocatable , public:: & + integer, dimension(:) ,allocatable :: & Nslip,& colinearSystem !< colinear system to the active slip system (only valid for fcc!) - logical, private :: & + logical :: & shortRangeStressCorrection, & !< flag indicating the use of the short range stress correction by a excess density gradient term probabilisticMultiplication @@ -160,13 +160,13 @@ module plastic_nonlocal end type tParameters - type, private :: tNonlocalMicrostructure + type :: tNonlocalMicrostructure real(pReal), allocatable, dimension(:,:) :: & tau_pass, & tau_Back end type tNonlocalMicrostructure - type, private :: tNonlocalState + type :: tNonlocalState real(pReal), pointer, dimension(:,:) :: & rho, & ! < all dislocations rhoSgl, & @@ -192,16 +192,16 @@ module plastic_nonlocal v_scr_neg end type tNonlocalState - type(tNonlocalState), allocatable, dimension(:), private :: & + type(tNonlocalState), allocatable, dimension(:) :: & deltaState, & dotState, & state - type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance) + type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance) - type(tNonlocalMicrostructure), dimension(:), allocatable, private :: microstructure + type(tNonlocalMicrostructure), dimension(:), allocatable :: microstructure - integer(kind(undefined_ID)), dimension(:,:), allocatable, private :: & + integer(kind(undefined_ID)), dimension(:,:), allocatable :: & plastic_nonlocal_outputID !< ID of each post result output public :: & @@ -1829,8 +1829,6 @@ subroutine plastic_nonlocal_updateCompatibility(orientation,i,e) ns, & ! number of active slip systems s1, & ! slip system index (me) s2 ! slip system index (my neighbor) - real(pReal), dimension(4) :: & - absoluteMisorientation ! absolute misorientation (without symmetry) between me and my neighbor real(pReal), dimension(2,totalNslip(phase_plasticityInstance(material_phaseAt(1,e))),& totalNslip(phase_plasticityInstance(material_phaseAt(1,e))),& nIPneighbors) :: & @@ -1841,7 +1839,7 @@ subroutine plastic_nonlocal_updateCompatibility(orientation,i,e) nThresholdValues logical, dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,e)))) :: & belowThreshold - type(rotation) :: rot + type(rotation) :: mis Nneighbors = nIPneighbors ph = material_phaseAt(1,e) @@ -1907,18 +1905,17 @@ subroutine plastic_nonlocal_updateCompatibility(orientation,i,e) !* Finally the smallest compatibility value is decreased until the sum is exactly equal to one. !* All values below the threshold are set to zero. else - rot = orientation(1,i,e)%misorientation(orientation(1,neighbor_i,neighbor_e)) - absoluteMisorientation = rot%asQuaternion() + mis = orientation(1,i,e)%misorientation(orientation(1,neighbor_i,neighbor_e)) mySlipSystems: do s1 = 1,ns neighborSlipSystems: do s2 = 1,ns my_compatibility(1,s2,s1,n) = math_inner(prm%slip_normal(1:3,s1), & - math_qRot(absoluteMisorientation, prm%slip_normal(1:3,s2))) & + mis%rotate(prm%slip_normal(1:3,s2))) & * abs(math_inner(prm%slip_direction(1:3,s1), & - math_qRot(absoluteMisorientation, prm%slip_direction(1:3,s2)))) + mis%rotate(prm%slip_direction(1:3,s2)))) my_compatibility(2,s2,s1,n) = abs(math_inner(prm%slip_normal(1:3,s1), & - math_qRot(absoluteMisorientation, prm%slip_normal(1:3,s2)))) & + mis%rotate(prm%slip_normal(1:3,s2)))) & * abs(math_inner(prm%slip_direction(1:3,s1), & - math_qRot(absoluteMisorientation, prm%slip_direction(1:3,s2)))) + mis%rotate(prm%slip_direction(1:3,s2)))) enddo neighborSlipSystems my_compatibilitySum = 0.0_pReal diff --git a/src/rotations.f90 b/src/rotations.f90 index e042cff21..a4a0bac88 100644 --- a/src/rotations.f90 +++ b/src/rotations.f90 @@ -64,6 +64,7 @@ module rotations procedure, public :: asRodrigues procedure, public :: asMatrix !------------------------------------------ + procedure, public :: fromQuaternion procedure, public :: fromEulers procedure, public :: fromAxisAngle procedure, public :: fromMatrix @@ -157,6 +158,18 @@ end function asHomochoric !--------------------------------------------------------------------------------------------------- ! Initialize rotation from different representations !--------------------------------------------------------------------------------------------------- +subroutine fromQuaternion(self,qu) + + class(rotation), intent(out) :: self + real(pReal), dimension(4), intent(in) :: qu + + if (dNeq(norm2(qu),1.0_pReal)) & + call IO_error(402,ext_msg='fromQuaternion') + + self%q = qu + +end subroutine fromQuaternion +!--------------------------------------------------------------------------------------------------- subroutine fromEulers(self,eu,degrees) class(rotation), intent(out) :: self