MOAB: Mesh Oriented datABase  (version 5.2.1)
imoab_testF.F
Go to the documentation of this file.
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_RegisterFortranApplication(appname, MPI_COMM_WORLD,
00085      & compid, pid)
00086 #else
00087       ierr = iMOAB_RegisterFortranApplication(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   
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines