July 2007 - moab-dev - mailman3.cels.anl.gov

r1224 - in MOAB/trunk: . test tools/qvdual
by tautges＠mcs.anl.gov 31 Jul '07

31 Jul '07

Author: tautges Date: 2007-07-31 07:32:25 -0500 (Tue, 31 Jul 2007) New Revision: 1224 Modified: MOAB/trunk/MeshTopoUtil.cpp MOAB/trunk/Tqdcfr.cpp MOAB/trunk/configure.in MOAB/trunk/test/Makefile.am MOAB/trunk/tools/qvdual/DrawDual.cpp MOAB/trunk/tools/qvdual/DrawDual.hpp MOAB/trunk/tools/qvdual/vtkMOABUtils.cxx MOAB/trunk/tools/qvdual/vtkMOABUtils.h Log: tools/qvdual/DrawDual.cpp/.hpp, vtkMOABUtils.cxx/.h, uiQVDual.ui.h: corrections to make qvdual work with more recent versions of vtk and graphviz configure.in: changing to make work with new versions of vtk and graphviz MeshTopoUtil.cpp: fixing typo in get_bridge_entities Tqdcfr.cpp: making various tag types integer rather than opaque test/Makefile.am: changing .cub files included with make dist Modified: MOAB/trunk/MeshTopoUtil.cpp =================================================================== --- MOAB/trunk/MeshTopoUtil.cpp 2007-07-26 06:43:08 UTC (rev 1223) +++ MOAB/trunk/MeshTopoUtil.cpp 2007-07-31 12:32:25 UTC (rev 1224) @@ -392,7 +392,7 @@ // get the vertices making up this sub-entity int num_bridge_verts = MBCN::VerticesPerEntity( MBCN::SubEntityType( from_type, bridge_dim, i ) ); MBCN::SubEntityVertexIndices( from_type, bridge_dim, i, bridge_indices ); - for (int j = 0; j < num_bridge_verts; ++i) + for (int j = 0; j < num_bridge_verts; ++j) bridge_verts[j]= connect[bridge_indices[j]]; //MBCN::SubEntityConn(connect, from_type, bridge_dim, i, &bridge_verts[0], num_bridge_verts); Modified: MOAB/trunk/Tqdcfr.cpp =================================================================== --- MOAB/trunk/Tqdcfr.cpp 2007-07-26 06:43:08 UTC (rev 1223) +++ MOAB/trunk/Tqdcfr.cpp 2007-07-31 12:32:25 UTC (rev 1224) @@ -150,7 +150,7 @@ MBErrorCode result = mdbImpl->tag_get_handle("cubIdTag", cubIdTag); if (MB_SUCCESS != result || MB_TAG_NOT_FOUND == result) { int default_val = -1; - result = mdbImpl->tag_create("cubIdTag", 4, MB_TAG_DENSE, + result = mdbImpl->tag_create("cubIdTag", sizeof(int), MB_TAG_DENSE, MB_TYPE_INTEGER, cubIdTag, &default_val); assert (MB_SUCCESS == result); } @@ -368,8 +368,8 @@ if (0 != nodeset_offset) { if (0 == nsTag) { int default_val = 0; - tmp_result = mdbImpl->tag_create(DIRICHLET_SET_TAG_NAME, 4, MB_TAG_SPARSE, - nsTag, &default_val); + tmp_result = mdbImpl->tag_create(DIRICHLET_SET_TAG_NAME, sizeof(int), MB_TAG_SPARSE, + MB_TYPE_INTEGER, nsTag, &default_val); if (MB_SUCCESS != tmp_result) result = tmp_result; } if (MB_SUCCESS == tmp_result) @@ -382,8 +382,8 @@ if (0 != sideset_offset) { if (0 == ssTag) { int default_val = 0; - tmp_result = mdbImpl->tag_create(NEUMANN_SET_TAG_NAME, 4, MB_TAG_SPARSE, - ssTag, &default_val); + tmp_result = mdbImpl->tag_create(NEUMANN_SET_TAG_NAME, sizeof(int), MB_TAG_SPARSE, + MB_TYPE_INTEGER, ssTag, &default_val); if (MB_SUCCESS != tmp_result) result = tmp_result; } if (MB_SUCCESS == tmp_result) @@ -513,7 +513,8 @@ FREADD(sideseth->numDF); MBTag distFactorTag; double *dum_val = NULL; - MBErrorCode result = mdbImpl->tag_create("distFactor", sizeof(double*), MB_TAG_SPARSE, distFactorTag, + MBErrorCode result = mdbImpl->tag_create("distFactor", sizeof(double*), MB_TAG_SPARSE, + distFactorTag, &dum_val); if (MB_SUCCESS != result && MB_ALREADY_ALLOCATED != result) return result; double *dist_data = new double[ss_dfs.size()]; @@ -573,7 +574,7 @@ int def_val = 1; MBTag sense_tag; tmp_result = mdbImpl->tag_create("SENSE", sizeof(int), - MB_TAG_SPARSE, sense_tag, &def_val); + MB_TAG_SPARSE, MB_TYPE_INTEGER, sense_tag, &def_val); if (tmp_result != MB_SUCCESS && tmp_result != MB_ALREADY_ALLOCATED) result = tmp_result; def_val = -1; tmp_result = mdbImpl->tag_set_data(sense_tag, &reverse_set, 1, &def_val); @@ -632,7 +633,7 @@ int def_val = 1; MBTag sense_tag; tmp_result = mdbImpl->tag_create("SENSE", sizeof(int), - MB_TAG_SPARSE, sense_tag, &def_val); + MB_TAG_SPARSE, MB_TYPE_INTEGER, sense_tag, &def_val); if (tmp_result != MB_SUCCESS && tmp_result != MB_ALREADY_ALLOCATED) result = tmp_result; def_val = -1; tmp_result = mdbImpl->tag_set_data(sense_tag, &reverse_set, 1, &def_val); @@ -682,7 +683,7 @@ // now do something with them... MBErrorCode result = mdbImpl->tag_create("Block_Attributes", blockh->attribOrder*sizeof(double), MB_TAG_SPARSE, - block_attribs, NULL); + MB_TYPE_DOUBLE, block_attribs, NULL); if (MB_SUCCESS != result && MB_ALREADY_ALLOCATED != result) return result; result = mdbImpl->tag_set_data(block_attribs, &(blockh->setHandle), 1, &(dbl_buf[0])); @@ -1030,7 +1031,7 @@ if (MB_SUCCESS != result || 0 == fixedFlagTag) { int dum_val = 0; result = mdbImpl->tag_create("NodeFixed", sizeof(int), MB_TAG_SPARSE, - fixedFlagTag, &dum_val); + MB_TYPE_INTEGER, fixedFlagTag, &dum_val); if (MB_SUCCESS != result) return result; } } @@ -1496,8 +1497,8 @@ int has_mid_nodes[] = {0, 0, 0, 0}; static MBTag hasMidNodesTag = 0; if (0 == hasMidNodesTag) { - result = instance->mdbImpl->tag_create(HAS_MID_NODES_TAG_NAME, 4*sizeof(int), MB_TAG_SPARSE, hasMidNodesTag, - has_mid_nodes); + result = instance->mdbImpl->tag_create(HAS_MID_NODES_TAG_NAME, 4*sizeof(int), MB_TAG_SPARSE, + MB_TYPE_INTEGER, hasMidNodesTag, has_mid_nodes); if (MB_SUCCESS != result && MB_ALREADY_ALLOCATED != result) return result; } @@ -1669,16 +1670,16 @@ int default_val = -1; MBErrorCode result; - result = instance->mdbImpl->tag_create(GLOBAL_ID_TAG_NAME, 4, MB_TAG_DENSE, + result = instance->mdbImpl->tag_create(GLOBAL_ID_TAG_NAME, sizeof(int), MB_TAG_DENSE, MB_TYPE_INTEGER, instance->globalIdTag, &default_val); if (MB_SUCCESS != result && MB_ALREADY_ALLOCATED != result) return result; if (feModelHeader.geomArray.numEntities > 0) { - result = instance->mdbImpl->tag_create(GEOM_DIMENSION_TAG_NAME, 4, MB_TAG_SPARSE, - instance->geomTag, &default_val); + result = instance->mdbImpl->tag_create(GEOM_DIMENSION_TAG_NAME, sizeof(int), MB_TAG_SPARSE, + MB_TYPE_INTEGER, instance->geomTag, &default_val); if (MB_SUCCESS != result && MB_ALREADY_ALLOCATED != result) return result; - result = instance->mdbImpl->tag_create("UNIQUE_ID", 4, MB_TAG_SPARSE, + result = instance->mdbImpl->tag_create("UNIQUE_ID", sizeof(int), MB_TAG_SPARSE, MB_TYPE_INTEGER, instance->uniqueIdTag, &default_val); if (MB_SUCCESS != result && MB_ALREADY_ALLOCATED != result) return result; @@ -1700,8 +1701,8 @@ } if (feModelHeader.blockArray.numEntities > 0) { - result = instance->mdbImpl->tag_create(MATERIAL_SET_TAG_NAME, 4, MB_TAG_SPARSE, - instance->blockTag, &default_val); + result = instance->mdbImpl->tag_create(MATERIAL_SET_TAG_NAME, sizeof(int), MB_TAG_SPARSE, + MB_TYPE_INTEGER, instance->blockTag, &default_val); if (MB_SUCCESS != result && MB_ALREADY_ALLOCATED != result) return result; result = Tqdcfr::BlockHeader::read_info_header(data_version, modelOffset, @@ -1712,8 +1713,8 @@ print_block_headers("Block headers:", feBlockH, feModelHeader.blockArray.numEntities); } if (feModelHeader.nodesetArray.numEntities > 0) { - result = instance->mdbImpl->tag_create(DIRICHLET_SET_TAG_NAME, 4, MB_TAG_SPARSE, - instance->nsTag, &default_val); + result = instance->mdbImpl->tag_create(DIRICHLET_SET_TAG_NAME, sizeof(int), MB_TAG_SPARSE, + MB_TYPE_INTEGER, instance->nsTag, &default_val); if (MB_SUCCESS != result && MB_ALREADY_ALLOCATED != result) return result; result = Tqdcfr::NodesetHeader::read_info_header(modelOffset, @@ -1724,8 +1725,8 @@ print_nodeset_headers("Nodeset headers:", feNodeSetH, feModelHeader.nodesetArray.numEntities); } if (feModelHeader.sidesetArray.numEntities > 0) { - result = instance->mdbImpl->tag_create(NEUMANN_SET_TAG_NAME, 4, MB_TAG_SPARSE, - instance->ssTag, &default_val); + result = instance->mdbImpl->tag_create(NEUMANN_SET_TAG_NAME, sizeof(int), MB_TAG_SPARSE, + MB_TYPE_INTEGER, instance->ssTag, &default_val); if (MB_SUCCESS != result && MB_ALREADY_ALLOCATED != result) return result; result = Tqdcfr::SidesetHeader::read_info_header(modelOffset, Modified: MOAB/trunk/configure.in =================================================================== --- MOAB/trunk/configure.in 2007-07-26 06:43:08 UTC (rev 1223) +++ MOAB/trunk/configure.in 2007-07-31 12:32:25 UTC (rev 1224) @@ -354,7 +354,7 @@ GRAPH_DIR= for dir in /usr /usr/local; do if test "x" = "x$GRAPH_DIR"; then - AC_CHECK_FILE([${dir}/include/graphviz/dotneato.h],[GRAPH_DIR="$dir"]) + AC_CHECK_FILE([${dir}/include/graphviz/gvc.h],[GRAPH_DIR="$dir"]) fi done AC_MSG_CHECKING([for GraphViz include dir]) @@ -377,7 +377,7 @@ GRAPH_INC_DIR= for dir in ${GRAPH_DIR}/include/graphviz ${GRAPH_DIR}/include ${GRAPH_DIR}; do if test "x" = "x$GRAPH_INC_DIR"; then - AC_CHECK_FILE([${dir}/dotneato.h],[GRAPH_INC_DIR="-I$dir"]) + AC_CHECK_FILE([${dir}/gvc.h],[GRAPH_INC_DIR="-I$dir"]) fi done AC_MSG_CHECKING([for GraphViz include dir]) @@ -405,16 +405,15 @@ ;; esac - GRAPH_LIBS="$GRAPH_LIB_DIR -ldotneato -lgraph" + GRAPH_LIBS="$GRAPH_LIB_DIR -lgvc -lgraph" old_CXXFLAGS="$CXXFLAGS" old_CPPFLAGS="$CPPFLAGS" CXXFLAGS="$GRAPH_INC_DIR" CPPFLAGS="$GRAPH_INC_DIR" old_LIBS="$LIBS" LIBS=$GRAPH_LIBS - AC_CHECK_HEADERS([dotneato.h],[],[AC_MSG_WARN([Missing GraphViz header]); GRAPHVIZ_MISSING=yes]) - AC_CHECK_LIB([graph],[aginitlib],[],[AC_MSG_WARN([Cannot find GraphViz library: -lgraph]); GRAPHVIZ_MISSING=yes]) - AC_CHECK_LIB([dotneato],[gvContext],[],[AC_MSG_WARN([Cannot find GraphViz library: -ldotneato]); GRAPHVIZ_MISSING=yes],[-lgraph]) + AC_CHECK_HEADERS([gvc.h],[],[AC_MSG_WARN([Missing GraphViz header]); GRAPHVIZ_MISSING=yes]) + AC_CHECK_LIB([gvc],[gvContext],[],[AC_MSG_WARN([Cannot find GraphViz library: -lgraph]); GRAPHVIZ_MISSING=yes]) LIBS="$old_LIBS" CXXFLAGS="$old_CXXFLAGS" CPPFLAGS="$old_CPPFLAGS" @@ -450,7 +449,7 @@ xyes|x) VTK_LIB_DIR= VTK_INC_DIR= - for dir in /usr/include/vtk /usr/local/include/vtk /include/vtk; do + for dir in /usr/include/vtk /usr/local/include/vtk /include/vtk /usr/include/vtk-5.0; do if test "x" = "x$VTK_INC_DIR"; then AC_CHECK_FILE([${dir}/vtkRenderer.h],[VTK_INC_DIR="-I$dir"]) fi @@ -461,11 +460,24 @@ else VTK_INCLUDES="$VTK_INC_DIR" fi + old_LIBS="$LIBS" + for dir in /usr/lib/vtk /usr/lib /usr/bin /usr; do + if test "x" = "x$VTK_LIB_DIR"; then + LIBS="$old_LIBS -L$dir" + unset ac_cv_lib_vtkCommon + unset ac_cv_lib_vtkCommon_main + AC_HAVE_LIBRARY([vtkCommon],[VTK_LIB_DIR=$dir]) + fi + done + if test "x" = "x$VTK_LIB_DIR"; then + VTK_MISSING=yes + fi + LIBS="$old_LIBS" ;; x?*) VTK_LIB_DIR= VTK_INC_DIR= - for dir in ${VTK_DIR}/include/vtk ${VTK_DIR}/include ${VTK_DIR}; do + for dir in ${VTK_DIR}/include/vtk ${VTK_DIR}/include/vtk/Rendering ${VTK_DIR}/include ${VTK_DIR}; do if test "x" = "x$VTK_INC_DIR"; then AC_CHECK_FILE([${dir}/vtkRenderer.h],[VTK_INC_DIR="-I$dir"]) fi @@ -490,9 +502,9 @@ for dir in ${VTK_DIR}/lib/vtk ${VTK_DIR}/lib ${VTK_DIR}/bin ${VTK_DIR}; do if test "x" = "x$VTK_LIB_DIR"; then LIBS="$old_LIBS -L$dir" - unset ac_cv_lib_vtkexpat - unset ac_cv_lib_vtkexpat_main - AC_HAVE_LIBRARY([vtkexpat],[VTK_LIB_DIR=$dir]) + unset ac_cv_lib_vtkCommon + unset ac_cv_lib_vtkCommon_main + AC_HAVE_LIBRARY([vtkCommon],[VTK_LIB_DIR=$dir]) fi done if test "x" = "x$VTK_LIB_DIR"; then @@ -549,7 +561,7 @@ SNL_CHECK_LIB_LIST([vtkexoIIc], [vtk_libs],[$vtk_libs_arg],[-L$VTK_LIB_DIR]) SNL_CHECK_LIB_LIST([vtkIO], [vtk_libs],[$vtk_libs_arg],[-L$VTK_LIB_DIR]) SNL_CHECK_LIB_LIST([vtkRendering], [vtk_libs],[$vtk_libs_arg],[-L$VTK_LIB_DIR]) - SNL_CHECK_LIB_LIST([vtkParallel], [vtk_libs],[$vtk_libs_arg],[-L$VTK_LIB_DIR]) + SNL_CHECK_LIB_LIST([vtkParallel], [vtk_libs],[$vtk_libs_arg],[-L$VTK_LIB_DIR]) SNL_CHECK_LIB_LIST([vtkHybrid], [vtk_libs],[$vtk_libs_arg],[-L$VTK_LIB_DIR]) # Make sure we found at least vtkCommon and vtkRendering Modified: MOAB/trunk/test/Makefile.am =================================================================== --- MOAB/trunk/test/Makefile.am 2007-07-26 06:43:08 UTC (rev 1223) +++ MOAB/trunk/test/Makefile.am 2007-07-31 12:32:25 UTC (rev 1224) @@ -10,7 +10,9 @@ 4k-tri-plane.vtk \ 4tet.g \ 8hex.g \ - block.cub \ + brick_cubit10.1.cub \ + brick_cubit10.2.cub \ + brick_cubit10.cub \ cell1.gen.gz \ cell2.gen.gz \ mb_big_test.g.gz \ Modified: MOAB/trunk/tools/qvdual/DrawDual.cpp =================================================================== --- MOAB/trunk/tools/qvdual/DrawDual.cpp 2007-07-26 06:43:08 UTC (rev 1223) +++ MOAB/trunk/tools/qvdual/DrawDual.cpp 2007-07-31 12:32:25 UTC (rev 1224) @@ -44,11 +44,11 @@ extern "C" { -#include "dotneato.h" +#include "gvc.h" // extern GVC_t *gvContext(); } -const bool my_debug = false; +const bool my_debug = true; const int SHEET_WINDOW_SIZE = 500; @@ -227,7 +227,7 @@ else { // more than one - traverse, choosing the chord if any first actors->InitTraversal(); - while (tmp_actor = actors->GetNextItem()) { + while ((tmp_actor = actors->GetNextItem())) { MBEntityHandle this_set = vtkMOABUtils::propSetMap[tmp_actor]; if (0 == this_set) continue; // get whether it's a dual surface or dual curve @@ -443,7 +443,7 @@ agwrite(this_gw.gvizGraph, stdout); } // neato_init_graph(this_gw.gvizGraph); - neato_layout(this_gw.gvizGraph); + gvLayout(gvContext(), this_gw.gvizGraph, "neato"); // adjustNodes(this_gw.gvizGraph); // spline_edges(this_gw.gvizGraph); // dotneato_postprocess(this_gw.gvizGraph, neato_nodesize); @@ -1320,8 +1320,8 @@ endloop = loop_verts.end(); if ((firstv == endloop && lastv != endloop) || (firstv != endloop && lastv == endloop)) { - agxset(this_gved, asym_weight->index, "10"); - agxset(this_gved, asym_len->index, "0.1"); + // agxset(this_gved, asym_weight->index, "10"); + // agxset(this_gved, asym_len->index, "0.1"); } } @@ -1789,8 +1789,8 @@ { // get the ids of these verts, equal to entity ids of their primal entities static std::vector<MBEntityHandle> primals; - primals.reserve(ents.size()); - ids.reserve(ents.size()); + primals.resize(ents.size()); + ids.resize(ents.size()); MBErrorCode result = MBI->tag_get_data(dualEntityTagHandle, ents, &primals[0]); if (MB_SUCCESS != result) { for (unsigned int i = 0; i < ents.size(); i++) ids[i] = 0; @@ -1798,6 +1798,9 @@ result = MBI->tag_get_data(vtkMOABUtils::globalId_tag(), &primals[0], ents.size(), &ids[0]); + if (MB_SUCCESS != result && MB_TAG_NOT_FOUND != result) + std::cerr << "tag_get_data returned non-zero result in get_primal_ids." << std::endl; + for (unsigned int i = 0; i < ents.size(); i++) { if (0 == ids[i]) ids[i] = MBI->id_from_handle(primals[i]); } @@ -1838,7 +1841,7 @@ MBRange saved_sets; // get all actors, check sets in propSetMap; save set, remove actor from map - while (tmp_actor = acoll->GetNextItem()) { + while ((tmp_actor = acoll->GetNextItem())) { std::map<vtkProp*,MBEntityHandle>::iterator mit = vtkMOABUtils::propSetMap.find(tmp_actor); if (mit == vtkMOABUtils::propSetMap.end()) continue; Modified: MOAB/trunk/tools/qvdual/DrawDual.hpp =================================================================== --- MOAB/trunk/tools/qvdual/DrawDual.hpp 2007-07-26 06:43:08 UTC (rev 1223) +++ MOAB/trunk/tools/qvdual/DrawDual.hpp 2007-07-31 12:32:25 UTC (rev 1224) @@ -4,6 +4,7 @@ class Agnode_t; class Agedge_t; class Agraph_t; + class QVTKWidget; #include "MBInterface.hpp" Modified: MOAB/trunk/tools/qvdual/vtkMOABUtils.cxx =================================================================== --- MOAB/trunk/tools/qvdual/vtkMOABUtils.cxx 2007-07-26 06:43:08 UTC (rev 1223) +++ MOAB/trunk/tools/qvdual/vtkMOABUtils.cxx 2007-07-31 12:32:25 UTC (rev 1224) @@ -123,7 +123,7 @@ MBTag vtkMOABUtils::globalIdTag = 0; MBTag vtkMOABUtils::categoryTag = 0; -bool debug = false; +bool vtkMOABUtils::debug = false; //vtkStandardNewMacro(vtkMOABUtils); Modified: MOAB/trunk/tools/qvdual/vtkMOABUtils.h =================================================================== --- MOAB/trunk/tools/qvdual/vtkMOABUtils.h 2007-07-26 06:43:08 UTC (rev 1223) +++ MOAB/trunk/tools/qvdual/vtkMOABUtils.h 2007-07-31 12:32:25 UTC (rev 1224) @@ -242,6 +242,9 @@ //! assign global ids, in preparation for writing mesh static void assign_global_ids(); + + //! debug parameter + static bool debug; private: //! private constructor so nobody can construct one

