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

  8:       DM :: dm
  9:       DMLabel :: label
 10:       Vec :: u
 11:       PetscViewer :: viewer
 12:       PetscSection :: section
 13:       PetscInt :: dim,numFields,numBC
 14:       PetscInt :: i,val
 15:       DMLabel, pointer :: nolabel(:) => NULL()
 16:       PetscInt, target, dimension(3) ::  numComp
 17:       PetscInt, pointer :: pNumComp(:)
 18:       PetscInt, target, dimension(12) ::  numDof
 19:       PetscInt, pointer :: pNumDof(:)
 20:       PetscInt, target, dimension(1) ::  bcField
 21:       PetscInt, pointer :: pBcField(:)
 22:       PetscInt, parameter :: zero = 0, one = 1, two = 2, eight = 8
 23:       PetscMPIInt :: size
 24:       IS, target, dimension(1) ::   bcCompIS
 25:       IS, target, dimension(1) ::   bcPointIS
 26:       IS, pointer :: pBcCompIS(:)
 27:       IS, pointer :: pBcPointIS(:)
 28:       PetscErrorCode :: ierr

 30:       PetscCallA(PetscInitialize(ierr))
 31:       PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD, size, ierr))
 32: !     Create a mesh
 33:       PetscCallA(DMCreate(PETSC_COMM_WORLD, dm, ierr))
 34:       PetscCallA(DMSetType(dm, DMPLEX, ierr))
 35:       PetscCallA(DMSetFromOptions(dm, ierr))
 36:       PetscCallA(DMViewFromOptions(dm, PETSC_NULL_VEC, '-dm_view', ierr))
 37:       PetscCallA(DMGetDimension(dm, dim, ierr))
 38: !     Create a scalar field u, a vector field v, and a surface vector field w
 39:       numFields  = 3
 40:       numComp(1) = 1
 41:       numComp(2) = dim
 42:       numComp(3) = dim-1
 43:       pNumComp => numComp
 44:       do i = 1, numFields*(dim+1)
 45:          numDof(i) = 0
 46:       end do
 47: !     Let u be defined on vertices
 48:       numDof(0*(dim+1)+1)     = 1
 49: !     Let v be defined on cells
 50:       numDof(1*(dim+1)+dim+1) = dim
 51: !     Let v be defined on faces
 52:       numDof(2*(dim+1)+dim)   = dim-1
 53:       pNumDof => numDof
 54: !     Setup boundary conditions
 55:       numBC = 1
 56: !     Test label retrieval
 57:       PetscCallA(DMGetLabel(dm, 'marker', label, ierr))
 58:       PetscCallA(DMLabelGetValue(label, zero, val, ierr))
 59:       PetscCheckA(size .ne. 1 .or. val .eq. -1,PETSC_COMM_SELF,PETSC_ERR_PLIB,'Error in library')
 60:       PetscCallA(DMLabelGetValue(label, eight, val, ierr))
 61:       PetscCheckA(size .ne. 1 .or. val .eq. 1,PETSC_COMM_SELF,PETSC_ERR_PLIB,'Error in library')
 62: !     Prescribe a Dirichlet condition on u on the boundary
 63: !       Label "marker" is made by the mesh creation routine
 64:       bcField(1) = 0
 65:       pBcField => bcField
 66:       PetscCallA(ISCreateStride(PETSC_COMM_WORLD, one, zero, one, bcCompIS(1), ierr))
 67:       pBcCompIS => bcCompIS
 68:       PetscCallA(DMGetStratumIS(dm, 'marker', one, bcPointIS(1),ierr))
 69:       pBcPointIS => bcPointIS
 70: !     Create a PetscSection with this data layout
 71:       PetscCallA(DMSetNumFields(dm, numFields,ierr))
 72:       PetscCallA(DMPlexCreateSection(dm,nolabel,pNumComp,pNumDof,numBC,pBcField,pBcCompIS,pBcPointIS,PETSC_NULL_IS,section,ierr))
 73:       PetscCallA(ISDestroy(bcCompIS(1), ierr))
 74:       PetscCallA(ISDestroy(bcPointIS(1), ierr))
 75: !     Name the Field variables
 76:       PetscCallA(PetscSectionSetFieldName(section, zero, 'u', ierr))
 77:       PetscCallA(PetscSectionSetFieldName(section, one,  'v', ierr))
 78:       PetscCallA(PetscSectionSetFieldName(section, two,  'w', ierr))
 79:       if (size .eq. 1) then
 80:         PetscCallA(PetscSectionView(section, PETSC_VIEWER_STDOUT_WORLD, ierr))
 81:       endif
 82: !     Tell the DM to use this data layout
 83:       PetscCallA(DMSetLocalSection(dm, section, ierr))
 84: !     Create a Vec with this layout and view it
 85:       PetscCallA(DMGetGlobalVector(dm, u, ierr))
 86:       PetscCallA(PetscViewerCreate(PETSC_COMM_WORLD, viewer, ierr))
 87:       PetscCallA(PetscViewerSetType(viewer, PETSCVIEWERVTK, ierr))
 88:       PetscCallA(PetscViewerFileSetName(viewer, 'sol.vtu', ierr))
 89:       PetscCallA(VecView(u, viewer, ierr))
 90:       PetscCallA(PetscViewerDestroy(viewer, ierr))
 91:       PetscCallA(DMRestoreGlobalVector(dm, u, ierr))
 92: !     Cleanup
 93:       PetscCallA(PetscSectionDestroy(section, ierr))
 94:       PetscCallA(DMDestroy(dm, ierr))

 96:       PetscCallA(PetscFinalize(ierr))
 97:       end program DMPlexTestField

 99: !/*TEST
100: !  build:
101: !    requires: defined(PETSC_USING_F90FREEFORM)
102: !
103: !  test:
104: !    suffix: 0
105: !    requires: triangle
106: !    args: -info :~sys,mat:
107: !
108: !  test:
109: !    suffix: 0_2
110: !    requires: triangle
111: !    nsize: 2
112: !    args: -info :~sys,mat,partitioner:
113: !
114: !  test:
115: !    suffix: 1
116: !    requires: ctetgen
117: !    args: -dm_plex_dim 3 -info :~sys,mat:
118: !
119: !TEST*/