MOAB: Mesh Oriented datABase  (version 5.5.0)
elem_eval_test.cpp
Go to the documentation of this file.
1 /**
2  * \file elem_eval_test.cpp
3  *
4  * \brief test ElemEvaluator and the various element types in MOAB
5  *
6  */
7 #include "moab/Core.hpp"
8 #include "moab/ReadUtilIface.hpp"
9 #include "moab/ElemEvaluator.hpp"
10 #include "moab/LocalDiscretization/LinearTri.hpp"
11 #include "moab/LocalDiscretization/LinearQuad.hpp"
12 #include "moab/LocalDiscretization/LinearHex.hpp"
13 #include "moab/LocalDiscretization/LinearTet.hpp"
14 #include "moab/LocalDiscretization/QuadraticHex.hpp"
15 #include "moab/CartVect.hpp"
16 #include "TestUtil.hpp"
17 
18 using namespace moab;
19 
20 void test_linear_tri();
21 void test_linear_quad();
22 void test_linear_hex();
23 void test_linear_tet();
24 void test_quadratic_hex();
25 void test_normal_linear_tri();
26 void test_normal_linear_quad();
27 void test_normal_linear_tet();
28 void test_normal_linear_hex();
29 ErrorCode create_mesh( Core& mb, EntityType type );
30 
31 CartVect hex_verts[] = {
32  // corners
33  CartVect( -1, -1, -1 ), CartVect( 1, -1, -1 ), CartVect( 1, 1, -1 ), CartVect( -1, 1, -1 ), CartVect( -1, -1, 1 ),
34  CartVect( 1, -1, 1 ), CartVect( 1, 1, 1 ), CartVect( -1, 1, 1 ),
35  // mid-edge (bottom, middle, top)
36  CartVect( 0, -1, -1 ), CartVect( 1, 0, -1 ), CartVect( 0, 1, -1 ), CartVect( -1, 0, -1 ), CartVect( -1, -1, 0 ),
37  CartVect( 1, -1, 0 ), CartVect( 1, 1, 0 ), CartVect( -1, 1, 0 ), CartVect( 0, -1, 1 ), CartVect( 1, 0, 1 ),
38  CartVect( 0, 1, 1 ), CartVect( -1, 0, 1 ),
39  // mid-face (middle, bottom, top)
40  CartVect( 0, -1, 0 ), CartVect( 1, 0, 0 ), CartVect( 0, 1, 0 ), CartVect( -1, 0, 0 ), CartVect( 0, 0, -1 ),
41  CartVect( 0, 0, 1 ),
42  // mid-element
43  CartVect( 0, 0, 0 ) };
44 
45 const double EPS1 = 1.0e-6;
46 
47 void test_eval( ElemEvaluator& ee, bool test_integrate )
48 {
49 
50  CartVect params, posn, params2;
51  int is_inside;
52  Matrix3 jacob;
53  ErrorCode rval;
54  int ent_dim = ee.get_moab()->dimension_from_handle( ee.get_ent_handle() );
55 
56  for( params[0] = -1; params[0] <= 1; params[0] += 0.2 )
57  {
58  for( params[1] = -1; params[1] <= 1; params[1] += 0.2 )
59  {
60  for( params[2] = -1; params[2] <= 1; params[2] += 0.2 )
61  {
62 
63  // forward/reverse evaluation should get back to the same point, within tol
64  if( !ee.inside( params.array(), EPS1 ) ) continue;
65 
66  rval = ee.eval( params.array(), posn.array() );CHECK_ERR( rval );
67  rval = ee.reverse_eval( posn.array(), EPS1, EPS1, params2.array(), &is_inside );CHECK_ERR( rval );
68  CHECK_REAL_EQUAL( 0.0, params[0] - params2[0], 3 * EPS1 );
69  if( ent_dim > 1 ) CHECK_REAL_EQUAL( 0.0, params[1] - params2[1], 3 * EPS1 );
70  if( ent_dim > 2 ) CHECK_REAL_EQUAL( 0.0, params[2] - params2[2], 3 * EPS1 );
71 
72  // jacobian should be >= 0
73  rval = ee.jacobian( params.array(), jacob.array() );CHECK_ERR( rval );
74  CHECK( jacob.determinant() >= 0.0 );
75  }
76  }
77  }
78 
79  if( !test_integrate ) return;
80 
81  // tag equal to coordinates should integrate to avg position, test that
82  Tag tag;
83  // make a temporary tag and set it on vertices to the vertex positions
84  rval = ee.get_moab()->tag_get_handle( NULL, 3, MB_TYPE_DOUBLE, tag, MB_TAG_DENSE | MB_TAG_CREAT );CHECK_ERR( rval );
85  rval = ee.get_moab()->tag_set_data( tag, ee.get_vert_handles(), ee.get_num_verts(), ee.get_vert_pos() );CHECK_ERR( rval );
86  // set that temporary tag on the evaluator so that's what gets integrated
87  rval = ee.set_tag_handle( tag, 0 );CHECK_ERR( rval );
88 
89  CartVect integral, avg;
90  rval = ee.integrate( integral.array() );CHECK_ERR( rval );
91 
92  // now integrate a const 1-valued function, using direct call to the integrate function
93  std::vector< double > one( ee.get_num_verts(), 1.0 );
94  double measure;
95  EvalSet es;
96  EntityHandle eh = ee.get_ent_handle();
97  rval = EvalSet::get_eval_set( ee.get_moab(), eh, es );CHECK_ERR( rval );
98  rval = ( *es.integrateFcn )( &one[0], ee.get_vert_pos(), ee.get_num_verts(),
99  ee.get_moab()->dimension_from_handle( eh ), 1, ee.get_work_space(), &measure );CHECK_ERR( rval );
100  if( measure ) integral /= measure;
101 
102  // check against avg of entity's vertices' positions
103  rval = ee.get_moab()->get_coords( &eh, 1, avg.array() );CHECK_ERR( rval );
104  CHECK_REAL_EQUAL( 0.0, ( avg - integral ).length(), EPS1 );
105 
106  // delete the temporary tag
107  rval = ee.get_moab()->tag_delete( tag );CHECK_ERR( rval );
108  // set the ee's tag back to coords
109  rval = ee.set_tag_handle( 0, 0 );CHECK_ERR( rval );
110 }
111 
112 void test_evals( ElemEvaluator& ee,
113  bool test_integrate,
114  EntityHandle* ents,
115  int num_ents,
116  Tag onetag,
117  double total_vol )
118 {
119  for( int i = 0; i < num_ents; i++ )
120  {
121  ee.set_ent_handle( ents[i] );
122  test_eval( ee, false );
123  }
124 
125  if( !test_integrate ) return;
126  if( !num_ents || !ents ) CHECK_ERR( MB_FAILURE );
127 
128  ErrorCode rval = ee.set_tag_handle( onetag, 0 );CHECK_ERR( rval );
129 
130  double tot_vol = 0.0;
131  for( int i = 0; i < num_ents; i++ )
132  {
133  double tmp_vol;
134  ee.set_ent_handle( ents[i] );
135  rval = ee.integrate( &tmp_vol );CHECK_ERR( rval );
136  tot_vol += tmp_vol;
137  }
138  CHECK_REAL_EQUAL( total_vol, tot_vol, EPS1 );
139 }
140 
141 int main()
142 {
143  int failures = 0;
144 
145  failures += RUN_TEST( test_linear_tri ); // currently failing linear tri, bad formulation, working on it...
146  failures += RUN_TEST( test_linear_quad );
147  failures += RUN_TEST( test_linear_hex );
148  failures += RUN_TEST( test_quadratic_hex );
149  failures += RUN_TEST( test_linear_tet );
150  failures += RUN_TEST( test_normal_linear_tri );
151  failures += RUN_TEST( test_normal_linear_quad );
152  failures += RUN_TEST( test_normal_linear_tet );
153  failures += RUN_TEST( test_normal_linear_hex );
154 
155  return failures;
156 }
157 
158 void test_linear_tri()
159 {
160  Core mb;
161  Range verts;
162  double tri_verts[] = { -1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, -1.0 };
163 
164  ErrorCode rval = mb.create_vertices( tri_verts, 3, verts );CHECK_ERR( rval );
165  EntityHandle tri;
166  std::vector< EntityHandle > connect;
167  std::copy( verts.begin(), verts.end(), std::back_inserter( connect ) );
168  rval = mb.create_element( MBTRI, &connect[0], 3, tri );CHECK_ERR( rval );
169 
170  ElemEvaluator ee( &mb, tri, 0 );
171  ee.set_tag_handle( 0, 0 );
172  ee.set_eval_set( MBTRI, LinearTri::eval_set() );
173 
174  test_eval( ee, true );
175 }
176 
177 void test_linear_quad()
178 {
179  Core mb;
180  Range verts;
181  ErrorCode rval = mb.create_vertices( (double*)hex_verts, 4, verts );CHECK_ERR( rval );
182  EntityHandle quad;
183  std::vector< EntityHandle > connect;
184  std::copy( verts.begin(), verts.end(), std::back_inserter( connect ) );
185  rval = mb.create_element( MBQUAD, &connect[0], 4, quad );CHECK_ERR( rval );
186 
187  ElemEvaluator ee( &mb, quad, 0 );
188  ee.set_tag_handle( 0, 0 );
189  ee.set_eval_set( MBQUAD, LinearQuad::eval_set() );
190 
191  test_eval( ee, true );
192 }
193 
194 void test_linear_hex()
195 {
196  Core mb;
197  Range verts;
198  ErrorCode rval = mb.create_vertices( (double*)hex_verts, 8, verts );CHECK_ERR( rval );
199  EntityHandle hex;
200  std::vector< EntityHandle > connect;
201  std::copy( verts.begin(), verts.end(), std::back_inserter( connect ) );
202  rval = mb.create_element( MBHEX, &connect[0], 8, hex );CHECK_ERR( rval );
203 
204  ElemEvaluator ee( &mb, hex, 0 );
205  ee.set_tag_handle( 0, 0 );
206  ee.set_eval_set( MBHEX, LinearHex::eval_set() );
207 
208  test_eval( ee, true );
209 }
210 
211 void test_quadratic_hex()
212 {
213  Core mb;
214  Range verts;
215  ErrorCode rval = mb.create_vertices( (double*)hex_verts, 27, verts );CHECK_ERR( rval );
216  EntityHandle hex;
217  std::vector< EntityHandle > connect;
218  std::copy( verts.begin(), verts.end(), std::back_inserter( connect ) );
219  rval = mb.create_element( MBHEX, &connect[0], 27, hex );CHECK_ERR( rval );
220 
221  ElemEvaluator ee( &mb, hex, 0 );
222  ee.set_tag_handle( 0, 0 );
223  ee.set_eval_set( MBHEX, QuadraticHex::eval_set() );
224 
225  test_eval( ee, false );
226 }
227 
228 void test_linear_tet()
229 {
230  Core mb;
231  Range verts;
232  ErrorCode rval = mb.create_vertices( (double*)hex_verts[0].array(), 8, verts );CHECK_ERR( rval );
233  EntityHandle starth = 1, *conn;
234  int conn_inds[] = { 1, 6, 4, 5, 1, 4, 6, 3, 0, 1, 3, 4, 1, 2, 3, 6, 3, 4, 6, 7 };
235  ReadUtilIface* ru;
236  rval = mb.query_interface( ru );CHECK_ERR( rval );
237  rval = ru->get_element_connect( 5, 4, MBTET, 1, starth, conn );CHECK_ERR( rval );
238  for( unsigned int i = 0; i < 20; i++ )
239  conn[i] = verts[conn_inds[i]];
240  EntityHandle tets[5];
241  for( unsigned int i = 0; i < 5; i++ )
242  tets[i] = starth + i;
243  ElemEvaluator ee( &mb, 0, 0 );
244  ee.set_tag_handle( 0, 0 );
245  ee.set_eval_set( MBTET, LinearTet::eval_set() );
246 
247  // make a tag whose value is one on all vertices
248  Tag tag;
249  rval = mb.tag_get_handle( NULL, 1, MB_TYPE_DOUBLE, tag, MB_TAG_DENSE | MB_TAG_CREAT );CHECK_ERR( rval );
250  std::vector< double > vals( verts.size(), 1.0 );
251  rval = mb.tag_set_data( tag, verts, &vals[0] );CHECK_ERR( rval );
252 
253  test_evals( ee, true, tets, 5, tag, 8.0 );
254 }
255 
256 void test_normal_linear_tri()
257 {
258  ErrorCode error;
259  Core mb;
260 
261  error = create_mesh( mb, MBTRI );CHECK_ERR( error );
262  Range faces;
263  error = mb.get_entities_by_dimension( 0, 2, faces );CHECK_ERR( error );
264 
265  ElemEvaluator ee( &mb, 0, 0 );
266  ee.set_eval_set( MBTRI, LinearTri::eval_set() );
267 
268  double nrms[8][3];
269 
270  for( int i = 0; i < 4; i++ )
271  {
272  ee.set_ent_handle( faces[i] );
273  ee.get_normal( 1, 0, nrms[2 * i] );
274  ee.get_normal( 1, 2, nrms[2 * i + 1] );
275  }
276 
277  for( int i = 0; i < 4; i++ )
278  {
279  if( i == 3 )
280  {
281  double val = nrms[7][0] * nrms[0][0] + nrms[7][1] * nrms[0][1] + nrms[7][2] * nrms[0][2];
282  CHECK_EQUAL( val, -1 );
283  }
284  else
285  {
286  double val = nrms[2 * i + 1][0] * nrms[2 * ( i + 1 )][0] + nrms[2 * i + 1][1] * nrms[2 * ( i + 1 )][1] +
287  nrms[2 * i + 1][2] * nrms[2 * ( i + 1 )][2];
288  CHECK_EQUAL( val, -1 );
289  }
290  }
291 }
292 
293 void test_normal_linear_quad()
294 {
295  ErrorCode error;
296  Core mb;
297 
298  error = create_mesh( mb, MBQUAD );CHECK_ERR( error );
299  Range faces;
300  error = mb.get_entities_by_dimension( 0, 2, faces );CHECK_ERR( error );
301 
302  ElemEvaluator ee( &mb, 0, 0 );
303  ee.set_eval_set( MBQUAD, LinearQuad::eval_set() );
304 
305  double nrms[8][3];
306 
307  for( int i = 0; i < 4; i++ )
308  {
309  ee.set_ent_handle( faces[i] );
310  ee.get_normal( 1, 0, nrms[2 * i] );
311  ee.get_normal( 1, 3, nrms[2 * i + 1] );
312  }
313 
314  for( int i = 0; i < 4; i++ )
315  {
316  if( i == 3 )
317  {
318  double val = nrms[7][0] * nrms[0][0] + nrms[7][1] * nrms[0][1] + nrms[7][2] * nrms[0][2];
319  CHECK_EQUAL( val, -1 );
320  }
321  else
322  {
323  double val = nrms[2 * i + 1][0] * nrms[2 * ( i + 1 )][0] + nrms[2 * i + 1][1] * nrms[2 * ( i + 1 )][1] +
324  nrms[2 * i + 1][2] * nrms[2 * ( i + 1 )][2];
325  CHECK_EQUAL( val, -1 );
326  }
327  }
328 }
329 
330 void test_normal_linear_tet()
331 {
332  ErrorCode error;
333  Core mb;
334 
335  error = create_mesh( mb, MBTET );CHECK_ERR( error );
336  Range cells;
337  error = mb.get_entities_by_dimension( 0, 3, cells );CHECK_ERR( error );
338 
339  ElemEvaluator ee( &mb, 0, 0 );
340  ee.set_eval_set( MBTET, LinearTet::eval_set() );
341 
342  double nrms[8][3];
343 
344  for( int i = 0; i < 4; i++ )
345  {
346  ee.set_ent_handle( cells[i] );
347  ee.get_normal( 2, 0, nrms[2 * i] );
348  ee.get_normal( 2, 2, nrms[2 * i + 1] );
349  }
350 
351  for( int i = 0; i < 4; i++ )
352  {
353  if( i == 3 )
354  {
355  double val = nrms[7][0] * nrms[0][0] + nrms[7][1] * nrms[0][1] + nrms[7][2] * nrms[0][2];
356  CHECK_EQUAL( val, -1 );
357  }
358  else
359  {
360  double val = nrms[2 * i + 1][0] * nrms[2 * ( i + 1 )][0] + nrms[2 * i + 1][1] * nrms[2 * ( i + 1 )][1] +
361  nrms[2 * i + 1][2] * nrms[2 * ( i + 1 )][2];
362  CHECK_EQUAL( val, -1 );
363  }
364  }
365 }
366 
367 void test_normal_linear_hex()
368 {
369  ErrorCode error;
370  Core mb;
371 
372  error = create_mesh( mb, MBHEX );CHECK_ERR( error );
373  Range cells;
374  error = mb.get_entities_by_dimension( 0, 3, cells );CHECK_ERR( error );
375 
376  ElemEvaluator ee( &mb, 0, 0 );
377  ee.set_eval_set( MBHEX, LinearHex::eval_set() );
378 
379  double nrms[8][3];
380 
381  for( int i = 0; i < 4; i++ )
382  {
383  ee.set_ent_handle( cells[i] );
384  ee.get_normal( 2, 0, nrms[2 * i] );
385  ee.get_normal( 2, 3, nrms[2 * i + 1] );
386  }
387 
388  for( int i = 0; i < 4; i++ )
389  {
390  if( i == 3 )
391  {
392  double val = nrms[7][0] * nrms[0][0] + nrms[7][1] * nrms[0][1] + nrms[7][2] * nrms[0][2];
393  CHECK_EQUAL( val, -1 );
394  }
395  else
396  {
397  double val = nrms[2 * i + 1][0] * nrms[2 * ( i + 1 )][0] + nrms[2 * i + 1][1] * nrms[2 * ( i + 1 )][1] +
398  nrms[2 * i + 1][2] * nrms[2 * ( i + 1 )][2];
399  CHECK_EQUAL( val, -1 );
400  }
401  }
402 }
403 
404 ErrorCode create_mesh( Core& mb, EntityType type )
405 {
406  ErrorCode error;
407  if( type == MBTRI )
408  {
409  const double coords[] = { 0, 0, 0, 1, 0, 0, 0, 1, 0, -1, 0, 0, 0, -1, 0 };
410  const size_t num_vtx = sizeof( coords ) / sizeof( double ) / 3;
411 
412  const int conn[] = { 0, 1, 2, 0, 2, 3, 0, 3, 4, 0, 4, 1 };
413  const size_t num_elems = sizeof( conn ) / sizeof( int ) / 3;
414 
415  EntityHandle verts[num_vtx], faces[num_elems];
416  for( size_t i = 0; i < num_vtx; ++i )
417  {
418  error = mb.create_vertex( coords + 3 * i, verts[i] );CHECK_ERR( error );
419  }
420 
421  for( size_t i = 0; i < num_elems; ++i )
422  {
423  EntityHandle c[3];
424  for( int j = 0; j < 3; j++ )
425  c[j] = verts[conn[3 * i + j]];
426 
427  error = mb.create_element( MBTRI, c, 3, faces[i] );CHECK_ERR( error );
428  }
429  }
430  else if( type == MBQUAD )
431  {
432  const double coords[] = { 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 1, 0, -1, 1,
433  0, -1, 0, 0, -1, -1, 0, 0, -1, 0, 1, -1, 0 };
434  const size_t num_vtx = sizeof( coords ) / sizeof( double ) / 3;
435 
436  const int conn[] = { 0, 1, 2, 3, 0, 3, 4, 5, 0, 5, 6, 7, 0, 7, 8, 1 };
437  const size_t num_elems = sizeof( conn ) / sizeof( int ) / 4;
438 
439  EntityHandle verts[num_vtx], faces[num_elems];
440  for( size_t i = 0; i < num_vtx; ++i )
441  {
442  error = mb.create_vertex( coords + 3 * i, verts[i] );CHECK_ERR( error );
443  }
444 
445  for( size_t i = 0; i < num_elems; ++i )
446  {
447  EntityHandle c[4];
448  for( int j = 0; j < 4; j++ )
449  c[j] = verts[conn[4 * i + j]];
450 
451  error = mb.create_element( MBQUAD, c, 4, faces[i] );CHECK_ERR( error );
452  }
453  }
454  else if( type == MBTET )
455  {
456  const double coords[] = { 0, 0, 0, 1, 0, 0, 0, 1, 0, -1, 0, 0, 0, -1, 0, 0, 0, 1 };
457  const size_t num_vtx = sizeof( coords ) / sizeof( double ) / 3;
458 
459  const int conn[] = { 0, 1, 2, 5, 0, 2, 3, 5, 0, 3, 4, 5, 0, 4, 1, 5 };
460  const size_t num_elems = sizeof( conn ) / sizeof( int ) / 4;
461 
462  EntityHandle verts[num_vtx], cells[num_elems];
463  for( size_t i = 0; i < num_vtx; ++i )
464  {
465  error = mb.create_vertex( coords + 3 * i, verts[i] );CHECK_ERR( error );
466  }
467 
468  for( size_t i = 0; i < num_elems; ++i )
469  {
470  EntityHandle c[4];
471  for( int j = 0; j < 4; j++ )
472  c[j] = verts[conn[4 * i + j]];
473 
474  error = mb.create_element( MBTET, c, 4, cells[i] );CHECK_ERR( error );
475  }
476  }
477  else if( type == MBHEX )
478  {
479  const double coords[] = { 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 1, 0, -1, 1, 0, -1, 0, 0,
480  -1, -1, 0, 0, -1, 0, 1, -1, 0, 0, 0, 1, 1, 0, 1, 1, 1, 1,
481  0, 1, 1, -1, 1, 1, -1, 0, 1, -1, -1, 1, 0, -1, 1, 1, -1, 1 };
482  const size_t num_vtx = sizeof( coords ) / sizeof( double ) / 3;
483 
484  const int conn[] = { 0, 1, 2, 3, 9, 10, 11, 12, 0, 3, 4, 5, 9, 12, 13, 14,
485  0, 5, 6, 7, 9, 14, 15, 16, 0, 7, 8, 1, 9, 16, 17, 10 };
486  const size_t num_elems = sizeof( conn ) / sizeof( int ) / 8;
487 
488  EntityHandle verts[num_vtx], cells[num_elems];
489  for( size_t i = 0; i < num_vtx; ++i )
490  {
491  error = mb.create_vertex( coords + 3 * i, verts[i] );CHECK_ERR( error );
492  }
493 
494  for( size_t i = 0; i < num_elems; ++i )
495  {
496  EntityHandle c[8];
497  for( int j = 0; j < 8; j++ )
498  c[j] = verts[conn[8 * i + j]];
499 
500  error = mb.create_element( MBHEX, c, 8, cells[i] );CHECK_ERR( error );
501  }
502  }
503 
504  return MB_SUCCESS;
505 }