 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
1 !> @example DirectAccessNoHolesF90.F90
2 !! \brief Use direct access to MOAB data to avoid calling through API, in Fortran90 \n
3 !!
4 !! This example creates a 1d row of quad elements, such that all quad and vertex handles
5 !! are contiguous in the handle space and in the database. Then it shows how to get access
6 !! to pointers to MOAB-native data for vertex coordinates, quad connectivity, and tag storage
7 !! (vertex to quad adjacency lists aren't accessible from Fortran because they are std::vector's).
8 !! This allows applications to access this data directly
9 !! without going through MOAB's API. In cases where the mesh is not changing (or only mesh
10 !! vertices are moving), this can save significant execution time in applications.
11 !!
12 !! Using direct access (or MOAB in general) from Fortran is complicated in two specific ways:
13 !! 1) There is no way to use arrays with specified dimension starting/ending values with ISO_C_BINDING;
14 !! therefore, all arrays passed back from MOAB/iMesh must use 1-based indices; this makes it a bit
15 !! more difficult to traverse the direct arrays. In this example, I explicitly use indices that
16 !! look like my_array(1+v_ind...) to remind users of this.
17 !! 2) Arithmetic on handles is complicated by the fact that Fortran has no unsigned integer type. I get
18 !! around this by assigning to a C-typed variable first, before the handle arithmetic. This seems to
19 !! work fine. C-typed variables are part of the Fortran95 standard.
20 !!
21 !! ----------------------
22 !! | | | |
23 !! | | | | ...
24 !! | | | |
25 !! ----------------------
26 !!
27 !! -# Initialize MOAB \n
28 !! -# Create a quad mesh, as depicted above
29 !! -# Create 2 dense tags (tag1, tag2) for avg position to assign to quads, and # verts per quad (tag3)
30 !! -# Get connectivity, coordinate, tag1 iterators
31 !! -# Iterate through quads, computing midpoint based on vertex positions, set on quad-based tag1
32 !! -# Iterate through vertices, summing positions into tag2 on connected quads and incrementing vertex count
33 !! -# Iterate through quads, normalizing tag2 by vertex count and comparing values of tag1 and tag2
34 !!
35 !! <b>To compile</b>: \n
36 !! make DirectAccessNoHolesF90 \n
37 !! <b>To run</b>: ./DirectAccessNoHolesF90 [-nquads <# quads>]\n
38 !!
39 !! VSM: Note that IBM xlf compilers do not like dynamic sizing functions for
40 !! C_F_POINTER calls. This leads to internal compiler failure.
41 !! http://www-01.ibm.com/support/docview.wss?uid=swg1IV49428
42 !!
43 #define CHECK(a) if (a .ne. iBase_SUCCESS) call exit(a)
44
45 program directaccessnoholesf90
46 use, intrinsic :: iso_c_binding
47 implicit none
48 #include "iMesh_f.h"
49
50 ! declarations
51 imesh_instance :: imesh
52 ibase_entitysethandle root_set
53 integer :: nverts, tmpflag
54 integer ierr, nquads, nquads_tmp
55 ibase_taghandle tag1h, tag2h
56 TYPE(c_ptr) verts_ptr, quads_ptr, conn_ptr, x_ptr, y_ptr, z_ptr, tag1_ptr, tag2_ptr
57 TYPE(c_ptr) tmp_quads_ptr, offsets_ptr ! pointers that are equivalence'd to arrays using c_f_pointer
58 real(c_double), pointer :: x(:), y(:), z(:), tag1(:), tag2(:) ! arrays equivalence'd to pointers
59 integer, pointer :: offsets(:) ! arrays equivalence'd to pointers
60 ibase_entityhandle, pointer :: quads(:), verts(:), conn(:), tmp_quads(:) ! arrays equivalence'd to pointers
61 ibase_entityarriterator viter, qiter
62 integer(C_SIZE_T) :: startv, startq, v_ind, e_ind ! needed to do handle arithmetic
63 integer count, vpere, e, i, j, v, offsets_size, tmp_quads_size, n_dis
64 character*120 opt
65
66 ! initialize interface and get root set
67 call imesh_newmesh("", imesh, ierr)
68 check(ierr)
69 call imesh_getrootset(%val(imesh), root_set, ierr)
70 check(ierr)
71
72 ! create mesh
73 nquads_tmp = 1000
74 call create_mesh_no_holes(imesh, nquads_tmp, ierr)
75 check(ierr)
76
77 ! get all vertices and quads as arrays
78 nverts = 0
79 call imesh_getentities(%val(imesh), %val(root_set), %val(0), %val(ibase_vertex), &
80 verts_ptr, nverts, nverts, ierr)
81 check(ierr)
82 call c_f_pointer(verts_ptr, verts, [nverts])
83 nquads = 0
84 call imesh_getentities(%val(imesh), %val(root_set), %val(2), %val(imesh_quadrilateral), &
85 quads_ptr, nquads, nquads, ierr)
86 check(ierr)
87 call c_f_pointer(quads_ptr, quads, [nquads])
88
89 ! First vertex/quad is at start of vertex/quad list, and is offset for vertex/quad index calculation
90 startv = verts(1)
91 startq = quads(1)
92 call imesh_freememory(%val(imesh), quads_ptr)
93
94 ! create tag1 (element-based avg), tag2 (vertex-based avg)
95 opt = 'moab:TAG_STORAGE_TYPE=DENSE moab:TAG_DEFAULT_valUE=0'
96 call imesh_createtagwithoptions(%val(imesh), 'tag1', opt, %val(3), %val(ibase_double), &
97 tag1h, ierr, %val(4), %val(56))
98 check(ierr)
99 call imesh_createtagwithoptions(%val(imesh), 'tag2', opt, %val(3), %val(ibase_double), &
100 tag2h, ierr, %val(4), %val(56))
101 check(ierr)
102
103 ! Get iterator to verts, then pointer to coordinates
104 call imesh_initentarriter(%val(imesh), %val(root_set), %val(ibase_vertex), %val(imesh_all_topologies), &
105 %val(nverts), %val(0), viter, ierr)
106 check(ierr)
107 call imesh_coordsiterate(%val(imesh), %val(viter), x_ptr, y_ptr, z_ptr, count, ierr)
108 check(ierr)
109 check(count-nverts) ! check that all verts are in one contiguous chunk
110 call c_f_pointer(x_ptr, x, [nverts])
111 call c_f_pointer(y_ptr, y, [nverts])
112 call c_f_pointer(z_ptr, z, [nverts])
113
114 ! Get iterator to quads, then pointers to connectivity and tags
115 call imesh_initentarriter(%val(imesh), %val(root_set), %val(ibase_face), %val(imesh_all_topologies), &
116 %val(nquads), %val(0), qiter, ierr)
117 check(ierr)
118 call imesh_connectiterate(%val(imesh), %val(qiter), conn_ptr, vpere, count, ierr)
119 check(ierr)
120 call c_f_pointer(conn_ptr, conn, [vpere*nquads])
121 call imesh_tagiterate(%val(imesh), %val(tag1h), %val(qiter), tag1_ptr, count, ierr)
122 check(ierr)
123 call c_f_pointer(tag1_ptr, tag1, [3*nquads])
124 call imesh_tagiterate(%val(imesh), %val(tag2h), %val(qiter), tag2_ptr, count, ierr)
125 check(ierr)
126 call c_f_pointer(tag2_ptr, tag2, [3*nquads])
127
128 ! iterate over elements, computing tag1 from coords positions;
129 ! use (1+... in array indices for 1-based indexing
130 do i = 0, nquads-1
131 tag1(1+3*i+0) = 0.0
132 tag1(1+3*i+1) = 0.0
133 tag1(1+3*i+2) = 0.0 ! initialize position for this element
134 do j = 0, vpere-1 ! loop over vertices in this quad
135 v_ind = conn(1+vpere*i+j) ! assign to v_ind to allow handle arithmetic
136 v_ind = v_ind - startv ! vert index is just the offset from start vertex
137 tag1(1+3*i+0) = tag1(1+3*i+0) + x(1+v_ind)
138 tag1(1+3*i+1) = tag1(1+3*i+1) + y(1+v_ind) ! sum vertex positions into tag1...
139 tag1(1+3*i+2) = tag1(1+3*i+2) + z(1+v_ind)
140 end do
141 tag1(1+3*i+0) = tag1(1+3*i+0) / vpere;
142 tag1(1+3*i+1) = tag1(1+3*i+1) / vpere ! then normalize
143 tag1(1+3*i+2) = tag1(1+3*i+2) / vpere;
144 end do ! loop over elements in chunk
145
146 ! Iterate through vertices, summing positions into tag2 on connected elements;
147 ! get adjacencies all at once,
148 ! could also get them per vertex (time/memory tradeoff)
149 tmp_quads_size = 0
150 offsets_size = 0
151 tmpflag = imesh_quadrilateral
152 call imesh_getentarradj(%val(imesh), verts, %val(nverts), %val(tmpflag), &
153 tmp_quads_ptr, tmp_quads_size, tmp_quads_size, &
154 offsets_ptr, offsets_size, offsets_size, ierr)
155 check(ierr)
156
157 call c_f_pointer(tmp_quads_ptr, tmp_quads, [tmp_quads_size])
158 call c_f_pointer(offsets_ptr, offsets, [offsets_size])
159 call imesh_freememory(%val(imesh), verts_ptr)
160 do v = 0, nverts-1
161 do e = 1+offsets(1+v), 1+offsets(1+v+1)-1 ! 1+ because offsets are in terms of 0-based arrays
162 e_ind = tmp_quads(e) ! assign to e_ind to allow handle arithmetic
163 e_ind = e_ind - startq ! e_ind is the quad handle minus the starting quad handle
164 tag2(1+3*e_ind+0) = tag2(1+3*e_ind+0) + x(1+v)
165 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
166 tag2(1+3*e_ind+2) = tag2(1+3*e_ind+2) + z(1+v)
167 end do
168 end do
169 call imesh_freememory(%val(imesh), tmp_quads_ptr)
170 call imesh_freememory(%val(imesh), offsets_ptr)
171
172 ! Normalize tag2 by vertex count (vpere);
173 ! loop over elements using same approach as before
174 ! At the same time, compare values of tag1 and tag2
175 n_dis = 0
176 do e = 0, nquads-1
177 do j = 0, 2
178 tag2(1+3*e+j) = tag2(1+3*e+j) / vpere ! normalize by # verts
179 end do
180 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
181 print *, "Tag1, tag2 disagree for element ", e
182 print *, "tag1: ", tag1(1+3*e), tag1(1+3*e+1), tag1(1+3*e+2)
183 print *, "tag2: ", tag2(1+3*e), tag2(1+3*e+1), tag2(1+3*e+2)
184 n_dis = n_dis + 1
185 endif
186 end do
187 if (n_dis .eq. 0) then
188 print *, 'All tags agree, success!'
189 endif
190
191 ! Ok, we are done, shut down MOAB
192 call imesh_dtor(%val(imesh), ierr)
193 check(ierr)
194 end program directaccessnoholesf90
195
196 subroutine create_mesh_no_holes(imesh, nquads, ierr)
197 use, intrinsic :: iso_c_binding
198 implicit none
199 #include "iMesh_f.h"
200
201 imesh_instance imesh
202 integer nquads, ierr
203
204 real(C_DOUBLE), pointer :: coords(:,:)
205 TYPE(c_ptr) tmp_ents_ptr, stat_ptr
206 ibase_entityhandle, pointer :: connect(:), tmp_ents(:)
207 integer nverts, tmp_size, stat_size, i
208
209 ! first make the mesh, a 1d array of quads with left hand side x = elem_num;
210 ! vertices are numbered in layers
211 ! create verts, num is 2(nquads+1) because they're in a 1d row
212 nverts = 2*(nquads+1)
213 allocate(coords(0:2, 0:nverts-1))
214 do i = 0, nquads-1
215 coords(0,2*i) = i
216 coords(0,2*i+1) = i ! x values are all i
217 coords(1,2*i) = 0.0
218 coords(1,2*i+1) = 1.0 ! y coords
219 coords(2,2*i) = 0.0
220 coords(2,2*i+1) = 0.0 ! z values, all zero (2d mesh)
221 end do
222 ! last two vertices
223 coords(0, nverts-2) = nquads
224 coords(0, nverts-1) = nquads
225 coords(1, nverts-2) = 0.0
226 coords(1, nverts-1) = 1.0 ! y coords
227 coords(2, nverts-2) = 0.0
228 coords(2, nverts-1) = 0.0 ! z values, all zero (2d mesh)
229 tmp_size = 0
230 !!!!!! VSM: This is the culprit call that screws up IBM xlf compiler
231 call imesh_createvtxarr(%val(imesh), %val(nverts), %val(ibase_interleaved), &
232 coords, %val(3*nverts), tmp_ents_ptr, tmp_size, tmp_size, ierr)
233 check(ierr)
234 call c_f_pointer(tmp_ents_ptr, tmp_ents, [nverts])
235 deallocate(coords)
236 allocate(connect(0:4*nquads-1))
237 do i = 0, nquads-1
238 connect(4*i+0) = tmp_ents(1+2*i)
239 connect(4*i+1) = tmp_ents(1+2*i+2)
240 connect(4*i+2) = tmp_ents(1+2*i+3)
241 connect(4*i+3) = tmp_ents(1+2*i+1)
242 end do
243 stat_size = 0
244 stat_ptr = c_null_ptr
245 ! re-use tmp_ents here;
246 ! we know it'll always be larger than nquads so it's ok
247 call imesh_createentarr(%val(imesh), %val(imesh_quadrilateral), connect, %val(4*nquads), &
248 tmp_ents_ptr, tmp_size, tmp_size, stat_ptr, stat_size, stat_size, ierr)
249 check(ierr)
250
251 ierr = ibase_success
252
253 call imesh_freememory(%val(imesh), tmp_ents_ptr)
254 call imesh_freememory(%val(imesh), stat_ptr)
255 deallocate(connect)
256
257 return
258 end subroutine create_mesh_no_holes
259
1.9.1.