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