moab
|
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.
---------------------- | | | | | | | | ... | | | | ----------------------
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