1:       program DMPlexTestField
  2:       implicit none
  3: #include "petsc/finclude/petsc.h90"
  4:       DM :: dm
  5:       DMLabel :: label
  6:       Vec :: u
  7:       PetscViewer :: viewer
  8:       PetscSection :: section
  9:       PetscInt :: dim,numFields,numBC,i,val
 10:       PetscInt, target, dimension(3) ::                                 &
 11:      &     numComp
 12:       PetscInt, pointer :: pNumComp(:)
 13:       PetscInt, target, dimension(12) ::                                &
 14:      &     numDof
 15:       PetscInt, pointer :: pNumDof(:)
 16:       PetscInt, target, dimension(1) ::                                 &
 17:      &     bcField
 18:       PetscInt, pointer :: pBcField(:)
 19:       IS, target, dimension(1) ::                                       &
 20:      &     bcCompIS
 21:       IS, target, dimension(1) ::                                       &
 22:      &     bcPointIS
 23:       IS, pointer :: pBcCompIS(:)
 24:       IS, pointer :: pBcPointIS(:)
 25:       PetscBool :: interpolate
 26:       PetscErrorCode :: ierr

 28:       call PetscInitialize(PETSC_NULL_CHARACTER, ierr)
 29:       CHKERRQ(ierr)
 30:       dim = 2
 31:       call PetscOptionsGetInt(PETSC_NULL_CHARACTER, "-dim", dim,
 32:      $     PETSC_NULL_BOOL, ierr)
 33:       CHKERRQ(ierr)
 34:       interpolate = PETSC_TRUE
 35: !     Create a mesh
 36:       call DMPlexCreateBoxMesh(PETSC_COMM_WORLD, dim, interpolate, dm,
 37:      $     ierr)
 38:       CHKERRQ(ierr)
 39: !     Create a scalar field u, a vector field v, and a surface vector field w
 40:       numFields  = 3
 41:       numComp(1) = 1
 42:       numComp(2) = dim
 43:       numComp(3) = dim-1
 44:       pNumComp => numComp
 45:       do i = 1, numFields*(dim+1)
 46:          numDof(i) = 0
 47:       end do
 48: !     Let u be defined on vertices
 49:       numDof(0*(dim+1)+1)     = 1
 50: !     Let v be defined on cells
 51:       numDof(1*(dim+1)+dim+1) = dim
 52: !     Let v be defined on faces
 53:       numDof(2*(dim+1)+dim)   = dim-1
 54:       pNumDof => numDof
 55: !     Setup boundary conditions
 56:       numBC = 1
 57: !     Test label retrieval
 58:       call DMPlexGetLabel(dm, "marker", label, ierr)
 59:       CHKERRQ(ierr)
 60:       call DMLabelGetValue(label, 0, val, ierr)
 61:       CHKERRQ(ierr)
 62:       if (val .ne. -1) then
 63:         CHKERRQ(1)
 64:       endif
 65:       call DMLabelGetValue(label, 8, val, ierr)
 66:       CHKERRQ(ierr)
 67:       if (val .ne. 1) then
 68:         CHKERRQ(1)
 69:       endif
 70: !     Prescribe a Dirichlet condition on u on the boundary
 71: !       Label "marker" is made by the mesh creation routine
 72:       bcField(1) = 0
 73:       pBcField => bcField
 74:       call ISCreateStride(PETSC_COMM_WORLD, 1, 0, 1, bcCompIS(1), ierr)
 75:       CHKERRQ(ierr)
 76:       pBcCompIS => bcCompIS
 77:       call DMPlexGetStratumIS(dm, "marker", 1, bcPointIS(1),
 78:      $     ierr)
 79:       CHKERRQ(ierr)
 80:       pBcPointIS => bcPointIS
 81: !     Create a PetscSection with this data layout
 82:       call DMPlexCreateSection(dm, dim, numFields, pNumComp,
 83:      $     pNumDof, numBC, pBcField, pBcCompIS, pBcPointIS,
 84:      $     PETSC_NULL_OBJECT, section, ierr)
 85:       CHKERRQ(ierr)
 86:       call ISDestroy(bcCompIS(1), ierr)
 87:       CHKERRQ(ierr)
 88:       call ISDestroy(bcPointIS(1), ierr)
 89:       CHKERRQ(ierr)
 90: !     Name the Field variables
 91:       call PetscSectionSetFieldName(section, 0, "u", ierr)
 92:       CHKERRQ(ierr)
 93:       call PetscSectionSetFieldName(section, 1, "v", ierr)
 94:       CHKERRQ(ierr)
 95:       call PetscSectionSetFieldName(section, 2, "w", ierr)
 96:       CHKERRQ(ierr)
 97:       call PetscSectionView(section, PETSC_VIEWER_STDOUT_WORLD, ierr)
 98:       CHKERRQ(ierr)
 99: !     Tell the DM to use this data layout
100:       call DMSetDefaultSection(dm, section, ierr)
101:       CHKERRQ(ierr)
102: !     Create a Vec with this layout and view it
103:       call DMGetGlobalVector(dm, u, ierr)
104:       CHKERRQ(ierr)
105:       call PetscViewerCreate(PETSC_COMM_WORLD, viewer, ierr)
106:       CHKERRQ(ierr)
107:       call PetscViewerSetType(viewer, PETSCVIEWERVTK, ierr)
108:       CHKERRQ(ierr)
109:       call PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_VTK, ierr)
110:       CHKERRQ(ierr)
111:       call PetscViewerFileSetName(viewer, 'sol.vtk', ierr)
112:       CHKERRQ(ierr)
113:       call VecView(u, viewer, ierr)
114:       CHKERRQ(ierr)
115:       call PetscViewerDestroy(viewer, ierr)
116:       CHKERRQ(ierr)
117:       call DMRestoreGlobalVector(dm, u, ierr)
118:       CHKERRQ(ierr)
119: !     Cleanup
120:       call PetscSectionDestroy(section, ierr)
121:       CHKERRQ(ierr)
122:       call DMDestroy(dm, ierr)
123:       CHKERRQ(ierr)

125:       call PetscFinalize(ierr)
126:       end program DMPlexTestField