[MOAB-dev] r1193 - in MOAB/trunk: . tools/dagmc
kraftche at mcs.anl.gov
kraftche at mcs.anl.gov
Tue Jul 10 17:59:37 CDT 2007
Author: kraftche
Date: 2007-07-10 17:59:36 -0500 (Tue, 10 Jul 2007)
New Revision: 1193
Added:
MOAB/trunk/tools/dagmc/test_geom.cc
Modified:
MOAB/trunk/configure.in
MOAB/trunk/tools/dagmc/DagMC.cpp
MOAB/trunk/tools/dagmc/DagMC.hpp
MOAB/trunk/tools/dagmc/Makefile.am
Log:
o Fix building of DagMC w/out CGM
o Fix function w/ no return val if building w/out CGM
o Make DagMC construct OBB tree if one doesn't already exist,
regardless of input file type.
o Fix broken point-in-volume test: clear static lists in each call
o Use smarter algorithm for ponit-in-volume test to avoid
slow solid-angle algorithm when input location is closest
to a vertex.
o Add unit tests for geometry functions in DagMC
NOTE: Unit test for DagMC::fire_ray fails due to a known issue.
Modified: MOAB/trunk/configure.in
===================================================================
--- MOAB/trunk/configure.in 2007-07-10 17:31:22 UTC (rev 1192)
+++ MOAB/trunk/configure.in 2007-07-10 22:59:36 UTC (rev 1193)
@@ -664,8 +664,8 @@
AC_SUBST(CGM_DIR)
AC_SUBST(CGM_CONFIG_OPTIONS)
fi
+AM_CONDITIONAL( HAVE_CGM, [test "x$CGM_MISSING" = "xno"] )
-
################################################################################
# iMesh Babel
################################################################################
@@ -758,9 +758,6 @@
AC_MSG_WARN([Count not find VTK. Build will FAIL!!! Try \"--with-VTK=DIR\"])
fi
fi
-if test "x$CGM_MISSING" = "xyes"; then
- AC_MSG_WARN([cgm2moab/dagmc build will fail! CGM not found!])
-fi
if test "x$HAVE_HDF5" = "xno"; then
AC_MSG_WARN([
Modified: MOAB/trunk/tools/dagmc/DagMC.cpp
===================================================================
--- MOAB/trunk/tools/dagmc/DagMC.cpp 2007-07-10 17:31:22 UTC (rev 1192)
+++ MOAB/trunk/tools/dagmc/DagMC.cpp 2007-07-10 22:59:36 UTC (rev 1193)
@@ -366,10 +366,8 @@
// point is inside or outside of the volume. If the point is closest to
// an edge joining two facets, then the point is inside the volume if
// that edge is at a concavity wrt the volume. If the point is closest
-// to a vertex in the facetting, then if the normals of the facets are
-// sufficiently close to each other as simple inside/outside test is
-// used for each facet. If not, the (very) slow spherical-area method
-// is used.
+// to a vertex in the facetting, it relies on the fact that the closest
+// vertex cannot be a saddle: it must be a strict concavity or convexity.
MBErrorCode DagMC::point_in_volume( MBEntityHandle volume,
double x, double y, double z,
int& result )
@@ -382,6 +380,8 @@
// Get closest triangles in volume
static std::vector<MBEntityHandle> tris, surfs;
+ tris.clear();
+ surfs.clear();
const MBCartVect point(x,y,z);
MBErrorCode rval = obbTree.closest_to_location( point.array(), root, epsilon, tris, &surfs );
if (MB_SUCCESS != rval) return rval;
@@ -508,45 +508,79 @@
return MB_SUCCESS;
}
- // Otherwise test triangles about each vertex that the point was closest to
-
- // First get unique list of vertices
- std::vector<MBEntityHandle> vertices( tris.size() ), uvertices;
- for (unsigned i = 0; i < tris.size(); ++i) {
- assert( (unsigned)(topo[i]) < 3);
- moab_instance()->get_connectivity( tris[1], conn, len, true );
- vertices.push_back( conn[topo[i]] );
+ // Otherwise pick a vertex that the point was closest to remove
+ // any triangles from lists that are not adjacent to that vertex.
+ assert( topo[0] < 3 ); // if here, must be closest to vertex
+ rval = moab_instance()->get_connectivity( tris[0], conn, len );
+ if (MB_SUCCESS != rval || len != 3)
+ return MB_FAILURE;
+ const MBEntityHandle closest_vert = conn[topo[0]];
+ unsigned w = 1;
+ for (unsigned r = 1; r < tris.size(); ++r) {
+ assert( topo[r] < 3 ); // if here, must be closest to vertex
+ moab_instance()->get_connectivity( tris[r], conn, len );
+ if (conn[topo[r]] == closest_vert) {
+ tris[w] = tris[r];
+ surfs[w] = surfs[r];
+ normals[w] = normals[r];
+ ++w;
+ }
}
- uvertices = vertices;
- std::sort( uvertices.begin(), uvertices.end() );
- uvertices.erase( std::unique( uvertices.begin(), uvertices.end() ), uvertices.end() );
+ tris.resize(w);
+ surfs.resize(w);
+ normals.resize(w);
- // Now for each vertex, test if we are inside or outside all the triangles
- for (unsigned i = 0; i < uvertices.size(); ++i) {
- bool inside_all = true, outside_all = true;
+ // use algorithm from:
+ // (referance)
+ // to determine inside vs. outside.
+
+ // first find single triangle for which all other triangles are to
+ // one side.
+ bool some_above, some_below;
+ MBCartVect vertcoords, closestcoords;
+ moab_instance()->get_coords( &closest_vert, 1, closestcoords.array() );
+ for (unsigned i = 0; i < tris.size(); ++i) {
+ some_above = some_below = false;
for (unsigned j = 0; j < tris.size(); ++j) {
- if (vertices[j] != uvertices[i])
+ if (j == i)
continue;
-
- const double dot = normals[j] % (closest[j] - point);
- if (dot < -std::numeric_limits<double>::epsilon())
- inside_all = false;
- if (dot > std::numeric_limits<double>::epsilon())
- outside_all = false;
+
+ moab_instance()->get_connectivity( tris[j], conn, len );
+ for (int k = 0; k < len; ++k) {
+ moab_instance()->get_coords( conn+k, 1, vertcoords.array() );
+ const double dot = normals[i] % (closestcoords - vertcoords);
+ if (dot * dot / (normals[i] % normals[i]) > epsilon * epsilon) {
+ if (dot < 0.0)
+ some_above = true;
+ else
+ some_below = true;
+ }
+ }
}
- if (inside_all) {
- result = 1;
+
+ // All triangles are roughly co-planar: if we're inside of
+ // one then we're inside of all.
+ if (!some_above && !some_below) {
+ result = (normals[i] % (closestcoords - point)) > 0.0;
return MB_SUCCESS;
}
- else if (outside_all) {
- result = 0;
+ // All other triangles to one side of this triangle:
+ else if (some_above != some_below) {
+ // If all other triangles are above this one (some_above == true),
+ // then the vertex is a concavity so the input point is inside.
+ // If all other triangles are below this one (some_above == false),
+ // then the vertex is a convexity and the input point must be outside.
+ result = some_above;
return MB_SUCCESS;
}
}
- // If this far, then closest to only vertices in the facetting
- // and cannot determine inside/outside at any vertex.
- // Use the (very) slow method
+ // If we got here, then something's wrong somewhere.
+ // We appear to be closest to a "saddle" vertex in the
+ // triangulation. That shoudn't be possible (must be
+ // closer to at least one of the adjacent triagles than
+ // to the saddle vertex.)
+ assert(false /*shouldn't be here*/);
return point_in_volume_slow( volume, x, y, z, result );
}
@@ -787,7 +821,8 @@
return rval;
// Build OBB trees for everything, but only if we only read geometry
- if (is_geom) {
+ // Changed to build obb tree if tree does not already exist. -- JK
+ if (!have_obb_tree()) {
rval = build_obbs(surfs, vols);
if (MB_SUCCESS != rval) {
std::cerr << "Failed to build obb." << std::endl;
@@ -815,6 +850,15 @@
return MB_SUCCESS;
}
+bool DagMC::have_obb_tree()
+{
+ MBRange entities;
+ MBErrorCode rval = mbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET,
+ &obbTag, 0, 1,
+ entities );
+ return MB_SUCCESS == rval && !entities.empty();
+}
+
MBEntityHandle DagMC::entity_by_id( int dimension, int id )
{
assert(0 <= dimension && 3 >= dimension);
@@ -1319,6 +1363,8 @@
return MB_SUCCESS;
#endif
+#else
+ return MB_FAILURE;
#endif
}
Modified: MOAB/trunk/tools/dagmc/DagMC.hpp
===================================================================
--- MOAB/trunk/tools/dagmc/DagMC.hpp 2007-07-10 17:31:22 UTC (rev 1192)
+++ MOAB/trunk/tools/dagmc/DagMC.hpp 2007-07-10 22:59:36 UTC (rev 1193)
@@ -135,6 +135,8 @@
double &len);
private:
+ bool have_obb_tree();
+
MBTag get_tag( const char* name, int size, MBTagType store, MBDataType type,
bool create_if_missing = true);
Modified: MOAB/trunk/tools/dagmc/Makefile.am
===================================================================
--- MOAB/trunk/tools/dagmc/Makefile.am 2007-07-10 17:31:22 UTC (rev 1192)
+++ MOAB/trunk/tools/dagmc/Makefile.am 2007-07-10 22:59:36 UTC (rev 1193)
@@ -3,26 +3,28 @@
@CGM_CONFIG_OPTIONS@
DEFS = $(DEFINES)
-INCLUDES += -I$(top_srcdir) -I$(top_builddir) $(CGM_INCLUDES)
-CPPFLAGS += $(CGM_CPPFLAGS)
+INCLUDES += -I$(top_srcdir) -I$(top_builddir)
+
+lib_LTLIBRARIES = libdagmc.la
+libdagmc_la_SOURCES = DagMC.cpp
+libdagmc_la_includedir = $(includedir)
+libdagmc_la_include_HEADERS = cgm2moab.hpp cubfile.h DagMC.hpp
+
+if HAVE_CGM
+INCLUDES += $(CGM_INCLUDES)
+CPPFLAGS += $(CGM_CPPFLAGS) -DCGM
CXXFLAGS += $(CGM_CXXFLAGS)
+libdagmc_la_SOURCES += cgm2moab.cc
+
bin_PROGRAMS = cgm2moab
cgm2moab_SOURCES = main.cc cubfile.c
cgm2moab_LDFLAGS = $(CGM_LIBS_LTFLAGS)
-
-libdagmc_la_SOURCES = \
- cgm2moab.cc DagMC.cpp
-
-libMOAB_la_include_HEADERS = cgm2moab.hpp cubfile.h DagMC.hpp
-
-lib_LTLIBRARIES = libdagmc.la
-
-libMOAB_la_includedir = $(includedir)
-
-LDADD = libdagmc.la $(top_builddir)/libMOAB.la $(CGM_LIBS_LINK)
+cgm2moab_LDADD = libdagmc.la $(top_builddir)/libMOAB.la $(CGM_LIBS_LINK)
cgm2moab_DEPENDENCIES = libdagmc.la $(top_builddir)/libMOAB.la
dist_man1_MANS = cgm2moab.man
+endif
+
# Automake doesn't seem to have a directory defined for
# platform-dependent data (or include) files. So put
@@ -31,3 +33,9 @@
cfgdir = $(libdir)
cfg_DATA = dagmc.make
+TESTS = test_geom
+check_PROGRAMS = $(TESTS)
+test_geom_SOURCES = test_geom.cc
+test_geom_LDADD = libdagmc.la $(top_builddir)/libMOAB.la
+
+
Added: MOAB/trunk/tools/dagmc/test_geom.cc
===================================================================
--- MOAB/trunk/tools/dagmc/test_geom.cc (rev 0)
+++ MOAB/trunk/tools/dagmc/test_geom.cc 2007-07-10 22:59:36 UTC (rev 1193)
@@ -0,0 +1,489 @@
+#include "MBInterface.hpp"
+#include "MBCore.hpp"
+#include "DagMC.hpp"
+#include "MBTagConventions.hpp"
+
+#include <vector>
+#include <iostream>
+#include <math.h>
+#include <limits>
+
+#define CHKERR if (MB_SUCCESS != rval) return rval
+
+const double ROOT2 = 1.4142135623730951;
+
+// Create file containing geometry for 2x2x2 cube
+// centered at origin and having a convex +Z face
+// (center of face at origin).
+MBErrorCode write_geometry( const char* output_file_name );
+
+MBErrorCode test_ray_fire( DagMC& );
+
+MBErrorCode test_point_in_volume( DagMC& );
+
+MBErrorCode test_measure_volume( DagMC& );
+
+MBErrorCode test_measure_area( DagMC& );
+
+MBErrorCode test_surface_sense( DagMC& );
+
+
+MBErrorCode write_geometry( const char* output_file_name )
+{
+ MBErrorCode rval;
+ MBCore moab_instance;
+ MBInterface& moab = moab_instance;
+
+ // Define a 2x2x2 cube centered at orgin
+ // with concavity in +Z face.
+ const double coords[] = {
+ 1, -1, -1,
+ 1, 1, -1,
+ -1, 1, -1,
+ -1, -1, -1,
+ 1, -1, 1,
+ 1, 1, 1,
+ -1, 1, 1,
+ -1, -1, 1,
+ 0, 0, 0 };
+ const int connectivity[] = {
+ 0, 3, 1, 3, 2, 1, // -Z
+ 0, 1, 4, 5, 4, 1, // +X
+ 1, 2, 6, 6, 5, 1, // +Y
+ 6, 2, 3, 7, 6, 3, // -X
+ 0, 4, 3, 7, 3, 4, // -Y
+ 4, 5, 8, 5, 6, 8, // +Z
+ 6, 7, 8, 7, 4, 8 // +Z
+ };
+ const unsigned tris_per_surf[] = { 2, 2, 2, 2, 2, 4 };
+
+ // Create the geometry
+ const unsigned num_verts = sizeof(coords) / (3*sizeof(double));
+ const unsigned num_tris = sizeof(connectivity) / (3*sizeof(int));
+ const unsigned num_surfs = sizeof(tris_per_surf) / sizeof(unsigned);
+ MBEntityHandle verts[num_verts], tris[num_tris], surfs[num_surfs];
+ for (unsigned i = 0; i < num_verts; ++i) {
+ rval = moab.create_vertex( coords + 3*i, verts[i] );
+ CHKERR;
+ }
+ for (unsigned i = 0; i < num_tris; ++i) {
+ const MBEntityHandle conn[] = { verts[connectivity[3*i ]],
+ verts[connectivity[3*i+1]],
+ verts[connectivity[3*i+2]] };
+ rval = moab.create_element( MBTRI, conn, 3, tris[i] );
+ CHKERR;
+ }
+
+ // create CAD topology
+ MBEntityHandle* tri_iter = tris;
+ for (unsigned i = 0; i < num_surfs; ++i) {
+ rval = moab.create_meshset( MESHSET_SET, surfs[i] );
+ CHKERR;
+ rval = moab.add_entities( surfs[i], tri_iter, tris_per_surf[i] );
+ CHKERR;
+ tri_iter += tris_per_surf[i];
+ }
+
+ MBTag dim_tag, id_tag, sense_tag;
+ rval = moab.tag_create( GEOM_DIMENSION_TAG_NAME,
+ sizeof(int),
+ MB_TAG_SPARSE,
+ MB_TYPE_INTEGER,
+ dim_tag, 0, true );
+ CHKERR;
+ rval = moab.tag_create( GLOBAL_ID_TAG_NAME,
+ sizeof(int),
+ MB_TAG_DENSE,
+ MB_TYPE_INTEGER,
+ id_tag, 0, true );
+ CHKERR;
+ rval = moab.tag_create( "GEOM_SENSE_2",
+ 2*sizeof(MBEntityHandle),
+ MB_TAG_SPARSE,
+ MB_TYPE_HANDLE,
+ sense_tag, 0, true );
+ CHKERR;
+
+ std::vector<int> dims( num_surfs, 2 );
+ rval = moab.tag_set_data( dim_tag, surfs, num_surfs, &dims[0] );
+ CHKERR;
+ std::vector<int> ids( num_surfs );
+ for (size_t i = 0; i < ids.size(); ++i) ids[i] = i+1;
+ rval = moab.tag_set_data( id_tag, surfs, num_surfs, &ids[0] );
+ CHKERR;
+
+ MBEntityHandle volume;
+ rval = moab.create_meshset( MESHSET_SET, volume );
+ CHKERR;
+ for (unsigned i = 0; i < num_surfs; ++i) {
+ rval = moab.add_parent_child( volume, surfs[i] );
+ CHKERR;
+ }
+
+ std::vector<MBEntityHandle> senses( 2*num_surfs, 0 );
+ for (size_t i = 0; i < senses.size(); i += 2)
+ senses[i] = volume;
+ rval = moab.tag_set_data( sense_tag, surfs, num_surfs, &senses[0] );
+ CHKERR;
+
+ const int three = 3;
+ const int one = 1;
+ rval = moab.tag_set_data( dim_tag, &volume, 1, &three );
+ CHKERR;
+ rval = moab.tag_set_data( id_tag, &volume, 1, &one );
+ CHKERR;
+
+ rval = moab.write_mesh( output_file_name );
+ CHKERR;
+
+ return MB_SUCCESS;
+}
+
+static bool run_test( std::string name, int argc, char* argv[] )
+{
+ if (argc == 1)
+ return true;
+ for (int i = 1; i < argc; ++i)
+ if (name == argv[i])
+ return true;
+ return false;
+}
+
+#define RUN_TEST(A) do { \
+ if (run_test( #A, argc, argv )) { \
+ std::cout << #A << "... " << std::endl; \
+ if (MB_SUCCESS != A ( dagmc ) ) { \
+ ++errors; \
+ } \
+ } \
+} while(false)
+
+int main( int argc, char* argv[] )
+{
+ MBErrorCode rval;
+ const char* filename = "test_geom.h5m";
+
+ rval = write_geometry( filename );
+ if (MB_SUCCESS != rval) {
+ remove( filename );
+ std::cerr << "Failed to create input file: " << filename << std::endl;
+ return 1;
+ }
+
+ DagMC& dagmc = *DagMC::instance();
+ rval = dagmc.load_file_and_init( filename, strlen(filename), 0, 0 );
+ remove( filename );
+ if (MB_SUCCESS != rval) {
+ std::cerr << "Failed to initialize DagMC." << std::endl;
+ return 2;
+ }
+
+ int errors = 0;
+ RUN_TEST( test_ray_fire );
+ RUN_TEST( test_point_in_volume );
+ RUN_TEST( test_measure_volume );
+ RUN_TEST( test_measure_area );
+ RUN_TEST( test_surface_sense );
+
+ return errors;
+}
+
+MBErrorCode test_surface_sense( DagMC& dagmc )
+{
+ MBErrorCode rval;
+ MBInterface& moab = *dagmc.moab_instance();
+
+ MBTag dim_tag = dagmc.geom_tag();
+ MBRange surfs, vols;
+ const int two = 2, three = 3;
+ const void* ptrs[] = { &two, &three };
+ rval = moab.get_entities_by_type_and_tag( 0, MBENTITYSET, &dim_tag, ptrs, 1, surfs );
+ CHKERR;
+ rval = moab.get_entities_by_type_and_tag( 0, MBENTITYSET, &dim_tag, ptrs+1, 1, vols );
+ CHKERR;
+
+ if (vols.size() != 1) {
+ std::cerr << "ERROR: Expected 1 volume in input, found " << vols.size() << std::endl;
+ return MB_FAILURE;
+ }
+ if (surfs.size() != 6) {
+ std::cerr << "ERROR: Expected 6 surfaces in input, found " << surfs.size() << std::endl;
+ return MB_FAILURE;
+ }
+
+ for (MBRange::iterator i = surfs.begin(); i != surfs.end(); ++i) {
+ int sense = 0;
+ rval = dagmc.surface_sense( vols.front(), 1, &*i, &sense );
+ if (MB_SUCCESS != rval || sense != 1) {
+ std::cerr << "ERROR: Expected 1 for surface sense, got " << sense << std::endl;
+ return MB_FAILURE;
+ }
+ }
+
+ return MB_SUCCESS;
+}
+
+MBErrorCode test_measure_volume( DagMC& dagmc )
+{
+ MBErrorCode rval;
+ MBInterface& moab = *dagmc.moab_instance();
+
+ MBTag dim_tag = dagmc.geom_tag();
+ MBRange vols;
+ const int three = 3;
+ const void* ptr = &three;
+ rval = moab.get_entities_by_type_and_tag( 0, MBENTITYSET, &dim_tag, &ptr, 1, vols );
+ CHKERR;
+
+ if (vols.size() != 1) {
+ std::cerr << "ERROR: Expected 1 volume in input, found " << vols.size() << std::endl;
+ return MB_FAILURE;
+ }
+
+ // input file is 2x2x2 cube with convexity in +Z face that touches the origin.
+ // expected volume is 8 (2x2x2) less the volume of the pyrimid concavity
+ double result;
+ const double vol = 2*2*2 - 1*4./3;
+
+ rval = dagmc.measure_volume( vols.front(), result );
+ CHKERR;
+ if (fabs(result - vol) > std::numeric_limits<double>::epsilon()) {
+ std::cerr << "ERROR: Expected " << vol << " as measure of volume, got " << result << std::endl;
+ return MB_FAILURE;
+ }
+
+ return MB_SUCCESS;
+}
+
+MBErrorCode test_measure_area( DagMC& dagmc )
+{
+ MBErrorCode rval;
+ MBInterface& moab = *dagmc.moab_instance();
+
+ MBTag dim_tag = dagmc.geom_tag();
+ MBRange surfs;
+ const int two = 2;
+ const void* ptr = &two;
+ rval = moab.get_entities_by_type_and_tag( 0, MBENTITYSET, &dim_tag, &ptr, 1, surfs );
+ CHKERR;
+
+ if (surfs.size() != 6) {
+ std::cerr << "ERROR: Expected 6 surfaces in input, found " << surfs.size() << std::endl;
+ return MB_FAILURE;
+ }
+
+ int ids[6];
+ rval = moab.tag_get_data( dagmc.id_tag(), surfs, ids );
+ CHKERR;
+
+ // expect area of 4 for all faces except face 6.
+ // face 6 should have area == 4*sqrt(2)
+ MBRange::iterator iter = surfs.begin();
+ for (unsigned i = 0; i < 6; ++i, ++iter) {
+ double expected = 4.0;
+ if (ids[i] == 6)
+ expected *= ROOT2;
+
+ double result;
+
+ rval = dagmc.measure_area( *iter, result );
+ CHKERR;
+ if (fabs(result - expected) > std::numeric_limits<double>::epsilon()) {
+ std::cerr << "ERROR: Expected area of surface " << ids[i] << " to be "
+ << expected << ". Got " << result << std::endl;
+ return MB_FAILURE;
+ }
+ }
+
+ return MB_SUCCESS;
+}
+
+
+struct ray_fire {
+ int prev_surf;
+ double origin[3], direction[3];
+ int hit_surf;
+ double distance;
+};
+
+MBErrorCode test_ray_fire( DagMC& dagmc )
+{
+ const struct ray_fire tests[] = {
+ /* src origin direction dest dist */
+ { 1, { 0.0, 0.0, -1. }, { -1.0/ROOT2, 0.0, 1.0/ROOT2 }, 4, ROOT2 },
+ { 1, { 0.0, 0.0, -1. }, { 1.0/ROOT2, 0.0, 1.0/ROOT2 }, 2, ROOT2 },
+ { 1, { 0.0, 0.0, -1. }, { 0.0, 1.0/ROOT2, 1.0/ROOT2 }, 3, ROOT2 },
+ { 1, { 0.5, 0.5, -1. }, { 0.0, 0.0, 1.0 }, 6, 3 },
+ { 2, { 1.0, 0.0, 0.5 }, { -1.0, 0.0, 0.0 }, 6, 1 } };
+
+ MBErrorCode rval;
+ MBInterface& moab = *dagmc.moab_instance();
+
+ MBTag dim_tag = dagmc.geom_tag();
+ MBRange surfs, vols;
+ const int two = 2, three = 3;
+ const void* ptrs[] = { &two, &three };
+ rval = moab.get_entities_by_type_and_tag( 0, MBENTITYSET, &dim_tag, ptrs, 1, surfs );
+ CHKERR;
+ rval = moab.get_entities_by_type_and_tag( 0, MBENTITYSET, &dim_tag, ptrs+1, 1, vols );
+ CHKERR;
+
+ if (vols.size() != 1) {
+ std::cerr << "ERROR: Expected 1 volume in input, found " << vols.size() << std::endl;
+ return MB_FAILURE;
+ }
+ if (surfs.size() != 6) {
+ std::cerr << "ERROR: Expected 6 surfaces in input, found " << surfs.size() << std::endl;
+ return MB_FAILURE;
+ }
+
+ int ids[6];
+ rval = moab.tag_get_data( dagmc.id_tag(), surfs, ids );
+ CHKERR;
+ MBEntityHandle surf[6];
+ std::copy( surfs.begin(), surfs.end(), surf );
+
+ const int num_test = sizeof(tests) / sizeof(tests[0]);
+ for (int i = 0; i < num_test; ++i) {
+ int* ptr = std::find( ids, ids+6, tests[i].prev_surf );
+ int idx = ptr - ids;
+ if (idx >= 6) {
+ std::cerr << "Surface " << tests[i].prev_surf << " not found." << std::endl;
+ return MB_FAILURE;
+ }
+ const MBEntityHandle src_surf = surf[idx];
+
+ ptr = std::find( ids, ids+6, tests[i].hit_surf );
+ idx = ptr - ids;
+ if (idx >= 6) {
+ std::cerr << "Surface " << tests[i].hit_surf << " not found." << std::endl;
+ return MB_FAILURE;
+ }
+ const MBEntityHandle hit_surf = surf[idx];
+
+ double dist;
+ MBEntityHandle result;
+ rval = dagmc.ray_fire( vols.front(),
+ src_surf,
+ 2,
+ tests[i].direction[0],
+ tests[i].direction[1],
+ tests[i].direction[2],
+ tests[i].origin[0],
+ tests[i].origin[1],
+ tests[i].origin[2],
+ HUGE_VAL,
+ dist, result );
+
+ if (result != hit_surf || fabs(dist - tests[i].distance) > 1e-6) {
+ MBEntityHandle *p = std::find( surf, surf+6, result );
+ idx = p - surf;
+ int id = idx > 5 ? 0 : ids[idx];
+
+ std::cerr << "Rayfire test failed for " << std::endl
+ << "\t ray from (" << tests[i].origin[0]
+ << ", " << tests[i].origin[1] << ", "
+ << tests[i].origin[2] << ") going ["
+ << tests[i].direction[0] << ", "
+ << tests[i].direction[1] << ", "
+ << tests[i].direction[2] << "]" << std::endl
+ << "\t Beginning on surface " << tests[i].prev_surf << std::endl
+ << "\t Expected to hit surface " << tests[i].hit_surf << " after "
+ << tests[i].distance << " units." << std::endl
+ << "\t Actually hit surface " << id << " after " << dist << " units."
+ << std::endl;
+ return MB_FAILURE;
+ }
+ }
+
+ return MB_SUCCESS;
+}
+
+struct PointInVol { double coords[3]; int result; };
+
+MBErrorCode test_point_in_volume( DagMC& dagmc )
+{
+ const char* const NAME_ARR[] = { "Boundary", "Outside", "Inside" };
+ const char* const* names = NAME_ARR + 1;
+ const int INSIDE = 1, OUTSIDE = 0, BOUNDARY = -1;
+ const struct PointInVol tests[] = {
+ { { 0.0, 0.0, 0.5 }, OUTSIDE },
+ { { 0.0, 0.0,-0.5 }, INSIDE },
+ { { 0.7, 0.0, 0.0 }, INSIDE },
+ { {-0.7, 0.0, 0.0 }, INSIDE },
+ { { 0.0,-0.7, 0.0 }, INSIDE },
+ { { 0.0,-0.7, 0.0 }, INSIDE },
+ { { 1.1, 1.1, 1.1 }, OUTSIDE },
+ { {-1.1, 1.1, 1.1 }, OUTSIDE },
+ { {-1.1,-1.1, 1.1 }, OUTSIDE },
+ { { 1.1,-1.1, 1.1 }, OUTSIDE },
+ { { 1.1, 1.1,-1.1 }, OUTSIDE },
+ { {-1.1, 1.1,-1.1 }, OUTSIDE },
+ { {-1.1,-1.1,-1.1 }, OUTSIDE },
+ { { 1.1,-1.1,-1.1 }, OUTSIDE },
+ { { 1.0, 0.0, 0.0 }, BOUNDARY},
+ { {-1.0, 0.0, 0.0 }, BOUNDARY},
+ { { 0.0, 1.0, 0.0 }, BOUNDARY},
+ { { 0.0,-1.0, 0.0 }, BOUNDARY},
+ { { 0.0, 0.0, 0.0 }, BOUNDARY},
+ { { 0.0, 0.0,-1.0 }, BOUNDARY} };
+ const int num_test = sizeof(tests) / sizeof(tests[0]);
+
+ MBErrorCode rval;
+ MBInterface& moab = *dagmc.moab_instance();
+
+ MBTag dim_tag = dagmc.geom_tag();
+ MBRange vols;
+ const int three = 3;
+ const void* ptr = &three;
+ rval = moab.get_entities_by_type_and_tag( 0, MBENTITYSET, &dim_tag, &ptr, 1, vols );
+ CHKERR;
+ if (vols.size() != 1) {
+ std::cerr << "ERROR: Expected 1 volume in input, found " << vols.size() << std::endl;
+ return MB_FAILURE;
+ }
+ const MBEntityHandle vol = vols.front();
+
+ for (int i = 0; i < num_test; ++i) {
+ int result;
+
+ rval = dagmc.point_in_volume( vol,
+ tests[i].coords[0],
+ tests[i].coords[1],
+ tests[i].coords[2],
+ result );
+ CHKERR;
+ if (result != tests[i].result) {
+ std::cerr << "ERROR testing point_in_volume[" << i << "]:" << std::endl
+ << "\tExpected " << names[tests[i].result]
+ << " for (" << tests[i].coords[0] << ", "
+ << tests[i].coords[1] << ", " << tests[i].coords[2]
+ << "). Got " << names[result] << std::endl;
+ return MB_FAILURE;
+ }
+
+ // point_in_volume_slow doesn't to boundary.
+ if (tests[i].result == BOUNDARY)
+ continue;
+
+ rval = dagmc.point_in_volume_slow( vol,
+ tests[i].coords[0],
+ tests[i].coords[1],
+ tests[i].coords[2],
+ result );
+ CHKERR;
+
+ if (result != tests[i].result) {
+ std::cerr << "ERROR testing point_in_volume_slow[" << i << "]:" << std::endl
+ << "\tExpected " << names[tests[i].result]
+ << " for (" << tests[i].coords[0] << ", "
+ << tests[i].coords[1] << ", " << tests[i].coords[2]
+ << "). Got " << names[result] << std::endl;
+ return MB_FAILURE;
+ }
+ }
+
+ return MB_SUCCESS;
+}
+
+
More information about the moab-dev
mailing list