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