imeshp_test.F90

Go to the documentation of this file.

! push parallel mesh into moab, F90 version, test
!
! This program shows how to push a mesh into MOAB in parallel from Fortran90, with sufficient
! information to resolve boundary sharing and exchange a layer of ghost information.
!
! After resolving the sharing, mesh is saved to a file, then read back again, in parallel
! the mesh is just a set of 6 quads that form a cube; it is distributed to at most 6 processors
! by default, this test is run on 2 processors
#define ERROR(rval) if (0 .ne. rval) call exit(1)
 
program imeshp_test
 
 use iso_c_binding
 implicit none
 
#include "moab/MOABConfig.h"
#include "mpif.h"
#include "iMeshP_f.h"
 
 ! declarations
 ! imesh is the instance handle
 imesh_instance imesh
 ! second instance, to check loading of mesh saved from first instance
 imesh_instance imesh2
 
 ! NUMV, NUME, NVPERE are the hardwired here; these are for the whole mesh,
 ! local mesh determined later
 integer numv, nume, nvpere
 parameter(numv = 8) ! # vertices in whole mesh
 parameter(nume = 6) ! # elements in whole mesh
 parameter(nvpere = 4) ! # vertices per element
 ! ents, verts will be arrays storing vertex/entity handles
 ibase_entityhandle, pointer :: ents(:), verts(:)
 ibase_entitysethandle root_set, root_set2, new_set
 TYPE(c_ptr) :: vertsptr, entsptr
 ! storage for vertex positions, element connectivity indices, global vertex ids
 real*8 coords(0:3*numv-1)
 integer iconn(0:4*nume-1), gids(0:numv-1)
 !
 ! local variables
 integer lgids(0:numv-1), lconn(0:4*nume-1)
 real*8 lcoords(0:3*numv-1)
 integer lnumv, lvids(0:numv-1), gvids(0:numv-1)
 integer lvpe, ltp ! lvpe = # vertices per entity, ltp = element type
 integer ic, ie, iv, istart, iend, ierr, indv, lnume, rank, sz
 integer iv2, ie2 ! after loading
 
 ! local variables for parallel runs; index2 for second imesh instance, used for loading from file
 imeshp_partitionhandle imeshp, imeshp2
 imeshp_parthandle part, part2
 ibase_handle_t mpi_comm_c
