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