Actual source code: ex2f90.F

petsc-3.4.2 2013-07-02
  1:       program main
  2:       implicit none
  3: !
  4: #include <finclude/petsc.h90>
  5: !
  6:       DM dm
  7:       PetscInt, target, dimension(3) :: EC
  8:       PetscInt, target, dimension(2) :: VE
  9:       PetscInt, pointer :: pEC(:), pVE(:)
 10:       PetscInt, pointer :: nClosure(:)
 11:       PetscInt, pointer :: nJoin(:)
 12:       PetscInt, pointer :: nMeet(:)
 13:       PetscInt       dim, cell, size
 14:       PetscInt i0,i1,i2,i3,i4,i5,i6,i7
 15:       PetscInt i8,i9,i10,i11
 16:       PetscErrorCode ierr

 18:       i0 = 0
 19:       i1 = 1
 20:       i2 = 2
 21:       i3 = 3
 22:       i4 = 4
 23:       i5 = 5
 24:       i6 = 6
 25:       i7 = 7
 26:       i8 = 8
 27:       i9 = 9
 28:       i10 = 10
 29:       i11 = 11

 31:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
 32:       CHKERRQ(ierr)
 33:       call DMPlexCreate(PETSC_COMM_WORLD, dm, ierr)
 34:       CHKERRQ(ierr)
 35:       dim = 2
 36:       call DMPlexSetDimension(dm, dim, ierr)
 37:       CHKERRQ(ierr)

 39: ! Make Doublet Mesh from Fig 2 of Flexible Representation of Computational Meshes,
 40: ! except indexing is from 0 instead of 1 and we obey the new restrictions on
 41: ! numbering: cells, vertices, faces, edges
 42:       call DMPlexSetChart(dm, i0, i11, ierr)
 43:       CHKERRQ(ierr)
 44: !     cells
 45:       call DMPlexSetConeSize(dm, i0, i3, ierr)
 46:       CHKERRQ(ierr)
 47:       call DMPlexSetConeSize(dm, i1, i3, ierr)
 48:       CHKERRQ(ierr)
 49: !     edges
 50:       call DMPlexSetConeSize(dm,  i6, i2, ierr)
 51:       CHKERRQ(ierr)
 52:       call DMPlexSetConeSize(dm,  i7, i2, ierr)
 53:       CHKERRQ(ierr)
 54:       call DMPlexSetConeSize(dm,  i8, i2, ierr)
 55:       CHKERRQ(ierr)
 56:       call DMPlexSetConeSize(dm,  i9, i2, ierr)
 57:       CHKERRQ(ierr)
 58:       call DMPlexSetConeSize(dm, i10, i2, ierr)
 59:       CHKERRQ(ierr)

 61:       call DMSetUp(dm, ierr)
 62:       CHKERRQ(ierr)

 64:       EC(1) = 6
 65:       EC(2) = 7
 66:       EC(3) = 8
 67:       pEC => EC
 68:       call DMPlexSetCone(dm, i0, pEC, ierr)
 69:       CHKERRQ(ierr)
 70:       EC(1) = 7
 71:       EC(2) = 9
 72:       EC(3) = 10
 73:       pEC => EC
 74:       call DMPlexSetCone(dm, i1 , pEC, ierr)
 75:       CHKERRQ(ierr)

 77:       VE(1) = 2
 78:       VE(2) = 3
 79:       pVE => VE
 80:       call DMPlexSetCone(dm, i6 , pVE, ierr)
 81:       CHKERRQ(ierr)
 82:       VE(1) = 3
 83:       VE(2) = 4
 84:       pVE => VE
 85:       call DMPlexSetCone(dm, i7 , pVE, ierr)
 86:       CHKERRQ(ierr)
 87:       VE(1) = 4
 88:       VE(2) = 2
 89:       pVE => VE
 90:       call DMPlexSetCone(dm, i8 , pVE, ierr)
 91:       CHKERRQ(ierr)
 92:       VE(1) = 3
 93:       VE(2) = 5
 94:       pVE => VE
 95:       call DMPlexSetCone(dm, i9 , pVE, ierr)
 96:       CHKERRQ(ierr)
 97:       VE(1) = 5
 98:       VE(2) = 4
 99:       pVE => VE
100:       call DMPlexSetCone(dm, i10 , pVE, ierr)
101:       CHKERRQ(ierr)

103:       call DMPlexSymmetrize(dm,ierr)
104:       CHKERRQ(ierr)
105:       call DMPlexStratify(dm,ierr)
106:       CHKERRQ(ierr)
107:       call DMView(dm, PETSC_VIEWER_STDOUT_WORLD, ierr)

109: !     Test Closure
110:       do cell = 0,1
111:          call DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE,          &
112:      &        nClosure, ierr)
113:          CHKERRQ(ierr)
114:          write(*,*) nClosure
115:          call DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE,      &
116:      &        nClosure, ierr)
117:       end do

119: !     Test Join
120:       size  = 2
121:       VE(1) = 6
122:       VE(2) = 7
123:       pVE => VE
124:       call DMPlexGetJoin(dm, size, pVE, nJoin, ierr)
125:       write(*,*) 'Join of',pVE,'is',nJoin
126:       call DMPlexRestoreJoin(dm, size, pVE, nJoin, ierr)
127:       size  = 2
128:       VE(1) = 9
129:       VE(2) = 7
130:       pVE => VE
131:       call DMPlexGetJoin(dm, size, pVE, nJoin, ierr)
132:       write(*,*) 'Join of',pVE,'is',nJoin
133:       call DMPlexRestoreJoin(dm, size, pVE, nJoin, ierr)
134: !     Test Full Join
135:       size  = 3
136:       EC(1) = 3
137:       EC(2) = 4
138:       EC(3) = 5
139:       pEC => EC
140:       call DMPlexGetFullJoin(dm, size, pEC, nJoin, ierr)
141:       write(*,*) 'Full Join of',pEC,'is',nJoin
142:       call DMPlexRestoreJoin(dm, size, pEC, nJoin, ierr)
143: !     Test Meet
144:       size  = 2
145:       VE(1) = 0
146:       VE(2) = 1
147:       pVE => VE
148:       call DMPlexGetMeet(dm, size, pVE, nMeet, ierr)
149:       write(*,*) 'Meet of',pVE,'is',nMeet
150:       call DMPlexRestoreMeet(dm, size, pVE, nMeet, ierr)
151:       size  = 2
152:       VE(1) = 6
153:       VE(2) = 7
154:       pVE => VE
155:       call DMPlexGetMeet(dm, size, pVE, nMeet, ierr)
156:       write(*,*) 'Meet of',pVE,'is',nMeet
157:       call DMPlexRestoreMeet(dm, size, pVE, nMeet, ierr)

159:       call DMDestroy(dm, ierr)
160:       CHKERRQ(ierr)
161:       call PetscFinalize(ierr)
162:       CHKERRQ(ierr)
163:       end