Actual source code: ex98f90.F90

  1: program ex98f90
  2: #include "petsc/finclude/petsc.h"
  3:     use petsc
  4:     implicit none

  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)

 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

 24:     PetscCallA(PetscInitialize(ierr))
 25:     PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD,numProc,ierr))

 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 ')

 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))

 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))

 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))

 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

 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

 80:     PetscCallA(PetscSectionSetUp(section,ierr))

 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))

102:     PetscCallA(PetscSectionDestroy(section,ierr))
103:     PetscCallA(DMDestroy(dm,ierr))

105:     PetscCallA(PetscFinalize(ierr))
106: end program ex98f90

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*/