moab
DirectAccessNoHolesF90.F90

Use direct access to MOAB data to avoid calling through API, in Fortran90
This example creates a 1d row of quad elements, such that all quad and vertex handles are contiguous in the handle space and in the database. Then it shows how to get access to pointers to MOAB-native data for vertex coordinates, quad connectivity, and tag storage (vertex to quad adjacency lists aren't accessible from Fortran because they are std::vector's). This allows applications to access this data directly without going through MOAB's API. In cases where the mesh is not changing (or only mesh vertices are moving), this can save significant execution time in applications.

Using direct access (or MOAB in general) from Fortran is complicated in two specific ways: 1) There is no way to use arrays with specified dimension starting/ending values with ISO_C_BINDING; therefore, all arrays passed back from MOAB/iMesh must use 1-based indices; this makes it a bit more difficult to traverse the direct arrays. In this example, I explicitly use indices that look like my_array(1+v_ind...) to remind users of this. 2) Arithmetic on handles is complicated by the fact that Fortran has no unsigned integer type. I get around this by assigning to a C-typed variable first, before the handle arithmetic. This seems to work fine. C-typed variables are part of the Fortran95 standard.

---------------------- | | | | | | | | ... | | | | ----------------------

  1. Initialize MOAB
  2. Create a quad mesh, as depicted above
  3. Create 2 dense tags (tag1, tag2) for avg position to assign to quads, and # verts per quad (tag3)
  4. Get connectivity, coordinate, tag1 iterators
  5. Iterate through quads, computing midpoint based on vertex positions, set on quad-based tag1
  6. Iterate through vertices, summing positions into tag2 on connected quads and incrementing vertex count
  7. Iterate through quads, normalizing tag2 by vertex count and comparing values of tag1 and tag2

