Actual source code: ex1f90.F

petsc-master 2016-07-26
Report Typos and Errors
  1:       program DMPlexTestField
  2:       implicit none
  3: #include "petsc/finclude/petsc.h90"
  4:       DM :: dm
  5:       DMLabel :: label
  6:       Vec :: u
  7:       PetscViewer :: viewer
  8:       PetscSection :: section
  9:       PetscInt :: dim,numCells,numFields,numBC
 10:       PetscInt :: i,val
 11:       PetscInt, target, dimension(3) ::                                 &
 12:      &     numComp
 13:       PetscInt, pointer :: pNumComp(:)
 14:       PetscInt, target, dimension(12) ::                                &
 15:      &     numDof
 16:       PetscInt, pointer :: pNumDof(:)
 17:       PetscInt, target, dimension(1) ::                                 &
 18:      &     bcField
 19:       PetscInt, pointer :: pBcField(:)
 20:       IS, target, dimension(1) ::                                       &
 21:      &     bcCompIS
 22:       IS, target, dimension(1) ::                                       &
 23:      &     bcPointIS
 24:       IS, pointer :: pBcCompIS(:)
 25:       IS, pointer :: pBcPointIS(:)
 26:       PetscBool :: interpolate
 27:       PetscErrorCode :: ierr

 29:       call PetscInitialize(PETSC_NULL_CHARACTER, ierr)
 30:       CHKERRQ(ierr)
 31:       dim = 2
 32:       call PetscOptionsGetInt(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER,   &
 33:      &                        '-dim', dim,PETSC_NULL_BOOL, ierr)
 34:       CHKERRQ(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,         &
 43:      &     interpolate, dm, ierr)
 44:       CHKERRQ(ierr)
 45: !     Create a scalar field u, a vector field v, and a surface vector field w
 46:       numFields  = 3
 47:       numComp(1) = 1
 48:       numComp(2) = dim
 49:       numComp(3) = dim-1
 50:       pNumComp => numComp
 51:       do i = 1, numFields*(dim+1)
 52:          numDof(i) = 0
 53:       end do
 54: !     Let u be defined on vertices
 55:       numDof(0*(dim+1)+1)     = 1
 56: !     Let v be defined on cells
 57:       numDof(1*(dim+1)+dim+1) = dim
 58: !     Let v be defined on faces
 59:       numDof(2*(dim+1)+dim)   = dim-1
 60:       pNumDof => numDof
 61: !     Setup boundary conditions
 62:       numBC = 1
 63: !     Test label retrieval
 64:       call DMGetLabel(dm, 'marker', label, ierr)
 65:       CHKERRQ(ierr)
 66:       call DMLabelGetValue(label, 0, val, ierr)
 67:       CHKERRQ(ierr)
 68:       if (val .ne. -1) then
 69:         CHKERRQ(1)
 70:       endif
 71:       call DMLabelGetValue(label, 8, val, ierr)
 72:       CHKERRQ(ierr)
 73:       if (val .ne. 1) then
 74:         CHKERRQ(1)
 75:       endif
 76: !     Prescribe a Dirichlet condition on u on the boundary
 77: !       Label "marker" is made by the mesh creation routine
 78:       bcField(1) = 0
 79:       pBcField => bcField
 80:       call ISCreateStride(PETSC_COMM_WORLD, 1, 0, 1, bcCompIS(1), ierr)
 81:       CHKERRQ(ierr)
 82:       pBcCompIS => bcCompIS
 83:       call DMGetStratumIS(dm, 'marker', 1, bcPointIS(1),
 84:      &    ierr)
 85:       CHKERRQ(ierr)
 86:       pBcPointIS => bcPointIS
 87: !     Create a PetscSection with this data layout
 88:       call DMPlexCreateSection(dm, dim, numFields, pNumComp,
 89:      &     pNumDof, numBC, pBcField, pBcCompIS, pBcPointIS,
 90:      &     PETSC_NULL_OBJECT, section, ierr)
 91:       CHKERRQ(ierr)
 92:       call ISDestroy(bcCompIS(1), ierr)
 93:       CHKERRQ(ierr)
 94:       call ISDestroy(bcPointIS(1), ierr)
 95:       CHKERRQ(ierr)
 96: !     Name the Field variables
 97:       call PetscSectionSetFieldName(section, 0, 'u', ierr)
 98:       CHKERRQ(ierr)
 99:       call PetscSectionSetFieldName(section, 1, 'v', ierr)
100:       CHKERRQ(ierr)
101:       call PetscSectionSetFieldName(section, 2, 'w', ierr)
102:       CHKERRQ(ierr)
103:       call PetscSectionView(section, PETSC_VIEWER_STDOUT_WORLD, ierr)
104:       CHKERRQ(ierr)
105: !     Tell the DM to use this data layout
106:       call DMSetDefaultSection(dm, section, ierr)
107:       CHKERRQ(ierr)
108: !     Create a Vec with this layout and view it
109:       call DMGetGlobalVector(dm, u, ierr)
110:       CHKERRQ(ierr)
111:       call PetscViewerCreate(PETSC_COMM_WORLD, viewer, ierr)
112:       CHKERRQ(ierr)
113:       call PetscViewerSetType(viewer, PETSCVIEWERVTK, ierr)
114:       CHKERRQ(ierr)
115:       call PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_VTK, ierr)
116:       CHKERRQ(ierr)
117:       call PetscViewerFileSetName(viewer, 'sol.vtk', ierr)
118:       CHKERRQ(ierr)
119:       call VecView(u, viewer, ierr)
120:       CHKERRQ(ierr)
121:       call PetscViewerDestroy(viewer, ierr)
122:       CHKERRQ(ierr)
123:       call DMRestoreGlobalVector(dm, u, ierr)
124:       CHKERRQ(ierr)
125: !     Cleanup
126:       call PetscSectionDestroy(section, ierr)
127:       CHKERRQ(ierr)
128:       call DMDestroy(dm, ierr)
129:       CHKERRQ(ierr)

131:       call PetscFinalize(ierr)
132:       end program DMPlexTestField