1c3871b17SNoam T! Use DMPlexGetClosureIndices to check for shared node DOF 2c3871b17SNoam T! The mesh consists of two tetrehadra, sharing a (triangular) face, hence 3c3871b17SNoam T! the number of shared DOF equals 3 nodes x 3 dof/node = 9 4c3871b17SNoam T! 5c3871b17SNoam T#include <petsc/finclude/petscdmplex.h> 6c3871b17SNoam T#include <petsc/finclude/petscdm.h> 7c3871b17SNoam T#include <petsc/finclude/petscsys.h> 8c3871b17SNoam T#include <petsc/finclude/petsc.h> 9c5e229c2SMartin Diehlprogram main 10c3871b17SNoam T 11*02c639afSMartin Diehl use petscdm 12*02c639afSMartin Diehl use petscdmplex 13c3871b17SNoam T 14c3871b17SNoam T implicit none 15c3871b17SNoam T 16c3871b17SNoam T DM :: dm, cdm 17c3871b17SNoam T PetscInt :: cStart, cEnd 18c3871b17SNoam T PetscInt :: cdim, nIdx, idx, cnt, Nf 19c3871b17SNoam T PetscInt, parameter :: sharedNodes = 3, zero = 0 20c3871b17SNoam T PetscSection :: gS 21c3871b17SNoam T PetscErrorCode :: ierr 22c3871b17SNoam T 23c3871b17SNoam T PetscInt, allocatable :: idxMatrix(:, :), offsets(:) 24c3871b17SNoam T PetscInt, pointer, dimension(:) :: indices 25c3871b17SNoam T 26c3871b17SNoam T PetscCallA(PetscInitialize(ierr)) 27c3871b17SNoam T 28c3871b17SNoam T PetscCallA(DMCreate(PETSC_COMM_WORLD, dm, ierr)) 29c3871b17SNoam T PetscCallA(DMSetType(dm, DMPLEX, ierr)) 30c3871b17SNoam T PetscCallA(DMSetFromOptions(dm, ierr)) 31c3871b17SNoam T 32c3871b17SNoam T PetscCallA(DMGetCoordinateDM(dm, cdm, ierr)) 33c3871b17SNoam T PetscCallA(DMGetCoordinateDim(cdm, cdim, ierr)) 34c3871b17SNoam T PetscCallA(DMGetGlobalSection(cdm, gS, ierr)) 35c3871b17SNoam T 36c3871b17SNoam T PetscCallA(DMPlexGetHeightStratum(dm, zero, cStart, cEnd, ierr)) 37c3871b17SNoam T PetscCallA(PetscSectionGetNumFields(gS, Nf, ierr)) 38c3871b17SNoam T allocate (offsets(Nf + 1), source=zero) 39c3871b17SNoam T 40c3871b17SNoam T ! Indices per cell 41c3871b17SNoam T ! cell 0 (cStart) 42c3871b17SNoam T PetscCallA(DMPlexGetClosureIndices(cdm, gS, gS, cStart, PETSC_TRUE, nIdx, indices, offsets, PETSC_NULL_SCALAR_POINTER, ierr)) 43c3871b17SNoam T allocate (idxMatrix(nIdx, cEnd - cStart)) 44c3871b17SNoam T idxMatrix(1:nIdx, cStart + 1) = indices 45c3871b17SNoam T ! Check size and content of output field offsets array 46c3871b17SNoam T PetscCheckA(size(offsets) == (Nf + 1) .and. offsets(1) == zero, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Wrong field offsets") 47c3871b17SNoam T PetscCallA(DMPlexRestoreClosureIndices(cdm, gS, gS, cStart, PETSC_TRUE, nIdx, indices, offsets, PETSC_NULL_SCALAR_POINTER, ierr)) 48c3871b17SNoam T ! cell 1 (cEnd - 1) 49c3871b17SNoam T PetscCallA(DMPlexGetClosureIndices(cdm, gS, gS, cEnd - 1, PETSC_TRUE, nIdx, indices, offsets, PETSC_NULL_SCALAR_POINTER, ierr)) 50c3871b17SNoam T idxMatrix(1:nIdx, cEnd - Cstart) = indices 51c3871b17SNoam T PetscCallA(DMPlexRestoreClosureIndices(cdm, gS, gS, cEnd - 1, PETSC_TRUE, nIdx, indices, offsets, PETSC_NULL_SCALAR_POINTER, ierr)) 52c3871b17SNoam T 53c3871b17SNoam T ! Check number of shared indices between cell 0 and cell 1 54c3871b17SNoam T cnt = 0 55c3871b17SNoam T do idx = 1, nIdx 56c3871b17SNoam T cnt = cnt + count(idxMatrix(idx, 1) == idxMatrix(1:nIdx, cEnd)) 57c3871b17SNoam T end do 58c3871b17SNoam T PetscCheckA(cnt == sharedNodes*cdim, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Wrong DOF indices") 59c3871b17SNoam T 60c3871b17SNoam T ! Cleanup 61c3871b17SNoam T deallocate (offsets) 62c3871b17SNoam T deallocate (idxMatrix) 63c3871b17SNoam T PetscCallA(DMDestroy(dm, ierr)) 64c3871b17SNoam T PetscCallA(PetscFinalize(ierr)) 65c3871b17SNoam T 66c3871b17SNoam Tend program main 67c3871b17SNoam T 68c3871b17SNoam T! /*TEST 69c3871b17SNoam T! 70c3871b17SNoam T! test: 71c3871b17SNoam T! args : -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/doublet-tet.msh 72c3871b17SNoam T! output_file: output/empty.out 73c3871b17SNoam T! 74c3871b17SNoam T! TEST*/ 75