1 program DMPlexTestField 2#include "petsc/finclude/petscdmplex.h" 3#include "petsc/finclude/petscdmlabel.h" 4 use petscdmplex 5 use petscsys 6 implicit none 7 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 29 30 PetscCallA(PetscInitialize(ierr)) 31 PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD, size, ierr)) 32! Create a mesh 33 PetscCallA(DMCreate(PETSC_COMM_WORLD, dm, ierr)) 34 PetscCallA(DMSetType(dm, DMPLEX, ierr)) 35 PetscCallA(DMSetFromOptions(dm, ierr)) 36 PetscCallA(DMViewFromOptions(dm, PETSC_NULL_VEC, '-dm_view', ierr)) 37 PetscCallA(DMGetDimension(dm, dim, ierr)) 38! Create a scalar field u, a vector field v, and a surface vector field w 39 numFields = 3 40 numComp(1) = 1 41 numComp(2) = dim 42 numComp(3) = dim-1 43 pNumComp => numComp 44 do i = 1, numFields*(dim+1) 45 numDof(i) = 0 46 end do 47! Let u be defined on vertices 48 numDof(0*(dim+1)+1) = 1 49! Let v be defined on cells 50 numDof(1*(dim+1)+dim+1) = dim 51! Let v be defined on faces 52 numDof(2*(dim+1)+dim) = dim-1 53 pNumDof => numDof 54! Setup boundary conditions 55 numBC = 1 56! Test label retrieval 57 PetscCallA(DMGetLabel(dm, 'marker', label, ierr)) 58 PetscCallA(DMLabelGetValue(label, zero, val, ierr)) 59 PetscCheckA(size .ne. 1 .or. val .eq. -1,PETSC_COMM_SELF,PETSC_ERR_PLIB,'Error in library') 60 PetscCallA(DMLabelGetValue(label, eight, val, ierr)) 61 PetscCheckA(size .ne. 1 .or. val .eq. 1,PETSC_COMM_SELF,PETSC_ERR_PLIB,'Error in library') 62! Prescribe a Dirichlet condition on u on the boundary 63! Label "marker" is made by the mesh creation routine 64 bcField(1) = 0 65 pBcField => bcField 66 PetscCallA(ISCreateStride(PETSC_COMM_WORLD, one, zero, one, bcCompIS(1), ierr)) 67 pBcCompIS => bcCompIS 68 PetscCallA(DMGetStratumIS(dm, 'marker', one, bcPointIS(1),ierr)) 69 pBcPointIS => bcPointIS 70! Create a PetscSection with this data layout 71 PetscCallA(DMSetNumFields(dm, numFields,ierr)) 72 PetscCallA(DMPlexCreateSection(dm,nolabel,pNumComp,pNumDof,numBC,pBcField,pBcCompIS,pBcPointIS,PETSC_NULL_IS,section,ierr)) 73 PetscCallA(ISDestroy(bcCompIS(1), ierr)) 74 PetscCallA(ISDestroy(bcPointIS(1), ierr)) 75! Name the Field variables 76 PetscCallA(PetscSectionSetFieldName(section, zero, 'u', ierr)) 77 PetscCallA(PetscSectionSetFieldName(section, one, 'v', ierr)) 78 PetscCallA(PetscSectionSetFieldName(section, two, 'w', ierr)) 79 if (size .eq. 1) then 80 PetscCallA(PetscSectionView(section, PETSC_VIEWER_STDOUT_WORLD, ierr)) 81 endif 82! Tell the DM to use this data layout 83 PetscCallA(DMSetLocalSection(dm, section, ierr)) 84! Create a Vec with this layout and view it 85 PetscCallA(DMGetGlobalVector(dm, u, ierr)) 86 PetscCallA(PetscViewerCreate(PETSC_COMM_WORLD, viewer, ierr)) 87 PetscCallA(PetscViewerSetType(viewer, PETSCVIEWERVTK, ierr)) 88 PetscCallA(PetscViewerFileSetName(viewer, 'sol.vtu', ierr)) 89 PetscCallA(VecView(u, viewer, ierr)) 90 PetscCallA(PetscViewerDestroy(viewer, ierr)) 91 PetscCallA(DMRestoreGlobalVector(dm, u, ierr)) 92! Cleanup 93 PetscCallA(PetscSectionDestroy(section, ierr)) 94 PetscCallA(DMDestroy(dm, ierr)) 95 96 PetscCallA(PetscFinalize(ierr)) 97 end program DMPlexTestField 98 99!/*TEST 100! build: 101! requires: defined(PETSC_USING_F90FREEFORM) 102! 103! test: 104! suffix: 0 105! requires: triangle 106! args: -info :~sys,mat: 107! 108! test: 109! suffix: 0_2 110! requires: triangle 111! nsize: 2 112! args: -info :~sys,mat,partitioner: 113! 114! test: 115! suffix: 1 116! requires: ctetgen 117! args: -dm_plex_dim 3 -info :~sys,mat: 118! 119!TEST*/ 120