imoab_testF.F90

Go to the documentation of this file.

 
SUBROUTINE errorout(ierr, message)
integer ierr
character*(*) message
if (ierr.ne.0) then
 print *, message
 call exit (1)
end if
return
end
 
program fdriver
 
use iso_c_binding
use imoab
 
#include "moab/MOABConfig.h"
#ifdef MOAB_HAVE_MPI
include 'mpif.h'
#endif
 
#ifndef MOAB_MESH_DIR
#error Specify MOAB_MESH_DIR path
#endif
 
 integer ierr, num_procs, my_id
 integer pid
 integer compid
 character :: appname*32
 character :: filename*1024
 character :: fname*1024
 character :: readopts*1024
 integer ngv, nge, ndim, nparts
 integer nghlay
 integer nverts(3), nelem(3), nblocks(3), nsbc(3), ndbc(3)
 ! large enough work arrays
 integer iwork(100000)
 double precision dwork(100000)
 ! indices in work arrays for vertex ids, ranks, coordinates
 integer vid, vra, vco
 !
 ! indices for free memory (index in integer or double work arrays)
 integer ifree, dfree
 ! size of coordinates array
 integer nco
 integer bid
 ! for some tags in the file
 character :: tagname1*128, tagname2*128
 integer tagtype(2), enttype(2), num_co
 integer tagindex(2)
 integer ntsync
 
 integer era, beid, eid
 ! vertices per element, number of elements in block
 integer vpere, nebl, blockid
 ! iWORK(eCO) start for connectivity
 integer sizeconn, eco
 ! IWORK(egID) , IWORK(elID) starts for global el ID, local elem ID
 integer egid, elid, eown
 integer itag, dtag
 
 ! indices for surface BC element, reference surf BC, value
 integer isbc, irbc, ivbc
 
 character outfile*1024, wopts*1024
 my_id = 0
 fname = 'unittest/io/p8ex1.h5m'//c_null_char
 
#ifdef MOAB_HAVE_MPI
 call mpi_init ( ierr )
 call errorout(ierr, 'fail to initialize MPI')
 
 ! find out MY process ID, and how many processes were started.
 
 call mpi_comm_rank (mpi_comm_world, my_id, ierr)
 call errorout(ierr, 'fail to get MPI rank')
 
 call mpi_comm_size (mpi_comm_world, num_procs, ierr)
 call errorout(ierr, 'fail to get MPI size')
 
 if (my_id .eq. 0) then
 print *, " I'm process ", my_id, " out of ", num_procs, " processes."
 end if
#endif
 
 ierr = imoab_initialize()
 call errorout(ierr, 'fail to initialize iMOAB')
 
 appname = 'PROTEUS'//c_null_char
 ! give a unique external id; here we just use 8?
 compid = 8
 ! number of ghost layers needed by application
 nghlay = 0
 
 
#ifdef MOAB_HAVE_MPI
 ierr = imoab_registerapplication(appname, mpi_comm_world, compid, pid)
#else
 ierr = imoab_registerapplication(appname, compid, pid)
#endif
 call errorout(ierr, 'fail to register application')
 
 filename = moab_mesh_dir//fname
 ierr = imoab_readheaderinfo( filename, ngv, nge, ndim, nparts)
 call errorout(ierr, 'fail to read header info')
 
 if (0.eq.my_id) then
 print *, filename, ' has ', nparts, ' parts in partition', ngv, ' vertices ', nge, ' elements of dimension ', ndim
 endif
#ifdef MOAB_HAVE_MPI
 nghlay = 1
 readopts='PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;'// &
 'PARALLEL_RESOLVE_SHARED_ENTS'//c_null_char
#else
 readopts=c_null_char
