223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

Merge branch 'only-use-rotation-class' into 'development' 

Only use rotation class

See merge request damask/DAMASK!110
This commit is contained in:
Vitesh 2019-12-07 17:33:34 +01:00
parent 0f4f415c5b 04492eeb0a
commit 3e269f0419
8 changed files with 94 additions and 159 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

24
src/grid/DAMASK_grid.f90
View File

@ -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)

12
src/grid/grid_mech_FEM.f90
View File

@ -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)

20
src/grid/grid_mech_spectral_basic.f90
View File

@ -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
!--------------------------------------------------------------------------------------------------

27
src/grid/grid_mech_spectral_polarisation.f90
View File

@ -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

32
src/grid/spectral_utilities.f90
View File

@ -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)

76
src/math.f90
View File

@ -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

49
src/plastic_nonlocal.f90
View File

@ -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

13
src/rotations.f90
View File

@ -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