! integer MPI_COMM_WORLD
 
 
 ! vertex positions, cube; side 2
 ! (first index varying fastest)
 data coords / &
 -1., -1., -1, 1., -1., -1., 1., 1., -1., -1., 1., -1., &
 -1., -1., 1, 1., -1., 1., 1., 1., 1., -1., 1., 1. /
 
 ! quad index numbering, each quad ccw, sides then bottom then top
 data iconn / &
 0, 1, 5, 4, &
 1, 2, 6, 5, &
 2, 3, 7, 6, &
 3, 0, 4, 7, &
 0, 3, 2, 1, &
 4, 5, 6, 7 /
 
 data lvpe /4/ ! quads in this example
 data ltp / imesh_quadrilateral / ! from iBase_f.h
 
 ! initialize global vertex ids
 do iv = 0, numv-1
 lgids(iv) = iv+1
 end do
 
 ! init the parallel partition
 call mpi_init(ierr)
 call mpi_comm_size(mpi_comm_world, sz, ierr)
 call mpi_comm_rank(mpi_comm_world, rank, ierr)
 error(ierr)
 ! compute starting/ending element numbers
 lnume = nume / sz
 istart = rank * lnume
 iend = istart + lnume - 1
 if (rank .eq. sz-1) then
 iend = nume-1
 lnume = iend - istart + 1
 endif
 
 ! for my elements, figure out which vertices I use and accumulate local indices and coords
 ! lvids stores the local 0-based index for each vertex; -1 means vertex i isn't used locally
 ! also build up connectivity indices for local elements, in lconn
 do iv = 0, numv-1
 lvids(iv) = -1
 end do
 lnumv = -1
 do ie = istart, iend
 do iv = 0, lvpe-1
 indv = iconn(lvpe*ie + iv)
 if (lvids(indv) .eq. -1) then
 lnumv = lnumv + 1 ! increment local # verts
 do ic = 0, 2 ! cache local coords
 lcoords(3*lnumv+ic) = coords(3*indv+ic)
 end do
 lvids(indv) = lnumv
 gvids(lnumv) = 1+indv
 end if
 lconn(lvpe*(ie-istart)+iv) = lvids(indv)
 end do ! do iv
 end do ! do ie
 
 lnumv = lnumv + 1
 
 ! now create the mesh; this also initializes parallel sharing and ghost exchange
 imesh = 0
 imeshp = 0
 call create_mesh(imesh, imeshp, part, mpi_comm_world, lnumv, lnume, gvids, lvpe, ltp, lcoords, lconn, &
 vertsptr, entsptr, ierr)
 error(ierr)
 call c_f_pointer(vertsptr, verts, [lnumv])
 call c_f_pointer(entsptr, ents, [lnume])
 
 ! this will save the mesh in parallel
 call imeshp_saveall(%VAL(imesh), %VAL(imeshp), %VAL(part), "test2.h5m", " moab:PARALLEL=WRITE_PART ", ierr)
 error(ierr)
 call imesh_getrootset(%VAL(imesh), root_set, ierr)
 error(ierr)
 call imeshp_getnumoftypeall(%VAL(imesh), %VAL(imeshp), %VAL(root_set), %VAL(ibase_vertex), iv, ierr)
 error(ierr)
 call imeshp_getnumoftypeall(%VAL(imesh), %VAL(imeshp), %VAL(root_set), %VAL(ibase_face), ie, ierr)
 error(ierr)
 
 ! load the saved mesh in a new instance, in parallel
 call imesh_newmesh("MOAB", imesh2, ierr)
 error(ierr)
 call imesh_getrootset(%VAL(imesh2), root_set2, ierr)
 error(ierr)
 call imeshp_getcommunicator(%VAL(imesh2), mpi_comm_world, mpi_comm_c, ierr)
 error(ierr)
 call imeshp_createpartitionall(%VAL(imesh2), %VAL(mpi_comm_c), imeshp2, ierr)
 error(ierr)
 call imeshp_createpart(%VAL(imesh2), %VAL(imeshp2), part2, ierr)
 error(ierr)
 iv2 = 0
 ie2 = 0
 ! see if load works in a new file set
 call imesh_createentset(%VAL(imesh2), %VAL(0), new_set, ierr)
 error(ierr)
 call imeshp_loadall( %VAL(imesh2), %VAL(imeshp2), %VAL(new_set), "test2.h5m", &
 " moab:PARALLEL=READ_PART moab:PARALLEL_RESOLVE_SHARED_ENTS moab:PARTITION=PARALLEL_PARTITION ", &
 ierr)
 error(ierr)
 
 call imeshp_getnumoftypeall(%VAL(imesh2), %VAL(imeshp2), %VAL(root_set2), %VAL(ibase_vertex), iv2, ierr)
 error(ierr)
 call imeshp_getnumoftypeall(%VAL(imesh2), %VAL(imeshp2), %VAL(root_set2), %VAL(ibase_face), ie2, ierr)
 error(ierr)
 
 if ( (iv.ne.iv2 ) .or. (ie.ne.ie2) ) then
 write(0, *) "inconsistent number of elements and vertices"
 write(0, *) " saved: " , iv, ie, " loaded: " , iv2, ie2
 call exit(1)
 endif
 
 call imesh_dtor(%VAL(imesh), ierr);
 error(ierr);
 call imesh_dtor(%VAL(imesh2), ierr);
 error(ierr);
 
 call mpi_finalize(ierr)
 stop
end program imeshp_test
 
subroutine create_mesh( &
 ! interfaces
 imesh, imeshp, part, &
 ! input
 comm, numv, nume, vgids, nvpe, tp, posn, iconn, &
 ! output
 vertsPtr, entsPtr, ierr)
 !
 ! create a mesh with numv vertices and nume elements, with elements of type tp
 ! vertices have positions in posn (3 coordinates each, interleaved xyzxyz...), indexed from 0
 ! elements have nvpe vertices per entity, with connectivity indices stored in iconn, referencing
 ! vertices using 0-based indices; vertex and entity handles are output in arrays passed in
 !
 ! if imesh/imeshp are 0, imesh/imeshp are initialized in this subroutine
 !
 
 use iso_c_binding
 implicit none
 
