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