1program ex98f90 2#include "petsc/finclude/petsc.h" 3 use petsc 4 implicit none 5 6 ! Get the fortran kind associated with PetscInt and PetscReal so that we can use literal constants. 7 PetscInt :: dummyPetscInt 8 PetscReal :: dummyPetscreal 9 integer,parameter :: kPI = kind(dummyPetscInt) 10 integer,parameter :: kPR = kind(dummyPetscReal) 11 12 type(tDM) :: dm,pdm 13 type(tPetscSection) :: section 14 character(len=PETSC_MAX_PATH_LEN) :: ifilename,iobuffer 15 PetscInt :: sdim,s,pStart,pEnd,p,numVS,numPoints 16 PetscInt,dimension(:),pointer :: constraints 17 type(tIS) :: setIS,pointIS 18 PetscInt,dimension(:),pointer :: setID,pointID 19 PetscErrorCode :: ierr 20 PetscBool :: flg 21 PetscMPIInt :: numProc 22 MPI_Comm :: comm 23 24 PetscCallA(PetscInitialize(ierr)) 25 PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD,numProc,ierr)) 26 27 PetscCallA(PetscOptionsGetString(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-i',ifilename,flg,ierr)) 28 PetscCheckA(flg,PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,'missing input file name -i <input file name>') 29 30 PetscCallA(DMPlexCreateFromFile(PETSC_COMM_WORLD,ifilename,PETSC_NULL_CHARACTER,PETSC_TRUE,dm,ierr)) 31 PetscCallA(DMPlexDistributeSetDefault(dm,PETSC_FALSE,ierr)) 32 PetscCallA(DMSetFromOptions(dm,ierr)) 33 34 if (numproc > 1) then 35 PetscCallA(DMPlexDistribute(dm,0_kPI,PETSC_NULL_SF,pdm,ierr)) 36 PetscCallA(DMDestroy(dm,ierr)) 37 dm = pdm 38 end if 39 PetscCallA(DMViewFromOptions(dm,PETSC_NULL_OPTIONS,'-dm_view',ierr)) 40 41 PetscCallA(DMGetDimension(dm,sdim,ierr)) 42 PetscCallA(PetscObjectGetComm(dm,comm,ierr)) 43 PetscCallA(PetscSectionCreate(comm,section,ierr)) 44 PetscCallA(PetscSectionSetNumFields(section,1_kPI,ierr)) 45 PetscCallA(PetscSectionSetFieldName(section,0_kPI,'U',ierr)) 46 PetscCallA(PetscSectionSetFieldComponents(section,0_kPI,sdim,ierr)) 47 PetscCallA(DMPlexGetChart(dm,pStart,pEnd,ierr)) 48 PetscCallA(PetscSectionSetChart(section,pStart,pEnd,ierr)) 49 50 ! initialize the section storage for a P1 field 51 PetscCallA(DMPlexGetDepthStratum(dm,0_kPI,pStart,pEnd,ierr)) 52 do p = pStart,pEnd-1 53 PetscCallA(PetscSectionSetDof(section,p,sdim,ierr)) 54 PetscCallA(PetscSectionSetFieldDof(section,p,0_kPI,sdim,ierr)) 55 end do 56 57 ! add constraints at all vertices belonging to a vertex set: 58 ! first pass is to reserve space 59 PetscCallA(DMGetLabelSize(dm,'Vertex Sets',numVS,ierr)) 60 write(iobuffer,'("# Vertex set: ",i3,"\n")' ) numVS 61 PetscCallA(PetscPrintf(PETSC_COMM_WORLD,iobuffer,ierr)) 62 PetscCallA(DMGetLabelIdIS(dm,'Vertex Sets',setIS,ierr)) 63 PetscCallA(ISGetIndicesF90(setIS,setID,ierr)) 64 do s = 1,numVS 65 PetscCallA(DMGetStratumIS(dm,'Vertex Sets',setID(s),pointIS,ierr)) 66 PetscCallA(DMGetStratumSize(dm,'Vertex Sets',setID(s),numPoints,ierr)) 67 write(iobuffer,'("set ",i3," size ",i3,"\n")' ) s,numPoints 68 PetscCallA(PetscPrintf(PETSC_COMM_WORLD,iobuffer,ierr)) 69 PetscCallA(ISGetIndicesF90(pointIS,pointID,ierr)) 70 do p = 1,numPoints 71 write(iobuffer,'(" point ",i3,"\n")' ) pointID(p) 72 PetscCallA(PetscPrintf(PETSC_COMM_WORLD,iobuffer,ierr)) 73 PetscCallA(PetscSectionSetConstraintDof(section,pointID(p),1_kPI,ierr)) 74 PetscCallA(PetscSectionSetFieldConstraintDof(section,pointID(p),0_kPI,1_kPI,ierr)) 75 end do 76 PetscCallA(ISRestoreIndicesF90(pointIS,pointID,ierr)) 77 PetscCallA(ISDestroy(pointIS,ierr)) 78 end do 79 80 PetscCallA(PetscSectionSetUp(section,ierr)) 81 82 ! add constraints at all vertices belonging to a vertex set: 83 ! second pass is to assign constraints to a specific component / dof 84 allocate(constraints(1)) 85 do s = 1,numVS 86 PetscCallA(DMGetStratumIS(dm,'Vertex Sets',setID(s),pointIS,ierr)) 87 PetscCallA(DMGetStratumSize(dm,'Vertex Sets',setID(s),numPoints,ierr)) 88 PetscCallA(ISGetIndicesF90(pointIS,pointID,ierr)) 89 do p = 1,numPoints 90 constraints(1) = mod(setID(s),sdim) 91 PetscCallA(PetscSectionSetConstraintIndicesF90(section,pointID(p),constraints,ierr)) 92 PetscCallA(PetscSectionSetFieldConstraintIndicesF90(section,pointID(p),0_kPI,constraints,ierr)) 93 end do 94 PetscCallA(ISRestoreIndicesF90(pointIS,pointID,ierr)) 95 PetscCallA(ISDestroy(pointIS,ierr)) 96 end do 97 deallocate(constraints) 98 PetscCallA(ISRestoreIndicesF90(setIS,setID,ierr)) 99 PetscCallA(ISDestroy(setIS,ierr)) 100 PetscCallA(PetscObjectViewFromOptions(section,PETSC_NULL_SECTION,'-dm_section_view',ierr)) 101 102 PetscCallA(PetscSectionDestroy(section,ierr)) 103 PetscCallA(DMDestroy(dm,ierr)) 104 105 PetscCallA(PetscFinalize(ierr)) 106end program ex98f90 107 108! /*TEST 109! build: 110! requires: exodusii pnetcdf !complex 111! testset: 112! args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/SquareFaceSet.exo -dm_view -dm_section_view 113! nsize: 1 114! test: 115! suffix: 0 116! args: 117! TEST*/ 118