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