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

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

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

 34:   PetscInt :: i

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

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

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

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

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

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

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

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

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

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