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

Merge branch '38-introduce-rudimentary-PETSc-based-FEM-solver' of magit1.mpie.de:damask/DAMASK into 38-introduce-rudimentary-PETSc-based-FEM-solver

This commit is contained in:
Jaeyong Jung 2018-09-20 09:57:30 +02:00
parent 3beab556a2 cc262ae198
commit abe5c1d825
10 changed files with 380 additions and 498 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
VERSION
View File

@ -1 +1 @@
v2.0.2-540-gce6e6679 v2.0.2-565-g59043d58

11
examples/ConfigFiles/Crystallite_All.config
View File

@ -5,15 +5,10 @@
(output) orientation # quaternion (output) orientation # quaternion
(output) eulerangles # orientation as Bunge triple in degree (output) eulerangles # orientation as Bunge triple in degree
(output) grainrotation # deviation from initial orientation as axis (1-3) and angle in degree (4) in crystal reference coordinates (output) grainrotation # deviation from initial orientation as axis (1-3) and angle in degree (4) in crystal reference coordinates
(output) grainrotationx # deviation from initial orientation as angle in degrees around sample reference x axis (output) f # deformation gradient tensor
(output) grainrotationy # deviation from initial orientation as angle in degrees around sample reference y axis
(output) grainrotationz # deviation from initial orientation as angle in degrees around sample reference z axis
(output) f # deformation gradient tensor; synonyms: "defgrad"
(output) fe # elastic deformation gradient tensor (output) fe # elastic deformation gradient tensor
(output) fp # plastic deformation gradient tensor (output) fp # plastic deformation gradient tensor
(output) e # total strain as Green-Lagrange tensor (output) p # first Piola-Kichhoff stress tensor
(output) ee # elastic strain as Green-Lagrange tensor (output) s # second Piola-Kichhoff stress tensor
(output) p # first Piola-Kichhoff stress tensor; synonyms: "firstpiola", "1stpiola"
(output) s # second Piola-Kichhoff stress tensor; synonyms: "tstar", "secondpiola", "2ndpiola"
(output) lp # plastic velocity gradient tensor (output) lp # plastic velocity gradient tensor
(output) elasmatrix # elastic stiffness matrix (output) elasmatrix # elastic stiffness matrix

2
examples/SpectralMethod/Polycrystal/material.config
View File

@ -18,8 +18,6 @@ mech none
(output) f # deformation gradient tensor; synonyms: "defgrad" (output) f # deformation gradient tensor; synonyms: "defgrad"
(output) fe # elastic deformation gradient tensor (output) fe # elastic deformation gradient tensor
(output) fp # plastic deformation gradient tensor (output) fp # plastic deformation gradient tensor
(output) e # total strain as Green-Lagrange tensor
(output) ee # elastic strain as Green-Lagrange tensor
(output) p # first Piola-Kichhoff stress tensor; synonyms: "firstpiola", "1stpiola" (output) p # first Piola-Kichhoff stress tensor; synonyms: "firstpiola", "1stpiola"
(output) lp # plastic velocity gradient tensor (output) lp # plastic velocity gradient tensor

22
src/DAMASK_FEM.f90
View File

@ -42,9 +42,9 @@ program DAMASK_FEM
restartWrite, & restartWrite, &
restartInc restartInc
use numerics, only: & use numerics, only: &
worldrank, &
maxCutBack, & maxCutBack, &
stagItMax, & stagItMax
worldrank
use mesh, only: & use mesh, only: &
mesh_Nboundaries, & mesh_Nboundaries, &
mesh_boundaries, & mesh_boundaries, &
@ -73,11 +73,7 @@ program DAMASK_FEM
COMPONENT_SOLUTE_CVaH_ID, & COMPONENT_SOLUTE_CVaH_ID, &
COMPONENT_SOLUTE_CVaHPOT_ID, & COMPONENT_SOLUTE_CVaHPOT_ID, &
COMPONENT_MGTWIN_PHI_ID, & COMPONENT_MGTWIN_PHI_ID, &
FIELD_MECH_label, & FIELD_MECH_label
FIELD_THERMAL_label, &
FIELD_DAMAGE_label, &
FIELD_SOLUTE_label, &
FIELD_MGTWIN_label
use FEM_mech use FEM_mech
implicit none implicit none
@ -405,18 +401,6 @@ program DAMASK_FEM
case(FIELD_MECH_ID) case(FIELD_MECH_ID)
write(6,'(2x,a)') 'Field '//trim(FIELD_MECH_label) write(6,'(2x,a)') 'Field '//trim(FIELD_MECH_label)
case(FIELD_THERMAL_ID)
write(6,'(2x,a)') 'Field '//trim(FIELD_THERMAL_label)
case(FIELD_DAMAGE_ID)
write(6,'(2x,a)') 'Field '//trim(FIELD_DAMAGE_label)
case(FIELD_MGTWIN_ID)
write(6,'(2x,a)') 'Field '//trim(FIELD_MGTWIN_label)
case(FIELD_SOLUTE_ID)
write(6,'(2x,a)') 'Field '//trim(FIELD_SOLUTE_label)
end select end select
do faceSet = 1_pInt, mesh_Nboundaries do faceSet = 1_pInt, mesh_Nboundaries
do component = 1_pInt, loadCases(currentLoadCase)%fieldBC(field)%nComponents do component = 1_pInt, loadCases(currentLoadCase)%fieldBC(field)%nComponents

22
src/FEM_mech.f90
View File

