Actual source code: ex1f90.F90

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

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

108: !/*TEST
109: !  build:
110: !    requires: defined(PETSC_USING_F90FREEFORM)
111: !
112: !  test:
113: !    suffix: 0
114: !    requires: triangle
115: !    args: -info :~sys,mat:
116: !
117: !  test:
118: !    suffix: 0_2
119: !    requires: triangle
120: !    nsize: 2
121: !    args: -info :~sys,mat,partitioner:
122: !
123: !  test:
124: !    suffix: 1
125: !    requires: ctetgen
126: !    args: -dm_plex_dim 3 -info :~sys,mat:
127: !
128: !TEST*/