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

this structure seems cleaner, i.e., BC per face

This commit is contained in:
Sharan Roongta 2023-12-17 00:04:55 +01:00
parent e4ad5a94a5
commit 0798548f0f
3 changed files with 37 additions and 64 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

44
src/mesh/DAMASK_mesh.f90
View File

@ -27,8 +27,7 @@ program DAMASK_mesh
integer :: N = 0, & !< number of increments
f_out = 1 !< frequency of result writes
logical :: estimate_rate = .true. !< follow trajectory of former loadcase
integer, allocatable, dimension(:) :: tag
type(tMechBC) :: mechBC
type(tMechBC), allocatable, dimension(:) :: mechBC
end type tLoadCase
@ -82,7 +81,6 @@ program DAMASK_mesh
type(tSolutionState), allocatable, dimension(:) :: solres
PetscInt :: faceSet, currentFaceSet, dimPlex
PetscErrorCode :: err_PETSc
integer(kind(COMPONENT_UNDEFINED_ID)) :: ID
external :: &
quit
character(len=:), allocatable :: &
@ -124,25 +122,16 @@ program DAMASK_mesh
allocate(loadCases(load_steps%length))
do l = 1, load_steps%length
load_step => load_steps%get_dict(l)
step_bc => load_step%get_dict('boundary_conditions')
step_mech => step_bc%get_list('mechanical')
loadCases(l)%mechBC%nComponents = dimPlex !< X, Y (, Z) displacements
allocate(loadCases(l)%mechBC%componentBC(dimPlex))
do component = 1, dimPlex
select case (component)
case (1)
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.)
allocate(loadCases(l)%mechBC(mesh_Nboundaries))
loadCases(l)%mechBC(:)%nComponents = dimPlex !< X, Y (, Z) displacements
do faceSet = 1, mesh_Nboundaries
allocate(loadCases(l)%mechBC(faceSet)%Value(dimPlex), source = 0.0_pREAL)
allocate(loadCases(l)%mechBC(faceSet)%Mask(dimPlex), source = .false.)
end do
do m = 1, step_mech%length
@ -152,11 +141,11 @@ program DAMASK_mesh
if (mesh_boundaries(faceSet) == mech_BC%get_asInt('tag')) currentFaceSet = faceSet
end do
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
mech_u => mech_BC%get_list('dot_u')
if (mech_u%get_asStr(component) /= 'x') then
loadCases(l)%mechBC%componentBC(component)%Mask(currentFaceSet) = .true.
loadCases(l)%mechBC%componentBC(component)%Value(currentFaceSet) = mech_u%get_asReal(component)
loadCases(l)%mechBC(currentFaceSet)%Mask(component) = .true.
loadCases(l)%mechBC(currentFaceSet)%Value(component) = mech_u%get_asReal(component)
end if
end do
end do
@ -183,12 +172,12 @@ program DAMASK_mesh
print'(2x,a)', 'Field '//trim(FIELD_MECH_label)
do faceSet = 1, mesh_Nboundaries
do component = 1, loadCases(l)%mechBC%nComponents
if (loadCases(l)%mechBC%componentBC(component)%Mask(faceSet)) &
do component = 1, dimPlex
if (loadCases(l)%mechBC(faceSet)%Mask(component)) &
print'(a,i2,a,i2,a,f12.7)', &
' Face ', mesh_boundaries(faceSet), &
' Component ', component, &
' Value ', loadCases(l)%mechBC%componentBC(component)%Value(faceSet)
' Value ', loadCases(l)%mechBC(faceSet)%Value(component)
end do
end do
print'(2x,a,f12.6)', 'time: ', loadCases(l)%t
@ -203,7 +192,7 @@ program DAMASK_mesh
!--------------------------------------------------------------------------------------------------
! doing initialization depending on active solvers
call FEM_Utilities_init(num_mesh)
call FEM_mechanical_init(loadCases(1)%mechBC,num_mesh)
call FEM_mechanical_init(loadCases(1)%mechBC(:),num_mesh)
call config_numerics_deallocate()
if (worldrank == 0) then
@ -247,14 +236,14 @@ program DAMASK_mesh
'-',stepFraction, '/', subStepFactor**cutBackLevel
flush(IO_STDOUT)
call FEM_mechanical_forward(guess,Delta_t,Delta_t_prev,loadCases(l)%mechBC)
call FEM_mechanical_forward(guess,Delta_t,Delta_t_prev,loadCases(l)%mechBC(:))
!--------------------------------------------------------------------------------------------------
! solve fields
stagIter = 0
stagIterate = .true.
do while (stagIterate)
solres(1) = FEM_mechanical_solution(incInfo,Delta_t,Delta_t_prev,loadCases(l)%mechBC)
solres(1) = FEM_mechanical_solution(incInfo,Delta_t,Delta_t_prev,loadCases(l)%mechBC(:))
if (.not. solres(1)%converged) exit
stagIter = stagIter + 1
@ -315,4 +304,5 @@ program DAMASK_mesh
call quit(0) ! no complains ;)
end program DAMASK_mesh

