Actual source code: ex2f90.F90

petsc-3.13.4 2020-08-01
Report Typos and Errors
  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:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
 31:       if (ierr .ne. 0) then
 32:         print*,'Unable to initialize PETSc'
 33:         stop
 34:       endif

 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)

 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)

 55:       call DMSetUp(dm, ierr);CHKERRA(ierr)

 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)

 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)

 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)

 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

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)

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*/