1!  push parallel mesh into moab, F90 version, test
2!
3! This program shows how to push a mesh into MOAB in parallel from Fortran90, with sufficient
4! information to resolve boundary sharing and exchange a layer of ghost information.
5!
6! After resolving the sharing, mesh is saved to a file, then read back again, in parallel
7!  the mesh is just a set of 6 quads that form a cube; it is distributed to at most 6 processors
8!  by default, this test is run on 2 processors
9#define ERROR(rval) if (0 .ne. rval) call exit(1)
10
11program imeshp_test
12
13  use ISO_C_BINDING
14  implicit none
15
16#include "moab/MOABConfig.h"
17#include "mpif.h"
18#include "iMeshP_f.h"
19
20  ! declarations
21  ! imesh is the instance handle
22  iMesh_Instance imesh
23  ! second instance, to check loading of mesh saved from first instance
24  iMesh_Instance imesh2
25
26  ! NUMV, NUME, NVPERE are the hardwired here; these are for the whole mesh,
27  ! local mesh determined later
28  integer NUMV, NUME, NVPERE
29  parameter (NUMV = 8)   ! # vertices in whole mesh
30  parameter (NUME = 6)   ! # elements in whole mesh
31  parameter (NVPERE = 4) ! # vertices per element
32  ! ents, verts will be arrays storing vertex/entity handles
33  iBase_EntityHandle, pointer :: ents(:), verts(:)
34  iBase_EntitySetHandle root_set, root_set2, new_set
35  TYPE(C_PTR) :: vertsPtr, entsPtr
36  ! storage for vertex positions, element connectivity indices, global vertex ids
37  real*8 coords(0:3*NUMV-1)
38  integer iconn(0:4*NUME-1), gids(0:NUMV-1)
39  !
40  ! local variables
41  integer lgids(0:NUMV-1), lconn(0:4*NUME-1)
42  real*8 lcoords(0:3*NUMV-1)
43  integer lnumv, lvids(0:NUMV-1), gvids(0:NUMV-1)
44  integer lvpe, ltp ! lvpe = # vertices per entity, ltp = element type
45  integer ic, ie, iv, istart, iend, ierr, indv, lnume, rank, sz
46  integer iv2, ie2 !  after loading
47
48  ! local variables for parallel runs; index2 for second imesh instance, used for loading from file
49  iMeshP_PartitionHandle imeshp, imeshp2
50  iMeshP_PartHandle part, part2
51  IBASE_HANDLE_T mpi_comm_c
52!    integer MPI_COMM_WORLD
53
54
55  ! vertex positions, cube; side 2
56  ! (first index varying fastest)
57  data coords / &
58       -1., -1., -1,  1., -1., -1.,  1., 1., -1.,  -1., 1., -1., &
59       -1., -1.,  1,  1., -1.,  1.,  1., 1.,  1.,  -1., 1.,  1. /
60
61  ! quad index numbering, each quad ccw, sides then bottom then top
62  data iconn / &
63       0, 1, 5, 4,  &
64       1, 2, 6, 5,  &
65       2, 3, 7, 6,  &
66       3, 0, 4, 7,  &
67       0, 3, 2, 1,  &
68       4, 5, 6, 7 /
69
70  data lvpe /4/ ! quads in this example
71  data ltp / iMesh_QUADRILATERAL / ! from iBase_f.h
72
73  ! initialize global vertex ids
74  do iv = 0, NUMV-1
75     lgids(iv) = iv+1
76  end do
77
78  ! init the parallel partition
79  call MPI_INIT(ierr)
80  call MPI_COMM_SIZE(MPI_COMM_WORLD, sz, ierr)
81  call MPI_COMM_RANK(MPI_COMM_WORLD, rank, ierr)
82  ERROR(ierr)
83  ! compute starting/ending element numbers
84  lnume = NUME / sz
85  istart = rank * lnume
86  iend = istart + lnume - 1
87  if (rank .eq. sz-1) then
88     iend = NUME-1
89     lnume = iend - istart + 1
90  endif
91
92  ! for my elements, figure out which vertices I use and accumulate local indices and coords
93  ! lvids stores the local 0-based index for each vertex; -1 means vertex i isn't used locally
94  ! also build up connectivity indices for local elements, in lconn
95  do iv = 0, NUMV-1
96     lvids(iv) = -1
97  end do
98  lnumv = -1
99  do ie = istart, iend
100     do iv = 0, lvpe-1
101        indv = iconn(lvpe*ie + iv)
102        if (lvids(indv) .eq. -1) then
103           lnumv = lnumv + 1 ! increment local # verts
104           do ic = 0, 2 ! cache local coords
105              lcoords(3*lnumv+ic) = coords(3*indv+ic)
106           end do
107           lvids(indv) = lnumv
108           gvids(lnumv) = 1+indv
109        end if
110        lconn(lvpe*(ie-istart)+iv) = lvids(indv)
111     end do  ! do iv
112  end do  ! do ie
113
114  lnumv = lnumv + 1
115
116  ! now create the mesh; this also initializes parallel sharing and ghost exchange
117  imesh = 0
118  imeshp = 0
119  call create_mesh(imesh, imeshp, part, MPI_COMM_WORLD, lnumv, lnume, gvids, lvpe, ltp, lcoords, lconn, &
120       vertsPtr, entsPtr, ierr)
121  ERROR(ierr)
122  call c_f_pointer(vertsPtr, verts, [lnumv])
123  call c_f_pointer(entsPtr, ents, [lnume])
124
125 ! this will save the mesh in parallel
126  call iMeshP_saveAll(%VAL(imesh), %VAL(imeshp), %VAL(part), "test2.h5m", " moab:PARALLEL=WRITE_PART ", ierr)
127  ERROR(ierr)
128  call iMesh_getRootSet(%VAL(imesh), root_set, ierr)
129  ERROR(ierr)
130  call iMeshP_getNumOfTypeAll(%VAL(imesh), %VAL(imeshp), %VAL(root_set), %VAL(iBase_VERTEX), iv, ierr)
131  ERROR(ierr)
132  call iMeshP_getNumOfTypeAll(%VAL(imesh), %VAL(imeshp), %VAL(root_set), %VAL(iBase_FACE), ie, ierr)
133  ERROR(ierr)
134
135  ! load the saved mesh in a new instance, in parallel
136  call iMesh_newMesh("MOAB", imesh2, ierr)
137  ERROR(ierr)
138  call iMesh_getRootSet(%VAL(imesh2), root_set2, ierr)
139  ERROR(ierr)
140  call iMeshP_getCommunicator(%VAL(imesh2), MPI_COMM_WORLD, mpi_comm_c, ierr)
141  ERROR(ierr)
142  call iMeshP_createPartitionAll(%VAL(imesh2), %VAL(mpi_comm_c), imeshp2, ierr)
143  ERROR(ierr)
144  call iMeshP_createPart(%VAL(imesh2), %VAL(imeshp2), part2, ierr)
145  ERROR(ierr)
146  iv2 = 0
147  ie2 = 0
148 ! see if load works in a new file set
149  call iMesh_createEntSet(%VAL(imesh2), %VAL(0), new_set, ierr)
150  ERROR(ierr)
151  call iMeshP_loadAll( %VAL(imesh2), %VAL(imeshp2), %VAL(new_set), "test2.h5m", &
152   " moab:PARALLEL=READ_PART moab:PARALLEL_RESOLVE_SHARED_ENTS moab:PARTITION=PARALLEL_PARTITION ", &
153   ierr)
154  ERROR(ierr)
155
156  call iMeshP_getNumOfTypeAll(%VAL(imesh2), %VAL(imeshp2), %VAL(root_set2), %VAL(iBase_VERTEX), iv2, ierr)
157  ERROR(ierr)
158  call iMeshP_getNumOfTypeAll(%VAL(imesh2), %VAL(imeshp2), %VAL(root_set2), %VAL(iBase_FACE), ie2, ierr)
159  ERROR(ierr)
160
161  if ( (iv.ne.iv2 ) .or. (ie.ne.ie2) ) then
162     write(0, *) "inconsistent number of elements and vertices"
163     write(0, *) " saved: " , iv, ie, " loaded: " , iv2, ie2
164     call exit(1)
165  endif
166
167  call iMesh_dtor(%VAL(imesh), ierr);
168  ERROR(ierr);
169  call iMesh_dtor(%VAL(imesh2), ierr);
170  ERROR(ierr);
171
172  call MPI_FINALIZE(ierr)
173  stop
174end program imeshp_test
175
176subroutine create_mesh( &
177  !     interfaces
178     imesh, imeshp, part, &
179  !     input
180     comm, numv, nume, vgids, nvpe, tp, posn, iconn, &
181  !     output
182     vertsPtr, entsPtr, ierr)
183  !
184  ! create a mesh with numv vertices and nume elements, with elements of type tp
185  ! vertices have positions in posn (3 coordinates each, interleaved xyzxyz...), indexed from 0
186  ! elements have nvpe vertices per entity, with connectivity indices stored in iconn, referencing
187  ! vertices using 0-based indices; vertex and entity handles are output in arrays passed in
188  !
189  ! if imesh/imeshp are 0, imesh/imeshp are initialized in this subroutine
190  !
191
192  use ISO_C_BINDING
193  implicit none
194
195#ifdef MOAB_HAVE_MPI
196#  include "iMeshP_f.h"
197#  include "mpif.h"
198#else
199#  include "iMesh_f.h"
200#endif
201
202  ! subroutine arguments
203  iMesh_Instance imesh
204  TYPE(C_PTR) :: vertsPtr, entsPtr
205  integer numv, nume, nvpe, vgids(0:*), iconn(0:*), ierr, tp
206  real*8 posn(0:*)
207
208  iMeshP_PartitionHandle imeshp
209  integer comm
210
211
212  ! local variables
213  integer comm_sz, comm_rank, numa, numo, iv, ie
214  TYPE(C_PTR) :: statsPtr
215  integer, allocatable, target :: stats(:)
216  iBase_TagHandle tagh
217  integer i
218  iBase_EntityHandle, pointer :: verts(:), ents(:)
219  iBase_EntityHandle, allocatable :: conn(:)
220  iBase_EntitySetHandle root_set
221  iBase_EntitySetHandle file_set
222
223  IBASE_HANDLE_T mpi_comm_c
224  TYPE(C_PTR) :: partsPtr
225  iMeshP_PartHandle, pointer :: parts(:)
226  iMeshP_PartHandle part
227  integer partsa, partso
228  character (len=10) filename
229
230  ! create the Mesh instance
231  if (imesh .eq. 0) then
232     call iMesh_newMesh("MOAB", imesh, ierr)
233  end if
234
235
236  if (imeshp .eq. 0) then
237     call iMeshP_getCommunicator(%VAL(imesh), MPI_COMM_WORLD, mpi_comm_c, ierr)
238     ERROR(ierr)
239     call iMeshP_createPartitionAll(%VAL(imesh), %VAL(mpi_comm_c), imeshp, ierr)
240     ERROR(ierr)
241     call iMeshP_createPart(%VAL(imesh), %VAL(imeshp), part, ierr)
242     ERROR(ierr)
243  else
244     partsa = 0
245     call iMeshP_getLocalParts(%VAL(imesh), %VAL(imeshp), partsPtr, partsa, partso, ierr)
246     ERROR(ierr)
247     call c_f_pointer(partsPtr, parts, [partso])
248     part = parts(1)
249  end if
250  call MPI_COMM_RANK(comm, comm_rank, ierr)
251  ERROR(ierr)
252  call MPI_COMM_SIZE(comm, comm_sz, ierr)
253  ERROR(ierr)
254
255  ! create the vertices, all in one call
256  numa = 0
257  call iMesh_createVtxArr(%VAL(imesh), %VAL(numv), %VAL(iBase_INTERLEAVED), posn, %VAL(3*numv), &
258       vertsPtr, numa, numo, ierr)
259  ERROR(ierr)
260
261  ! fill in the connectivity array, based on indexing from iconn
262  allocate (conn(0:nvpe*nume-1))
263  call c_f_pointer(vertsPtr, verts, [numv])
264  do i = 0, nvpe*nume-1
265     conn(i) = verts(1+iconn(i))
266  end do
267  ! create the elements
268  numa = 0
269  allocate(stats(0:nume-1))
270  statsPtr = C_LOC(stats(0))
271  call iMesh_createEntArr(%VAL(imesh), %VAL(tp), conn, %VAL(nvpe*nume), &
272       entsPtr, numa, numo, statsPtr, numa, numo, ierr)
273  deallocate(stats)
274  deallocate(conn)
275
276  ! take care of parallel stuff
277
278  ! add entities to part, using iMesh
279  call c_f_pointer(entsPtr, ents, [numo])
280  call iMesh_addEntArrToSet(%VAL(imesh), ents, %VAL(numo), %VAL(part), ierr)
281  ERROR(ierr)
282  ! set global ids on vertices, needed for sharing between procs
283  call iMesh_getTagHandle(%VAL(imesh), "GLOBAL_ID", tagh, ierr, %VAL(9))
284  if (iBase_SUCCESS .ne. ierr) then
285     ! didn't get handle, need to create the tag
286     call iMesh_createTag(%VAL(imesh), "GLOBAL_ID", %VAL(iBase_INTEGER), tagh, ierr)
287     ERROR(ierr)
288  end if
289  call iMesh_setIntArrData(%VAL(imesh), verts, %VAL(numv), %VAL(tagh), vgids, %VAL(numv), ierr)
290  ERROR(ierr)
291  ! now resolve shared verts and exchange ghost cells
292  call iMeshP_syncMeshAll(%VAL(imesh), %VAL(imeshp), ierr)
293  ERROR(ierr)
294
295  call iMeshP_createGhostEntsAll(%VAL(imesh), %VAL(imeshp), %VAL(2), %VAL(1), %VAL(1), %VAL(0), ierr)
296  ERROR(ierr)
297
298  call iMesh_freeMemory(%VAL(imesh), vertsPtr);
299  call iMesh_freeMemory(%VAL(imesh), entsPtr);
300
301  return
302end subroutine create_mesh
303