22 if( MB_SUCCESS != rval ) return rval
24 const double ROOT2 = 1.4142135623730951;
64 const double coords[] = { 1, -1, -1, 1, 1, -1, -1, 1, -1, -1, -1, -1, 1, -1,
65 1, 1, 1, 1, -1, 1, 1, -1, -1, 1, 0, 0, 0 };
66 const int connectivity[] = {
75 const unsigned tris_per_surf[] = { 2, 2, 2, 2, 2, 4 };
78 const unsigned num_verts =
sizeof( coords ) / ( 3 *
sizeof(
double ) );
79 const unsigned num_tris =
sizeof( connectivity ) / ( 3 *
sizeof(
int ) );
80 const unsigned num_surfs =
sizeof( tris_per_surf ) /
sizeof(
unsigned );
81 EntityHandle verts[num_verts], tris[num_tris], surfs[num_surfs];
82 for(
unsigned i = 0; i < num_verts; ++i )
84 rval =
moab->create_vertex( coords + 3 * i, verts[i] );
CHKERR;
86 for(
unsigned i = 0; i < num_tris; ++i )
88 const EntityHandle conn[] = { verts[connectivity[3 * i]], verts[connectivity[3 * i + 1]],
89 verts[connectivity[3 * i + 2]] };
95 for(
unsigned i = 0; i < num_surfs; ++i )
98 rval =
moab->add_entities( surfs[i], tri_iter, tris_per_surf[i] );
CHKERR;
99 tri_iter += tris_per_surf[i];
107 std::vector< int > dims( num_surfs, 2 );
109 std::vector< int > ids( num_surfs );
110 for(
size_t i = 0; i < ids.size(); ++i )
116 for(
unsigned i = 0; i < num_surfs; ++i )
118 rval =
moab->add_parent_child( volume, surfs[i] );
CHKERR;
121 std::vector< EntityHandle > senses( 2 * num_surfs, 0 );
122 for(
size_t i = 0; i < senses.size(); i += 2 )
124 rval =
moab->tag_set_data( sense_tag, surfs, num_surfs, &senses[0] );
CHKERR;
131 rval =
moab->write_mesh( output_file_name );
CHKERR;
143 const double coords[] = { 1, -1, -1, 1, 1, -1, 0, 1, -1, 0, -1, -1, 1, -1, 1, 1, 1, 1, 0, 1, 1, 0, -1, 1,
145 0.01, -1, -1, 0.01, 1, -1, -1, 1, -1, -1, -1, -1, 0.01, -1, 1, 0.01, 1, 1, -1, 1, 1, -1,
147 const int connectivity[] = { 0, 3, 1, 3, 2, 1,
155 const unsigned tris_per_surf = 2;
156 const unsigned num_cubes = 2;
157 const unsigned num_verts =
sizeof( coords ) / ( 3 *
sizeof(
double ) );
158 const unsigned num_tris = num_cubes *
sizeof( connectivity ) / ( 3 *
sizeof(
int ) );
159 const unsigned num_surfs = num_tris / tris_per_surf;
160 EntityHandle verts[num_verts], tris[num_tris], surfs[num_surfs];
162 for(
unsigned i = 0; i < num_verts; ++i )
164 rval =
moab->create_vertex( coords + 3 * i, verts[i] );
CHKERR;
167 for(
unsigned i = 0; i < num_tris / 2; ++i )
169 const EntityHandle conn[] = { verts[connectivity[3 * i]], verts[connectivity[3 * i + 1]],
170 verts[connectivity[3 * i + 2]] };
174 for(
unsigned i = 0; i < num_tris / 2; ++i )
176 const EntityHandle conn[] = { verts[8 + connectivity[3 * i]], verts[8 + connectivity[3 * i + 1]],
177 verts[8 + connectivity[3 * i + 2]] };
178 rval =
moab->create_element(
MBTRI, conn, 3, tris[num_tris / 2 + i] );
CHKERR;
183 for(
unsigned i = 0; i < num_surfs; ++i )
186 rval =
moab->add_entities( surfs[i], tri_iter, tris_per_surf );
CHKERR;
187 tri_iter += tris_per_surf;
195 std::vector< int > dims( num_surfs, 2 );
197 std::vector< int > ids( num_surfs );
198 for(
size_t i = 0; i < ids.size(); ++i )
204 for(
unsigned i = 0; i < num_surfs; ++i )
206 rval =
moab->add_parent_child( volume, surfs[i] );
CHKERR;
209 std::vector< EntityHandle > senses( 2 * num_surfs, 0 );
210 for(
size_t i = 0; i < senses.size(); i += 2 )
212 rval =
moab->tag_set_data( sense_tag, surfs, num_surfs, &senses[0] );
CHKERR;
219 rval =
moab->write_mesh( output_file_name );
CHKERR;
224 #define RUN_TEST( A ) \
227 std::cout << #A << "... " << std::endl; \
228 if( MB_SUCCESS != A( gqt ) ) \
234 int main(
int argc,
char* argv[] )
237 const char*
filename =
"test_geom.h5m";
240 int fail = MPI_Init( &argc, &argv );
281 std::cout <<
"-------------------------------------" << std::endl;
282 std::cout <<
"Re-running tests with MBI constructor" << std::endl;
283 std::cout <<
"-------------------------------------" << std::endl;
314 fail = MPI_Finalize();
336 double overlap_thickness = 0.1;
352 double overlap_thickness = 3.0;
371 const int two = 2, three = 3;
372 const void* ptrs[] = { &two, &three };
376 if( vols.
size() != 2 )
378 std::cerr <<
"ERROR: Expected 2 volumes in input, found " << vols.
size() << std::endl;
381 if( surfs.
size() != 6 )
383 std::cerr <<
"ERROR: Expected 6 surfaces in input, found " << surfs.
size() << std::endl;
393 std::cerr <<
"ERROR: Expected 1 for surface sense, got " << sense << std::endl;
408 const int two = 2, three = 3;
409 const void* ptrs[] = { &two, &three };
413 if( vols.
size() != 2 )
415 std::cerr <<
"ERROR: Expected 2 volumes in input, found " << vols.
size() << std::endl;
418 if( surfs.
size() != 12 )
420 std::cerr <<
"ERROR: Expected 12 surfaces in input, found " << surfs.
size() << std::endl;
430 std::cerr <<
"ERROR: Expected 1 for surface sense, got " << sense << std::endl;
445 const void* ptr = &three;
448 if( vols.
size() != 2 )
450 std::cerr <<
"ERROR: Expected 2 volumes in input, found " << vols.
size() << std::endl;
457 const double vol = 2 * 2 * 2 - 1 * 4. / 3;
460 if( fabs( result - vol ) > 10 * std::numeric_limits< double >::epsilon() )
462 std::cerr <<
"ERROR: Expected " << vol <<
" as measure of volume, got " << result << std::endl;
476 const void* ptr = &three;
479 if( vols.
size() != 2 )
481 std::cerr <<
"ERROR: Expected 2 volumes in input, found " << vols.
size() << std::endl;
487 const double vol = ( 1 + 1.01 ) * 2 * 2;
490 if( fabs( result - vol ) > 2 * std::numeric_limits< double >::epsilon() )
492 std::cerr <<
"ERROR: Expected " << vol <<
" as measure of volume, got " << result << std::endl;
507 const void* ptr = &two;
510 if( surfs.
size() != 6 )
512 std::cerr <<
"ERROR: Expected 6 surfaces in input, found " << surfs.
size() << std::endl;
522 for(
unsigned i = 0; i < 6; ++i, ++iter )
524 double expected = 4.0;
525 if( ids[i] == 6 ) expected *=
ROOT2;
530 if( fabs( result - expected ) > std::numeric_limits< double >::epsilon() )
532 std::cerr <<
"ERROR: Expected area of surface " << ids[i] <<
" to be " << expected <<
". Got " << result
549 const void* ptr = &two;
552 const unsigned num_surfs = 12;
553 if( surfs.
size() != num_surfs )
555 std::cerr <<
"ERROR: Expected " << num_surfs <<
" surfaces in input, found " << surfs.
size() << std::endl;
562 const double x_area = 2 * 2;
563 const double yz_area0 = 2 * 1;
564 const double yz_area1 = 2 * 1.01;
566 for(
unsigned i = 0; i < num_surfs; ++i, ++iter )
568 double expected, result;
569 if( 0 == i || 2 == i || 4 == i || 5 == i )
571 else if( 1 == i || 3 == i || 7 == i || 9 == i )
573 else if( 6 == i || 8 == i || 10 == i || 11 == i )
577 if( fabs( result - expected ) > std::numeric_limits< double >::epsilon() )
579 std::cerr <<
"ERROR: Expected area of surface " << ids[i] <<
" to be " << expected <<
". Got " << result
603 { 1, { 0.0, 0.0, -1. }, { -1.0 /
ROOT2, 0.0, 1.0 /
ROOT2 }, 4,
ROOT2 },
609 { 1, { 0.5, 0.5, -1. }, { 0.0, 0.0, 1.0 }, 6, 1.5 },
611 { 2, { 1.0, 0.0, 0.5 }, { -1.0, 0.0, 0.0 }, 6, 0.5 },
613 { 2, { 1.0, 0.0, 0.0 }, { -1.0, 0.0, 0.0 }, 4, 2.0 },
615 { 1, { 0.0, 0.0, -1. }, { 0.0, 0.0, 1.0 }, 6, 1.0 },
617 { 2, { 1.0, 0.0, 0.5 }, { -1.0 /
ROOT2, 1.0 /
ROOT2, 0.0 }, 3,
ROOT2 } };
624 const int two = 2, three = 3;
625 const void* ptrs[] = { &two, &three };
629 if( vols.
size() != 2 )
631 std::cerr <<
"ERROR: Expected 2 volumes in input, found " << vols.
size() << std::endl;
634 if( surfs.
size() != 6 )
636 std::cerr <<
"ERROR: Expected 6 surfaces in input, found " << surfs.
size() << std::endl;
643 std::copy( surfs.
begin(), surfs.
end(), surf );
645 const int num_test =
sizeof( tests ) /
sizeof( tests[0] );
646 for(
int i = 0; i < num_test; ++i )
648 int* ptr = std::find( ids, ids + 6, tests[i].
prev_surf );
652 std::cerr <<
"Surface " << tests[i].
prev_surf <<
" not found." << std::endl;
657 ptr = std::find( ids, ids + 6, tests[i].
hit_surf );
661 std::cerr <<
"Surface " << tests[i].
hit_surf <<
" not found." << std::endl;
675 int id = idx > 5 ? 0 : ids[idx];
677 std::cerr <<
"Rayfire test failed for " << std::endl
678 <<
"\t ray from (" << tests[i].
origin[0] <<
", " << tests[i].
origin[1] <<
", "
680 <<
", " << tests[i].
direction[2] <<
"]" << std::endl
681 <<
"\t Beginning on surface " << tests[i].
prev_surf << std::endl
682 <<
"\t Expected to hit surface " << tests[i].
hit_surf <<
" after " << tests[i].
distance
683 <<
" units." << std::endl
684 <<
"\t Actually hit surface " <<
id <<
" after " << dist <<
" units." << std::endl;
690 std::vector< std::pair< int, GeomQueryTool::RayHistory* > > boundary_tests;
691 boundary_tests.push_back( std::make_pair( 1, &history ) );
692 boundary_tests.push_back( std::make_pair( 0, &history ) );
696 for(
unsigned int bt = 0; bt < boundary_tests.size(); ++bt )
699 int expected = boundary_tests[bt].first;
705 if( expected == 1 ) uvw = -uvw;
707 int boundary_result = -1;
711 if( boundary_result != expected )
713 std::cerr <<
"DagMC::test_volume_boundary failed (" << ( ( expected == 0 ) ?
"+" :
"-" ) <<
" dir,"
714 << ( ( h ) ?
"+" :
"-" ) <<
" history, i=" << i <<
")" << std::endl;
741 { 4, { 1.0, 0.0, 0.0 }, { 1.0, 0.0, 0.0 }, 2, 0.0 },
742 { 2, { 1.0, 0.0, 0.0 }, { -1.0, 0.0, 0.0 }, 4, 1.0 },
744 { 4, { 0.5, 0.0, 0.0 }, { 1.0, 0.0, 0.0 }, 2, 0.5 },
745 { 2, { 0.5, 0.0, 0.0 }, { -1.0, 0.0, 0.0 }, 4, 0.5 },
747 { 4, { 0.01, 0.0, 0.0 }, { 1.0, 0.0, 0.0 }, 8, 0.0 },
748 { 2, { 0.01, 0.0, 0.0 }, { -1.0, 0.0, 0.0 }, 4, 0.01 },
750 { 10, { 0.005, 0.0, 0.0 }, { 1.0, 0.0, 0.0 }, 8, 0.005 },
751 { 2, { 0.005, 0.0, 0.0 }, { -1.0, 0.0, 0.0 }, 4, 0.005 },
753 { 10, { 0.0, 0.0, 0.0 }, { 1.0, 0.0, 0.0 }, 8, 0.01 },
754 { 8, { 0.0, 0.0, 0.0 }, { -1.0, 0.0, 0.0 }, 4, 0.0 },
756 { 10, { -0.5, 0.0, 0.0 }, { 1.0, 0.0, 0.0 }, 8, 0.51 },
757 { 8, { -0.5, 0.0, 0.0 }, { -1.0, 0.0, 0.0 }, 10, 0.5 },
759 { 10, { -1.0, 0.0, 0.0 }, { 1.0, 0.0, 0.0 }, 8, 1.01 },
760 { 8, { -1.0, 0.0, 0.0 }, { -1.0, 0.0, 0.0 }, 10, 0.0 } };
767 const int two = 2, three = 3;
768 const void* ptrs[] = { &two, &three };
772 if( vols.
size() != 2 )
774 std::cerr <<
"ERROR: Expected 2 volumes in input, found " << vols.
size() << std::endl;
777 const unsigned num_surf = 12;
778 if( surfs.
size() != num_surf )
780 std::cerr <<
"ERROR: Expected " << num_surf <<
" surfaces in input, found " << surfs.
size() << std::endl;
787 std::copy( surfs.
begin(), surfs.
end(), surf );
789 const int num_test =
sizeof( tests ) /
sizeof( tests[0] );
790 for(
int i = 0; i < num_test; ++i )
792 int* ptr = std::find( ids, ids + num_surf, tests[i].
prev_surf );
793 unsigned idx = ptr - ids;
794 if( idx >= num_surf )
796 std::cerr <<
"Surface " << tests[i].
prev_surf <<
" not found." << std::endl;
801 ptr = std::find( ids, ids + num_surf, tests[i].
hit_surf );
803 if( idx >= num_surf )
805 std::cerr <<
"Surface " << tests[i].
hit_surf <<
" not found." << std::endl;
818 int id = idx > num_surf - 1 ? 0 : ids[idx];
820 std::cerr <<
"Rayfire test failed for " << std::endl
821 <<
"\t ray from (" << tests[i].
origin[0] <<
", " << tests[i].
origin[1] <<
", "
823 <<
", " << tests[i].
direction[2] <<
"]" << std::endl
824 <<
"\t Beginning on surface " << tests[i].
prev_surf << std::endl
825 <<
"\t Expected to hit surface " << tests[i].
hit_surf <<
" after " << tests[i].
distance
826 <<
" units." << std::endl
827 <<
"\t Actually hit surface " <<
id <<
" after " << dist <<
" units." << std::endl;
844 const char*
const NAME_ARR[] = {
"Boundary",
"Outside",
"Inside" };
845 const char*
const* names = NAME_ARR + 1;
846 const int INSIDE = 1, OUTSIDE = 0, BOUNDARY = -1;
847 const struct PointInVol tests[] = { { { 0.0, 0.0, 0.5 }, OUTSIDE, { 0.0, 0.0, 0.0 } },
848 { { 0.0, 0.0, -0.5 }, INSIDE, { 0.0, 0.0, 0.0 } },
849 { { 0.7, 0.0, 0.0 }, INSIDE, { 0.0, 0.0, 0.0 } },
850 { { -0.7, 0.0, 0.0 }, INSIDE, { 0.0, 0.0, 0.0 } },
851 { { 0.0, -0.7, 0.0 }, INSIDE, { 0.0, 0.0, 0.0 } },
852 { { 0.0, -0.7, 0.0 }, INSIDE, { 0.0, 0.0, 0.0 } },
853 { { 1.1, 1.1, 1.1 }, OUTSIDE, { 0.0, 0.0, 0.0 } },
854 { { -1.1, 1.1, 1.1 }, OUTSIDE, { 0.0, 0.0, 0.0 } },
855 { { -1.1, -1.1, 1.1 }, OUTSIDE, { 0.0, 0.0, 0.0 } },
856 { { 1.1, -1.1, 1.1 }, OUTSIDE, { 0.0, 0.0, 0.0 } },
857 { { 1.1, 1.1, -1.1 }, OUTSIDE, { 0.0, 0.0, 0.0 } },
858 { { -1.1, 1.1, -1.1 }, OUTSIDE, { 0.0, 0.0, 0.0 } },
859 { { -1.1, -1.1, -1.1 }, OUTSIDE, { 0.0, 0.0, 0.0 } },
860 { { 1.1, -1.1, -1.1 }, OUTSIDE, { 0.0, 0.0, 0.0 } },
862 { { 0.5, 0.0, 0.0 }, INSIDE, { -1.0, 0.0, 0.0 } },
863 { { 0.5, 0.0, 0.0 }, INSIDE, { 1.0, 0.0, 0.0 } },
864 { { 0.0, 0.0, 2.0 }, OUTSIDE, { 0.0, 0.0, -1.0 } },
865 { { 0.5, 0.0, -0.5 }, INSIDE, { -1.0, 0.0, 0.0 } },
866 { { 0.5, -0.5, -2.0 }, OUTSIDE, { 0.0, 0.0, 1.0 } } };
874 const int num_test =
sizeof( tests ) /
sizeof( tests[0] );
883 const void* ptr = &three;
885 if( vols.
size() != 2 )
887 std::cerr <<
"ERROR: Expected 2 volumes in input, found " << vols.
size() << std::endl;
892 for(
int i = 0; i < num_test; ++i )
898 std::cerr <<
"ERROR testing point_in_volume[" << i <<
"]:" << std::endl
899 <<
"\tExpected " << names[tests[i].
result] <<
" for (" << tests[i].
coords[0] <<
", "
900 << tests[i].
coords[1] <<
", " << tests[i].
coords[2] <<
"). Got " << names[
result] << std::endl;
905 if( tests[i].
result == BOUNDARY )
continue;
911 std::cerr <<
"ERROR testing point_in_volume_slow[" << i <<
"]:" << std::endl
912 <<
"\tExpected " << names[tests[i].
result] <<
" for (" << tests[i].
coords[0] <<
", "
913 << tests[i].
coords[1] <<
", " << tests[i].
coords[2] <<
"). Got " << names[
result] << std::endl;
928 const void* ptr = &three;
931 if( vols.
size() != 2 )
933 std::cerr <<
"ERROR: Expected 2 volumes in input, found " << vols.
size() << std::endl;
941 if(
result != 9.0 )
return MB_FAILURE;
946 if(
result != 9.0 )
return MB_FAILURE;
947 if( surface == 1 )
return MB_FAILURE;
954 const char*
const NAME_ARR[] = {
"Boundary",
"Outside",
"Inside" };
955 const char*
const* names = NAME_ARR + 1;
956 const int INSIDE = 1, OUTSIDE = 0, BOUNDARY = -1;
957 const struct PointInVol tests[] = { { { 0.5, 0.0, 0.5 }, INSIDE, { 0.0, 0.0, 0.0 } },
958 { { 0.5, 0.0, -0.5 }, INSIDE, { 0.0, 0.0, 0.0 } },
959 { { 0.5, 0.0, 0.0 }, INSIDE, { 0.0, 0.0, 0.0 } },
960 { { -0.5, 0.0, 0.0 }, INSIDE, { 0.0, 0.0, 0.0 } },
961 { { 0.5, 0.5, 0.0 }, INSIDE, { 0.0, 0.0, 0.0 } },
962 { { 0.5, -0.5, 0.0 }, INSIDE, { 0.0, 0.0, 0.0 } },
963 { { 1.1, 1.1, 1.1 }, OUTSIDE, { 0.0, 0.0, 0.0 } },
964 { { -1.1, 1.1, 1.1 }, OUTSIDE, { 0.0, 0.0, 0.0 } },
965 { { -1.1, -1.1, 1.1 }, OUTSIDE, { 0.0, 0.0, 0.0 } },
966 { { 1.1, -1.1, 1.1 }, OUTSIDE, { 0.0, 0.0, 0.0 } },
967 { { 1.1, 1.1, -1.1 }, OUTSIDE, { 0.0, 0.0, 0.0 } },
968 { { -1.1, 1.1, -1.1 }, OUTSIDE, { 0.0, 0.0, 0.0 } },
969 { { -1.1, -1.1, -1.1 }, OUTSIDE, { 0.0, 0.0, 0.0 } },
970 { { 1.1, -1.1, -1.1 }, OUTSIDE, { 0.0, 0.0, 0.0 } },
972 { { 0.5, 0.0, 0.0 }, INSIDE, { -1.0, 0.0, 0.0 } },
973 { { 0.5, 0.0, 0.0 }, INSIDE, { 1.0, 0.0, 0.0 } },
974 { { 0.0, 0.0, 2.0 }, OUTSIDE, { 0.0, 0.0, -1.0 } },
975 { { 0.5, 0.0, -0.5 }, INSIDE, { -1.0, 0.0, 0.0 } },
976 { { 0.5, -0.5, -2.0 }, OUTSIDE, { 0.0, 0.0, 1.0 } },
978 { { 0.005, 0.0, 0.0 }, INSIDE, { -1.0, 0.0, 0.0 } },
979 { { 0.005, 0.0, 0.0 }, INSIDE, { 1.0, 0.0, 0.0 } },
980 { { 0.005, 0.0, 2.0 }, OUTSIDE, { 0.0, 0.0, -1.0 } },
981 { { 0.005, 0.0, -0.5 }, INSIDE, { -1.0, 0.0, 0.0 } },
982 { { 0.005, -0.5, -2.0 }, OUTSIDE, { 0.0, 0.0, 1.0 } } };
984 const int num_test =
sizeof( tests ) /
sizeof( tests[0] );
992 const void* ptr = &three;
994 if( vols.
size() != 2 )
996 std::cerr <<
"ERROR: Expected 2 volumes in input, found " << vols.
size() << std::endl;
1000 for(
int i = 0; i < num_test; ++i )
1006 std::cerr <<
"ERROR testing point_in_volume[" << i <<
"]:" << std::endl
1007 <<
"\tExpected " << names[tests[i].
result] <<
" for (" << tests[i].
coords[0] <<
", "
1008 << tests[i].
coords[1] <<
", " << tests[i].
coords[2] <<
"). Got " << names[
result] << std::endl;
1013 if( tests[i].
result == BOUNDARY )
continue;
1019 std::cerr <<
"ERROR testing point_in_volume_slow[" << i <<
"]:" << std::endl
1020 <<
"\tExpected " << names[tests[i].
result] <<
" for (" << tests[i].
coords[0] <<
", "
1021 << tests[i].
coords[1] <<
", " << tests[i].
coords[2] <<
"). Got " << names[
result] << std::endl;
1042 Range surfs, explicit_vols;
1043 const int two = 2, three = 3;
1044 const void* ptrs[] = { &two, &three };
1050 if( explicit_vols.
size() != 2 )
1052 std::cerr <<
"ERROR: Expected 2 explicit volumes in input, found " << explicit_vols.
size() << std::endl;
1055 const unsigned num_surf = 12;
1056 if( surfs.
size() != num_surf )
1058 std::cerr <<
"ERROR: Expected " << num_surf <<
" surfaces in input, found " << surfs.
size() << std::endl;
1064 double point[] = { -0.9, 0, 0 };
1065 const double dir[] = { 1, 0, 0 };
1068 const int INSIDE = 1;
1072 std::cerr <<
"ERROR: particle not inside explicit volume" << std::endl;
1081 if( next_surf != surfs[7] || fabs( dist - 0.91 ) > 1e-6 )
1083 std::cerr <<
"ERROR: failed on advance 1" << std::endl;
1086 for(
unsigned i = 0; i < 3; i++ )
1087 point[i] += dist *
dir[i];
1096 if( next_surf != surfs[3] || fabs( dist - 0.0 ) > 1e-6 )
1098 std::cerr <<
"ERROR: failed on advance 2" << std::endl;
1101 for(
unsigned i = 0; i < 3; i++ )
1102 point[i] += dist *
dir[i];
1110 if( next_surf != surfs[1] || fabs( dist - 0.99 ) > 1e-6 )
1112 std::cerr <<
"ERROR: failed on advance 3" << std::endl;
1115 for(
unsigned i = 0; i < 3; i++ )
1116 point[i] += dist *
dir[i];