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: #include <petsc/finclude/petscdmplex.h>
  6: #include <petsc/finclude/petscdm.h>
  7: #include <petsc/finclude/petscsys.h>
  8: #include <petsc/finclude/petsc.h>
  9: program main

 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
 20:   PetscSection :: gS
 21:   PetscErrorCode :: ierr
 22:   PetscInt, allocatable :: idxMatrix(:, :), offsets(:)
 23:   PetscInt, pointer, dimension(:) :: indices

 25:   PetscCallA(PetscInitialize(ierr))

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

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

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

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

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

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

 65: end program main

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