#endif
 
 ! now let us load the mesh in parallel
 ierr = imoab_loadmesh(pid, filename, readopts, nghlay)
 call errorout(ierr, 'fail to read file in parallel')
 
 ! number of vertices/elements/blocks/sidesets in the mesh
 ierr = imoab_getmeshinfo(pid, nverts, nelem, nblocks, nsbc, ndbc)
 call errorout(ierr, 'fail to get mesh info')
 
 vid = 1
 ierr = imoab_getvertexid(pid, nverts(3), iwork(vid) )
 call errorout(ierr, 'failed to get vertex id info')
 
 vra = vid + nverts(3)
 ierr = imoab_getvertexownership(pid, nverts(3), iwork(vra) )
 call errorout(ierr, 'failed to get vertex owner ranks')
 
 ifree = vra + nverts(3)
 
 ! double * coords = (double*) malloc(3*nverts[2]*sizeof(double));
 vco = 1
 nco = 3 * nverts(3)
 
 ierr = imoab_getvisibleverticescoordinates(pid, nco, dwork(vco))
 print *, 'nCO = ', nco
 print *, 'ierr = ', ierr
 call errorout(ierr, 'failed to get coordinates')
 dfree = vco + 3 * nverts(3)
 
 bid = ifree
 ierr = imoab_getblockid(pid, nblocks(3), iwork(bid))
 call errorout(ierr, 'failed to get block info')
 ifree = ifree + nblocks(3)
 
 ! the 2 tags used in this example exist in the file, already
 ! first tag, INTFIELD is on vertices, integer
 ! second tag DFIELD is on elements, double
 tagtype(1)=0 !dense, int
 tagtype(2)=1 !dense, double
 enttype(1)=0 ! on verts
 enttype(2)=1 ! on elem
 num_co = 1
 tagname1 ='INTFIELD'//c_null_char
 ierr = imoab_definetagstorage(pid, tagname1, tagtype(1), num_co, tagindex(1) )
 call errorout(ierr, 'failed to get tag INTFIELD')
 
 tagname2 ='DFIELD'//c_null_char
 ierr = imoab_definetagstorage(pid, tagname2, tagtype(2), num_co, tagindex(2) )
 call errorout(ierr, 'failed to get tag DFIELD')
 
 ! synchronize one of the tags only, just to see what happens
 ! put in the sync array just first tag index (INTFIELD)
 ntsync =2
 
 ierr = imoab_synchronizetags(pid, ntsync, tagindex, tagtype )
 call errorout(ierr, 'failed to sync tag INTFIELD')
 
 ! start printing some information, retrieved from each task
 do irk=0, num_procs-1
 if (irk .eq. my_id) then
 
 ! printf some of the block info */
 print *, 'on rank ', my_id, ' there are '
 print *, nverts(3), ' visible vertices of which ',nverts(1), ' local ', nverts(2), ' ghost'
 print *, nblocks(3), ' visible blocks'
 print *, nsbc(3), ' visible Neumann BCs'
 print *, ndbc(3), ' visible dirichlet BCs'
 
 print *, 'on rank ', my_id, ' vertex info:'
 do i=1,nverts(3)
 write(*, 100) i, iwork(vra+i-1), iwork( vid+i-1), dwork(vco+3*i-3), dwork(vco+3*i-2), dwork(vco+3*i-1)
100 FORMAT(' vertex local id ', i3, ' rank ID', i3, ' global ID:', i3, ' coords:', 3f11.3)
 enddo
 
 eid = ifree
 beid = eid + nelem(3)
 era = beid + nelem(3)
 ierr = imoab_getvisibleelementsinfo(pid, nelem(3),iwork(eid), iwork(era), iwork(beid) )
 call errorout(ierr, 'failed to get all elem info')
 ifree = era + nelem(3)
 do i=1, nelem(3)
 write(*, 101) iwork(eid+i-1), iwork(era+i-1), iwork(beid+i-1)
