 |
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 #include "SphereDecomp.hpp"
2 #include "moab/MeshTopoUtil.hpp"
3 #include "moab/Range.hpp"
4 #include "moab/CN.hpp"
5 #include <cmath>
6 #include <cassert>
7 #include <iostream>
8
9 #define RR \
10 if( MB_SUCCESS != result ) return result
11
12 const char* SUBDIV_VERTICES_TAG_NAME = "subdiv_vertices";
13
14 using namespace moab;
15
16 SphereDecomp::SphereDecomp( Interface* impl )
17 {
18 mbImpl = impl;
19 }
20
21 ErrorCode SphereDecomp::build_sphere_mesh( const char* sphere_radii_tag_name, EntityHandle* hex_set )
22 {
23 ErrorCode result = mbImpl->tag_get_handle( sphere_radii_tag_name, 1, MB_TYPE_DOUBLE, sphereRadiiTag );RR;
24
25 // need to make sure all interior edges and faces are created
26 Range all_verts;
27 result = mbImpl->get_entities_by_type( 0, MBVERTEX, all_verts );RR;
28 MeshTopoUtil mtu( mbImpl );
29 result = mtu.construct_aentities( all_verts );RR;
30
31 // create tag to hold vertices
32 result = mbImpl->tag_get_handle( SUBDIV_VERTICES_TAG_NAME, 9, MB_TYPE_HANDLE, subdivVerticesTag,
33 MB_TAG_DENSE | MB_TAG_EXCL );RR;
34
35 // compute nodal positions for each dimension element
36 result = compute_nodes( 1 );RR;
37 result = compute_nodes( 2 );RR;
38 result = compute_nodes( 3 );RR;
39
40 // build hex elements
41 std::vector< EntityHandle > sphere_hexes, interstic_hexes;
42 result = build_hexes( sphere_hexes, interstic_hexes );RR;
43
44 result = mbImpl->tag_delete( subdivVerticesTag );RR;
45
46 if( NULL != hex_set )
47 {
48 if( 0 == *hex_set )
49 {
50 EntityHandle this_set;
51 // make a new set
52 result = mbImpl->create_meshset( MESHSET_SET, this_set );RR;
53 *hex_set = this_set;
54 }
55
56 // save all the hexes to this set
57 result = mbImpl->add_entities( *hex_set, &sphere_hexes[0], sphere_hexes.size() );RR;
58 result = mbImpl->add_entities( *hex_set, &interstic_hexes[0], interstic_hexes.size() );RR;
59 }
60
61 return result;
62 }
63
64 ErrorCode SphereDecomp::compute_nodes( const int dim )
65 {
66 // get facets of that dimension
67 Range these_ents;
68 const EntityType the_types[4] = { MBVERTEX, MBEDGE, MBTRI, MBTET };
69
70 ErrorCode result = mbImpl->get_entities_by_dimension( 0, dim, these_ents );RR;
71 assert( mbImpl->type_from_handle( *these_ents.begin() ) == the_types[dim] &&
72 mbImpl->type_from_handle( *these_ents.rbegin() ) == the_types[dim] );
73
74 EntityHandle subdiv_vertices[9];
75 MeshTopoUtil mtu( mbImpl );
76 double avg_pos[3], vert_pos[12], new_vert_pos[12], new_new_vert_pos[3];
77 double radii[4], unitv[3];
78 int num_verts = CN::VerticesPerEntity( the_types[dim] );
79
80 for( Range::iterator rit = these_ents.begin(); rit != these_ents.end(); ++rit )
81 {
82
83 // get vertices
84 const EntityHandle* connect;
85 int num_connect;
86 result = mbImpl->get_connectivity( *rit, connect, num_connect );RR;
87
88 // compute center
89 result = mtu.get_average_position( connect, num_connect, avg_pos );RR;
90
91 // create center vertex
92 result = mbImpl->create_vertex( avg_pos, subdiv_vertices[num_verts] );RR;
93
94 // get coords of other vertices
95 result = mbImpl->get_coords( connect, num_connect, vert_pos );RR;
96
97 // get radii associated with each vertex
98 result = mbImpl->tag_get_data( sphereRadiiTag, connect, num_connect, radii );RR;
99
100 // compute subdiv vertex position for each vertex
101 for( int i = 0; i < num_verts; i++ )
102 {
103 for( int j = 0; j < 3; j++ )
104 unitv[j] = avg_pos[j] - vert_pos[3 * i + j];
105 double vlength = sqrt( unitv[0] * unitv[0] + unitv[1] * unitv[1] + unitv[2] * unitv[2] );
106 if( vlength < radii[i] )
107 {
108 std::cout << "Radius too large at vertex " << i << std::endl;
109 result = MB_FAILURE;
110 continue;
111 }
112
113 for( int j = 0; j < 3; j++ )
114 unitv[j] /= vlength;
115
116 for( int j = 0; j < 3; j++ )
117 new_vert_pos[3 * i + j] = vert_pos[3 * i + j] + radii[i] * unitv[j];
118
119 // create vertex at this position
120 ErrorCode tmp_result = mbImpl->create_vertex( &new_vert_pos[3 * i], subdiv_vertices[i] );
121 if( MB_SUCCESS != tmp_result ) result = tmp_result;
122 }
123
124 if( MB_SUCCESS != result ) return result;
125
126 // compute subdiv vertex positions for vertices inside spheres; just mid-pt between
127 // previous subdiv vertex and corner vertex
128 for( int i = 0; i < num_verts; i++ )
129 {
130 for( int j = 0; j < 3; j++ )
131 new_new_vert_pos[j] = .5 * ( vert_pos[3 * i + j] + new_vert_pos[3 * i + j] );
132
133 result = mbImpl->create_vertex( new_new_vert_pos, subdiv_vertices[num_verts + 1 + i] );
134 }
135
136 // set the tag
137 result = mbImpl->tag_set_data( subdivVerticesTag, &( *rit ), 1, subdiv_vertices );RR;
138 }
139
140 return result;
141 }
142
143 ErrorCode SphereDecomp::build_hexes( std::vector< EntityHandle >& sphere_hexes,
144 std::vector< EntityHandle >& interstic_hexes )
145 {
146 // build hexes inside each tet element separately
147 Range tets;
148 ErrorCode result = mbImpl->get_entities_by_type( 0, MBTET, tets );RR;
149
150 for( Range::iterator vit = tets.begin(); vit != tets.end(); ++vit )
151 {
152 result = subdivide_tet( *vit, sphere_hexes, interstic_hexes );RR;
153 }
154
155 return MB_SUCCESS;
156 }
157
158 ErrorCode SphereDecomp::subdivide_tet( EntityHandle tet,
159 std::vector< EntityHandle >& sphere_hexes,
160 std::vector< EntityHandle >& interstic_hexes )
161 {
162 // 99: (#subdiv_verts/entity=9) * (#edges=6 + #faces=4 + 1=tet)
163 EntityHandle subdiv_verts[99];
164
165 // get tet connectivity
166 std::vector< EntityHandle > tet_conn;
167 ErrorCode result = mbImpl->get_connectivity( &tet, 1, tet_conn );RR;
168
169 for( int dim = 1; dim <= 3; dim++ )
170 {
171 // get entities of this dimension
172 std::vector< EntityHandle > ents;
173 if( dim != 3 )
174 {
175 result = mbImpl->get_adjacencies( &tet, 1, dim, false, ents );RR;
176 }
177 else
178 ents.push_back( tet );
179
180 // for each, get subdiv verts & put into vector
181 for( std::vector< EntityHandle >::iterator vit = ents.begin(); vit != ents.end(); ++vit )
182 {
183 result = retrieve_subdiv_verts( tet, *vit, &tet_conn[0], dim, subdiv_verts );RR;
184 }
185 }
186
187 // ok, subdiv_verts are in canonical order; now create the hexes, using pre-computed templates
188
189 // Templates are specified in terms of the vertices making up each hex; vertices are specified
190 // by specifying the facet index and type they resolve, and the index of that vertex in that
191 // facet's subdivision vertices list.
192
193 // Each facet is subdivided into:
194 // - a mid vertex
195 // - one vertex for each corner vertex on the facet (located on a line between the mid vertex
196 // and
197 // the corresponding corner vertex, a distance equal to the sphere radius away from the corner
198 // vertex)
199 // - one vertex midway between each corner vertex and the corresponding "sphere surface" vertex
200 // For edges, tris and tets this gives 5, 7 and 9 subdivision vertices, respectively.
201 // Subdivision vertices appear in the list in the order: sphere surface vertices, mid vertex,
202 // sphere interior vertices. In each of those sub lists, vertices are listed in the canonical
203 // order of the corresponding corner vertices for that facet.
204
205 // Subdivision vertices for facetes are indexed by listing the facet type they resolve (EDGE,
206 // FACE, TET), the index of that facet (integer = 0..5, 0..3, 0 for edges, tris, tet, resp), and
207 // subdivision index (AINDEX..EINDEX for edges, AINDEX..GINDEX for tris, AINDEX..IINDEX for
208 // tets).
209
210 // Subdivision vertices for all facets of a tet are stored in one subdivision vertex vector, in
211 // order of increasing facet dimension and index (index varies fastest). The ESV, FSV, and TSV
212 // macros are used to compute the indices into that vector for various parameters. The CV macro
213 // is used to index into the tet connectivity vector.
214
215 // Subdivision templates for splitting the tet into 28 hexes were derived by hand, and are
216 // listed below (using the indexing scheme described above).
217
218 #define EDGE 0
219 #define FACE 1
220 #define TET 2
221 #define AINDEX 0
222 #define BINDEX 1
223 #define CINDEX 2
224 #define DINDEX 3
225 #define EINDEX 4
226 #define FINDEX 5
227 #define GINDEX 6
228 #define HINDEX 7
229 #define IINDEX 8
230 #define V0INDEX 0
231 #define V1INDEX 1
232 #define V2INDEX 2
233 #define V3INDEX 3
234 #define CV( a ) tet_conn[a]
235 #define ESV( a, b ) subdiv_verts[(a)*9 + ( b )]
236 #define FSV( a, b ) subdiv_verts[54 + (a)*9 + ( b )]
237 #define TSV( a, b ) subdiv_verts[90 + (a)*9 + ( b )]
238
239 EntityHandle this_connect[8], this_hex;
240
241 // first, interstices hexes, three per vertex/spherical surface
242 // V0:
243 int i = 0;
244 this_connect[i++] = ESV( 0, AINDEX );
245 this_connect[i++] = ESV( 0, CINDEX );
246 this_connect[i++] = FSV( 3, DINDEX );
247 this_connect[i++] = FSV( 3, AINDEX );
248 this_connect[i++] = FSV( 0, AINDEX );
249 this_connect[i++] = FSV( 0, DINDEX );
250 this_connect[i++] = TSV( 0, EINDEX );
251 this_connect[i++] = TSV( 0, AINDEX );
252 result = mbImpl->create_element( MBHEX, this_connect, 8, this_hex );RR;
253 interstic_hexes.push_back( this_hex );
254
255 i = 0;
256 this_connect[i++] = FSV( 0, AINDEX );
257 this_connect[i++] = FSV( 0, DINDEX );
258 this_connect[i++] = TSV( 0, EINDEX );
259 this_connect[i++] = TSV( 0, AINDEX );
260 this_connect[i++] = ESV( 3, AINDEX );
261 this_connect[i++] = ESV( 3, CINDEX );
262 this_connect[i++] = FSV( 2, DINDEX );
263 this_connect[i++] = FSV( 2, AINDEX );
264 result = mbImpl->create_element( MBHEX, this_connect, 8, this_hex );RR;
265 interstic_hexes.push_back( this_hex );
266
267 i = 0;
268 this_connect[i++] = FSV( 3, AINDEX );
269 this_connect[i++] = FSV( 3, DINDEX );
270 this_connect[i++] = ESV( 2, CINDEX );
271 this_connect[i++] = ESV( 2, BINDEX );
272 this_connect[i++] = TSV( 0, AINDEX );
273 this_connect[i++] = TSV( 0, EINDEX );
274 this_connect[i++] = FSV( 2, DINDEX );
275 this_connect[i++] = FSV( 2, AINDEX );
276 result = mbImpl->create_element( MBHEX, this_connect, 8, this_hex );RR;
277 interstic_hexes.push_back( this_hex );
278
279 // V1:
280 i = 0;
281 this_connect[i++] = ESV( 0, CINDEX );
282 this_connect[i++] = ESV( 0, BINDEX );
283 this_connect[i++] = FSV( 3, CINDEX );
284 this_connect[i++] = FSV( 3, DINDEX );
285 this_connect[i++] = FSV( 0, DINDEX );
286 this_connect[i++] = FSV( 0, BINDEX );
287 this_connect[i++] = TSV( 0, BINDEX );
288 this_connect[i++] = TSV( 0, EINDEX );
289 result = mbImpl->create_element( MBHEX, this_connect, 8, this_hex );RR;
290 interstic_hexes.push_back( this_hex );
291
292 i = 0;
293 this_connect[i++] = FSV( 0, DINDEX );
294 this_connect[i++] = FSV( 0, BINDEX );
295 this_connect[i++] = TSV( 0, BINDEX );
296 this_connect[i++] = TSV( 0, EINDEX );
297 this_connect[i++] = ESV( 4, CINDEX );
298 this_connect[i++] = ESV( 4, AINDEX );
299 this_connect[i++] = FSV( 1, AINDEX );
300 this_connect[i++] = FSV( 1, DINDEX );
301 result = mbImpl->create_element( MBHEX, this_connect, 8, this_hex );RR;
302 interstic_hexes.push_back( this_hex );
303
304 i = 0;
305 this_connect[i++] = FSV( 1, DINDEX );
306 this_connect[i++] = FSV( 1, AINDEX );
307 this_connect[i++] = TSV( 0, BINDEX );
308 this_connect[i++] = TSV( 0, EINDEX );
309 this_connect[i++] = ESV( 1, CINDEX );
310 this_connect[i++] = ESV( 1, AINDEX );
311 this_connect[i++] = FSV( 3, CINDEX );
312 this_connect[i++] = FSV( 3, DINDEX );
313 result = mbImpl->create_element( MBHEX, this_connect, 8, this_hex );RR;
314 interstic_hexes.push_back( this_hex );
315
316 // V2:
317 i = 0;
318 this_connect[i++] = FSV( 3, DINDEX );
319 this_connect[i++] = ESV( 1, CINDEX );
320 this_connect[i++] = ESV( 1, BINDEX );
321 this_connect[i++] = FSV( 3, BINDEX );
322 this_connect[i++] = TSV( 0, EINDEX );
323 this_connect[i++] = FSV( 1, DINDEX );
324 this_connect[i++] = FSV( 1, BINDEX );
325 this_connect[i++] = TSV( 0, CINDEX );
326 result = mbImpl->create_element( MBHEX, this_connect, 8, this_hex );RR;
327 interstic_hexes.push_back( this_hex );
328
329 i = 0;
330 this_connect[i++] = TSV( 0, EINDEX );
331 this_connect[i++] = FSV( 1, DINDEX );
332 this_connect[i++] = FSV( 1, BINDEX );
333 this_connect[i++] = TSV( 0, CINDEX );
334 this_connect[i++] = FSV( 2, DINDEX );
335 this_connect[i++] = ESV( 5, CINDEX );
336 this_connect[i++] = ESV( 5, AINDEX );
337 this_connect[i++] = FSV( 2, CINDEX );
338 result = mbImpl->create_element( MBHEX, this_connect, 8, this_hex );RR;
339 interstic_hexes.push_back( this_hex );
340
341 i = 0;
342 this_connect[i++] = TSV( 0, CINDEX );
343 this_connect[i++] = FSV( 2, CINDEX );
344 this_connect[i++] = ESV( 2, AINDEX );
345 this_connect[i++] = FSV( 3, BINDEX );
346 this_connect[i++] = TSV( 0, EINDEX );
347 this_connect[i++] = FSV( 2, DINDEX );
348 this_connect[i++] = ESV( 2, CINDEX );
349 this_connect[i++] = FSV( 3, DINDEX );
350 result = mbImpl->create_element( MBHEX, this_connect, 8, this_hex );RR;
351 interstic_hexes.push_back( this_hex );
352
353 // V3:
354 i = 0;
355 this_connect[i++] = TSV( 0, EINDEX );
356 this_connect[i++] = FSV( 1, DINDEX );
357 this_connect[i++] = ESV( 5, CINDEX );
358 this_connect[i++] = FSV( 2, DINDEX );
359 this_connect[i++] = TSV( 0, DINDEX );
360 this_connect[i++] = FSV( 1, CINDEX );
361 this_connect[i++] = ESV( 5, BINDEX );
362 this_connect[i++] = FSV( 2, BINDEX );
363 result = mbImpl->create_element( MBHEX, this_connect, 8, this_hex );RR;
364 interstic_hexes.push_back( this_hex );
365
366 i = 0;
367 this_connect[i++] = FSV( 0, DINDEX );
368 this_connect[i++] = ESV( 4, CINDEX );
369 this_connect[i++] = FSV( 1, DINDEX );
370 this_connect[i++] = TSV( 0, EINDEX );
371 this_connect[i++] = FSV( 0, CINDEX );
372 this_connect[i++] = ESV( 4, BINDEX );
373 this_connect[i++] = FSV( 1, CINDEX );
374 this_connect[i++] = TSV( 0, DINDEX );
375 result = mbImpl->create_element( MBHEX, this_connect, 8, this_hex );RR;
376 interstic_hexes.push_back( this_hex );
377
378 i = 0;
379 this_connect[i++] = ESV( 3, CINDEX );
380 this_connect[i++] = FSV( 0, DINDEX );
381 this_connect[i++] = TSV( 0, EINDEX );
382 this_connect[i++] = FSV( 2, DINDEX );
383 this_connect[i++] = ESV( 3, BINDEX );
384 this_connect[i++] = FSV( 0, CINDEX );
385 this_connect[i++] = TSV( 0, DINDEX );
386 this_connect[i++] = FSV( 2, BINDEX );
387 result = mbImpl->create_element( MBHEX, this_connect, 8, this_hex );RR;
388 interstic_hexes.push_back( this_hex );
389
390 // now, the sphere interiors, four hexes per vertex sphere
391
392 // V0:
393 i = 0;
394 this_connect[i++] = CV( V0INDEX );
395 this_connect[i++] = ESV( 0, DINDEX );
396 this_connect[i++] = FSV( 3, EINDEX );
397 this_connect[i++] = ESV( 2, EINDEX );
398 this_connect[i++] = ESV( 3, DINDEX );
399 this_connect[i++] = FSV( 0, EINDEX );
400 this_connect[i++] = TSV( 0, FINDEX );
401 this_connect[i++] = FSV( 2, EINDEX );
402 result = mbImpl->create_element( MBHEX, this_connect, 8, this_hex );RR;
403 sphere_hexes.push_back( this_hex );
404
405 i = 0;
406 this_connect[i++] = ESV( 0, DINDEX );
407 this_connect[i++] = ESV( 0, AINDEX );
408 this_connect[i++] = FSV( 3, AINDEX );
409 this_connect[i++] = FSV( 3, EINDEX );
410 this_connect[i++] = FSV( 0, EINDEX );
411 this_connect[i++] = FSV( 0, AINDEX );
412 this_connect[i++] = TSV( 0, AINDEX );
413 this_connect[i++] = TSV( 0, FINDEX );
414 result = mbImpl->create_element( MBHEX, this_connect, 8, this_hex );RR;
415 sphere_hexes.push_back( this_hex );
416
417 i = 0;
418 this_connect[i++] = FSV( 3, EINDEX );
419 this_connect[i++] = FSV( 3, AINDEX );
420 this_connect[i++] = ESV( 2, BINDEX );
421 this_connect[i++] = ESV( 2, EINDEX );
422 this_connect[i++] = TSV( 0, FINDEX );
423 this_connect[i++] = TSV( 0, AINDEX );
424 this_connect[i++] = FSV( 2, AINDEX );
425 this_connect[i++] = FSV( 2, EINDEX );
426 result = mbImpl->create_element( MBHEX, this_connect, 8, this_hex );RR;
427 sphere_hexes.push_back( this_hex );
428
429 i = 0;
430 this_connect[i++] = TSV( 0, FINDEX );
431 this_connect[i++] = TSV( 0, AINDEX );
432 this_connect[i++] = FSV( 2, AINDEX );
433 this_connect[i++] = FSV( 2, EINDEX );
434 this_connect[i++] = FSV( 0, EINDEX );
435 this_connect[i++] = FSV( 0, AINDEX );
436 this_connect[i++] = ESV( 3, AINDEX );
437 this_connect[i++] = ESV( 3, DINDEX );
438 result = mbImpl->create_element( MBHEX, this_connect, 8, this_hex );RR;
439 sphere_hexes.push_back( this_hex );
440
441 // V1:
442 i = 0;
443 this_connect[i++] = CV( V1INDEX );
444 this_connect[i++] = ESV( 1, DINDEX );
445 this_connect[i++] = FSV( 3, GINDEX );
446 this_connect[i++] = ESV( 0, EINDEX );
447 this_connect[i++] = ESV( 4, DINDEX );
448 this_connect[i++] = FSV( 1, EINDEX );
449 this_connect[i++] = TSV( 0, GINDEX );
450 this_connect[i++] = FSV( 0, FINDEX );
451 result = mbImpl->create_element( MBHEX, this_connect, 8, this_hex );RR;
452 sphere_hexes.push_back( this_hex );
453
454 i = 0;
455 this_connect[i++] = FSV( 3, GINDEX );
456 this_connect[i++] = ESV( 1, DINDEX );
457 this_connect[i++] = ESV( 1, AINDEX );
458 this_connect[i++] = FSV( 3, CINDEX );
459 this_connect[i++] = TSV( 0, GINDEX );
460 this_connect[i++] = FSV( 1, EINDEX );
461 this_connect[i++] = FSV( 1, AINDEX );
462 this_connect[i++] = TSV( 0, BINDEX );
463 result = mbImpl->create_element( MBHEX, this_connect, 8, this_hex );RR;
464 sphere_hexes.push_back( this_hex );
465
466 i = 0;
467 this_connect[i++] = TSV( 0, GINDEX );
468 this_connect[i++] = FSV( 1, EINDEX );
469 this_connect[i++] = FSV( 1, AINDEX );
470 this_connect[i++] = TSV( 0, BINDEX );
471 this_connect[i++] = FSV( 0, FINDEX );
472 this_connect[i++] = ESV( 4, DINDEX );
473 this_connect[i++] = ESV( 4, AINDEX );
474 this_connect[i++] = FSV( 0, BINDEX );
475 result = mbImpl->create_element( MBHEX, this_connect, 8, this_hex );RR;
476 sphere_hexes.push_back( this_hex );
477
478 i = 0;
479 this_connect[i++] = ESV( 0, BINDEX );
480 this_connect[i++] = ESV( 0, EINDEX );
481 this_connect[i++] = FSV( 3, GINDEX );
482 this_connect[i++] = FSV( 3, CINDEX );
483 this_connect[i++] = FSV( 0, BINDEX );
484 this_connect[i++] = FSV( 0, FINDEX );
485 this_connect[i++] = TSV( 0, GINDEX );
486 this_connect[i++] = TSV( 0, BINDEX );
487 result = mbImpl->create_element( MBHEX, this_connect, 8, this_hex );RR;
488 sphere_hexes.push_back( this_hex );
489
490 // V2:
491 i = 0;
492 this_connect[i++] = ESV( 1, BINDEX );
493 this_connect[i++] = ESV( 1, EINDEX );
494 this_connect[i++] = FSV( 3, FINDEX );
495 this_connect[i++] = FSV( 3, BINDEX );
496 this_connect[i++] = FSV( 1, BINDEX );
497 this_connect[i++] = FSV( 1, FINDEX );
498 this_connect[i++] = TSV( 0, HINDEX );
499 this_connect[i++] = TSV( 0, CINDEX );
500 result = mbImpl->create_element( MBHEX, this_connect, 8, this_hex );RR;
501 sphere_hexes.push_back( this_hex );
502
503 i = 0;
504 this_connect[i++] = FSV( 3, FINDEX );
505 this_connect[i++] = ESV( 1, EINDEX );
506 this_connect[i++] = CV( V2INDEX );
507 this_connect[i++] = ESV( 2, DINDEX );
508 this_connect[i++] = TSV( 0, HINDEX );
509 this_connect[i++] = FSV( 1, FINDEX );
510 this_connect[i++] = ESV( 5, DINDEX );
511 this_connect[i++] = FSV( 2, GINDEX );
512 result = mbImpl->create_element( MBHEX, this_connect, 8, this_hex );RR;
513 sphere_hexes.push_back( this_hex );
514
515 i = 0;
516 this_connect[i++] = TSV( 0, HINDEX );
517 this_connect[i++] = FSV( 1, FINDEX );
518 this_connect[i++] = ESV( 5, DINDEX );
519 this_connect[i++] = FSV( 2, GINDEX );
520 this_connect[i++] = TSV( 0, CINDEX );
521 this_connect[i++] = FSV( 1, BINDEX );
522 this_connect[i++] = ESV( 5, AINDEX );
523 this_connect[i++] = FSV( 2, CINDEX );
524 result = mbImpl->create_element( MBHEX, this_connect, 8, this_hex );RR;
525 sphere_hexes.push_back( this_hex );
526
527 i = 0;
528 this_connect[i++] = FSV( 3, BINDEX );
529 this_connect[i++] = FSV( 3, FINDEX );
530 this_connect[i++] = ESV( 2, DINDEX );
531 this_connect[i++] = ESV( 2, AINDEX );
532 this_connect[i++] = TSV( 0, CINDEX );
533 this_connect[i++] = TSV( 0, HINDEX );
534 this_connect[i++] = FSV( 2, GINDEX );
535 this_connect[i++] = FSV( 2, CINDEX );
536 result = mbImpl->create_element( MBHEX, this_connect, 8, this_hex );RR;
537 sphere_hexes.push_back( this_hex );
538
539 // V3:
540 i = 0;
541 this_connect[i++] = FSV( 0, CINDEX );
542 this_connect[i++] = ESV( 4, BINDEX );
543 this_connect[i++] = FSV( 1, CINDEX );
544 this_connect[i++] = TSV( 0, DINDEX );
545 this_connect[i++] = FSV( 0, GINDEX );
546 this_connect[i++] = ESV( 4, EINDEX );
547 this_connect[i++] = FSV( 1, GINDEX );
548 this_connect[i++] = TSV( 0, IINDEX );
549 result = mbImpl->create_element( MBHEX, this_connect, 8, this_hex );RR;
550 sphere_hexes.push_back( this_hex );
551
552 i = 0;
553 this_connect[i++] = ESV( 3, BINDEX );
554 this_connect[i++] = FSV( 0, CINDEX );
555 this_connect[i++] = TSV( 0, DINDEX );
556 this_connect[i++] = FSV( 2, BINDEX );
557 this_connect[i++] = ESV( 3, EINDEX );
558 this_connect[i++] = FSV( 0, GINDEX );
559 this_connect[i++] = TSV( 0, IINDEX );
560 this_connect[i++] = FSV( 2, FINDEX );
561 result = mbImpl->create_element( MBHEX, this_connect, 8, this_hex );RR;
562 sphere_hexes.push_back( this_hex );
563
564 i = 0;
565 this_connect[i++] = TSV( 0, DINDEX );
566 this_connect[i++] = FSV( 1, CINDEX );
567 this_connect[i++] = ESV( 5, BINDEX );
568 this_connect[i++] = FSV( 2, BINDEX );
569 this_connect[i++] = TSV( 0, IINDEX );
570 this_connect[i++] = FSV( 1, GINDEX );
571 this_connect[i++] = ESV( 5, EINDEX );
572 this_connect[i++] = FSV( 2, FINDEX );
573 result = mbImpl->create_element( MBHEX, this_connect, 8, this_hex );RR;
574 sphere_hexes.push_back( this_hex );
575
576 i = 0;
577 this_connect[i++] = FSV( 0, GINDEX );
578 this_connect[i++] = ESV( 4, EINDEX );
579 this_connect[i++] = FSV( 1, GINDEX );
580 this_connect[i++] = TSV( 0, IINDEX );
581 this_connect[i++] = ESV( 3, EINDEX );
582 this_connect[i++] = CV( V3INDEX );
583 this_connect[i++] = ESV( 5, EINDEX );
584 this_connect[i++] = FSV( 2, FINDEX );
585 result = mbImpl->create_element( MBHEX, this_connect, 8, this_hex );RR;
586 sphere_hexes.push_back( this_hex );
587
588 return result;
589 }
590
591 ErrorCode SphereDecomp::retrieve_subdiv_verts( EntityHandle tet,
592 EntityHandle this_ent,
593 const EntityHandle* tet_conn,
594 const int dim,
595 EntityHandle* subdiv_verts )
596 {
597 // get the subdiv verts for this entity
598 ErrorCode result;
599
600 // if it's a tet, just put them on the end & return
601 if( tet == this_ent )
602 {
603 result = mbImpl->tag_get_data( subdivVerticesTag, &this_ent, 1, &subdiv_verts[90] );
604 return MB_SUCCESS;
605 }
606
607 // if it's a sub-entity, need to find index, relative orientation, and offset
608 // get connectivity of sub-entity
609 std::vector< EntityHandle > this_conn;
610 result = mbImpl->get_connectivity( &this_ent, 1, this_conn );RR;
611
612 // get relative orientation
613 std::vector< int > conn_tet_indices( this_conn.size() );
614 for( size_t i = 0; i < this_conn.size(); ++i )
615 conn_tet_indices[i] = std::find( tet_conn, tet_conn + 4, this_conn[i] ) - tet_conn;
616 int sense, side_no, offset;
617 int success = CN::SideNumber( MBTET, &conn_tet_indices[0], this_conn.size(), dim, side_no, sense, offset );
618 if( -1 == success ) return MB_FAILURE;
619
620 // start of this entity's subdiv_verts; edges go first, then preceding sides, then this one;
621 // this assumes 6 edges/tet
622 EntityHandle* subdiv_start = &subdiv_verts[( ( dim - 1 ) * 6 + side_no ) * 9];
623
624 // get subdiv_verts and put them into proper place
625 result = mbImpl->tag_get_data( subdivVerticesTag, &this_ent, 1, subdiv_start );
626
627 // could probably do this more elegantly, but isn't worth it
628 #define SWITCH( a, b ) \
629 { \
630 EntityHandle tmp_handle = a; \
631 ( a ) = b; \
632 ( b ) = tmp_handle; \
633 }
634 switch( dim )
635 {
636 case 1:
637 if( offset != 0 || sense == -1 )
638 {
639 SWITCH( subdiv_start[0], subdiv_start[1] );
640 SWITCH( subdiv_start[3], subdiv_start[4] );
641 }
642 break;
643 case 2:
644 // rotate first
645 if( 0 != offset )
646 {
647 std::rotate( subdiv_start, subdiv_start + offset, subdiv_start + 3 );
648 std::rotate( subdiv_start + 4, subdiv_start + 4 + offset, subdiv_start + 7 );
649 }
650 // now flip, if necessary
651 if( -1 == sense )
652 {
653 SWITCH( subdiv_start[1], subdiv_start[2] );
654 SWITCH( subdiv_start[5], subdiv_start[6] );
655 }
656 break;
657 default:
658 return MB_FAILURE;
659 }
660
661 // ok, we're done
662 return MB_SUCCESS;
663 }
1.9.1.