diff --git a/src/mesh/DAMASK_mesh.f90 b/src/mesh/DAMASK_mesh.f90 index b8eb64d7d..6b6400d62 100644 --- a/src/mesh/DAMASK_mesh.f90 +++ b/src/mesh/DAMASK_mesh.f90 @@ -25,10 +25,9 @@ program DAMASK_mesh type :: tLoadCase real(pREAL) :: t = 0.0_pREAL !< length of increment 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 - integer, allocatable, dimension(:) :: tag - type(tMechBC) :: mechBC + type(tMechBC), allocatable, dimension(:) :: mechBC end type tLoadCase @@ -40,20 +39,16 @@ program DAMASK_mesh t = 0.0_pREAL, & !< elapsed time t_0 = 0.0_pREAL, & !< begin of interval Delta_t = 0.0_pREAL, & !< current time interval - Delta_t_prev = 0.0_pREAL, & !< previous time interval - t_remaining = 0.0_pREAL !< remaining time of current load case + Delta_t_prev = 0.0_pREAL !< previous time interval logical :: & guess, & !< guess along former trajectory stagIterate integer :: & l, & - i, & m, & errorID, & cutBackLevel = 0, & !< cut back level \f$ t = \frac{t_{inc}}{2^l} \f$ stepFraction = 0, & !< fraction of current time interval - currentLoadcase = 0, & !< current load case - currentFace = 0, & inc, & !< current increment in current load case totalIncsCounter = 0, & !< total # of increments statUnit = 0, & !< file unit for statistics output @@ -82,7 +77,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 +118,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 +137,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 +168,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 @@ -230,7 +215,6 @@ program DAMASK_mesh stepFraction = 0 ! fraction scaled by stepFactor**cutLevel subStepLooping: do while (stepFraction < subStepFactor**cutBackLevel) - t_remaining = loadCases(l)%t + t_0 - t t = t + Delta_t ! forward target time stepFraction = stepFraction + 1 ! count step @@ -315,4 +299,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 4b6621b75..220852b89 100644 --- a/src/mesh/FEM_utilities.f90 +++ b/src/mesh/FEM_utilities.f90 @@ -39,14 +39,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 @@ -55,15 +47,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 @@ -74,11 +61,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 5540c6ea2..13ae970ae 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 @@ -112,7 +112,7 @@ subroutine FEM_mechanical_init(mechBC,num_mesh) PetscInt, dimension(:), pointer :: pNumComp, pNumDof, pBcField, pBcPoint PetscInt :: numBC, bcSize, nc, & - field, faceSet, topologDim, nNodalPoints, & + component, faceSet, topologDim, nNodalPoints, & cellStart, cellEnd, cell, basis IS :: bcPoint @@ -208,17 +208,17 @@ 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 component = 1, dimPlex + if (mechBC(faceSet)%Mask(component)) 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 component = 1, dimPlex + if (mechBC(faceSet)%Mask(component)) then 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) call DMGetStratumSize(mechanical_mesh,'Face Sets',mesh_boundaries(faceSet),bcSize,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 @@ -380,7 +380,7 @@ subroutine FEM_mechanical_formResidual(dm_local,xx_local,f_local,dummy,err_PETSc real(pREAL), dimension(cellDof), target :: f_scal PetscReal :: IcellJMat(dimPlex,dimPlex) PetscReal, dimension(:),pointer :: pV0, pCellJ, pInvcellJ, basisField, basisFieldDer - PetscInt :: cellStart, cellEnd, cell, field, face, & + PetscInt :: cellStart, cellEnd, cell, component, face, & qPt, basis, comp, cidx, & numFields, & bcSize,m,i @@ -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 component = 1_pPETSCINT, dimPlex + if (params%mechBC(face)%Mask(component)) 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) + call utilities_projectBCValues(x_local,section,0_pPETSCINT,component-1,bcPoints, & + 0.0_pREAL,params%mechBC(face)%Value(component),params%Delta_t) call ISDestroy(bcPoints,err_PETSc) CHKERRQ(err_PETSc) 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) :: K_eA, K_eB - PetscInt :: cellStart, cellEnd, cell, field, face, & + PetscInt :: cellStart, cellEnd, cell, component, face, & qPt, basis, comp, cidx,bcSize, m, i IS :: bcPoints @@ -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 component = 1, dimPlex + if (params%mechBC(face)%Mask(component)) 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) + call utilities_projectBCValues(x_local,section,0_pPETSCINT,component-1,bcPoints, & + 0.0_pREAL,params%mechBC(face)%Value(component),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, & @@ -675,7 +675,7 @@ subroutine FEM_mechanical_forward(guess,Delta_t,Delta_t_prev,mechBC) logical, intent(in) :: & guess - PetscInt :: field, face, bcSize + PetscInt :: component, face, bcSize DM :: dm_local Vec :: x_local PetscSection :: section @@ -700,14 +700,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 component = 1, dimPlex + if (mechBC(face)%Mask(component)) 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) + call utilities_projectBCValues(solution_local,section,0_pPETSCINT,component-1,bcPoints, & + 0.0_pREAL,mechBC(face)%Value(component),Delta_t_prev) call ISDestroy(bcPoints,err_PETSc) CHKERRQ(err_PETSc) end if