1: #include <petsc/finclude/petscdmplex.h>
  2: program main
  3:   use petscdm
  4:   use petscdmplex
  5:   implicit none

  7:   DM dm
  8:   PetscInt, target, dimension(3) :: EC
  9:   PetscInt, target, dimension(2) :: VE
 10:   PetscInt, pointer, dimension(:) :: pEC, pVE, nClosure, nJoin, nMeet
 11:   PetscInt, parameter :: dim = 2
 12:   PetscInt cell, size, nC
 13:   PetscErrorCode ierr

 15:   PetscCallA(PetscInitialize(ierr))

 17:   PetscCallA(DMPlexCreate(PETSC_COMM_WORLD, dm, ierr))
 18:   PetscCallA(PetscObjectSetName(dm, 'Mesh', ierr))
 19:   PetscCallA(DMSetDimension(dm, dim, ierr))

 21: ! Make Doublet Mesh from Fig 2 of Flexible Representation of Computational Meshes,
 22: ! except indexing is from 0 instead of 1 and we obey the new restrictions on
 23: ! numbering: cells, vertices, faces, edges
 24:   PetscCallA(DMPlexSetChart(dm, 0_PETSC_INT_KIND, 11_PETSC_INT_KIND, ierr))
 25: !     cells
 26:   PetscCallA(DMPlexSetConeSize(dm, 0_PETSC_INT_KIND, 3_PETSC_INT_KIND, ierr))
 27:   PetscCallA(DMPlexSetConeSize(dm, 1_PETSC_INT_KIND, 3_PETSC_INT_KIND, ierr))
 28: !     edges
 29:   PetscCallA(DMPlexSetConeSize(dm, 6_PETSC_INT_KIND, 2_PETSC_INT_KIND, ierr))
 30:   PetscCallA(DMPlexSetConeSize(dm, 7_PETSC_INT_KIND, 2_PETSC_INT_KIND, ierr))
 31:   PetscCallA(DMPlexSetConeSize(dm, 8_PETSC_INT_KIND, 2_PETSC_INT_KIND, ierr))
 32:   PetscCallA(DMPlexSetConeSize(dm, 9_PETSC_INT_KIND, 2_PETSC_INT_KIND, ierr))
 33:   PetscCallA(DMPlexSetConeSize(dm, 10_PETSC_INT_KIND, 2_PETSC_INT_KIND, ierr))

 35:   PetscCallA(DMSetUp(dm, ierr))

 37:   EC = [6, 7, 8]
 38:   pEC => EC
 39:   PetscCallA(DMPlexSetCone(dm, 0_PETSC_INT_KIND, pEC, ierr))
 40:   EC = [7, 9, 10]
 41:   pEC => EC
 42:   PetscCallA(DMPlexSetCone(dm, 1_PETSC_INT_KIND, pEC, ierr))

 44:   VE = [2, 3]
 45:   pVE => VE
 46:   PetscCallA(DMPlexSetCone(dm, 6_PETSC_INT_KIND, pVE, ierr))
 47:   VE = [3, 4]
 48:   pVE => VE
 49:   PetscCallA(DMPlexSetCone(dm, 7_PETSC_INT_KIND, pVE, ierr))
 50:   VE = [4, 2]
 51:   pVE => VE
 52:   PetscCallA(DMPlexSetCone(dm, 8_PETSC_INT_KIND, pVE, ierr))
 53:   VE = [3, 5]
 54:   pVE => VE
 55:   PetscCallA(DMPlexSetCone(dm, 9_PETSC_INT_KIND, pVE, ierr))
 56:   VE = [5, 4]
 57:   pVE => VE
 58:   PetscCallA(DMPlexSetCone(dm, 10_PETSC_INT_KIND, pVE, ierr))

 60:   PetscCallA(DMPlexSymmetrize(dm, ierr))
 61:   PetscCallA(DMPlexStratify(dm, ierr))
 62:   PetscCallA(DMView(dm, PETSC_VIEWER_STDOUT_WORLD, ierr))

 64: ! Test Closure
 65:   do cell = 0, 1
 66:     PetscCallA(DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, nC, nClosure, ierr))
 67: !   Different Fortran compilers print a different number of columns
 68: !   per row producing different outputs in the test runs hence
 69: !   do not print the nClosure
 70:     write (*, 1000) 'nClosure ', nClosure
 71: 1000 format(a, 30i4)
 72:     PetscCallA(DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, nC, nClosure, ierr))
 73:   end do

 75: ! Test Join
 76:   size = 2
 77:   VE = [6, 7]
 78:   pVE => VE
 79:   PetscCallA(DMPlexGetJoin(dm, size, pVE, PETSC_NULL_INTEGER, nJoin, ierr))
 80:   write (*, 1001) 'Join of', pVE
 81:   write (*, 1002) '  is', nJoin
 82:   PetscCallA(DMPlexRestoreJoin(dm, size, pVE, PETSC_NULL_INTEGER, nJoin, ierr))
 83:   size = 2
 84:   VE = [9, 7]
 85:   pVE => VE
 86:   PetscCallA(DMPlexGetJoin(dm, size, pVE, PETSC_NULL_INTEGER, nJoin, ierr))
 87:   write (*, 1001) 'Join of', pVE
 88: 1001 format(a, 10i5)
 89:   write (*, 1002) '  is', nJoin
 90: 1002 format(a, 10i5)
 91:   PetscCallA(DMPlexRestoreJoin(dm, size, pVE, PETSC_NULL_INTEGER, nJoin, ierr))
 92: ! Test Full Join
 93:   size = 3
 94:   EC = [3, 4, 5]
 95:   pEC => EC
 96:   PetscCallA(DMPlexGetFullJoin(dm, size, pEC, PETSC_NULL_INTEGER, nJoin, ierr))
 97:   write (*, 1001) 'Full Join of', pEC
 98:   write (*, 1002) '  is', nJoin
 99:   PetscCallA(DMPlexRestoreJoin(dm, size, pEC, PETSC_NULL_INTEGER, nJoin, ierr))
100: ! Test Meet
101:   size = 2
102:   VE = [0, 1]
103:   pVE => VE
104:   PetscCallA(DMPlexGetMeet(dm, size, pVE, PETSC_NULL_INTEGER, nMeet, ierr))
105:   write (*, 1001) 'Meet of', pVE
106:   write (*, 1002) '  is', nMeet
107:   PetscCallA(DMPlexRestoreMeet(dm, size, pVE, PETSC_NULL_INTEGER, nMeet, ierr))
108:   size = 2
109:   VE = [6, 7]
110:   pVE => VE
111:   PetscCallA(DMPlexGetMeet(dm, size, pVE, PETSC_NULL_INTEGER, nMeet, ierr))
112:   write (*, 1001) 'Meet of', pVE
113:   write (*, 1002) '  is', nMeet
114:   PetscCallA(DMPlexRestoreMeet(dm, size, pVE, PETSC_NULL_INTEGER, nMeet, ierr))

116:   PetscCallA(DMDestroy(dm, ierr))
117:   PetscCallA(PetscFinalize(ierr))
118: end
119: !
120: !/*TEST
121: !
122: !   test:
123: !     suffix: 0
124: !
125: !TEST*/