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