MOAB: Mesh Oriented datABase  (version 5.4.1)
DirectAccessNoHolesF90.F90
Go to the documentation of this file.
00001 !> @example DirectAccessNoHolesF90.F90
00002 !! \brief Use direct access to MOAB data to avoid calling through API, in Fortran90 \n
00003 !!
00004 !! This example creates a 1d row of quad elements, such that all quad and vertex handles
00005 !! are contiguous in the handle space and in the database.  Then it shows how to get access
00006 !! to pointers to MOAB-native data for vertex coordinates, quad connectivity, and tag storage
00007 !! (vertex to quad adjacency lists aren't accessible from Fortran because they are std::vector's).
00008 !! This allows applications to access this data directly
00009 !! without going through MOAB's API.  In cases where the mesh is not changing (or only mesh
00010 !! vertices are moving), this can save significant execution time in applications.
00011 !!
00012 !! Using direct access (or MOAB in general) from Fortran is complicated in two specific ways:
00013 !! 1) There is no way to use arrays with specified dimension starting/ending values with ISO_C_BINDING;
00014 !!    therefore, all arrays passed back from MOAB/iMesh must use 1-based indices; this makes it a bit
00015 !!    more difficult to traverse the direct arrays.  In this example, I explicitly use indices that
00016 !!    look like my_array(1+v_ind...) to remind users of this.
00017 !! 2) Arithmetic on handles is complicated by the fact that Fortran has no unsigned integer type.  I get
00018 !!    around this by assigning to a C-typed variable first, before the handle arithmetic.  This seems to
00019 !!    work fine.  C-typed variables are part of the Fortran95 standard.
00020 !!
00021 !!  ----------------------
00022 !!  |      |      |      |
00023 !!  |      |      |      | ...
00024 !!  |      |      |      |
00025 !!  ----------------------
00026 !!
00027 !!    -#  Initialize MOAB \n
00028 !!    -#  Create a quad mesh, as depicted above
00029 !!    -#  Create 2 dense tags (tag1, tag2) for avg position to assign to quads, and # verts per quad (tag3)
00030 !!    -#  Get connectivity, coordinate, tag1 iterators
00031 !!    -#  Iterate through quads, computing midpoint based on vertex positions, set on quad-based tag1
00032 !!    -#  Iterate through vertices, summing positions into tag2 on connected quads and incrementing vertex count
00033 !!    -#  Iterate through quads, normalizing tag2 by vertex count and comparing values of tag1 and tag2
00034 !!
00035 !! <b>To compile</b>: \n
00036 !!    make DirectAccessNoHolesF90 MOAB_DIR=<installdir>  \n
00037 !! <b>To run</b>: ./DirectAccessNoHolesF90 [-nquads <# quads>]\n
00038 !!
00039 !! VSM: Note that IBM xlf compilers do not like dynamic sizing functions for
00040 !!      C_F_POINTER calls. This leads to internal compiler failure.
00041 !!      http://www-01.ibm.com/support/docview.wss?uid=swg1IV49428
00042 !!
00043 #define CHECK(a)   if (a .ne. iBase_SUCCESS) call exit(a)
00044 
00045 program DirectAccessNoHolesF90
00046   use, intrinsic :: ISO_C_BINDING
00047   implicit none
00048 #include "iMesh_f.h"
00049 
00050   ! declarations
00051   iMesh_Instance :: imesh
00052   iBase_EntitySetHandle root_set
00053   integer :: nverts, tmpflag
00054   integer ierr, nquads, nquads_tmp
00055   iBase_TagHandle tag1h, tag2h
00056   TYPE(C_PTR) verts_ptr, quads_ptr, conn_ptr, x_ptr, y_ptr, z_ptr, tag1_ptr, tag2_ptr
00057   TYPE(C_PTR) tmp_quads_ptr, offsets_ptr  ! pointers that are equivalence'd to arrays using c_f_pointer
00058   real(C_DOUBLE), pointer :: x(:), y(:), z(:), tag1(:), tag2(:) ! arrays equivalence'd to pointers
00059   integer, pointer :: offsets(:)                       ! arrays equivalence'd to pointers
00060   iBase_EntityHandle, pointer :: quads(:), verts(:), conn(:), tmp_quads(:)  ! arrays equivalence'd to pointers
00061   iBase_EntityArrIterator viter, qiter
00062   integer(C_SIZE_T) :: startv, startq, v_ind, e_ind  ! needed to do handle arithmetic
00063   integer count, vpere, e, i, j, v, offsets_size, tmp_quads_size, n_dis
00064   character*120 opt
00065 
00066   ! initialize interface and get root set
00067   call iMesh_newMesh("", imesh, ierr)
00068   CHECK(ierr)
00069   call iMesh_getRootSet(%val(imesh), root_set, ierr)
00070   CHECK(ierr)
00071 
00072   ! create mesh
00073   nquads_tmp = 1000
00074   call create_mesh_no_holes(imesh, nquads_tmp, ierr)
00075   CHECK(ierr)
00076 
00077   ! get all vertices and quads as arrays
00078   nverts = 0
00079   call iMesh_getEntities(%val(imesh), %val(root_set), %val(0), %val(iBase_VERTEX), &
00080        verts_ptr, nverts, nverts, ierr)
00081   CHECK(ierr)
00082   call c_f_pointer(verts_ptr, verts, [nverts])
00083   nquads = 0
00084   call iMesh_getEntities(%val(imesh), %val(root_set), %val(2), %val(iMesh_QUADRILATERAL), &
00085        quads_ptr, nquads, nquads, ierr)
00086   CHECK(ierr)
00087   call c_f_pointer(quads_ptr, quads, [nquads])
00088 
00089   ! First vertex/quad is at start of vertex/quad list, and is offset for vertex/quad index calculation
00090   startv = verts(1)
00091   startq = quads(1)
00092   call iMesh_freeMemory(%val(imesh), quads_ptr)
00093 
00094   ! create tag1 (element-based avg), tag2 (vertex-based avg)
00095   opt = 'moab:TAG_STORAGE_TYPE=DENSE moab:TAG_DEFAULT_valUE=0'
00096   call iMesh_createTagWithOptions(%val(imesh), 'tag1', opt, %val(3), %val(iBase_DOUBLE), &
00097        tag1h, ierr, %val(4), %val(56))
00098   CHECK(ierr)
00099   call iMesh_createTagWithOptions(%val(imesh), 'tag2', opt, %val(3), %val(iBase_DOUBLE), &
00100        tag2h, ierr, %val(4), %val(56))
00101   CHECK(ierr)
00102 
00103   ! Get iterator to verts, then pointer to coordinates
00104   call iMesh_initEntArrIter(%val(imesh), %val(root_set), %val(iBase_VERTEX), %val(iMesh_ALL_TOPOLOGIES), &
00105        %val(nverts), %val(0), viter, ierr)
00106   CHECK(ierr)
00107   call iMesh_coordsIterate(%val(imesh), %val(viter), x_ptr, y_ptr, z_ptr, count, ierr)
00108   CHECK(ierr)
00109   CHECK(count-nverts) ! check that all verts are in one contiguous chunk
00110   call c_f_pointer(x_ptr, x, [nverts])
00111   call c_f_pointer(y_ptr, y, [nverts])
00112   call c_f_pointer(z_ptr, z, [nverts])
00113 
00114   ! Get iterator to quads, then pointers to connectivity and tags
00115   call iMesh_initEntArrIter(%val(imesh), %val(root_set), %val(iBase_FACE), %val(iMesh_ALL_TOPOLOGIES), &
00116        %val(nquads), %val(0), qiter, ierr)
00117   CHECK(ierr)
00118   call iMesh_connectIterate(%val(imesh), %val(qiter), conn_ptr, vpere, count, ierr)
00119   CHECK(ierr)
00120   call c_f_pointer(conn_ptr, conn, [vpere*nquads])
00121   call iMesh_tagIterate(%val(imesh), %val(tag1h), %val(qiter), tag1_ptr, count, ierr)
00122   CHECK(ierr)
00123   call c_f_pointer(tag1_ptr, tag1, [3*nquads])
00124   call iMesh_tagIterate(%val(imesh), %val(tag2h), %val(qiter), tag2_ptr, count, ierr)
00125   CHECK(ierr)
00126   call c_f_pointer(tag2_ptr, tag2, [3*nquads])
00127   
00128   ! iterate over elements, computing tag1 from coords positions;
00129   ! use (1+... in array indices for 1-based indexing
00130   do i = 0, nquads-1
00131      tag1(1+3*i+0) = 0.0
00132      tag1(1+3*i+1) = 0.0
00133      tag1(1+3*i+2) = 0.0 ! initialize position for this element
00134      do j = 0, vpere-1 ! loop over vertices in this quad
00135         v_ind = conn(1+vpere*i+j) ! assign to v_ind to allow handle arithmetic
00136         v_ind = v_ind - startv ! vert index is just the offset from start vertex
00137         tag1(1+3*i+0) = tag1(1+3*i+0) + x(1+v_ind)
00138         tag1(1+3*i+1) = tag1(1+3*i+1) + y(1+v_ind) ! sum vertex positions into tag1...
00139         tag1(1+3*i+2) = tag1(1+3*i+2) + z(1+v_ind)
00140      end do
00141      tag1(1+3*i+0) = tag1(1+3*i+0) / vpere;
00142      tag1(1+3*i+1) = tag1(1+3*i+1) / vpere  ! then normalize
00143      tag1(1+3*i+2) = tag1(1+3*i+2) / vpere;
00144   end do ! loop over elements in chunk
00145     
00146   ! Iterate through vertices, summing positions into tag2 on connected elements;
00147   ! get adjacencies all at once,
00148   ! could also get them per vertex (time/memory tradeoff)
00149   tmp_quads_size = 0
00150   offsets_size = 0
00151   tmpflag = iMesh_QUADRILATERAL
00152   call iMesh_getEntArrAdj(%val(imesh), verts, %val(nverts), %val(tmpflag), &
00153        tmp_quads_ptr, tmp_quads_size, tmp_quads_size, &
00154        offsets_ptr, offsets_size, offsets_size, ierr)
00155   CHECK(ierr)
00156 
00157   call c_f_pointer(tmp_quads_ptr, tmp_quads, [tmp_quads_size])
00158   call c_f_pointer(offsets_ptr, offsets, [offsets_size])
00159   call iMesh_freeMemory(%val(imesh), verts_ptr)
00160   do v = 0, nverts-1
00161      do e = 1+offsets(1+v), 1+offsets(1+v+1)-1  ! 1+ because offsets are in terms of 0-based arrays
00162         e_ind = tmp_quads(e) ! assign to e_ind to allow handle arithmetic
00163         e_ind = e_ind - startq ! e_ind is the quad handle minus the starting quad handle
00164         tag2(1+3*e_ind+0) = tag2(1+3*e_ind+0) + x(1+v)
00165         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
00166         tag2(1+3*e_ind+2) = tag2(1+3*e_ind+2) + z(1+v)
00167      end do
00168   end do
00169   call iMesh_freeMemory(%val(imesh), tmp_quads_ptr)
00170   call iMesh_freeMemory(%val(imesh), offsets_ptr)
00171 
00172   ! Normalize tag2 by vertex count (vpere); 
00173   ! loop over elements using same approach as before
00174   ! At the same time, compare values of tag1 and tag2
00175   n_dis = 0
00176   do e = 0, nquads-1
00177      do j = 0, 2
00178         tag2(1+3*e+j) = tag2(1+3*e+j) / vpere  ! normalize by # verts
00179      end do
00180      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*e+2)) then
00181         print *, "Tag1, tag2 disagree for element ", e
00182         print *, "tag1: ", tag1(1+3*e), tag1(1+3*e+1), tag1(1+3*e+2)
00183         print *, "tag2: ", tag2(1+3*e), tag2(1+3*e+1), tag2(1+3*e+2)
00184         n_dis = n_dis + 1
00185      endif
00186   end do
00187   if (n_dis .eq. 0) then
00188     print *, 'All tags agree, success!'
00189   endif
00190 
00191     ! Ok, we are done, shut down MOAB
00192   call iMesh_dtor(%val(imesh), ierr)
00193   CHECK(ierr)
00194 end program DirectAccessNoHolesF90
00195 
00196 subroutine create_mesh_no_holes(imesh, nquads, ierr) 
00197   use, intrinsic :: ISO_C_BINDING
00198   implicit none
00199 #include "iMesh_f.h"
00200 
00201   iMesh_Instance imesh
00202   integer nquads, ierr
00203 
00204   real(C_DOUBLE), pointer :: coords(:,:)
00205   TYPE(C_PTR) tmp_ents_ptr, stat_ptr
00206   iBase_EntityHandle, pointer :: connect(:), tmp_ents(:)
00207   integer nverts, tmp_size, stat_size, i
00208 
00209   ! first make the mesh, a 1d array of quads with left hand side x = elem_num;
00210   ! vertices are numbered in layers
00211   ! create verts, num is 2(nquads+1) because they're in a 1d row
00212   nverts = 2*(nquads+1)
00213   allocate(coords(0:2, 0:nverts-1))
00214   do i = 0, nquads-1
00215      coords(0,2*i) = i
00216      coords(0,2*i+1) = i ! x values are all i
00217      coords(1,2*i) = 0.0
00218      coords(1,2*i+1) = 1.0 ! y coords
00219      coords(2,2*i) = 0.0
00220      coords(2,2*i+1) = 0.0 ! z values, all zero (2d mesh)
00221   end do
00222   ! last two vertices
00223   coords(0, nverts-2) = nquads
00224   coords(0, nverts-1) = nquads
00225   coords(1, nverts-2) = 0.0
00226   coords(1, nverts-1) = 1.0 ! y coords
00227   coords(2, nverts-2) = 0.0
00228   coords(2, nverts-1) = 0.0 ! z values, all zero (2d mesh)
00229   tmp_size = 0
00230   !!!!!! VSM: This is the culprit call that screws up IBM xlf compiler
00231   call iMesh_createVtxArr(%val(imesh), %val(nverts), %val(iBase_INTERLEAVED), &
00232        coords, %val(3*nverts), tmp_ents_ptr, tmp_size, tmp_size, ierr)
00233   CHECK(ierr)
00234   call c_f_pointer(tmp_ents_ptr, tmp_ents, [nverts])
00235   deallocate(coords)
00236   allocate(connect(0:4*nquads-1))
00237   do i = 0, nquads-1
00238      connect(4*i+0) = tmp_ents(1+2*i)
00239      connect(4*i+1) = tmp_ents(1+2*i+2)
00240      connect(4*i+2) = tmp_ents(1+2*i+3)
00241      connect(4*i+3) = tmp_ents(1+2*i+1)
00242   end do
00243   stat_size = 0
00244   stat_ptr = C_NULL_PTR
00245   ! re-use tmp_ents here;
00246   ! we know it'll always be larger than nquads so it's ok
00247   call iMesh_createEntArr(%val(imesh), %val(iMesh_QUADRILATERAL), connect, %val(4*nquads), &
00248        tmp_ents_ptr, tmp_size, tmp_size, stat_ptr, stat_size, stat_size, ierr)
00249   CHECK(ierr)
00250 
00251   ierr = iBase_SUCCESS
00252 
00253   call iMesh_freeMemory(%val(imesh), tmp_ents_ptr)
00254   call iMesh_freeMemory(%val(imesh), stat_ptr)
00255   deallocate(connect)
00256 
00257   return
00258 end subroutine create_mesh_no_holes
00259 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines