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