MOAB: Mesh Oriented datABase  (version 5.2.1)
imeshp_test.F90
Go to the documentation of this file.
00001 !  push parallel mesh into moab, F90 version, test
00002 !
00003 ! This program shows how to push a mesh into MOAB in parallel from Fortran90, with sufficient
00004 ! information to resolve boundary sharing and exchange a layer of ghost information.
00005 !
00006 ! After resolving the sharing, mesh is saved to a file, then read back again, in parallel
00007 !  the mesh is just a set of 6 quads that form a cube; it is distributed to at most 6 processors
00008 !  by default, this test is run on 2 processors
00009 #define ERROR(rval) if (0 .ne. rval) call exit(1)
00010 
00011 program imeshp_test
00012 
00013   use ISO_C_BINDING
00014   implicit none
00015 
00016 #include "moab/MOABConfig.h"
00017 #include "mpif.h"
00018 #include "iMeshP_f.h"
00019 
00020   ! declarations
00021   ! imesh is the instance handle
00022   iMesh_Instance imesh
00023   ! second instance, to check loading of mesh saved from first instance
00024   iMesh_Instance imesh2
00025 
00026   ! NUMV, NUME, NVPERE are the hardwired here; these are for the whole mesh,
00027   ! local mesh determined later
00028   integer NUMV, NUME, NVPERE
00029   parameter (NUMV = 8)   ! # vertices in whole mesh
00030   parameter (NUME = 6)   ! # elements in whole mesh
00031   parameter (NVPERE = 4) ! # vertices per element
00032   ! ents, verts will be arrays storing vertex/entity handles
00033   iBase_EntityHandle, pointer :: ents(:), verts(:)
00034   iBase_EntitySetHandle root_set, root_set2, new_set
00035   TYPE(C_PTR) :: vertsPtr, entsPtr
00036   ! storage for vertex positions, element connectivity indices, global vertex ids
00037   real*8 coords(0:3*NUMV-1)
00038   integer iconn(0:4*NUME-1), gids(0:NUMV-1)
00039   !
00040   ! local variables
00041   integer lgids(0:NUMV-1), lconn(0:4*NUME-1)
00042   real*8 lcoords(0:3*NUMV-1)
00043   integer lnumv, lvids(0:NUMV-1), gvids(0:NUMV-1)
00044   integer lvpe, ltp ! lvpe = # vertices per entity, ltp = element type
00045   integer ic, ie, iv, istart, iend, ierr, indv, lnume, rank, sz
00046   integer iv2, ie2 !  after loading
00047 
00048   ! local variables for parallel runs; index2 for second imesh instance, used for loading from file
00049   iMeshP_PartitionHandle imeshp, imeshp2
00050   iMeshP_PartHandle part, part2
00051   IBASE_HANDLE_T mpi_comm_c
00052 !    integer MPI_COMM_WORLD
00053 
00054 
00055   ! vertex positions, cube; side 2
00056   ! (first index varying fastest)
00057   data coords / &
00058        -1., -1., -1,  1., -1., -1.,  1., 1., -1.,  -1., 1., -1., &
00059        -1., -1.,  1,  1., -1.,  1.,  1., 1.,  1.,  -1., 1.,  1. /
00060 
00061   ! quad index numbering, each quad ccw, sides then bottom then top
00062   data iconn / &
00063        0, 1, 5, 4,  &
00064        1, 2, 6, 5,  &
00065        2, 3, 7, 6,  &
00066        3, 0, 4, 7,  &
00067        0, 3, 2, 1,  &
00068        4, 5, 6, 7 /
00069 
00070   data lvpe /4/ ! quads in this example
00071   data ltp / iMesh_QUADRILATERAL / ! from iBase_f.h
00072 
00073   ! initialize global vertex ids
00074   do iv = 0, NUMV-1
00075      lgids(iv) = iv+1
00076   end do
00077 
00078   ! init the parallel partition
00079   call MPI_INIT(ierr)
00080   call MPI_COMM_SIZE(MPI_COMM_WORLD, sz, ierr)
00081   call MPI_COMM_RANK(MPI_COMM_WORLD, rank, ierr)
00082   ERROR(ierr)
00083   ! compute starting/ending element numbers
00084   lnume = NUME / sz
00085   istart = rank * lnume
00086   iend = istart + lnume - 1
00087   if (rank .eq. sz-1) then
00088      iend = NUME-1
00089      lnume = iend - istart + 1
00090   endif
00091 
00092   ! for my elements, figure out which vertices I use and accumulate local indices and coords
00093   ! lvids stores the local 0-based index for each vertex; -1 means vertex i isn't used locally
00094   ! also build up connectivity indices for local elements, in lconn
00095   do iv = 0, NUMV-1
00096      lvids(iv) = -1
00097   end do
00098   lnumv = -1
00099   do ie = istart, iend
00100      do iv = 0, lvpe-1
00101         indv = iconn(lvpe*ie + iv)
00102         if (lvids(indv) .eq. -1) then
00103            lnumv = lnumv + 1 ! increment local # verts
00104            do ic = 0, 2 ! cache local coords
00105               lcoords(3*lnumv+ic) = coords(3*indv+ic)
00106            end do
00107            lvids(indv) = lnumv
00108            gvids(lnumv) = 1+indv
00109         end if
00110         lconn(lvpe*(ie-istart)+iv) = lvids(indv)
00111      end do  ! do iv
00112   end do  ! do ie
00113 
00114   lnumv = lnumv + 1
00115 
00116   ! now create the mesh; this also initializes parallel sharing and ghost exchange
00117   imesh = 0
00118   imeshp = 0
00119   call create_mesh(imesh, imeshp, part, MPI_COMM_WORLD, lnumv, lnume, gvids, lvpe, ltp, lcoords, lconn, &
00120        vertsPtr, entsPtr, ierr)
00121   ERROR(ierr)
00122   call c_f_pointer(vertsPtr, verts, [lnumv])
00123   call c_f_pointer(entsPtr, ents, [lnume])
00124 
00125  ! this will save the mesh in parallel
00126   call iMeshP_saveAll(%VAL(imesh), %VAL(imeshp), %VAL(part), "test2.h5m", " moab:PARALLEL=WRITE_PART ", ierr)
00127   ERROR(ierr)
00128   call iMesh_getRootSet(%VAL(imesh), root_set, ierr)
00129   ERROR(ierr)
00130   call iMeshP_getNumOfTypeAll(%VAL(imesh), %VAL(imeshp), %VAL(root_set), %VAL(iBase_VERTEX), iv, ierr)
00131   ERROR(ierr)
00132   call iMeshP_getNumOfTypeAll(%VAL(imesh), %VAL(imeshp), %VAL(root_set), %VAL(iBase_FACE), ie, ierr)
00133   ERROR(ierr)
00134 
00135   ! load the saved mesh in a new instance, in parallel
00136   call iMesh_newMesh("MOAB", imesh2, ierr)
00137   ERROR(ierr)
00138   call iMesh_getRootSet(%VAL(imesh2), root_set2, ierr)
00139   ERROR(ierr)
00140   call iMeshP_getCommunicator(%VAL(imesh2), MPI_COMM_WORLD, mpi_comm_c, ierr)
00141   ERROR(ierr)
00142   call iMeshP_createPartitionAll(%VAL(imesh2), %VAL(mpi_comm_c), imeshp2, ierr)
00143   ERROR(ierr)
00144   call iMeshP_createPart(%VAL(imesh2), %VAL(imeshp2), part2, ierr)
00145   ERROR(ierr)
00146   iv2 = 0
00147   ie2 = 0
00148  ! see if load works in a new file set
00149   call iMesh_createEntSet(%VAL(imesh2), %VAL(0), new_set, ierr)
00150   ERROR(ierr)
00151   call iMeshP_loadAll( %VAL(imesh2), %VAL(imeshp2), %VAL(new_set), "test2.h5m", &
00152    " moab:PARALLEL=READ_PART moab:PARALLEL_RESOLVE_SHARED_ENTS moab:PARTITION=PARALLEL_PARTITION ", &
00153    ierr)
00154   ERROR(ierr)
00155 
00156   call iMeshP_getNumOfTypeAll(%VAL(imesh2), %VAL(imeshp2), %VAL(root_set2), %VAL(iBase_VERTEX), iv2, ierr)
00157   ERROR(ierr)
00158   call iMeshP_getNumOfTypeAll(%VAL(imesh2), %VAL(imeshp2), %VAL(root_set2), %VAL(iBase_FACE), ie2, ierr)
00159   ERROR(ierr)
00160 
00161   if ( (iv.ne.iv2 ) .or. (ie.ne.ie2) ) then
00162      write(0, *) "inconsistent number of elements and vertices"
00163      write(0, *) " saved: " , iv, ie, " loaded: " , iv2, ie2
00164      call exit(1)
00165   endif
00166 
00167   call iMesh_dtor(%VAL(imesh), ierr);
00168   ERROR(ierr);
00169   call iMesh_dtor(%VAL(imesh2), ierr);
00170   ERROR(ierr);
00171 
00172   call MPI_FINALIZE(ierr)
00173   stop
00174 end program imeshp_test
00175 
00176 subroutine create_mesh( &
00177   !     interfaces
00178      imesh, imeshp, part, &
00179   !     input
00180      comm, numv, nume, vgids, nvpe, tp, posn, iconn, &
00181   !     output
00182      vertsPtr, entsPtr, ierr)
00183   !
00184   ! create a mesh with numv vertices and nume elements, with elements of type tp
00185   ! vertices have positions in posn (3 coordinates each, interleaved xyzxyz...), indexed from 0
00186   ! elements have nvpe vertices per entity, with connectivity indices stored in iconn, referencing
00187   ! vertices using 0-based indices; vertex and entity handles are output in arrays passed in
00188   !
00189   ! if imesh/imeshp are 0, imesh/imeshp are initialized in this subroutine
00190   !
00191 
00192   use ISO_C_BINDING
00193   implicit none
00194 
00195 #ifdef MOAB_HAVE_MPI
00196 #  include "iMeshP_f.h"
00197 #  include "mpif.h"
00198 #else
00199 #  include "iMesh_f.h"
00200 #endif
00201 
00202   ! subroutine arguments
00203   iMesh_Instance imesh
00204   TYPE(C_PTR) :: vertsPtr, entsPtr
00205   integer numv, nume, nvpe, vgids(0:*), iconn(0:*), ierr, tp
00206   real*8 posn(0:*)
00207 
00208   iMeshP_PartitionHandle imeshp
00209   integer comm
00210 
00211 
00212   ! local variables
00213   integer comm_sz, comm_rank, numa, numo, iv, ie
00214   TYPE(C_PTR) :: statsPtr
00215   integer, allocatable, target :: stats(:)
00216   iBase_TagHandle tagh
00217   integer i
00218   iBase_EntityHandle, pointer :: verts(:), ents(:)
00219   iBase_EntityHandle, allocatable :: conn(:)
00220   iBase_EntitySetHandle root_set
00221   iBase_EntitySetHandle file_set
00222 
00223   IBASE_HANDLE_T mpi_comm_c
00224   TYPE(C_PTR) :: partsPtr
00225   iMeshP_PartHandle, pointer :: parts(:)
00226   iMeshP_PartHandle part
00227   integer partsa, partso
00228   character (len=10) filename
00229 
00230   ! create the Mesh instance
00231   if (imesh .eq. 0) then
00232      call iMesh_newMesh("MOAB", imesh, ierr)
00233   end if
00234 
00235 
00236   if (imeshp .eq. 0) then
00237      call iMeshP_getCommunicator(%VAL(imesh), MPI_COMM_WORLD, mpi_comm_c, ierr)
00238      ERROR(ierr)
00239      call iMeshP_createPartitionAll(%VAL(imesh), %VAL(mpi_comm_c), imeshp, ierr)
00240      ERROR(ierr)
00241      call iMeshP_createPart(%VAL(imesh), %VAL(imeshp), part, ierr)
00242      ERROR(ierr)
00243   else
00244      partsa = 0
00245      call iMeshP_getLocalParts(%VAL(imesh), %VAL(imeshp), partsPtr, partsa, partso, ierr)
00246      ERROR(ierr)
00247      call c_f_pointer(partsPtr, parts, [partso])
00248      part = parts(1)
00249   end if
00250   call MPI_COMM_RANK(comm, comm_rank, ierr)
00251   ERROR(ierr)
00252   call MPI_COMM_SIZE(comm, comm_sz, ierr)
00253   ERROR(ierr)
00254 
00255   ! create the vertices, all in one call
00256   numa = 0
00257   call iMesh_createVtxArr(%VAL(imesh), %VAL(numv), %VAL(iBase_INTERLEAVED), posn, %VAL(3*numv), &
00258        vertsPtr, numa, numo, ierr)
00259   ERROR(ierr)
00260 
00261   ! fill in the connectivity array, based on indexing from iconn
00262   allocate (conn(0:nvpe*nume-1))
00263   call c_f_pointer(vertsPtr, verts, [numv])
00264   do i = 0, nvpe*nume-1
00265      conn(i) = verts(1+iconn(i))
00266   end do
00267   ! create the elements
00268   numa = 0
00269   allocate(stats(0:nume-1))
00270   statsPtr = C_LOC(stats(0))
00271   call iMesh_createEntArr(%VAL(imesh), %VAL(tp), conn, %VAL(nvpe*nume), &
00272        entsPtr, numa, numo, statsPtr, numa, numo, ierr)
00273   deallocate(stats)
00274   deallocate(conn)
00275 
00276   ! take care of parallel stuff
00277 
00278   ! add entities to part, using iMesh
00279   call c_f_pointer(entsPtr, ents, [numo])
00280   call iMesh_addEntArrToSet(%VAL(imesh), ents, %VAL(numo), %VAL(part), ierr)
00281   ERROR(ierr)
00282   ! set global ids on vertices, needed for sharing between procs
00283   call iMesh_getTagHandle(%VAL(imesh), "GLOBAL_ID", tagh, ierr, %VAL(9))
00284   if (iBase_SUCCESS .ne. ierr) then
00285      ! didn't get handle, need to create the tag
00286      call iMesh_createTag(%VAL(imesh), "GLOBAL_ID", %VAL(iBase_INTEGER), tagh, ierr)
00287      ERROR(ierr)
00288   end if
00289   call iMesh_setIntArrData(%VAL(imesh), verts, %VAL(numv), %VAL(tagh), vgids, %VAL(numv), ierr)
00290   ERROR(ierr)
00291   ! now resolve shared verts and exchange ghost cells
00292   call iMeshP_syncMeshAll(%VAL(imesh), %VAL(imeshp), ierr)
00293   ERROR(ierr)
00294 
00295   call iMeshP_createGhostEntsAll(%VAL(imesh), %VAL(imeshp), %VAL(2), %VAL(1), %VAL(1), %VAL(0), ierr)
00296   ERROR(ierr)
00297 
00298   call iMesh_freeMemory(%VAL(imesh), vertsPtr);
00299   call iMesh_freeMemory(%VAL(imesh), entsPtr);
00300 
00301   return
00302 end subroutine create_mesh
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines