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

grouping variables for better readability

This commit is contained in:
Martin Diehl 2020-01-21 06:10:19 +01:00
parent db96ee0fc2
commit 843676fb10
3 changed files with 69 additions and 58 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

107
src/mesh/FEM_mech.f90
View File

@ -67,28 +67,32 @@ contains
subroutine FEM_mech_init(fieldBC) subroutine FEM_mech_init(fieldBC)
type(tFieldBC), intent(in) :: fieldBC type(tFieldBC), intent(in) :: fieldBC
DM :: mech_mesh DM :: mech_mesh
PetscFE :: mechFE PetscFE :: mechFE
PetscQuadrature :: mechQuad, functional PetscQuadrature :: mechQuad, functional
PetscDS :: mechDS PetscDS :: mechDS
PetscDualSpace :: mechDualSpace PetscDualSpace :: mechDualSpace
DMLabel :: BCLabel DMLabel :: BCLabel
PetscInt, dimension(:), pointer :: pNumComp, pNumDof, pBcField, pBcPoint PetscInt, dimension(:), pointer :: pNumComp, pNumDof, pBcField, pBcPoint
PetscInt :: numBC, bcSize, nc PetscInt :: numBC, bcSize, nc, &
field, faceSet, topologDim, nNodalPoints, &
cellStart, cellEnd, cell, basis
IS :: bcPoint IS :: bcPoint
IS, pointer :: pBcComps(:), pBcPoints(:) IS, dimension(:), pointer :: pBcComps, pBcPoints
PetscSection :: section PetscSection :: section
PetscInt :: field, faceSet, topologDim, nNodalPoints
PetscReal, dimension(:), pointer :: qPointsP, qWeightsP, & PetscReal, dimension(:), pointer :: qPointsP, qWeightsP, &
nodalPointsP, nodalWeightsP nodalPointsP, nodalWeightsP,pV0, pCellJ, pInvcellJ
PetscReal, allocatable, target :: nodalPoints(:), nodalWeights(:)
PetscScalar, pointer :: px_scal(:)
PetscScalar, allocatable, target :: x_scal(:)
PetscReal :: detJ PetscReal :: detJ
PetscReal, allocatable, target :: cellJMat(:,:) PetscReal, allocatable, target :: cellJMat(:,:)
PetscReal, pointer :: pV0(:), pCellJ(:), pInvcellJ(:)
PetscInt :: cellStart, cellEnd, cell, basis PetscScalar, pointer :: px_scal(:)
character(len=7), parameter :: prefix = 'mechFE_' PetscScalar, allocatable, target :: x_scal(:)
character(len=*), parameter :: prefix = 'mechFE_'
PetscErrorCode :: ierr PetscErrorCode :: ierr
write(6,'(/,a)') ' <<<+- FEM_mech init -+>>>' write(6,'(/,a)') ' <<<+- FEM_mech init -+>>>'
@ -196,13 +200,11 @@ subroutine FEM_mech_init(fieldBC)
call VecSet(solution ,0.0,ierr); CHKERRQ(ierr) call VecSet(solution ,0.0,ierr); CHKERRQ(ierr)
call VecSet(solution_rate ,0.0,ierr); CHKERRQ(ierr) call VecSet(solution_rate ,0.0,ierr); CHKERRQ(ierr)
allocate(x_scal(cellDof)) allocate(x_scal(cellDof))
allocate(nodalPoints (dimPlex)) allocate(nodalWeightsP(1))
allocate(nodalWeights(1)) allocate(nodalPointsP(dimPlex))
nodalPointsP => nodalPoints
nodalWeightsP => nodalWeights
allocate(pv0(dimPlex)) allocate(pv0(dimPlex))
allocate(pcellJ(dimPlex*dimPlex)) allocate(pcellJ(dimPlex**2))
allocate(pinvcellJ(dimPlex*dimPlex)) allocate(pinvcellJ(dimPlex**2))
allocate(cellJMat(dimPlex,dimPlex)) allocate(cellJMat(dimPlex,dimPlex))
call DMGetSection(mech_mesh,section,ierr); CHKERRQ(ierr) call DMGetSection(mech_mesh,section,ierr); CHKERRQ(ierr)
call DMGetDS(mech_mesh,mechDS,ierr); CHKERRQ(ierr) call DMGetDS(mech_mesh,mechDS,ierr); CHKERRQ(ierr)
@ -212,7 +214,7 @@ subroutine FEM_mech_init(fieldBC)
call DMPlexGetHeightStratum(mech_mesh,0,cellStart,cellEnd,ierr) call DMPlexGetHeightStratum(mech_mesh,0,cellStart,cellEnd,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
do cell = cellStart, cellEnd-1 !< loop over all elements do cell = cellStart, cellEnd-1 !< loop over all elements
x_scal = 0.0 x_scal = 0.0_pReal
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])
@ -221,7 +223,7 @@ subroutine FEM_mech_init(fieldBC)
CHKERRQ(ierr) CHKERRQ(ierr)
call PetscQuadratureGetData(functional,dimPlex,nc,nNodalPoints,nodalPointsP,nodalWeightsP,ierr) call PetscQuadratureGetData(functional,dimPlex,nc,nNodalPoints,nodalPointsP,nodalWeightsP,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
x_scal(basis+1:basis+dimPlex) = pV0 + matmul(transpose(cellJMat),nodalPointsP + 1.0) x_scal(basis+1:basis+dimPlex) = pV0 + matmul(transpose(cellJMat),nodalPointsP + 1.0_pReal)
enddo enddo
px_scal => x_scal px_scal => x_scal
call DMPlexVecSetClosure(mech_mesh,section,solution_local,cell,px_scal,INSERT_ALL_VALUES,ierr) call DMPlexVecSetClosure(mech_mesh,section,solution_local,cell,px_scal,INSERT_ALL_VALUES,ierr)
@ -283,6 +285,9 @@ end function FEM_mech_solution
subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr) subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr)
DM :: dm_local DM :: dm_local
PetscObject,intent(in) :: dummy
PetscErrorCode :: ierr
PetscDS :: prob PetscDS :: prob
Vec :: x_local, f_local, xx_local Vec :: x_local, f_local, xx_local
PetscSection :: section PetscSection :: section
@ -294,10 +299,10 @@ subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr)
qPt, basis, comp, cidx qPt, basis, comp, cidx
PetscReal :: detFAvg PetscReal :: detFAvg
PetscReal :: BMat(dimPlex*dimPlex,cellDof) PetscReal :: BMat(dimPlex*dimPlex,cellDof)
PetscObject,intent(in) :: dummy
PetscInt :: bcSize PetscInt :: bcSize
IS :: bcPoints IS :: bcPoints
PetscErrorCode :: ierr
allocate(pV0(dimPlex)) allocate(pV0(dimPlex))
allocate(pcellJ(dimPlex**2)) allocate(pcellJ(dimPlex**2))
@ -316,7 +321,7 @@ subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr)
call DMGetStratumIS(dm_local,'Face Sets',mesh_boundaries(face),bcPoints,ierr) call DMGetStratumIS(dm_local,'Face Sets',mesh_boundaries(face),bcPoints,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call utilities_projectBCValues(x_local,section,0,field-1,bcPoints, & call utilities_projectBCValues(x_local,section,0,field-1,bcPoints, &
0.0,params%fieldBC%componentBC(field)%Value(face),params%timeinc) 0.0_pReal,params%fieldBC%componentBC(field)%Value(face),params%timeinc)
call ISDestroy(bcPoints,ierr); CHKERRQ(ierr) call ISDestroy(bcPoints,ierr); CHKERRQ(ierr)
endif endif
endif endif
@ -403,29 +408,36 @@ end subroutine FEM_mech_formResidual
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr) subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr)
DM :: dm_local
DM :: dm_local
Mat :: Jac_pre, Jac
PetscObject, intent(in) :: dummy
PetscErrorCode :: ierr
PetscDS :: prob PetscDS :: prob
Vec :: x_local, xx_local Vec :: x_local, xx_local
Mat :: Jac_pre, Jac
PetscSection :: section, gSection PetscSection :: section, gSection
PetscReal, dimension(1, cellDof) :: MatB
PetscReal, dimension(dimPlex**2,cellDof) :: BMat, BMatAvg, MatA
PetscReal, dimension(3,3) :: F, FAvg, FInv
PetscReal :: detJ PetscReal :: detJ
PetscReal, dimension(:), pointer :: basisField, basisFieldDer, & PetscReal, dimension(:), pointer :: basisField, basisFieldDer, &
pV0, pCellJ, pInvcellJ pV0, pCellJ, pInvcellJ
PetscInt :: cellStart, cellEnd, cell, field, face, &
qPt, basis, comp, cidx,bcSize PetscScalar, dimension(:), pointer :: pK_e, x_scal
PetscScalar,dimension(cellDOF,cellDOF), target :: K_e, & PetscScalar,dimension(cellDOF,cellDOF), target :: K_e, &
K_eA , & K_eA , &
K_eB K_eB
PetscScalar, target :: K_eVec(cellDof*cellDof) PetscScalar,(cellDof**2) ,target :: K_eVec
PetscReal :: BMat (dimPlex*dimPlex,cellDof), &
BMatAvg(dimPlex*dimPlex,cellDof), & PetscInt :: cellStart, cellEnd, cell, field, face, &
MatA (dimPlex*dimPlex,cellDof), & qPt, basis, comp, cidx,bcSize
MatB (1 ,cellDof)
PetscScalar, dimension(:), pointer :: pK_e, x_scal
PetscReal, dimension(3,3) :: F, FAvg, FInv
PetscObject, intent(in) :: dummy
IS :: bcPoints IS :: bcPoints
PetscErrorCode :: ierr
allocate(pV0(dimPlex)) allocate(pV0(dimPlex))
allocate(pcellJ(dimPlex**2)) allocate(pcellJ(dimPlex**2))
@ -440,7 +452,7 @@ subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr)
call DMGetGlobalSection(dm_local,gSection,ierr); CHKERRQ(ierr) call DMGetGlobalSection(dm_local,gSection,ierr); CHKERRQ(ierr)
call DMGetLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr) call DMGetLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr)
call VecWAXPY(x_local,1.0,xx_local,solution_local,ierr); CHKERRQ(ierr) call VecWAXPY(x_local,1.0_pReal,xx_local,solution_local,ierr); CHKERRQ(ierr)
do field = 1, dimPlex; do face = 1, mesh_Nboundaries do field = 1, dimPlex; do face = 1, mesh_Nboundaries
if (params%fieldBC%componentBC(field)%Mask(face)) then if (params%fieldBC%componentBC(field)%Mask(face)) then
call DMGetStratumSize(dm_local,'Face Sets',mesh_boundaries(face),bcSize,ierr) call DMGetStratumSize(dm_local,'Face Sets',mesh_boundaries(face),bcSize,ierr)
@ -448,7 +460,7 @@ subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr)
call DMGetStratumIS(dm_local,'Face Sets',mesh_boundaries(face),bcPoints,ierr) call DMGetStratumIS(dm_local,'Face Sets',mesh_boundaries(face),bcPoints,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call utilities_projectBCValues(x_local,section,0,field-1,bcPoints, & call utilities_projectBCValues(x_local,section,0,field-1,bcPoints, &
0.0,params%fieldBC%componentBC(field)%Value(face),params%timeinc) 0.0_pReal,params%fieldBC%componentBC(field)%Value(face),params%timeinc)
call ISDestroy(bcPoints,ierr); CHKERRQ(ierr) call ISDestroy(bcPoints,ierr); CHKERRQ(ierr)
endif endif
endif endif
@ -501,7 +513,6 @@ subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr)
(matmul(matmul(transpose(BMatAvg), & (matmul(matmul(transpose(BMatAvg), &
reshape(FInv(1:dimPlex,1:dimPlex),shape=[dimPlex*dimPlex,1],order=[2,1])),MatB) + & reshape(FInv(1:dimPlex,1:dimPlex),shape=[dimPlex*dimPlex,1],order=[2,1])),MatB) + &
K_eB)/real(dimPlex) K_eB)/real(dimPlex)
else else
K_e = K_eA K_e = K_eA
endif endif
@ -536,18 +547,18 @@ subroutine FEM_mech_forward(guess,timeinc,timeinc_old,fieldBC)
type(tFieldBC), intent(in) :: & type(tFieldBC), intent(in) :: &
fieldBC fieldBC
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
timeinc_old, & timeinc_old, &
timeinc timeinc
logical, intent(in) :: & logical, intent(in) :: &
guess guess
PetscInt :: field, face
DM :: dm_local PetscInt :: field, face, bcSize
Vec :: x_local DM :: dm_local
PetscSection :: section Vec :: x_local
PetscInt :: bcSize PetscSection :: section
IS :: bcPoints IS :: bcPoints
PetscErrorCode :: ierr PetscErrorCode :: ierr
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! forward last inc ! forward last inc
@ -557,7 +568,7 @@ subroutine FEM_mech_forward(guess,timeinc,timeinc_old,fieldBC)
call SNESGetDM(mech_snes,dm_local,ierr); CHKERRQ(ierr) !< retrieve mesh info from mech_snes into dm_local call SNESGetDM(mech_snes,dm_local,ierr); CHKERRQ(ierr) !< retrieve mesh info from mech_snes into dm_local
call DMGetSection(dm_local,section,ierr); CHKERRQ(ierr) call DMGetSection(dm_local,section,ierr); CHKERRQ(ierr)
call DMGetLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr) call DMGetLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr)
call VecSet(x_local,0.0,ierr); CHKERRQ(ierr) call VecSet(x_local,0.00_pReal,ierr); CHKERRQ(ierr)
call DMGlobalToLocalBegin(dm_local,solution,INSERT_VALUES,x_local,ierr) !< retrieve my partition of global solution vector call DMGlobalToLocalBegin(dm_local,solution,INSERT_VALUES,x_local,ierr) !< retrieve my partition of global solution vector
CHKERRQ(ierr) CHKERRQ(ierr)
call DMGlobalToLocalEnd(dm_local,solution,INSERT_VALUES,x_local,ierr) call DMGlobalToLocalEnd(dm_local,solution,INSERT_VALUES,x_local,ierr)
@ -570,7 +581,7 @@ subroutine FEM_mech_forward(guess,timeinc,timeinc_old,fieldBC)
call DMGetStratumIS(dm_local,'Face Sets',mesh_boundaries(face),bcPoints,ierr) call DMGetStratumIS(dm_local,'Face Sets',mesh_boundaries(face),bcPoints,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call utilities_projectBCValues(solution_local,section,0,field-1,bcPoints, & call utilities_projectBCValues(solution_local,section,0,field-1,bcPoints, &
0.0,fieldBC%componentBC(field)%Value(face),timeinc_old) 0.0_pReal,fieldBC%componentBC(field)%Value(face),timeinc_old)
call ISDestroy(bcPoints,ierr); CHKERRQ(ierr) call ISDestroy(bcPoints,ierr); CHKERRQ(ierr)
endif endif
endif endif

16
src/mesh/FEM_utilities.f90
View File

@ -23,7 +23,7 @@ module FEM_utilities
private private
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
logical, public :: cutBack = .false. !< cut back of BVP solver in case convergence is not achieved or a material point is terminally ill logical, public :: cutBack = .false. !< cut back of BVP solver in case convergence is not achieved or a material point is terminally ill
integer, public, parameter :: maxFields = 6 integer, public, parameter :: maxFields = 6
integer, public :: nActiveFields = 0 integer, public :: nActiveFields = 0
@ -56,15 +56,15 @@ module FEM_utilities
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! derived types ! derived types
type, public :: tSolutionState !< return type of solution from FEM solver variants type, public :: tSolutionState !< return type of solution from FEM solver variants
logical :: converged = .true. logical :: converged = .true.
logical :: stagConverged = .true. logical :: stagConverged = .true.
integer :: iterationsNeeded = 0 integer :: iterationsNeeded = 0
end type tSolutionState end type tSolutionState
type, public :: tComponentBC type, public :: tComponentBC
integer(kind(COMPONENT_UNDEFINED_ID)) :: ID integer(kind(COMPONENT_UNDEFINED_ID)) :: ID
real(pReal), allocatable :: Value(:) real(pReal), allocatable, dimension(:) :: Value
logical, allocatable :: Mask(:) logical, allocatable, dimension(:) :: Mask
end type tComponentBC end type tComponentBC
type, public :: tFieldBC type, public :: tFieldBC
@ -79,8 +79,8 @@ module FEM_utilities
outputfrequency = 1, & !< frequency of result writes outputfrequency = 1, & !< frequency of result writes
logscale = 0 !< linear/logarithmic time inc flag logscale = 0 !< linear/logarithmic time inc flag
logical :: followFormerTrajectory = .true. !< follow trajectory of former loadcase logical :: followFormerTrajectory = .true. !< follow trajectory of former loadcase
integer, allocatable :: faceID(:) integer, allocatable, dimension(:) :: faceID
type(tFieldBC), allocatable :: fieldBC(:) type(tFieldBC), allocatable, dimension(:) :: fieldBC
end type tLoadCase end type tLoadCase
public :: & public :: &

4
src/mesh/FEM_zoo.f90
View File

@ -37,7 +37,7 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine FEM_Zoo_init subroutine FEM_Zoo_init
write(6,'(/,a)') ' <<<+- FEM_Zoo init -+>>>' write(6,'(/,a)') ' <<<+- FEM_Zoo init -+>>>'; flush(6)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! 2D linear ! 2D linear
@ -53,7 +53,7 @@ subroutine FEM_Zoo_init
! 2D quadratic ! 2D quadratic
FEM_Zoo_nQuadrature(2,2) = 3 FEM_Zoo_nQuadrature(2,2) = 3
allocate(FEM_Zoo_QuadratureWeights(2,2)%p(3)) allocate(FEM_Zoo_QuadratureWeights(2,2)%p(3))
FEM_Zoo_QuadratureWeights(2,2)%p(1:3) = 1.0_pReal/3.0_pReal FEM_Zoo_QuadratureWeights(2,2)%p(1:3) = 1.0_pReal/3.0_pReal
allocate(FEM_Zoo_QuadraturePoints (2,2)%p(6)) allocate(FEM_Zoo_QuadraturePoints (2,2)%p(6))