1 0

r1223 - MOAB/trunk
by dcthomp＠mcs.anl.gov 26 Jul '07

26 Jul '07

Author: dcthomp Date: 2007-07-26 01:43:08 -0500 (Thu, 26 Jul 2007) New Revision: 1223 Modified: MOAB/trunk/configure.in Log: ENH: Recent VTK builds have some extra libraries that must be linked. Modified: MOAB/trunk/configure.in =================================================================== --- MOAB/trunk/configure.in 2007-07-23 19:20:35 UTC (rev 1222) +++ MOAB/trunk/configure.in 2007-07-26 06:43:08 UTC (rev 1223) @@ -531,6 +531,7 @@ SNL_CHECK_LIB_LIST([GL], [vtk_x_libs],[$x_libs_arg],[-L$x_libraries]) SNL_CHECK_LIB_LIST([GLU], [vtk_x_libs],[$x_libs_arg],[-L$x_libraries]) vtk_libs_arg="-L$x_libraries $vtk_x_libs $X_EXTRA_LIBS -lX11 $X_PRE_LIBS $vtk_sys_libs" + SNL_CHECK_LIB_LIST([vtksys], [vtk_libs],[$vtk_libs_arg],[-L$VTK_LIB_DIR]) SNL_CHECK_LIB_LIST([vtkCommon], [vtk_libs],[$vtk_libs_arg],[-L$VTK_LIB_DIR]) SNL_CHECK_LIB_LIST([vtkDICOMParser],[vtk_libs],[$vtk_libs_arg],[-L$VTK_LIB_DIR]) SNL_CHECK_LIB_LIST([vtkfreetype], [vtk_libs],[$vtk_libs_arg],[-L$VTK_LIB_DIR]) @@ -542,10 +543,14 @@ SNL_CHECK_LIB_LIST([vtkftgl], [vtk_libs],[$vtk_libs_arg],[-L$VTK_LIB_DIR]) SNL_CHECK_LIB_LIST([vtkFiltering], [vtk_libs],[$vtk_libs_arg],[-L$VTK_LIB_DIR]) SNL_CHECK_LIB_LIST([vtkImaging], [vtk_libs],[$vtk_libs_arg],[-L$VTK_LIB_DIR]) + SNL_CHECK_LIB_LIST([verdict112], [vtk_libs],[$vtk_libs_arg],[-L$VTK_LIB_DIR]) SNL_CHECK_LIB_LIST([vtkGraphics], [vtk_libs],[$vtk_libs_arg],[-L$VTK_LIB_DIR]) + SNL_CHECK_LIB_LIST([vtkNetCDF], [vtk_libs],[$vtk_libs_arg],[-L$VTK_LIB_DIR]) + SNL_CHECK_LIB_LIST([vtkexoIIc], [vtk_libs],[$vtk_libs_arg],[-L$VTK_LIB_DIR]) SNL_CHECK_LIB_LIST([vtkIO], [vtk_libs],[$vtk_libs_arg],[-L$VTK_LIB_DIR]) SNL_CHECK_LIB_LIST([vtkRendering], [vtk_libs],[$vtk_libs_arg],[-L$VTK_LIB_DIR]) SNL_CHECK_LIB_LIST([vtkParallel], [vtk_libs],[$vtk_libs_arg],[-L$VTK_LIB_DIR]) + SNL_CHECK_LIB_LIST([vtkHybrid], [vtk_libs],[$vtk_libs_arg],[-L$VTK_LIB_DIR]) # Make sure we found at least vtkCommon and vtkRendering if ! echo "x $vtk_libs" | grep lvtkCommon >/dev/null; then

