1: #include <petsc/finclude/petscdmplex.h>
  2: #include <petsc/finclude/petscdmlabel.h>
  3: program DMPlexTestField
  4:   use petscdm
  5:   use petscdmplex
  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:   PetscInt, target, dimension(3) ::  numComp
 16:   PetscInt, pointer :: pNumComp(:)
 17:   PetscInt, target, dimension(12) ::  numDof
 18:   PetscInt, pointer :: pNumDof(:)
 19:   PetscInt, target, dimension(1) ::  bcField
 20:   PetscInt, pointer :: pBcField(:)
 21:   PetscInt, parameter :: zero = 0, one = 1, two = 2, eight = 8
 22:   PetscMPIInt :: size
 23:   IS, target, dimension(1) ::   bcCompIS
 24:   IS, target, dimension(1) ::   bcPointIS
 25:   IS, pointer :: pBcCompIS(:)
 26:   IS, pointer :: pBcPointIS(:)
 27:   PetscErrorCode :: ierr

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

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

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