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:     if (.not. flg) then
 29:         SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"missing input file name -i <input file name>")
 30:     end if

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

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

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

 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

 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

 82:     PetscCallA(PetscSectionSetUp(section,ierr))

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

104:     PetscCallA(PetscSectionDestroy(section,ierr))
105:     PetscCallA(DMDestroy(dm,ierr))

107:     PetscCallA(PetscFinalize(ierr))
108: end program ex98f90

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