1 0

r1222 - in MOAB/trunk: . test/obb tools/dagmc
by kraftche＠mcs.anl.gov 23 Jul '07

23 Jul '07

Author: kraftche Date: 2007-07-23 14:20:35 -0500 (Mon, 23 Jul 2007) New Revision: 1222 Modified: MOAB/trunk/MBOrientedBoxTreeTool.cpp MOAB/trunk/MBOrientedBoxTreeTool.hpp MOAB/trunk/test/obb/obb_test.cpp MOAB/trunk/tools/dagmc/DagMC.cpp Log: Change algorithm for closest-point and point-in-volume. o Old closest point function did: - find closest point on boundary - calculate distance d between input point and closest point - return all points within d+epsilon of input point o Old closest point function very susceptable to rounding error if input point is far away (d+epsilon == d for large values of d). o New closest point function just returns the closest single point o New addtional function to intersect sphere with triangles o New point-in-volume test does: - get closest point - get any triangles within epsilon of closest point by doing sphere-intersect-trinagles where center == closest point and radius == epsilon. This method is better because o eliminates potential source of rounding error o greatly simplifies point-in-volume because the result set of triangles are all topologically adjacent - no need to search larger set of tris from old closest-point method for a subset that are adjacent This method might be slower because of the second tree traversal. However, it may not be because the closest-point query is a little faster and the sphere-intersect-triangles query is much quicker than either closest-point query. Fixes bug from Paul (tube_1e-6.h5m). Modified: MOAB/trunk/MBOrientedBoxTreeTool.cpp =================================================================== --- MOAB/trunk/MBOrientedBoxTreeTool.cpp 2007-07-20 18:24:11 UTC (rev 1221) +++ MOAB/trunk/MBOrientedBoxTreeTool.cpp 2007-07-23 19:20:35 UTC (rev 1222) @@ -480,7 +480,107 @@ return MB_SUCCESS; } +/********************** General Sphere/Triangel Intersection ***************/ +MBErrorCode MBOrientedBoxTreeTool::sphere_intersect_triangles( + const double* center_v, + double radius, + MBEntityHandle tree_root, + std::vector<MBEntityHandle>& facets_out, + std::vector<MBEntityHandle>* sets_out ) +{ + const double radsqr = radius * radius; + MBOrientedBox b; + MBErrorCode rval; + MBRange sets; + const MBCartVect center(center_v); + MBCartVect closest, coords[3]; + const MBEntityHandle* conn; + int num_conn; +#ifndef MB_OBB_USE_VECTOR_QUERIES + MBRange tris; + MBRange::const_iterator t; +#else + std::vector<MBEntityHandle> tris; + std::vector<MBEntityHandle>::const_iterator t; +#endif + + std::vector<MBEntityHandle> stack, children; + stack.reserve(60); + stack.push_back(tree_root); + stack.push_back(0); + while (!stack.empty()) { + MBEntityHandle surf = stack.back(); stack.pop_back(); + MBEntityHandle node = stack.back(); stack.pop_back(); + if (!surf && sets_out) { + rval = get_moab_instance()->get_entities_by_type( node, MBENTITYSET, sets ); + if (!sets.empty()) + surf = sets.front(); + } + + // check if sphere intersects box + rval = box( node, b ); + if (MB_SUCCESS != rval) + return rval; + b.closest_location_in_box( center, closest ); + closest -= center; + if ((closest % closest) > radsqr) + continue; + + // push child boxes on stack + children.clear(); + rval = instance->get_child_meshsets( node, children ); + if (MB_SUCCESS != rval) + return rval; + if (!children.empty()) { + assert(children.size() == 2); + stack.push_back( children[0] ); + stack.push_back( surf ); + stack.push_back( children[1] ); + stack.push_back( surf ); + continue; + } + + // if leaf, intersect sphere with triangles +#ifndef MB_OBB_USE_VECTOR_QUERIES +# ifdef MB_OBB_USE_TYPE_QUERIES + rval = get_moab_instance()->get_entities_by_type( node, MBTRI, tris ); +# else + rval = get_moab_instance()->get_entities_by_handle( node, tris ); +# endif +#else + rval = get_moab_instance()->get_entities_by_handle( node, tris ); +#endif + if (MB_SUCCESS != rval) + return rval; + + for (t = tris.begin(); t != tris.end(); ++t) { + rval = get_moab_instance()->get_connectivity( *t, conn, num_conn, true ); + if (MB_SUCCESS != rval) + return rval; + if (num_conn != 3) + continue; + + rval = get_moab_instance()->get_coords( conn, num_conn, coords[0].array() ); + if (MB_SUCCESS != rval) + return rval; + + MBGeomUtil::closest_location_on_tri( center, coords, closest ); + closest -= center; + if ((closest % closest) <= radsqr && + std::find(facets_out.begin(),facets_out.end(),*t) == facets_out.end()) { + facets_out.push_back( *t ); + if (sets_out) + sets_out->push_back( surf ); + } + } + } + + return MB_SUCCESS; +} + + + /********************** General Ray/Tree and Ray/Triangel Intersection ***************/ @@ -843,53 +943,134 @@ /********************** Closest Point code ***************/ +struct MBOBBTreeCPFrame { + MBOBBTreeCPFrame( double d, MBEntityHandle n, MBEntityHandle s ) + : dist_sqr(d), node(n), mset(s) {} + double dist_sqr; + MBEntityHandle node; + MBEntityHandle mset; +}; + MBErrorCode MBOrientedBoxTreeTool::closest_to_location( const double* point, MBEntityHandle root, - double tolerance, double* point_out, MBEntityHandle& facet_out, MBEntityHandle* set_out ) { - std::vector<MBEntityHandle> facets(1), sets; - MBErrorCode rval = closest_to_location( point, root, tolerance, facets, set_out ? &sets : 0 ); - if (MB_SUCCESS != rval) - return rval; - if (facets.empty()) - return MB_ENTITY_NOT_FOUND; - facet_out = facets.front(); - if (set_out) - *set_out = sets.front(); + MBErrorCode rval; + const MBCartVect loc( point ); + double smallest_dist_sqr = std::numeric_limits<double>::max(); - const MBEntityHandle* conn; - int len; - rval = instance->get_connectivity( facet_out, conn, len, true ); - if (MB_SUCCESS != rval) - return rval; + MBEntityHandle current_set = 0; + MBRange sets; + std::vector<MBEntityHandle> children(2); + std::vector<double> coords; + std::vector<MBOBBTreeCPFrame> stack; + + stack.push_back( MBOBBTreeCPFrame(0.0, root, current_set) ); - std::vector<MBCartVect> coords( len ); - rval = instance->get_coords( conn, len, coords[0].array() ); - if (MB_SUCCESS != rval) - return rval; + while( !stack.empty() ) { + + // pop from top of stack + MBEntityHandle node = stack.back().node; + double dist_sqr = stack.back().dist_sqr; + current_set = stack.back().mset; + stack.pop_back(); + + // If current best result is closer than the box, skip this tree node. + if (dist_sqr > smallest_dist_sqr) + continue; + + // Check if this node has a set associated with it + if (set_out && !current_set) { + sets.clear(); + rval = instance->get_entities_by_type( node, MBENTITYSET, sets ); + if (MB_SUCCESS != rval) + return rval; + if (!sets.empty()) { + if (sets.size() != 1) + return MB_MULTIPLE_ENTITIES_FOUND; + current_set = sets.front(); + } + } + + // Get child boxes + children.clear(); + rval = instance->get_child_meshsets( node, children ); + if (MB_SUCCESS != rval) + return rval; + + // if not a leaf node + if (!children.empty()) { + if (children.size() != 2) + return MB_MULTIPLE_ENTITIES_FOUND; + + // get boxes + MBOrientedBox box1, box2; + rval = box( children[0], box1 ); + if (MB_SUCCESS != rval) return rval; + rval = box( children[1], box2 ); + if (MB_SUCCESS != rval) return rval; + + // get distance from each box + MBCartVect pt1, pt2; + box1.closest_location_in_box( loc, pt1 ); + box2.closest_location_in_box( loc, pt2 ); + pt1 -= loc; + pt2 -= loc; + const double dsqr1 = pt1 % pt1; + const double dsqr2 = pt2 % pt2; + + // push children on tree such that closer one is on top + if (dsqr1 < dsqr2) { + stack.push_back( MBOBBTreeCPFrame(dsqr2, children[1], current_set ) ); + stack.push_back( MBOBBTreeCPFrame(dsqr1, children[0], current_set ) ); + } + else { + stack.push_back( MBOBBTreeCPFrame(dsqr1, children[0], current_set ) ); + stack.push_back( MBOBBTreeCPFrame(dsqr2, children[1], current_set ) ); + } + } + else { // LEAF NODE + MBRange facets; + rval = instance->get_entities_by_dimension( node, 2, facets ); + if (MB_SUCCESS != rval) + return rval; + + const MBEntityHandle* conn; + int len; + MBCartVect tmp, diff; + for (MBRange::iterator i = facets.begin(); i != facets.end(); ++i) { + rval = instance->get_connectivity( *i, conn, len, true ); + if (MB_SUCCESS != rval) + return rval; + + coords.resize(3*len); + rval = instance->get_coords( conn, len, &coords[0] ); + if (MB_SUCCESS != rval) + return rval; + + if (len == 3) + MBGeomUtil::closest_location_on_tri( loc, (MBCartVect*)(&coords[0]), tmp ); + else + MBGeomUtil::closest_location_on_polygon( loc, (MBCartVect*)(&coords[0]), len, tmp ); + + diff = tmp - loc; + dist_sqr = diff % diff; + if (dist_sqr < smallest_dist_sqr) { + smallest_dist_sqr = dist_sqr; + facet_out = *i; + tmp.get( point_out ); + if (set_out) + *set_out = current_set; + } + } + } // LEAF NODE + } - if (len == 3) - MBGeomUtil::closest_location_on_tri( *reinterpret_cast<const MBCartVect*>(point), - &coords[0], - *reinterpret_cast<MBCartVect*>(point_out) ); - else - MBGeomUtil::closest_location_on_polygon( *reinterpret_cast<const MBCartVect*>(point), - &coords[0], len, - *reinterpret_cast<MBCartVect*>(point_out) ); return MB_SUCCESS; } - -struct MBOBBTreeCPFrame { - MBOBBTreeCPFrame( double d, MBEntityHandle n, MBEntityHandle s ) - : dist_sqr(d), node(n), mset(s) {} - double dist_sqr; - MBEntityHandle node; - MBEntityHandle mset; -}; MBErrorCode MBOrientedBoxTreeTool::closest_to_location( const double* point, MBEntityHandle root, Modified: MOAB/trunk/MBOrientedBoxTreeTool.hpp =================================================================== --- MOAB/trunk/MBOrientedBoxTreeTool.hpp 2007-07-20 18:24:11 UTC (rev 1221) +++ MOAB/trunk/MBOrientedBoxTreeTool.hpp 2007-07-23 19:20:35 UTC (rev 1222) @@ -203,7 +203,6 @@ */ MBErrorCode closest_to_location( const double* point, MBEntityHandle tree_root, - double tolerance, double* point_out, MBEntityHandle& facet_out, MBEntityHandle* set_out = 0); @@ -221,6 +220,23 @@ std::vector<MBEntityHandle>& facets_out, std::vector<MBEntityHandle>* sets_out = 0 ); + /**\brief Find facets intersected by a sphere + * + * Find facets intersected by a spherical volume. + *\param center Sphere center + *\param radius Sphere radius + *\param tree_root Root of OBB tree + *\param facets_out List of triangles intersecting sphere + *\param sets_out If not null, sets owning facets are appended to this + * list in an order corresponding to the entries in + * facets_out. + */ + MBErrorCode sphere_intersect_triangles( const double* center, + double radius, + MBEntityHandle tree_root, + std::vector<MBEntityHandle>& facets_out, + std::vector<MBEntityHandle>* sets_out = 0 ); + /**\brief Get oriented box at node in tree * * Get the oriented box for a node in an oriented bounding box tree. Modified: MOAB/trunk/test/obb/obb_test.cpp =================================================================== --- MOAB/trunk/test/obb/obb_test.cpp 2007-07-20 18:24:11 UTC (rev 1221) +++ MOAB/trunk/test/obb/obb_test.cpp 2007-07-23 19:20:35 UTC (rev 1222) @@ -1295,7 +1295,6 @@ // find closest point usnig tree rval = tool.closest_to_location( points[i].array(), root_set, - tolerance, t_result.array(), t_tri, set_ptr ); Modified: MOAB/trunk/tools/dagmc/DagMC.cpp =================================================================== --- MOAB/trunk/tools/dagmc/DagMC.cpp 2007-07-20 18:24:11 UTC (rev 1221) +++ MOAB/trunk/tools/dagmc/DagMC.cpp 2007-07-23 19:20:35 UTC (rev 1222) @@ -371,164 +371,60 @@ double x, double y, double z, int& result ) { + MBErrorCode rval; const double epsilon = moabMCNPTolerance; // Get OBB Tree for volume assert(volume - setOffset < rootSets.size()); MBEntityHandle root = rootSets[volume - setOffset]; - // Get closest triangles in volume + // Get closest point in triangulation + const MBCartVect point(x,y,z); + MBEntityHandle facet; + MBCartVect closest, diff; + rval = obbTree.closest_to_location( point.array(), root, closest.array(), facet ); + if (MB_SUCCESS != rval) return rval; + + // Check for on-boundary case + diff = closest - point; + if (diff%diff <= epsilon*epsilon) { + result = -1; + return MB_SUCCESS; + } + + // Get triangles at closest point 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 ); + rval = obbTree.sphere_intersect_triangles( closest.array(), epsilon, root, tris, &surfs ); if (MB_SUCCESS != rval) return rval; - // Get sense of each surface wrt volume - std::vector<int> senses( surfs.size() ); - rval = surface_sense( volume, surfs.size(), &surfs[0], &senses[0] ); - if (MB_SUCCESS != rval) return rval; - - // Get closest point on each triangle, and group triangles - // by intersection type - std::vector<int> edge_tris, vertex_tris; - - // Get corrected normal, closest point and topolgocal closest for each triangle - std::vector<MBCartVect> normals( tris.size() ), closest( tris.size() ); - std::vector<int> topo( tris.size() ); + // One-triangle case : determine if point is above or below triangle const MBEntityHandle* conn; - int len; - MBCartVect coords[3]; - for (unsigned i = 0; i < tris.size(); ++i) { - // volume tree shouldn't contain non-manifold surfaces - assert(senses[i]); - - rval = moab_instance()->get_connectivity( tris[i], conn, len ); + int len, sense; + MBCartVect coords[3], normal; + if (tris.size() == 1) { + rval = moab_instance()->get_connectivity( tris.front(), conn, len ); if (MB_SUCCESS != rval || len != 3) return MB_FAILURE; rval = moab_instance()->get_coords( conn, len, coords[0].array() ); - if (MB_SUCCESS != rval) - return rval; + if (MB_SUCCESS != rval) return rval; - MBGeomUtil::closest_location_on_tri( point, coords, epsilon, closest[i], topo[i] ); + rval = surface_sense( volume, 1, &surfs.front(), &sense ); + if (MB_SUCCESS != rval) return rval; + + // get triangle normal coords[1] -= coords[0]; coords[2] -= coords[0]; - normals[i] = senses[i] * coords[1] * coords[2]; - } - - // If within tolerance of triangle, then return that point is on boundary - if ((point - closest[0]).length() < epsilon) { - result = -1; + normal = sense * coords[1] * coords[2]; + // compare relative sense + result = (normal % (point - closest) < 0.0); return MB_SUCCESS; } - // if only one, then return - if (tris.size() == 1) { - // Ignore the possibility of being in the plane of the - // triangle. If parallel to one triangle, should have - // at least one other triangle. - result = (normals[0] % (closest[0] - point) > 0.0); - return MB_SUCCESS; - } + // many triangle case : determine if closest point is convexity or concavity - // If we found any triangle for which the closest point was in - // the interior, just use that one - std::vector<int>::iterator ii = std::find( topo.begin(), topo.end(), 6 ); - if (ii != topo.end()) { - int idx = ii - topo.begin(); - result = (normals[idx] % (closest[idx] - point) > 0.0); - return MB_SUCCESS; - } - // Otherwise if we found any triangle for which the - // closest point was at an edge, find it's neighbor and use the pair - for (ii = topo.begin(); ii != topo.end() && *ii < 3; ++ii); - if (ii != topo.end()) // found one - { - unsigned idx1 = ii - topo.end(); - int idx2 = -1; - MBEntityHandle edge[2]; - moab_instance()->get_connectivity( tris[idx1], conn, len, true ); - edge[0] = conn[*ii - 3]; - edge[1] = conn[(*ii - 2) % 3]; - for (unsigned i = 0; i < tris.size(); ++i) { - if (i == idx1) - continue; - moab_instance()->get_connectivity( tris[i], conn, len, true ); - int j; - for (j = 0; j < 3; ++j) - if (conn[j] == edge[0]) - if (conn[(j+1)%3] == edge[1] || conn[(j+2)%3] == edge[1]) - break; - if (j < 3) { - idx2 = i; - break; - } - } - - if (idx2 < 0) { // if no adjacent triangle - assert(false); - result = (normals[idx1] % (closest[idx1] - point) > 0.0); - return MB_SUCCESS; - } - - // Check if inside or outside of both facets first. - // This avoids potential rounding/accuracy issues - // when testing for convexity/concavity if the facets - // are near co-planer, and handles the special case - // where both facets are exactly coplanar. - const double dot1 = normals[idx1] % (closest[idx1] - point); - const double dot2 = normals[idx2] % (closest[idx2] - point); - const bool inside1 = dot1 > -std::numeric_limits<double>::epsilon(); - const bool inside2 = dot2 > -std::numeric_limits<double>::epsilon(); - const bool outside1 = dot1 < std::numeric_limits<double>::epsilon(); - const bool outside2 = dot2 < std::numeric_limits<double>::epsilon(); - if (inside1 && inside2) - result = 1; - else if (outside1 && outside2) - result = 0; - else { - // Edge is at a concavity if the projection of the normal of - // triangle A into the plane of triangle B points into triangle - // B. It is a convexity if it points out of triangle B (from the - // shared edge). This test can be reduced to checking the sign of - // dot product of the normal of triangle A with any vector pointing - // into triangle B from the shared edge. - - // Get triangle coordinates - moab_instance()->get_connectivity( tris[idx2], conn, len, true ); - moab_instance()->get_coords( conn, len, coords[0].array() ); - // Get a vector from edge end point to opposite vertex - const MBCartVect vect2 = coords[(topo[idx2] - 1) % 3] - coords[topo[idx2]]; - result = (vect2 % normals[idx1] >= 0.0); // true if a concavity - result = inside1 || inside2; - } - return MB_SUCCESS; - } - - // 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; - } - } - tris.resize(w); - surfs.resize(w); - normals.resize(w); - // use algorithm from: // (referance) // to determine inside vs. outside. @@ -536,20 +432,40 @@ // 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; + + rval = moab_instance()->get_connectivity( tris[i], conn, len ); + if (MB_SUCCESS != rval || len != 3) + return MB_FAILURE; + rval = moab_instance()->get_coords( conn, len, coords[0].array() ); + if (MB_SUCCESS != rval) return rval; + + rval = surface_sense( volume, 1, &surfs.front(), &sense ); + if (MB_SUCCESS != rval) return rval; + + // get triangle normal + coords[1] -= coords[0]; + coords[2] -= coords[0]; + normal = sense * coords[1] * coords[2]; + closest = coords[0]; + const double norm_len_sqr = normal % normal; + for (unsigned j = 0; j < tris.size(); ++j) { if (j == i) continue; - moab_instance()->get_connectivity( tris[j], conn, len ); + rval = moab_instance()->get_connectivity( tris[j], conn, len ); + if (MB_SUCCESS != rval || len != 3) + return MB_FAILURE; + rval = moab_instance()->get_coords( conn, len, coords[0].array() ); + if (MB_SUCCESS != rval) + return rval; + 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) + const double dot = normal % (coords[k] - closest); + if (dot * dot / norm_len_sqr > epsilon * epsilon) { + if (dot > 0.0) some_above = true; else some_below = true; @@ -560,7 +476,7 @@ // 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; + result = (normal % (closest - point)) > 0.0; return MB_SUCCESS; } // All other triangles to one side of this triangle: @@ -586,9 +502,6 @@ // detemine distance to nearest surface MBErrorCode DagMC::closest_to_location( MBEntityHandle volume, double* coords, double& result) { - - const double epsilon = moabMCNPTolerance; - // Get OBB Tree for volume assert(volume - setOffset < rootSets.size()); MBEntityHandle root = rootSets[volume - setOffset]; @@ -597,7 +510,7 @@ const MBCartVect point(coords); MBCartVect nearest; MBEntityHandle facet_out; - MBErrorCode rval = obbTree.closest_to_location( point.array(), root, epsilon, nearest.array(), facet_out ); + MBErrorCode rval = obbTree.closest_to_location( point.array(), root, nearest.array(), facet_out ); if (MB_SUCCESS != rval) return rval; // calculate distance between point and nearest facet

1 0

r1221 - MOAB/trunk/tools/dagmc
by tautges＠mcs.anl.gov 20 Jul '07

20 Jul '07

Author: tautges Date: 2007-07-20 13:24:11 -0500 (Fri, 20 Jul 2007) New Revision: 1221 Modified: MOAB/trunk/tools/dagmc/DagMC.cpp Log: Initializing normal to 0.0, from Paul Wilson. Modified: MOAB/trunk/tools/dagmc/DagMC.cpp =================================================================== --- MOAB/trunk/tools/dagmc/DagMC.cpp 2007-07-19 20:38:44 UTC (rev 1220) +++ MOAB/trunk/tools/dagmc/DagMC.cpp 2007-07-20 18:24:11 UTC (rev 1221) @@ -1285,7 +1285,7 @@ MBErrorCode rval = obbTree.closest_to_location( in_pt, root, tolerance(), facets ); assert(MB_SUCCESS == rval); - MBCartVect coords[3], normal; + MBCartVect coords[3], normal(0.0); const MBEntityHandle *conn; int len; for (unsigned i = 0; i < facets.size(); ++i) {

1 0

r1216 - MOAB/trunk
by kraftche＠mcs.anl.gov 18 Jul '07

18 Jul '07

Author: kraftche Date: 2007-07-18 17:55:21 -0500 (Wed, 18 Jul 2007) New Revision: 1216 Modified: MOAB/trunk/ReadHDF5.cpp Log: Fix really stupid memory leak: I added a wrapper class around raw arrays to take care of deleting the arrays, but forgot to give the class a destructor to do the actual deletion Modified: MOAB/trunk/ReadHDF5.cpp =================================================================== --- MOAB/trunk/ReadHDF5.cpp 2007-07-18 19:58:03 UTC (rev 1215) +++ MOAB/trunk/ReadHDF5.cpp 2007-07-18 22:55:21 UTC (rev 1216) @@ -510,6 +510,7 @@ public: auto_array( size_t s ) : data(new T[s]) {} auto_array() : data(0) {} + ~auto_array() { delete [] data; } T* get () { return data; } const T* get () const { return data; } T& operator[]( size_t i ) { return data[i]; }

1 0

r1215 - MOAB/trunk
by kraftche＠mcs.anl.gov 18 Jul '07

18 Jul '07

Author: kraftche Date: 2007-07-18 14:58:03 -0500 (Wed, 18 Jul 2007) New Revision: 1215 Modified: MOAB/trunk/VtkTest.cpp Log: add unit test for writing subset of mesh to vtk file Modified: MOAB/trunk/VtkTest.cpp =================================================================== --- MOAB/trunk/VtkTest.cpp 2007-07-18 19:23:46 UTC (rev 1214) +++ MOAB/trunk/VtkTest.cpp 2007-07-18 19:58:03 UTC (rev 1215) @@ -94,6 +94,8 @@ DECLARE_TEST(tensor_attrib_float); DECLARE_TEST(tensor_attrib_double); +DECLARE_TEST(subset); + int main( int argc, char* argv[] ) { int *test_indices = (int*)malloc( sizeof(int) * num_tests ); @@ -1246,3 +1248,100 @@ return check_tag_data( file, mb_type, 9 ); } + +bool test_subset( ) +{ + MBCore moab_inst; + MBInterface& moab = moab_inst; + MBErrorCode rval; + + // create 9 nodes in grid pattern + MBEntityHandle verts[9]; + const double coords[][3] = { { 0, 0, 0 }, + { 1, 0, 0 }, + { 2, 0, 0 }, + { 0, 1, 0 }, + { 1, 1, 0 }, + { 2, 1, 0 }, + { 0, 2, 0 }, + { 1, 2, 0 }, + { 2, 2, 0 } }; + for (unsigned i = 0; i < 9; ++i) { + rval = moab.create_vertex(coords[i], verts[i]); + assert(MB_SUCCESS == rval); + } + + // create 4 quad elements in grid pattern + const int conn[][4] = { { 0, 1, 4, 3 }, + { 1, 2, 5, 4 }, + { 3, 4, 7, 6 }, + { 4, 5, 8, 7 } }; + MBEntityHandle econn[4], elems[4]; + for (unsigned i = 0; i < 4; ++i) { + for (unsigned j = 0; j < 4; ++j) + econn[j] = verts[conn[i][j]]; + rval = moab.create_element( MBQUAD, econn, 4, elems[i] ); + assert(MB_SUCCESS == rval); + } + + // create 3 meshsets + MBEntityHandle sets[3]; + for (unsigned i = 0;i < 3; ++i) { + rval = moab.create_meshset( 0, sets[i] ); + assert(MB_SUCCESS == rval); + } + + // add element 3 to set 0 + rval = moab.add_entities( sets[0], elems+3, 1 ); + assert(MB_SUCCESS == rval); + // add node 2 to set 1 + rval = moab.add_entities( sets[1], verts+2, 1 ); + assert(MB_SUCCESS == rval); + // add element 2 and 3 to set 2 + rval = moab.add_entities( sets[2], elems+2, 2 ); + assert(MB_SUCCESS == rval); + + // make set 2 a child of set 1 + rval = moab.add_child_meshset( sets[1], sets[2] ); + assert(MB_SUCCESS == rval); + // put set 1 in set 0 + rval = moab.add_entities( sets[0], sets+1, 1 ); + assert(MB_SUCCESS == rval); + + // write sets[0] to vtk file + rval = moab.write_mesh( "tmp_file.vtk", sets, 1 ); + MBCHECK(rval); + + // read data back in + moab.delete_mesh(); + rval = moab.load_mesh( "tmp_file.vtk" ); + remove( "tmp_file.vtk" ); + MBCHECK(rval); + + // writer should have written all three sets, + // so the resulting mesh should be elems[2], elems[3], + // and verts[2] + MBRange new_elems, new_verts; + rval = moab.get_entities_by_type( 0, MBQUAD, new_elems ); + MBCHECK(rval); + CHECK( new_elems.size() == 2 ); + rval = moab.get_entities_by_type( 0, MBVERTEX, new_verts ); + MBCHECK(rval); + CHECK( new_verts.size() == 7 ); + + // vertex not in element closure should have coords of 2,0,0 + MBRange elem_verts; + rval = moab.get_adjacencies( new_elems, 0, false, elem_verts, MBInterface::UNION ); + MBCHECK(rval); + CHECK(elem_verts.size() == 6); + MBRange free_verts( new_verts.subtract(elem_verts ) ); + CHECK(free_verts.size() == 1 ); + double vcoords[3]; + rval = moab.get_coords( free_verts, vcoords ); + CHECK( vcoords[0] == 2 ); + CHECK( vcoords[1] == 0 ); + CHECK( vcoords[2] == 0 ); + + return true; +} +

1 0

r1214 - MOAB/trunk
by kraftche＠mcs.anl.gov 18 Jul '07

18 Jul '07

Author: kraftche Date: 2007-07-18 14:23:46 -0500 (Wed, 18 Jul 2007) New Revision: 1214 Modified: MOAB/trunk/WriteVtk.cpp Log: fix inefficiency in previous commit : closure of elements code belongs outside entityset loop Modified: MOAB/trunk/WriteVtk.cpp =================================================================== --- MOAB/trunk/WriteVtk.cpp 2007-07-18 19:14:54 UTC (rev 1213) +++ MOAB/trunk/WriteVtk.cpp 2007-07-18 19:23:46 UTC (rev 1214) @@ -171,22 +171,22 @@ nodes.merge( node_i, elem_i ); elems.merge( elem_i, set_i ); std::copy( set_i, a.end(), std::back_inserter(sets) ); - - for (MBRange::const_iterator e = elems.begin(); e != elems.end(); ++e) - { - const MBEntityHandle* conn; - int conn_len; - rval = mbImpl->get_connectivity( *e, conn, conn_len ); - if (MB_SUCCESS != rval) return rval; - - for (int i = 0; i < conn_len; ++i) - nodes.insert( conn[i] ); - } a.clear(); rval = mbImpl->get_child_meshsets( set, a ); std::copy( a.begin(), a.end(), std::back_inserter(sets) ); } + + for (MBRange::const_iterator e = elems.begin(); e != elems.end(); ++e) + { + const MBEntityHandle* conn; + int conn_len; + rval = mbImpl->get_connectivity( *e, conn, conn_len ); + if (MB_SUCCESS != rval) return rval; + + for (int i = 0; i < conn_len; ++i) + nodes.insert( conn[i] ); + } } if (nodes.empty())

1 0

r1213 - MOAB/trunk
by kraftche＠mcs.anl.gov 18 Jul '07

18 Jul '07

Author: kraftche Date: 2007-07-18 14:14:54 -0500 (Wed, 18 Jul 2007) New Revision: 1213 Modified: MOAB/trunk/WriteVtk.cpp Log: When WriteVtk is passed one or more entity sets to write, write all entities in the closure of those sets (contained and child sets), not just entities contained directly Modified: MOAB/trunk/WriteVtk.cpp =================================================================== --- MOAB/trunk/WriteVtk.cpp 2007-07-18 19:13:45 UTC (rev 1212) +++ MOAB/trunk/WriteVtk.cpp 2007-07-18 19:14:54 UTC (rev 1213) @@ -30,6 +30,9 @@ #include <iostream> #include <stdio.h> #include <assert.h> +#include <vector> +#include <set> +#include <iterator> #include "MBInterface.hpp" #include "MBRange.hpp" @@ -145,9 +148,18 @@ } else { - for (int i = 0; i < num_sets; ++i) - { - MBEntityHandle set = set_list[i]; + std::set<MBEntityHandle> visited; + std::vector<MBEntityHandle> sets; + sets.reserve(num_sets); + std::copy( set_list, set_list + num_sets, std::back_inserter(sets) ); + while (!sets.empty()) { + // get next set + MBEntityHandle set = sets.back(); + sets.pop_back(); + // skip sets we've already done + if (!visited.insert(set).second) + continue; + MBRange a; rval = mbImpl->get_entities_by_handle( set, a ); if (MB_SUCCESS != rval) return rval; @@ -158,17 +170,22 @@ set_i = a.lower_bound( elem_i, a.end(), CREATE_HANDLE( MBENTITYSET, 0, e ) ); nodes.merge( node_i, elem_i ); elems.merge( elem_i, set_i ); - } + std::copy( set_i, a.end(), std::back_inserter(sets) ); - for (MBRange::const_iterator e = elems.begin(); e != elems.end(); ++e) - { - const MBEntityHandle* conn; - int conn_len; - rval = mbImpl->get_connectivity( *e, conn, conn_len ); - if (MB_SUCCESS != rval) return rval; + for (MBRange::const_iterator e = elems.begin(); e != elems.end(); ++e) + { + const MBEntityHandle* conn; + int conn_len; + rval = mbImpl->get_connectivity( *e, conn, conn_len ); + if (MB_SUCCESS != rval) return rval; + + for (int i = 0; i < conn_len; ++i) + nodes.insert( conn[i] ); + } - for (int i = 0; i < conn_len; ++i) - nodes.insert( conn[i] ); + a.clear(); + rval = mbImpl->get_child_meshsets( set, a ); + std::copy( a.begin(), a.end(), std::back_inserter(sets) ); } }

1 0

r1212 - MOAB/trunk
by kraftche＠mcs.anl.gov 18 Jul '07

18 Jul '07

Author: kraftche Date: 2007-07-18 14:13:45 -0500 (Wed, 18 Jul 2007) New Revision: 1212 Modified: MOAB/trunk/MBCore.cpp Log: Fix bug in function for blocked coordinate access that I added on May 15, 2007 Modified: MOAB/trunk/MBCore.cpp =================================================================== --- MOAB/trunk/MBCore.cpp 2007-07-18 18:25:05 UTC (rev 1211) +++ MOAB/trunk/MBCore.cpp 2007-07-18 19:13:45 UTC (rev 1212) @@ -574,7 +574,7 @@ if (*iter > end_ent) // trim block to those in the sequence iter -= *iter - end_ent; - const MBEntityID offset = start_ent - range_start; + const MBEntityID offset = range_start - start_ent; const MBEntityID count = *iter - range_start + 1; ++iter; // done with iter for this iteration, advance to next handle

1 0

r1210 - MOAB/trunk/tools/dagmc
by kraftche＠mcs.anl.gov 17 Jul '07

17 Jul '07

Author: kraftche Date: 2007-07-17 16:22:01 -0500 (Tue, 17 Jul 2007) New Revision: 1210 Modified: MOAB/trunk/tools/dagmc/DagMC.cpp MOAB/trunk/tools/dagmc/Makefile.am MOAB/trunk/tools/dagmc/cgm2moab.cc MOAB/trunk/tools/dagmc/main.cc Log: o Change DagMC and other tools to use new method for initializing ACIS in CGM - now works when using libicubitNN.so w/out CGM ACIS source or headers. o Fix build of test_geom w/ CGM support enabled Modified: MOAB/trunk/tools/dagmc/DagMC.cpp =================================================================== --- MOAB/trunk/tools/dagmc/DagMC.cpp 2007-07-17 21:20:35 UTC (rev 1209) +++ MOAB/trunk/tools/dagmc/DagMC.cpp 2007-07-17 21:22:01 UTC (rev 1210) @@ -6,9 +6,8 @@ #include "MBGeomUtil.hpp" #ifdef CGM +#include "InitCGMA.hpp" #include "CGMApp.hpp" -#include "AcisQueryEngine.hpp" -#include "AcisModifyEngine.hpp" #include "CubitDefines.h" #include "GeometryQueryTool.hpp" #include "CubitVector.hpp" @@ -762,9 +761,8 @@ #ifdef CGM // initialize cgm - CGMApp::instance()->startup( 0, NULL ); - AcisQueryEngine::instance(); - AcisModifyEngine::instance(); + InitCGMA::initialize_cgma(); + InitCGMA::initialize_engine("ACIS"); CGMApp::instance()->attrib_manager()->set_all_auto_read_flags(true); CGMApp::instance()->attrib_manager()->set_all_auto_actuate_flags(true); CubitStatus s = GeometryQueryTool::instance()->import_solid_model( scfile.c_str(), "ACIS_SAT"); Modified: MOAB/trunk/tools/dagmc/Makefile.am =================================================================== --- MOAB/trunk/tools/dagmc/Makefile.am 2007-07-17 21:20:35 UTC (rev 1209) +++ MOAB/trunk/tools/dagmc/Makefile.am 2007-07-17 21:22:01 UTC (rev 1210) @@ -36,6 +36,6 @@ TESTS = test_geom check_PROGRAMS = $(TESTS) test_geom_SOURCES = test_geom.cc -test_geom_LDADD = libdagmc.la $(top_builddir)/libMOAB.la +test_geom_LDADD = libdagmc.la $(top_builddir)/libMOAB.la $(CGM_LIBS_LTFLAGS) $(CGM_LIBS_LINK) Modified: MOAB/trunk/tools/dagmc/cgm2moab.cc =================================================================== --- MOAB/trunk/tools/dagmc/cgm2moab.cc 2007-07-17 21:20:35 UTC (rev 1209) +++ MOAB/trunk/tools/dagmc/cgm2moab.cc 2007-07-17 21:22:01 UTC (rev 1210) @@ -1,6 +1,3 @@ -#include "AcisQueryEngine.hpp" -#include "AcisModifyEngine.hpp" -#include "CGMApp.hpp" #include "GeometryQueryTool.hpp" #include "GeometryQueryEngine.hpp" #include "ModelQueryEngine.hpp" Modified: MOAB/trunk/tools/dagmc/main.cc =================================================================== --- MOAB/trunk/tools/dagmc/main.cc 2007-07-17 21:20:35 UTC (rev 1209) +++ MOAB/trunk/tools/dagmc/main.cc 2007-07-17 21:22:01 UTC (rev 1210) @@ -1,6 +1,5 @@ #include "GeometryQueryTool.hpp" -#include "AcisQueryEngine.hpp" -#include "AcisModifyEngine.hpp" +#include "InitCGMA.hpp" #include "CGMApp.hpp" #include "MBCore.hpp" #include "cgm2moab.hpp" @@ -114,9 +113,8 @@ // Initialize CGM - CGMApp::instance()->startup( 0, NULL ); - AcisQueryEngine::instance(); - AcisModifyEngine::instance(); + InitCGMA::initialize_cgma(); + InitCGMA::initialize_engine("ACIS"); if (actuate_attribs) { CGMApp::instance()->attrib_manager()->set_all_auto_read_flags( actuate_attribs ); CGMApp::instance()->attrib_manager()->set_all_auto_actuate_flags( actuate_attribs );

1 0