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

use rotation class

This commit is contained in:
Martin Diehl 2019-12-02 20:23:50 +01:00
parent 8a9d3f8d6d
commit f5292019e5
5 changed files with 23 additions and 59 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
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)
@ -255,7 +254,7 @@ function grid_mech_FEM_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation
!--------------------------------------------------------------------------------------------------
! update stiffness (and gamma operator)
S = utilities_maskedCompliance(rotation_BC%asMatrix(),stress_BC%maskLogical,C_volAvg)
S = utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg)
!--------------------------------------------------------------------------------------------------
! set module wide available data
params%stress_mask = stress_BC%maskFloat
@ -507,7 +506,7 @@ subroutine formResidual(da_local,x_local, &
! evaluate constitutive response
call Utilities_constitutiveResponse(P_current,&
P_av,C_volAvg,devNull, &
F,params%timeinc,params%rotation_BC%asMatrix())
F,params%timeinc,params%rotation_BC)
call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr)
!--------------------------------------------------------------------------------------------------

7
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
@ -225,7 +224,7 @@ function grid_mech_spectral_basic_solution(incInfoIn,timeinc,timeinc_old,stress_
!--------------------------------------------------------------------------------------------------
! update stiffness (and gamma operator)
S = utilities_maskedCompliance(rotation_BC%asMatrix(),stress_BC%maskLogical,C_volAvg)
S = utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg)
if(num%update_gamma) call utilities_updateGamma(C_minMaxAvg)
!--------------------------------------------------------------------------------------------------
@ -457,7 +456,7 @@ subroutine formResidual(in, F, &
! evaluate constitutive response
call utilities_constitutiveResponse(residuum, & ! "residuum" gets field of first PK stress (to save memory)
P_av,C_volAvg,C_minMaxAvg, &
F,params%timeinc,params%rotation_BC%asMatrix())
F,params%timeinc,params%rotation_BC)
call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr)
!--------------------------------------------------------------------------------------------------

7
src/grid/grid_mech_spectral_polarisation.f90
View File

@ -187,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
@ -241,7 +240,7 @@ function grid_mech_spectral_polarisation_solution(incInfoIn,timeinc,timeinc_old,
!--------------------------------------------------------------------------------------------------
! update stiffness (and gamma operator)
S = utilities_maskedCompliance(rotation_BC%asMatrix(),stress_BC%maskLogical,C_volAvg)
S = utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg)
if (num%update_gamma) then
call utilities_updateGamma(C_minMaxAvg)
C_scale = C_minMaxAvg
@ -546,7 +545,7 @@ subroutine formResidual(in, FandF_tau, &
! evaluate constitutive response
call utilities_constitutiveResponse(residual_F, & ! "residuum" gets field of first PK stress (to save memory)
P_av,C_volAvg,C_minMaxAvg, &
F - residual_F_tau/polarBeta,params%timeinc,params%rotation_BC%asMatrix())
F - residual_F_tau/polarBeta,params%timeinc,params%rotation_BC)
call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr)
!--------------------------------------------------------------------------------------------------

18
src/grid/spectral_utilities.f90
View File

@ -687,9 +687,10 @@ 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
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
@ -707,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 ............................................'
@ -836,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(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
type(rotation), intent(in), optional :: rotation_BC !< rotation of load frame
integer :: &
@ -863,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)

35
src/math.f90
View File

@ -1301,41 +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 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