1: !   Use DMPlexGetClosureIndices to check for shared node DOF
  2: !   The mesh consists of two tetrehadra, sharing a (triangular) face, hence
  3: !   the number of shared DOF equals 3 nodes x 3 dof/node = 9
  4: !
  5: program main
  6: #include <petsc/finclude/petscdmplex.h>
  7: #include <petsc/finclude/petscdm.h>
  8: #include <petsc/finclude/petscsys.h>
  9: #include <petsc/finclude/petsc.h>

 11:   use PETScDM
 12:   use PETScDMplex

 14:   implicit none

 16:   DM :: dm, cdm
 17:   PetscInt :: cStart, cEnd
 18:   PetscInt :: cdim, nIdx, idx, cnt, Nf
 19:   PetscInt, parameter :: sharedNodes = 3, zero = 0
 20:   PetscSection :: gS
 21:   PetscErrorCode :: ierr

 23:   PetscInt, allocatable :: idxMatrix(:, :), offsets(:)
 24:   PetscInt, pointer, dimension(:) :: indices

 26:   PetscCallA(PetscInitialize(ierr))

 28:   PetscCallA(DMCreate(PETSC_COMM_WORLD, dm, ierr))
 29:   PetscCallA(DMSetType(dm, DMPLEX, ierr))
 30:   PetscCallA(DMSetFromOptions(dm, ierr))

 32:   PetscCallA(DMGetCoordinateDM(dm, cdm, ierr))
 33:   PetscCallA(DMGetCoordinateDim(cdm, cdim, ierr))
 34:   PetscCallA(DMGetGlobalSection(cdm, gS, ierr))

 36:   PetscCallA(DMPlexGetHeightStratum(dm, zero, cStart, cEnd, ierr))
 37:   PetscCallA(PetscSectionGetNumFields(gS, Nf, ierr))
 38:   allocate (offsets(Nf + 1), source=zero)

 40:   ! Indices per cell
 41:   ! cell 0 (cStart)
 42:   PetscCallA(DMPlexGetClosureIndices(cdm, gS, gS, cStart, PETSC_TRUE, nIdx, indices, offsets, PETSC_NULL_SCALAR_POINTER, ierr))
 43:   allocate (idxMatrix(nIdx, cEnd - cStart))
 44:   idxMatrix(1:nIdx, cStart + 1) = indices
 45:   ! Check size and content of output field offsets array
 46:   PetscCheckA(size(offsets) == (Nf + 1) .and. offsets(1) == zero, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Wrong field offsets")
 47:   PetscCallA(DMPlexRestoreClosureIndices(cdm, gS, gS, cStart, PETSC_TRUE, nIdx, indices, offsets, PETSC_NULL_SCALAR_POINTER, ierr))
 48:   ! cell 1 (cEnd - 1)
 49:   PetscCallA(DMPlexGetClosureIndices(cdm, gS, gS, cEnd - 1, PETSC_TRUE, nIdx, indices, offsets, PETSC_NULL_SCALAR_POINTER, ierr))
 50:   idxMatrix(1:nIdx, cEnd - Cstart) = indices
 51:   PetscCallA(DMPlexRestoreClosureIndices(cdm, gS, gS, cEnd - 1, PETSC_TRUE, nIdx, indices, offsets, PETSC_NULL_SCALAR_POINTER, ierr))

 53:   ! Check number of shared indices between cell 0 and cell 1
 54:   cnt = 0
 55:   do idx = 1, nIdx
 56:     cnt = cnt + count(idxMatrix(idx, 1) == idxMatrix(1:nIdx, cEnd))
 57:   end do
 58:   PetscCheckA(cnt == sharedNodes*cdim, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Wrong DOF indices")

 60:   ! Cleanup
 61:   deallocate (offsets)
 62:   deallocate (idxMatrix)
 63:   PetscCallA(DMDestroy(dm, ierr))
 64:   PetscCallA(PetscFinalize(ierr))

 66: end program main

 68: ! /*TEST
 69: !
 70: ! test:
 71: !   args : -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/doublet-tet.msh
 72: !   output_file: output/empty.out
 73: !
 74: ! TEST*/