Mesh Oriented datABase  (version 5.5.0)
An array-based unstructured mesh library
SphereDecomp.cpp
Go to the documentation of this file.
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 }