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