moab
|
00001 ! PushParMeshIntoMoabF90: push parallel mesh into moab, F90 version 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 ! To successfully link this example, you need to specify FCFLAGS that include: 00006 ! a) -DUSE_MPI, and 00007 ! b) flags required to link Fortran90 MPI programs with the C++ compiler; these flags 00008 ! can often be found on your system by inspecting the output of 'mpif90 -show' 00009 ! For example, using gcc, the link line looks like: 00010 ! make MOAB_DIR=<moab install dir> FCFLAGS="-DUSE_MPI -I/usr/lib/openmpi/include -pthread -I/usr/lib/openmpi/lib -L/usr/lib/openmpi/lib -lmpi_f90 -lmpi_f77 -lmpi -lopen-rte -lopen-pal -ldl -Wl,--export-dynamic -lnsl -lutil -lm -ldl" PushParMeshIntoMoabF90 00011 ! 00012 ! Usage: PushParMeshIntoMoab 00013 00014 program PushParMeshIntoMoab 00015 00016 use ISO_C_BINDING 00017 implicit none 00018 00019 #include "mpif.h" 00020 00021 #ifdef USE_MPI 00022 # include "iMeshP_f.h" 00023 #else 00024 # include "iMesh_f.h" 00025 #endif 00026 00027 ! declarations 00028 ! imesh is the instance handle 00029 iMesh_Instance imesh 00030 ! NUMV, NUME, NVPERE are the hardwired here; these are for the whole mesh, 00031 ! local mesh determined later 00032 integer NUMV, NUME, NVPERE 00033 parameter (NUMV = 8) ! # vertices in whole mesh 00034 parameter (NUME = 6) ! # elements in whole mesh 00035 parameter (NVPERE = 4) ! # vertices per element 00036 ! ents, verts will be arrays storing vertex/entity handles 00037 iBase_EntityHandle, pointer :: ents, verts 00038 iBase_EntitySetHandle root_set 00039 TYPE(C_PTR) :: vertsPtr, entsPtr 00040 ! storage for vertex positions, element connectivity indices, global vertex ids 00041 real*8 coords(0:3*NUMV-1) 00042 integer iconn(0:4*NUME-1), gids(0:NUMV-1) 00043 ! 00044 ! local variables 00045 integer lgids(0:NUMV-1), lconn(0:4*NUME-1) 00046 real*8 lcoords(0:3*NUMV-1) 00047 integer lnumv, lvids(0:NUMV-1), gvids(0:NUMV-1) 00048 integer lvpe, ltp ! lvpe = # vertices per entity, ltp = element type 00049 integer ic, ie, iv, istart, iend, ierr, indv, lnume, rank, sz 00050 00051 #ifdef USE_MPI 00052 ! local variables for parallel runs 00053 iMeshP_PartitionHandle imeshp 00054 ! integer MPI_COMM_WORLD 00055 #endif 00056 00057 ! vertex positions, latlon coords, (lat, lon, lev), fortran ordering 00058 ! (first index varying fastest) 00059 data coords / & 00060 0.0, -45.0, 0.0, 90.0, -45.0, 0.0, 180.0, -45.0, 0.0, 270.0, -45.0, 0.0, & 00061 0.0, 45.0, 0.0, 90.0, 45.0, 0.0, 180.0, 45.0, 0.0, 270.0, 45.0, 0.0 / 00062 00063 ! quad index numbering, each quad ccw, sides then bottom then top 00064 data iconn / & 00065 0, 1, 5, 4, & 00066 1, 2, 6, 5, & 00067 2, 3, 7, 6, & 00068 3, 0, 4, 7, & 00069 0, 3, 2, 1, & 00070 4, 5, 6, 7 / 00071 00072 data lvpe /4/ ! quads in this example 00073 data ltp / iMesh_QUADRILATERAL / ! from iBase_f.h 00074 00075 ! initialize global vertex ids 00076 do iv = 0, NUMV-1 00077 lgids(iv) = iv+1 00078 end do 00079 00080 #ifdef USE_MPI 00081 ! init the parallel partition 00082 call MPI_INIT(ierr) 00083 call MPI_COMM_SIZE(MPI_COMM_WORLD, sz, ierr) 00084 call MPI_COMM_RANK(MPI_COMM_WORLD, rank, ierr) 00085 ! compute starting/ending element numbers 00086 lnume = NUME / sz 00087 istart = rank * lnume 00088 iend = istart + lnume - 1 00089 if (rank .eq. sz-1) then 00090 iend = NUME-1 00091 lnume = iend - istart + 1 00092 endif 00093 #else 00094 ! set the starting/ending element numbers 00095 istart = 0 00096 iend = NUME-1 00097 lnume = NUME 00098 #endif 00099 00100 ! for my elements, figure out which vertices I use and accumulate local indices and coords 00101 ! lvids stores the local 0-based index for each vertex; -1 means vertex i isn't used locally 00102 ! also build up connectivity indices for local elements, in lconn 00103 do iv = 0, NUMV-1 00104 lvids(iv) = -1 00105 end do 00106 lnumv = -1 00107 do ie = istart, iend 00108 do iv = 0, lvpe-1 00109 indv = iconn(lvpe*ie + iv) 00110 if (lvids(indv) .eq. -1) then 00111 lnumv = lnumv + 1 ! increment local # verts 00112 do ic = 0, 2 ! cache local coords 00113 lcoords(3*lnumv+ic) = coords(3*indv+ic) 00114 end do 00115 lvids(indv) = lnumv 00116 gvids(lnumv) = 1+indv 00117 end if 00118 lconn(lvpe*(ie-istart)+iv) = lvids(indv) 00119 end do ! do iv 00120 end do ! do ie 00121 00122 lnumv = lnumv + 1 00123 00124 ! now create the mesh; this also initializes parallel sharing and ghost exchange 00125 imesh = 0 00126 imeshp = 0 00127 call create_mesh(imesh, imeshp, MPI_COMM_WORLD, lnumv, lnume, gvids, lvpe, ltp, lcoords, lconn, & 00128 vertsPtr, entsPtr, ierr) 00129 call c_f_pointer(vertsPtr, verts, [lnumv]) 00130 call c_f_pointer(entsPtr, ents, [lnume]) 00131 00132 ! get/report number of vertices, elements 00133 call iMesh_getRootSet(%VAL(imesh), root_set, ierr) 00134 iv = 0 00135 ie = 0 00136 #ifdef USE_MPI 00137 call iMeshP_getNumOfTypeAll(%VAL(imesh), %VAL(imeshp), %VAL(root_set), %VAL(iBase_VERTEX), iv, ierr) 00138 call iMeshP_getNumOfTypeAll(%VAL(imesh), %VAL(imeshp), %VAL(root_set), %VAL(iBase_FACE), ie, ierr) 00139 if (rank .eq. 0) then 00140 write(0,*) "Number of vertices = ", iv 00141 write(0,*) "Number of entities = ", ie 00142 endif 00143 #else 00144 call iMesh_getNumOfTypeAll(%VAL(imesh), %VAL(root_set), %VAL(iBase_VERTEX), iv, ierr) 00145 call iMesh_getNumOfTypeAll(%VAL(imesh), %VAL(root_set), %VAL(iBase_FACE), ie, ierr) 00146 write(0,*) "Number of vertices = ", iv 00147 write(0,*) "Number of entities = ", ie 00148 #endif 00149 00150 ! from here, can use verts and ents as (1-based) arrays of entity handles for input to other iMesh functions 00151 00152 call MPI_FINALIZE(ierr) 00153 stop 00154 end program PushParMeshIntoMoab 00155 00156 subroutine create_mesh( & 00157 ! interfaces 00158 imesh, imeshp, & 00159 ! input 00160 comm, numv, nume, vgids, nvpe, tp, posn, iconn, & 00161 ! output 00162 vertsPtr, entsPtr, ierr) 00163 ! 00164 ! create a mesh with numv vertices and nume elements, with elements of type tp 00165 ! vertices have positions in posn (3 coordinates each, interleaved xyzxyz...), indexed from 0 00166 ! elements have nvpe vertices per entity, with connectivity indices stored in iconn, referencing 00167 ! vertices using 0-based indices; vertex and entity handles are output in arrays passed in 00168 ! 00169 ! if imesh/imeshp are 0, imesh/imeshp are initialized in this subroutine 00170 ! 00171 00172 use ISO_C_BINDING 00173 implicit none 00174 00175 #ifdef USE_MPI 00176 # include "iMeshP_f.h" 00177 # include "mpif.h" 00178 #else 00179 # include "iMesh_f.h" 00180 #endif 00181 00182 ! subroutine arguments 00183 iMesh_Instance imesh 00184 TYPE(C_PTR) :: vertsPtr, entsPtr 00185 integer numv, nume, nvpe, vgids(0:*), iconn(0:*), ierr, tp 00186 real*8 posn(0:*) 00187 #ifdef USE_MPI 00188 iMeshP_PartitionHandle imeshp 00189 integer comm 00190 #endif 00191 00192 ! local variables 00193 integer comm_sz, comm_rank, numa, numo, iv, ie 00194 TYPE(C_PTR) :: statsPtr 00195 integer, allocatable, target :: stats(:) 00196 iBase_TagHandle tagh 00197 integer i 00198 iBase_EntityHandle, pointer :: verts(:), ents(:) 00199 iBase_EntityHandle, allocatable :: conn(:) 00200 iBase_EntitySetHandle root_set 00201 #ifdef USE_MPI 00202 IBASE_HANDLE_T mpi_comm_c 00203 TYPE(C_PTR) :: partsPtr 00204 iMeshP_PartHandle, pointer :: parts(:) 00205 iMeshP_PartHandle part 00206 integer partsa, partso 00207 #endif 00208 00209 ! create the Mesh instance 00210 if (imesh .eq. 0) then 00211 call iMesh_newMesh("MOAB", imesh, ierr) 00212 end if 00213 00214 #ifdef USE_MPI 00215 if (imeshp .eq. 0) then 00216 call iMeshP_getCommunicator(%VAL(imesh), MPI_COMM_WORLD, mpi_comm_c, ierr) 00217 call iMeshP_createPartitionAll(%VAL(imesh), %VAL(mpi_comm_c), imeshp, ierr) 00218 call iMeshP_createPart(%VAL(imesh), %VAL(imeshp), part, ierr) 00219 else 00220 partsa = 0 00221 call iMeshP_getLocalParts(%VAL(imesh), %VAL(imeshp), partsPtr, partsa, partso, ierr) 00222 call c_f_pointer(partsPtr, parts, [partso]) 00223 part = parts(1) 00224 end if 00225 call MPI_COMM_RANK(comm, comm_rank, ierr) 00226 call MPI_COMM_SIZE(comm, comm_sz, ierr) 00227 #endif 00228 00229 ! create the vertices, all in one call 00230 numa = 0 00231 call iMesh_createVtxArr(%VAL(imesh), %VAL(numv), %VAL(iBase_INTERLEAVED), posn, %VAL(3*numv), & 00232 vertsPtr, numa, numo, ierr) 00233 00234 ! fill in the connectivity array, based on indexing from iconn 00235 allocate (conn(0:nvpe*nume-1)) 00236 call c_f_pointer(vertsPtr, verts, [numv]) 00237 do i = 0, nvpe*nume-1 00238 conn(i) = verts(1+iconn(i)) 00239 end do 00240 ! create the elements 00241 numa = 0 00242 allocate(stats(0:nume-1)) 00243 statsPtr = C_LOC(stats(0)) 00244 call iMesh_createEntArr(%VAL(imesh), %VAL(tp), conn, %VAL(nvpe*nume), & 00245 entsPtr, numa, numo, statsPtr, numa, numo, ierr) 00246 deallocate(stats) 00247 deallocate(conn) 00248 00249 #ifdef USE_MPI 00250 ! take care of parallel stuff 00251 00252 ! add entities to part, using iMesh 00253 call c_f_pointer(entsPtr, ents, [numo]) 00254 call iMesh_addEntArrToSet(%VAL(imesh), ents, %VAL(numo), %VAL(part), ierr) 00255 ! set global ids on vertices, needed for sharing between procs 00256 call iMesh_getTagHandle(%VAL(imesh), "GLOBAL_ID", tagh, ierr, %VAL(9)) 00257 if (iBase_SUCCESS .ne. ierr) then 00258 ! didn't get handle, need to create the tag 00259 call iMesh_createTag(%VAL(imesh), "GLOBAL_ID", %VAL(iBase_INTEGER), tagh, ierr, %VAL(9)) 00260 end if 00261 call iMesh_setIntArrData(%VAL(imesh), verts, %VAL(numv), %VAL(tagh), vgids, %VAL(numv), ierr) 00262 00263 ! now resolve shared verts and exchange ghost cells 00264 call iMeshP_syncMeshAll(%VAL(imesh), %VAL(imeshp), ierr) 00265 call iMesh_getRootSet(%VAL(imesh), root_set, ierr) 00266 call iMeshP_getNumOfTypeAll(%VAL(imesh), %VAL(imeshp), %VAL(root_set), %VAL(iBase_VERTEX), iv, ierr) 00267 call iMeshP_getNumOfTypeAll(%VAL(imesh), %VAL(imeshp), %VAL(root_set), %VAL(iBase_FACE), ie, ierr) 00268 if (comm_rank .eq. 0) then 00269 write(0,*) "After syncMeshAll:" 00270 write(0,*) " Number of vertices = ", iv 00271 write(0,*) " Number of entities = ", ie 00272 endif 00273 00274 call iMeshP_createGhostEntsAll(%VAL(imesh), %VAL(imeshp), %VAL(2), %VAL(1), %VAL(1), %VAL(0), ierr) 00275 if (comm_rank .eq. 0) then 00276 write(0,*) "After createGhostEntsAll:" 00277 write(0,*) " Number of vertices = ", iv 00278 write(0,*) " Number of entities = ", ie 00279 endif 00280 #endif 00281 00282 return 00283 end subroutine create_mesh