Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
ElementTest.cpp
Go to the documentation of this file.
1 #include "TestUtil.hpp"
2 #include "ElemUtil.hpp"
3 #include "moab/Core.hpp"
4 #include "moab/Range.hpp"
5 #include <iostream>
6 
7 #ifndef MESHDIR
8 #error Specify MESHDIR to compile test
9 #endif
10 
11 using namespace moab;
12 
13 void test_tet();
14 void test_hex();
15 void test_spectral_hex();
16 void test_spectral_quad();
17 void test_spherical_quad();
18 void test_linear_tri();
19 void test_spherical_tri();
20 
21 int main()
22 {
23  int rval = 0;
24  rval += RUN_TEST( test_tet );
25  rval += RUN_TEST( test_hex );
26  rval += RUN_TEST( test_spectral_hex );
27  rval += RUN_TEST( test_spectral_quad );
28  rval += RUN_TEST( test_spherical_quad );
29  rval += RUN_TEST( test_linear_tri );
30  rval += RUN_TEST( test_spherical_tri );
31  return rval;
32 }
33 
34 void test_tet()
35 {
37 } // test_tet()
38 
39 void test_hex()
40 {
41  double positions[] = { 236.80706050970281, -139.55422526228017, 193.27999999999997,
42  236.47511729348639, -141.33020962638582, 193.27999999999997,
43  237.8457938295229, -142.57076074835663, 193.27999999999997,
44  239.12702305519684, -139.96608933577852, 193.27999999999997,
45  236.80841321361444, -139.55341321335499, 202.654,
46  236.47655014713746, -141.32980272396816, 202.654,
47  237.8477913707564, -142.57047282187165, 202.654,
48  239.12865103844533, -139.96531051891105, 202.654 };
49  CartVect x( 235.96518686964933, -142.43503000077749, 188.19999999999987 );
50  std::vector< CartVect > vertices;
51  for( int i = 0; i < 8; i++ )
52  vertices.push_back( CartVect( positions + 3 * i ) );
53 
54  moab::Element::LinearHex hex( vertices );
55  double tol( 0.0001 );
56  if( hex.inside_box( x, tol ) )
57  {
58  CartVect nat_par = hex.ievaluate( x, 0.0001 );
59  std::cout << nat_par << "\n";
60  }
61 
62  double positions2[] = { 49.890500000000024, -20.376134375374882, 312.72000000000003,
63  52.015875000000044, -19.149048546996006, 312.72000000000003,
64  48.430375821458099, -18.548796774572125, 312.72000000000003,
65  47.717616239031223, -21.191360829777231, 312.72000000000003,
66  49.890500000000024, -20.376134375374882, 322.88,
67  52.015875000000044, -19.149048546996006, 322.88,
68  48.429930354643275, -18.52828610485021, 322.88,
69  47.720552036968819, -21.167591146685712, 322.88 };
70 
71  CartVect x2( 51.469000000000015, -20.145482942833631, 317.80000000000001 );
72 
73  vertices.clear();
74  for( int i = 0; i < 8; i++ )
75  vertices.push_back( CartVect( positions2 + 3 * i ) );
76  moab::Element::LinearHex hex2( vertices );
77  if( hex2.inside_box( x2, tol ) )
78  {
79  try
80  {
81  CartVect nat_par = hex.ievaluate( x, 0.0001 );
82  std::cout << nat_par << "\n";
83  }
85  {
86  // nothing
87  }
88  }
89 
90 } // test_hex()
91 
93 {
94  // first load a model that has spectral elements
95  moab::Core* mb = new moab::Core();
96  std::string meshFile = STRINGIFY( MESHDIR ) "/spectral.h5m";
97  moab::ErrorCode rval = mb->load_mesh( meshFile.c_str() );
98  if( moab::MB_SUCCESS != rval )
99  {
100  std::cout << "Problems reading file " << meshFile << "." << std::endl;
101  delete mb;
102  return;
103  }
104 
105  // get the ent set with SEM_DIMS tag
106  moab::Range spectral_sets;
107  moab::Tag sem_tag;
108  rval = mb->tag_get_handle( "SEM_DIMS", 3, moab::MB_TYPE_INTEGER, sem_tag );
109  if( moab::MB_SUCCESS != rval )
110  {
111  std::cout << "can't find tag, no spectral set\n";
112  delete mb;
113  return;
114  }
115  rval = mb->get_entities_by_type_and_tag( 0, moab::MBENTITYSET, &sem_tag, NULL, 1, spectral_sets );
116  if( moab::MB_SUCCESS != rval || spectral_sets.empty() )
117  {
118  std::cout << "can't get sem set\n";
119  delete mb;
120  return;
121  }
122 
123  moab::Range ents;
124 
125  int sem_dims[3];
126  moab::EntityHandle firstSemSet = spectral_sets[0];
127  rval = mb->tag_get_data( sem_tag, &firstSemSet, 1, (void*)sem_dims );
128  if( moab::MB_SUCCESS != rval )
129  {
130  delete mb;
131  return;
132  }
133 
134  rval = mb->get_entities_by_dimension( firstSemSet, 3, ents );
135  if( moab::MB_SUCCESS != rval )
136  {
137  delete mb;
138  return;
139  }
140  std::cout << "Found " << ents.size() << " " << 3 << "-dimensional entities:" << std::endl;
141 
142  if( sem_dims[0] != sem_dims[1] || sem_dims[0] != sem_dims[2] )
143  {
144  std::cout << " dimensions are different. bail out\n";
145  delete mb;
146  return;
147  }
148 
149  // get the SEM_X ...tags
150  moab::Tag xm1Tag, ym1Tag, zm1Tag;
151  int ntot = sem_dims[0] * sem_dims[1] * sem_dims[2];
152  rval = mb->tag_get_handle( "SEM_X", ntot, moab::MB_TYPE_DOUBLE, xm1Tag );
153  if( moab::MB_SUCCESS != rval )
154  {
155  std::cout << "can't get xm1tag \n";
156  delete mb;
157  return;
158  }
159  rval = mb->tag_get_handle( "SEM_Y", ntot, moab::MB_TYPE_DOUBLE, ym1Tag );
160  if( moab::MB_SUCCESS != rval )
161  {
162  std::cout << "can't get ym1tag \n";
163  delete mb;
164  return;
165  }
166  rval = mb->tag_get_handle( "SEM_Z", ntot, moab::MB_TYPE_DOUBLE, zm1Tag );
167  if( moab::MB_SUCCESS != rval )
168  {
169  std::cout << "can't get zm1tag \n";
170  delete mb;
171  return;
172  }
173  moab::Tag velTag;
174 
175  rval = mb->tag_get_handle( "VX", ntot, moab::MB_TYPE_DOUBLE, velTag );
176  if( moab::MB_SUCCESS != rval )
177  {
178  std::cout << "can't get veltag \n";
179  delete mb;
180  return;
181  }
182  moab::Element::SpectralHex specHex( sem_dims[0] );
183 
184  // compute the data for some elements
185  for( moab::Range::iterator rit = ents.begin(); rit != ents.end(); ++rit )
186  {
187  // get the tag pointers to the internal storage for xm1, to not copy the values
188  moab::EntityHandle eh = *rit;
189  const double* xval;
190  const double* yval;
191  const double* zval;
192  rval = mb->tag_get_by_ptr( xm1Tag, &eh, 1, (const void**)&xval );
193  if( moab::MB_SUCCESS != rval )
194  {
195  std::cout << "can't get xm1 values \n";
196  delete mb;
197  return;
198  }
199  rval = mb->tag_get_by_ptr( ym1Tag, &eh, 1, (const void**)&yval );
200  if( moab::MB_SUCCESS != rval )
201  {
202  std::cout << "can't get ym1 values \n";
203  delete mb;
204  return;
205  }
206  rval = mb->tag_get_by_ptr( zm1Tag, &eh, 1, (const void**)&zval );
207  if( moab::MB_SUCCESS != rval )
208  {
209  std::cout << "can't get zm1 values \n";
210  delete mb;
211  return;
212  }
213  if( rit == ents.begin() )
214  {
215  std::cout << " xm1 for first element: \n";
216  for( int i = 0; i < ntot; i++ )
217  std::cout << " " << xval[i];
218  std::cout << "\n";
219  }
220  specHex.set_gl_points( (double*)xval, (double*)yval, (double*)zval );
221  // first evaluate a point, then inverse it to see if we get the same thing
222  moab::CartVect rst( 0.1, -0.1, 0.5 );
223  moab::CartVect pos = specHex.evaluate( rst );
224  moab::CartVect inverse = specHex.ievaluate( pos, 0.0001 );
225  std::cout << "difference" << rst - inverse << "\n";
226  Matrix3 jac = specHex.jacobian( rst );
227  std::cout << "jacobian: \n" << jac << " \n determinant: " << jac.determinant() << "\n";
228  // evaluate vx at some point
229  const double* vx;
230  rval = mb->tag_get_by_ptr( velTag, &eh, 1, (const void**)&vx );
231  if( moab::MB_SUCCESS != rval )
232  {
233  std::cout << "can't get vel values \n";
234  delete mb;
235  return;
236  }
237  double vel1 = specHex.evaluate_scalar_field( rst, vx );
238  std::cout << "velocity: " << vel1 << "\n";
239  // compute integral over vx:
240  double integral = specHex.integrate_scalar_field( vx );
241  std::cout << "integral over vx: " << integral << "\n";
242  }
243  std::cout << "success...\n";
244 
245  delete mb;
246 }
247 
249 {
250  // first load a model that has spectral elements
251  moab::Core* mb = new moab::Core();
252  // use the grid on Sphere from mbcslam folder
253  std::string meshFile = STRINGIFY( MESHDIR ) "/mbcslam/eulerHomme.vtk";
254  moab::ErrorCode rval = mb->load_mesh( meshFile.c_str() );
255  if( moab::MB_SUCCESS != rval )
256  {
257  std::cout << "Problems reading file " << meshFile << "." << std::endl;
258  delete mb;
259  return;
260  }
261 
262  // for each element, compute the gl points and project them on sphere
263  // the radius is the same as the one from intersection test on sphere
264  // double R = 6. * sqrt(3.) / 2; // input
265 
266  moab::Range ents;
267 
268  rval = mb->get_entities_by_type( 0, moab::MBQUAD, ents ); // get all quads
269  if( moab::MB_SUCCESS != rval )
270  {
271  delete mb;
272  return;
273  }
274  std::cout << "Found " << ents.size() << " " << 2 << "-dimensional entities:" << std::endl;
275 
276  //
277  int NP = 4; // test this....
278  moab::Element::SpectralQuad specQuad( NP );
279 
280  // compute the gl points for some elements
281  for( moab::Range::iterator rit = ents.begin(); rit != ents.end(); ++rit )
282  {
283 
284  const moab::EntityHandle* conn4 = NULL;
285  int num_nodes = 0;
286  rval = mb->get_connectivity( *rit, conn4, num_nodes );
287  if( moab::MB_SUCCESS != rval )
288  {
289  std::cout << "can't get connectivity for quad \n";
290  delete mb;
291  return;
292  }
293  assert( num_nodes == 4 );
294 
295  std::vector< CartVect > verts( 4 );
296  rval = mb->get_coords( conn4, num_nodes, &( verts[0][0] ) );
297  if( moab::MB_SUCCESS != rval )
298  {
299  std::cout << "can't get coords for quad \n";
300  delete mb;
301  return;
302  }
303 
304  specQuad.set_vertices( verts );
305  specQuad.compute_gl_positions();
306  // do something with the gl positions, project them on a sphere, and create another mesh?
307  if( rit == ents.begin() )
308  {
309  std::cout << " gl points for first element: \n";
310  int size;
311  double* xyz[3];
312  specQuad.get_gl_points( xyz[0], xyz[1], xyz[2], size );
313  for( int i = 0; i < size; i++ )
314  std::cout << xyz[0][i] << " " << xyz[1][i] << " " << xyz[2][i] << "\n";
315  }
316 
317  // project them on a sphere, and create another mesh with it?
318  }
319  std::cout << "success...\n";
320 
321  delete mb;
322 }
324 {
325  // example from one coupler test, run like this
326  // ./mbcoupler_test -meshes sphere_16p.h5m mpas_p8.h5m -itag vertex_field -meth 4 -outfile
327  // dd.h5m method 4 is spherical
328  double positions[] = { -0.88882388032987436, -0.069951956448441419, 0.45287838714646161, -0.88226455385461389,
329  -0.13973697758043971, 0.4495362433757738, -0.84497006020160348, -0.13383011007602069,
330  0.51779831884618843, -0.85072691325794214, -0.066953660115039074, 0.52132612293631853 };
331  CartVect x( -0.85408569769999998, -0.12391301439999999, 0.50515659540000002 );
332  std::vector< CartVect > vertices;
333  for( int i = 0; i < 4; i++ )
334  vertices.push_back( CartVect( positions + 3 * i ) );
335 
336  moab::Element::SphericalQuad squad( vertices );
337  double tol( 0.0001 );
338  if( squad.inside_box( x, tol ) )
339  {
340  CartVect nat_par = squad.ievaluate( x, 0.000001 );
341  std::cout << nat_par << "\n";
342  }
343  std::cout << "success...\n";
344 }
345 
347 {
348  double positions[] = { 0, 0, 0, 2, 0, 0, 0, 3, 0 };
349  CartVect x( 1, 0.5, 0 );
350  std::vector< CartVect > vertices;
351  for( int i = 0; i < 3; i++ )
352  vertices.push_back( CartVect( positions + 3 * i ) );
353 
354  moab::Element::LinearTri tri( vertices );
355  double tol( 0.0001 );
356  if( tri.inside_box( x, tol ) )
357  {
358  CartVect nat_par = tri.ievaluate( x );
359  std::cout << x << " :" << nat_par << "\n";
360  }
361 
362  x = CartVect( 0, 2, 0 );
363  if( tri.inside_box( x, tol ) )
364  {
365  CartVect nat_par = tri.ievaluate( x );
366  std::cout << x << " :" << nat_par << "\n";
367  }
368 
369  x = CartVect( 1, 0, 0.5 );
370  if( tri.inside_box( x, tol ) )
371  {
372  CartVect nat_par = tri.ievaluate( x );
373  std::cout << x << " :" << nat_par << "\n";
374  }
375 
376  double positions2[] = { -0.866026, -0.500001, 0., 0.866026, -0.500001, 0., 0.000000, 100.000000, 0. };
377  x = CartVect( 0, 0, 0 );
378  std::vector< CartVect > vertices2;
379  for( int i = 0; i < 3; i++ )
380  vertices2.push_back( CartVect( positions2 + 3 * i ) );
381 
382  moab::Element::LinearTri tri2( vertices2 );
383 
384  if( tri2.inside_box( x, tol ) )
385  {
386  CartVect nat_par = tri2.ievaluate( x );
387  std::cout << x << " :" << nat_par << "\n";
388  }
389 
390  std::cout << "vertices2 " << vertices2[0] << " " << vertices2[1] << " " << vertices2[2] << "\n";
391 
392  x = CartVect( -0.866026, -0.500001, 0. );
393  std::cout << x << " :" << tri2.ievaluate( x ) << "\n";
394 
395  x = CartVect( +0.866026, -0.500001, 0. );
396  std::cout << x << " :" << tri2.ievaluate( x ) << "\n";
397  x = CartVect( 0.000000, 100.000000, 0. );
398  std::cout << x << " :" << tri2.ievaluate( x ) << "\n";
399 
400  std::cout << "success...\n";
401 }
402 
404 {
405  // example from one coupler test, run like this
406  // ./mbcoupler_test -meshes tri_fl_8p.h5m mpas_p8.h5m -itag vertex_field -meth 4 -outfile
407  // oo.h5m -eps 1.e-9 method 4 is spherical
408  double positions[] = { -0.86339258282987197, -0.17004443185241255, 0.47501383044112816,
409  -0.80777478326268271, -0.15172299908552511, 0.5696314870803928,
410  -0.8655618847392077, -0.061613422011313854, 0.49699739427361828 };
411  CartVect x( -0.85408569769999998, -0.12391301439999999, 0.50515659540000002 );
412  std::vector< CartVect > vertices;
413  for( int i = 0; i < 3; i++ )
414  vertices.push_back( CartVect( positions + 3 * i ) );
415 
416  moab::Element::SphericalTri sphtri( vertices );
417  double tol( 0.000001 );
418  if( sphtri.inside_box( x, tol ) )
419  {
420  CartVect nat_par = sphtri.ievaluate( x, 0.000001 );
421  std::cout << nat_par << "\n";
422  }
423  std::cout << "success...\n";
424 }