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

Merge branch 'mesh-clean' into 'development' 

Mesh clean

See merge request damask/DAMASK!889
This commit is contained in:
Javier Velo 2024-01-01 18:43:09 +00:00
parent 14bafc62db 42dd33d14c
commit f57924e9fe
3 changed files with 44 additions and 76 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

47
src/mesh/DAMASK_mesh.f90
View File

@ -25,10 +25,9 @@ program DAMASK_mesh
type :: tLoadCase type :: tLoadCase
real(pREAL) :: t = 0.0_pREAL !< length of increment real(pREAL) :: t = 0.0_pREAL !< length of increment
integer :: N = 0, & !< number of increments integer :: N = 0, & !< number of increments
f_out = 1 !< frequency of result writes f_out = 1 !< frequency of result writes
logical :: estimate_rate = .true. !< follow trajectory of former loadcase logical :: estimate_rate = .true. !< follow trajectory of former loadcase
integer, allocatable, dimension(:) :: tag type(tMechBC), allocatable, dimension(:) :: mechBC
type(tMechBC) :: mechBC
end type tLoadCase end type tLoadCase
@ -40,20 +39,16 @@ program DAMASK_mesh
t = 0.0_pREAL, & !< elapsed time t = 0.0_pREAL, & !< elapsed time
t_0 = 0.0_pREAL, & !< begin of interval t_0 = 0.0_pREAL, & !< begin of interval
Delta_t = 0.0_pREAL, & !< current time interval Delta_t = 0.0_pREAL, & !< current time interval
Delta_t_prev = 0.0_pREAL, & !< previous time interval Delta_t_prev = 0.0_pREAL !< previous time interval
t_remaining = 0.0_pREAL !< remaining time of current load case
logical :: & logical :: &
guess, & !< guess along former trajectory guess, & !< guess along former trajectory
stagIterate stagIterate
integer :: & integer :: &
l, & l, &
i, &
m, & m, &
errorID, & errorID, &
cutBackLevel = 0, & !< cut back level \f$ t = \frac{t_{inc}}{2^l} \f$ cutBackLevel = 0, & !< cut back level \f$ t = \frac{t_{inc}}{2^l} \f$
stepFraction = 0, & !< fraction of current time interval stepFraction = 0, & !< fraction of current time interval
currentLoadcase = 0, & !< current load case
currentFace = 0, &
inc, & !< current increment in current load case inc, & !< current increment in current load case
totalIncsCounter = 0, & !< total # of increments totalIncsCounter = 0, & !< total # of increments
statUnit = 0, & !< file unit for statistics output statUnit = 0, & !< file unit for statistics output
@ -82,7 +77,6 @@ program DAMASK_mesh
type(tSolutionState), allocatable, dimension(:) :: solres type(tSolutionState), allocatable, dimension(:) :: solres
PetscInt :: faceSet, currentFaceSet, dimPlex PetscInt :: faceSet, currentFaceSet, dimPlex
PetscErrorCode :: err_PETSc PetscErrorCode :: err_PETSc
integer(kind(COMPONENT_UNDEFINED_ID)) :: ID
external :: & external :: &
quit quit
character(len=:), allocatable :: & character(len=:), allocatable :: &
@ -124,25 +118,16 @@ program DAMASK_mesh
allocate(loadCases(load_steps%length)) allocate(loadCases(load_steps%length))
do l = 1, load_steps%length do l = 1, load_steps%length
load_step => load_steps%get_dict(l) load_step => load_steps%get_dict(l)
step_bc => load_step%get_dict('boundary_conditions') step_bc => load_step%get_dict('boundary_conditions')
step_mech => step_bc%get_list('mechanical') step_mech => step_bc%get_list('mechanical')
loadCases(l)%mechBC%nComponents = dimPlex !< X, Y (, Z) displacements allocate(loadCases(l)%mechBC(mesh_Nboundaries))
allocate(loadCases(l)%mechBC%componentBC(dimPlex)) loadCases(l)%mechBC(:)%nComponents = dimPlex !< X, Y (, Z) displacements
do component = 1, dimPlex do faceSet = 1, mesh_Nboundaries
select case (component) allocate(loadCases(l)%mechBC(faceSet)%Value(dimPlex), source = 0.0_pREAL)
case (1) allocate(loadCases(l)%mechBC(faceSet)%Mask(dimPlex), source = .false.)
loadCases(l)%mechBC%componentBC(component)%ID = COMPONENT_MECH_X_ID
case (2)
loadCases(l)%mechBC%componentBC(component)%ID = COMPONENT_MECH_Y_ID
case (3)
loadCases(l)%mechBC%componentBC(component)%ID = COMPONENT_MECH_Z_ID
end select
end do
do component = 1, dimPlex
allocate(loadCases(l)%mechBC%componentBC(component)%Value(mesh_Nboundaries), source = 0.0_pREAL)
allocate(loadCases(l)%mechBC%componentBC(component)%Mask (mesh_Nboundaries), source = .false.)
end do end do
do m = 1, step_mech%length do m = 1, step_mech%length
@ -152,11 +137,11 @@ program DAMASK_mesh
if (mesh_boundaries(faceSet) == mech_BC%get_asInt('tag')) currentFaceSet = faceSet if (mesh_boundaries(faceSet) == mech_BC%get_asInt('tag')) currentFaceSet = faceSet
end do end do
if (currentFaceSet < 0) call IO_error(error_ID = 837, ext_msg = 'invalid BC') if (currentFaceSet < 0) call IO_error(error_ID = 837, ext_msg = 'invalid BC')
mech_u => mech_BC%get_list('dot_u')
do component = 1, dimPlex do component = 1, dimPlex
mech_u => mech_BC%get_list('dot_u')
if (mech_u%get_asStr(component) /= 'x') then if (mech_u%get_asStr(component) /= 'x') then
loadCases(l)%mechBC%componentBC(component)%Mask(currentFaceSet) = .true. loadCases(l)%mechBC(currentFaceSet)%Mask(component) = .true.
loadCases(l)%mechBC%componentBC(component)%Value(currentFaceSet) = mech_u%get_asReal(component) loadCases(l)%mechBC(currentFaceSet)%Value(component) = mech_u%get_asReal(component)
end if end if
end do end do
end do end do
@ -183,12 +168,12 @@ program DAMASK_mesh
print'(2x,a)', 'Field '//trim(FIELD_MECH_label) print'(2x,a)', 'Field '//trim(FIELD_MECH_label)
do faceSet = 1, mesh_Nboundaries do faceSet = 1, mesh_Nboundaries
do component = 1, loadCases(l)%mechBC%nComponents do component = 1, dimPlex
if (loadCases(l)%mechBC%componentBC(component)%Mask(faceSet)) & if (loadCases(l)%mechBC(faceSet)%Mask(component)) &
print'(a,i2,a,i2,a,f12.7)', & print'(a,i2,a,i2,a,f12.7)', &
' Face ', mesh_boundaries(faceSet), & ' Face ', mesh_boundaries(faceSet), &
' Component ', component, & ' Component ', component, &
' Value ', loadCases(l)%mechBC%componentBC(component)%Value(faceSet) ' Value ', loadCases(l)%mechBC(faceSet)%Value(component)
end do end do
end do end do
print'(2x,a,f12.6)', 'time: ', loadCases(l)%t print'(2x,a,f12.6)', 'time: ', loadCases(l)%t
@ -230,7 +215,6 @@ program DAMASK_mesh
stepFraction = 0 ! fraction scaled by stepFactor**cutLevel stepFraction = 0 ! fraction scaled by stepFactor**cutLevel
subStepLooping: do while (stepFraction < subStepFactor**cutBackLevel) subStepLooping: do while (stepFraction < subStepFactor**cutBackLevel)
t_remaining = loadCases(l)%t + t_0 - t
t = t + Delta_t ! forward target time t = t + Delta_t ! forward target time
stepFraction = stepFraction + 1 ! count step stepFraction = stepFraction + 1 ! count step
@ -315,4 +299,5 @@ program DAMASK_mesh
call quit(0) ! no complains ;) call quit(0) ! no complains ;)
end program DAMASK_mesh end program DAMASK_mesh

