Actual source code: ex2f90.F90

  1:       program main
  2: #include <petsc/finclude/petscdmplex.h>
  3:       use petscdmplex
  4:       use petscsys
  5:       implicit none

  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

 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

 30:       PetscCallA(PetscInitialize(ierr))

 32:       PetscCallA(DMPlexCreate(PETSC_COMM_WORLD, dm, ierr))
 33:       PetscCallA(PetscObjectSetName(dm, 'Mesh', ierr))
 34:       dim = 2
 35:       PetscCallA(DMSetDimension(dm, dim, ierr))

 37: ! Make Doublet Mesh from Fig 2 of Flexible Representation of Computational Meshes,
 38: ! except indexing is from 0 instead of 1 and we obey the new restrictions on
 39: ! numbering: cells, vertices, faces, edges
 40:       PetscCallA(DMPlexSetChart(dm, i0, i11, ierr))
 41: !     cells
 42:       PetscCallA(DMPlexSetConeSize(dm, i0, i3, ierr))
 43:       PetscCallA(DMPlexSetConeSize(dm, i1, i3, ierr))
 44: !     edges
 45:       PetscCallA(DMPlexSetConeSize(dm,  i6, i2, ierr))
 46:       PetscCallA(DMPlexSetConeSize(dm,  i7, i2, ierr))
 47:       PetscCallA(DMPlexSetConeSize(dm,  i8, i2, ierr))
 48:       PetscCallA(DMPlexSetConeSize(dm,  i9, i2, ierr))
 49:       PetscCallA(DMPlexSetConeSize(dm, i10, i2, ierr))

 51:       PetscCallA(DMSetUp(dm, ierr))

 53:       EC(1) = 6
 54:       EC(2) = 7
 55:       EC(3) = 8
 56:       pEC => EC
 57:       PetscCallA(DMPlexSetCone(dm, i0, pEC, ierr))
 58:       EC(1) = 7
 59:       EC(2) = 9
 60:       EC(3) = 10
 61:       pEC => EC
 62:       PetscCallA(DMPlexSetCone(dm, i1 , pEC, ierr))

 64:       VE(1) = 2
 65:       VE(2) = 3
 66:       pVE => VE
 67:       PetscCallA(DMPlexSetCone(dm, i6 , pVE, ierr))
 68:       VE(1) = 3
 69:       VE(2) = 4
 70:       pVE => VE
 71:       PetscCallA(DMPlexSetCone(dm, i7 , pVE, ierr))
 72:       VE(1) = 4
 73:       VE(2) = 2
 74:       pVE => VE
 75:       PetscCallA(DMPlexSetCone(dm, i8 , pVE, ierr))
 76:       VE(1) = 3
 77:       VE(2) = 5
 78:       pVE => VE
 79:       PetscCallA(DMPlexSetCone(dm, i9 , pVE, ierr))
 80:       VE(1) = 5
 81:       VE(2) = 4
 82:       pVE => VE
 83:       PetscCallA(DMPlexSetCone(dm, i10 , pVE, ierr))

 85:       PetscCallA(DMPlexSymmetrize(dm,ierr))
 86:       PetscCallA(DMPlexStratify(dm,ierr))
 87:       PetscCallA(DMView(dm, PETSC_VIEWER_STDOUT_WORLD, ierr))

 89: !     Test Closure
 90:       do cell = 0,1
 91:          PetscCallA(DMPlexGetTransitiveClosure(dm,cell,PETSC_TRUE,nClosure,ierr))
 92: !     Different Fortran compilers print a different number of columns
 93: !     per row producing different outputs in the test runs hence
 94: !     do not print the nClosure
 95:         write(*,1000) 'nClosure ',nClosure
 96:  1000   format (a,30i4)
 97:         PetscCallA(DMPlexRestoreTransitiveClosure(dm,cell,PETSC_TRUE,nClosure,ierr))
 98:       end do

100: !     Test Join
101:       size  = 2
102:       VE(1) = 6
103:       VE(2) = 7
104:       pVE => VE
105:       PetscCallA(DMPlexGetJoin(dm, size, pVE, nJoin, ierr))
106:       write(*,1001) 'Join of',pVE
107:       write(*,1002) '  is',nJoin
108:       PetscCallA(DMPlexRestoreJoin(dm, size, pVE, nJoin, ierr))
109:       size  = 2
110:       VE(1) = 9
111:       VE(2) = 7
112:       pVE => VE
113:       PetscCallA(DMPlexGetJoin(dm, size, pVE, nJoin, ierr))
114:       write(*,1001) 'Join of',pVE
115:  1001 format (a,10i5)
116:        write(*,1002) '  is',nJoin
117:  1002  format (a,10i5)
118:      PetscCallA(DMPlexRestoreJoin(dm, size, pVE, nJoin, ierr))
119: !     Test Full Join
120:       size  = 3
121:       EC(1) = 3
122:       EC(2) = 4
123:       EC(3) = 5
124:       pEC => EC
125:       PetscCallA(DMPlexGetFullJoin(dm, size, pEC, nJoin, ierr))
126:       write(*,1001) 'Full Join of',pEC
127:       write(*,1002) '  is',nJoin
128:       PetscCallA(DMPlexRestoreJoin(dm, size, pEC, nJoin, ierr))
129: !     Test Meet
130:       size  = 2
131:       VE(1) = 0
132:       VE(2) = 1
133:       pVE => VE
134:       PetscCallA(DMPlexGetMeet(dm, size, pVE, nMeet, ierr))
135:       write(*,1001) 'Meet of',pVE
136:       write(*,1002) '  is',nMeet
137:       PetscCallA(DMPlexRestoreMeet(dm, size, pVE, nMeet, ierr))
138:       size  = 2
139:       VE(1) = 6
140:       VE(2) = 7
141:       pVE => VE
142:       PetscCallA(DMPlexGetMeet(dm, size, pVE, nMeet, ierr))
143:       write(*,1001) 'Meet of',pVE
144:       write(*,1002) '  is',nMeet
145:       PetscCallA(DMPlexRestoreMeet(dm, size, pVE, nMeet, ierr))

147:       PetscCallA(DMDestroy(dm, ierr))
148:       PetscCallA(PetscFinalize(ierr))
149:       end
150: !
151: !/*TEST
152: !
153: !   test:
154: !     suffix: 0
155: !
156: !TEST*/