Actual source code: ex1f90.F90
petsc-3.14.6 2021-03-30
1: program DMPlexTestField
2: #include "petsc/finclude/petscdmplex.h"
3: #include "petsc/finclude/petscdmlabel.h"
4: use petscdmplex
5: use petscsys
6: implicit none
8: DM :: dm
9: DMLabel :: label
10: Vec :: u
11: PetscViewer :: viewer
12: PetscSection :: section
13: PetscInt :: dim,numFields,numBC
14: PetscInt :: i,val
15: DMLabel, pointer :: nolabel(:) => NULL()
16: PetscInt, target, dimension(3) :: numComp
17: PetscInt, pointer :: pNumComp(:)
18: PetscInt, target, dimension(12) :: numDof
19: PetscInt, pointer :: pNumDof(:)
20: PetscInt, target, dimension(1) :: bcField
21: PetscInt, pointer :: pBcField(:)
22: PetscInt, parameter :: zero = 0, one = 1, two = 2, eight = 8
23: IS, target, dimension(1) :: bcCompIS
24: IS, target, dimension(1) :: bcPointIS
25: IS, pointer :: pBcCompIS(:)
26: IS, pointer :: pBcPointIS(:)
27: PetscBool :: interpolate,flg
28: PetscErrorCode :: ierr
30: call PetscInitialize(PETSC_NULL_CHARACTER, ierr)
31: if (ierr .ne. 0) then
32: print*,'Unable to initialize PETSc'
33: stop
34: endif
35: dim = 2
36: call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-dim', dim,flg, ierr);CHKERRA(ierr)
37: interpolate = PETSC_TRUE
38: ! Create a mesh
39: call DMPlexCreateBoxMesh(PETSC_COMM_WORLD, dim, PETSC_TRUE, PETSC_NULL_INTEGER, PETSC_NULL_REAL, PETSC_NULL_REAL, PETSC_NULL_INTEGER, interpolate, dm, ierr);CHKERRA(ierr)
40: ! Create a scalar field u, a vector field v, and a surface vector field w
41: numFields = 3
42: numComp(1) = 1
43: numComp(2) = dim
44: numComp(3) = dim-1
45: pNumComp => numComp
46: do i = 1, numFields*(dim+1)
47: numDof(i) = 0
48: end do
49: ! Let u be defined on vertices
50: numDof(0*(dim+1)+1) = 1
51: ! Let v be defined on cells
52: numDof(1*(dim+1)+dim+1) = dim
53: ! Let v be defined on faces
54: numDof(2*(dim+1)+dim) = dim-1
55: pNumDof => numDof
56: ! Setup boundary conditions
57: numBC = 1
58: ! Test label retrieval
59: call DMGetLabel(dm, 'marker', label, ierr);CHKERRA(ierr)
60: call DMLabelGetValue(label, zero, val, ierr);CHKERRA(ierr)
61: if (val .ne. -1) then
62: SETERRA(PETSC_COMM_SELF,PETSC_ERR_PLIB,'Error in library')
63: endif
64: call DMLabelGetValue(label, eight, val, ierr);CHKERRA(ierr)
65: if (val .ne. 1) then
66: SETERRA(PETSC_COMM_SELF,PETSC_ERR_PLIB,'Error in library')
67: endif
68: ! Prescribe a Dirichlet condition on u on the boundary
69: ! Label "marker" is made by the mesh creation routine
70: bcField(1) = 0
71: pBcField => bcField
72: call ISCreateStride(PETSC_COMM_WORLD, one, zero, one, bcCompIS(1), ierr);CHKERRA(ierr)
73: pBcCompIS => bcCompIS
74: call DMGetStratumIS(dm, 'marker', one, bcPointIS(1),ierr);CHKERRA(ierr)
75: pBcPointIS => bcPointIS
76: ! Create a PetscSection with this data layout
77: call DMSetNumFields(dm, numFields,ierr);CHKERRA(ierr)
78: call DMPlexCreateSection(dm,nolabel,pNumComp,pNumDof,numBC,pBcField,pBcCompIS,pBcPointIS,PETSC_NULL_IS,section,ierr)
79: CHKERRA(ierr)
80: call ISDestroy(bcCompIS(1), ierr);CHKERRA(ierr)
81: call ISDestroy(bcPointIS(1), ierr);CHKERRA(ierr)
82: ! Name the Field variables
83: call PetscSectionSetFieldName(section, zero, 'u', ierr);CHKERRA(ierr)
84: call PetscSectionSetFieldName(section, one, 'v', ierr);CHKERRA(ierr)
85: call PetscSectionSetFieldName(section, two, 'w', ierr);CHKERRA(ierr)
86: call PetscSectionView(section, PETSC_VIEWER_STDOUT_WORLD, ierr);CHKERRA(ierr)
87: ! Tell the DM to use this data layout
88: call DMSetLocalSection(dm, section, ierr);CHKERRA(ierr)
89: ! Create a Vec with this layout and view it
90: call DMGetGlobalVector(dm, u, ierr);CHKERRA(ierr)
91: call PetscViewerCreate(PETSC_COMM_WORLD, viewer, ierr);CHKERRA(ierr)
92: call PetscViewerSetType(viewer, PETSCVIEWERVTK, ierr);CHKERRA(ierr)
93: call PetscViewerFileSetName(viewer, 'sol.vtu', ierr);CHKERRA(ierr)
94: call VecView(u, viewer, ierr);CHKERRA(ierr)
95: call PetscViewerDestroy(viewer, ierr);CHKERRA(ierr)
96: call DMRestoreGlobalVector(dm, u, ierr);CHKERRA(ierr)
97: ! Cleanup
98: call PetscSectionDestroy(section, ierr);CHKERRA(ierr)
99: call DMDestroy(dm, ierr);CHKERRA(ierr)
101: call PetscFinalize(ierr)
102: end program DMPlexTestField
104: !/*TEST
105: ! build:
106: ! requires: define(PETSC_USING_F90FREEFORM)
107: !
108: ! test:
109: ! suffix: 0
110: ! requires: triangle
111: ! args: -info :~sys:
112: !
113: ! test:
114: ! suffix: 1
115: ! requires: ctetgen
116: ! args: -dim 3 -info :~sys:
117: !
118: !TEST*/