23
src/mesh/FEM_utilities.f90
View File

@ -39,14 +39,6 @@ module FEM_utilities
character(len=*), parameter, public :: & character(len=*), parameter, public :: &
FIELD_MECH_label = 'mechanical' FIELD_MECH_label = 'mechanical'
enum, bind(c); enumerator :: &
COMPONENT_UNDEFINED_ID, &
COMPONENT_MECH_X_ID, &
COMPONENT_MECH_Y_ID, &
COMPONENT_MECH_Z_ID
end enum
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! 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
@ -55,15 +47,10 @@ module FEM_utilities
PetscInt :: iterationsNeeded = 0_pPETSCINT PetscInt :: iterationsNeeded = 0_pPETSCINT
end type tSolutionState end type tSolutionState
type, public :: tComponentBC
integer(kind(COMPONENT_UNDEFINED_ID)) :: ID
real(pREAL), allocatable, dimension(:) :: Value
logical, allocatable, dimension(:) :: Mask
end type tComponentBC
type, public :: tMechBC type, public :: tMechBC
integer :: nComponents = 0 integer :: nComponents = 0
type(tComponentBC), allocatable, dimension(:) :: componentBC real(pREAL), allocatable, dimension(:) :: Value
logical, allocatable, dimension(:) :: Mask
end type tMechBC end type tMechBC
external :: & ! ToDo: write interfaces external :: & ! ToDo: write interfaces
@ -74,11 +61,7 @@ module FEM_utilities
public :: & public :: &
FEM_utilities_init, & FEM_utilities_init, &
utilities_constitutiveResponse, & utilities_constitutiveResponse, &
utilities_projectBCValues, & utilities_projectBCValues
COMPONENT_UNDEFINED_ID, &
COMPONENT_MECH_X_ID, &
COMPONENT_MECH_Y_ID, &
COMPONENT_MECH_Z_ID
contains contains

