1#include "petsc/finclude/petsc.h" 2program ex98f90 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 MPIU_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_OBJECT, '-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(ISGetIndices(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(ISGetIndices(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(ISRestoreIndices(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(ISGetIndices(pointIS, pointID, ierr)) 89 do p = 1, numPoints 90 constraints(1) = mod(setID(s), sdim) 91 PetscCallA(PetscSectionSetConstraintIndices(section, pointID(p), constraints, ierr)) 92 PetscCallA(PetscSectionSetFieldConstraintIndices(section, pointID(p), 0_kPI, constraints, ierr)) 93 end do 94 PetscCallA(ISRestoreIndices(pointIS, pointID, ierr)) 95 PetscCallA(ISDestroy(pointIS, ierr)) 96 end do 97 deallocate (constraints) 98 PetscCallA(ISRestoreIndices(setIS, setID, ierr)) 99 PetscCallA(ISDestroy(setIS, ierr)) 100 PetscCallA(PetscObjectViewFromOptions(section, PETSC_NULL_OBJECT, '-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