Actual source code: ex1f90.F90

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

      7:   DM :: dm
      8:   DMLabel :: label
      9:   Vec :: u
     10:   PetscViewer :: viewer
     11:   PetscSection :: section
     12:   PetscInt :: dim, 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, parameter :: zero = 0, one = 1, two = 2, eight = 8
     21:   PetscMPIInt :: size
     22:   IS, target, dimension(1) ::   bcCompIS
     23:   IS, target, dimension(1) ::   bcPointIS
     24:   IS, pointer :: pBcCompIS(:)
     25:   IS, pointer :: pBcPointIS(:)
     26:   PetscErrorCode :: ierr

     28:   PetscCallA(PetscInitialize(ierr))
     29:   PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD, size, ierr))
     30: !     Create a mesh
     31:   PetscCallA(DMCreate(PETSC_COMM_WORLD, dm, ierr))
     32:   PetscCallA(DMSetType(dm, DMPLEX, ierr))
     33:   PetscCallA(DMSetFromOptions(dm, ierr))
     34:   PetscCallA(DMViewFromOptions(dm, PETSC_NULL_OBJECT, '-dm_view', ierr))
     35:   PetscCallA(DMGetDimension(dm, dim, 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:   PetscCallA(DMGetLabel(dm, 'marker', label, ierr))
     56:   PetscCallA(DMLabelGetValue(label, zero, val, ierr))
     57:   PetscCheckA(size /= 1 .or. val == -1, PETSC_COMM_SELF, PETSC_ERR_PLIB, 'Error in library')
     58:   PetscCallA(DMLabelGetValue(label, eight, val, ierr))
     59:   PetscCheckA(size /= 1 .or. val == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, 'Error in library')
     60: !     Prescribe a Dirichlet condition on u on the boundary
     61: !       Label "marker" is made by the mesh creation routine
     62:   bcField(1) = 0
     63:   pBcField => bcField
     64:   PetscCallA(ISCreateStride(PETSC_COMM_WORLD, one, zero, one, bcCompIS(1), ierr))
     65:   pBcCompIS => bcCompIS
     66:   PetscCallA(DMGetStratumIS(dm, 'marker', one, bcPointIS(1), ierr))
     67:   pBcPointIS => bcPointIS
     68: !     Create a PetscSection with this data layout
     69:   PetscCallA(DMSetNumFields(dm, numFields, ierr))
     70:   PetscCallA(DMPlexCreateSection(dm, PETSC_NULL_DMLABEL_ARRAY, pNumComp, pNumDof, numBC, pBcField, pBcCompIS, pBcPointIS, PETSC_NULL_IS, section, ierr))
     71:   PetscCallA(ISDestroy(bcCompIS(1), ierr))
     72:   PetscCallA(ISDestroy(bcPointIS(1), ierr))
     73: !     Name the Field variables
     74:   PetscCallA(PetscSectionSetFieldName(section, zero, 'u', ierr))
     75:   PetscCallA(PetscSectionSetFieldName(section, one, 'v', ierr))
     76:   PetscCallA(PetscSectionSetFieldName(section, two, 'w', ierr))
     77:   if (size == 1) then
     78:     PetscCallA(PetscSectionView(section, PETSC_VIEWER_STDOUT_WORLD, ierr))
     79:   end if
     80: !     Tell the DM to use this data layout
     81:   PetscCallA(DMSetLocalSection(dm, section, ierr))
     82: !     Create a Vec with this layout and view it
     83:   PetscCallA(DMGetGlobalVector(dm, u, ierr))
     84:   PetscCallA(PetscViewerCreate(PETSC_COMM_WORLD, viewer, ierr))
     85:   PetscCallA(PetscViewerSetType(viewer, PETSCVIEWERVTK, ierr))
     86:   PetscCallA(PetscViewerFileSetName(viewer, 'sol.vtu', ierr))
     87:   PetscCallA(VecView(u, viewer, ierr))
     88:   PetscCallA(PetscViewerDestroy(viewer, ierr))
     89:   PetscCallA(DMRestoreGlobalVector(dm, u, ierr))
     90: !     Cleanup
     91:   PetscCallA(PetscSectionDestroy(section, ierr))
     92:   PetscCallA(DMDestroy(dm, ierr))

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

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