50
src/mesh/mesh_mech_FEM.f90
View File

@ -36,7 +36,7 @@ module mesh_mechanical_FEM
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! derived types ! derived types
type tSolutionParams type tSolutionParams
type(tMechBC) :: mechBC type(tMechBC),allocatable, dimension(:) :: mechBC
real(pREAL) :: Delta_t real(pREAL) :: Delta_t
end type tSolutionParams end type tSolutionParams
@ -99,7 +99,7 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine FEM_mechanical_init(mechBC,num_mesh) subroutine FEM_mechanical_init(mechBC,num_mesh)
type(tMechBC), intent(in) :: mechBC type(tMechBC), dimension(:), intent(in):: mechBC
type(tDict), pointer, intent(in) :: num_mesh type(tDict), pointer, intent(in) :: num_mesh
DM :: mechanical_mesh DM :: mechanical_mesh
@ -112,7 +112,7 @@ subroutine FEM_mechanical_init(mechBC,num_mesh)
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, & component, faceSet, topologDim, nNodalPoints, &
cellStart, cellEnd, cell, basis cellStart, cellEnd, cell, basis
IS :: bcPoint IS :: bcPoint
@ -208,17 +208,17 @@ subroutine FEM_mechanical_init(mechBC,num_mesh)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
end do end do
numBC = 0 numBC = 0
do field = 1, dimPlex; do faceSet = 1, mesh_Nboundaries do faceSet = 1, mesh_Nboundaries; do component = 1, dimPlex
if (mechBC%componentBC(field)%Mask(faceSet)) numBC = numBC + 1 if (mechBC(faceSet)%Mask(component)) numBC = numBC + 1
end do; end do end do; end do
allocate(pbcField(numBC), source=0_pPETSCINT) allocate(pbcField(numBC), source=0_pPETSCINT)
allocate(pbcComps(numBC)) allocate(pbcComps(numBC))
allocate(pbcPoints(numBC)) allocate(pbcPoints(numBC))
numBC = 0 numBC = 0
do field = 1, dimPlex; do faceSet = 1, mesh_Nboundaries do faceSet = 1, mesh_Nboundaries; do component = 1, dimPlex
if (mechBC%componentBC(field)%Mask(faceSet)) then if (mechBC(faceSet)%Mask(component)) then
numBC = numBC + 1 numBC = numBC + 1
call ISCreateGeneral(PETSC_COMM_WORLD,1_pPETSCINT,[field-1],PETSC_COPY_VALUES,pbcComps(numBC),err_PETSc) call ISCreateGeneral(PETSC_COMM_WORLD,1_pPETSCINT,[component-1],PETSC_COPY_VALUES,pbcComps(numBC),err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call DMGetStratumSize(mechanical_mesh,'Face Sets',mesh_boundaries(faceSet),bcSize,err_PETSc) call DMGetStratumSize(mechanical_mesh,'Face Sets',mesh_boundaries(faceSet),bcSize,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
@ -327,7 +327,7 @@ type(tSolutionState) function FEM_mechanical_solution( &
real(pREAL), intent(in) :: & real(pREAL), intent(in) :: &
Delta_t, & !< increment in time for current solution Delta_t, & !< increment in time for current solution
Delta_t_prev !< increment in time of last increment Delta_t_prev !< increment in time of last increment
type(tMechBC), intent(in) :: & type(tMechBC), dimension(:),intent(in) :: &
mechBC mechBC
character(len=*), intent(in) :: & character(len=*), intent(in) :: &
incInfoIn incInfoIn
@ -380,7 +380,7 @@ subroutine FEM_mechanical_formResidual(dm_local,xx_local,f_local,dummy,err_PETSc
real(pREAL), dimension(cellDof), target :: f_scal real(pREAL), dimension(cellDof), target :: f_scal
PetscReal :: IcellJMat(dimPlex,dimPlex) PetscReal :: IcellJMat(dimPlex,dimPlex)
PetscReal, dimension(:),pointer :: pV0, pCellJ, pInvcellJ, basisField, basisFieldDer PetscReal, dimension(:),pointer :: pV0, pCellJ, pInvcellJ, basisField, basisFieldDer
PetscInt :: cellStart, cellEnd, cell, field, face, & PetscInt :: cellStart, cellEnd, cell, component, face, &
qPt, basis, comp, cidx, & qPt, basis, comp, cidx, &
numFields, & numFields, &
bcSize,m,i bcSize,m,i
@ -406,14 +406,14 @@ subroutine FEM_mechanical_formResidual(dm_local,xx_local,f_local,dummy,err_PETSc
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call VecWAXPY(x_local,1.0_pREAL,xx_local,solution_local,err_PETSc) call VecWAXPY(x_local,1.0_pREAL,xx_local,solution_local,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
do field = 1_pPETSCINT, dimPlex; do face = 1_pPETSCINT, mesh_Nboundaries do face = 1_pPETSCINT, mesh_Nboundaries; do component = 1_pPETSCINT, dimPlex
if (params%mechBC%componentBC(field)%Mask(face)) then if (params%mechBC(face)%Mask(component)) then
call DMGetStratumSize(dm_local,'Face Sets',mesh_boundaries(face),bcSize,err_PETSc) call DMGetStratumSize(dm_local,'Face Sets',mesh_boundaries(face),bcSize,err_PETSc)
if (bcSize > 0) then if (bcSize > 0) then
call DMGetStratumIS(dm_local,'Face Sets',mesh_boundaries(face),bcPoints,err_PETSc) call DMGetStratumIS(dm_local,'Face Sets',mesh_boundaries(face),bcPoints,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call utilities_projectBCValues(x_local,section,0_pPETSCINT,field-1,bcPoints, & call utilities_projectBCValues(x_local,section,0_pPETSCINT,component-1,bcPoints, &
0.0_pREAL,params%mechBC%componentBC(field)%Value(face),params%Delta_t) 0.0_pREAL,params%mechBC(face)%Value(component),params%Delta_t)
call ISDestroy(bcPoints,err_PETSc) call ISDestroy(bcPoints,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
end if end if
@ -529,7 +529,7 @@ subroutine FEM_mechanical_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,err_P
real(pREAL),dimension(cellDOF,cellDOF), target :: K_e real(pREAL),dimension(cellDOF,cellDOF), target :: K_e
real(pREAL),dimension(cellDOF,cellDOF) :: K_eA, K_eB real(pREAL),dimension(cellDOF,cellDOF) :: K_eA, K_eB
PetscInt :: cellStart, cellEnd, cell, field, face, & PetscInt :: cellStart, cellEnd, cell, component, face, &
qPt, basis, comp, cidx,bcSize, m, i qPt, basis, comp, cidx,bcSize, m, i
IS :: bcPoints IS :: bcPoints
@ -556,14 +556,14 @@ subroutine FEM_mechanical_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,err_P
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call VecWAXPY(x_local,1.0_pREAL,xx_local,solution_local,err_PETSc) call VecWAXPY(x_local,1.0_pREAL,xx_local,solution_local,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
do field = 1, dimPlex; do face = 1, mesh_Nboundaries do face = 1, mesh_Nboundaries; do component = 1, dimPlex
if (params%mechBC%componentBC(field)%Mask(face)) then if (params%mechBC(face)%Mask(component)) then
call DMGetStratumSize(dm_local,'Face Sets',mesh_boundaries(face),bcSize,err_PETSc) call DMGetStratumSize(dm_local,'Face Sets',mesh_boundaries(face),bcSize,err_PETSc)
if (bcSize > 0) then if (bcSize > 0) then
call DMGetStratumIS(dm_local,'Face Sets',mesh_boundaries(face),bcPoints,err_PETSc) call DMGetStratumIS(dm_local,'Face Sets',mesh_boundaries(face),bcPoints,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call utilities_projectBCValues(x_local,section,0_pPETSCINT,field-1,bcPoints, & call utilities_projectBCValues(x_local,section,0_pPETSCINT,component-1,bcPoints, &
0.0_pREAL,params%mechBC%componentBC(field)%Value(face),params%Delta_t) 0.0_pREAL,params%mechBC(face)%Value(component),params%Delta_t)
call ISDestroy(bcPoints,err_PETSc) call ISDestroy(bcPoints,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
end if end if
@ -667,7 +667,7 @@ end subroutine FEM_mechanical_formJacobian
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine FEM_mechanical_forward(guess,Delta_t,Delta_t_prev,mechBC) subroutine FEM_mechanical_forward(guess,Delta_t,Delta_t_prev,mechBC)
type(tMechBC), intent(in) :: & type(tMechBC), dimension(:), intent(in) :: &
mechBC mechBC
real(pREAL), intent(in) :: & real(pREAL), intent(in) :: &
Delta_t_prev, & Delta_t_prev, &
@ -675,7 +675,7 @@ subroutine FEM_mechanical_forward(guess,Delta_t,Delta_t_prev,mechBC)
logical, intent(in) :: & logical, intent(in) :: &
guess guess
PetscInt :: field, face, bcSize PetscInt :: component, face, bcSize
DM :: dm_local DM :: dm_local
Vec :: x_local Vec :: x_local
PetscSection :: section PetscSection :: section
@ -700,14 +700,14 @@ subroutine FEM_mechanical_forward(guess,Delta_t,Delta_t_prev,mechBC)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call VecAXPY(solution_local,1.0_pREAL,x_local,err_PETSc) call VecAXPY(solution_local,1.0_pREAL,x_local,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
do field = 1, dimPlex; do face = 1, mesh_Nboundaries do face = 1, mesh_Nboundaries; do component = 1, dimPlex
if (mechBC%componentBC(field)%Mask(face)) then if (mechBC(face)%Mask(component)) then
call DMGetStratumSize(dm_local,'Face Sets',mesh_boundaries(face),bcSize,err_PETSc) call DMGetStratumSize(dm_local,'Face Sets',mesh_boundaries(face),bcSize,err_PETSc)
if (bcSize > 0) then if (bcSize > 0) then
call DMGetStratumIS(dm_local,'Face Sets',mesh_boundaries(face),bcPoints,err_PETSc) call DMGetStratumIS(dm_local,'Face Sets',mesh_boundaries(face),bcPoints,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call utilities_projectBCValues(solution_local,section,0_pPETSCINT,field-1,bcPoints, & call utilities_projectBCValues(solution_local,section,0_pPETSCINT,component-1,bcPoints, &
0.0_pREAL,mechBC%componentBC(field)%Value(face),Delta_t_prev) 0.0_pREAL,mechBC(face)%Value(component),Delta_t_prev)
call ISDestroy(bcPoints,err_PETSc) call ISDestroy(bcPoints,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
end if end if