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