From 843676fb10a3291d7817cf021fe4f43225a9d309 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 21 Jan 2020 06:10:19 +0100 Subject: [PATCH] grouping variables for better readability --- src/mesh/FEM_mech.f90 | 107 ++++++++++++++++++++----------------- src/mesh/FEM_utilities.f90 | 16 +++--- src/mesh/FEM_zoo.f90 | 4 +- 3 files changed, 69 insertions(+), 58 deletions(-) diff --git a/src/mesh/FEM_mech.f90 b/src/mesh/FEM_mech.f90 index 3a19f67c7..e56cb7151 100644 --- a/src/mesh/FEM_mech.f90 +++ b/src/mesh/FEM_mech.f90 @@ -67,28 +67,32 @@ contains subroutine FEM_mech_init(fieldBC) type(tFieldBC), intent(in) :: fieldBC + DM :: mech_mesh PetscFE :: mechFE PetscQuadrature :: mechQuad, functional PetscDS :: mechDS PetscDualSpace :: mechDualSpace DMLabel :: BCLabel + 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, pointer :: pBcComps(:), pBcPoints(:) + IS, dimension(:), pointer :: pBcComps, pBcPoints PetscSection :: section - PetscInt :: field, faceSet, topologDim, nNodalPoints + PetscReal, dimension(:), pointer :: qPointsP, qWeightsP, & - nodalPointsP, nodalWeightsP - PetscReal, allocatable, target :: nodalPoints(:), nodalWeights(:) - PetscScalar, pointer :: px_scal(:) - PetscScalar, allocatable, target :: x_scal(:) + nodalPointsP, nodalWeightsP,pV0, pCellJ, pInvcellJ PetscReal :: detJ PetscReal, allocatable, target :: cellJMat(:,:) - PetscReal, pointer :: pV0(:), pCellJ(:), pInvcellJ(:) - PetscInt :: cellStart, cellEnd, cell, basis - character(len=7), parameter :: prefix = 'mechFE_' + + PetscScalar, pointer :: px_scal(:) + PetscScalar, allocatable, target :: x_scal(:) + + character(len=*), parameter :: prefix = 'mechFE_' PetscErrorCode :: ierr 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_rate ,0.0,ierr); CHKERRQ(ierr) allocate(x_scal(cellDof)) - allocate(nodalPoints (dimPlex)) - allocate(nodalWeights(1)) - nodalPointsP => nodalPoints - nodalWeightsP => nodalWeights + allocate(nodalWeightsP(1)) + allocate(nodalPointsP(dimPlex)) allocate(pv0(dimPlex)) - allocate(pcellJ(dimPlex*dimPlex)) - allocate(pinvcellJ(dimPlex*dimPlex)) + allocate(pcellJ(dimPlex**2)) + allocate(pinvcellJ(dimPlex**2)) allocate(cellJMat(dimPlex,dimPlex)) call DMGetSection(mech_mesh,section,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) CHKERRQ(ierr) 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) CHKERRQ(ierr) cellJMat = reshape(pCellJ,shape=[dimPlex,dimPlex]) @@ -221,7 +223,7 @@ subroutine FEM_mech_init(fieldBC) CHKERRQ(ierr) call PetscQuadratureGetData(functional,dimPlex,nc,nNodalPoints,nodalPointsP,nodalWeightsP,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 px_scal => x_scal 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) DM :: dm_local + PetscObject,intent(in) :: dummy + PetscErrorCode :: ierr + PetscDS :: prob Vec :: x_local, f_local, xx_local PetscSection :: section @@ -294,10 +299,10 @@ subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr) qPt, basis, comp, cidx PetscReal :: detFAvg PetscReal :: BMat(dimPlex*dimPlex,cellDof) - PetscObject,intent(in) :: dummy + PetscInt :: bcSize IS :: bcPoints - PetscErrorCode :: ierr + allocate(pV0(dimPlex)) 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) CHKERRQ(ierr) 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) endif endif @@ -403,29 +408,36 @@ end subroutine FEM_mech_formResidual !-------------------------------------------------------------------------------------------------- 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 Vec :: x_local, xx_local - Mat :: Jac_pre, Jac + 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, dimension(:), pointer :: basisField, basisFieldDer, & 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, & K_eA , & K_eB - PetscScalar, target :: K_eVec(cellDof*cellDof) - PetscReal :: BMat (dimPlex*dimPlex,cellDof), & - BMatAvg(dimPlex*dimPlex,cellDof), & - MatA (dimPlex*dimPlex,cellDof), & - MatB (1 ,cellDof) - PetscScalar, dimension(:), pointer :: pK_e, x_scal - PetscReal, dimension(3,3) :: F, FAvg, FInv - PetscObject, intent(in) :: dummy + PetscScalar,(cellDof**2) ,target :: K_eVec + + PetscInt :: cellStart, cellEnd, cell, field, face, & + qPt, basis, comp, cidx,bcSize + IS :: bcPoints - PetscErrorCode :: ierr + allocate(pV0(dimPlex)) 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 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 if (params%fieldBC%componentBC(field)%Mask(face)) then 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) CHKERRQ(ierr) 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) endif endif @@ -501,7 +513,6 @@ subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr) (matmul(matmul(transpose(BMatAvg), & reshape(FInv(1:dimPlex,1:dimPlex),shape=[dimPlex*dimPlex,1],order=[2,1])),MatB) + & K_eB)/real(dimPlex) - else K_e = K_eA endif @@ -536,18 +547,18 @@ subroutine FEM_mech_forward(guess,timeinc,timeinc_old,fieldBC) type(tFieldBC), intent(in) :: & fieldBC - real(pReal), intent(in) :: & + real(pReal), intent(in) :: & timeinc_old, & timeinc - logical, intent(in) :: & + logical, intent(in) :: & guess - PetscInt :: field, face - DM :: dm_local - Vec :: x_local - PetscSection :: section - PetscInt :: bcSize - IS :: bcPoints - PetscErrorCode :: ierr + + PetscInt :: field, face, bcSize + DM :: dm_local + Vec :: x_local + PetscSection :: section + IS :: bcPoints + PetscErrorCode :: ierr !-------------------------------------------------------------------------------------------------- ! 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 DMGetSection(dm_local,section,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 CHKERRQ(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) CHKERRQ(ierr) 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) endif endif diff --git a/src/mesh/FEM_utilities.f90 b/src/mesh/FEM_utilities.f90 index 1303f0df8..34cdd7214 100644 --- a/src/mesh/FEM_utilities.f90 +++ b/src/mesh/FEM_utilities.f90 @@ -23,7 +23,7 @@ module FEM_utilities 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 :: nActiveFields = 0 @@ -56,15 +56,15 @@ module FEM_utilities !-------------------------------------------------------------------------------------------------- ! derived types type, public :: tSolutionState !< return type of solution from FEM solver variants - logical :: converged = .true. - logical :: stagConverged = .true. - integer :: iterationsNeeded = 0 + logical :: converged = .true. + logical :: stagConverged = .true. + integer :: iterationsNeeded = 0 end type tSolutionState type, public :: tComponentBC integer(kind(COMPONENT_UNDEFINED_ID)) :: ID - real(pReal), allocatable :: Value(:) - logical, allocatable :: Mask(:) + real(pReal), allocatable, dimension(:) :: Value + logical, allocatable, dimension(:) :: Mask end type tComponentBC type, public :: tFieldBC @@ -79,8 +79,8 @@ module FEM_utilities outputfrequency = 1, & !< frequency of result writes logscale = 0 !< linear/logarithmic time inc flag logical :: followFormerTrajectory = .true. !< follow trajectory of former loadcase - integer, allocatable :: faceID(:) - type(tFieldBC), allocatable :: fieldBC(:) + integer, allocatable, dimension(:) :: faceID + type(tFieldBC), allocatable, dimension(:) :: fieldBC end type tLoadCase public :: & diff --git a/src/mesh/FEM_zoo.f90 b/src/mesh/FEM_zoo.f90 index 5ef37e462..5ab1c3102 100644 --- a/src/mesh/FEM_zoo.f90 +++ b/src/mesh/FEM_zoo.f90 @@ -37,7 +37,7 @@ contains !-------------------------------------------------------------------------------------------------- subroutine FEM_Zoo_init - write(6,'(/,a)') ' <<<+- FEM_Zoo init -+>>>' + write(6,'(/,a)') ' <<<+- FEM_Zoo init -+>>>'; flush(6) !-------------------------------------------------------------------------------------------------- ! 2D linear @@ -53,7 +53,7 @@ subroutine FEM_Zoo_init ! 2D quadratic 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 allocate(FEM_Zoo_QuadraturePoints (2,2)%p(6))