1c4762a1bSJed Brown program DMPlexTestField 2c4762a1bSJed Brown#include "petsc/finclude/petscdmplex.h" 3c4762a1bSJed Brown#include "petsc/finclude/petscdmlabel.h" 4c4762a1bSJed Brown use petscdmplex 5c4762a1bSJed Brown use petscsys 6c4762a1bSJed Brown implicit none 7c4762a1bSJed Brown 8c4762a1bSJed Brown DM :: dm 9c4762a1bSJed Brown DMLabel :: label 10c4762a1bSJed Brown Vec :: u 11c4762a1bSJed Brown PetscViewer :: viewer 12c4762a1bSJed Brown PetscSection :: section 13c4762a1bSJed Brown PetscInt :: dim,numFields,numBC 14c4762a1bSJed Brown PetscInt :: i,val 15c4762a1bSJed Brown DMLabel, pointer :: nolabel(:) => NULL() 16c4762a1bSJed Brown PetscInt, target, dimension(3) :: numComp 17c4762a1bSJed Brown PetscInt, pointer :: pNumComp(:) 18c4762a1bSJed Brown PetscInt, target, dimension(12) :: numDof 19c4762a1bSJed Brown PetscInt, pointer :: pNumDof(:) 20c4762a1bSJed Brown PetscInt, target, dimension(1) :: bcField 21c4762a1bSJed Brown PetscInt, pointer :: pBcField(:) 22c4762a1bSJed Brown PetscInt, parameter :: zero = 0, one = 1, two = 2, eight = 8 23c4762a1bSJed Brown IS, target, dimension(1) :: bcCompIS 24c4762a1bSJed Brown IS, target, dimension(1) :: bcPointIS 25c4762a1bSJed Brown IS, pointer :: pBcCompIS(:) 26c4762a1bSJed Brown IS, pointer :: pBcPointIS(:) 27c4762a1bSJed Brown PetscBool :: interpolate,flg 28c4762a1bSJed Brown PetscErrorCode :: ierr 29c4762a1bSJed Brown 30c4762a1bSJed Brown call PetscInitialize(PETSC_NULL_CHARACTER, ierr) 31c4762a1bSJed Brown if (ierr .ne. 0) then 32c4762a1bSJed Brown print*,'Unable to initialize PETSc' 33c4762a1bSJed Brown stop 34c4762a1bSJed Brown endif 35c4762a1bSJed Brown dim = 2 36c4762a1bSJed Brown call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-dim', dim,flg, ierr);CHKERRA(ierr) 37c4762a1bSJed Brown interpolate = PETSC_TRUE 38c4762a1bSJed Brown! Create a mesh 39c4762a1bSJed Brown call DMPlexCreateBoxMesh(PETSC_COMM_WORLD, dim, PETSC_TRUE, PETSC_NULL_INTEGER, PETSC_NULL_REAL, PETSC_NULL_REAL, PETSC_NULL_INTEGER, interpolate, dm, ierr);CHKERRA(ierr) 40c4762a1bSJed Brown! Create a scalar field u, a vector field v, and a surface vector field w 41c4762a1bSJed Brown numFields = 3 42c4762a1bSJed Brown numComp(1) = 1 43c4762a1bSJed Brown numComp(2) = dim 44c4762a1bSJed Brown numComp(3) = dim-1 45c4762a1bSJed Brown pNumComp => numComp 46c4762a1bSJed Brown do i = 1, numFields*(dim+1) 47c4762a1bSJed Brown numDof(i) = 0 48c4762a1bSJed Brown end do 49c4762a1bSJed Brown! Let u be defined on vertices 50c4762a1bSJed Brown numDof(0*(dim+1)+1) = 1 51c4762a1bSJed Brown! Let v be defined on cells 52c4762a1bSJed Brown numDof(1*(dim+1)+dim+1) = dim 53c4762a1bSJed Brown! Let v be defined on faces 54c4762a1bSJed Brown numDof(2*(dim+1)+dim) = dim-1 55c4762a1bSJed Brown pNumDof => numDof 56c4762a1bSJed Brown! Setup boundary conditions 57c4762a1bSJed Brown numBC = 1 58c4762a1bSJed Brown! Test label retrieval 59c4762a1bSJed Brown call DMGetLabel(dm, 'marker', label, ierr);CHKERRA(ierr) 60c4762a1bSJed Brown call DMLabelGetValue(label, zero, val, ierr);CHKERRA(ierr) 61c4762a1bSJed Brown if (val .ne. -1) then 62c4762a1bSJed Brown SETERRA(PETSC_COMM_SELF,PETSC_ERR_PLIB,'Error in library') 63c4762a1bSJed Brown endif 64c4762a1bSJed Brown call DMLabelGetValue(label, eight, val, ierr);CHKERRA(ierr) 65c4762a1bSJed Brown if (val .ne. 1) then 66c4762a1bSJed Brown SETERRA(PETSC_COMM_SELF,PETSC_ERR_PLIB,'Error in library') 67c4762a1bSJed Brown endif 68c4762a1bSJed Brown! Prescribe a Dirichlet condition on u on the boundary 69c4762a1bSJed Brown! Label "marker" is made by the mesh creation routine 70c4762a1bSJed Brown bcField(1) = 0 71c4762a1bSJed Brown pBcField => bcField 72c4762a1bSJed Brown call ISCreateStride(PETSC_COMM_WORLD, one, zero, one, bcCompIS(1), ierr);CHKERRA(ierr) 73c4762a1bSJed Brown pBcCompIS => bcCompIS 74c4762a1bSJed Brown call DMGetStratumIS(dm, 'marker', one, bcPointIS(1),ierr);CHKERRA(ierr) 75c4762a1bSJed Brown pBcPointIS => bcPointIS 76c4762a1bSJed Brown! Create a PetscSection with this data layout 77c4762a1bSJed Brown call DMSetNumFields(dm, numFields,ierr);CHKERRA(ierr) 78c4762a1bSJed Brown call DMPlexCreateSection(dm,nolabel,pNumComp,pNumDof,numBC,pBcField,pBcCompIS,pBcPointIS,PETSC_NULL_IS,section,ierr) 79c4762a1bSJed Brown CHKERRA(ierr) 80c4762a1bSJed Brown call ISDestroy(bcCompIS(1), ierr);CHKERRA(ierr) 81c4762a1bSJed Brown call ISDestroy(bcPointIS(1), ierr);CHKERRA(ierr) 82c4762a1bSJed Brown! Name the Field variables 83c4762a1bSJed Brown call PetscSectionSetFieldName(section, zero, 'u', ierr);CHKERRA(ierr) 84c4762a1bSJed Brown call PetscSectionSetFieldName(section, one, 'v', ierr);CHKERRA(ierr) 85c4762a1bSJed Brown call PetscSectionSetFieldName(section, two, 'w', ierr);CHKERRA(ierr) 86c4762a1bSJed Brown call PetscSectionView(section, PETSC_VIEWER_STDOUT_WORLD, ierr);CHKERRA(ierr) 87c4762a1bSJed Brown! Tell the DM to use this data layout 88c4762a1bSJed Brown call DMSetLocalSection(dm, section, ierr);CHKERRA(ierr) 89c4762a1bSJed Brown! Create a Vec with this layout and view it 90c4762a1bSJed Brown call DMGetGlobalVector(dm, u, ierr);CHKERRA(ierr) 91c4762a1bSJed Brown call PetscViewerCreate(PETSC_COMM_WORLD, viewer, ierr);CHKERRA(ierr) 92c4762a1bSJed Brown call PetscViewerSetType(viewer, PETSCVIEWERVTK, ierr);CHKERRA(ierr) 93*8ec8862eSJed Brown call PetscViewerFileSetName(viewer, 'sol.vtu', ierr);CHKERRA(ierr) 94c4762a1bSJed Brown call VecView(u, viewer, ierr);CHKERRA(ierr) 95c4762a1bSJed Brown call PetscViewerDestroy(viewer, ierr);CHKERRA(ierr) 96c4762a1bSJed Brown call DMRestoreGlobalVector(dm, u, ierr);CHKERRA(ierr) 97c4762a1bSJed Brown! Cleanup 98c4762a1bSJed Brown call PetscSectionDestroy(section, ierr);CHKERRA(ierr) 99c4762a1bSJed Brown call DMDestroy(dm, ierr);CHKERRA(ierr) 100c4762a1bSJed Brown 101c4762a1bSJed Brown call PetscFinalize(ierr) 102c4762a1bSJed Brown end program DMPlexTestField 103c4762a1bSJed Brown 104c4762a1bSJed Brown!/*TEST 105c4762a1bSJed Brown! build: 106c4762a1bSJed Brown! requires: define(PETSC_USING_F90FREEFORM) 107c4762a1bSJed Brown! 108c4762a1bSJed Brown! test: 109c4762a1bSJed Brown! suffix: 0 110c4762a1bSJed Brown! requires: triangle 111c20d7725SJed Brown! args: -info :~sys: 112c4762a1bSJed Brown! 113c4762a1bSJed Brown! test: 114c4762a1bSJed Brown! suffix: 1 115c4762a1bSJed Brown! requires: ctetgen 116c20d7725SJed Brown! args: -dim 3 -info :~sys: 117c4762a1bSJed Brown! 118c4762a1bSJed Brown!TEST*/ 119