#ifdef MOAB_HAVE_MPI
# include "iMeshP_f.h"
# include "mpif.h"
#else
# include "iMesh_f.h"
#endif
 
 ! subroutine arguments
 imesh_instance imesh
 TYPE(c_ptr) :: vertsPtr, entsPtr
 integer numv, nume, nvpe, vgids(0:*), iconn(0:*), ierr, tp
 real*8 posn(0:*)
 
 imeshp_partitionhandle imeshp
 integer comm
 
 
 ! local variables
 integer comm_sz, comm_rank, numa, numo, iv, ie
 TYPE(c_ptr) :: statsPtr
 integer, allocatable, target :: stats(:)
 ibase_taghandle tagh
 integer i
 ibase_entityhandle, pointer :: verts(:), ents(:)
 ibase_entityhandle, allocatable :: conn(:)
 ibase_entitysethandle root_set
 ibase_entitysethandle file_set
 
 ibase_handle_t mpi_comm_c
 TYPE(c_ptr) :: partsPtr
 imeshp_parthandle, pointer :: parts(:)
 imeshp_parthandle part
 integer partsa, partso
 character (len=10) filename
 
 ! create the Mesh instance
 if (imesh .eq. 0) then
 call imesh_newmesh("MOAB", imesh, ierr)
 end if
 
 
 if (imeshp .eq. 0) then
 call imeshp_getcommunicator(%VAL(imesh), mpi_comm_world, mpi_comm_c, ierr)
 error(ierr)
 call imeshp_createpartitionall(%VAL(imesh), %VAL(mpi_comm_c), imeshp, ierr)
 error(ierr)
 call imeshp_createpart(%VAL(imesh), %VAL(imeshp), part, ierr)
 error(ierr)
 else
 partsa = 0
 call imeshp_getlocalparts(%VAL(imesh), %VAL(imeshp), partsptr, partsa, partso, ierr)
 error(ierr)
 call c_f_pointer(partsptr, parts, [partso])
 part = parts(1)
 end if
 call mpi_comm_rank(comm, comm_rank, ierr)
 error(ierr)
 call mpi_comm_size(comm, comm_sz, ierr)
 error(ierr)
 
 ! create the vertices, all in one call
 numa = 0
 call imesh_createvtxarr(%VAL(imesh), %VAL(numv), %VAL(ibase_interleaved), posn, %VAL(3*numv), &
 vertsptr, numa, numo, ierr)
 error(ierr)
 
 ! fill in the connectivity array, based on indexing from iconn
 allocate (conn(0:nvpe*nume-1))
 call c_f_pointer(vertsptr, verts, [numv])
 do i = 0, nvpe*nume-1
 conn(i) = verts(1+iconn(i))
 end do
 ! create the elements
 numa = 0
 allocate(stats(0:nume-1))
 statsptr = c_loc(stats(0))
 call imesh_createentarr(%VAL(imesh), %VAL(tp), conn, %VAL(nvpe*nume), &
 entsptr, numa, numo, statsptr, numa, numo, ierr)
 deallocate(stats)
 deallocate(conn)
 
 ! take care of parallel stuff
 
 ! add entities to part, using iMesh
 call c_f_pointer(entsptr, ents, [numo])
 call imesh_addentarrtoset(%VAL(imesh), ents, %VAL(numo), %VAL(part), ierr)
 error(ierr)
 ! set global ids on vertices, needed for sharing between procs
 call imesh_gettaghandle(%VAL(imesh), "GLOBAL_ID", tagh, ierr, %VAL(9))
 if (ibase_success .ne. ierr) then
 ! didn't get handle, need to create the tag
 call imesh_createtag(%VAL(imesh), "GLOBAL_ID", %VAL(ibase_integer), tagh, ierr)
 error(ierr)
 end if
 call imesh_setintarrdata(%VAL(imesh), verts, %VAL(numv), %VAL(tagh), vgids, %VAL(numv), ierr)
 error(ierr)
 ! now resolve shared verts and exchange ghost cells
 call imeshp_syncmeshall(%VAL(imesh), %VAL(imeshp), ierr)
 error(ierr)
 
 call imeshp_createghostentsall(%VAL(imesh), %VAL(imeshp), %VAL(2), %VAL(1), %VAL(1), %VAL(0), ierr)
 error(ierr)
 
 call imesh_freememory(%VAL(imesh), vertsptr);
 call imesh_freememory(%VAL(imesh), entsptr);
 
 return
end subroutine create_mesh