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:       use petscsys
  6:       use petscdmplex
  7: #include <petsc/finclude/petscsys.h>
  8: #include <petsc/finclude/petscdmplex.h>
  9:       implicit none
 10:       DM :: dm, dmi
 11:       PetscFV :: fvm
 12:       PetscInt, parameter :: dim = 3

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

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

 35:       PetscInt :: i

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

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

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

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

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

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

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

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

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

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