Actual source code: ex2f90.F
petsc-3.5.4 2015-05-23
1: program main
2: implicit none
3: !
4: #include <finclude/petsc.h90>
5: !
6: DM dm
7: PetscInt, target, dimension(3) :: EC
8: PetscInt, target, dimension(2) :: VE
9: PetscInt, pointer :: pEC(:), pVE(:)
10: PetscInt, pointer :: nClosure(:)
11: PetscInt, pointer :: nJoin(:)
12: PetscInt, pointer :: nMeet(:)
13: PetscInt dim, cell, size
14: PetscInt i0,i1,i2,i3,i4,i5,i6,i7
15: PetscInt i8,i9,i10,i11
16: PetscErrorCode ierr
18: i0 = 0
19: i1 = 1
20: i2 = 2
21: i3 = 3
22: i4 = 4
23: i5 = 5
24: i6 = 6
25: i7 = 7
26: i8 = 8
27: i9 = 9
28: i10 = 10
29: i11 = 11
31: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
32: CHKERRQ(ierr)
33: call DMPlexCreate(PETSC_COMM_WORLD, dm, ierr)
34: CHKERRQ(ierr)
35: dim = 2
36: call DMPlexSetDimension(dm, dim, ierr)
37: CHKERRQ(ierr)
39: ! Make Doublet Mesh from Fig 2 of Flexible Representation of Computational Meshes,
40: ! except indexing is from 0 instead of 1 and we obey the new restrictions on
41: ! numbering: cells, vertices, faces, edges
42: call DMPlexSetChart(dm, i0, i11, ierr)
43: CHKERRQ(ierr)
44: ! cells
45: call DMPlexSetConeSize(dm, i0, i3, ierr)
46: CHKERRQ(ierr)
47: call DMPlexSetConeSize(dm, i1, i3, ierr)
48: CHKERRQ(ierr)
49: ! edges
50: call DMPlexSetConeSize(dm, i6, i2, ierr)
51: CHKERRQ(ierr)
52: call DMPlexSetConeSize(dm, i7, i2, ierr)
53: CHKERRQ(ierr)
54: call DMPlexSetConeSize(dm, i8, i2, ierr)
55: CHKERRQ(ierr)
56: call DMPlexSetConeSize(dm, i9, i2, ierr)
57: CHKERRQ(ierr)
58: call DMPlexSetConeSize(dm, i10, i2, ierr)
59: CHKERRQ(ierr)
61: call DMSetUp(dm, ierr)
62: CHKERRQ(ierr)
64: EC(1) = 6
65: EC(2) = 7
66: EC(3) = 8
67: pEC => EC
68: call DMPlexSetCone(dm, i0, pEC, ierr)
69: CHKERRQ(ierr)
70: EC(1) = 7
71: EC(2) = 9
72: EC(3) = 10
73: pEC => EC
74: call DMPlexSetCone(dm, i1 , pEC, ierr)
75: CHKERRQ(ierr)
77: VE(1) = 2
78: VE(2) = 3
79: pVE => VE
80: call DMPlexSetCone(dm, i6 , pVE, ierr)
81: CHKERRQ(ierr)
82: VE(1) = 3
83: VE(2) = 4
84: pVE => VE
85: call DMPlexSetCone(dm, i7 , pVE, ierr)
86: CHKERRQ(ierr)
87: VE(1) = 4
88: VE(2) = 2
89: pVE => VE
90: call DMPlexSetCone(dm, i8 , pVE, ierr)
91: CHKERRQ(ierr)
92: VE(1) = 3
93: VE(2) = 5
94: pVE => VE
95: call DMPlexSetCone(dm, i9 , pVE, ierr)
96: CHKERRQ(ierr)
97: VE(1) = 5
98: VE(2) = 4
99: pVE => VE
100: call DMPlexSetCone(dm, i10 , pVE, ierr)
101: CHKERRQ(ierr)
103: call DMPlexSymmetrize(dm,ierr)
104: CHKERRQ(ierr)
105: call DMPlexStratify(dm,ierr)
106: CHKERRQ(ierr)
107: call DMView(dm, PETSC_VIEWER_STDOUT_WORLD, ierr)
109: ! Test Closure
110: do cell = 0,1
111: call DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &
112: & nClosure, ierr)
113: CHKERRQ(ierr)
114: write(*,*) nClosure
115: call DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &
116: & nClosure, ierr)
117: end do
119: ! Test Join
120: size = 2
121: VE(1) = 6
122: VE(2) = 7
123: pVE => VE
124: call DMPlexGetJoin(dm, size, pVE, nJoin, ierr)
125: write(*,*) 'Join of',pVE,'is',nJoin
126: call DMPlexRestoreJoin(dm, size, pVE, nJoin, ierr)
127: size = 2
128: VE(1) = 9
129: VE(2) = 7
130: pVE => VE
131: call DMPlexGetJoin(dm, size, pVE, nJoin, ierr)
132: write(*,*) 'Join of',pVE,'is',nJoin
133: call DMPlexRestoreJoin(dm, size, pVE, nJoin, ierr)
134: ! Test Full Join
135: size = 3
136: EC(1) = 3
137: EC(2) = 4
138: EC(3) = 5
139: pEC => EC
140: call DMPlexGetFullJoin(dm, size, pEC, nJoin, ierr)
141: write(*,*) 'Full Join of',pEC,'is',nJoin
142: call DMPlexRestoreJoin(dm, size, pEC, nJoin, ierr)
143: ! Test Meet
144: size = 2
145: VE(1) = 0
146: VE(2) = 1
147: pVE => VE
148: call DMPlexGetMeet(dm, size, pVE, nMeet, ierr)
149: write(*,*) 'Meet of',pVE,'is',nMeet
150: call DMPlexRestoreMeet(dm, size, pVE, nMeet, ierr)
151: size = 2
152: VE(1) = 6
153: VE(2) = 7
154: pVE => VE
155: call DMPlexGetMeet(dm, size, pVE, nMeet, ierr)
156: write(*,*) 'Meet of',pVE,'is',nMeet
157: call DMPlexRestoreMeet(dm, size, pVE, nMeet, ierr)
159: call DMDestroy(dm, ierr)
160: CHKERRQ(ierr)
161: call PetscFinalize(ierr)
162: CHKERRQ(ierr)
163: end