101 FORMAT( ' global ID ', i5, ' rank: ', i3, ' block ID: ', i4)
 enddo
 
 do i=1,nblocks(3)
 
 print *,' block index:', i, ' block ID ', iwork(bid+i-1)
 blockid = iwork(bid+i-1)
 ierr = imoab_getblockinfo(pid, blockid , vpere, nebl)
 call errorout(ierr, 'failed to elem block info')
 print *, ' has' , nebl, ' elements with ', vpere, 'verts'
 
 sizeconn = nebl * vpere
 eco = ifree
 ierr = imoab_getblockelementconnectivities(pid, blockid, sizeconn, iwork(eco) )
 print *, ierr
 call errorout(ierr, 'failed to get block elem connectivity')
 
 ifree = ifree + sizeconn
 eown = ifree
 ierr = imoab_getelementownership(pid, blockid, nebl, iwork(eown))
 call errorout(ierr, 'failed to get block elem ownership')
 ifree = ifree+nebl
 
 egid = ifree
 elid = ifree + nebl
 ierr = imoab_getelementid(pid, blockid, nebl, iwork(egid), iwork(elid) )
 call errorout(ierr, 'failed to get block elem IDs')
 ifree = elid + nebl
 
 do j=1, nebl
 write (*, 102) j, iwork(eown+j-1),iwork(egid+j-1), iwork(elid+j-1), (iwork(eco-1+(j-1)*vpere+k), k=1,vpere)
102 FORMAT(' elem ', i3, ' owned by', i3, ' gid:', i3, ' lid:', i3, ' : ', 10i5)
 enddo
 enddo
 
 ! query int tag values on vertices
 itag= ifree
 
 ierr = imoab_getinttagstorage(pid, tagname1, nverts(3), enttype(1), iwork(itag) )
 call errorout(ierr, 'failed to get INTFIELD tag')
 ifree = itag + nverts(3)
 print * , 'INTFIELD tag values'
 write(*, 103) (iwork(itag+k-1), k=1,nverts(3) )
103 FORMAT (10i8)
 
 dtag = dfree
 ! query double tag values on elements
 
 ierr = imoab_getdoubletagstorage(pid, tagname2, nelem(3), enttype(2), dwork(dtag) )
 call errorout(ierr, 'failed to get DFIELD tag')
 dfree = dtag + nelem(3)
 print *, 'DFIELD tag values: (not exchanged) '
 write (*, 104) (dwork(dtag+k-1), k=1,nelem(3) )
104 FORMAT ( 10f7.2 )
 
 ! query surface BCs
 isbc = ifree
 irbc = isbc + nsbc(3)
 ivbc = irbc + nsbc(3)
 
 ierr = imoab_getpointertosurfacebc(pid, nsbc(3), iwork(isbc),iwork(irbc), iwork(ivbc))
 call errorout(ierr, 'failed to get surf boundary conditions')
 ifree = ivbc + nsbc(3)
 print * , 'Surface boundary conditions '
 write (*, 105) (iwork(isbc+k-1),iwork(irbc+k-1), iwork(ivbc+k-1), k=1, nsbc(3))
105 FORMAT (' elem localID: ', i3, ' side:', i1, ' val:', i4)
 
 ! query vertex BCs
 ivebc = ifree
 ivabc = ivebc + ndbc(3)
 
 ierr = imoab_getpointertovertexbc(pid, ndbc(3), iwork(ivebc),iwork(ivabc))
 call errorout(ierr, 'failed to get vertex boundary conditions')
 ifree = ivabc + ndbc(3)
 
 print *, ' Vertex boundary conditions:'
 write (*, 106) (iwork(ivebc+k-1),iwork(ivabc+k-1), k=1, ndbc(3))
106 FORMAT (' vertex: ', i3, ' BC:', i6 )
 
 endif
 
#ifdef MOAB_HAVE_MPI
 call mpi_barrier(mpi_comm_world, ierr)
#endif
 call errorout(ierr, 'fail at barrier')
enddo
 
 outfile = 'fnew2.h5m'//c_null_char
#ifdef MOAB_HAVE_MPI
 wopts = 'PARALLEL=WRITE_PART'//c_null_char
#else
 wopts=c_null_char
#endif
 
 ! write out the mesh file to disk
 ierr = imoab_writemesh(pid, outfile, wopts)
 call errorout(ierr, 'fail to write the mesh file')
 
 ! all done. de-register and finalize
 ierr = imoab_deregisterapplication(pid)
 call errorout(ierr, 'fail to deregister application')
 
 ierr = imoab_finalize()
 call errorout(ierr, 'fail to finalize iMOAB')
 
#ifdef MOAB_HAVE_MPI
 call mpi_finalize ( ierr )
#endif
 call errorout(ierr, 'fail to finalize MPI')
 
 stop
end