1:       program DMPlexTestField
  2: #include "petsc/finclude/petscdmplex.h"
  3: #include "petsc/finclude/petscdmlabel.h"
  4:       use petscdmplex
  5:       implicit none

  7:       DM :: dm
  8:       DMLabel :: label
  9:       Vec :: u
 10:       PetscViewer :: viewer
 11:       PetscSection :: section
 12:       PetscInt :: dim,numCells,numFields,numBC
 13:       PetscInt :: i,val
 14:       PetscInt, target, dimension(3) ::  numComp
 15:       PetscInt, pointer :: pNumComp(:)
 16:       PetscInt, target, dimension(12) ::  numDof
 17:       PetscInt, pointer :: pNumDof(:)
 18:       PetscInt, target, dimension(1) ::  bcField
 19:       PetscInt, pointer :: pBcField(:)
 20:       PetscInt :: zero,eight
 21:       IS, target, dimension(1) ::   bcCompIS
 22:       IS, target, dimension(1) ::   bcPointIS
 23:       IS, pointer :: pBcCompIS(:)
 24:       IS, pointer :: pBcPointIS(:)
 25:       PetscBool :: interpolate
 26:       PetscErrorCode :: ierr

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

106:       call PetscFinalize(ierr)
107:       end program DMPlexTestField

109: !/*TEST
110: !  build:
111: !    requires: define(PETSC_USING_F90FREEFORM)
112: !
113: !  test:
114: !    suffix: 0
115: !    requires: triangle
116: !
117: !  test:
118: !    suffix: 1
119: !    requires: ctetgen
120: !    args: -dim 3
121: !
122: !TEST*/