This commit is contained in:
Martin Diehl 2020-01-20 23:46:48 +01:00
parent 80ae2f5f6b
commit db96ee0fc2
1 changed files with 9 additions and 14 deletions

View File

@ -197,21 +197,17 @@ end subroutine mesh_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine mesh_FEM_build_ipVolumes(dimPlex) subroutine mesh_FEM_build_ipVolumes(dimPlex)
PetscInt :: dimPlex PetscInt,intent(in):: dimPlex
PetscReal :: vol PetscReal :: vol
PetscReal, target :: cent(dimPlex), norm(dimPlex) PetscReal, pointer,dimension(:) :: pCent, pNorm
PetscReal, pointer :: pCent(:), pNorm(:)
PetscInt :: cellStart, cellEnd, cell PetscInt :: cellStart, cellEnd, cell
PetscErrorCode :: ierr PetscErrorCode :: ierr
if (.not. allocated(mesh_ipVolume)) then allocate(mesh_ipVolume(mesh_maxNips,mesh_NcpElems),source=0.0_pReal)
allocate(mesh_ipVolume(mesh_maxNips,mesh_NcpElems))
mesh_ipVolume = 0.0_pReal
endif
call DMPlexGetHeightStratum(geomMesh,0,cellStart,cellEnd,ierr); CHKERRQ(ierr) call DMPlexGetHeightStratum(geomMesh,0,cellStart,cellEnd,ierr); CHKERRQ(ierr)
pCent => cent allocate(pCent(dimPlex))
pNorm => norm allocate(pNorm(dimPlex))
do cell = cellStart, cellEnd-1 do cell = cellStart, cellEnd-1
call DMPlexComputeCellGeometryFVM(geomMesh,cell,vol,pCent,pNorm,ierr) call DMPlexComputeCellGeometryFVM(geomMesh,cell,vol,pCent,pNorm,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
@ -229,8 +225,7 @@ subroutine mesh_FEM_build_ipCoordinates(dimPlex,qPoints)
PetscInt, intent(in) :: dimPlex PetscInt, intent(in) :: dimPlex
PetscReal, intent(in) :: qPoints(mesh_maxNips*dimPlex) PetscReal, intent(in) :: qPoints(mesh_maxNips*dimPlex)
PetscReal, target :: v0(dimPlex), cellJ(dimPlex*dimPlex), invcellJ(dimPlex*dimPlex) PetscReal, pointer,dimension(:) :: pV0, pCellJ, pInvcellJ
PetscReal, pointer :: pV0(:), pCellJ(:), pInvcellJ(:)
PetscReal :: detJ PetscReal :: detJ
PetscInt :: cellStart, cellEnd, cell, qPt, dirI, dirJ, qOffset PetscInt :: cellStart, cellEnd, cell, qPt, dirI, dirJ, qOffset
PetscErrorCode :: ierr PetscErrorCode :: ierr
@ -238,9 +233,9 @@ subroutine mesh_FEM_build_ipCoordinates(dimPlex,qPoints)
allocate(mesh_ipCoordinates(3,mesh_maxNips,mesh_NcpElems),source=0.0_pReal) allocate(mesh_ipCoordinates(3,mesh_maxNips,mesh_NcpElems),source=0.0_pReal)
pV0 => v0 allocate(pV0(dimPlex))
pCellJ => cellJ allocatE(pCellJ(dimPlex**2))
pInvcellJ => invcellJ allocatE(pinvCellJ(dimPlex**2))
call DMPlexGetHeightStratum(geomMesh,0,cellStart,cellEnd,ierr); CHKERRQ(ierr) call DMPlexGetHeightStratum(geomMesh,0,cellStart,cellEnd,ierr); CHKERRQ(ierr)
do cell = cellStart, cellEnd-1 !< loop over all elements do cell = cellStart, cellEnd-1 !< loop over all elements
call DMPlexComputeCellGeometryAffineFEM(geomMesh,cell,pV0,pCellJ,pInvcellJ,detJ,ierr) call DMPlexComputeCellGeometryAffineFEM(geomMesh,cell,pV0,pCellJ,pInvcellJ,detJ,ierr)