1: ! setting up 3-D DMPlex using DMPlexCreateFromDAG(), DMPlexInterpolate() and
  2: ! DMPlexComputeCellGeometryFVM()
  3: ! Contributed by Adrian Croucher <a.croucher@auckland.ac.nz>
  4: #include <petsc/finclude/petscdmplex.h>
  5: program main
  6:   use petscdmplex
  7:   implicit none
  8:   DM :: dm, dmi
  9:   PetscFV :: fvm
 10:   PetscInt, parameter :: dim = 3

 12:   PetscInt :: depth = 1
 13:   PetscErrorCode :: ierr
 14:   PetscInt, dimension(2) :: numPoints
 15:   PetscInt, dimension(14) :: coneSize
 16:   PetscInt, dimension(16) :: cones
 17:   PetscInt, dimension(16) :: coneOrientations
 18:   PetscScalar, dimension(36) :: vertexCoords
 19:   PetscReal ::  vol = 0.
 20:   PetscReal, target, dimension(dim)  :: centroid
 21:   PetscReal, target, dimension(dim)  :: normal
 22:   PetscReal, pointer :: pcentroid(:)
 23:   PetscReal, pointer :: pnormal(:)

 25:   PetscReal, target, dimension(dim)  :: v0
 26:   PetscReal, target, dimension(dim*dim)  :: J
 27:   PetscReal, target, dimension(dim*dim)  :: invJ
 28:   PetscReal, pointer :: pv0(:)
 29:   PetscReal, pointer :: pJ(:)
 30:   PetscReal, pointer :: pinvJ(:)
 31:   PetscReal :: detJ

 33:   PetscInt :: i

 35:   numPoints = [12, 2]
 36:   coneSize = [8, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
 37:   cones = [2, 5, 4, 3, 6, 7, 8, 9, 3, 4, 11, 10, 7, 12, 13, 8]
 38:   coneOrientations = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
 39:   vertexCoords = [                                                  &
 40:  &  -0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, -0.5, 1.0, 0.0,           &
 41:  &  -0.5, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, 1.0, -0.5, 1.0, 1.0,           &
 42:  &   0.5, 0.0, 0.0, 0.5, 1.0, 0.0, 0.5, 0.0, 1.0, 0.5, 1.0, 1.0]
 43:   pcentroid => centroid
 44:   pnormal => normal

 46:   pv0 => v0
 47:   pJ => J
 48:   pinvJ => invJ

 50:   PetscCallA(PetscInitialize(ierr))
 51:   PetscCallA(DMPlexCreate(PETSC_COMM_WORLD, dm, ierr))
 52:   PetscCallA(PetscObjectSetName(dm, 'testplex', ierr))
 53:   PetscCallA(DMSetDimension(dm, dim, ierr))

 55:   PetscCallA(DMPlexCreateFromDAG(dm, depth, numPoints, coneSize, cones, coneOrientations, vertexCoords, ierr))

 57:   PetscCallA(DMPlexInterpolate(dm, dmi, ierr))
 58:   PetscCallA(DMPlexCopyCoordinates(dm, dmi, ierr))
 59:   PetscCallA(DMDestroy(dm, ierr))
 60:   dm = dmi

 62:   PetscCallA(DMView(dm, PETSC_VIEWER_STDOUT_WORLD, ierr))

 64:   do i = 0, 1
 65:     PetscCallA(DMPlexComputeCellGeometryFVM(dm, i, vol, pcentroid, pnormal, ierr))
 66:     write (*, '(a, i2, a, f8.4, a, 3(f8.4, 1x))') 'cell: ', i, ' volume: ', vol, ' centroid: ', pcentroid(1), pcentroid(2), pcentroid(3)
 67:     PetscCallA(DMPlexComputeCellGeometryAffineFEM(dm, i, pv0, pJ, pinvJ, detJ, ierr))
 68:   end do

 70:   PetscCallA(PetscFVCreate(PETSC_COMM_WORLD, fvm, ierr))
 71:   PetscCallA(PetscFVSetUp(fvm, ierr))
 72:   PetscCallA(PetscFVDestroy(fvm, ierr))

 74:   PetscCallA(DMDestroy(dm, ierr))
 75:   PetscCallA(PetscFinalize(ierr))
 76: end program main

 78: !/*TEST
 79: !
 80: !   test:
 81: !      suffix: 0
 82: !
 83: !TEST*/