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 23*2e3d3ef9SMatthew G. Knepley PetscMPIInt :: size 24c4762a1bSJed Brown IS, target, dimension(1) :: bcCompIS 25c4762a1bSJed Brown IS, target, dimension(1) :: bcPointIS 26c4762a1bSJed Brown IS, pointer :: pBcCompIS(:) 27c4762a1bSJed Brown IS, pointer :: pBcPointIS(:) 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 35*2e3d3ef9SMatthew G. Knepley call MPI_Comm_size(PETSC_COMM_WORLD, size, ierr);CHKERRA(ierr) 36c4762a1bSJed Brown! Create a mesh 3730602db0SMatthew G. Knepley call DMCreate(PETSC_COMM_WORLD, dm, ierr);CHKERRA(ierr) 3830602db0SMatthew G. Knepley call DMSetType(dm, DMPLEX, ierr);CHKERRA(ierr) 3930602db0SMatthew G. Knepley call DMSetFromOptions(dm, ierr);CHKERRA(ierr) 40*2e3d3ef9SMatthew G. Knepley call DMViewFromOptions(dm, PETSC_NULL_VEC, '-dm_view', ierr);CHKERRA(ierr) 4130602db0SMatthew G. Knepley call DMGetDimension(dm, dim, ierr);CHKERRA(ierr) 42c4762a1bSJed Brown! Create a scalar field u, a vector field v, and a surface vector field w 43c4762a1bSJed Brown numFields = 3 44c4762a1bSJed Brown numComp(1) = 1 45c4762a1bSJed Brown numComp(2) = dim 46c4762a1bSJed Brown numComp(3) = dim-1 47c4762a1bSJed Brown pNumComp => numComp 48c4762a1bSJed Brown do i = 1, numFields*(dim+1) 49c4762a1bSJed Brown numDof(i) = 0 50c4762a1bSJed Brown end do 51c4762a1bSJed Brown! Let u be defined on vertices 52c4762a1bSJed Brown numDof(0*(dim+1)+1) = 1 53c4762a1bSJed Brown! Let v be defined on cells 54c4762a1bSJed Brown numDof(1*(dim+1)+dim+1) = dim 55c4762a1bSJed Brown! Let v be defined on faces 56c4762a1bSJed Brown numDof(2*(dim+1)+dim) = dim-1 57c4762a1bSJed Brown pNumDof => numDof 58c4762a1bSJed Brown! Setup boundary conditions 59c4762a1bSJed Brown numBC = 1 60c4762a1bSJed Brown! Test label retrieval 61c4762a1bSJed Brown call DMGetLabel(dm, 'marker', label, ierr);CHKERRA(ierr) 62c4762a1bSJed Brown call DMLabelGetValue(label, zero, val, ierr);CHKERRA(ierr) 63*2e3d3ef9SMatthew G. Knepley if (size .eq. 1 .and. val .ne. -1) then 64c4762a1bSJed Brown SETERRA(PETSC_COMM_SELF,PETSC_ERR_PLIB,'Error in library') 65c4762a1bSJed Brown endif 66c4762a1bSJed Brown call DMLabelGetValue(label, eight, val, ierr);CHKERRA(ierr) 67*2e3d3ef9SMatthew G. Knepley if (size .eq. 1 .and. val .ne. 1) then 68c4762a1bSJed Brown SETERRA(PETSC_COMM_SELF,PETSC_ERR_PLIB,'Error in library') 69c4762a1bSJed Brown endif 70c4762a1bSJed Brown! Prescribe a Dirichlet condition on u on the boundary 71c4762a1bSJed Brown! Label "marker" is made by the mesh creation routine 72c4762a1bSJed Brown bcField(1) = 0 73c4762a1bSJed Brown pBcField => bcField 74c4762a1bSJed Brown call ISCreateStride(PETSC_COMM_WORLD, one, zero, one, bcCompIS(1), ierr);CHKERRA(ierr) 75c4762a1bSJed Brown pBcCompIS => bcCompIS 76c4762a1bSJed Brown call DMGetStratumIS(dm, 'marker', one, bcPointIS(1),ierr);CHKERRA(ierr) 77c4762a1bSJed Brown pBcPointIS => bcPointIS 78c4762a1bSJed Brown! Create a PetscSection with this data layout 79c4762a1bSJed Brown call DMSetNumFields(dm, numFields,ierr);CHKERRA(ierr) 80c4762a1bSJed Brown call DMPlexCreateSection(dm,nolabel,pNumComp,pNumDof,numBC,pBcField,pBcCompIS,pBcPointIS,PETSC_NULL_IS,section,ierr) 81c4762a1bSJed Brown CHKERRA(ierr) 82c4762a1bSJed Brown call ISDestroy(bcCompIS(1), ierr);CHKERRA(ierr) 83c4762a1bSJed Brown call ISDestroy(bcPointIS(1), ierr);CHKERRA(ierr) 84c4762a1bSJed Brown! Name the Field variables 85c4762a1bSJed Brown call PetscSectionSetFieldName(section, zero, 'u', ierr);CHKERRA(ierr) 86c4762a1bSJed Brown call PetscSectionSetFieldName(section, one, 'v', ierr);CHKERRA(ierr) 87c4762a1bSJed Brown call PetscSectionSetFieldName(section, two, 'w', ierr);CHKERRA(ierr) 88*2e3d3ef9SMatthew G. Knepley if (size .eq. 1) then 89c4762a1bSJed Brown call PetscSectionView(section, PETSC_VIEWER_STDOUT_WORLD, ierr);CHKERRA(ierr) 90*2e3d3ef9SMatthew G. Knepley endif 91c4762a1bSJed Brown! Tell the DM to use this data layout 92c4762a1bSJed Brown call DMSetLocalSection(dm, section, ierr);CHKERRA(ierr) 93c4762a1bSJed Brown! Create a Vec with this layout and view it 94c4762a1bSJed Brown call DMGetGlobalVector(dm, u, ierr);CHKERRA(ierr) 95c4762a1bSJed Brown call PetscViewerCreate(PETSC_COMM_WORLD, viewer, ierr);CHKERRA(ierr) 96c4762a1bSJed Brown call PetscViewerSetType(viewer, PETSCVIEWERVTK, ierr);CHKERRA(ierr) 978ec8862eSJed Brown call PetscViewerFileSetName(viewer, 'sol.vtu', ierr);CHKERRA(ierr) 98c4762a1bSJed Brown call VecView(u, viewer, ierr);CHKERRA(ierr) 99c4762a1bSJed Brown call PetscViewerDestroy(viewer, ierr);CHKERRA(ierr) 100c4762a1bSJed Brown call DMRestoreGlobalVector(dm, u, ierr);CHKERRA(ierr) 101c4762a1bSJed Brown! Cleanup 102c4762a1bSJed Brown call PetscSectionDestroy(section, ierr);CHKERRA(ierr) 103c4762a1bSJed Brown call DMDestroy(dm, ierr);CHKERRA(ierr) 104c4762a1bSJed Brown 105c4762a1bSJed Brown call PetscFinalize(ierr) 106c4762a1bSJed Brown end program DMPlexTestField 107c4762a1bSJed Brown 108c4762a1bSJed Brown!/*TEST 109c4762a1bSJed Brown! build: 110dfd57a17SPierre Jolivet! requires: defined(PETSC_USING_F90FREEFORM) 111c4762a1bSJed Brown! 112c4762a1bSJed Brown! test: 113c4762a1bSJed Brown! suffix: 0 114c4762a1bSJed Brown! requires: triangle 11530602db0SMatthew G. Knepley! args: -info :~sys,mat: 116c4762a1bSJed Brown! 117c4762a1bSJed Brown! test: 118*2e3d3ef9SMatthew G. Knepley! suffix: 0_2 119*2e3d3ef9SMatthew G. Knepley! requires: triangle 120*2e3d3ef9SMatthew G. Knepley! nsize: 2 121*2e3d3ef9SMatthew G. Knepley! args: -info :~sys,mat,partitioner: 122*2e3d3ef9SMatthew G. Knepley! 123*2e3d3ef9SMatthew G. Knepley! test: 124c4762a1bSJed Brown! suffix: 1 125c4762a1bSJed Brown! requires: ctetgen 12630602db0SMatthew G. Knepley! args: -dm_plex_dim 3 -info :~sys,mat: 127c4762a1bSJed Brown! 128c4762a1bSJed Brown!TEST*/ 129