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