@ -5,7 +5,6 @@
!> @brief FEM PETSc solver !> @brief FEM PETSc solver
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module FEM_mech module FEM_mech
use PETSC
#include <petsc/finclude/petscdmplex.h> #include <petsc/finclude/petscdmplex.h>
#include <petsc/finclude/petscdm.h> #include <petsc/finclude/petscdm.h>
#include <petsc/finclude/petscdmda.h> #include <petsc/finclude/petscdmda.h>
@ -15,7 +14,7 @@ use PETSC
use PETScsnes use PETScsnes
use PETScDM use PETScDM
use PETScDMplex use PETScDMplex
use PETSC use PETScDT
use prec, only: & use prec, only: &
pInt, & pInt, &
pReal pReal
@ -25,9 +24,6 @@ use PETSC
tSolutionState, & tSolutionState, &
tFieldBC, & tFieldBC, &
tComponentBC tComponentBC
use numerics, only: &
worldrank, &
worldsize
use mesh, only: & use mesh, only: &
mesh_Nboundaries, & mesh_Nboundaries, &
mesh_boundaries mesh_boundaries
@ -66,7 +62,7 @@ use PETSC
FEM_mech_solution ,& FEM_mech_solution ,&
FEM_mech_forward, & FEM_mech_forward, &
FEM_mech_destroy FEM_mech_destroy
external :: PETScerrorf
contains contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -82,7 +78,6 @@ subroutine FEM_mech_init(fieldBC)
use mesh, only: & use mesh, only: &
geomMesh geomMesh
use numerics, only: & use numerics, only: &
worldrank, &
itmax, & itmax, &
integrationOrder integrationOrder
use FEM_Zoo, only: & use FEM_Zoo, only: &
@ -106,8 +101,8 @@ subroutine FEM_mech_init(fieldBC)
IS, pointer :: pBcComps(:), pBcPoints(:) IS, pointer :: pBcComps(:), pBcPoints(:)
PetscSection :: section PetscSection :: section
PetscInt :: field, faceSet, topologDim, nNodalPoints PetscInt :: field, faceSet, topologDim, nNodalPoints
PetscReal, pointer :: qPointsP(:), qWeightsP(:), & PetscReal, dimension(:) , pointer :: qPointsP, qWeightsP, &
nodalPointsP(:), nodalWeightsP(:) nodalPointsP, nodalWeightsP
PetscReal, allocatable, target :: nodalPoints(:), nodalWeights(:) PetscReal, allocatable, target :: nodalPoints(:), nodalWeights(:)
PetscScalar, pointer :: px_scal(:) PetscScalar, pointer :: px_scal(:)
PetscScalar, allocatable, target :: x_scal(:) PetscScalar, allocatable, target :: x_scal(:)
@ -117,6 +112,7 @@ subroutine FEM_mech_init(fieldBC)
PetscInt :: cellStart, cellEnd, cell, basis PetscInt :: cellStart, cellEnd, cell, basis
character(len=7) :: prefix = 'mechFE_' character(len=7) :: prefix = 'mechFE_'
PetscErrorCode :: ierr PetscErrorCode :: ierr
PetscReal, allocatable, target, dimension(:) :: qWeights
write(6,'(/,a)') ' <<<+- FEM_mech init -+>>>' write(6,'(/,a)') ' <<<+- FEM_mech init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp() write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
@ -129,8 +125,6 @@ subroutine FEM_mech_init(fieldBC)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! Setup FEM mech discretization ! Setup FEM mech discretization
allocate(qPoints(dimPlex*FEM_Zoo_nQuadrature(dimPlex,integrationOrder)))
allocate(qWeights(FEM_Zoo_nQuadrature(dimPlex,integrationOrder)))
qPoints = FEM_Zoo_QuadraturePoints(dimPlex,integrationOrder)%p qPoints = FEM_Zoo_QuadraturePoints(dimPlex,integrationOrder)%p
qWeights = FEM_Zoo_QuadratureWeights(dimPlex,integrationOrder)%p qWeights = FEM_Zoo_QuadratureWeights(dimPlex,integrationOrder)%p
nQuadrature = FEM_Zoo_nQuadrature(dimPlex,integrationOrder) nQuadrature = FEM_Zoo_nQuadrature(dimPlex,integrationOrder)
@ -250,12 +244,17 @@ subroutine FEM_mech_init(fieldBC)
call PetscFEGetDualSpace(mechFE,mechDualSpace,ierr); CHKERRQ(ierr) call PetscFEGetDualSpace(mechFE,mechDualSpace,ierr); CHKERRQ(ierr)
call DMPlexGetHeightStratum(mech_mesh,0,cellStart,cellEnd,ierr) call DMPlexGetHeightStratum(mech_mesh,0,cellStart,cellEnd,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
write(6,*) 'cellDof', cellDof;flush(6)
write(6,*) 'cell start and end-1',cellStart,cellEnd-1;flush(6)
do cell = cellStart, cellEnd-1 !< loop over all elements do cell = cellStart, cellEnd-1 !< loop over all elements
write(6,*) 'cell',cell;flush(6)
x_scal = 0.0 x_scal = 0.0
call DMPlexComputeCellGeometryAffineFEM(mech_mesh,cell,pV0,pCellJ,pInvcellJ,detJ,ierr) call DMPlexComputeCellGeometryAffineFEM(mech_mesh,cell,pV0,pCellJ,pInvcellJ,detJ,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
cellJMat = reshape(pCellJ,shape=[dimPlex,dimPlex]) cellJMat = reshape(pCellJ,shape=[dimPlex,dimPlex])
do basis = 0, nBasis-1 do basis = 0, nBasis-1
write(6,*) 'nBasis-1',nBasis-1;flush(6)
write(6,*) 'basis',basis;flush(6)
call PetscDualSpaceGetFunctional(mechDualSpace,basis,functional,ierr) call PetscDualSpaceGetFunctional(mechDualSpace,basis,functional,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call PetscQuadratureGetData(functional,dimPlex,nc,nNodalPoints,nodalPointsP,nodalWeightsP,ierr) call PetscQuadratureGetData(functional,dimPlex,nc,nNodalPoints,nodalPointsP,nodalWeightsP,ierr)
@ -677,6 +676,7 @@ end subroutine FEM_mech_forward
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine FEM_mech_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr) subroutine FEM_mech_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr)
use numerics, only: & use numerics, only: &
worldrank, &
err_struct_tolAbs, & err_struct_tolAbs, &
err_struct_tolRel err_struct_tolRel
use IO, only: & use IO, only: &

42
src/FEM_utilities.f90
View File

@ -23,20 +23,16 @@ use PETScis
! grid related information information ! grid related information information
real(pReal), public :: wgt !< weighting factor 1/Nelems real(pReal), public :: wgt !< weighting factor 1/Nelems
real(pReal), public :: wgtDof !< weighting factor 1/Nelems real(pReal), public :: wgtDof !< weighting factor 1/Nelems
real(pReal), public :: C_volAvg(3,3,3,3)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! output data ! output data
Vec, public :: coordinatesVec Vec, public :: coordinatesVec
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! field labels information ! field labels information
character(len=*), parameter, public :: & character(len=*), parameter, public :: &
FIELD_MECH_label = 'mechanical', & FIELD_MECH_label = 'mechanical'
FIELD_THERMAL_label = 'thermal', &
FIELD_DAMAGE_label = 'damage', & integer(pInt), parameter :: structOrder = 2_pInt
FIELD_SOLUTE_label = 'solute', &
FIELD_MGTWIN_label = 'mgtwin'
enum, bind(c) enum, bind(c)
enumerator :: FIELD_UNDEFINED_ID, & enumerator :: FIELD_UNDEFINED_ID, &
@ -121,23 +117,12 @@ use PETScis
utilities_indexActiveSet, & utilities_indexActiveSet, &
utilities_destroy, & utilities_destroy, &
FIELD_MECH_ID, & FIELD_MECH_ID, &
FIELD_THERMAL_ID, &
FIELD_DAMAGE_ID, &
FIELD_SOLUTE_ID, &
FIELD_MGTWIN_ID, &
COMPONENT_MECH_X_ID, & COMPONENT_MECH_X_ID, &
COMPONENT_MECH_Y_ID, & COMPONENT_MECH_Y_ID, &
COMPONENT_MECH_Z_ID, & COMPONENT_MECH_Z_ID, &
COMPONENT_THERMAL_T_ID, & COMPONENT_THERMAL_T_ID
COMPONENT_DAMAGE_PHI_ID, &
COMPONENT_SOLUTE_CV_ID, &
COMPONENT_SOLUTE_CVPOT_ID, &
COMPONENT_SOLUTE_CH_ID, &
COMPONENT_SOLUTE_CHPOT_ID, &
COMPONENT_SOLUTE_CVaH_ID, &
COMPONENT_SOLUTE_CVaHPOT_ID, &
COMPONENT_MGTWIN_PHI_ID
external :: PETScErrorF
contains contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -176,14 +161,11 @@ subroutine utilities_init()
implicit none implicit none
character(len=1024) :: petsc_optionsPhysics, grainStr character(len=1024) :: petsc_optionsPhysics
integer(pInt) :: dimPlex integer(pInt) :: dimPlex
integer(pInt) :: headerID = 205_pInt integer(pInt) :: headerID = 205_pInt
PetscInt, dimension(:), pointer :: points PetscInt, allocatable :: nEntities(:), nOutputCells(:), nOutputNodes(:)
PetscInt, allocatable :: nEntities(:), nOutputCells(:), nOutputNodes(:), mappingCells(:) PetscInt :: dim
PetscInt :: cellStart, cellEnd, cell, ip, dim, ctr, qPt
PetscInt, allocatable :: connectivity(:,:)
Vec :: connectivityVec
PetscErrorCode :: ierr PetscErrorCode :: ierr
write(6,'(/,a)') ' <<<+- DAMASK_FEM_utilities init -+>>>' write(6,'(/,a)') ' <<<+- DAMASK_FEM_utilities init -+>>>'
@ -207,7 +189,7 @@ subroutine utilities_init()
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_defaultOptions),ierr) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_defaultOptions),ierr)
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_options),ierr) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_options),ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
!write(petsc_optionsPhysics,'(a,i0)') '-mechFE_petscspace_order ' , structOrder write(petsc_optionsPhysics,'(a,i0)') '-mechFE_petscspace_order ' , structOrder
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_optionsPhysics),ierr) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_optionsPhysics),ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
@ -249,8 +231,6 @@ subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData)
use debug, only: & use debug, only: &
debug_reset, & debug_reset, &
debug_info debug_info
use numerics, only: &
worldrank
use math, only: & use math, only: &
math_transpose33, & math_transpose33, &
math_rotate_forward33, & math_rotate_forward33, &
@ -258,10 +238,8 @@ subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData)
use FEsolving, only: & use FEsolving, only: &
restartWrite restartWrite
use homogenization, only: & use homogenization, only: &
materialpoint_F0, &
materialpoint_F, & materialpoint_F, &
materialpoint_P, & materialpoint_P, &
materialpoint_dPdF, &
materialpoint_stressAndItsTangent materialpoint_stressAndItsTangent
use mesh, only: & use mesh, only: &
mesh_NcpElems mesh_NcpElems
@ -314,9 +292,7 @@ subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData)
cutBack = .false. ! reset cutBack status cutBack = .false. ! reset cutBack status
P_av = sum(sum(materialpoint_P,dim=4),dim=3) * wgt ! average of P P_av = sum(sum(materialpoint_P,dim=4),dim=3) * wgt ! average of P
C_volAvg = sum(sum(materialpoint_dPdF,dim=6),dim=5) * wgt
call MPI_Allreduce(MPI_IN_PLACE,P_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,P_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
call MPI_Allreduce(MPI_IN_PLACE,C_volAvg,81,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD, ierr)
end subroutine utilities_constitutiveResponse end subroutine utilities_constitutiveResponse

247
src/constitutive.f90
View File

@ -19,16 +19,16 @@ module constitutive
constitutive_init, & constitutive_init, &
constitutive_homogenizedC, & constitutive_homogenizedC, &
constitutive_microstructure, & constitutive_microstructure, &
constitutive_LpAndItsTangent, & constitutive_LpAndItsTangents, &
constitutive_LiAndItsTangent, & constitutive_LiAndItsTangents, &
constitutive_initialFi, & constitutive_initialFi, &
constitutive_TandItsTangent, & constitutive_SandItsTangents, &
constitutive_collectDotState, & constitutive_collectDotState, &
constitutive_collectDeltaState, & constitutive_collectDeltaState, &
constitutive_postResults constitutive_postResults
private :: & private :: &
constitutive_hooke_TandItsTangent constitutive_hooke_SandItsTangents
contains contains
@ -346,6 +346,7 @@ end subroutine constitutive_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief returns the homogenize elasticity matrix !> @brief returns the homogenize elasticity matrix
!> ToDo: homogenizedC66 would be more consistent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function constitutive_homogenizedC(ipc,ip,el) function constitutive_homogenizedC(ipc,ip,el)
use prec, only: & use prec, only: &
@ -430,7 +431,7 @@ end subroutine constitutive_microstructure
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the velocity gradient !> @brief contains the constitutive equation for calculating the velocity gradient
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar3333, dLp_dFi3333, Tstar_v, Fi, ipc, ip, el) subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, S6, Fi, ipc, ip, el)
use prec, only: & use prec, only: &
pReal pReal
use math, only: & use math, only: &
@ -470,20 +471,21 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar3333, dLp_dFi3333, Tstar_v
ip, & !< integration point ip, & !< integration point
el !< element el !< element
real(pReal), intent(in), dimension(6) :: & real(pReal), intent(in), dimension(6) :: &
Tstar_v !< 2nd Piola-Kirchhoff stress S6 !< 2nd Piola-Kirchhoff stress (vector notation)
real(pReal), intent(in), dimension(3,3) :: & real(pReal), intent(in), dimension(3,3) :: &
Fi !< intermediate deformation gradient Fi !< intermediate deformation gradient
real(pReal), intent(out), dimension(3,3) :: & real(pReal), intent(out), dimension(3,3) :: &
Lp !< plastic velocity gradient Lp !< plastic velocity gradient
real(pReal), intent(out), dimension(3,3,3,3) :: & real(pReal), intent(out), dimension(3,3,3,3) :: &
dLp_dTstar3333, & !< derivative of Lp with respect to Tstar (4th-order tensor) dLp_dS, &
dLp_dFi3333 !< derivative of Lp with respect to Fi (4th-order tensor) dLp_dFi !< derivative of Lp with respect to Fi
real(pReal), dimension(6) :: & real(pReal), dimension(3,3,3,3) :: &
Mstar_v !< Mandel stress work conjugate with Lp dLp_dMp !< derivative of Lp with respect to Mandel stress
real(pReal), dimension(9,9) :: & real(pReal), dimension(9,9) :: &
dLp_dMstar !< derivative of Lp with respect to Mstar (4th-order tensor) dLp_dMp99 !< derivative of Lp with respect to Mstar (matrix notation)
real(pReal), dimension(3,3) :: & real(pReal), dimension(3,3) :: &
temp_33 Mp, & !< Mandel stress work conjugate with Lp
S !< 2nd Piola-Kirchhoff stress
integer(pInt) :: & integer(pInt) :: &
ho, & !< homogenization ho, & !< homogenization
tme !< thermal member position tme !< thermal member position
@ -493,47 +495,57 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar3333, dLp_dFi3333, Tstar_v
ho = material_homog(ip,el) ho = material_homog(ip,el)
tme = thermalMapping(ho)%p(ip,el) tme = thermalMapping(ho)%p(ip,el)
Mstar_v = math_Mandel33to6(math_mul33x33(math_mul33x33(transpose(Fi),Fi),math_Mandel6to33(Tstar_v))) S = math_Mandel6to33(S6)
Mp = math_mul33x33(math_mul33x33(transpose(Fi),Fi),S)
plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el))) plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
case (PLASTICITY_NONE_ID) plasticityType case (PLASTICITY_NONE_ID) plasticityType
Lp = 0.0_pReal Lp = 0.0_pReal
dLp_dMstar = 0.0_pReal dLp_dMp = 0.0_pReal
case (PLASTICITY_ISOTROPIC_ID) plasticityType case (PLASTICITY_ISOTROPIC_ID) plasticityType
call plastic_isotropic_LpAndItsTangent (Lp,dLp_dMstar,Mstar_v,ipc,ip,el) call plastic_isotropic_LpAndItsTangent (Lp,dLp_dMp99, math_Mandel33to6(Mp),ipc,ip,el)
dLp_dMp = math_Plain99to3333(dLp_dMp99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget
case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType
call plastic_phenopowerlaw_LpAndItsTangent (Lp,dLp_dMstar,Mstar_v,ipc,ip,el) call plastic_phenopowerlaw_LpAndItsTangent (Lp,dLp_dMp99, math_Mandel33to6(Mp),ipc,ip,el)
dLp_dMp = math_Plain99to3333(dLp_dMp99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget
case (PLASTICITY_KINEHARDENING_ID) plasticityType case (PLASTICITY_KINEHARDENING_ID) plasticityType
call plastic_kinehardening_LpAndItsTangent (Lp,dLp_dMstar,Mstar_v,ipc,ip,el) call plastic_kinehardening_LpAndItsTangent (Lp,dLp_dMp99, math_Mandel33to6(Mp),ipc,ip,el)
dLp_dMp = math_Plain99to3333(dLp_dMp99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget
case (PLASTICITY_NONLOCAL_ID) plasticityType case (PLASTICITY_NONLOCAL_ID) plasticityType
call plastic_nonlocal_LpAndItsTangent (Lp,dLp_dMstar,Mstar_v, & call plastic_nonlocal_LpAndItsTangent (Lp,dLp_dMp99, math_Mandel33to6(Mp), &
temperature(ho)%p(tme),ip,el) temperature(ho)%p(tme),ip,el)
dLp_dMp = math_Plain99to3333(dLp_dMp99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget
case (PLASTICITY_DISLOTWIN_ID) plasticityType case (PLASTICITY_DISLOTWIN_ID) plasticityType
call plastic_dislotwin_LpAndItsTangent (Lp,dLp_dMstar,Mstar_v, & call plastic_dislotwin_LpAndItsTangent (Lp,dLp_dMp99, math_Mandel33to6(Mp), &
temperature(ho)%p(tme),ipc,ip,el) temperature(ho)%p(tme),ipc,ip,el)
dLp_dMp = math_Plain99to3333(dLp_dMp99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget
case (PLASTICITY_DISLOUCLA_ID) plasticityType case (PLASTICITY_DISLOUCLA_ID) plasticityType
call plastic_disloucla_LpAndItsTangent (Lp,dLp_dMstar,Mstar_v, & call plastic_disloucla_LpAndItsTangent (Lp,dLp_dMp99, math_Mandel33to6(Mp), &
temperature(ho)%p(tme), ipc,ip,el) temperature(ho)%p(tme), ipc,ip,el)
dLp_dMp = math_Plain99to3333(dLp_dMp99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget
end select plasticityType end select plasticityType
dLp_dTstar3333 = math_Plain99to3333(dLp_dMstar) do concurrent(i = 1_pInt:3_pInt, j = 1_pInt:3_pInt)
temp_33 = math_mul33x33(Fi,math_Mandel6to33(Tstar_v)) dLp_dFi(i,j,1:3,1:3) = math_mul33x33(math_mul33x33(Fi,S),transpose(dLp_dMp(i,j,1:3,1:3))) + &
forall(i = 1_pInt:3_pInt, j = 1_pInt:3_pInt) & math_mul33x33(math_mul33x33(Fi,dLp_dMp(i,j,1:3,1:3)),S)
dLp_dFi3333(i,j,1:3,1:3) = math_mul33x33(temp_33,transpose(dLp_dTstar3333(i,j,1:3,1:3))) + & dLp_dS(i,j,1:3,1:3) = math_mul33x33(math_mul33x33(transpose(Fi),Fi),dLp_dMp(i,j,1:3,1:3)) ! ToDo: @PS: why not: dLp_dMp:(FiT Fi)
math_mul33x33(math_mul33x33(Fi,dLp_dTstar3333(i,j,1:3,1:3)),math_Mandel6to33(Tstar_v)) enddo
temp_33 = math_mul33x33(transpose(Fi),Fi) end subroutine constitutive_LpAndItsTangents
forall(i = 1_pInt:3_pInt, j = 1_pInt:3_pInt) &
dLp_dTstar3333(i,j,1:3,1:3) = math_mul33x33(temp_33,dLp_dTstar3333(i,j,1:3,1:3))
end subroutine constitutive_LpAndItsTangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the velocity gradient !> @brief contains the constitutive equation for calculating the velocity gradient
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_LiAndItsTangent(Li, dLi_dTstar3333, dLi_dFi3333, Tstar_v, Fi, ipc, ip, el) subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, S6, Fi, ipc, ip, el)
use prec, only: & use prec, only: &
pReal pReal
use math, only: & use math, only: &
@ -571,18 +583,18 @@ subroutine constitutive_LiAndItsTangent(Li, dLi_dTstar3333, dLi_dFi3333, Tstar_v
ip, & !< integration point ip, & !< integration point
el !< element el !< element
real(pReal), intent(in), dimension(6) :: & real(pReal), intent(in), dimension(6) :: &
Tstar_v !< 2nd Piola-Kirchhoff stress S6 !< 2nd Piola-Kirchhoff stress (vector notation)
real(pReal), intent(in), dimension(3,3) :: & real(pReal), intent(in), dimension(3,3) :: &
Fi !< intermediate deformation gradient Fi !< intermediate deformation gradient
real(pReal), intent(out), dimension(3,3) :: & real(pReal), intent(out), dimension(3,3) :: &
Li !< intermediate velocity gradient Li !< intermediate velocity gradient
real(pReal), intent(out), dimension(3,3,3,3) :: & real(pReal), intent(out), dimension(3,3,3,3) :: &
dLi_dTstar3333, & !< derivative of Li with respect to Tstar (4th-order tensor) dLi_dS, & !< derivative of Li with respect to S
dLi_dFi3333 dLi_dFi
real(pReal), dimension(3,3) :: & real(pReal), dimension(3,3) :: &
my_Li !< intermediate velocity gradient my_Li !< intermediate velocity gradient
real(pReal), dimension(3,3,3,3) :: & real(pReal), dimension(3,3,3,3) :: &
my_dLi_dTstar my_dLi_dS
real(pReal), dimension(3,3) :: & real(pReal), dimension(3,3) :: &
FiInv, & FiInv, &
temp_33 temp_33
@ -594,51 +606,51 @@ subroutine constitutive_LiAndItsTangent(Li, dLi_dTstar3333, dLi_dFi3333, Tstar_v
i, j i, j
Li = 0.0_pReal Li = 0.0_pReal
dLi_dTstar3333 = 0.0_pReal dLi_dS = 0.0_pReal
dLi_dFi3333 = 0.0_pReal dLi_dFi = 0.0_pReal
plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el))) plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
case (PLASTICITY_isotropic_ID) plasticityType case (PLASTICITY_isotropic_ID) plasticityType
call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dTstar, Tstar_v, ipc, ip, el) call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dS, S6, ipc, ip, el)
case default plasticityType case default plasticityType
my_Li = 0.0_pReal my_Li = 0.0_pReal
my_dLi_dTstar = 0.0_pReal my_dLi_dS = 0.0_pReal
end select plasticityType end select plasticityType
Li = Li + my_Li Li = Li + my_Li
dLi_dTstar3333 = dLi_dTstar3333 + my_dLi_dTstar dLi_dS = dLi_dS + my_dLi_dS
KinematicsLoop: do k = 1_pInt, phase_Nkinematics(material_phase(ipc,ip,el)) KinematicsLoop: do k = 1_pInt, phase_Nkinematics(material_phase(ipc,ip,el))
kinematicsType: select case (phase_kinematics(k,material_phase(ipc,ip,el))) kinematicsType: select case (phase_kinematics(k,material_phase(ipc,ip,el)))
case (KINEMATICS_cleavage_opening_ID) kinematicsType case (KINEMATICS_cleavage_opening_ID) kinematicsType
call kinematics_cleavage_opening_LiAndItsTangent(my_Li, my_dLi_dTstar, Tstar_v, ipc, ip, el) call kinematics_cleavage_opening_LiAndItsTangent(my_Li, my_dLi_dS, S6, ipc, ip, el)
case (KINEMATICS_slipplane_opening_ID) kinematicsType case (KINEMATICS_slipplane_opening_ID) kinematicsType
call kinematics_slipplane_opening_LiAndItsTangent(my_Li, my_dLi_dTstar, Tstar_v, ipc, ip, el) call kinematics_slipplane_opening_LiAndItsTangent(my_Li, my_dLi_dS, S6, ipc, ip, el)
case (KINEMATICS_thermal_expansion_ID) kinematicsType case (KINEMATICS_thermal_expansion_ID) kinematicsType
call kinematics_thermal_expansion_LiAndItsTangent(my_Li, my_dLi_dTstar, ipc, ip, el) call kinematics_thermal_expansion_LiAndItsTangent(my_Li, my_dLi_dS, ipc, ip, el)
case (KINEMATICS_vacancy_strain_ID) kinematicsType case (KINEMATICS_vacancy_strain_ID) kinematicsType
call kinematics_vacancy_strain_LiAndItsTangent(my_Li, my_dLi_dTstar, ipc, ip, el) call kinematics_vacancy_strain_LiAndItsTangent(my_Li, my_dLi_dS, ipc, ip, el)
case (KINEMATICS_hydrogen_strain_ID) kinematicsType case (KINEMATICS_hydrogen_strain_ID) kinematicsType
call kinematics_hydrogen_strain_LiAndItsTangent(my_Li, my_dLi_dTstar, ipc, ip, el) call kinematics_hydrogen_strain_LiAndItsTangent(my_Li, my_dLi_dS, ipc, ip, el)
case default kinematicsType case default kinematicsType
my_Li = 0.0_pReal my_Li = 0.0_pReal
my_dLi_dTstar = 0.0_pReal my_dLi_dS = 0.0_pReal
end select kinematicsType end select kinematicsType
Li = Li + my_Li Li = Li + my_Li
dLi_dTstar3333 = dLi_dTstar3333 + my_dLi_dTstar dLi_dS = dLi_dS + my_dLi_dS
enddo KinematicsLoop enddo KinematicsLoop
FiInv = math_inv33(Fi) FiInv = math_inv33(Fi)
detFi = math_det33(Fi) detFi = math_det33(Fi)
Li = math_mul33x33(math_mul33x33(Fi,Li),FiInv)*detFi !< push forward to intermediate configuration Li = math_mul33x33(math_mul33x33(Fi,Li),FiInv)*detFi !< push forward to intermediate configuration
temp_33 = math_mul33x33(FiInv,Li) temp_33 = math_mul33x33(FiInv,Li)
forall(i = 1_pInt:3_pInt, j = 1_pInt:3_pInt) do concurrent(i = 1_pInt:3_pInt, j = 1_pInt:3_pInt)
dLi_dTstar3333(1:3,1:3,i,j) = math_mul33x33(math_mul33x33(Fi,dLi_dTstar3333(1:3,1:3,i,j)),FiInv)*detFi dLi_dS(1:3,1:3,i,j) = math_mul33x33(math_mul33x33(Fi,dLi_dS(1:3,1:3,i,j)),FiInv)*detFi
dLi_dFi3333 (1:3,1:3,i,j) = dLi_dFi3333(1:3,1:3,i,j) + Li*FiInv(j,i) dLi_dFi(1:3,1:3,i,j) = dLi_dFi(1:3,1:3,i,j) + Li*FiInv(j,i)
dLi_dFi3333 (1:3,i,1:3,j) = dLi_dFi3333(1:3,i,1:3,j) + math_I3*temp_33(j,i) + Li*FiInv(j,i) dLi_dFi(1:3,i,1:3,j) = dLi_dFi(1:3,i,1:3,j) + math_I3*temp_33(j,i) + Li*FiInv(j,i)
end forall end do
end subroutine constitutive_LiAndItsTangent end subroutine constitutive_LiAndItsTangents
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -696,10 +708,10 @@ end function constitutive_initialFi
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to !> @brief returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to
!> the elastic deformation gradient depending on the selected elastic law (so far no case switch !> the elastic/intermediate deformation gradients depending on the selected elastic law
!! because only Hooke is implemented !! (so far no case switch because only Hooke is implemented)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_TandItsTangent(T, dT_dFe, dT_dFi, Fe, Fi, ipc, ip, el) subroutine constitutive_SandItsTangents(S, dS_dFe, dS_dFi, Fe, Fi, ipc, ip, el)
use prec, only: & use prec, only: &
pReal pReal
@ -712,30 +724,28 @@ subroutine constitutive_TandItsTangent(T, dT_dFe, dT_dFi, Fe, Fi, ipc, ip, el)
Fe, & !< elastic deformation gradient Fe, & !< elastic deformation gradient
Fi !< intermediate deformation gradient Fi !< intermediate deformation gradient
real(pReal), intent(out), dimension(3,3) :: & real(pReal), intent(out), dimension(3,3) :: &
T !< 2nd Piola-Kirchhoff stress tensor S !< 2nd Piola-Kirchhoff stress tensor
real(pReal), intent(out), dimension(3,3,3,3) :: & real(pReal), intent(out), dimension(3,3,3,3) :: &
dT_dFe, & !< derivative of 2nd P-K stress with respect to elastic deformation gradient dS_dFe, & !< derivative of 2nd P-K stress with respect to elastic deformation gradient
dT_dFi !< derivative of 2nd P-K stress with respect to intermediate deformation gradient dS_dFi !< derivative of 2nd P-K stress with respect to intermediate deformation gradient
call constitutive_hooke_TandItsTangent(T, dT_dFe, dT_dFi, Fe, Fi, ipc, ip, el) call constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, Fe, Fi, ipc, ip, el)
end subroutine constitutive_TandItsTangent end subroutine constitutive_SandItsTangents
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to !> @brief returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to
!> the elastic deformation gradient using hookes law !> the elastic and intermeidate deformation gradients using Hookes law
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_hooke_TandItsTangent(T, dT_dFe, dT_dFi, Fe, Fi, ipc, ip, el) subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, Fe, Fi, ipc, ip, el)
use prec, only: & use prec, only: &
pReal pReal
use math, only : & use math, only : &
math_mul3x3, &
math_mul33x33, & math_mul33x33, &
math_mul3333xx33, & math_mul3333xx33, &
math_Mandel66to3333, & math_Mandel66to3333, &
math_trace33, &
math_I3 math_I3
use material, only: & use material, only: &
material_phase, & material_phase, &
@ -758,10 +768,10 @@ subroutine constitutive_hooke_TandItsTangent(T, dT_dFe, dT_dFi, Fe, Fi, ipc, ip,
Fe, & !< elastic deformation gradient Fe, & !< elastic deformation gradient
Fi !< intermediate deformation gradient Fi !< intermediate deformation gradient
real(pReal), intent(out), dimension(3,3) :: & real(pReal), intent(out), dimension(3,3) :: &
T !< 2nd Piola-Kirchhoff stress tensor in lattice configuration S !< 2nd Piola-Kirchhoff stress tensor in lattice configuration
real(pReal), intent(out), dimension(3,3,3,3) :: & real(pReal), intent(out), dimension(3,3,3,3) :: &
dT_dFe, & !< derivative of 2nd P-K stress with respect to elastic deformation gradient dS_dFe, & !< derivative of 2nd P-K stress with respect to elastic deformation gradient
dT_dFi !< derivative of 2nd P-K stress with respect to intermediate deformation gradient dS_dFi !< derivative of 2nd P-K stress with respect to intermediate deformation gradient
real(pReal), dimension(3,3) :: E real(pReal), dimension(3,3) :: E
real(pReal), dimension(3,3,3,3) :: C real(pReal), dimension(3,3,3,3) :: C
integer(pInt) :: & integer(pInt) :: &
@ -784,22 +794,22 @@ subroutine constitutive_hooke_TandItsTangent(T, dT_dFe, dT_dFi, Fe, Fi, ipc, ip,
enddo DegradationLoop enddo DegradationLoop
E = 0.5_pReal*(math_mul33x33(transpose(Fe),Fe)-math_I3) !< Green-Lagrange strain in unloaded configuration E = 0.5_pReal*(math_mul33x33(transpose(Fe),Fe)-math_I3) !< Green-Lagrange strain in unloaded configuration
T = math_mul3333xx33(C,math_mul33x33(math_mul33x33(transpose(Fi),E),Fi)) !< 2PK stress in lattice configuration in work conjugate with GL strain pulled back to lattice configuration S = math_mul3333xx33(C,math_mul33x33(math_mul33x33(transpose(Fi),E),Fi)) !< 2PK stress in lattice configuration in work conjugate with GL strain pulled back to lattice configuration
dT_dFe = 0.0_pReal dS_dFe = 0.0_pReal
forall (i=1_pInt:3_pInt, j=1_pInt:3_pInt) forall (i=1_pInt:3_pInt, j=1_pInt:3_pInt)
dT_dFe(i,j,1:3,1:3) = & dS_dFe(i,j,1:3,1:3) = &
math_mul33x33(Fe,math_mul33x33(math_mul33x33(Fi,C(i,j,1:3,1:3)),transpose(Fi))) !< dT_ij/dFe_kl = C_ijmn * Fi_lm * Fi_on * Fe_ko math_mul33x33(Fe,math_mul33x33(math_mul33x33(Fi,C(i,j,1:3,1:3)),transpose(Fi))) !< dS_ij/dFe_kl = C_ijmn * Fi_lm * Fi_on * Fe_ko
dT_dFi(i,j,1:3,1:3) = 2.0_pReal*math_mul33x33(math_mul33x33(E,Fi),C(i,j,1:3,1:3)) !< dT_ij/dFi_kl = C_ijln * E_km * Fe_mn dS_dFi(i,j,1:3,1:3) = 2.0_pReal*math_mul33x33(math_mul33x33(E,Fi),C(i,j,1:3,1:3)) !< dS_ij/dFi_kl = C_ijln * E_km * Fe_mn
end forall end forall
end subroutine constitutive_hooke_TandItsTangent end subroutine constitutive_hooke_SandItsTangents
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the rate of change of microstructure !> @brief contains the constitutive equation for calculating the rate of change of microstructure
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, subdt, subfracArray,ipc, ip, el) subroutine constitutive_collectDotState(S6, FeArray, Fi, FpArray, subdt, subfracArray,ipc, ip, el)
use prec, only: & use prec, only: &
pReal, & pReal, &
pLongInt pLongInt
@ -807,6 +817,10 @@ subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, subdt, subfra
debug_level, & debug_level, &
debug_constitutive, & debug_constitutive, &
debug_levelBasic debug_levelBasic
use math, only: &
math_Mandel6to33, &
math_Mandel33to6, &
math_mul33x33
use mesh, only: & use mesh, only: &
mesh_NcpElems, & mesh_NcpElems, &
mesh_maxNips mesh_maxNips
@ -863,8 +877,12 @@ subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, subdt, subfra
real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
FeArray, & !< elastic deformation gradient FeArray, & !< elastic deformation gradient
FpArray !< plastic deformation gradient FpArray !< plastic deformation gradient
real(pReal), intent(in), dimension(3,3) :: &
Fi !< intermediate deformation gradient
real(pReal), intent(in), dimension(6) :: & real(pReal), intent(in), dimension(6) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor (Mandel) S6 !< 2nd Piola Kirchhoff stress (vector notation)
real(pReal), dimension(3,3) :: &
Mstar
integer(pInt) :: & integer(pInt) :: &
ho, & !< homogenization ho, & !< homogenization
tme, & !< thermal member position tme, & !< thermal member position
@ -873,35 +891,50 @@ subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, subdt, subfra
ho = material_homog( ip,el) ho = material_homog( ip,el)
tme = thermalMapping(ho)%p(ip,el) tme = thermalMapping(ho)%p(ip,el)
Mstar = math_mul33x33(math_mul33x33(transpose(Fi),Fi),math_Mandel6to33(S6))
plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el))) plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
case (PLASTICITY_ISOTROPIC_ID) plasticityType case (PLASTICITY_ISOTROPIC_ID) plasticityType
call plastic_isotropic_dotState (Tstar_v,ipc,ip,el) call plastic_isotropic_dotState (math_Mandel33to6(Mstar),ipc,ip,el)
case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType
call plastic_phenopowerlaw_dotState(Tstar_v,ipc,ip,el) call plastic_phenopowerlaw_dotState(math_Mandel33to6(Mstar),ipc,ip,el)
case (PLASTICITY_KINEHARDENING_ID) plasticityType case (PLASTICITY_KINEHARDENING_ID) plasticityType
call plastic_kinehardening_dotState(Tstar_v,ipc,ip,el) call plastic_kinehardening_dotState(math_Mandel33to6(Mstar),ipc,ip,el)
case (PLASTICITY_DISLOTWIN_ID) plasticityType case (PLASTICITY_DISLOTWIN_ID) plasticityType
call plastic_dislotwin_dotState (Tstar_v,temperature(ho)%p(tme), & call plastic_dislotwin_dotState (math_Mandel33to6(Mstar),temperature(ho)%p(tme), &
ipc,ip,el) ipc,ip,el)
case (PLASTICITY_DISLOUCLA_ID) plasticityType case (PLASTICITY_DISLOUCLA_ID) plasticityType
call plastic_disloucla_dotState (Tstar_v,temperature(ho)%p(tme), & call plastic_disloucla_dotState (math_Mandel33to6(Mstar),temperature(ho)%p(tme), &
ipc,ip,el) ipc,ip,el)
case (PLASTICITY_NONLOCAL_ID) plasticityType case (PLASTICITY_NONLOCAL_ID) plasticityType
call plastic_nonlocal_dotState (Tstar_v,FeArray,FpArray,temperature(ho)%p(tme), & call plastic_nonlocal_dotState (math_Mandel33to6(Mstar),FeArray,FpArray,temperature(ho)%p(tme), &
subdt,subfracArray,ip,el) subdt,subfracArray,ip,el)
end select plasticityType end select plasticityType
SourceLoop: do s = 1_pInt, phase_Nsources(material_phase(ipc,ip,el)) SourceLoop: do s = 1_pInt, phase_Nsources(material_phase(ipc,ip,el))
sourceType: select case (phase_source(s,material_phase(ipc,ip,el))) sourceType: select case (phase_source(s,material_phase(ipc,ip,el)))
case (SOURCE_damage_anisoBrittle_ID) sourceType case (SOURCE_damage_anisoBrittle_ID) sourceType
call source_damage_anisoBrittle_dotState (Tstar_v, ipc, ip, el) call source_damage_anisoBrittle_dotState (S6, ipc, ip, el) !< correct stress?
case (SOURCE_damage_isoDuctile_ID) sourceType case (SOURCE_damage_isoDuctile_ID) sourceType
call source_damage_isoDuctile_dotState ( ipc, ip, el) call source_damage_isoDuctile_dotState ( ipc, ip, el)
case (SOURCE_damage_anisoDuctile_ID) sourceType case (SOURCE_damage_anisoDuctile_ID) sourceType
call source_damage_anisoDuctile_dotState ( ipc, ip, el) call source_damage_anisoDuctile_dotState ( ipc, ip, el)
case (SOURCE_thermal_externalheat_ID) sourceType case (SOURCE_thermal_externalheat_ID) sourceType
call source_thermal_externalheat_dotState( ipc, ip, el) call source_thermal_externalheat_dotState( ipc, ip, el)
end select sourceType end select sourceType
enddo SourceLoop enddo SourceLoop
end subroutine constitutive_collectDotState end subroutine constitutive_collectDotState
@ -910,7 +943,7 @@ end subroutine constitutive_collectDotState
!> @brief for constitutive models having an instantaneous change of state !> @brief for constitutive models having an instantaneous change of state
!> will return false if delta state is not needed/supported by the constitutive model !> will return false if delta state is not needed/supported by the constitutive model
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_collectDeltaState(Tstar_v, Fe, ipc, ip, el) subroutine constitutive_collectDeltaState(S6, Fe, Fi, ipc, ip, el)
use prec, only: & use prec, only: &
pReal, & pReal, &
pLongInt pLongInt
@ -918,6 +951,10 @@ subroutine constitutive_collectDeltaState(Tstar_v, Fe, ipc, ip, el)
debug_level, & debug_level, &
debug_constitutive, & debug_constitutive, &
debug_levelBasic debug_levelBasic
use math, only: &
math_Mandel6to33, &
math_Mandel33to6, &
math_mul33x33
use material, only: & use material, only: &
phase_plasticity, & phase_plasticity, &
phase_source, & phase_source, &
@ -945,29 +982,43 @@ subroutine constitutive_collectDeltaState(Tstar_v, Fe, ipc, ip, el)
ip, & !< integration point ip, & !< integration point
el !< element el !< element
real(pReal), intent(in), dimension(6) :: & real(pReal), intent(in), dimension(6) :: &
Tstar_v !< 2nd Piola-Kirchhoff stress S6 !< 2nd Piola Kirchhoff stress (vector notation)
real(pReal), intent(in), dimension(3,3) :: & real(pReal), intent(in), dimension(3,3) :: &
Fe !< elastic deformation gradient Fe, & !< elastic deformation gradient
Fi !< intermediate deformation gradient
real(pReal), dimension(3,3) :: &
Mstar
integer(pInt) :: & integer(pInt) :: &
s !< counter in source loop s !< counter in source loop
Mstar = math_mul33x33(math_mul33x33(transpose(Fi),Fi),math_Mandel6to33(S6))
plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el))) plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
case (PLASTICITY_KINEHARDENING_ID) plasticityType case (PLASTICITY_KINEHARDENING_ID) plasticityType
call plastic_kinehardening_deltaState(Tstar_v,ipc,ip,el) call plastic_kinehardening_deltaState(math_Mandel33to6(Mstar),ipc,ip,el)
case (PLASTICITY_NONLOCAL_ID) plasticityType case (PLASTICITY_NONLOCAL_ID) plasticityType
call plastic_nonlocal_deltaState(Tstar_v,ip,el) call plastic_nonlocal_deltaState(math_Mandel33to6(Mstar),ip,el)
end select plasticityType end select plasticityType
SourceLoop: do s = 1_pInt, phase_Nsources(material_phase(ipc,ip,el)) sourceLoop: do s = 1_pInt, phase_Nsources(material_phase(ipc,ip,el))
sourceType: select case (phase_source(s,material_phase(ipc,ip,el))) sourceType: select case (phase_source(s,material_phase(ipc,ip,el)))
case (SOURCE_damage_isoBrittle_ID) sourceType case (SOURCE_damage_isoBrittle_ID) sourceType
call source_damage_isoBrittle_deltaState (constitutive_homogenizedC(ipc,ip,el), Fe, & call source_damage_isoBrittle_deltaState (constitutive_homogenizedC(ipc,ip,el), Fe, &
ipc, ip, el) ipc, ip, el)
case (SOURCE_vacancy_irradiation_ID) sourceType case (SOURCE_vacancy_irradiation_ID) sourceType
call source_vacancy_irradiation_deltaState(ipc, ip, el) call source_vacancy_irradiation_deltaState(ipc, ip, el)
case (SOURCE_vacancy_thermalfluc_ID) sourceType case (SOURCE_vacancy_thermalfluc_ID) sourceType
call source_vacancy_thermalfluc_deltaState(ipc, ip, el) call source_vacancy_thermalfluc_deltaState(ipc, ip, el)
end select sourceType end select sourceType
enddo SourceLoop enddo SourceLoop
end subroutine constitutive_collectDeltaState end subroutine constitutive_collectDeltaState
@ -976,7 +1027,7 @@ end subroutine constitutive_collectDeltaState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief returns array of constitutive results !> @brief returns array of constitutive results
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function constitutive_postResults(Tstar_v, FeArray, ipc, ip, el) function constitutive_postResults(S6, FeArray, ipc, ip, el)
use prec, only: & use prec, only: &
pReal pReal
use mesh, only: & use mesh, only: &
@ -1036,7 +1087,7 @@ function constitutive_postResults(Tstar_v, FeArray, ipc, ip, el)
real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
FeArray !< elastic deformation gradient FeArray !< elastic deformation gradient
real(pReal), intent(in), dimension(6) :: & real(pReal), intent(in), dimension(6) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor (Mandel) S6 !< 2nd Piola Kirchhoff stress (vector notation)
integer(pInt) :: & integer(pInt) :: &
startPos, endPos startPos, endPos
integer(pInt) :: & integer(pInt) :: &
@ -1054,22 +1105,22 @@ function constitutive_postResults(Tstar_v, FeArray, ipc, ip, el)
plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el))) plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
case (PLASTICITY_ISOTROPIC_ID) plasticityType case (PLASTICITY_ISOTROPIC_ID) plasticityType
constitutive_postResults(startPos:endPos) = plastic_isotropic_postResults(Tstar_v,ipc,ip,el) constitutive_postResults(startPos:endPos) = plastic_isotropic_postResults(S6,ipc,ip,el)
case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType
constitutive_postResults(startPos:endPos) = & constitutive_postResults(startPos:endPos) = &
plastic_phenopowerlaw_postResults(Tstar_v,ipc,ip,el) plastic_phenopowerlaw_postResults(S6,ipc,ip,el)
case (PLASTICITY_KINEHARDENING_ID) plasticityType case (PLASTICITY_KINEHARDENING_ID) plasticityType
constitutive_postResults(startPos:endPos) = & constitutive_postResults(startPos:endPos) = &
plastic_kinehardening_postResults(Tstar_v,ipc,ip,el) plastic_kinehardening_postResults(S6,ipc,ip,el)
case (PLASTICITY_DISLOTWIN_ID) plasticityType case (PLASTICITY_DISLOTWIN_ID) plasticityType
constitutive_postResults(startPos:endPos) = & constitutive_postResults(startPos:endPos) = &
plastic_dislotwin_postResults(Tstar_v,temperature(ho)%p(tme),ipc,ip,el) plastic_dislotwin_postResults(S6,temperature(ho)%p(tme),ipc,ip,el)
case (PLASTICITY_DISLOUCLA_ID) plasticityType case (PLASTICITY_DISLOUCLA_ID) plasticityType
constitutive_postResults(startPos:endPos) = & constitutive_postResults(startPos:endPos) = &
plastic_disloucla_postResults(Tstar_v,temperature(ho)%p(tme),ipc,ip,el) plastic_disloucla_postResults(S6,temperature(ho)%p(tme),ipc,ip,el)
case (PLASTICITY_NONLOCAL_ID) plasticityType case (PLASTICITY_NONLOCAL_ID) plasticityType
constitutive_postResults(startPos:endPos) = & constitutive_postResults(startPos:endPos) = &
plastic_nonlocal_postResults (Tstar_v,FeArray,ip,el) plastic_nonlocal_postResults (S6,FeArray,ip,el)
end select plasticityType end select plasticityType
SourceLoop: do s = 1_pInt, phase_Nsources(material_phase(ipc,ip,el)) SourceLoop: do s = 1_pInt, phase_Nsources(material_phase(ipc,ip,el))

