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 if (size .eq. 1 .and. val .ne. -1) then 60 SETERRA(PETSC_COMM_SELF,PETSC_ERR_PLIB,'Error in library') 61 endif 62 PetscCallA(DMLabelGetValue(label, eight, val, 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! Prescribe a Dirichlet condition on u on the boundary 67! Label "marker" is made by the mesh creation routine 68 bcField(1) = 0 69 pBcField => bcField 70 PetscCallA(ISCreateStride(PETSC_COMM_WORLD, one, zero, one, bcCompIS(1), ierr)) 71 pBcCompIS => bcCompIS 72 PetscCallA(DMGetStratumIS(dm, 'marker', one, bcPointIS(1),ierr)) 73 pBcPointIS => bcPointIS 74! Create a PetscSection with this data layout 75 PetscCallA(DMSetNumFields(dm, numFields,ierr)) 76 PetscCallA(DMPlexCreateSection(dm,nolabel,pNumComp,pNumDof,numBC,pBcField,pBcCompIS,pBcPointIS,PETSC_NULL_IS,section,ierr)) 77 PetscCallA(ISDestroy(bcCompIS(1), ierr)) 78 PetscCallA(ISDestroy(bcPointIS(1), ierr)) 79! Name the Field variables 80 PetscCallA(PetscSectionSetFieldName(section, zero, 'u', ierr)) 81 PetscCallA(PetscSectionSetFieldName(section, one, 'v', ierr)) 82 PetscCallA(PetscSectionSetFieldName(section, two, 'w', ierr)) 83 if (size .eq. 1) then 84 PetscCallA(PetscSectionView(section, PETSC_VIEWER_STDOUT_WORLD, ierr)) 85 endif 86! Tell the DM to use this data layout 87 PetscCallA(DMSetLocalSection(dm, section, ierr)) 88! Create a Vec with this layout and view it 89 PetscCallA(DMGetGlobalVector(dm, u, ierr)) 90 PetscCallA(PetscViewerCreate(PETSC_COMM_WORLD, viewer, ierr)) 91 PetscCallA(PetscViewerSetType(viewer, PETSCVIEWERVTK, ierr)) 92 PetscCallA(PetscViewerFileSetName(viewer, 'sol.vtu', ierr)) 93 PetscCallA(VecView(u, viewer, ierr)) 94 PetscCallA(PetscViewerDestroy(viewer, ierr)) 95 PetscCallA(DMRestoreGlobalVector(dm, u, ierr)) 96! Cleanup 97 PetscCallA(PetscSectionDestroy(section, ierr)) 98 PetscCallA(DMDestroy(dm, ierr)) 99 100 PetscCallA(PetscFinalize(ierr)) 101 end program DMPlexTestField 102 103!/*TEST 104! build: 105! requires: defined(PETSC_USING_F90FREEFORM) 106! 107! test: 108! suffix: 0 109! requires: triangle 110! args: -info :~sys,mat: 111! 112! test: 113! suffix: 0_2 114! requires: triangle 115! nsize: 2 116! args: -info :~sys,mat,partitioner: 117! 118! test: 119! suffix: 1 120! requires: ctetgen 121! args: -dm_plex_dim 3 -info :~sys,mat: 122! 123!TEST*/ 124