diff --git a/src/mesh/DAMASK_mesh.f90 b/src/mesh/DAMASK_mesh.f90 index b8eb64d7d..69688a039 100644 --- a/src/mesh/DAMASK_mesh.f90 +++ b/src/mesh/DAMASK_mesh.f90 @@ -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 diff --git a/src/mesh/FEM_utilities.f90 b/src/mesh/FEM_utilities.f90 index 15a2168a8..cc906a92b 100644 --- a/src/mesh/FEM_utilities.f90 +++ b/src/mesh/FEM_utilities.f90 @@ -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 diff --git a/src/mesh/mesh_mech_FEM.f90 b/src/mesh/mesh_mech_FEM.f90 index 824ebecdd..16cf94ae5 100644 --- a/src/mesh/mesh_mech_FEM.f90 +++ b/src/mesh/mesh_mech_FEM.f90 @@ -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