MOAB: Mesh Oriented datABase
(version 5.3.1)
|
00001 00002 SUBROUTINE errorout(ierr, message) 00003 integer ierr 00004 character*(*) message 00005 if (ierr.ne.0) then 00006 print *, message 00007 call exit (1) 00008 end if 00009 return 00010 end 00011 00012 #include "moab/MOABConfig.h" 00013 program fdriver 00014 #ifdef MOAB_HAVE_MPI 00015 include 'mpif.h' 00016 #endif 00017 integer ierr, num_procs, my_id 00018 integer pid 00019 integer compid 00020 character*10 appname 00021 character*200 filename 00022 character*14 fname 00023 character*150 readopts 00024 integer ngv, nge, ndim, nparts 00025 integer nghlay 00026 integer nverts(3), nelem(3), nblocks(3), nsbc(3), ndbc(3) 00027 C large enough work arrays 00028 integer iwork(100000) 00029 real*8 dwork(100000) 00030 C indices in work arrays for vertex ids, ranks, coordinates 00031 integer vID, vRA, vCO 00032 C 00033 C indices for free memory (index in integer or double work arrays) 00034 integer ifree, dfree 00035 C size of coordinates array 00036 integer nCO 00037 integer bID 00038 C for some tags in the file 00039 character*20 tagname1, tagname2 00040 integer tagtype(2), enttype(2), num_co 00041 integer tagindex(2) 00042 integer stags(2), itags(2), ntsync 00043 00044 integer eRA, beID, eID 00045 C vertice per element, number of elements in block 00046 integer vpere, nebl, blockID 00047 C iWORK(eCO) start for connectivity 00048 integer sizeconn, eCO 00049 C IWORK(egID) , IWORK(elID) starts for global el ID, local elem ID 00050 integer egID, elID, eOWN 00051 integer iTAG, dTAG 00052 00053 C indices for surface BC element, reference surf BC, value 00054 integer isBC, irBC, ivBC 00055 00056 character*100 outfile, wopts 00057 my_id = 0 00058 fname = '/io/p8ex1.h5m' // char(0) 00059 #ifdef MOAB_HAVE_MPI 00060 call MPI_INIT ( ierr ) 00061 call errorout(ierr, 'fail to initialize MPI') 00062 00063 c find out MY process ID, and how many processes were started. 00064 00065 call MPI_COMM_RANK (MPI_COMM_WORLD, my_id, ierr) 00066 call errorout(ierr, 'fail to get MPI rank') 00067 00068 call MPI_COMM_SIZE (MPI_COMM_WORLD, num_procs, ierr) 00069 call errorout(ierr, 'fail to get MPI size') 00070 00071 if (my_id .eq. 0) then 00072 print *, " I'm process ", my_id, " out of ", 00073 & num_procs, " processes." 00074 00075 end if 00076 #endif 00077 ierr = iMOAB_InitializeFortran() 00078 call errorout(ierr, 'fail to initialize iMOAB') 00079 00080 appname = 'PROTEUS'//CHAR(0) 00081 C give a unique external id; here we just use 8? 00082 compid = 8 00083 #ifdef MOAB_HAVE_MPI 00084 ierr = iMOAB_RegisterApplicationFortran(appname, MPI_COMM_WORLD, 00085 & compid, pid) 00086 #else 00087 ierr = iMOAB_RegisterApplicationFortran(appname, compid, 00088 & pid) 00089 #endif 00090 call errorout(ierr, 'fail to register application') 00091 C 00092 filename = MESHDIR2//fname 00093 ierr = iMOAB_ReadHeaderInfo ( filename, ngv, nge, ndim, nparts) 00094 call errorout(ierr, 'fail to read header info') 00095 00096 if (0.eq.my_id) then 00097 print *, filename, ' has ', nparts, ' parts in partition', 00098 & ngv, ' vertices ', nge, ' elements of dimension ', ndim 00099 endif 00100 #ifdef MOAB_HAVE_MPI 00101 readopts='PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;'// 00102 & 'PARALLEL_RESOLVE_SHARED_ENTS'//CHAR(0) 00103 #else 00104 readopts=CHAR(0) 00105 #endif 00106 00107 C number of ghost layers needed by application 00108 nghlay = 1 00109 00110 C now let us load the mesh in parallel 00111 ierr = iMOAB_LoadMesh(pid, filename, readopts, nghlay) 00112 call errorout(ierr, 'fail to read file in parallel') 00113 00114 C number of vertices/elements/blocks/sidesets in the mesh 00115 ierr = iMOAB_GetMeshInfo(pid, nverts, nelem, nblocks, nsbc, ndbc) 00116 call errorout(ierr, 'fail to get mesh info') 00117 00118 vID = 1 00119 ierr = iMOAB_GetVertexID(pid, nverts(3), IWORK(vID) ) 00120 call errorout(ierr, 'failed to get vertex id info') 00121 00122 vRA = vID + nverts(3) 00123 ierr = iMOAB_GetVertexOwnership(pid, nverts(3), IWORK(vRA) ) 00124 call errorout(ierr, 'failed to get vertex owner ranks') 00125 00126 ifree = vRA + nverts(3) 00127 00128 C double * coords = (double*) malloc(3*nverts[2]*sizeof(double)); 00129 vCO = 1 00130 nCO = 3 * nverts(3) 00131 00132 ierr = iMOAB_GetVisibleVerticesCoordinates(pid, nCO, DWORK(vCO)) 00133 call errorout(ierr, 'failed to get coordinates') 00134 dfree = vCO + 3 * nverts(3) 00135 00136 bID = ifree 00137 ierr = iMOAB_GetBlockID(pid, nblocks(3), IWORK(bID)) 00138 call errorout(ierr, 'failed to get block info') 00139 ifree = ifree + nblocks(3) 00140 00141 C the 2 tags used in this example exist in the file, already 00142 C first tag, INTFIELD is on vertices, integer 00143 C second tag DFIELD is on elements, double 00144 tagtype(1)=0 !dense, int 00145 tagtype(2)=1 !dense, double 00146 enttype(1)=0 ! on verts 00147 enttype(2)=1 ! on elem 00148 num_co = 1 00149 tagname1 ='INTFIELD'//CHAR(0) 00150 ierr = iMOAB_DefineTagStorage(pid, tagname1, tagtype(1), num_co, 00151 & tagindex(1) ) 00152 call errorout(ierr, 'failed to get tag INTFIELD') 00153 00154 tagname2 ='DFIELD'//CHAR(0) 00155 ierr = iMOAB_DefineTagStorage(pid, tagname2, tagtype(2), num_co, 00156 & tagindex(2) ) 00157 call errorout(ierr, 'failed to get tag DFIELD') 00158 00159 C synchronize one of the tags only, just to see what happens 00160 C put in the sync array just first tag index (INTFIELD) 00161 ntsync =1 00162 stags(1) = tagIndex(1) 00163 itags(1) = tagType(1) 00164 00165 ierr = iMOAB_SynchronizeTags(pid, ntsync, stags, itags ) 00166 call errorout(ierr, 'failed to sync tag INTFIELD') 00167 00168 C start printing some information, retrieved from each task 00169 do irk=0, num_procs-1 00170 if (irk .eq. my_id) then 00171 00172 C printf some of the block info */ 00173 print *, 'on rank ', my_id, ' there are ' 00174 print *, nverts(3), ' visible vertices of which ',nverts(1), 00175 & ' local ', nverts(2), ' ghost' 00176 print *, nblocks(3), ' visible blocks' 00177 print *, nsbc(3), ' visible Neumann BCs' 00178 print *, ndbc(3), ' visible dirichlet BCs' 00179 00180 print *, 'on rank ', my_id, ' vertex info:' 00181 do i=1,nverts(3) 00182 write(*, 100) i, IWORK(vRA+i-1), IWORK( vID+i-1), 00183 & DWORK(vCO+3*i-3), DWORK(vCO+3*i-2), DWORK(vCO+3*i-1) 00184 100 FORMAT(' vertex local id ', I3, ' rank ID', I3, ' global ID:' 00185 & , I3, ' coords:', 3F11.3) 00186 enddo 00187 eID = ifree 00188 beID = eID + nelem(3) 00189 eRA = beID + nelem(3) 00190 ierr = iMOAB_GetVisibleElementsInfo(pid, nelem(3),IWORK(eID), 00191 & IWORK(eRA), IWORK(beID) ) 00192 call errorout(ierr, 'failed to get all elem info') 00193 ifree = eRA + nelem(3) 00194 do i=1, nelem(3) 00195 write(*, 101) IWORK(eID+i-1), IWORK(eRA+i-1), IWORK(beID+i-1) 00196 101 FORMAT( ' global ID ', I5, ' rank: ', 00197 & I3, ' block ID: ' I4) 00198 enddo 00199 00200 do i=1,nblocks(3) 00201 00202 print *,' block index:', i, ' block ID ', IWORK(bID+i-1) 00203 blockID = IWORK(bID+i-1) 00204 ierr = iMOAB_GetBlockInfo(pid, blockID , vpere, nebl) 00205 call errorout(ierr, 'failed to elem block info') 00206 print *, ' has' , nebl, ' elements with ', vpere, 'verts' 00207 00208 sizeconn = nebl * vpere 00209 eCO = ifree 00210 ierr = iMOAB_GetBlockElementConnectivities(pid, blockID, 00211 & sizeconn, IWORK(eCO) ) 00212 call errorout(ierr, 'failed to get block elem connectivity') 00213 00214 ifree = ifree + sizeconn 00215 eOWN = ifree 00216 ierr = iMOAB_GetElementOwnership(pid, blockID, nebl, 00217 & IWORK(eOWN)) 00218 call errorout(ierr, 'failed to get block elem ownership') 00219 ifree = ifree+nebl 00220 00221 egID = ifree 00222 elID = ifree + nebl 00223 ierr = iMOAB_GetElementID(pid, blockID, nebl, 00224 & IWORK(egID), IWORK(elID) ) 00225 call errorout(ierr, 'failed to get block elem IDs') 00226 ifree = elID + nebl 00227 00228 do j=1, nebl 00229 write (*, 102) j, IWORK(eOWN+j-1),IWORK(egID+j-1), 00230 & IWORK(elID+j-1), (IWORK(eCO-1+(j-1)*vpere+k), k=1,vpere) 00231 102 FORMAT(' elem ', I3, ' owned by', I3, ' gid:', I3, ' lid:', 00232 & I3, ' : ', 10I5) 00233 enddo 00234 enddo 00235 00236 C query int tag values on vertices 00237 iTAG= ifree 00238 00239 ierr = iMOAB_GetIntTagStorage(pid, tagname1, nverts(3), 00240 & enttype(1), IWORK(iTAG) ) 00241 call errorout(ierr, 'failed to get INTFIELD tag') 00242 ifree = iTAG + nverts(3) 00243 print * , 'INTFIELD tag values' 00244 write(*, 103) (IWORK(iTAG+k-1), k=1,nverts(3) ) 00245 103 FORMAT (10I8) 00246 00247 dTAG = dfree 00248 C query double tag values on elements 00249 00250 ierr = iMOAB_GetDoubleTagStorage(pid, tagname2, nelem(3), 00251 & entType(2), DWORK(dTAG) ) 00252 call errorout(ierr, 'failed to get DFIELD tag') 00253 dfree = dTAG + nelem(3) 00254 print *, 'DFIELD tag values: (not exchanged) ' 00255 write (*, 104) (DWORK(dTAG+k-1), k=1,nelem(3) ) 00256 104 FORMAT ( 10F7.2 ) 00257 00258 C query surface BCs 00259 isBC = ifree 00260 irBC = isBC + nsbc(3) 00261 ivBC = irBC + nsbc(3) 00262 00263 ierr = iMOAB_GetPointerToSurfaceBC(pid, nsbc(3), IWORK(isBC), 00264 & IWORK(irBC), IWORK(ivBC)) 00265 call errorout(ierr, 'failed to get surf boundary conditions') 00266 ifree = ivBC + nsbc(3) 00267 print * , 'Surface boundary conditions ' 00268 write (*, 105) (IWORK(isBC+k-1), 00269 & IWORK(irBC+k-1), IWORK(ivBC+k-1), k=1, nsbc(3)) 00270 105 FORMAT (' elem localID: ', I3, ' side:', I1, ' val:', I4) 00271 00272 C query vertex BCs 00273 iveBC = ifree 00274 ivaBC = iveBC + ndbc(3) 00275 00276 ierr = iMOAB_GetPointerToVertexBC(pid, ndbc(3), IWORK(iveBC), 00277 & IWORK(ivaBC)) 00278 call errorout(ierr, 'failed to get vertex boundary conditions') 00279 ifree = ivaBC + ndbc(3) 00280 00281 print *, ' Vertex boundary conditions:' 00282 write (*, 106) (IWORK(iveBC+k-1), 00283 & IWORK(ivaBC+k-1), k=1, ndbc(3)) 00284 106 FORMAT (' vertex: ', I3, ' BC:', I6 ) 00285 00286 endif 00287 #ifdef MOAB_HAVE_MPI 00288 call MPI_Barrier(MPI_COMM_WORLD, ierr) 00289 #endif 00290 call errorout(ierr, 'fail at barrier') 00291 00292 enddo 00293 00294 outfile = 'fnew2.h5m'//CHAR(0) 00295 #ifdef MOAB_HAVE_MPI 00296 wopts = 'PARALLEL=WRITE_PART'//CHAR(0) 00297 #else 00298 wopts=CHAR(0) 00299 #endif 00300 C write out the mesh file to disk 00301 ierr = iMOAB_WriteMesh(pid, outfile, wopts) 00302 call errorout(ierr, 'fail to write the mesh file') 00303 00304 C all done. de-register and finalize 00305 ierr = iMOAB_DeregisterApplication(pid) 00306 call errorout(ierr, 'fail to deregister application') 00307 00308 ierr = iMOAB_Finalize() 00309 call errorout(ierr, 'fail to finalize iMOAB') 00310 00311 #ifdef MOAB_HAVE_MPI 00312 call MPI_FINALIZE ( ierr ) 00313 #endif 00314 call errorout(ierr, 'fail to finalize MPI') 00315 00316 stop 00317 end