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 if (.not. flg) then 29 SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"missing input file name -i <input file name>") 30 end if 31 32 PetscCallA(DMPlexCreateFromFile(PETSC_COMM_WORLD,ifilename,PETSC_NULL_CHARACTER,PETSC_TRUE,dm,ierr)) 33 PetscCallA(DMPlexDistributeSetDefault(dm,PETSC_FALSE,ierr)) 34 PetscCallA(DMSetFromOptions(dm,ierr)) 35 36 if (numproc > 1) then 37 PetscCallA(DMPlexDistribute(dm,0_kPI,PETSC_NULL_SF,pdm,ierr)) 38 PetscCallA(DMDestroy(dm,ierr)) 39 dm = pdm 40 end if 41 PetscCallA(DMViewFromOptions(dm,PETSC_NULL_OPTIONS,"-dm_view",ierr)) 42 43 PetscCallA(DMGetDimension(dm,sdim,ierr)) 44 PetscCallA(PetscObjectGetComm(dm,comm,ierr)) 45 PetscCallA(PetscSectionCreate(comm,section,ierr)) 46 PetscCallA(PetscSectionSetNumFields(section,1_kPI,ierr)) 47 PetscCallA(PetscSectionSetFieldName(section,0_kPI,"U",ierr)) 48 PetscCallA(PetscSectionSetFieldComponents(section,0_kPI,sdim,ierr)) 49 PetscCallA(DMPlexGetChart(dm,pStart,pEnd,ierr)) 50 PetscCallA(PetscSectionSetChart(section,pStart,pEnd,ierr)) 51 52 ! initialize the section storage for a P1 field 53 PetscCallA(DMPlexGetDepthStratum(dm,0_kPI,pStart,pEnd,ierr)) 54 do p = pStart,pEnd-1 55 PetscCallA(PetscSectionSetDof(section,p,sdim,ierr)) 56 PetscCallA(PetscSectionSetFieldDof(section,p,0_kPI,sdim,ierr)) 57 end do 58 59 ! add constraints at all vertices belonging to a vertex set: 60 ! first pass is to reserve space 61 PetscCallA(DMGetLabelSize(dm,"Vertex Sets",numVS,ierr)) 62 write(iobuffer,'("# Vertex set: ",i3,"\n")' ) numVS 63 PetscCallA(PetscPrintf(PETSC_COMM_WORLD,iobuffer,ierr)) 64 PetscCallA(DMGetLabelIdIS(dm,"Vertex Sets",setIS,ierr)) 65 PetscCallA(ISGetIndicesF90(setIS,setID,ierr)) 66 do s = 1,numVS 67 PetscCallA(DMGetStratumIS(dm,"Vertex Sets",setID(s),pointIS,ierr)) 68 PetscCallA(DMGetStratumSize(dm,"Vertex Sets",setID(s),numPoints,ierr)) 69 write(iobuffer,'("set ",i3," size ",i3,"\n")' ) s,numPoints 70 PetscCallA(PetscPrintf(PETSC_COMM_WORLD,iobuffer,ierr)) 71 PetscCallA(ISGetIndicesF90(pointIS,pointID,ierr)) 72 do p = 1,numPoints 73 write(iobuffer,'(" point ",i3,"\n")' ) pointID(p) 74 PetscCallA(PetscPrintf(PETSC_COMM_WORLD,iobuffer,ierr)) 75 PetscCallA(PetscSectionSetConstraintDof(section,pointID(p),1_kPI,ierr)) 76 PetscCallA(PetscSectionSetFieldConstraintDof(section,pointID(p),0_kPI,1_kPI,ierr)) 77 end do 78 PetscCallA(ISRestoreIndicesF90(pointIS,pointID,ierr)) 79 PetscCallA(ISDestroy(pointIS,ierr)) 80 end do 81 82 PetscCallA(PetscSectionSetUp(section,ierr)) 83 84 ! add constraints at all vertices belonging to a vertex set: 85 ! second pass is to assign constraints to a specific component / dof 86 allocate(constraints(1)) 87 do s = 1,numVS 88 PetscCallA(DMGetStratumIS(dm,"Vertex Sets",setID(s),pointIS,ierr)) 89 PetscCallA(DMGetStratumSize(dm,"Vertex Sets",setID(s),numPoints,ierr)) 90 PetscCallA(ISGetIndicesF90(pointIS,pointID,ierr)) 91 do p = 1,numPoints 92 constraints(1) = mod(setID(s),sdim) 93 PetscCallA(PetscSectionSetConstraintIndicesF90(section,pointID(p),constraints,ierr)) 94 PetscCallA(PetscSectionSetFieldConstraintIndicesF90(section,pointID(p),0_kPI,constraints,ierr)) 95 end do 96 PetscCallA(ISRestoreIndicesF90(pointIS,pointID,ierr)) 97 PetscCallA(ISDestroy(pointIS,ierr)) 98 end do 99 deallocate(constraints) 100 PetscCallA(ISRestoreIndicesF90(setIS,setID,ierr)) 101 PetscCallA(ISDestroy(setIS,ierr)) 102 PetscCallA(PetscObjectViewFromOptions(section,PETSC_NULL_SECTION,"-dm_section_view",ierr)) 103 104 PetscCallA(PetscSectionDestroy(section,ierr)) 105 PetscCallA(DMDestroy(dm,ierr)) 106 107 PetscCallA(PetscFinalize(ierr)) 108end program ex98f90 109 110! /*TEST 111! build: 112! requires: exodusii pnetcdf !complex 113! testset: 114! args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/SquareFaceSet.exo -dm_view -dm_section_view 115! nsize: 1 116! test: 117! suffix: 0 118! args: 119! TEST*/ 120