 |
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 "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 {
36 moab::Element::LinearTet tet;
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 }
84 catch( Element::Map::EvaluationError )
85 {
86 // nothing
87 }
88 }
89
90 } // test_hex()
91
92 void test_spectral_hex()
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
248 void test_spectral_quad()
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 }
323 void test_spherical_quad()
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
346 void test_linear_tri()
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
403 void test_spherical_tri()
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 }
1.9.1.