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> 9program main 10 11 use petscdm 12 use petscdmplex 13 14 implicit none 15 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 22 23 PetscInt, allocatable :: idxMatrix(:, :), offsets(:) 24 PetscInt, pointer, dimension(:) :: indices 25 26 PetscCallA(PetscInitialize(ierr)) 27 28 PetscCallA(DMCreate(PETSC_COMM_WORLD, dm, ierr)) 29 PetscCallA(DMSetType(dm, DMPLEX, ierr)) 30 PetscCallA(DMSetFromOptions(dm, ierr)) 31 32 PetscCallA(DMGetCoordinateDM(dm, cdm, ierr)) 33 PetscCallA(DMGetCoordinateDim(cdm, cdim, ierr)) 34 PetscCallA(DMGetGlobalSection(cdm, gS, ierr)) 35 36 PetscCallA(DMPlexGetHeightStratum(dm, zero, cStart, cEnd, ierr)) 37 PetscCallA(PetscSectionGetNumFields(gS, Nf, ierr)) 38 allocate (offsets(Nf + 1), source=zero) 39 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)) 52 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") 59 60 ! Cleanup 61 deallocate (offsets) 62 deallocate (idxMatrix) 63 PetscCallA(DMDestroy(dm, ierr)) 64 PetscCallA(PetscFinalize(ierr)) 65 66end program main 67 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*/ 75