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