1c4762a1bSJed Brown program main 2c4762a1bSJed Brown#include <petsc/finclude/petscdmplex.h> 3c4762a1bSJed Brown use petscdmplex 4c4762a1bSJed Brown use petscsys 5c4762a1bSJed Brown implicit none 6c4762a1bSJed Brown 7c4762a1bSJed Brown DM dm 8c4762a1bSJed Brown PetscInt, target, dimension(3) :: EC 9c4762a1bSJed Brown PetscInt, target, dimension(2) :: VE 10c4762a1bSJed Brown PetscInt, pointer :: pEC(:), pVE(:) 11c4762a1bSJed Brown PetscInt, pointer :: nClosure(:) 12c4762a1bSJed Brown PetscInt, pointer :: nJoin(:) 13c4762a1bSJed Brown PetscInt, pointer :: nMeet(:) 14c4762a1bSJed Brown PetscInt dim, cell, size 15*d33816bfSBarry Smith PetscInt i0,i1,i2,i3,i6,i7 16c4762a1bSJed Brown PetscInt i8,i9,i10,i11 17c4762a1bSJed Brown PetscErrorCode ierr 18c4762a1bSJed Brown 19c4762a1bSJed Brown i0 = 0 20c4762a1bSJed Brown i1 = 1 21c4762a1bSJed Brown i2 = 2 22c4762a1bSJed Brown i3 = 3 23c4762a1bSJed Brown i6 = 6 24c4762a1bSJed Brown i7 = 7 25c4762a1bSJed Brown i8 = 8 26c4762a1bSJed Brown i9 = 9 27c4762a1bSJed Brown i10 = 10 28c4762a1bSJed Brown i11 = 11 29c4762a1bSJed Brown 30c4762a1bSJed Brown call PetscInitialize(PETSC_NULL_CHARACTER,ierr) 31c4762a1bSJed Brown if (ierr .ne. 0) then 32c4762a1bSJed Brown print*,'Unable to initialize PETSc' 33c4762a1bSJed Brown stop 34c4762a1bSJed Brown endif 35c4762a1bSJed Brown 36c4762a1bSJed Brown call DMPlexCreate(PETSC_COMM_WORLD, dm, ierr);CHKERRA(ierr) 37c4762a1bSJed Brown call PetscObjectSetName(dm, 'Mesh', ierr);CHKERRA(ierr) 38c4762a1bSJed Brown dim = 2 39c4762a1bSJed Brown call DMSetDimension(dm, dim, ierr);CHKERRA(ierr) 40c4762a1bSJed Brown 41c4762a1bSJed Brown! Make Doublet Mesh from Fig 2 of Flexible Representation of Computational Meshes, 42c4762a1bSJed Brown! except indexing is from 0 instead of 1 and we obey the new restrictions on 43c4762a1bSJed Brown! numbering: cells, vertices, faces, edges 44c4762a1bSJed Brown call DMPlexSetChart(dm, i0, i11, ierr);CHKERRA(ierr) 45c4762a1bSJed Brown! cells 46c4762a1bSJed Brown call DMPlexSetConeSize(dm, i0, i3, ierr);CHKERRA(ierr) 47c4762a1bSJed Brown call DMPlexSetConeSize(dm, i1, i3, ierr);CHKERRA(ierr) 48c4762a1bSJed Brown! edges 49c4762a1bSJed Brown call DMPlexSetConeSize(dm, i6, i2, ierr);CHKERRA(ierr) 50c4762a1bSJed Brown call DMPlexSetConeSize(dm, i7, i2, ierr);CHKERRA(ierr) 51c4762a1bSJed Brown call DMPlexSetConeSize(dm, i8, i2, ierr);CHKERRA(ierr) 52c4762a1bSJed Brown call DMPlexSetConeSize(dm, i9, i2, ierr);CHKERRA(ierr) 53c4762a1bSJed Brown call DMPlexSetConeSize(dm, i10, i2, ierr);CHKERRA(ierr) 54c4762a1bSJed Brown 55c4762a1bSJed Brown call DMSetUp(dm, ierr);CHKERRA(ierr) 56c4762a1bSJed Brown 57c4762a1bSJed Brown EC(1) = 6 58c4762a1bSJed Brown EC(2) = 7 59c4762a1bSJed Brown EC(3) = 8 60c4762a1bSJed Brown pEC => EC 61c4762a1bSJed Brown call DMPlexSetCone(dm, i0, pEC, ierr);CHKERRA(ierr) 62c4762a1bSJed Brown EC(1) = 7 63c4762a1bSJed Brown EC(2) = 9 64c4762a1bSJed Brown EC(3) = 10 65c4762a1bSJed Brown pEC => EC 66c4762a1bSJed Brown call DMPlexSetCone(dm, i1 , pEC, ierr);CHKERRA(ierr) 67c4762a1bSJed Brown 68c4762a1bSJed Brown VE(1) = 2 69c4762a1bSJed Brown VE(2) = 3 70c4762a1bSJed Brown pVE => VE 71c4762a1bSJed Brown call DMPlexSetCone(dm, i6 , pVE, ierr);CHKERRA(ierr) 72c4762a1bSJed Brown VE(1) = 3 73c4762a1bSJed Brown VE(2) = 4 74c4762a1bSJed Brown pVE => VE 75c4762a1bSJed Brown call DMPlexSetCone(dm, i7 , pVE, ierr);CHKERRA(ierr) 76c4762a1bSJed Brown VE(1) = 4 77c4762a1bSJed Brown VE(2) = 2 78c4762a1bSJed Brown pVE => VE 79c4762a1bSJed Brown call DMPlexSetCone(dm, i8 , pVE, ierr);CHKERRA(ierr) 80c4762a1bSJed Brown VE(1) = 3 81c4762a1bSJed Brown VE(2) = 5 82c4762a1bSJed Brown pVE => VE 83c4762a1bSJed Brown call DMPlexSetCone(dm, i9 , pVE, ierr);CHKERRA(ierr) 84c4762a1bSJed Brown VE(1) = 5 85c4762a1bSJed Brown VE(2) = 4 86c4762a1bSJed Brown pVE => VE 87c4762a1bSJed Brown call DMPlexSetCone(dm, i10 , pVE, ierr);CHKERRA(ierr) 88c4762a1bSJed Brown 89c4762a1bSJed Brown call DMPlexSymmetrize(dm,ierr);CHKERRA(ierr) 90c4762a1bSJed Brown call DMPlexStratify(dm,ierr);CHKERRA(ierr) 91c4762a1bSJed Brown call DMView(dm, PETSC_VIEWER_STDOUT_WORLD, ierr);CHKERRA(ierr) 92c4762a1bSJed Brown 93c4762a1bSJed Brown! Test Closure 94c4762a1bSJed Brown do cell = 0,1 95c4762a1bSJed Brown call DMPlexGetTransitiveClosure(dm,cell,PETSC_TRUE,nClosure,ierr);CHKERRA(ierr) 96c4762a1bSJed Brown! Different Fortran compilers print a different number of columns 97c4762a1bSJed Brown! per row producing different outputs in the test runs hence 98c4762a1bSJed Brown! do not print the nClosure 99c4762a1bSJed Brown write(*,1000) 'nClosure ',nClosure 100c4762a1bSJed Brown 1000 format (a,30i4) 101c4762a1bSJed Brown call DMPlexRestoreTransitiveClosure(dm,cell,PETSC_TRUE,nClosure,ierr);CHKERRA(ierr) 102c4762a1bSJed Brown end do 103c4762a1bSJed Brown 104c4762a1bSJed Brown! Test Join 105c4762a1bSJed Brown size = 2 106c4762a1bSJed Brown VE(1) = 6 107c4762a1bSJed Brown VE(2) = 7 108c4762a1bSJed Brown pVE => VE 109c4762a1bSJed Brown call DMPlexGetJoin(dm, size, pVE, nJoin, ierr);CHKERRA(ierr) 110c4762a1bSJed Brown write(*,1001) 'Join of',pVE 111c4762a1bSJed Brown write(*,1002) ' is',nJoin 112c4762a1bSJed Brown call DMPlexRestoreJoin(dm, size, pVE, nJoin, ierr);CHKERRA(ierr) 113c4762a1bSJed Brown size = 2 114c4762a1bSJed Brown VE(1) = 9 115c4762a1bSJed Brown VE(2) = 7 116c4762a1bSJed Brown pVE => VE 117c4762a1bSJed Brown call DMPlexGetJoin(dm, size, pVE, nJoin, ierr);CHKERRA(ierr) 118c4762a1bSJed Brown write(*,1001) 'Join of',pVE 119c4762a1bSJed Brown 1001 format (a,10i5) 120c4762a1bSJed Brown write(*,1002) ' is',nJoin 121c4762a1bSJed Brown 1002 format (a,10i5) 122c4762a1bSJed Brown call DMPlexRestoreJoin(dm, size, pVE, nJoin, ierr);CHKERRA(ierr) 123c4762a1bSJed Brown! Test Full Join 124c4762a1bSJed Brown size = 3 125c4762a1bSJed Brown EC(1) = 3 126c4762a1bSJed Brown EC(2) = 4 127c4762a1bSJed Brown EC(3) = 5 128c4762a1bSJed Brown pEC => EC 129c4762a1bSJed Brown call DMPlexGetFullJoin(dm, size, pEC, nJoin, ierr);CHKERRA(ierr) 130c4762a1bSJed Brown write(*,1001) 'Full Join of',pEC 131c4762a1bSJed Brown write(*,1002) ' is',nJoin 132c4762a1bSJed Brown call DMPlexRestoreJoin(dm, size, pEC, nJoin, ierr);CHKERRA(ierr) 133c4762a1bSJed Brown! Test Meet 134c4762a1bSJed Brown size = 2 135c4762a1bSJed Brown VE(1) = 0 136c4762a1bSJed Brown VE(2) = 1 137c4762a1bSJed Brown pVE => VE 138c4762a1bSJed Brown call DMPlexGetMeet(dm, size, pVE, nMeet, ierr);CHKERRA(ierr) 139c4762a1bSJed Brown write(*,1001) 'Meet of',pVE 140c4762a1bSJed Brown write(*,1002) ' is',nMeet 141c4762a1bSJed Brown call DMPlexRestoreMeet(dm, size, pVE, nMeet, ierr);CHKERRA(ierr) 142c4762a1bSJed Brown size = 2 143c4762a1bSJed Brown VE(1) = 6 144c4762a1bSJed Brown VE(2) = 7 145c4762a1bSJed Brown pVE => VE 146c4762a1bSJed Brown call DMPlexGetMeet(dm, size, pVE, nMeet, ierr);CHKERRA(ierr) 147c4762a1bSJed Brown write(*,1001) 'Meet of',pVE 148c4762a1bSJed Brown write(*,1002) ' is',nMeet 149c4762a1bSJed Brown call DMPlexRestoreMeet(dm, size, pVE, nMeet, ierr);CHKERRA(ierr) 150c4762a1bSJed Brown 151c4762a1bSJed Brown call DMDestroy(dm, ierr);CHKERRA(ierr) 152c4762a1bSJed Brown call PetscFinalize(ierr) 153c4762a1bSJed Brown end 154c4762a1bSJed Brown! 155c4762a1bSJed Brown!/*TEST 156c4762a1bSJed Brown! 157c4762a1bSJed Brown! test: 158c4762a1bSJed Brown! suffix: 0 159c4762a1bSJed Brown! 160c4762a1bSJed Brown!TEST*/ 161