23
src/mesh/FEM_utilities.f90
View File

@ -38,14 +38,6 @@ module FEM_utilities
character(len=*), parameter, public :: &
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
type, public :: tSolutionState !< return type of solution from FEM solver variants
@ -54,15 +46,10 @@ module FEM_utilities
PetscInt :: iterationsNeeded = 0_pPETSCINT
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
integer :: nComponents = 0
type(tComponentBC), allocatable, dimension(:) :: componentBC
real(pREAL), allocatable, dimension(:) :: Value
logical, allocatable, dimension(:) :: Mask
end type tMechBC
external :: & ! ToDo: write interfaces
@ -73,11 +60,7 @@ module FEM_utilities
public :: &
FEM_utilities_init, &
utilities_constitutiveResponse, &
utilities_projectBCValues, &
COMPONENT_UNDEFINED_ID, &
COMPONENT_MECH_X_ID, &
COMPONENT_MECH_Y_ID, &
COMPONENT_MECH_Z_ID
utilities_projectBCValues
contains

34
src/mesh/mesh_mech_FEM.f90
View File

@ -36,7 +36,7 @@ module mesh_mechanical_FEM
!--------------------------------------------------------------------------------------------------
! derived types
type tSolutionParams
type(tMechBC) :: mechBC
type(tMechBC),allocatable, dimension(:) :: mechBC
real(pREAL) :: Delta_t
end type tSolutionParams
@ -99,7 +99,7 @@ contains
!--------------------------------------------------------------------------------------------------
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
DM :: mechanical_mesh
@ -208,15 +208,15 @@ subroutine FEM_mechanical_init(mechBC,num_mesh)
CHKERRQ(err_PETSc)
end do
numBC = 0
do field = 1, dimPlex; do faceSet = 1, mesh_Nboundaries
if (mechBC%componentBC(field)%Mask(faceSet)) numBC = numBC + 1
do faceSet = 1, mesh_Nboundaries; do field = 1, dimPlex
if (mechBC(faceSet)%Mask(field)) numBC = numBC + 1
end do; end do
allocate(pbcField(numBC), source=0_pPETSCINT)
allocate(pbcComps(numBC))
allocate(pbcPoints(numBC))
numBC = 0
do field = 1, dimPlex; do faceSet = 1, mesh_Nboundaries
if (mechBC%componentBC(field)%Mask(faceSet)) then
do faceSet = 1, mesh_Nboundaries; do field = 1, dimPlex
if (mechBC(faceSet)%Mask(field)) then
numBC = numBC + 1
call ISCreateGeneral(PETSC_COMM_WORLD,1_pPETSCINT,[field-1],PETSC_COPY_VALUES,pbcComps(numBC),err_PETSc)
CHKERRQ(err_PETSc)
@ -327,7 +327,7 @@ type(tSolutionState) function FEM_mechanical_solution( &
real(pREAL), intent(in) :: &
Delta_t, & !< increment in time for current solution
Delta_t_prev !< increment in time of last increment
type(tMechBC), intent(in) :: &
type(tMechBC), dimension(:),intent(in) :: &
mechBC
character(len=*), intent(in) :: &
incInfoIn
@ -406,14 +406,14 @@ subroutine FEM_mechanical_formResidual(dm_local,xx_local,f_local,dummy,err_PETSc
CHKERRQ(err_PETSc)
call VecWAXPY(x_local,1.0_pREAL,xx_local,solution_local,err_PETSc)
CHKERRQ(err_PETSc)
do field = 1_pPETSCINT, dimPlex; do face = 1_pPETSCINT, mesh_Nboundaries
if (params%mechBC%componentBC(field)%Mask(face)) then
do face = 1_pPETSCINT, mesh_Nboundaries; do field = 1_pPETSCINT, dimPlex
if (params%mechBC(face)%Mask(field)) then
call DMGetStratumSize(dm_local,'Face Sets',mesh_boundaries(face),bcSize,err_PETSc)
if (bcSize > 0) then
call DMGetStratumIS(dm_local,'Face Sets',mesh_boundaries(face),bcPoints,err_PETSc)
CHKERRQ(err_PETSc)
call utilities_projectBCValues(x_local,section,0_pPETSCINT,field-1,bcPoints, &
0.0_pREAL,params%mechBC%componentBC(field)%Value(face),params%Delta_t)
0.0_pREAL,params%mechBC(face)%Value(field),params%Delta_t)
call ISDestroy(bcPoints,err_PETSc)
CHKERRQ(err_PETSc)
end if
@ -556,14 +556,14 @@ subroutine FEM_mechanical_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,err_P
CHKERRQ(err_PETSc)
call VecWAXPY(x_local,1.0_pREAL,xx_local,solution_local,err_PETSc)
CHKERRQ(err_PETSc)
do field = 1, dimPlex; do face = 1, mesh_Nboundaries
if (params%mechBC%componentBC(field)%Mask(face)) then
do face = 1, mesh_Nboundaries; do field = 1, dimPlex
if (params%mechBC(face)%Mask(field)) then
call DMGetStratumSize(dm_local,'Face Sets',mesh_boundaries(face),bcSize,err_PETSc)
if (bcSize > 0) then
call DMGetStratumIS(dm_local,'Face Sets',mesh_boundaries(face),bcPoints,err_PETSc)
CHKERRQ(err_PETSc)
call utilities_projectBCValues(x_local,section,0_pPETSCINT,field-1,bcPoints, &
0.0_pREAL,params%mechBC%componentBC(field)%Value(face),params%Delta_t)
0.0_pREAL,params%mechBC(face)%Value(field),params%Delta_t)
call ISDestroy(bcPoints,err_PETSc)
CHKERRQ(err_PETSc)
end if
@ -667,7 +667,7 @@ end subroutine FEM_mechanical_formJacobian
!--------------------------------------------------------------------------------------------------
subroutine FEM_mechanical_forward(guess,Delta_t,Delta_t_prev,mechBC)
type(tMechBC), intent(in) :: &
type(tMechBC), dimension(:), intent(in) :: &
mechBC
real(pREAL), intent(in) :: &
Delta_t_prev, &
@ -701,14 +701,14 @@ subroutine FEM_mechanical_forward(guess,Delta_t,Delta_t_prev,mechBC)
CHKERRQ(err_PETSc)
call VecAXPY(solution_local,1.0_pREAL,x_local,err_PETSc)
CHKERRQ(err_PETSc)
do field = 1, dimPlex; do face = 1, mesh_Nboundaries
if (mechBC%componentBC(field)%Mask(face)) then
do face = 1, mesh_Nboundaries; do field = 1, dimPlex
if (mechBC(face)%Mask(field)) then
call DMGetStratumSize(dm_local,'Face Sets',mesh_boundaries(face),bcSize,err_PETSc)
if (bcSize > 0) then
call DMGetStratumIS(dm_local,'Face Sets',mesh_boundaries(face),bcPoints,err_PETSc)
CHKERRQ(err_PETSc)
call utilities_projectBCValues(solution_local,section,0_pPETSCINT,field-1,bcPoints, &
0.0_pREAL,mechBC%componentBC(field)%Value(face),Delta_t_prev)
0.0_pREAL,mechBC(face)%Value(field),Delta_t_prev)
call ISDestroy(bcPoints,err_PETSc)
CHKERRQ(err_PETSc)
end if