478
src/crystallite.f90
View File

File diff suppressed because it is too large Load Diff

32
src/homogenization.f90
View File

@ -448,8 +448,6 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
subStepSizeHomog, & subStepSizeHomog, &
stepIncreaseHomog, & stepIncreaseHomog, &
nMPstate nMPstate
use math, only: &
math_transpose33
use FEsolving, only: & use FEsolving, only: &
FEsolving_execElem, & FEsolving_execElem, &
FEsolving_execIP, & FEsolving_execIP, &
@ -496,6 +494,7 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
crystallite_converged, & crystallite_converged, &
crystallite_stressAndItsTangent, & crystallite_stressAndItsTangent, &
crystallite_orientations crystallite_orientations
#ifdef DEBUG
use debug, only: & use debug, only: &
debug_level, & debug_level, &
debug_homogenization, & debug_homogenization, &
@ -504,6 +503,7 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
debug_levelSelective, & debug_levelSelective, &
debug_e, & debug_e, &
debug_i debug_i
#endif
implicit none implicit none
real(pReal), intent(in) :: dt !< time increment real(pReal), intent(in) :: dt !< time increment
@ -517,18 +517,16 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
mySource, & mySource, &
myNgrains myNgrains
!-------------------------------------------------------------------------------------------------- #ifdef DEBUG
! initialize to starting condition
if (iand(debug_level(debug_homogenization), debug_levelBasic) /= 0_pInt) then if (iand(debug_level(debug_homogenization), debug_levelBasic) /= 0_pInt) then
!$OMP CRITICAL (write2out)
write(6,'(/a,i5,1x,i2)') '<< HOMOG >> Material Point start at el ip ', debug_e, debug_i write(6,'(/a,i5,1x,i2)') '<< HOMOG >> Material Point start at el ip ', debug_e, debug_i
write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< HOMOG >> F0', & write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< HOMOG >> F0', &
math_transpose33(materialpoint_F0(1:3,1:3,debug_i,debug_e)) transpose(materialpoint_F0(1:3,1:3,debug_i,debug_e))
write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< HOMOG >> F', & write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< HOMOG >> F', &
math_transpose33(materialpoint_F(1:3,1:3,debug_i,debug_e)) transpose(materialpoint_F(1:3,1:3,debug_i,debug_e))
!$OMP END CRITICAL (write2out)
endif endif
#endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! initialize restoration points of ... ! initialize restoration points of ...
@ -608,10 +606,8 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
! calculate new subStep and new subFrac ! calculate new subStep and new subFrac
materialpoint_subFrac(i,e) = materialpoint_subFrac(i,e) + materialpoint_subStep(i,e) materialpoint_subFrac(i,e) = materialpoint_subFrac(i,e) + materialpoint_subStep(i,e)
!$OMP FLUSH(materialpoint_subFrac)
materialpoint_subStep(i,e) = min(1.0_pReal-materialpoint_subFrac(i,e), & materialpoint_subStep(i,e) = min(1.0_pReal-materialpoint_subFrac(i,e), &
stepIncreaseHomog*materialpoint_subStep(i,e)) ! introduce flexibility for step increase/acceleration stepIncreaseHomog*materialpoint_subStep(i,e)) ! introduce flexibility for step increase/acceleration
!$OMP FLUSH(materialpoint_subStep)
steppingNeeded: if (materialpoint_subStep(i,e) > subStepMinHomog) then steppingNeeded: if (materialpoint_subStep(i,e) > subStepMinHomog) then
@ -671,7 +667,6 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
hydrogenfluxState(mappingHomogenization(2,i,e))%subState0(:,mappingHomogenization(1,i,e)) = & hydrogenfluxState(mappingHomogenization(2,i,e))%subState0(:,mappingHomogenization(1,i,e)) = &
hydrogenfluxState(mappingHomogenization(2,i,e))%State( :,mappingHomogenization(1,i,e))! ...internal hydrogen transport state hydrogenfluxState(mappingHomogenization(2,i,e))%State( :,mappingHomogenization(1,i,e))! ...internal hydrogen transport state
materialpoint_subF0(1:3,1:3,i,e) = materialpoint_subF(1:3,1:3,i,e) ! ...def grad materialpoint_subF0(1:3,1:3,i,e) = materialpoint_subF(1:3,1:3,i,e) ! ...def grad
!$OMP FLUSH(materialpoint_subF0)
endif steppingNeeded endif steppingNeeded
else converged else converged
@ -689,7 +684,6 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
!$OMP END CRITICAL (setTerminallyIll) !$OMP END CRITICAL (setTerminallyIll)
else ! cutback makes sense else ! cutback makes sense
materialpoint_subStep(i,e) = subStepSizeHomog * materialpoint_subStep(i,e) ! crystallite had severe trouble, so do a significant cutback materialpoint_subStep(i,e) = subStepSizeHomog * materialpoint_subStep(i,e) ! crystallite had severe trouble, so do a significant cutback
!$OMP FLUSH(materialpoint_subStep)
#ifdef DEBUG #ifdef DEBUG
if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0_pInt & if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0_pInt &
@ -752,8 +746,9 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
if (materialpoint_subStep(i,e) > subStepMinHomog) then if (materialpoint_subStep(i,e) > subStepMinHomog) then
materialpoint_requested(i,e) = .true. materialpoint_requested(i,e) = .true.
materialpoint_subF(1:3,1:3,i,e) = materialpoint_subF0(1:3,1:3,i,e) + & materialpoint_subF(1:3,1:3,i,e) = materialpoint_subF0(1:3,1:3,i,e) &
materialpoint_subStep(i,e) * (materialpoint_F(1:3,1:3,i,e) - materialpoint_F0(1:3,1:3,i,e)) + materialpoint_subStep(i,e) * (materialpoint_F(1:3,1:3,i,e) &
- materialpoint_F0(1:3,1:3,i,e))
materialpoint_subdt(i,e) = materialpoint_subStep(i,e) * dt materialpoint_subdt(i,e) = materialpoint_subStep(i,e) * dt
materialpoint_doneAndHappy(1:2,i,e) = [.false.,.true.] materialpoint_doneAndHappy(1:2,i,e) = [.false.,.true.]
endif endif
@ -810,13 +805,6 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
materialpoint_doneAndHappy(1:2,i,e) = homogenization_updateState(i,e) materialpoint_doneAndHappy(1:2,i,e) = homogenization_updateState(i,e)
materialpoint_converged(i,e) = all(materialpoint_doneAndHappy(1:2,i,e)) ! converged if done and happy materialpoint_converged(i,e) = all(materialpoint_doneAndHappy(1:2,i,e)) ! converged if done and happy
endif endif
!$OMP FLUSH(materialpoint_converged)
if (materialpoint_converged(i,e)) then
if (iand(debug_level(debug_homogenization), debug_levelBasic) /= 0_pInt) then
!$OMP CRITICAL (distributionMPState)
!$OMP END CRITICAL (distributionMPState)
endif
endif
endif endif
enddo IpLooping3 enddo IpLooping3
enddo elementLooping3 enddo elementLooping3
@ -838,9 +826,7 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
enddo elementLooping4 enddo elementLooping4
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
else else
!$OMP CRITICAL (write2out)
write(6,'(/,a,/)') '<< HOMOG >> Material Point terminally ill' write(6,'(/,a,/)') '<< HOMOG >> Material Point terminally ill'
!$OMP END CRITICAL (write2out)
endif endif
end subroutine materialpoint_stressAndItsTangent end subroutine materialpoint_stressAndItsTangent

10
src/prec.f90
View File

@ -87,16 +87,6 @@ module prec
integer(pInt), pointer, dimension(:,:,:) :: p integer(pInt), pointer, dimension(:,:,:) :: p
end type end type
#ifdef FEM
type, public :: tOutputData
integer(pInt) :: &
sizeIpCells = 0_pInt , &
sizeResults = 0_pInt
real(pReal), allocatable, dimension(:,:) :: &
output !< output data
end type
#endif
public :: & public :: &
prec_init, & prec_init, &
dEq, & dEq, &