1:       program DMPlexTestField
  2:       implicit none
  3: #include "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:      &     bcPointIS
 21:       IS, pointer :: pBcPointIS(:)
 22:       PetscBool :: interpolate
 23:       PetscErrorCode :: ierr

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

117:       call PetscFinalize(ierr)
118:       end program DMPlexTestField