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