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