1: ! setting up DMPlex for finite elements
  2: ! Contributed by Pratheek Shanthraj <p.shanthraj@mpie.de>
  3: program main
  4: #include <petsc/finclude/petsc.h>
  5:   use petsc
  6:   implicit none
  7:   DM :: dm
  8:   PetscDS :: ds
  9:   PetscInt :: dim = 3, zero = 0
 10:   PetscBool :: simplex = PETSC_TRUE
 11:   PetscBool :: interpolate = PETSC_TRUE
 12:   PetscReal :: refinementLimit = 0.0
 13:   PetscErrorCode :: ierr
 14:   PetscTabulation, pointer :: tab(:)
 15:   PetscFE fe, rfe
 16:   PetscObject obj
 17:   PetscInt :: one = 1, mone = -1

 19:   PetscCallA(PetscInitialize(PETSC_NULL_CHARACTER, ierr))
 20:   PetscCallA(DMPlexCreateDoublet(PETSC_COMM_WORLD, dim, simplex, interpolate, refinementLimit, dm, ierr))
 21:   PetscCallA(PetscFECreateDefault(PETSC_COMM_WORLD, dim, one, simplex, 'name', mone, fe, ierr))
 22:   PetscCallA(PetscObjectSetName(fe, 'name', ierr))
 23:   PetscCallA(DMSetField(dm, zero, PETSC_NULL_DMLABEL, PetscObjectCast(fe), ierr))
 24:   PetscCallA(DMSetField(dm, one, PETSC_NULL_DMLABEL, PetscObjectCast(fe), ierr))

 26:   PetscCallA(DMSetUp(dm, ierr))
 27:   PetscCallA(DMCreateDS(dm, ierr))
 28:   PetscCallA(DMGetDS(dm, ds, ierr))
 29:   PetscCallA(PetscDSGetTabulation(ds, tab, ierr))
 30:   print *, tab(1)%ptr%T(1)%ptr
 31:   print *, tab(1)%ptr%T(2)%ptr
 32:   print *, tab(2)%ptr%T(1)%ptr
 33:   print *, tab(2)%ptr%T(2)%ptr
 34:   PetscCallA(PetscDSRestoreTabulation(ds, tab, ierr))

 36:   PetscCallA(PetscDSGetDiscretization(ds, zero, obj, ierr))
 37:   PetscObjectSpecificCast(rfe, obj)
 38:   PetscCallA(PetscFEDestroy(fe, ierr))
 39:   PetscCallA(DMDestroy(dm, ierr))
 40:   PetscCallA(PetscFinalize(ierr))
 41: end program main
 42: !/*TEST
 43: !
 44: !  test:
 45: !    nsize: 1
 46: !
 47: !TEST*/