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