Actual source code: ex1f90.F90
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: PetscMPIInt :: size
24: IS, target, dimension(1) :: bcCompIS
25: IS, target, dimension(1) :: bcPointIS
26: IS, pointer :: pBcCompIS(:)
27: IS, pointer :: pBcPointIS(:)
28: PetscErrorCode :: ierr
30: PetscCallA(PetscInitialize(ierr))
31: PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD, size, ierr))
32: ! Create a mesh
33: PetscCallA(DMCreate(PETSC_COMM_WORLD, dm, ierr))
34: PetscCallA(DMSetType(dm, DMPLEX, ierr))
35: PetscCallA(DMSetFromOptions(dm, ierr))
36: PetscCallA(DMViewFromOptions(dm, PETSC_NULL_VEC, '-dm_view', ierr))
37: PetscCallA(DMGetDimension(dm, dim, ierr))
38: ! Create a scalar field u, a vector field v, and a surface vector field w
39: numFields = 3
40: numComp(1) = 1
41: numComp(2) = dim
42: numComp(3) = dim-1
43: pNumComp => numComp
44: do i = 1, numFields*(dim+1)
45: numDof(i) = 0
46: end do
47: ! Let u be defined on vertices
48: numDof(0*(dim+1)+1) = 1
49: ! Let v be defined on cells
50: numDof(1*(dim+1)+dim+1) = dim
51: ! Let v be defined on faces
52: numDof(2*(dim+1)+dim) = dim-1
53: pNumDof => numDof
54: ! Setup boundary conditions
55: numBC = 1
56: ! Test label retrieval
57: PetscCallA(DMGetLabel(dm, 'marker', label, ierr))
58: PetscCallA(DMLabelGetValue(label, zero, val, ierr))
59: if (size .eq. 1 .and. val .ne. -1) then
60: SETERRA(PETSC_COMM_SELF,PETSC_ERR_PLIB,'Error in library')
61: endif
62: PetscCallA(DMLabelGetValue(label, eight, val, ierr))
63: if (size .eq. 1 .and. val .ne. 1) then
64: SETERRA(PETSC_COMM_SELF,PETSC_ERR_PLIB,'Error in library')
65: endif
66: ! Prescribe a Dirichlet condition on u on the boundary
67: ! Label "marker" is made by the mesh creation routine
68: bcField(1) = 0
69: pBcField => bcField
70: PetscCallA(ISCreateStride(PETSC_COMM_WORLD, one, zero, one, bcCompIS(1), ierr))
71: pBcCompIS => bcCompIS
72: PetscCallA(DMGetStratumIS(dm, 'marker', one, bcPointIS(1),ierr))
73: pBcPointIS => bcPointIS
74: ! Create a PetscSection with this data layout
75: PetscCallA(DMSetNumFields(dm, numFields,ierr))
76: PetscCallA(DMPlexCreateSection(dm,nolabel,pNumComp,pNumDof,numBC,pBcField,pBcCompIS,pBcPointIS,PETSC_NULL_IS,section,ierr))
77: PetscCallA(ISDestroy(bcCompIS(1), ierr))
78: PetscCallA(ISDestroy(bcPointIS(1), ierr))
79: ! Name the Field variables
80: PetscCallA(PetscSectionSetFieldName(section, zero, 'u', ierr))
81: PetscCallA(PetscSectionSetFieldName(section, one, 'v', ierr))
82: PetscCallA(PetscSectionSetFieldName(section, two, 'w', ierr))
83: if (size .eq. 1) then
84: PetscCallA(PetscSectionView(section, PETSC_VIEWER_STDOUT_WORLD, ierr))
85: endif
86: ! Tell the DM to use this data layout
87: PetscCallA(DMSetLocalSection(dm, section, ierr))
88: ! Create a Vec with this layout and view it
89: PetscCallA(DMGetGlobalVector(dm, u, ierr))
90: PetscCallA(PetscViewerCreate(PETSC_COMM_WORLD, viewer, ierr))
91: PetscCallA(PetscViewerSetType(viewer, PETSCVIEWERVTK, ierr))
92: PetscCallA(PetscViewerFileSetName(viewer, 'sol.vtu', ierr))
93: PetscCallA(VecView(u, viewer, ierr))
94: PetscCallA(PetscViewerDestroy(viewer, ierr))
95: PetscCallA(DMRestoreGlobalVector(dm, u, ierr))
96: ! Cleanup
97: PetscCallA(PetscSectionDestroy(section, ierr))
98: PetscCallA(DMDestroy(dm, ierr))
100: PetscCallA(PetscFinalize(ierr))
101: end program DMPlexTestField
103: !/*TEST
104: ! build:
105: ! requires: defined(PETSC_USING_F90FREEFORM)
106: !
107: ! test:
108: ! suffix: 0
109: ! requires: triangle
110: ! args: -info :~sys,mat:
111: !
112: ! test:
113: ! suffix: 0_2
114: ! requires: triangle
115: ! nsize: 2
116: ! args: -info :~sys,mat,partitioner:
117: !
118: ! test:
119: ! suffix: 1
120: ! requires: ctetgen
121: ! args: -dm_plex_dim 3 -info :~sys,mat:
122: !
123: !TEST*/