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, parameter :: depth = 1
 13:   PetscErrorCode :: ierr
 14:   PetscInt, dimension(2) :: numPoints = [12, 2]
 15:   PetscInt, dimension(14) :: coneSize = [8, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
 16:   PetscInt, dimension(16) :: cones = [2, 5, 4, 3, 6, 7, 8, 9, 3, 4, 11, 10, 7, 12, 13, 8]
 17:   PetscInt, dimension(16) :: coneOrientations = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
 18:   PetscScalar, dimension(36), parameter :: vertexCoords = [ &
 19:                                                           -0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, -0.5, 1.0, 0.0, &
 20:                                                           -0.5, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, 1.0, -0.5, 1.0, 1.0, &
 21:                                                           0.5, 0.0, 0.0, 0.5, 1.0, 0.0, 0.5, 0.0, 1.0, 0.5, 1.0, 1.0]
 22:   PetscReal :: vol
 23:   PetscReal, target, dimension(dim)  :: centroid
 24:   PetscReal, target, dimension(dim)  :: normal
 25:   PetscReal, pointer :: pcentroid(:)
 26:   PetscReal, pointer :: pnormal(:)

 28:   PetscReal, target, dimension(dim)  :: v0
 29:   PetscReal, target, dimension(dim*dim)  :: J
 30:   PetscReal, target, dimension(dim*dim)  :: invJ
 31:   PetscReal, pointer :: pv0(:)
 32:   PetscReal, pointer :: pJ(:)
 33:   PetscReal, pointer :: pinvJ(:)
 34:   PetscReal :: detJ

 36:   PetscInt :: i

 38:   pcentroid => centroid
 39:   pnormal => normal

 41:   pv0 => v0
 42:   pJ => J
 43:   pinvJ => invJ

 45:   PetscCallA(PetscInitialize(ierr))
 46:   PetscCallA(DMPlexCreate(PETSC_COMM_WORLD, dm, ierr))
 47:   PetscCallA(PetscObjectSetName(dm, 'testplex', ierr))
 48:   PetscCallA(DMSetDimension(dm, dim, ierr))

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

 52:   PetscCallA(DMPlexInterpolate(dm, dmi, ierr))
 53:   PetscCallA(DMPlexCopyCoordinates(dm, dmi, ierr))
 54:   PetscCallA(DMDestroy(dm, ierr))
 55:   dm = dmi

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

 59:   do i = 0, 1
 60:     PetscCallA(DMPlexComputeCellGeometryFVM(dm, i, vol, pcentroid, pnormal, ierr))
 61:     write (*, '(a, i2, a, f8.4, a, 3(f8.4, 1x))') 'cell: ', i, ' volume: ', vol, ' centroid: ', pcentroid(1), pcentroid(2), pcentroid(3)
 62:     PetscCallA(DMPlexComputeCellGeometryAffineFEM(dm, i, pv0, pJ, pinvJ, detJ, ierr))
 63:   end do

 65:   PetscCallA(PetscFVCreate(PETSC_COMM_WORLD, fvm, ierr))
 66:   PetscCallA(PetscFVSetUp(fvm, ierr))
 67:   PetscCallA(PetscFVDestroy(fvm, ierr))

 69:   PetscCallA(DMDestroy(dm, ierr))
 70:   PetscCallA(PetscFinalize(ierr))
 71: end program main

 73: !/*TEST
 74: !
 75: !   test:
 76: !      suffix: 0
 77: !
 78: !TEST*/