Actual source code: ex2f90.F90
1: program main
2: #include <petsc/finclude/petscdmplex.h>
3: use petscdmplex
4: use petscsys
5: implicit none
7: DM dm
8: PetscInt, target, dimension(3) :: EC
9: PetscInt, target, dimension(2) :: VE
10: PetscInt, pointer :: pEC(:), pVE(:)
11: PetscInt, pointer :: nClosure(:)
12: PetscInt, pointer :: nJoin(:)
13: PetscInt, pointer :: nMeet(:)
14: PetscInt dim, cell, size
15: PetscInt i0,i1,i2,i3,i6,i7
16: PetscInt i8,i9,i10,i11
17: PetscErrorCode ierr
19: i0 = 0
20: i1 = 1
21: i2 = 2
22: i3 = 3
23: i6 = 6
24: i7 = 7
25: i8 = 8
26: i9 = 9
27: i10 = 10
28: i11 = 11
30: PetscCallA(PetscInitialize(ierr))
32: PetscCallA(DMPlexCreate(PETSC_COMM_WORLD, dm, ierr))
33: PetscCallA(PetscObjectSetName(dm, 'Mesh', ierr))
34: dim = 2
35: PetscCallA(DMSetDimension(dm, dim, ierr))
37: ! Make Doublet Mesh from Fig 2 of Flexible Representation of Computational Meshes,
38: ! except indexing is from 0 instead of 1 and we obey the new restrictions on
39: ! numbering: cells, vertices, faces, edges
40: PetscCallA(DMPlexSetChart(dm, i0, i11, ierr))
41: ! cells
42: PetscCallA(DMPlexSetConeSize(dm, i0, i3, ierr))
43: PetscCallA(DMPlexSetConeSize(dm, i1, i3, ierr))
44: ! edges
45: PetscCallA(DMPlexSetConeSize(dm, i6, i2, ierr))
46: PetscCallA(DMPlexSetConeSize(dm, i7, i2, ierr))
47: PetscCallA(DMPlexSetConeSize(dm, i8, i2, ierr))
48: PetscCallA(DMPlexSetConeSize(dm, i9, i2, ierr))
49: PetscCallA(DMPlexSetConeSize(dm, i10, i2, ierr))
51: PetscCallA(DMSetUp(dm, ierr))
53: EC(1) = 6
54: EC(2) = 7
55: EC(3) = 8
56: pEC => EC
57: PetscCallA(DMPlexSetCone(dm, i0, pEC, ierr))
58: EC(1) = 7
59: EC(2) = 9
60: EC(3) = 10
61: pEC => EC
62: PetscCallA(DMPlexSetCone(dm, i1 , pEC, ierr))
64: VE(1) = 2
65: VE(2) = 3
66: pVE => VE
67: PetscCallA(DMPlexSetCone(dm, i6 , pVE, ierr))
68: VE(1) = 3
69: VE(2) = 4
70: pVE => VE
71: PetscCallA(DMPlexSetCone(dm, i7 , pVE, ierr))
72: VE(1) = 4
73: VE(2) = 2
74: pVE => VE
75: PetscCallA(DMPlexSetCone(dm, i8 , pVE, ierr))
76: VE(1) = 3
77: VE(2) = 5
78: pVE => VE
79: PetscCallA(DMPlexSetCone(dm, i9 , pVE, ierr))
80: VE(1) = 5
81: VE(2) = 4
82: pVE => VE
83: PetscCallA(DMPlexSetCone(dm, i10 , pVE, ierr))
85: PetscCallA(DMPlexSymmetrize(dm,ierr))
86: PetscCallA(DMPlexStratify(dm,ierr))
87: PetscCallA(DMView(dm, PETSC_VIEWER_STDOUT_WORLD, ierr))
89: ! Test Closure
90: do cell = 0,1
91: PetscCallA(DMPlexGetTransitiveClosure(dm,cell,PETSC_TRUE,nClosure,ierr))
92: ! Different Fortran compilers print a different number of columns
93: ! per row producing different outputs in the test runs hence
94: ! do not print the nClosure
95: write(*,1000) 'nClosure ',nClosure
96: 1000 format (a,30i4)
97: PetscCallA(DMPlexRestoreTransitiveClosure(dm,cell,PETSC_TRUE,nClosure,ierr))
98: end do
100: ! Test Join
101: size = 2
102: VE(1) = 6
103: VE(2) = 7
104: pVE => VE
105: PetscCallA(DMPlexGetJoin(dm, size, pVE, nJoin, ierr))
106: write(*,1001) 'Join of',pVE
107: write(*,1002) ' is',nJoin
108: PetscCallA(DMPlexRestoreJoin(dm, size, pVE, nJoin, ierr))
109: size = 2
110: VE(1) = 9
111: VE(2) = 7
112: pVE => VE
113: PetscCallA(DMPlexGetJoin(dm, size, pVE, nJoin, ierr))
114: write(*,1001) 'Join of',pVE
115: 1001 format (a,10i5)
116: write(*,1002) ' is',nJoin
117: 1002 format (a,10i5)
118: PetscCallA(DMPlexRestoreJoin(dm, size, pVE, nJoin, ierr))
119: ! Test Full Join
120: size = 3
121: EC(1) = 3
122: EC(2) = 4
123: EC(3) = 5
124: pEC => EC
125: PetscCallA(DMPlexGetFullJoin(dm, size, pEC, nJoin, ierr))
126: write(*,1001) 'Full Join of',pEC
127: write(*,1002) ' is',nJoin
128: PetscCallA(DMPlexRestoreJoin(dm, size, pEC, nJoin, ierr))
129: ! Test Meet
130: size = 2
131: VE(1) = 0
132: VE(2) = 1
133: pVE => VE
134: PetscCallA(DMPlexGetMeet(dm, size, pVE, nMeet, ierr))
135: write(*,1001) 'Meet of',pVE
136: write(*,1002) ' is',nMeet
137: PetscCallA(DMPlexRestoreMeet(dm, size, pVE, nMeet, ierr))
138: size = 2
139: VE(1) = 6
140: VE(2) = 7
141: pVE => VE
142: PetscCallA(DMPlexGetMeet(dm, size, pVE, nMeet, ierr))
143: write(*,1001) 'Meet of',pVE
144: write(*,1002) ' is',nMeet
145: PetscCallA(DMPlexRestoreMeet(dm, size, pVE, nMeet, ierr))
147: PetscCallA(DMDestroy(dm, ierr))
148: PetscCallA(PetscFinalize(ierr))
149: end
150: !
151: !/*TEST
152: !
153: ! test:
154: ! suffix: 0
155: !
156: !TEST*/