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*/