To compile:
make DirectAccessNoHolesF90 MOAB_DIR=<installdir>
To run: ./DirectAccessNoHolesF90 [-nquads <# quads>]

00001 
00039 
00040 #define CHECK(a) \
00041   if (a .ne. iBase_SUCCESS) call exit(a)
00042 
00043 module freem
00044   interface
00045      subroutine c_free(ptr) bind(C, name='free')
00046        use ISO_C_BINDING
00047        type(C_ptr), value, intent(in) :: ptr
00048      end subroutine c_free
00049   end interface
00050 end module freem
00051 
00052 program DirectAccessNoHolesF90
00053   use ISO_C_BINDING
00054   use freem
00055   implicit none
00056 #include "iMesh_f.h"
00057 
00058   ! declarations
00059   iMesh_Instance imesh
00060   iBase_EntitySetHandle root_set
00061   integer ierr, nquads, nquads_tmp, nverts
00062   iBase_TagHandle tag1h, tag2h, tag3h
00063   TYPE(C_PTR) verts_ptr, quads_ptr, conn_ptr, x_ptr, y_ptr, z_ptr, tag1_ptr, tag2_ptr, tag3_ptr, 
00064        tmp_quads_ptr, offsets_ptr  ! pointers that are equivalence'd to arrays using c_f_pointer
00065   real(C_DOUBLE), pointer :: x(:), y(:), z(:), tag1(:), tag2(:) ! arrays equivalence'd to pointers
00066   integer, pointer :: tag3(:), offsets(:)                       ! arrays equivalence'd to pointers
00067   iBase_EntityHandle, pointer :: quads(:), verts(:), conn(:), tmp_quads(:)  ! arrays equivalence'd to pointers
00068   iBase_EntityArrIterator viter, qiter
00069   integer(C_SIZE_T) :: startv, startq, v_ind, e_ind  ! needed to do handle arithmetic
00070   integer count, vpere, e, i, j, v, offsets_size, tmp_quads_size, n_dis
00071   character*120 opt
00072 
00073   ! initialize interface and get root set
00074   call iMesh_newMesh("", imesh, ierr); CHECK(ierr)
00075   call iMesh_getRootSet(%VAL(imesh), root_set, ierr); CHECK(ierr)
00076 
00077   ! create mesh
00078   nquads_tmp = 1000
00079   call create_mesh_no_holes(imesh, nquads_tmp, ierr); CHECK(ierr)
00080 
00081   ! get all vertices and quads as arrays
00082   nverts = 0
00083   call iMesh_getEntities(%VAL(imesh), %VAL(root_set), %VAL(0), %VAL(iBase_VERTEX), &
00084        verts_ptr, nverts, nverts, ierr); CHECK(ierr)
00085   call c_f_pointer(verts_ptr, verts, [nverts])
00086   nquads = 0
00087   call iMesh_getEntities(%VAL(imesh), %VAL(root_set), %VAL(2), %VAL(iMesh_QUADRILATERAL), &
00088        quads_ptr, nquads, nquads, ierr); CHECK(ierr)
00089   call c_f_pointer(quads_ptr, quads, [nquads])
00090 
00091   ! First vertex/quad is at start of vertex/quad list, and is offset for vertex/quad index calculation
00092   startv = verts(1); startq = quads(1)
00093   call c_free(quads_ptr)  ! free memory we allocated in MOAB
00094 
00095   ! create tag1 (element-based avg), tag2 (vertex-based avg)
00096   opt = 'moab:TAG_STORAGE_TYPE=DENSE moab:TAG_DEFAULT_VALUE=0'
00097   call iMesh_createTagWithOptions(%VAL(imesh), 'tag1', opt, %VAL(3), %VAL(iBase_DOUBLE), &
00098        tag1h, ierr, %VAL(4), %VAL(56)); CHECK(ierr)
00099   call iMesh_createTagWithOptions(%VAL(imesh), 'tag2', opt, %VAL(3), %VAL(iBase_DOUBLE), &
00100        tag2h, ierr, %VAL(4), %VAL(56)); CHECK(ierr)
00101 
00102   ! Get iterator to verts, then pointer to coordinates
00103   call iMesh_initEntArrIter(%VAL(imesh), %VAL(root_set), %VAL(iBase_VERTEX), %VAL(iMesh_ALL_TOPOLOGIES), &
00104        %VAL(nverts), %VAL(0), viter, ierr); CHECK(ierr)
00105   call iMesh_coordsIterate(%VAL(imesh), %VAL(viter), x_ptr, y_ptr, z_ptr, count, ierr); CHECK(ierr)
00106   CHECK(count-nverts) ! check that all verts are in one contiguous chunk
00107   call c_f_pointer(x_ptr, x, [nverts]); call c_f_pointer(y_ptr, y, [nverts]); call c_f_pointer(z_ptr, z, [nverts])
00108 
00109   ! Get iterator to quads, then pointers to connectivity and tags
00110   call iMesh_initEntArrIter(%VAL(imesh), %VAL(root_set), %VAL(iBase_FACE), %VAL(iMesh_ALL_TOPOLOGIES), &
00111        %VAL(nquads), %VAL(0), qiter, ierr); CHECK(ierr)
00112   call iMesh_connectIterate(%VAL(imesh), %VAL(qiter), conn_ptr, vpere, count, ierr); CHECK(ierr)
00113   call c_f_pointer(conn_ptr, conn, [vpere*nquads])
00114   call iMesh_tagIterate(%VAL(imesh), %VAL(tag1h), %VAL(qiter), tag1_ptr, count, ierr); CHECK(ierr)
00115   call c_f_pointer(tag1_ptr, tag1, [3*nquads])
00116   call iMesh_tagIterate(%VAL(imesh), %VAL(tag2h), %VAL(qiter), tag2_ptr, count, ierr); CHECK(ierr)
00117   call c_f_pointer(tag2_ptr, tag2, [3*nquads])
00118   
00119   ! iterate over elements, computing tag1 from coords positions; use (1+... in array indices for 1-based indexing
00120   do i = 0, nquads-1
00121      tag1(1+3*i+0) = 0.0; tag1(1+3*i+1) = 0.0; tag1(1+3*i+2) = 0.0 ! initialize position for this element
00122      do j = 0, vpere-1 ! loop over vertices in this quad
00123         v_ind = conn(1+vpere*i+j) ! assign to v_ind to allow handle arithmetic
00124         v_ind = v_ind - startv ! vert index is just the offset from start vertex
00125         tag1(1+3*i+0) = tag1(1+3*i+0) + x(1+v_ind)
00126         tag1(1+3*i+1) = tag1(1+3*i+1) + y(1+v_ind) ! sum vertex positions into tag1...
00127         tag1(1+3*i+2) = tag1(1+3*i+2) + z(1+v_ind)
00128      end do
00129      tag1(1+3*i+0) = tag1(1+3*i+0) / vpere;
00130      tag1(1+3*i+1) = tag1(1+3*i+1) / vpere; ! then normalize
00131      tag1(1+3*i+2) = tag1(1+3*i+2) / vpere;
00132   end do ! loop over elements in chunk
00133     
00134   ! Iterate through vertices, summing positions into tag2 on connected elements; get adjacencies all at once,
00135   ! could also get them per vertex (time/memory tradeoff)
00136   tmp_quads_size = 0
00137   offsets_size = 0
00138   call iMesh_getEntArrAdj(%VAL(imesh), verts, %VAL(nverts), %VAL(iMesh_QUADRILATERAL), &
00139        tmp_quads_ptr, tmp_quads_size, tmp_quads_size, &
00140        offsets_ptr, offsets_size, offsets_size, ierr); CHECK(ierr) 
00141   call c_f_pointer(tmp_quads_ptr, tmp_quads, [tmp_quads_size])
00142   call c_f_pointer(offsets_ptr, offsets, [offsets_size])
00143   call c_free(verts_ptr)  ! free memory we allocated in MOAB
00144   do v = 0, nverts-1
00145      do e = 1+offsets(1+v), 1+offsets(1+v+1)-1  ! 1+ because offsets are in terms of 0-based arrays
00146         e_ind = tmp_quads(e) ! assign to e_ind to allow handle arithmetic
00147         e_ind = e_ind - startq ! e_ind is the quad handle minus the starting quad handle
00148         tag2(1+3*e_ind+0) = tag2(1+3*e_ind+0) + x(1+v)
00149         tag2(1+3*e_ind+1) = tag2(1+3*e_ind+1) + y(1+v)   ! tag on each element is 3 doubles, x/y/z
00150         tag2(1+3*e_ind+2) = tag2(1+3*e_ind+2) + z(1+v)
00151      end do
00152   end do
00153   call c_free(tmp_quads_ptr)  ! free memory we allocated in MOAB
00154   call c_free(offsets_ptr)    ! free memory we allocated in MOAB
00155 
00156   ! Normalize tag2 by vertex count (vpere); loop over elements using same approach as before
00157   ! At the same time, compare values of tag1 and tag2
00158   n_dis = 0
00159   do e = 0, nquads-1
00160      do j = 0, 2
00161         tag2(1+3*e+j) = tag2(1+3*e+j) / vpere;  ! normalize by # verts
00162      end do
00163      if (tag1(1+3*e) .ne. tag2(1+3*e) .or. tag1(1+3*e+1) .ne. tag2(1+3*e+1) .or. tag1(1+3*e+2) .ne. tag2(1+3*i+2)) then
00164         print *, "Tag1, tag2 disagree for element ", e
00165         print *, "tag1: ", tag1(1+3*e), tag1(1+3*e+1), tag1(1+3*e+2)
00166         print *, "tag2: ", tag2(1+3*e), tag2(1+3*e+1), tag2(1+3*e+2)
00167         n_dis = n_dis + 1
00168      end if
00169   end do
00170   if (n_dis .eq. 0) print *, 'All tags agree, success!'
00171 
00172     ! Ok, we are done, shut down MOAB
00173   call iMesh_dtor(%VAL(imesh), ierr)
00174   return
00175 end program DirectAccessNoHolesF90
00176 
00177 subroutine create_mesh_no_holes(imesh, nquads, ierr) 
00178   use ISO_C_BINDING
00179   use freem
00180   implicit none
00181 #include "iMesh_f.h"
00182 
00183   iMesh_Instance imesh
00184   integer nquads, ierr
00185 
00186   real(C_DOUBLE), pointer :: coords(:,:)
00187   TYPE(C_PTR) connect_ptr, tmp_ents_ptr, stat_ptr
00188   iBase_EntityHandle, pointer :: connect(:), tmp_ents(:)
00189   integer, pointer :: stat(:)
00190   integer nverts, tmp_size, stat_size, i
00191 
00192   ! first make the mesh, a 1d array of quads with left hand side x = elem_num; vertices are numbered in layers
00193   ! create verts, num is 2(nquads+1) because they're in a 1d row
00194   nverts = 2*(nquads+1)
00195   allocate(coords(0:2, 0:nverts-1))
00196   do i = 0, nquads-1
00197      coords(0,2*i) = i; coords(0,2*i+1) = i ! x values are all i
00198      coords(1,2*i) = 0.0; coords(1,2*i+1) = 1.0 ! y coords
00199      coords(2,2*i) = 0.0; coords(2,2*i+1) = 0.0 ! z values, all zero (2d mesh)
00200   end do
00201   ! last two vertices
00202   coords(0, nverts-2) = nquads; coords(0, nverts-1) = nquads
00203   coords(1, nverts-2) = 0.0; coords(1, nverts-1) = 1.0 ! y coords
00204   coords(2, nverts-2) = 0.0; coords(2, nverts-1) = 0.0 ! z values, all zero (2d mesh)
00205   tmp_size = 0
00206   call iMesh_createVtxArr(%VAL(imesh), %VAL(nverts), %VAL(iBase_INTERLEAVED), &
00207        coords, %VAL(3*nverts), tmp_ents_ptr, tmp_size, tmp_size, ierr); CHECK(ierr)
00208   call c_f_pointer(tmp_ents_ptr, tmp_ents, [nverts])
00209   deallocate(coords)
00210   allocate(connect(0:4*nquads-1))
00211   do i = 0, nquads-1
00212      connect(4*i+0) = tmp_ents(1+2*i)
00213      connect(4*i+1) = tmp_ents(1+2*i+2)
00214      connect(4*i+2) = tmp_ents(1+2*i+3)
00215      connect(4*i+3) = tmp_ents(1+2*i+1)
00216   end do
00217   stat_size = 0
00218   stat_ptr = C_NULL_PTR
00219   ! re-use tmp_ents here; we know it'll always be larger than nquads so it's ok
00220   call iMesh_createEntArr(%VAL(imesh), %VAL(iMesh_QUADRILATERAL), connect, %VAL(4*nquads), &
00221        tmp_ents_ptr, tmp_size, tmp_size, stat_ptr, stat_size, stat_size, ierr); CHECK(ierr)
00222 
00223   ierr = iBase_SUCCESS
00224 
00225   call c_free(tmp_ents_ptr)
00226   call c_free(stat_ptr)
00227   deallocate(connect)
00228 
00229   return
00230 end subroutine create_mesh_no_holes
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines