moab
DirectAccessNoHolesF90.F90
Go to the documentation of this file.
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