[MOAB-dev] commit/MOAB: 14 new changesets
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Mon Feb 10 14:08:41 CST 2014
14 new commits in MOAB:
https://bitbucket.org/fathomteam/moab/commits/d371b7209156/
Changeset: d371b7209156
Branch: None
User: pshriwise
Date: 2014-01-06 19:10:37
Summary: Merged fathomteam/moab into master
Affected #: 287 files
diff --git a/.gitignore b/.gitignore
index 6d29b5e..e824a6b 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,99 +1,142 @@
-MOABConfig.cmake
-moab.files
-config/test-driver
-moab.config
-moab.includes
-moab.creator*
-*.ccm
-*.cub
+*~
+*.a
aclocal.m4
autom4te.cache/
-config.h
-config.h.in
-config.log
-config.lt
-config.status
+bin
+bin/*
+*.ccm
config/config.guess
config/config.sub
config/depcomp
+config.h
+config.h.in
config/install-sh
config/libtool.m4
+config.log
+config.lt
config/ltmain.sh
+config/lt~obsolete.m4
config/ltoptions.m4
config/ltsugar.m4
config/ltversion.m4
-config/lt~obsolete.m4
config/missing
+config.status
+config/test-driver
configure
+.cproject
+*.cub
+cubit??.jou
+.deps
doc/config.tex
doc/dev.dox
-doc/user.dox
doc/user/*
+doc/user.dox
examples/examples.make
+examples/DeformMeshRemap
+examples/DirectAccessNoHoles
+examples/DirectAccessWithHoles
+examples/FileRead
+examples/GeomSetHierarchy
+examples/GetEntities
+examples/*.h5m
+examples/HelloMOAB
+examples/HelloParMOAB
+examples/itaps/FindConnectF
+examples/itaps/ListSetsNTagsCXX
+examples/itaps/ListSetsNTagsF90
+examples/itaps/TagIterateC
+examples/itaps/TagIterateF
+examples/KDTree
+examples/ObbTree
+examples/PushParMeshIntoMoabF90
+examples/ReduceExchangeTags
+examples/SetsNTags
+examples/SkinMesh
+examples/StructuredMeshSimple
+examples/SurfArea
+examples/TestExodusII
+history??.jou
+hpctoolkit-*
+*.hpcstruct
+include
+include/*
itaps/iBase_f.h
itaps/igeom/FBiGeom-Defs.inc
+itaps/igeom/FBiGeom_protos.h
+itaps/igeom/testgeom
+itaps/igeom/testSmooth2
+itaps/igeom/testSmoothGeom
+itaps/imesh/FindAdjacencyF90
itaps/imesh/iMesh-Defs.inc
+itaps/imesh/iMesh_extensions_protos.h
itaps/imesh/iMeshP_extensions_protos.h
itaps/imesh/iMeshP_protos.h
-itaps/imesh/iMesh_extensions_protos.h
itaps/imesh/iMesh_protos.h
+itaps/imesh/MOAB_iMesh_extensions_tests
+itaps/imesh/MOAB_iMeshP_unit_tests
+itaps/imesh/MOAB_iMesh_unit_tests
+itaps/imesh/partest
+itaps/imesh/ScdMeshF77
+itaps/imesh/ScdMeshF90
+itaps/imesh/testc_cbind
+*.la
+*.la
+*.lai
+lib
+lib/*
+.libs
libtool
+*.lo
+*.log
+makefile
+Makefile
+*/Makefile
+*/**/Makefile
+*/**/Makefile.in
+*/Makefile.in
+Makefile.in
+moab.config
+MOABConfig.cmake
+moab.creator*
+moab.files
+moab.includes
moab.make
+*.o
+.project
+*.rej
+share/*
+share/doc/moab
+share/man/man1
+*.so
src/FCDefs.h
+src/io/mhdf/h5minfo
+src/io/mhdf/h5mvalidate
src/MBCN_protos.h
-src/MOAB_FCDefs.h
src/moab/EntityHandle.hpp
-src/moab/Version.h
+src/MOAB_FCDefs.h
src/moab/stamp-h2
src/moab/stamp-h3
+src/moab/Version.h
src/parallel/moab_mpi_config.h
src/parallel/stamp-h4
src/stamp-h5
stamp-h1
-tools/mbcoupler/tests/
-tools/mbzoltan/Config.moab
-tools/vtkMOABReader/CMakeLists.txt
-tools/vtkMOABReaderNew/CMakeLists.txt
-.deps
-Makefile.in
-Makefile
-*/Makefile.in
-*/Makefile
-*/**/Makefile.in
-*/**/Makefile
-.libs
-*.o
-*.log
-*.lo
-*.la
-*.a
-*.so
-include/*
-lib/*
-share/*
-bin/*
-*~
-examples/HelloMoabPar
-examples/TestExodusII
-itaps/igeom/FBiGeom_protos.h
-itaps/igeom/testSmooth2
-itaps/igeom/testSmoothGeom
-itaps/igeom/testgeom
-itaps/imesh/FindAdjacencyF90
-itaps/imesh/MOAB_iMeshP_unit_tests
-itaps/imesh/ScdMeshF77
-itaps/imesh/ScdMeshF90
-itaps/imesh/partest
test/adaptive_kd_tree_tests
test/bsp_tree_poly_test
test/bsp_tree_test
-test/*.gen
+test/CMakeLists.txt
test/coords_connect_iterate
test/cropvol_test
test/dual/dual_test
+test/elem_eval_test
test/file_options_test
+test/*.g
+test/*.g
+test/*.gen
+test/*.gen
test/geom_util_test
test/gttool_test
+test/gttool_test
test/h5file/dump_sets
test/h5file/h5legacy
test/h5file/h5partial
@@ -102,13 +145,14 @@ test/h5file/h5regression
test/h5file/h5sets_test
test/h5file/h5test
test/h5file/h5varlen
-test/*.g
test/homxform_test
+test/io/*.ccmg
+test/io/ccmio_test
test/io/cub_file_test
test/io/exodus_test
+test/io/*.g
test/io/gmsh_test
test/io/ideas_test
-test/io/*.g
test/io/nastran_test
test/io/read_cgm_test
test/io/read_nc
@@ -131,15 +175,18 @@ test/obb/obb_time
test/obb/obb_tree_tool
test/obb_test
test/oldinc/test_oldinc
+test/read_mpas_nc
test/parallel/*.h5m
-test/parallel/*.vtk
test/parallel/mbparallelcomm_test
test/parallel/mhdf_parallel
+test/parallel/mpastrvpart
test/parallel/par_coupler_test
test/parallel/par_intx_sph
test/parallel/parallel_hdf5_test
test/parallel/parallel_unit_tests
test/parallel/parallel_write_test
+test/parallel/par_coupler_test
+test/parallel/par_intx_sph
test/parallel/parmerge
test/parallel/partcheck
test/parallel/pcomm_serial
@@ -150,38 +197,68 @@ test/parallel/scdtest
test/parallel/structured3
test/parallel/uber_parallel_test
test/parallel/ucdtrvpart
+test/parallel/*.vtk
test/perf/adj_time
test/perf/perf
test/perf/perftool
+test/perf/point_in_elem
+test/perf/runtest
test/perf/seqperf
test/perf/tstt_perf_binding
+test/perf/point_location/elem_eval_time
+test/perf/point_location/point_location
test/range_test
test/reorder_test
test/scdseq_test
+test/scd_test_partn
test/seq_man_test
+test/spatial_locator_test
+test/test_boundbox
test/tag_test
test/test_adj
test/test_prog_opt
+test/TestRunner.hpp
test/var_len_test
test/var_len_test_no_template
test/xform_test
+tools/dagmc/dagmc_preproc
tools/dagmc/pt_vol_test
+tools/dagmc/quads_to_tris
tools/dagmc/ray_fire_test
tools/dagmc/test_geom
tools/dagmc/update_coords
+tools/hexmodops
+tools/mbconvert
tools/mbcoupler/*.h5m
+tools/mbcoupler/*.g
+tools/mbcoupler/tests/
tools/mbcslam/case1_test
tools/mbcslam/intersect1.h5m
-tools/mbcslam/intx.vtk
tools/mbcslam/intx1.vtk
tools/mbcslam/intx_in_plane_test
tools/mbcslam/intx_on_sphere_test
+tools/mbcslam/intx.vtk
tools/mbcslam/lagr.h5m
-tools/mbcslam/spec_visu_test
tools/mbcslam/spectral.vtk
+tools/mbcslam/spec_visu_test
tools/mbcslam/spherical_area_test
-.project
-.cproject
-examples/*.h5m
-examples/ReduceExchangeTags
-
+tools/mbdepth
+tools/mbgsets
+tools/mbmem
+tools/mbsize
+tools/mbskin
+tools/mbsurfplot
+tools/mbtagprop
+tools/mbzoltan/Config.moab
+tools/spheredecomp
+tools/vtkMOABReader/CMakeLists.txt
+tools/vtkMOABReaderNew/CMakeLists.txt
+tools/vtkMOABReaderNew/CMakeCache.txt
+tools/vtkMOABReaderNew/CMakeFiles/*
+tools/vtkMOABReaderNew/cmake_install.cmake
+tools/vtkMOABReaderNew/vtkMoabReaderPlugin.qrc
+tools/vtkMOABReaderNew/vtkMoabReaderPluginInit.cxx
+tools/vtkMOABReaderNew/vtkMoabReaderPlugin_Plugin.cxx
+tools/vtkMOABReaderNew/vtkMoabReaderPlugin_Plugin.h
+tools/vtkMOABReaderNew/vtkSMvtkMoabReaderPluginInstantiator.cxx
+tools/vtkMOABReaderNew/vtkSMvtkMoabReaderPluginInstantiator.h
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 3774046..6b18435 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,3 +1,9 @@
+# check if we are using MPI - reset compilers accordingly
+if ( MOAB_USE_MPI )
+ SET(CMAKE_CXX_COMPILER ${MPI_CXX_COMPILER})
+ SET(CMAKE_C_COMPILER ${MPI_C_COMPILER})
+endif ( MOAB_USE_MPI )
+
project( MOAB )
cmake_minimum_required( VERSION 2.4 )
@@ -100,7 +106,30 @@ project( MOAB )
# iMesh
option ( MOAB_BUILD_IMESH "Build the iMesh interface?" ON )
+ # check for MPI package
+ if ( MOAB_USE_MPI )
+ find_package( MPI )
+ # CMake FindMPI script is sorely lacking:
+ if ( MPI_LIBRARY AND MPI_INCLUDE_PATH )
+ set( MPI_FOUND 1 )
+ endif ( MPI_LIBRARY AND MPI_INCLUDE_PATH )
+
+ if ( MPI_FOUND )
+ set ( MOAB_DEFINES "${MOAB_DEFINES} -DUSE_MPI" )
+ endif ( MPI_FOUND )
+ endif ( MOAB_USE_MPI )
+
+ if ( MOAB_USE_NETCDF )
+ find_package( NetCDF )
+ endif ( MOAB_USE_NETCDF )
+
+ if ( MOAB_USE_HDF )
+ # HDF5
+ find_package( HDF5 )
+ endif ( MOAB_USE_HDF )
+
add_subdirectory( src )
+ add_subdirectory( itaps/imesh )
add_subdirectory( tools )
add_subdirectory( test )
-
+
diff --git a/Makefile.am b/Makefile.am
index 1c91622..b8c3c9e 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -48,7 +48,7 @@ install-data-hook:
echo -n "SET(MOAB_LIBRARIES \"" >> $(DESTDIR)$(cfgdir)/MOABConfig.cmake
# next line: get the post-substitution set of link commands, then filter out leading
# and trailing spaces, then terminate with close-quote and close-paren
- echo -n "@EXPORT_LDFLAGS@ -L$(MOAB_LIBDIR) -lMOAB @LIBS@ @PNETCDF_LIBS@ @NETCDF_LIBS@ @HDF5_LIBS@ @CCMIO_LIBS@ @CGM_LIBS@" | \
+ echo -n "@EXPORT_LDFLAGS@ -L$(libdir) -lMOAB @LIBS@ @PNETCDF_LIBS@ @NETCDF_LIBS@ @CGNS_LIBS@ @HDF5_LIBS@ @CCMIO_LIBS@ @CGM_LIBS@" | \
sed -e 's/^ *//;s/ *$$/\")\n/' >> $(DESTDIR)$(cfgdir)/MOABConfig.cmake
# Generate a file to be installed in $libdir containing the configuration
diff --git a/MeshFiles/unittest/125hex.g b/MeshFiles/unittest/125hex.g
index 6cf0433..2042583 100644
Binary files a/MeshFiles/unittest/125hex.g and b/MeshFiles/unittest/125hex.g differ
diff --git a/MeshFiles/unittest/16_unmerged_hex.h5m b/MeshFiles/unittest/16_unmerged_hex.h5m
new file mode 100644
index 0000000..2ebde86
Binary files /dev/null and b/MeshFiles/unittest/16_unmerged_hex.h5m differ
diff --git a/MeshFiles/unittest/1hex.g b/MeshFiles/unittest/1hex.g
index 3dcd658..0affc6e 100644
Binary files a/MeshFiles/unittest/1hex.g and b/MeshFiles/unittest/1hex.g differ
diff --git a/MeshFiles/unittest/1khex.g b/MeshFiles/unittest/1khex.g
index e6c9417..6b0a8fc 100644
Binary files a/MeshFiles/unittest/1khex.g and b/MeshFiles/unittest/1khex.g differ
diff --git a/MeshFiles/unittest/64bricks_512hex_256part.h5m b/MeshFiles/unittest/64bricks_512hex_256part.h5m
new file mode 100644
index 0000000..d7df13c
Binary files /dev/null and b/MeshFiles/unittest/64bricks_512hex_256part.h5m differ
diff --git a/MeshFiles/unittest/Makefile.am b/MeshFiles/unittest/Makefile.am
index f347315..a1784ee 100644
--- a/MeshFiles/unittest/Makefile.am
+++ b/MeshFiles/unittest/Makefile.am
@@ -1,6 +1,7 @@
SUBDIRS = io
EXTRA_DIST = 125hex.g \
+ 16_unmerged_hex.h5m \
1hex.g \
1khex.g \
1tet.g \
@@ -21,6 +22,8 @@ EXTRA_DIST = 125hex.g \
testtet.g \
ptest.cub \
64bricks_512hex.h5m \
+ 64bricks_1khex.h5m\
+ 64bricks_12ktet.h5m\
partBed.smf \
PB.h5m \
polyPB.txt \
@@ -28,5 +31,7 @@ EXTRA_DIST = 125hex.g \
quads.h5m \
SECrop2.h5m \
poly14.txt \
- BedCrop2.h5m
+ BedCrop2.h5m \
+ mpas_p8.h5m \
+ Homme_2pt.h5m
diff --git a/MeshFiles/unittest/io/2d_naca0012.cgns b/MeshFiles/unittest/io/2d_naca0012.cgns
new file mode 100644
index 0000000..4ce0452
Binary files /dev/null and b/MeshFiles/unittest/io/2d_naca0012.cgns differ
diff --git a/MeshFiles/unittest/io/Makefile.am b/MeshFiles/unittest/io/Makefile.am
index 8257e5f..dbce393 100644
--- a/MeshFiles/unittest/io/Makefile.am
+++ b/MeshFiles/unittest/io/Makefile.am
@@ -7,6 +7,7 @@ EXTRA_DIST = HommeMapping.nc \
cubtest12.cub \
cubtest.jou \
dum.sat \
+ dum.stp \
gmsh2.msh \
hex_2x2x2_ss.exo \
ho_test.g \
@@ -18,6 +19,8 @@ EXTRA_DIST = HommeMapping.nc \
test.nas \
test.unv \
sample.stl \
- three.smf
+ three.smf \
+ mpasx1.642.t.2.nc \
+ 2d_naca0012.cgns
diff --git a/MeshFiles/unittest/io/mpasx1.642.t.2.nc b/MeshFiles/unittest/io/mpasx1.642.t.2.nc
new file mode 100644
index 0000000..f055c49
Binary files /dev/null and b/MeshFiles/unittest/io/mpasx1.642.t.2.nc differ
diff --git a/MeshFiles/unittest/mpas_p8.h5m b/MeshFiles/unittest/mpas_p8.h5m
new file mode 100644
index 0000000..75662d2
Binary files /dev/null and b/MeshFiles/unittest/mpas_p8.h5m differ
diff --git a/MeshFiles/unittest/surfrandomtris-4part.h5m b/MeshFiles/unittest/surfrandomtris-4part.h5m
new file mode 100644
index 0000000..1ca984a
Binary files /dev/null and b/MeshFiles/unittest/surfrandomtris-4part.h5m differ
diff --git a/config/cgns.m4 b/config/cgns.m4
new file mode 100644
index 0000000..16b19e3
--- /dev/null
+++ b/config/cgns.m4
@@ -0,0 +1,92 @@
+#######################################################################################
+# Check for CGNS library
+# Sets HAVE_CGNS to 'yes' or 'no'
+# If HAVE_CGNS == yes, then exports:
+# CGNS_CPPFLAGS
+# CGNS_LDFLAGS
+# CGNS_LIBS
+#######################################################################################
+AC_DEFUN([FATHOM_CHECK_CGNS],[
+
+AC_MSG_CHECKING([if CGNS support is enabled])
+AC_ARG_WITH(cgns,
+[AC_HELP_STRING([--with-cgns<:@=DIR@:>@], [Specify CGNS library to use for CGNS file format])
+AC_HELP_STRING([--without-cgns], [Disable support for CGNS file format])],
+[CGNS_ARG=$withval
+DISTCHECK_CONFIGURE_FLAGS="$DISTCHECK_CONFIGURE_FLAGS --with-cgns=\"${withval}\""
+]
+, [CGNS_ARG=])
+if test "xno" != "x$CGNS_ARG"; then
+ AC_MSG_RESULT([yes])
+else
+ AC_MSG_RESULT([no])
+fi
+
+ # if CGNS support is not disabled
+HAVE_CGNS=no
+if test "xno" != "x$CGNS_ARG"; then
+ HAVE_CGNS=yes
+
+ # if a path is specified, update LIBS and INCLUDES accordingly
+ if test "xyes" != "x$CGNS_ARG" && test "x" != "x$CGNS_ARG"; then
+ if test -d "${CGNS_ARG}/lib"; then
+ CGNS_LDFLAGS="-L${CGNS_ARG}/lib"
+ elif test -d "${CGNS_ARG}"; then
+ CGNS_LDFLAGS="-L${CGNS_ARG}"
+ else
+ AC_MSG_ERROR("$CGNS_ARG is not a directory.")
+ fi
+ if test -d "${CGNS_ARG}/include"; then
+ CGNS_CPPFLAGS="-I${CGNS_ARG}/include"
+ elif test -d "${CGNS_ARG}/inc"; then
+ CGNS_CPPFLAGS="-I${CGNS_ARG}/inc"
+ else
+ CGNS_CPPFLAGS="-I${CGNS_ARG}"
+ fi
+ fi
+
+ old_CPPFLAGS="$CPPFLAGS"
+ CPPFLAGS="$CGNS_CPPFLAGS $CPPFLAGS"
+ old_LDFLAGS="$LDFLAGS"
+ LDFLAGS="$CGNS_LDFLAGS $HDF5_LDFLAGS $LDFLAGS"
+
+ # Check for C library
+ AC_LANG_PUSH([C])
+
+ AC_CHECK_HEADER([cgnslib.h], [],
+ [AC_MSG_WARN([[CGNS header not found.]]); HAVE_CGNS=no] )
+
+ # Check if cgns is usable by itself
+ AC_CHECK_LIB( [cgns], [cg_open], [CGNS_LIBS="-lcgns"], [
+ # Check if cgns is usable with HDF5
+ unset ac_cv_lib_cgns
+ unset ac_cv_lib_cgns_cg_open
+ # If we haven't already looked for HDF5 libraries, again now incase
+ # they're in the CGNS lib directory.
+ FATHOM_DETECT_HDF5_LIBS
+ LDFLAGS="$LDFLAGS $HDF5_LDFLAGS"
+ AC_CHECK_LIB( [cgns], [cg_open], [CGNS_LIBS="-lcgns -lhdf5_hl"], [
+ # Try one more time with HDF5 and libcurl
+ unset ac_cv_lib_cgns
+ unset ac_cv_lib_cgns_cg_open
+ AC_CHECK_LIB( [cgns], [cg_open], [CGNS_LIBS="-lcgns -lhdf5_hl -lcurl"],
+ [HAVE_CGNS=no], [-lhdf5_hl $HDF5_LIBS -lcurl] )],
+ [-lhdf5_hl $HDF5_LIBS] )],
+ )
+
+ CPPFLAGS="$old_CPPFLAGS"
+ LDFLAGS="$old_LDFLAGS"
+ AC_LANG_POP([C])
+
+ if test "x$HAVE_CGNS" = "xno"; then
+ if test "x$CGNS_ARG" != "x"; then
+ AC_MSG_ERROR("CGNS not found or not working")
+ else
+ AC_MSG_WARN("CGNS support disabled")
+ fi
+ CGNS_CPPFLAGS=
+ CGNS_LDFLAGS=
+ fi
+fi
+
+]) # FATHOM_HAVE_CGNS
diff --git a/config/compiler.m4 b/config/compiler.m4
index 6a7e2ed..ce83a5f 100644
--- a/config/compiler.m4
+++ b/config/compiler.m4
@@ -222,6 +222,14 @@ if test "xyes" = "x$enable_debug"; then
CFLAGS="$CFLAGS -g"
FCFLAGS="$FCFLAGS -g"
FFLAGS="$FFLAGS -g"
+ # Add -fstack-protector-all option for g++ in debug mode
+ if test "x$cxx_compiler" = "xGNU"; then
+ CXXFLAGS="$CXXFLAGS -fstack-protector-all"
+ fi
+ # Add -fstack-protector-all option for gcc in debug mode
+ if test "x$cc_compiler" = "xGNU"; then
+ CFLAGS="$CFLAGS -fstack-protector-all"
+ fi
fi
if test "xyes" = "x$enable_cxx_optimize"; then
CXXFLAGS="$CXXFLAGS -O2 -DNDEBUG"
diff --git a/configure.ac b/configure.ac
index 130e4a9..ea22a5e 100644
--- a/configure.ac
+++ b/configure.ac
@@ -446,26 +446,45 @@ AM_CPPFLAGS="$PNETCDF_CPPFLAGS $AM_CPPFLAGS"
EXPORT_LDFLAGS="$PNETCDF_LDFLAGS $EXPORT_LDFLAGS"
AC_SUBST(PNETCDF_LIBS)
+################################################################################
+# CGNS OPTIONS
+################################################################################
+
+old_LDFLAGS="$LDFLAGS"
+LDFLAGS="$LDFLAGS $HDF5_LDFLAGS"
+FATHOM_CHECK_CGNS
+LDFLAGS="$old_LDFLAGS"
+if test "xno" != "x$HAVE_CGNS"; then
+ AM_CPPFLAGS="$AM_CPPFLAGS -DCGNS_FILE"
+fi
+AM_CONDITIONAL(CGNS_FILE, [test "xno" != "x$HAVE_CGNS"])
+AM_CPPFLAGS="$CGNS_CPPFLAGS $AM_CPPFLAGS"
+EXPORT_LDFLAGS="$CGNS_LDFLAGS $EXPORT_LDFLAGS"
+AC_SUBST(CGNS_LIBS)
+
#################################################################################
# Documentation
#################################################################################
-AC_ARG_ENABLE([doxygen],
-[AC_HELP_STRING([[--enable-doxygen@<:@=DIR@:>@]],[Specify directory where Doxygen program is installed])
-AC_HELP_STRING([--disable-doxygen],[Do not generate API documentation (default)])],
- [ENABLE_DOXYGEN="$enableval"],[ENABLE_DOXYGEN=no] )
-if test "x$ENABLE_DOXYGEN" = "xyes"; then
- AC_PATH_PROG( [DOXYGEN], [doxygen], [no] )
-elif test "x$ENABLE_DOXYGEN" != "xno"; then
- AC_PATH_PROG( [DOXYGEN], [doxygen], [no], [$ENABLE_DOXYGEN] )
+AC_ARG_ENABLE([docs],
+[AC_HELP_STRING([--enable-docs],[indicate to check for doxygen installation])],
+ [ENABLE_DOCS=yes],[ENABLE_DOCS=no] )
+AC_ARG_WITH([doxygen],
+[AC_HELP_STRING([--with-doxygen=DIR], [Specify directory where Doxygen program is installed. By default, it is /usr/bin])],
+[WITH_DOXYGEN="$withval"], [WITH_DOXYGEN=no])
+
+if test "x$WITH_DOXYGEN" != "xno"; then
+ AC_PATH_PROGS( [DOXYGEN], [doxygen], [no],[$WITH_DOXYGEN] )
+else
+ AC_PATH_PROGS( [DOXYGEN], [doxygen], [no],[$PATH])
fi
-if test "x$ENABLE_DOXYGEN" != "xno"; then
+if test "x$ENABLE_DOCS" = "xyes"; then
if test "x$DOXYGEN" = "xno"; then
AC_MSG_ERROR("Doxygen executable not found.")
fi
fi
AC_SUBST([DOXYGEN])
-AM_CONDITIONAL([ENABLE_DOXYGEN],[test "x$ENABLE_DOXYGEN" != "xno"])
+AM_CONDITIONAL([ENABLE_DOCS],[test "x$ENABLE_DOCS" != "xno"])
################################################################################
@@ -554,8 +573,6 @@ AC_ARG_WITH([parmetis],
;;
esac] )
-
-
################################################################################
# Optional Libraries
################################################################################
@@ -620,7 +637,6 @@ AC_HELP_STRING([--disable-tools],[Disable all tools by default])],
[ENABLE_TOOLS=$enableval],[ENABLE_TOOLS=] )
# Individual tools
MB_OPTIONAL_TOOL([mbconvert], [yes])
-MB_OPTIONAL_TOOL([pointsearch], [no])
MB_OPTIONAL_TOOL([hexmodops], [yes])
MB_OPTIONAL_TOOL([qvdual], [no] )
MB_OPTIONAL_TOOL([vtkMOABReader],[${VTKMOAB_DEFAULT}] )
@@ -698,7 +714,7 @@ if test "xyes" = "x$WITH_ZOLTAN"; then
[#include <stdlib.h>])
CPPFLAGS="$old_CPPFLAGS"
AC_LANG_POP([C++])
-
+
# Build LDFLAGS for Zoltan
found=no
if test "x" != "x$ZOLTAN_DIR"; then
@@ -720,7 +736,7 @@ if test "xyes" = "x$WITH_ZOLTAN"; then
ZOLTAN_LIB_FLAGS="$ZOLTAN_LIB_FLAGS -L$ZOLTAN_DIR"
fi
fi
-
+
# Test for Zoltan library, possibly requiring parmetis
old_LDFLAGS="$LDFLAGS"
LDFLAGS="$LDFLAGS $ZOLTAN_LIB_FLAGS"
@@ -731,13 +747,17 @@ if test "xyes" = "x$WITH_ZOLTAN"; then
[AC_MSG_ERROR([zoltan library not found or not usable. Perhaps you need to specify --with-parmetis=@<:@DIR@:>@?])],
[-lparmetis -lmetis])])
LDFLAGS="$old_LDFLAGS"
-
+
+ # Define macro HAVE_ZOLTAN
+ AM_CPPFLAGS="$AM_CPPFLAGS -DHAVE_ZOLTAN"
fi
AC_SUBST(ZOLTAN_LIBS)
AC_SUBST(ZOLTAN_LIB_FLAGS)
AC_SUBST(ZOLTAN_INC_FLAGS)
+EXPORT_LDFLAGS="$ZOLTAN_LIB_FLAGS $EXPORT_LDFLAGS"
+
################################################################################
# VTK for qvdual and/or vtkMOABReader
################################################################################
@@ -1143,6 +1163,7 @@ AC_SUBST(CGM_CPPFLAGS)
AC_SUBST(CGM_LDFLAGS)
AC_SUBST(CGM_LTFLAGS)
AC_SUBST(CGM_LIBS)
+AC_SUBST(CGM_DIR)
AM_CONDITIONAL( HAVE_CGM, [test "x$CGM_MISSING" = "xno"] )
@@ -1180,6 +1201,7 @@ AC_CONFIG_FILES([Makefile
src/Makefile
src/io/Makefile
src/io/mhdf/Makefile
+ src/LocalDiscretization/Makefile
src/parallel/Makefile
src/oldinc/Makefile
test/Makefile
@@ -1187,6 +1209,7 @@ AC_CONFIG_FILES([Makefile
test/dual/Makefile
test/obb/Makefile
test/perf/Makefile
+ test/perf/point_location/Makefile
test/io/Makefile
test/parallel/Makefile
test/oldinc/Makefile
@@ -1209,7 +1232,6 @@ AC_CONFIG_FILES([Makefile
tools/vtkMOABReaderNew/CMakeLists.txt
doc/Makefile
doc/user.dox
- doc/dev.dox
doc/config.tex
MeshFiles/Makefile
MeshFiles/unittest/Makefile
diff --git a/doc/DG/figure1.jpg b/doc/DG/figure1.jpg
new file mode 100644
index 0000000..e102a97
Binary files /dev/null and b/doc/DG/figure1.jpg differ
diff --git a/doc/DG/figure2.jpg b/doc/DG/figure2.jpg
new file mode 100644
index 0000000..7e8b199
Binary files /dev/null and b/doc/DG/figure2.jpg differ
diff --git a/doc/DG/figure3.jpg b/doc/DG/figure3.jpg
new file mode 100644
index 0000000..8810bfc
Binary files /dev/null and b/doc/DG/figure3.jpg differ
diff --git a/doc/DG/moabDG.h b/doc/DG/moabDG.h
new file mode 100644
index 0000000..32d7fe1
--- /dev/null
+++ b/doc/DG/moabDG.h
@@ -0,0 +1,273 @@
+/*! \page developerguide Developer's Guide (MOAB 4.6)
+
+ \subpage dg-contents
+
+ \subpage dg-figures
+
+*/
+
+/*! \page dg-figures List of Figures
+
+ \ref figure1
+
+ \ref figure2
+
+ \ref figure3
+*/
+
+
+/*! \page dg-contents Table of Contents
+
+ \ref sequence
+
+ \ref manager
+
+ \ref s-mesh
+
+ \ref sets
+
+ \section sequence EntitySequence & SequenceData
+
+ \subsection figure1 Figure 1: EntitySequences For One SequenceData
+ \image html figure1.jpg
+
+ \ref dg-figures "List of Figures"
+
+The <I>SequenceData</I> class manages as set of arrays of per-entity values. Each
+<I>SequenceData</I> has a start and end handle denoting the block of entities for which
+the arrays contain data. The arrays managed by a <I>SequenceData</I> instance are
+divided into three groups:
+
+- Type-specific data (connectivity, coordinates, etc.): zero or more arrays.
+- Adjacency data: zero or one array.
+- Dense tag data: zero or more arrays.
+.
+
+The abstract <I>EntitySequence</I> class is a non-strict subset of a <I>SequenceData</I>.
+It contains a pointer to a <I>SequenceData</I> and the start and end handles to indi-
+cate the subset of the referenced <I>SequenceData</I>. The <I>EntitySequence</I> class is
+used to represent the regions of valid (or allocated) handles in a <I>SequenceData</I>.
+A <I>SequenceData</I> is expected to be referenced by one or more <I>EntitySequence</I>
+instances.
+
+Initial <I>EntitySequence</I> and <I>SequenceData</I> pairs are typically created in one
+of two configurations. When reading from a file, a <I>SequenceData</I> will be created
+to represent all of a single type of entity contained in a file. As all entries in the <I>SequenceData</I> correspond to valid handles (entities read from the file) a single
+<I>EntitySequence</I> instance corresponding to the entire <I>SequenceData</I> is initially
+created. The second configuration arises when allocating a single entity. If no
+entities have been allocated yet, a new <I>SequenceData</I> must be created to store
+the entity data. It is created with a constant size (e.g. 4k entities). The new
+<I>EntitySequence</I> corresponds to only the first entity in the <I>SequenceData</I>: the
+one allocated entity. As subsequent entities are allocated, the <I>EntitySequence</I>
+is extended to cover more of the corresponding <I>SequenceData</I>.
+
+Concrete subclasses of the <I>EntitySequence</I> class are responsible for rep-
+resenting specific types of entities using the array storage provided by the
+<I>SequenceData</I> class. They also handle allocating <I>SequenceData</I> instances with
+appropriate arrays for storing a particular type of entity. Each concrete subclass
+typically provides two constructors corresponding to the two initial allocation
+configurations described in the previous paragraph. <I>EntitySequence</I> imple-
+mentations also provide a split method, which is a type of factory method. It
+modifies the called sequence and creates a new sequence such that the range of
+entities represented by the original sequence is split.
+
+The <I>VertexSequence</I> class provides an <I>EntitySequence</I> for storing ver-
+tex data. It references a SequenceData containing three arrays of doubles
+for storing the blocked vertex coordinate data. The <I>ElementSequence</I> class
+extends the <I>EntitySequence</I> interface with element-specific functionality. The
+<I>UnstructuredElemSeq</I> class is the concrete implementation of <I>ElementSequence</I>
+used to represent unstructured elements, polygons, and polyhedra. <I>MeshSetSequence</I>
+is the <I>EntitySequence</I> used for storing entity sets.
+
+Each <I>EntitySequence</I> implementation also provides an implementation of
+the values per entity method. This value is used to determine if an exist-
+ing <I>SequenceData</I> that has available entities is suitable for storing a particular
+entity. For example, <I>UnstructuredElemSeq</I> returns the number of nodes per el-
+ement from values per entity. When allocating a new element with a specific
+number of nodes, this value is used to determine if that element may be stored
+in a specific <I>SequenceData</I>. For vertices, this value is always zero. This could
+be changed to the number of coordinates per vertex, allowing representation of
+mixed-dimension data. However, API changes would be required to utilize such
+a feature. Sequences for which the corresponding data cannot be used to store
+new entities (e.g. structured mesh discussed in a later section) will return -1 or
+some other invalid value.
+
+ \ref dg-contents "Top"
+
+ \section manager TypeSequenceManager & SequenceManager
+
+The <I>TypeSequenceManager</I> class maintains an organized set of <I>EntitySequence</I>
+instances and corresponding <I>SequenceData</I> instances. It is used to manage
+all such instances for entities of a single <I>EntityType</I>. <I>TypeSequenceManager</I>
+enforces the following four rules on its contained data:
+
+-# No two <I>SequenceData</I> instances may overlap.
+-# No two <I>EntitySequence</I> instances may overlap.
+-# Every <I>EntitySequence</I> must be a subset of a <I>SequenceData</I>.
+-# Any pair of <I>EntitySequence</I> instances referencing the same <I>SequenceData</I> must be separated by at least one unallocated handle.
+.
+
+ \subsection figure2 Figure 2: SequenceManager and Related Classes
+ \image html figure2.jpg
+
+ \ref dg-figures "List of Figures"
+
+The first three rules are required for the validity of the data model. The
+fourth rule avoids unnecessary inefficiency. It is implemented by merging such
+adjacent sequences. In some cases, other classes (e.g. <I>SequenceManager</I>) can
+modify an <I>EntitySequence</I> such that the fourth rule is violated. In such cases,
+the <I>TypeSequenceManager::notify</I> prepended or <I>TypeSequenceManager::notify</I> appended
+method must be called to maintain the integrity of the data<sup>1</sup>. The above rules
+(including the fourth) are assumed in many other methods of the <I>TypeSequenceManager</I>
+class, such that those methods will fail or behave unexpectedly if the managed
+data does not conform to the rules.
+
+<I>TypeSequenceManager</I> contains three principal data structures. The first is
+a <I>std::set</I> of <I>EntitySequence</I> pointers sorted using a custom comparison op-
+erator that queries the start and end handles of the referenced sequences. The
+comparison operation is defined as: <I>a->end_handle() < b->start_handle()</I>.
+This method of comparison has the advantage that a sequence corresponding to
+a specific handle can be located by searching the set for a “sequence” beginning
+and ending with the search value. The lower bound and find methods pro-
+vided by the library are guaranteed to return the sequence, if it exists. Using
+such a comparison operator will result in undefined behavior if the set contains
+overlapping sequences. This is acceptable, as rule two above prohibits such
+a configuration. However, some care must be taken in writing and modifying
+methods in <I>TypeSequenceManager</I> so as to avoid having overlapping sequences
+as a transitory state of some operation.
+
+The second important data member of <I>TypeSequenceManager</I> is a pointer
+to the last referenced <I>EntitySequence</I>. This “cached” value is used to speed up
+searches by entity handle. This pointer is never null unless the sequence is empty.
+This rule is maintained to avoid unnecessary branches in fast query paths. In
+cases where the last referenced sequence is deleted, <I>TypeSequenceManager</I> will
+typically assign an arbitrary sequence (e.g. the first one) to the last referenced
+pointer.
+
+The third data member of <I>TypeSequenceManager</I> is a <I>std::set</I> of <I>SequenceData</I>
+instances that are not completely covered by a <I>EntitySequence</I> instance<sup>2</sup>.
+This list is searched when allocating new handles. <I>TypeSequenceManager</I> also
+embeds in each <I>SequenceData</I> instance a reference to the first corresponding
+<I>EntitySequence</I> so that it may be located quickly from only the <I>SequenceData</I>
+pointer.
+
+The <I>SequenceManager</I> class contains an array of <I>TypeSequenceManager</I> in-
+stances, one for each <I>EntityType</I>. It also provides all type-specific operations
+such as allocating the correct <I>EntitySequence</I> subtype for a given <I>EntityType</I>.
+
+<sup>1</sup>This source of potential error can be eliminated with changes to the entity set representation. This is discussed in a later section.
+
+<sup>2</sup>Given rule four for the data managed by a <I>TypeSequenceManager</I>, any <I>SequenceData</I> for which all handles are allocated will be referenced by exactly one <I>EntitySequence</I>.
+
+ \ref dg-contents "Top"
+
+ \section s-mesh Structured Mesh
+
+Structured mesh storage is implemented using subclasses of <I>SequenceData</I>:
+<I>ScdElementData</I> and <I>ScdVertexData</I>. The <I>StructuredElementSeq</I> class is
+used to access the structured element connectivity. A standard <I>VertexSequence</I>
+instance is used to access the ScdVertexData because the vertex data storage
+is the same as for unstructured mesh.
+
+ \ref dg-contents "Top"
+
+ \section sets Entity Sets
+
+- MeshSetSequence
+
+The <I>MeshSetSequence</I> class is the same as most other subclasses of <I>EntitySequence</I>
+in that it utilizes SequenceData to store its data. A single array in the <I>SequenceData</I>
+is used to store instances of the MeshSet class, one per allocated <I>EntityHandle</I>.
+<I>SequenceData</I> allocates all of its managed arrays using malloc and free as
+simple arrays of bytes. <I>MeshSetSequence</I> does in-place construction and de-
+struction of <I>MeshSet</I> instances within that array. This is similar to what is
+done by <I>std::vector</I> and other container classes that may own more storage
+than is required at a given time for contained objects.
+
+- MeshSet
+
+ \subsection figure3 Figure 3: SequenceManager and Related Classes
+ \image html figure3.jpg
+
+ \ref dg-figures "List of Figures"
+
+The <I>MeshSet</I> class is used to represent a single entity set instance in MOAB.
+The class is optimized to minimize storage (further possible improvements in
+storage size are discussed later.)
+
+Figure 3 shows the memory layout of an instance of the <I>MeshSet</I> class.
+The flags member holds the set creation bit flags: <I>MESHSET_TRACK_OWNER</I>,
+<I>MESHSET_SET</I>, and <I>MESHSET_ORDERED</I>. The presence of the <I>MESHSET_TRACK_OWNER</I>
+indicates that reverse links from the contained entities back to the owning set
+should be maintained in the adjacency list of each entity. The <I>MESHSET_SET</I>
+and <I>MESHSET_ORDERED</I> bits are mutually exclusive, and as such most code only
+tests for the <I>MESHSET_ORDERED</I>, meaning that in practice the <I>MESHSET_SET</I> bit is
+ignored. <I>MESHSET_ORDERED</I> indicates that the set may contain duplicate handles
+and that the order that the handles are added to the set should be preserved.
+In practice, such sets are stored as a simple list of handles. <I>MESHSET_SET</I> (or in
+practice, the lack of <I>MESHSET_ORDERED</I>) indicates that the order of the handles
+need not be preserved and that the set may not contain duplicate handles. Such
+sets are stored in a sorted range-compacted format similar to that of the Range
+class.
+
+The memory for storing contents, parents, and children are each handled in
+the same way. The data in the class is composed of a 2-bit ‘size’ field and two
+values, where the two values may either be two handles or two pointers. The size
+bit-fields are grouped together to reduce the required amount of memory. If the
+numerical value of the 2-bit size field is 0 then the corresponding list is empty.
+If the 2-bit size field is either 1 or 2, then the contents of the corresponding list
+are stored directly in the corresponding two data fields of the MeshSet object.
+If the 2-bit size field has a value of 3 (11 binary), then the corresponding two
+data fields store the begin and end pointers of an external array of handles.
+The number of handles in the external array can be obtained by taking the
+difference of the start and end pointers. Note that unlike <I>std::vector</I>, we
+do not store both an allocated and used size. We store only the ‘used’ size
+and call std::realloc whenever the used size is modified, thus we rely on the
+std::malloc implementation in the standard C library to track ‘allocated’ size
+for us. In practice this performs well but does not return memory to the ‘system’
+when lists shrink (unless they shrink to zero). This overall scheme could exhibit
+poor performance if the size of one of the data lists in the set frequently changes
+between less than two and more than two handles, as this will result in frequent
+releasing and re-allocating of the memory for the corresponding array.
+
+If the <I>MESHSET_ORDERED</I> flag is not present, then the set contents list (parent
+and child lists are unaffected) is stored in a range-compacted format. In this
+format the number of handles stored in the array is always a multiple of two.
+Each consecutive pair of handles indicate the start and end, inclusive, of a range
+of handles contained in the set. All such handle range pairs are stored in sorted
+order and do not overlap. Nor is the end handle of one range ever one less than
+the start handle of the next. All such ‘adjacent’ range pairs are merged into a
+single pair. The code for insertion and removal of handles from range-formatted
+set content lists is fairly complex. The implementation will guarantee that a
+given call to insert entities into a range or remove entities from a range is never
+worse than O(ln n) + O(m + n), where ‘n’ is the number of handles to insert
+and ‘m’ is the number of handles already contained in the set. So it is generally
+much more efficient to build Ranges of handles to insert (and remove) and call
+MOAB to insert (or remove) the entire list at once rather than making may
+calls to insert (or remove) one or a few handles from the contents of a set.
+The set storage could probably be further minimized by allowing up to six
+handles in one of the lists to be elided. That is, as there are six potential ‘slots’
+in the MeshSet object then if two of the lists are empty it should be possible to store up to six values of the remaining list directly in the MeshSet object.
+However, the additional runtime cost of such complexity could easily outweigh
+any storage advantage. Further investigation into this has not been done because
+the primary motivation for the storage optimization was to support binary trees.
+
+Another possible optimization of storage would be to remove the <I>MeshSet</I>
+object entirely and instead store the data in a ‘blocked’ format. The corre-
+sponding <I>SequenceData</I> would contain four arrays: flags, parents, children, and
+contents instead of a single array of <I>MeshSet</I> objects. If this were done then
+no storage need ever be allocated for parent or child links if none of the sets
+in a <I>SequenceData</I> has parent or child links. The effectiveness of the storage
+reduction would depend greatly on how sets get grouped into <I>SequenceDatas</I>.
+This alternate storage scheme might also allow for better cache utilization as it
+would group like data together. It is often the case that application code that
+is querying the contents of one set will query the contents of many but never
+query the parents or children of any set. Or that an application will query only
+parent or child links of a set without every querying other set properties. The
+downside of this solution is that it makes the implementation a little less mod-
+ular and maintainable because the existing logic contained in the <I>MeshSet</I> class
+would need to be spread throughout the <I>MeshSetSequence</I> class.
+
+ \ref dg-contents "Top"
+*/
+
diff --git a/doc/DG/styleguide.h b/doc/DG/styleguide.h
new file mode 100644
index 0000000..3a13e45
--- /dev/null
+++ b/doc/DG/styleguide.h
@@ -0,0 +1,91 @@
+/*!
+ \page styleguide Coding Style Guide
+Code developed in %MOAB should follow the coding styles described here. Any deviations from this style
+guide will result in severe berating and other verbal abuse.
+
+\section dirstructure MOAB Directory Structure
+%MOAB source code is organized in the following directory structure: \n
+ - doc: Documentation is put here, along with the input file for Doxygen. Most %MOAB documentation is doxygen-processed.
+ - examples: Examples of %MOAB usage, both small and large. These programs are not meant to be used as unit tests, but
+ rather as further documentation on %MOAB usage.
+ - src : Mesh related source codes. It includes:
+ - io: %MOAB Input/Output classes.
+ - moab: %MOAB core classes.
+ - lotte: Computational Meshing basics.
+ - parallel: Parallel mesh computation, i/o data processing methods.
+ - test: All unit test programs should go below this directory.
+ Please put the unit tests into their related subdirectories based on the test's
+ purpose if possible.
+If you're designing a new class or other code for %MOAB and are not sure where to put it, try to find something similar
+and put it there. Otherwise, email the %MOAB email list for pointers. <em> In general, you should not need to create new
+subdirectories in the %MOAB source code, except when implementing a new algorithm with more than about 2 files.</em>
+
+\section sourcestyle Source Code Style and Best Practices
+%MOAB code should abide by the following general rules:
+ - Names:
+ - Class names should be in the CamelBack style, e.g. EdgeMesh or VertexMesher.
+ - Class member variables should be camelBack, e.g. EdgeMesh::schemeType; each member variable, e.g. int memberVariable,
+ should have set/get functions void member_variable(int newval) and int member_variable(), respectively.
+ - Enumeration values should be all captitalized, with underscores avoided if possible (the enumeration name indicates
+ the general purpose of the enumeration, so e.g. we use EQUAL, not EQUAL_MESH)
+ - Source code should not contain tabs or MS-DOS newlines; tabs and other indentations should be set to a width of 2 spaces.
+ For general tips on how to set your editor for this, see the %MOAB-dev discussion starting with <a href="https://lists.mcs.anl.gov/mailman/private/moab-dev/2011/000519.html">this message</a>.
+ - Each class header should be fully commented; that includes:
+ - A \\file comment block at the top of the file; DO NOT include things like Author and Date blocks; this stuff is available
+ from subversion if we really need to know.
+ - A \\class comment block, formatted like those in the %MOAB core classes. THE FIELDS AFTER THE CLASS NAME ARE VERY IMPORTANT,
+ as they tell developers how to include the class declaration in their code. This information goes into the "Detailed
+ Description" portion of the class documentation. This block should include any features or limitations of the class.
+ Eventually, we'll impose some standard boilerplate that each meshing class should use.
+ Until then, please keep this block to around a paragraph.
+ - Each function in both the public and private interfaces should be commented, INCLUDING ANY ARGUMENTS AND RETURN VALUES.
+ See the %MOAB classes for examples of how to format these comments. As a rule of thumb, your code should run through
+ Doxygen without generating any warnings; in fact, Doxygen is sometimes helpful at pointing out inconsistencies in your
+ class declaration.
+ - Developers should avoid using \#include in header files, as they propagate dependencies more widely than necessary. The only
+ cases where other includes are needed are to import the declaration for a parent class, and to declare types used as
+ non-pointer and non-reference function arguments. In most cases, a forward-declaration statement (e.g. 'class MKCore')
+ will suffice.
+ - Naming classes and other top-level constructs:
+ - No names should be added to the global namespace. Everything should be
+ in the MOAB namespace. An exception can be made for names with a static
+ scope declared in a .cpp file, but class member functions never have a
+ static scope.
+ - Names should be kept as private as possible. If declaring a struct or
+ utility class that is used internally by some other class, consider
+ defining it in the .cpp file of the main class or a separate header
+ only included in that .cpp file and using (if necessary) only forward
+ delcarations (e.g. \c struct \c Point3D;) in the header file used
+ by other code. If that is not possible, then consider nesting the
+ definitions such that the scope of the name is limited to that of the
+ class using it.
+ - Any names introduced into the top-level MOAB namespace should be
+ sufficiently unique to avoid conflicts with other code. If you must
+ introduce a class to the top-level MOAB namespace, don't choose
+ an overly genereric name like \c Point3D .
+ - Constants and Macros
+ - Don't use a pre-processor macro where a const variable or an inline or
+ template function will suffice.
+ There is absolutely benefit to the former over the later with modern
+ compilers. Further, using macros bypasses typechecking that the compiler
+ would otherwise do for you and if used in headers, introduce names into
+ the global rather than MOAB namespace.
+ - Don't define constants that are already provided by standard libraries.
+ For example, use \c M_PI as defined in \c math.h rather than defining
+ your own constant.
+\section commits Making Repository Commits
+As a general rule, developers should update frequently, and commit changes often. However, the repository should always remain
+in a state where the code can be compiled. Most of the time, the code should also successfully execute "make check" run from the
+top-level directory. If you commit code that violates this principal, it should be your first priority to return the repository
+code to a compilable state, and your second priority to make sure "make check" runs without errors.
+
+Commits to the repository should also come with a non-trivial, useful, non-verbose log message. Oftentimes the best way to generate
+this message is to run 'commit -a', and include a comment on
+each file that changed, then Ctrl+O to write out, followed by 'Enter' and Ctrl+X. Many times it is helpful to state that 'make check runs successfully' at the end of the log message.
+Although it would be possible and many software projects do it, we prefer not to force successful execution of the test suite
+before every commit. Developers should make every effort to avoid having to impose this constraint, by running a make check
+before every commit.
+
+Top: \ref index
+
+ */
diff --git a/doc/MOAB-UG.doc b/doc/MOAB-UG.doc
index 281d619..b7668b8 100644
Binary files a/doc/MOAB-UG.doc and b/doc/MOAB-UG.doc differ
This diff is so big that we needed to truncate the remainder.
https://bitbucket.org/fathomteam/moab/commits/d67a959548a4/
Changeset: d67a959548a4
Branch: None
User: pshriwise
Date: 2014-01-07 16:52:20
Summary: Merged fathomteam/moab into master
Affected #: 11 files
diff --git a/tools/mbcslam/CslamUtils.cpp b/tools/mbcslam/CslamUtils.cpp
index 6cca02d..e806d61 100644
--- a/tools/mbcslam/CslamUtils.cpp
+++ b/tools/mbcslam/CslamUtils.cpp
@@ -30,7 +30,7 @@ double area2D(double *a, double *b, double *c)
// (b-a)x(c-a) / 2
return ((b[0] - a[0]) * (c[1] - a[1]) - (b[1] - a[1]) * (c[0] - a[0])) / 2;
}
-int borderPointsOfXinY2(double * X, int nX, double * Y, int nY, double * P, int side[MAXEDGES])
+int borderPointsOfXinY2(double * X, int nX, double * Y, int nY, double * P, int side[MAXEDGES], double epsilon_area)
{
// 2 triangles, 3 corners, is the corner of X in Y?
// Y must have a positive area
@@ -54,7 +54,7 @@ int borderPointsOfXinY2(double * X, int nX, double * Y, int nY, double * P, int
double area2 = (B[0] - A[0]) * (C[1] - A[1])
- (C[0] - A[0]) * (B[1] - A[1]);
- if (area2 < 0.)
+ if (area2 < -epsilon_area)
{
inside = 0;
break;
diff --git a/tools/mbcslam/CslamUtils.hpp b/tools/mbcslam/CslamUtils.hpp
index 05775df..66f4e99 100644
--- a/tools/mbcslam/CslamUtils.hpp
+++ b/tools/mbcslam/CslamUtils.hpp
@@ -21,7 +21,7 @@ namespace moab
{
double dist2(double * a, double * b);
double area2D(double *a, double *b, double *c);
-int borderPointsOfXinY2(double * X, int nX, double * Y, int nY, double * P, int side[MAXEDGES]);
+int borderPointsOfXinY2(double * X, int nX, double * Y, int nY, double * P, int side[MAXEDGES], double epsilon_area);
int SortAndRemoveDoubles2(double * P, int & nP, double epsilon);
// the marks will show what edges of blue intersect the red
diff --git a/tools/mbcslam/Intx2Mesh.cpp b/tools/mbcslam/Intx2Mesh.cpp
index c1da3fd..625a442 100644
--- a/tools/mbcslam/Intx2Mesh.cpp
+++ b/tools/mbcslam/Intx2Mesh.cpp
@@ -18,7 +18,7 @@
namespace moab {
-Intx2Mesh::Intx2Mesh(Interface * mbimpl):mb(mbimpl), parcomm(NULL), myTree(NULL), remote_cells(NULL)
+Intx2Mesh::Intx2Mesh(Interface * mbimpl): mb(mbimpl), parcomm(NULL), remote_cells(NULL)
{
dbg_1=0;
box_error=0;
diff --git a/tools/mbcslam/Intx2Mesh.hpp b/tools/mbcslam/Intx2Mesh.hpp
index 2daa052..d54b0b5 100644
--- a/tools/mbcslam/Intx2Mesh.hpp
+++ b/tools/mbcslam/Intx2Mesh.hpp
@@ -34,7 +34,6 @@ namespace moab {
// forward declarations
class ParallelComm;
-class AdaptiveKDTree;
class TupleList;
class Intx2Mesh
@@ -68,7 +67,7 @@ public:
ErrorCode GetOrderedNeighbors(EntityHandle set, EntityHandle quad,
EntityHandle neighbors[MAXEDGES]);
- void SetErrorTolerance(double eps) { epsilon_1=eps;}
+ void SetErrorTolerance(double eps) { epsilon_1=eps; epsilon_area = eps*eps/2;}
//void SetEntityType (EntityType tp) { type=tp;}
@@ -97,7 +96,8 @@ public:
ErrorCode correct_intersection_points_positions();
- void enable_debug() {dbg_1=1;};
+ void enable_debug() {dbg_1 = 1;}
+ void disable_debug() {dbg_1 = 0;}
protected: // so it can be accessed in derived classes, InPlane and OnSphere
Interface * mb;
@@ -139,10 +139,10 @@ protected: // so it can be accessed in derived classes, InPlane and OnSphere
std::vector<std::vector<EntityHandle> *> extraNodesVec;
double epsilon_1;
+ double epsilon_area;
ParallelComm * parcomm;
- AdaptiveKDTree *myTree;
std::vector<double> allBoxes;
double box_error;
/* \brief Local root of the kdtree */
diff --git a/tools/mbcslam/Intx2MeshInPlane.cpp b/tools/mbcslam/Intx2MeshInPlane.cpp
index 9fb2a94..2ddea3a 100644
--- a/tools/mbcslam/Intx2MeshInPlane.cpp
+++ b/tools/mbcslam/Intx2MeshInPlane.cpp
@@ -106,7 +106,7 @@ int Intx2MeshInPlane::computeIntersectionBetweenRedAndBlue(EntityHandle red, Ent
}
int side[MAXEDGES] = { 0 };// this refers to what side? blue or red?
- int extraPoints = borderPointsOfXinY2(blueCoords2D, nsBlue, redCoords2D, nsRed, &(P[2 * nP]), side);
+ int extraPoints = borderPointsOfXinY2(blueCoords2D, nsBlue, redCoords2D, nsRed, &(P[2 * nP]), side, epsilon_area);
if (extraPoints >= 1)
{
for (int k = 0; k < nsBlue; k++)
@@ -133,7 +133,7 @@ int Intx2MeshInPlane::computeIntersectionBetweenRedAndBlue(EntityHandle red, Ent
}
nP += extraPoints;
- extraPoints = borderPointsOfXinY2(redCoords2D, nsRed, blueCoords2D, nsBlue, &(P[2 * nP]), side);
+ extraPoints = borderPointsOfXinY2(redCoords2D, nsRed, blueCoords2D, nsBlue, &(P[2 * nP]), side, epsilon_area);
if (extraPoints >= 1)
{
for (int k = 0; k < nsRed; k++)
diff --git a/tools/mbcslam/Intx2MeshOnSphere.cpp b/tools/mbcslam/Intx2MeshOnSphere.cpp
index a1d786b..ef4d219 100644
--- a/tools/mbcslam/Intx2MeshOnSphere.cpp
+++ b/tools/mbcslam/Intx2MeshOnSphere.cpp
@@ -124,7 +124,7 @@ int Intx2MeshOnSphere::computeIntersectionBetweenRedAndBlue(EntityHandle red, En
return 1; // some unforeseen error
int side[MAXEDGES] = { 0 };// this refers to what side? blue or red?
- int extraPoints = borderPointsOfXinY2(blueCoords2D, nsBlue, redCoords2D, nsRed, &(P[2 * nP]), side);
+ int extraPoints = borderPointsOfXinY2(blueCoords2D, nsBlue, redCoords2D, nsRed, &(P[2 * nP]), side, epsilon_area);
if (extraPoints >= 1)
{
for (int k = 0; k < nsBlue; k++)
@@ -144,7 +144,7 @@ int Intx2MeshOnSphere::computeIntersectionBetweenRedAndBlue(EntityHandle red, En
}
nP += extraPoints;
- extraPoints = borderPointsOfXinY2(redCoords2D, nsRed, blueCoords2D, nsBlue, &(P[2 * nP]), side);
+ extraPoints = borderPointsOfXinY2(redCoords2D, nsRed, blueCoords2D, nsBlue, &(P[2 * nP]), side, epsilon_area);
if (extraPoints >= 1)
{
for (int k = 0; k < nsRed; k++)
@@ -370,10 +370,11 @@ int Intx2MeshOnSphere::findNodes(EntityHandle red, int nsRed, EntityHandle blue,
std::stringstream fff;
fff << "file0" << count<< ".vtk";
- mb->write_mesh(fff.str().c_str(), &outSet, 1);
+ mb->write_mesh(fff.str().c_str(), &outSet, 1);
}
}
+ disable_debug();
delete[] foundIds;
foundIds = NULL;
return 0;
diff --git a/tools/mbcslam/README b/tools/mbcslam/README
new file mode 100644
index 0000000..3501d1c
--- /dev/null
+++ b/tools/mbcslam/README
@@ -0,0 +1,26 @@
+module for intersection + transport methods in moab
+
+contains several utilities and "make check" tests
+the core is the intersection class, Intx2Mesh, with 2 derived
+classes, Intx2MeshInPlane and Intx2MeshOnSphere
+it takes 2 mesh sets that overlap in plane or on a sphere, and
+creates an output set with convex polygons. Each convex polygon
+in the output set has 2 parents, each parent in one of the input sets
+"convex" is understood either on a gnomonic plane of a
+sphere or in a plane
+
+the input files can be visualized in visit or paraview, and for that we prefer
+conversion to vtk, if needed.
+True writing/reading works in parallel only with hdf5 format
+
+example for advection:
+create_dp ~/source/MOABsource/MeshFiles/unittest/mpas_p8.h5m HN16DP.h5m -t 0.1 -dt 0.05
+
+then
+intx_imesh
+
+mbconvert out.h5m out.vtk
+
+compare the TracerAverage plot on HN16DP.h5m with the plot on out.vtk
+
+
diff --git a/tools/mbcslam/advection.F90 b/tools/mbcslam/advection.F90
index fb102ba..ad92f88 100644
--- a/tools/mbcslam/advection.F90
+++ b/tools/mbcslam/advection.F90
@@ -26,7 +26,7 @@ program advection
use ISO_C_BINDING
implicit none
iMesh_Instance, INTENT(IN) , VALUE :: instance
- iBase_EntityHandle, INTENT(OUT) :: opEulerSet
+ iBase_EntityHandle, INTENT(IN), VALUE :: opEulerSet
integer(c_int) , INTENT (OUT) :: ierr
END SUBROUTINE update_tracer
END INTERFACE
@@ -35,18 +35,22 @@ program advection
! imesh is the instance handle
iMesh_Instance imesh
- ! ents, verts will be arrays storing vertex/entity handles
- iBase_EntityHandle, pointer :: ents, verts
+ ! cells will be storing 2d cells
+ TYPE(C_PTR) cells_ptr
+ iBase_EntityHandle, pointer :: cells(:)
+ INTEGER ents_alloc, ents_size
+
iBase_EntitySetHandle root_set
iBase_EntitySetHandle opEulerSet
CHARACTER (LEN=200) options
CHARACTER (LEN=200) filename
CHARACTER (LEN=200) optionswrite
CHARACTER (LEN=200) outname
- TYPE(C_PTR) :: vertsPtr, entsPtr
+ ! TYPE(C_PTR) :: vertsPtr, entsPtr
integer rank, sz, ierr
integer lenopt, lenname
+ integer isList
#ifdef USE_MPI
! local variables for parallel runs
@@ -95,6 +99,24 @@ program advection
%VAL(lenopt) ) !119) !strlen(options));
CHECK("fail to load mesh in parallel ")
+ isList = 0
+ call iMesh_createEntSet(%VAL(imesh), %VAL(isList), opEulerSet, ierr)
+ CHECK("fail to create euler set ")
+
+ ents_alloc = 0
+ ents_size = 0
+
+ call iMesh_getEntities(%VAL(imesh),%VAL(root_set), &
+ %VAL(iBase_FACE), &
+ %VAL(iMesh_ALL_TOPOLOGIES), cells_ptr, &
+ ents_alloc, ents_size, ierr);
+ CHECK("fail to get 2d entities ")
+
+ call c_f_pointer(cells_ptr, cells, [ents_size])
+
+ call iMesh_addEntArrToSet(%VAL(imesh), cells, %VAL(ents_size), &
+ %VAL(opEulerSet), ierr)
+
call update_tracer(imesh, opEulerSet, ierr)
CHECK("fail to update tracer ")
diff --git a/tools/mbcslam/intx_imesh.cpp b/tools/mbcslam/intx_imesh.cpp
index 0c49ac4..cb7b494 100644
--- a/tools/mbcslam/intx_imesh.cpp
+++ b/tools/mbcslam/intx_imesh.cpp
@@ -16,7 +16,7 @@
#ifdef __cplusplus
extern "C" {
#endif
- void update_tracer( iMesh_Instance instance, iBase_EntitySetHandle * opEulerSet, int * ierr);
+ void update_tracer( iMesh_Instance instance, iBase_EntitySetHandle eulerSet, int * ierr);
#ifdef __cplusplus
} // extern "C"
@@ -71,7 +71,23 @@ int main(int argc, char* argv[]){
printf("There's %d entity sets here on process rank %d \n", num_sets, rank);
iBase_EntitySetHandle euler_set;
- update_tracer( imesh, &euler_set, &ierr);
+
+ iMesh_createEntSet(imesh, 0, &euler_set, &ierr);
+ IMESH_ASSERT(ierr);
+
+ iBase_EntityHandle *cells = NULL;
+ int ents_alloc = 0;
+ int ents_size = 0;
+
+ iMesh_getEntities(imesh, root, iBase_FACE, iMesh_ALL_TOPOLOGIES, &cells,
+ &ents_alloc, &ents_size, &ierr);
+ IMESH_ASSERT(ierr);
+
+ iMesh_addEntArrToSet(imesh, cells, ents_size, euler_set, &ierr);
+
+ IMESH_ASSERT(ierr);
+
+ update_tracer( imesh, euler_set, &ierr);
IMESH_ASSERT(ierr);
// write everything
diff --git a/tools/mbcslam/intx_on_sphere_test.cpp b/tools/mbcslam/intx_on_sphere_test.cpp
index 350a39b..360aa28 100644
--- a/tools/mbcslam/intx_on_sphere_test.cpp
+++ b/tools/mbcslam/intx_on_sphere_test.cpp
@@ -69,10 +69,12 @@ int main(int argc, char* argv[])
return 1;
Intx2MeshOnSphere worker(mb);
- double radius= 6. * sqrt(3.) / 2; // input
- worker.SetErrorTolerance(radius*1.e-8);
+
+ worker.SetErrorTolerance(R*1.e-8);
//worker.SetEntityType(moab::MBQUAD);
- worker.SetRadius(radius);
+ worker.SetRadius(R);
+ //worker.enable_debug();
+
rval = worker.intersect_meshes(sf1, sf2, outputSet);
//compute total area with 2 methods
@@ -81,11 +83,15 @@ int main(int argc, char* argv[])
rval = mb->write_mesh(newFile, &outputSet, 1);
if (MB_SUCCESS != rval)
return 1;
- double area_method1 = area_on_sphere_lHuiller(mb, outputSet, radius);
- double area_method2 = area_on_sphere(mb, outputSet, radius);
+ double initial_area = area_on_sphere_lHuiller(mb, sf1, R);
+ double area_method1 = area_on_sphere_lHuiller(mb, outputSet, R);
+ double area_method2 = area_on_sphere(mb, outputSet, R);
+ std::cout << "initial area: " << initial_area << "\n";
std::cout<< " area with l'Huiller: " << area_method1 << " with Girard: " << area_method2<< "\n";
- std::cout << " relative difference: " << fabs(area_method1-area_method2)/area_method1 << "\n";
+ std::cout << " relative difference areas " << fabs(area_method1-area_method2)/area_method1 << "\n";
+ std::cout << " relative error " << fabs(area_method1-initial_area)/area_method1 << "\n";
+
return 0;
}
diff --git a/tools/mbcslam/wrap_intx.cpp b/tools/mbcslam/wrap_intx.cpp
index e38c897..b973acf 100644
--- a/tools/mbcslam/wrap_intx.cpp
+++ b/tools/mbcslam/wrap_intx.cpp
@@ -23,25 +23,14 @@ bool debug = false;
extern "C" {
#endif
-void update_tracer( iMesh_Instance instance, iBase_EntitySetHandle * opEulerSet, int * ierr)
+void update_tracer( iMesh_Instance instance, iBase_EntitySetHandle imesh_euler_set, int * ierr)
{
Range ents;
moab::Interface * mb =MOABI;
*ierr =1;
- ErrorCode rval = mb->get_entities_by_dimension(0, 2, ents);// root 0
- ERRORV(rval, "can't get all 2d entities from root");
- EntityHandle euler_set;
- //EntityHandle lagr_set;
-
- rval = mb->create_meshset(MESHSET_SET, euler_set);
- ERRORV(rval , "can't create arrival mesh set");
-
- *opEulerSet = (iBase_EntitySetHandle)euler_set;
-
- rval = mb->add_entities(euler_set, ents);
- ERRORV(rval , "can't add ents to arrival set");
+ EntityHandle euler_set = (EntityHandle) imesh_euler_set;
Intx2MeshOnSphere worker(mb);
worker.SetRadius(radius);
@@ -49,7 +38,8 @@ void update_tracer( iMesh_Instance instance, iBase_EntitySetHandle * opEulerSet,
worker.SetErrorTolerance(gtol);
EntityHandle covering_lagr_set;
- rval = mb->create_meshset(MESHSET_SET, covering_lagr_set);
+
+ ErrorCode rval = mb->create_meshset(MESHSET_SET, covering_lagr_set);
ERRORV(rval , "can't create covering set ");
// we need to update the correlation tag and remote tuples
https://bitbucket.org/fathomteam/moab/commits/eaab3e1d2b0e/
Changeset: eaab3e1d2b0e
Branch: None
User: pshriwise
Date: 2014-01-07 21:15:39
Summary: Merged fathomteam/moab into master
Affected #: 47 files
diff --git a/.gitignore b/.gitignore
index e824a6b..d842c2b 100644
--- a/.gitignore
+++ b/.gitignore
@@ -165,6 +165,7 @@ test/io/vtk_test
test/kd_tree_test
test/kd_tree_time
test/kd_tree_tool
+test/lloyd_smoother_test
test/mbcn_test
test/mbfacet_test
test/mbground_test
diff --git a/examples/DeformMeshRemap.cpp b/examples/DeformMeshRemap.cpp
new file mode 100644
index 0000000..5442d1d
--- /dev/null
+++ b/examples/DeformMeshRemap.cpp
@@ -0,0 +1,488 @@
+/** @example DeformMeshRemap.cpp
+ * Description: Account for mesh deformation of a solid due to structural mechanics\n
+ * In this example there are two meshes, a "master" and "slave" mesh. In the master mesh,
+ * the solid material is deformed, to mimic what happens when a solid heats up and deforms.
+ * The fluid mesh is smoothed to account for those deformations, and tags on the fluid are
+ * remapped to those new positions. Then mesh positions and state variables are transferred
+ * to the slave mesh, mimicing another mesh used by some other physics.
+ *
+ * To run: ./DeformMeshRemap [<master_meshfile><slave_meshfile>]\n
+ * (default values can run if users don't specify the mesh files)
+ */
+
+#include "moab/Core.hpp"
+#include "moab/Range.hpp"
+#include "moab/LloydSmoother.hpp"
+#include "moab/ProgOptions.hpp"
+#include "moab/BoundBox.hpp"
+#include "moab/SpatialLocator.hpp"
+#include "MBTagConventions.hpp"
+#include "DataCoupler.hpp"
+
+#ifdef USE_MPI
+# include "moab/ParallelComm.hpp"
+#endif
+
+#include <iostream>
+#include <set>
+#include <sstream>
+#include <assert.h>
+
+using namespace moab;
+using namespace std;
+
+#ifndef MESH_DIR
+#define MESH_DIR "."
+#endif
+
+ErrorCode read_file(string &fname, EntityHandle &seth,
+ Range &solids, Range &solid_elems, Range &fluids, Range &fluid_elems);
+void deform_func(const BoundBox &bbox, double *xold, double *xnew);
+ErrorCode deform_master(Range &fluid_elems, Range &solid_elems, Tag &xnew);
+ErrorCode smooth_master(int dim, Tag xnew, EntityHandle &master, Range &fluids);
+ErrorCode write_to_coords(Range &elems, Tag tagh);
+
+const int SOLID_SETNO = 100, FLUID_SETNO = 200;
+
+Interface *mb;
+#define RR(a) if (MB_SUCCESS != rval) {cout << a << endl; return MB_FAILURE;}
+
+const bool debug = true;
+
+class DeformMeshRemap
+{
+public:
+
+ //! enumerator for solid/fluid, master/slave
+ enum {MASTER=0, SLAVE, SOLID, FLUID};
+
+ //! constructor
+ //! if master is NULL, the MOAB part is run in serial;
+ //! if slave is NULL but the master isn't, the slave is copied from the master
+ //! Create communicators using moab::ParallelComm::get_pcomm
+ DeformMeshRemap(Interface *impl, ParallelComm *master = NULL, ParallelComm *slave = NULL);
+
+ //! destructor
+ ~DeformMeshRemap();
+
+ //! execute the deformed mesh process
+ ErrorCode execute();
+
+ //! add a set number
+ ErrorCode add_set_no(int fluid_or_solid, int set_no);
+
+ //! remove a set number
+ ErrorCode remove_set_no(int fluid_or_solid, int set_no);
+
+ //! get the set numbers
+ ErrorCode get_set_nos(int fluid_or_solid, std::set<int> &set_nos) const;
+
+ //! get the xNew tag handle
+ inline Tag x_new() const {return xNew;}
+
+ //! get the tag name
+ std::string x_new_name() const {return xNewName;}
+
+ //! set the tag name
+ void x_new_name(const std::string &name) {xNewName = name;}
+
+ //! get/set the file name
+ std::string get_file_name(int m_or_s) const;
+
+ //! get/set the file name
+ void set_file_name(int m_or_s, const std::string &name);
+
+private:
+ //! apply a known deformation to the solid elements, putting the results in the xNew tag; also
+ //! write current coordinates to the xNew tag for fluid elements
+ ErrorCode deform_master(Range &fluid_elems, Range &solid_elems, const char *tag_name = NULL);
+
+ //! read a file and establish proper ranges
+ ErrorCode read_file(int m_or_s, string &fname, EntityHandle &seth);
+
+ //! write the input tag to the coordinates for the vertices in the input elems
+ ErrorCode write_to_coords(Range &elems, Tag tagh);
+
+ //! write the tag to the vertices, then save to the specified file
+ ErrorCode write_and_save(Range &ents, EntityHandle seth, Tag tagh, const char *filename);
+
+ //! moab interface
+ Interface *mbImpl;
+
+#ifdef USE_MPI
+ //! ParallelComm for master, slave meshes
+ ParallelComm *pcMaster, *pcSlave;
+#endif
+
+ //! material set numbers for fluid materials
+ std::set<int> fluidSetNos;
+
+ //! material set numbers for solid materials
+ std::set<int> solidSetNos;
+
+ //! sets defining master/slave meshes
+ EntityHandle masterSet, slaveSet;
+
+ //! sets in master/slave meshes
+ Range fluidSets[2], solidSets[2];
+
+ //! elements in master/slave meshes
+ Range fluidElems[2], solidElems[2];
+
+ //! filenames for master/slave meshes
+ std::string masterFileName, slaveFileName;
+
+ //! tag used for new positions
+ Tag xNew;
+
+ //! tag name used for new positions
+ std::string xNewName;
+};
+
+ //! add a set number
+inline ErrorCode DeformMeshRemap::add_set_no(int f_or_s, int set_no)
+{
+ std::set<int> *this_set;
+ switch (f_or_s) {
+ case FLUID:
+ this_set = &fluidSetNos; break;
+ case SOLID:
+ this_set = &solidSetNos; break;
+ default:
+ assert(false && "f_or_s should be FLUID or SOLID.");
+ return MB_FAILURE;
+ }
+
+ this_set->insert(set_no);
+
+ return MB_SUCCESS;
+}
+
+ //! remove a set number
+inline ErrorCode DeformMeshRemap::remove_set_no(int f_or_s, int set_no)
+{
+ std::set<int> *this_set;
+ switch (f_or_s) {
+ case FLUID:
+ this_set = &fluidSetNos; break;
+ case SOLID:
+ this_set = &solidSetNos; break;
+ default:
+ assert(false && "f_or_s should be FLUID or SOLID.");
+ return MB_FAILURE;
+ }
+ std::set<int>::iterator sit = this_set->find(set_no);
+ if (sit != this_set->end()) {
+ this_set->erase(*sit);
+ return MB_SUCCESS;
+ }
+
+ return MB_FAILURE;
+}
+
+ //! get the set numbers
+inline ErrorCode DeformMeshRemap::get_set_nos(int f_or_s, std::set<int> &set_nos) const
+{
+ const std::set<int> *this_set;
+ switch (f_or_s) {
+ case FLUID:
+ this_set = &fluidSetNos; break;
+ case SOLID:
+ this_set = &solidSetNos; break;
+ default:
+ assert(false && "f_or_s should be FLUID or SOLID.");
+ return MB_FAILURE;
+ }
+
+ set_nos = *this_set;
+
+ return MB_SUCCESS;
+}
+
+ErrorCode DeformMeshRemap::execute()
+{
+ // read master/slave files and get fluid/solid material sets
+ ErrorCode rval = read_file(MASTER, masterFileName, masterSet);
+ if (MB_SUCCESS != rval) return rval;
+
+ rval = read_file(SLAVE, slaveFileName, slaveSet);
+ if (MB_SUCCESS != rval) return rval;
+
+ Range src_elems = solidElems[MASTER];
+ src_elems.merge(fluidElems[MASTER]);
+ // locate slave vertices in master, orig coords; do this with a data coupler, so you can
+ // later interpolate
+ Range tgt_verts, tmp_range = solidElems[SLAVE];
+ tmp_range.merge(fluidElems[SLAVE]);
+ rval = mbImpl->get_adjacencies(tmp_range, 0, false, tgt_verts, Interface::UNION);
+ RR("Failed to get target verts.");
+
+
+ // initialize data coupler on source elements
+ DataCoupler dc_master(mbImpl, NULL, src_elems, 0);
+
+ // locate slave vertices, caching results in dc
+ rval = dc_master.locate_points(tgt_verts); RR("Point location of tgt verts failed.");
+ int num_located = dc_master.spatial_locator()->local_num_located();
+ if (num_located != (int)tgt_verts.size()) {
+ rval = MB_FAILURE;
+ std::cout << "Only " << num_located << " out of " << tgt_verts.size() << " target points successfully located." << std::endl;
+ return rval;
+ }
+
+ // deform the master's solid mesh, put results in a new tag
+ rval = deform_master(fluidElems[MASTER], solidElems[MASTER], "xnew"); RR("");
+
+ { // to isolate the lloyd smoother & delete when done
+
+ // smooth the master mesh
+ LloydSmoother ll(mbImpl, NULL, fluidElems[MASTER], xNew);
+ rval = ll.perform_smooth();
+ RR("Failed in lloyd smoothing.");
+ cout << "Lloyd smoothing required " << ll.num_its() << " iterations." << endl;
+ }
+
+ // map new locations to slave
+ // interpolate xNew to slave points
+ rval = dc_master.interpolate((int)DataCoupler::VOLUME, "xnew"); RR("Failed to interpolate target solution.");
+
+ // transfer xNew to coords, for master and slave
+ rval = write_to_coords(fluidElems[MASTER], xNew); RR("Failed writing tag to master fluid verts.");
+ rval = write_to_coords(tgt_verts, xNew); RR("Failed writing tag to slave verts.");
+
+ if (debug) {
+ std::string str;
+#ifdef USE_MPI
+ if (pcMaster && pcMaster->size() > 1)
+ str = "PARALLEL=WRITE_PART";
+#endif
+ rval = mbImpl->write_file("smoothed_master.vtk", NULL, str.c_str(), &masterSet, 1);
+#ifdef USE_MPI
+ str.clear();
+ if (pcSlave && pcSlave->size() > 1)
+ str = "PARALLEL=WRITE_PART";
+#endif
+ rval = mbImpl->write_file("slave_interp.vtk", NULL, str.c_str(), &slaveSet, 1);
+ }
+
+ return MB_SUCCESS;
+}
+
+std::string DeformMeshRemap::get_file_name(int m_or_s) const
+{
+ switch (m_or_s) {
+ case MASTER:
+ return masterFileName;
+ case SLAVE:
+ return slaveFileName;
+ default:
+ assert(false && "m_or_s should be MASTER or SLAVE.");
+ return std::string();
+ }
+}
+
+void DeformMeshRemap::set_file_name(int m_or_s, const std::string &name)
+{
+ switch (m_or_s) {
+ case MASTER:
+ masterFileName = name; break;
+ case SLAVE:
+ slaveFileName = name; break;
+ default:
+ assert(false && "m_or_s should be MASTER or SLAVE.");
+ }
+}
+
+DeformMeshRemap::DeformMeshRemap(Interface *impl, ParallelComm *master, ParallelComm *slave)
+ : mbImpl(impl), pcMaster(master), pcSlave(slave), masterSet(0), slaveSet(0), xNew(0), xNewName("xnew")
+{
+ if (!pcSlave && pcMaster)
+ pcSlave = pcMaster;
+}
+
+DeformMeshRemap::~DeformMeshRemap()
+{
+ // delete the tag
+ mbImpl->tag_delete(xNew);
+}
+
+int main(int argc, char **argv) {
+
+ ErrorCode rval;
+
+ ProgOptions po("Deformed mesh options");
+ po.addOpt<std::string> ("master,m", "Specify the master meshfile name" );
+ po.addOpt<std::string> ("slave,s", "Specify the slave meshfile name" );
+ po.parseCommandLine(argc, argv);
+ std::string foo;
+ string masterf, slavef;
+ if(!po.getOpt("master", &masterf))
+ masterf = string(MESH_DIR) + string("/rodquad.g");
+ if(!po.getOpt("slave", &slavef))
+ slavef = string(MESH_DIR) + string("/rodtri.g");
+
+ mb = new Core();
+
+ DeformMeshRemap *dfr;
+#ifdef USE_MPI
+ ParallelComm *pc = new ParallelComm(mb, MPI_COMM_WORLD);
+ dfr = new DeformMeshRemap(mb, pc);
+#else
+ dfr = new DeformMeshRemap(mb);
+#endif
+ dfr->set_file_name(DeformMeshRemap::MASTER, masterf);
+ dfr->set_file_name(DeformMeshRemap::SLAVE, slavef);
+ rval = dfr->add_set_no(DeformMeshRemap::SOLID, SOLID_SETNO); RR("Failed to add solid set no.");
+ rval = dfr->add_set_no(DeformMeshRemap::FLUID, FLUID_SETNO); RR("Failed to add fluid set no.");
+
+ rval = dfr->execute();
+
+ delete dfr;
+ delete mb;
+
+ return rval;
+}
+
+ErrorCode DeformMeshRemap::write_and_save(Range &ents, EntityHandle seth, Tag tagh, const char *filename)
+{
+ ErrorCode rval = write_to_coords(ents, tagh); RR("");
+ rval = mbImpl->write_file(filename, NULL, NULL, &seth, 1); RR("");
+ return rval;
+}
+
+ErrorCode DeformMeshRemap::write_to_coords(Range &elems, Tag tagh)
+{
+ // write the tag to coordinates
+ Range verts;
+ ErrorCode rval = mbImpl->get_adjacencies(elems, 0, false, verts, Interface::UNION);
+ RR("Failed to get adj vertices.");
+ std::vector<double> coords(3*verts.size());
+ rval = mbImpl->tag_get_data(tagh, verts, &coords[0]);
+ RR("Failed to get tag data.");
+ rval = mbImpl->set_coords(verts, &coords[0]);
+ RR("Failed to set coordinates.");
+ return MB_SUCCESS;
+}
+
+void deform_func(const BoundBox &bbox, double *xold, double *xnew)
+{
+/* Deformation function based on max delx and dely at top of rod
+ const double RODWIDTH = 0.2, RODHEIGHT = 0.5;
+ // function: origin is at middle base of rod, and is .5 high
+ // top of rod is (0,.55) on left and (.2,.6) on right
+ double delx = 0.5*RODWIDTH;
+
+ double xfrac = (xold[0] + .5*RODWIDTH)/RODWIDTH, yfrac = xold[1]/RODHEIGHT;
+ xnew[0] = xold[0] + yfrac * delx;
+ xnew[1] = xold[1] + yfrac * (1.0 + xfrac) * 0.05;
+*/
+
+/* Deformation function based on fraction of bounding box dimension in each direction */
+ double frac = 0.01; // taken from approximate relative deformation from LLNL Diablo of XX09 assys
+ CartVect *xo = reinterpret_cast<CartVect*>(xold), *xn = reinterpret_cast<CartVect*>(xnew);
+ CartVect disp = frac * (*xo - bbox.bMin);
+ *xn = *xo + disp;
+}
+
+ErrorCode DeformMeshRemap::deform_master(Range &fluid_elems, Range &solid_elems, const char *tag_name)
+{
+ // deform elements with an analytic function
+
+ // create the tag
+ ErrorCode rval = mbImpl->tag_get_handle((tag_name ? tag_name : ""), 3, MB_TYPE_DOUBLE, xNew, MB_TAG_CREAT|MB_TAG_DENSE);
+ RR("Failed to create xnew tag.");
+
+ // get all the vertices and coords in the fluid, set xnew to them
+ Range verts;
+ rval = mbImpl->get_adjacencies(fluid_elems, 0, false, verts, Interface::UNION);
+ RR("Failed to get vertices.");
+ std::vector<double> coords(3*verts.size(), 0.0);
+ rval = mbImpl->get_coords(verts, &coords[0]);
+ RR("Failed to get vertex coords.");
+ rval = mbImpl->tag_set_data(xNew, verts, &coords[0]);
+ RR("Failed to set xnew tag on fluid verts.");
+
+ // get the bounding box of the solid mesh
+ BoundBox bbox;
+ bbox.update(*mbImpl, solid_elems);
+
+ // get all the vertices and coords in the solid
+ verts.clear();
+ rval = mbImpl->get_adjacencies(solid_elems, 0, false, verts, Interface::UNION);
+ RR("Failed to get vertices.");
+ coords.resize(3*verts.size(), 0.0);
+ rval = mbImpl->get_coords(verts, &coords[0]);
+ RR("Failed to get vertex coords.");
+ unsigned int num_verts = verts.size();
+ for (unsigned int i = 0; i < num_verts; i++)
+ deform_func(bbox, &coords[3*i], &coords[3*i]);
+
+ // set the new tag to those coords
+ rval = mbImpl->tag_set_data(xNew, verts, &coords[0]);
+ RR("Failed to set tag data.");
+
+ return MB_SUCCESS;
+}
+
+ErrorCode DeformMeshRemap::read_file(int m_or_s, string &fname, EntityHandle &seth)
+{
+ // create meshset
+ ErrorCode rval = mbImpl->create_meshset(0, seth);
+ RR("Couldn't create master/slave set.");
+ ostringstream ostr;
+#ifdef USE_MPI
+ ParallelComm *pc = (m_or_s == MASTER ? pcMaster : pcSlave);
+ if (pc && pc->size() > 1) {
+ if (debug) ostr << "DEBUG_IO=1;CPUTIME;";
+ ostr << "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARALLEL_RESOLVE_SHARED_ENTS;"
+ << "PARALLEL_GHOSTS=2.0.1;PARALLEL_COMM=" << pc->get_id();
+ }
+#endif
+ rval = mbImpl->load_file(fname.c_str(), &seth, ostr.str().c_str());
+ RR("Couldn't load master/slave mesh.");
+
+ // get material sets for solid/fluid
+ Tag tagh;
+ rval = mbImpl->tag_get_handle(MATERIAL_SET_TAG_NAME, tagh); RR("Couldn't get material set tag name.");
+ for (std::set<int>::iterator sit = solidSetNos.begin(); sit != solidSetNos.end(); sit++) {
+ Range sets;
+ int set_no = *sit;
+ const void *setno_ptr = &set_no;
+ rval = mbImpl->get_entities_by_type_and_tag(seth, MBENTITYSET, &tagh, &setno_ptr, 1, sets);
+ if (sets.empty()) rval = MB_FAILURE;
+ RR("Couldn't get any solid sets.");
+ solidSets[m_or_s].merge(sets);
+ }
+
+ // get solid entities, and dimension
+ Range tmp_range;
+ for (Range::iterator rit = solidSets[m_or_s].begin(); rit != solidSets[m_or_s].end(); rit++) {
+ rval = mbImpl->get_entities_by_handle(*rit, tmp_range, true);
+ RR("Failed to get entities in solid.");
+ }
+ int dim = mbImpl->dimension_from_handle(*tmp_range.rbegin());
+ assert(dim > 0 && dim < 4);
+
+ solidElems[m_or_s] = tmp_range.subset_by_dimension(dim);
+
+ for (std::set<int>::iterator sit = fluidSetNos.begin(); sit != fluidSetNos.end(); sit++) {
+ Range sets;
+ int set_no = *sit;
+ const void *setno_ptr = &set_no;
+ rval = mbImpl->get_entities_by_type_and_tag(seth, MBENTITYSET, &tagh, &setno_ptr, 1, sets);
+ if (sets.empty()) rval = MB_FAILURE;
+ RR("Couldn't get any fluid sets.");
+ fluidSets[m_or_s].merge(sets);
+ }
+
+ // get fluid entities, and dimension
+ tmp_range.clear();
+ for (Range::iterator rit = fluidSets[m_or_s].begin(); rit != fluidSets[m_or_s].end(); rit++) {
+ rval = mbImpl->get_entities_by_handle(*rit, tmp_range, true);
+ RR("Failed to get entities in fluid.");
+ }
+
+ fluidElems[m_or_s] = tmp_range.subset_by_dimension(dim);
+
+ return MB_SUCCESS;
+}
diff --git a/examples/LloydRelaxation.cpp b/examples/LloydRelaxation.cpp
index 0cf9d1e..34002cd 100644
--- a/examples/LloydRelaxation.cpp
+++ b/examples/LloydRelaxation.cpp
@@ -14,8 +14,10 @@
* in the current directory (H5M format must be used since the file is written in parallel).
*/
-#include "moab/ParallelComm.hpp"
-#include "MBParallelConventions.h"
+#ifdef USE_MPI
+# include "moab/ParallelComm.hpp"
+# include "MBParallelConventions.h"
+#endif
#include "moab/Core.hpp"
#include "moab/Skinner.hpp"
#include "moab/CN.hpp"
@@ -30,7 +32,7 @@ string test_file_name = string(MESH_DIR) + string("/surfrandomtris-4part.h5m");
#define RC if (MB_SUCCESS != rval) return rval
-ErrorCode perform_lloyd_relaxation(ParallelComm *pc, Range &verts, Range &cells, Tag fixed,
+ErrorCode perform_lloyd_relaxation(Interface *mb, Range &verts, Range &cells, Tag fixed,
int num_its, int report_its);
int main(int argc, char **argv)
@@ -38,7 +40,9 @@ int main(int argc, char **argv)
int num_its = 10;
int report_its = 1;
+#ifdef USE_MPI
MPI_Init(&argc, &argv);
+#endif
// need option handling here for input filename
if (argc > 1){
@@ -49,10 +53,13 @@ int main(int argc, char **argv)
// get MOAB and ParallelComm instances
Interface *mb = new Core;
if (NULL == mb) return 1;
+ int nprocs = 1;
+
+#ifdef USE_MPI
// get the ParallelComm instance
- ParallelComm* pcomm = new ParallelComm(mb, MPI_COMM_WORLD);
- int nprocs = pcomm->size();
-
+ ParallelComm *pcomm = new ParallelComm(mb, MPI_COMM_WORLD);
+ nprocs = pcomm->size();
+#endif
string options;
if (nprocs > 1) // if reading in parallel, need to tell it how
options = "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARALLEL_RESOLVE_SHARED_ENTS;PARALLEL_GHOSTS=2.0.1;DEBUG_IO=0;DEBUG_PIO=0";
@@ -76,34 +83,45 @@ int main(int argc, char **argv)
// ok to mark non-owned skin vertices too, I won't move those anyway
// use MOAB's skinner class to find the skin
Skinner skinner(mb);
- rval = skinner.find_skin(faces, true, skin_verts); RC; // 'true' param indicates we want vertices back, not faces
+ rval = skinner.find_skin(0, faces, true, skin_verts); RC; // 'true' param indicates we want vertices back, not faces
std::vector<int> fix_tag(skin_verts.size(), 1); // initialized to 1 to indicate fixed
rval = mb->tag_set_data(fixed, skin_verts, &fix_tag[0]); RC;
// now perform the Lloyd relaxation
- rval = perform_lloyd_relaxation(pcomm, verts, faces, fixed, num_its, report_its); RC;
+ rval = perform_lloyd_relaxation(mb, verts, faces, fixed, num_its, report_its); RC;
// delete fixed tag, since we created it here
rval = mb->tag_delete(fixed); RC;
// output file, using parallel write
- rval = mb->write_file("lloydfinal.h5m", NULL, "PARALLEL=WRITE_PART"); RC;
+
+#ifdef USE_MPI
+ options = "PARALLEL=WRITE_PART";
+#endif
+
+ rval = mb->write_file("lloydfinal.h5m", NULL, options.c_str()); RC;
// delete MOAB instance
delete mb;
+#ifdef USE_MPI
MPI_Finalize();
+#endif
return 0;
}
-ErrorCode perform_lloyd_relaxation(ParallelComm *pcomm, Range &verts, Range &faces, Tag fixed,
+ErrorCode perform_lloyd_relaxation(Interface *mb, Range &verts, Range &faces, Tag fixed,
int num_its, int report_its)
{
ErrorCode rval;
- Interface *mb = pcomm->get_moab();
- int nprocs = pcomm->size();
+ int nprocs = 1;
+
+#ifdef USE_MPI
+ ParallelComm *pcomm = ParallelComm::get_pcomm(mb, 0);
+ nprocs = pcomm->size();
+#endif
// perform Lloyd relaxation:
// 1. setup: set vertex centroids from vertex coords; filter to owned verts; get fixed tags
@@ -120,8 +138,10 @@ ErrorCode perform_lloyd_relaxation(ParallelComm *pcomm, Range &verts, Range &fac
// filter verts down to owned ones and get fixed tag for them
Range owned_verts, shared_owned_verts;
if (nprocs > 1) {
+#ifdef USE_MPI
rval = pcomm->filter_pstatus(verts, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, &owned_verts);
if (rval != MB_SUCCESS) return rval;
+#endif
}
else
owned_verts = verts;
@@ -133,11 +153,13 @@ ErrorCode perform_lloyd_relaxation(ParallelComm *pcomm, Range &verts, Range &fac
vcentroids.resize(3*owned_verts.size());
rval = mb->tag_get_data(centroid, owned_verts, &vcentroids[0]); RC;
+#ifdef USE_MPI
// get shared owned verts, for exchanging tags
rval = pcomm->get_shared_entities(-1, shared_owned_verts, 0, false, true); RC;
// workaround: if no shared owned verts, put a non-shared one in the list, to prevent exchanging tags
// for all shared entities
if (shared_owned_verts.empty()) shared_owned_verts.insert(*verts.begin());
+#endif
// some declarations for later iterations
std::vector<double> fcentroids(3*faces.size()); // fcentroids for face centroids
@@ -195,16 +217,22 @@ ErrorCode perform_lloyd_relaxation(ParallelComm *pcomm, Range &verts, Range &fac
// 2c. exchange tags on owned verts
if (nprocs > 1) {
+#ifdef USE_MPI
rval = pcomm->exchange_tags(centroid, shared_owned_verts); RC;
+#endif
}
if (!(nit%report_its)) {
// global reduce for maximum delta, then report it
double global_max = mxdelta;
+ int myrank = 0;
+#ifdef USE_MPI
if (nprocs > 1)
MPI_Reduce(&mxdelta, &global_max, 1, MPI_DOUBLE, MPI_MAX, 0, pcomm->comm());
- if (1 == nprocs || !pcomm->rank())
+ myrank = pcomm->rank();
+#endif
+ if (1 == nprocs || !myrank)
cout << "Max delta = " << global_max << endl;
}
}
diff --git a/examples/makefile b/examples/makefile
index 5ae5af0..a6e3d1f 100644
--- a/examples/makefile
+++ b/examples/makefile
@@ -4,7 +4,7 @@ include ${MOAB_DIR}/lib/iMesh-Defs.inc
.SUFFIXES: .o .cpp .F90
-# MESH_DIR is the top-level MOAB source directory, used to locate mesh files that come with MOAB source
+# MESH_DIR is the directory containing mesh files that come with MOAB source
MESH_DIR="../MeshFiles/unittest"
EXAMPLES = HelloMOAB GetEntities SetsNTags StructuredMeshSimple DirectAccessWithHoles DirectAccessNoHoles
@@ -14,45 +14,48 @@ EXOIIEXAMPLES = TestExodusII
default: ${EXAMPLES}
-HelloMOAB : HelloMOAB.o
+HelloMOAB : HelloMOAB.o ${MOAB_LIBDIR}/libMOAB.la
${MOAB_CXX} -o $@ $< ${MOAB_LIBS_LINK}
-GetEntities: GetEntities.o
+GetEntities: GetEntities.o ${MOAB_LIBDIR}/libMOAB.la
${MOAB_CXX} -o $@ $< ${MOAB_LIBS_LINK}
-SetsNTags: SetsNTags.o
+SetsNTags: SetsNTags.o ${MOAB_LIBDIR}/libMOAB.la
${MOAB_CXX} -o $@ $< ${MOAB_LIBS_LINK}
-LloydRelaxation: LloydRelaxation.o
+LloydRelaxation: LloydRelaxation.o ${MOAB_LIBDIR}/libMOAB.la
${MOAB_CXX} -o $@ $< ${MOAB_LIBS_LINK}
-StructuredMeshSimple : StructuredMeshSimple.o
+StructuredMeshSimple : StructuredMeshSimple.o ${MOAB_LIBDIR}/libMOAB.la
${MOAB_CXX} -o $@ $< ${MOAB_LIBS_LINK}
-DirectAccessWithHoles: DirectAccessWithHoles.o
+DirectAccessWithHoles: DirectAccessWithHoles.o ${MOAB_LIBDIR}/libMOAB.la
${MOAB_CXX} -o $@ $< ${MOAB_LIBS_LINK}
-DirectAccessNoHoles: DirectAccessNoHoles.o
+DirectAccessNoHoles: DirectAccessNoHoles.o ${MOAB_LIBDIR}/libMOAB.la
${MOAB_CXX} -o $@ $< ${MOAB_LIBS_LINK}
-DirectAccessNoHolesF90: DirectAccessNoHolesF90.o
+DirectAccessNoHolesF90: DirectAccessNoHolesF90.o ${MOAB_LIBDIR}/libMOAB.la
${MOAB_CXX} -o $@ $< ${IMESH_LIBS}
-ReduceExchangeTags : ReduceExchangeTags.o
+ReduceExchangeTags : ReduceExchangeTags.o ${MOAB_LIBDIR}/libMOAB.la
${MOAB_CXX} -o $@ $< ${MOAB_LIBS_LINK}
-HelloParMOAB: HelloParMOAB.o
+HelloParMOAB: HelloParMOAB.o ${MOAB_LIBDIR}/libMOAB.la
${MOAB_CXX} -o $@ $< ${MOAB_LIBS_LINK}
-TestExodusII: TestExodusII.o
+TestExodusII: TestExodusII.o ${MOAB_LIBDIR}/libMOAB.la
${MOAB_CXX} -o $@ $< ${MOAB_LIBS_LINK}
-point_in_elem_search: point_in_elem_search.o
+point_in_elem_search: point_in_elem_search.o ${MOAB_LIBDIR}/libMOAB.la
${MOAB_CXX} -o $@ $< ${MOAB_LIBS_LINK}
PushParMeshIntoMoabF90: PushParMeshIntoMoabF90.o
${MOAB_CXX} -o $@ $< ${IMESH_LIBS} -lgfortran -L/usr/lib/openmpi/lib -lmpi_f90 -lmpi_f77 -lmpi -lopen-rte -lopen-pal -ldl
+DeformMeshRemap: DeformMeshRemap.o ${MOAB_LIBDIR}/libMOAB.la
+ ${MOAB_CXX} -o $@ $< ${MOAB_LIBS_LINK} -lmbcoupler ${MOAB_LIBS_LINK}
+
clean:
rm -rf *.o ${EXAMPLES} ${PAREXAMPLES} ${EXOIIEXAMPLES}
diff --git a/src/AdaptiveKDTree.cpp b/src/AdaptiveKDTree.cpp
index f66acbc..1324c30 100644
--- a/src/AdaptiveKDTree.cpp
+++ b/src/AdaptiveKDTree.cpp
@@ -1303,7 +1303,8 @@ namespace moab {
ErrorCode AdaptiveKDTree::point_search(const double *point,
EntityHandle& leaf_out,
- double tol,
+ const double iter_tol,
+ const double inside_tol,
bool *multiple_leaves,
EntityHandle *start_node,
CartVect *params)
@@ -1322,7 +1323,7 @@ namespace moab {
treeStats.nodesVisited++;
ErrorCode rval = get_bounding_box(box, &node);
if (MB_SUCCESS != rval) return rval;
- if (!box.contains_point(point, tol)) return MB_SUCCESS;
+ if (!box.contains_point(point, iter_tol)) return MB_SUCCESS;
rval = moab()->get_child_meshsets( node, children );
if (MB_SUCCESS != rval)
@@ -1346,7 +1347,7 @@ namespace moab {
treeStats.leavesVisited++;
if (myEval && params) {
- rval = myEval->find_containing_entity(node, point, tol,
+ rval = myEval->find_containing_entity(node, point, iter_tol, inside_tol,
leaf_out, params->array(), &treeStats.leafObjectTests);
if (MB_SUCCESS != rval) return rval;
}
@@ -1358,7 +1359,8 @@ namespace moab {
ErrorCode AdaptiveKDTree::point_search(const double *point,
AdaptiveKDTreeIter& leaf_it,
- double tol,
+ const double iter_tol,
+ const double /*inside_tol*/,
bool *multiple_leaves,
EntityHandle *start_node)
{
@@ -1372,7 +1374,7 @@ namespace moab {
leaf_it.mBox[1] = boundBox.bMax;
// test that point is inside tree
- if (!boundBox.contains_point(point, tol)) {
+ if (!boundBox.contains_point(point, iter_tol)) {
treeStats.nodesVisited++;
return MB_ENTITY_NOT_FOUND;
}
@@ -1423,7 +1425,8 @@ namespace moab {
ErrorCode AdaptiveKDTree::distance_search(const double from_point[3],
const double distance,
std::vector<EntityHandle>& result_list,
- double tol,
+ const double iter_tol,
+ const double inside_tol,
std::vector<double> *result_dists,
std::vector<CartVect> *result_params,
EntityHandle *tree_root)
@@ -1445,7 +1448,7 @@ namespace moab {
// (zero if inside box)
BoundBox box;
rval = get_bounding_box(box);
- if (MB_SUCCESS == rval && !box.contains_point(from_point, tol)) {
+ if (MB_SUCCESS == rval && !box.contains_point(from_point, iter_tol)) {
treeStats.nodesVisited++;
return MB_SUCCESS;
}
@@ -1484,7 +1487,7 @@ namespace moab {
if (myEval && result_params) {
EntityHandle ent;
CartVect params;
- rval = myEval->find_containing_entity(node.handle, from_point, tol,
+ rval = myEval->find_containing_entity(node.handle, from_point, iter_tol, inside_tol,
ent, params.array(), &treeStats.leafObjectTests);
if (MB_SUCCESS != rval) return rval;
else if (ent) {
@@ -1802,7 +1805,7 @@ namespace moab {
// the same distance from the input point as the current closest
// point is.
CartVect diff = closest_pt - from;
- rval = distance_search(from_coords, sqrt(diff%diff), leaves, 0.0, NULL, NULL, &tree_root);
+ rval = distance_search(from_coords, sqrt(diff%diff), leaves, 1.0e-10, 1.0e-6, NULL, NULL, &tree_root);
if (MB_SUCCESS != rval) return rval;
// Check any close leaves to see if they contain triangles that
@@ -1834,7 +1837,7 @@ namespace moab {
// get leaves of tree that intersect sphere
assert(tree_root);
- rval = distance_search(center, radius, leaves, 0.0, NULL, NULL, &tree_root);
+ rval = distance_search(center, radius, leaves, 1.0e-10, 1.0e-6, NULL, NULL, &tree_root);
if (MB_SUCCESS != rval) return rval;
// search each leaf for triangles intersecting sphere
diff --git a/src/BVHTree.cpp b/src/BVHTree.cpp
index 8903070..b65eae4 100644
--- a/src/BVHTree.cpp
+++ b/src/BVHTree.cpp
@@ -435,7 +435,8 @@ namespace moab
ErrorCode BVHTree::find_point(const std::vector<double> &point,
const unsigned int &index,
- const double tol,
+ const double iter_tol,
+ const double inside_tol,
std::pair<EntityHandle, CartVect> &result)
{
if (index == 0) treeStats.numTraversals++;
@@ -453,7 +454,7 @@ namespace moab
for(Range::iterator i = entities.begin(); i != entities.end(); i++) {
treeStats.leafObjectTests++;
myEval->set_ent_handle(*i);
- myEval->reverse_eval(&point[0], tol, params.array(), &is_inside);
+ myEval->reverse_eval(&point[0], iter_tol, inside_tol, params.array(), &is_inside);
if (is_inside) {
result.first = *i;
result.second = params;
@@ -469,11 +470,11 @@ namespace moab
rval = mbImpl->get_child_meshsets(startSetHandle+index, children);
if (MB_SUCCESS != rval || children.size() != 2) return rval;
- if((node.Lmax+tol) < (node.Rmin-tol)){
- if(point[node.dim] <= (node.Lmax + tol))
- return find_point(point, children[0]-startSetHandle, tol, result);
- else if(point[ node.dim] >= (node.Rmin - tol))
- return find_point(point, children[1]-startSetHandle, tol, result);
+ if((node.Lmax+iter_tol) < (node.Rmin-iter_tol)){
+ if(point[node.dim] <= (node.Lmax + iter_tol))
+ return find_point(point, children[0]-startSetHandle, iter_tol, inside_tol, result);
+ else if(point[ node.dim] >= (node.Rmin - iter_tol))
+ return find_point(point, children[1]-startSetHandle, iter_tol, inside_tol, result);
result = std::make_pair(0, CartVect(&point[0])); //point lies in empty space.
return MB_SUCCESS;
}
@@ -482,11 +483,11 @@ namespace moab
//left of Rmin, you must be on the left
//we can't be sure about the boundaries since the boxes overlap
//this was a typo in the paper which caused pain.
- if(point[node.dim] < (node.Rmin - tol))
- return find_point(point, children[0]-startSetHandle, tol, result);
+ if(point[node.dim] < (node.Rmin - iter_tol))
+ return find_point(point, children[0]-startSetHandle, iter_tol, inside_tol, result);
//if you are on the right Lmax, you must be on the right
- else if(point[ node.dim] > (node.Lmax+tol))
- return find_point(point, children[1]-startSetHandle, tol, result);
+ else if(point[ node.dim] > (node.Lmax+iter_tol))
+ return find_point(point, children[1]-startSetHandle, iter_tol, inside_tol, result);
/* pg5 of paper
* However, instead of always traversing either subtree
@@ -501,22 +502,22 @@ namespace moab
//branch predicted..
//bool dir = (point[ node.dim] - node.Rmin) <=
// (node.Lmax - point[ node.dim]);
- find_point(point, children[0]-startSetHandle, tol, result);
+ find_point(point, children[0]-startSetHandle, iter_tol, inside_tol, result);
if(result.first == 0){
- return find_point(point, children[1]-startSetHandle, tol, result);
+ return find_point(point, children[1]-startSetHandle, iter_tol, inside_tol, result);
}
return MB_SUCCESS;
}
- EntityHandle BVHTree::bruteforce_find(const double *point, const double tol) {
+ EntityHandle BVHTree::bruteforce_find(const double *point, const double iter_tol, const double inside_tol) {
treeStats.numTraversals++;
CartVect params;
for(unsigned int i = 0; i < myTree.size(); i++) {
- if(myTree[i].dim != 3 || !myTree[i].box.contains_point(point, tol)) continue;
+ if(myTree[i].dim != 3 || !myTree[i].box.contains_point(point, iter_tol)) continue;
if (myEval) {
EntityHandle entity = 0;
treeStats.leavesVisited++;
- ErrorCode rval = myEval->find_containing_entity(startSetHandle+i, point, tol,
+ ErrorCode rval = myEval->find_containing_entity(startSetHandle+i, point, iter_tol, inside_tol,
entity, params.array(), &treeStats.leafObjectTests);
if (entity) return entity;
else if (MB_SUCCESS != rval) return 0;
@@ -542,7 +543,8 @@ namespace moab
ErrorCode BVHTree::point_search(const double *point,
EntityHandle& leaf_out,
- double tol,
+ const double iter_tol,
+ const double inside_tol,
bool *multiple_leaves,
EntityHandle *start_node,
CartVect *params)
@@ -572,7 +574,7 @@ namespace moab
// test box of this node
ErrorCode rval = get_bounding_box(box, &this_set);
if (MB_SUCCESS != rval) return rval;
- if (!box.contains_point(point, tol)) continue;
+ if (!box.contains_point(point, iter_tol)) continue;
// else if not a leaf, test children & put on list
else if (myTree[ind].dim != 3) {
@@ -581,7 +583,7 @@ namespace moab
continue;
}
else if (myTree[ind].dim == 3 && myEval && params) {
- rval = myEval->find_containing_entity(startSetHandle+ind, point, tol,
+ rval = myEval->find_containing_entity(startSetHandle+ind, point, iter_tol, inside_tol,
leaf_out, params->array(), &treeStats.leafObjectTests);
if (leaf_out || MB_SUCCESS != rval) return rval;
}
@@ -599,7 +601,8 @@ namespace moab
ErrorCode BVHTree::distance_search(const double from_point[3],
const double distance,
std::vector<EntityHandle>& result_list,
- double params_tol,
+ const double iter_tol,
+ const double inside_tol,
std::vector<double> *result_dists,
std::vector<CartVect> *result_params,
EntityHandle *tree_root)
@@ -653,7 +656,7 @@ namespace moab
if (myEval && result_params) {
EntityHandle ent;
CartVect params;
- rval = myEval->find_containing_entity(startSetHandle+ind, from_point, params_tol,
+ rval = myEval->find_containing_entity(startSetHandle+ind, from_point, iter_tol, inside_tol,
ent, params.array(), &treeStats.leafObjectTests);
if (MB_SUCCESS != rval) return rval;
else if (ent) {
diff --git a/src/LloydSmoother.cpp b/src/LloydSmoother.cpp
new file mode 100644
index 0000000..b1614d2
--- /dev/null
+++ b/src/LloydSmoother.cpp
@@ -0,0 +1,232 @@
+#include "moab/LloydSmoother.hpp"
+#include "moab/Error.hpp"
+#include "moab/Skinner.hpp"
+#include "moab/CN.hpp"
+#include "moab/CartVect.hpp"
+#include "moab/BoundBox.hpp"
+
+#ifdef USE_MPI
+#include "moab/ParallelComm.hpp"
+#include "MBParallelConventions.h"
+#endif
+
+#include <iostream>
+
+#define RR(a) if (MB_SUCCESS != rval) { \
+ errorHandler->set_last_error(a); \
+ return rval;}
+
+namespace moab {
+
+LloydSmoother::LloydSmoother(Interface *impl, ParallelComm *pc, Range &elms, Tag ctag,
+ Tag ftag, double at, double rt)
+ : mbImpl(impl), myPcomm(pc), myElems(elms), coordsTag(ctag), fixedTag(ftag), absTol(at), relTol(rt),
+ reportIts(0), numIts(0), iCreatedTag(false)
+{
+ ErrorCode rval = mbImpl->query_interface(errorHandler);
+ if (rval) throw rval;
+}
+
+LloydSmoother::~LloydSmoother()
+{
+ if (iCreatedTag && fixedTag) {
+ ErrorCode rval = mbImpl->tag_delete(fixedTag);
+ if (rval) throw rval;
+ }
+}
+
+ErrorCode LloydSmoother::perform_smooth()
+{
+ ErrorCode rval;
+
+ // first figure out tolerance to use
+ if (0 > absTol) {
+ // no tolerance set - get one relative to bounding box around elements
+ BoundBox bb;
+ rval = bb.update(*mbImpl, myElems);
+ RR("Failed to compute bounding box around elements.");
+ absTol = relTol * bb.diagonal_length();
+ }
+
+ // initialize if we need to
+ rval = initialize();
+ RR("Failed to initialize.");
+
+ // get all vertices
+ Range verts;
+ rval = mbImpl->get_adjacencies(myElems, 0, false, verts, Interface::UNION);
+ RR("Failed to get all vertices.");
+
+ // perform Lloyd relaxation:
+ // 1. setup: set vertex centroids from vertex coords; filter to owned verts; get fixed tags
+
+ // get all verts coords into tag; don't need to worry about filtering out fixed verts,
+ // we'll just be setting to their fixed coords
+ std::vector<double> vcentroids(3*verts.size());
+ if (!coordsTag) {
+ rval = mbImpl->get_coords(verts, &vcentroids[0]); RR("Failed to get vert coords.");
+ }
+ else {
+ rval = mbImpl->tag_get_data(coordsTag, verts, &vcentroids[0]); RR("Failed to get vert coords tag values.");
+ }
+
+ Tag centroid;
+ rval = mbImpl->tag_get_handle("", 3, MB_TYPE_DOUBLE, centroid, MB_TAG_CREAT | MB_TAG_DENSE);
+ RR("Couldn't get tag handle.");
+ rval = mbImpl->tag_set_data(centroid, verts, &vcentroids[0]); RR("Failed setting centroid tag.");
+
+ Range owned_verts, shared_owned_verts;
+#ifdef USE_MPI
+ // filter verts down to owned ones and get fixed tag for them
+ if (myPcomm && myPcomm->size() > 1) {
+ rval = myPcomm->filter_pstatus(verts, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, &owned_verts);
+ RR("Failed to filter on pstatus.");
+ // get shared owned verts, for exchanging tags
+ rval = myPcomm->filter_pstatus(owned_verts, PSTATUS_SHARED, PSTATUS_AND, -1, &shared_owned_verts);
+ RR("Failed to filter for shared owned.");
+ // workaround: if no shared owned verts, put a non-shared one in the list, to prevent exchanging tags
+ // for all shared entities
+ if (shared_owned_verts.empty()) shared_owned_verts.insert(*verts.begin());
+ }
+ else
+ owned_verts = verts;
+#else
+ owned_verts = verts;
+#endif
+
+ std::vector<unsigned char> fix_tag(owned_verts.size());
+ rval = mbImpl->tag_get_data(fixedTag, owned_verts, &fix_tag[0]); RR("Failed to get fixed tag.");
+
+ // now fill vcentroids array with positions of just owned vertices, since those are the ones
+ // we're actually computing
+ vcentroids.resize(3*owned_verts.size());
+ rval = mbImpl->tag_get_data(centroid, owned_verts, &vcentroids[0]); RR("Failed to get centroid tag.");
+
+ // some declarations for later iterations
+ std::vector<double> fcentroids(3*myElems.size()); // fcentroids for element centroids
+ std::vector<double> ctag(3*CN::MAX_NODES_PER_ELEMENT); // temporary coordinate storage for verts bounding an element
+ const EntityHandle *conn; // const ptr & size to elem connectivity
+ int nconn;
+ Range::iterator eit, vit; // for iterating over elems, verts
+ int e, v; // for indexing into centroid vectors
+ std::vector<EntityHandle> adj_elems; // used in vertex iteration
+
+ // 2. while !converged
+ double resid = DBL_MAX;
+ numIts = 0;
+ while (resid > absTol) {
+ numIts++;
+ resid = 0.0;
+
+ // 2a. foreach elem: centroid = sum(vertex centroids)/num_verts_in_cell
+ for (eit = myElems.begin(), e = 0; eit != myElems.end(); eit++, e++) {
+ // get verts for this elem
+ rval = mbImpl->get_connectivity(*eit, conn, nconn); RR("Failed to get connectivity.");
+ // get centroid tags for those verts
+ rval = mbImpl->tag_get_data(centroid, conn, nconn, &ctag[0]); RR("Failed to get centroid.");
+ fcentroids[3*e+0] = fcentroids[3*e+1] = fcentroids[3*e+2] = 0.0;
+ for (v = 0; v < nconn; v++) {
+ fcentroids[3*e+0] += ctag[3*v+0];
+ fcentroids[3*e+1] += ctag[3*v+1];
+ fcentroids[3*e+2] += ctag[3*v+2];
+ }
+ for (v = 0; v < 3; v++) fcentroids[3*e+v] /= nconn;
+ }
+ rval = mbImpl->tag_set_data(centroid, myElems, &fcentroids[0]); RR("Failed to set elem centroid.");
+
+ // 2b. foreach owned vertex:
+ for (vit = owned_verts.begin(), v = 0; vit != owned_verts.end(); vit++, v++) {
+ // if !fixed
+ if (fix_tag[v]) continue;
+ // vertex centroid = sum(cell centroids)/ncells
+ adj_elems.clear();
+ rval = mbImpl->get_adjacencies(&(*vit), 1, 2, false, adj_elems); RR("Failed getting adjs.");
+ rval = mbImpl->tag_get_data(centroid, &adj_elems[0], adj_elems.size(), &fcentroids[0]);
+ RR("Failed to get elem centroid.");
+ double vnew[] = {0.0, 0.0, 0.0};
+ for (e = 0; e < (int)adj_elems.size(); e++) {
+ vnew[0] += fcentroids[3*e+0];
+ vnew[1] += fcentroids[3*e+1];
+ vnew[2] += fcentroids[3*e+2];
+ }
+ for (e = 0; e < 3; e++) vnew[e] /= adj_elems.size();
+ double delta = (CartVect(vnew)-CartVect(&vcentroids[3*v])).length();
+ resid = (v ? std::max(resid, delta) : delta);
+ for (e = 0; e < 3; e++) vcentroids[3*v+e] = vnew[e];
+ }
+
+ // set the centroid tag; having them only in vcentroids array isn't enough, as vertex centroids are
+ // accessed randomly in loop over faces
+ rval = mbImpl->tag_set_data(centroid, owned_verts, &vcentroids[0]); RR("Failed to set vertex centroid.");
+
+#ifdef USE_MPI
+ // 2c. exchange tags on owned verts
+ if (myPcomm && myPcomm->size() > 1) {
+ rval = myPcomm->exchange_tags(centroid, shared_owned_verts); RR("Failed to exchange tags.");
+ }
+#endif
+
+ if (reportIts && !(numIts%reportIts)) {
+ double global_max = resid;
+#ifdef USE_MPI
+ // global reduce for maximum delta, then report it
+ if (myPcomm && myPcomm->size() > 1)
+ MPI_Reduce(&resid, &global_max, 1, MPI_DOUBLE, MPI_MAX, 0, myPcomm->comm());
+ if (!myPcomm || !myPcomm->rank())
+#endif
+ std::cout << "Max residual = " << global_max << std::endl;
+ }
+
+ } // end while
+
+ // write the tag back onto vertex coordinates
+ if (!coordsTag) {
+ rval = mbImpl->set_coords(owned_verts, &vcentroids[0]); RR("Failed to set vertex coords.");
+ }
+ else {
+ rval = mbImpl->tag_set_data(coordsTag, owned_verts, &vcentroids[0]);
+ RR("Failed to set vert coords tag values.");
+ }
+
+ return MB_SUCCESS;
+}
+
+ErrorCode LloydSmoother::initialize()
+{
+ ErrorCode rval = MB_SUCCESS;
+
+ // first, check for tag; if we don't have it, make one and mark skin as fixed
+ if (!fixedTag) {
+ unsigned char fixed = 0x0;
+ rval = mbImpl->tag_get_handle("", 1, MB_TYPE_OPAQUE, fixedTag, MB_TAG_DENSE|MB_TAG_CREAT, &fixed);
+ RR("Trouble making fixed tag.");
+ iCreatedTag = true;
+
+ // get the skin; get facets, because we might need to filter on shared entities
+ Skinner skinner(mbImpl);
+ Range skin, skin_verts;
+ rval = skinner.find_skin(0, myElems, false, skin);
+ RR("Unable to find skin.");
+
+#ifdef USE_MPI
+ // need to do a little extra if we're working in parallel
+ if (myPcomm) {
+ // filter out ghost and interface facets
+ rval = myPcomm->filter_pstatus(skin, PSTATUS_GHOST|PSTATUS_INTERFACE, PSTATUS_NOT);
+ RR("Failed to filter on shared status.");
+ }
+#endif
+ // get the vertices from those entities
+ rval = mbImpl->get_adjacencies(skin, 0, false, skin_verts, Interface::UNION);
+ RR("Trouble getting vertices.");
+
+ // mark them fixed
+ std::vector<unsigned char> marks(skin.size(), 0x1);
+ rval = mbImpl->tag_set_data(fixedTag, skin_verts, &marks[0]);
+ RR("Unable to set tag on skin.");
+ }
+
+ return MB_SUCCESS;
+}
+
+}
diff --git a/src/LocalDiscretization/ElemEvaluator.cpp b/src/LocalDiscretization/ElemEvaluator.cpp
index b132da9..4ca3943 100644
--- a/src/LocalDiscretization/ElemEvaluator.cpp
+++ b/src/LocalDiscretization/ElemEvaluator.cpp
@@ -6,6 +6,7 @@
// need to include eval set types here to support get_eval_set; alternative would be to have some
// type of registration, but we'd still need static registration for the built-in types
+#include "moab/LinearTri.hpp"
#include "moab/LinearQuad.hpp"
#include "moab/LinearTet.hpp"
#include "moab/LinearHex.hpp"
@@ -16,12 +17,12 @@
namespace moab {
ErrorCode EvalSet::evaluate_reverse(EvalFcn eval, JacobianFcn jacob, InsideFcn inside_f,
const double *posn, const double *verts, const int nverts,
- const int ndim, const double tol, double *work,
- double *params, bool *inside) {
+ const int ndim, const double iter_tol, const double inside_tol,
+ double *work, double *params, bool *inside) {
// TODO: should differentiate between epsilons used for
// Newton Raphson iteration, and epsilons used for curved boundary geometry errors
// right now, fix the tolerance used for NR
- const double error_tol_sqr = tol*tol;
+ const double error_tol_sqr = iter_tol*iter_tol;
CartVect *cvparams = reinterpret_cast<CartVect*>(params);
const CartVect *cvposn = reinterpret_cast<const CartVect*>(posn);
@@ -30,41 +31,48 @@ namespace moab {
CartVect new_pos;
// evaluate that first guess to get a new position
- ErrorCode rval = (*eval)(cvparams->array(), verts, ndim, ndim, work, new_pos.array());
+ ErrorCode rval = (*eval)(cvparams->array(), verts, ndim,
+ 3, // hardwire to num_tuples to 3 since the field is coords
+ work, new_pos.array());
if (MB_SUCCESS != rval) return rval;
// residual is diff between old and new pos; need to minimize that
CartVect res = new_pos - *cvposn;
Matrix3 J;
+ bool dum, *tmp_inside = (inside ? inside : &dum);
int iters=0;
// while |res| larger than tol
while (res % res > error_tol_sqr) {
if(++iters>10) {
- if (inside) {
- // if we haven't converged but we're outside, that's defined as success
- *inside = (*inside_f)(params, ndim, tol);
- if (!(*inside)) return MB_SUCCESS;
- }
- return MB_FAILURE;
+ // if we haven't converged but we're outside, that's defined as success
+ *tmp_inside = (*inside_f)(params, ndim, inside_tol);
+ if (!(*tmp_inside)) return MB_SUCCESS;
+ else return MB_FAILURE;
}
// get jacobian at current params
rval = (*jacob)(cvparams->array(), verts, nverts, ndim, work, J[0]);
double det = J.determinant();
- assert(det > std::numeric_limits<double>::epsilon());
+ if (det < std::numeric_limits<double>::epsilon()) {
+ *tmp_inside = (*inside_f)(params, ndim, inside_tol);
+ if (!(*tmp_inside)) return MB_SUCCESS;
+ else return MB_FAILURE;
+ }
// new params tries to eliminate residual
*cvparams -= J.inverse(1.0/det) * res;
// get the new forward-evaluated position, and its difference from the target pt
- rval = (*eval)(params, verts, ndim, ndim, work, new_pos.array());
+ rval = (*eval)(params, verts, ndim,
+ 3, // hardwire to num_tuples to 3 since the field is coords
+ work, new_pos.array());
if (MB_SUCCESS != rval) return rval;
res = new_pos - *cvposn;
}
if (inside)
- *inside = (*inside_f)(params, ndim, tol);
+ *inside = (*inside_f)(params, ndim, inside_tol);
return MB_SUCCESS;
}// Map::evaluate_reverse()
@@ -85,6 +93,7 @@ namespace moab {
case MBEDGE:
break;
case MBTRI:
+ if (LinearTri::compatible(tp, num_vertices, eval_set)) return MB_SUCCESS;
break;
case MBQUAD:
if (LinearQuad::compatible(tp, num_vertices, eval_set)) return MB_SUCCESS;
@@ -105,9 +114,9 @@ namespace moab {
return MB_NOT_IMPLEMENTED;
}
- ErrorCode ElemEvaluator::find_containing_entity(Range &entities, const double *point, double tol,
- EntityHandle &containing_ent, double *params,
- unsigned int *num_evals)
+ ErrorCode ElemEvaluator::find_containing_entity(Range &entities, const double *point, const double iter_tol,
+ const double inside_tol, EntityHandle &containing_ent,
+ double *params, unsigned int *num_evals)
{
bool is_inside;
ErrorCode rval = MB_SUCCESS;
@@ -116,7 +125,7 @@ namespace moab {
for(i = entities.begin(); i != entities.end(); i++) {
nevals++;
set_ent_handle(*i);
- rval = reverse_eval(point, tol, params, &is_inside);
+ rval = reverse_eval(point, iter_tol, inside_tol, params, &is_inside);
if (MB_SUCCESS != rval) return rval;
if (is_inside) break;
}
diff --git a/src/LocalDiscretization/LinearHex.cpp b/src/LocalDiscretization/LinearHex.cpp
index 14de2d9..5e730af 100644
--- a/src/LocalDiscretization/LinearHex.cpp
+++ b/src/LocalDiscretization/LinearHex.cpp
@@ -96,10 +96,12 @@ namespace moab
ErrorCode LinearHex::reverseEvalFcn(EvalFcn eval, JacobianFcn jacob, InsideFcn ins,
const double *posn, const double *verts, const int nverts, const int ndim,
- const double tol, double *work, double *params, bool *is_inside)
+ const double iter_tol, const double inside_tol, double *work,
+ double *params, bool *is_inside)
{
assert(posn && verts);
- return EvalSet::evaluate_reverse(eval, jacob, ins, posn, verts, nverts, ndim, tol, work, params, is_inside);
+ return EvalSet::evaluate_reverse(eval, jacob, ins, posn, verts, nverts, ndim, iter_tol, inside_tol, work,
+ params, is_inside);
}
bool LinearHex::insideFcn(const double *params, const int ndim, const double tol)
diff --git a/src/LocalDiscretization/LinearQuad.cpp b/src/LocalDiscretization/LinearQuad.cpp
index ddd4918..63c9959 100644
--- a/src/LocalDiscretization/LinearQuad.cpp
+++ b/src/LocalDiscretization/LinearQuad.cpp
@@ -16,23 +16,23 @@ namespace moab
*/
const double LinearQuad::gauss[1][2] = { { 2.0, 0.0 } };
- ErrorCode LinearQuad::jacobianFcn(const double *params, const double *verts, const int /*nverts*/, const int ndim,
+ ErrorCode LinearQuad::jacobianFcn(const double *params, const double *verts, const int /*nverts*/, const int /*ndim*/,
double *, double *result)
{
Matrix3 *J = reinterpret_cast<Matrix3*>(result);
*J = Matrix3(0.0);
for (unsigned i = 0; i < 4; ++i) {
- const double params_p = 1 + params[0]*corner[i][0];
+ const double xi_p = 1 + params[0]*corner[i][0];
const double eta_p = 1 + params[1]*corner[i][1];
- const double dNi_dparams = corner[i][0] * eta_p;
- const double dNi_deta = corner[i][1] * params_p;
- (*J)(0,0) += dNi_dparams * verts[i*ndim+0];
- (*J)(1,0) += dNi_dparams * verts[i*ndim+1];
- (*J)(0,1) += dNi_deta * verts[i*ndim+0];
- (*J)(1,1) += dNi_deta * verts[i*ndim+1];
+ const double dNi_dxi = corner[i][0] * eta_p;
+ const double dNi_deta = corner[i][1] * xi_p;
+ (*J)(0,0) += dNi_dxi * verts[i*3+0];
+ (*J)(1,0) += dNi_dxi * verts[i*3+1];
+ (*J)(0,1) += dNi_deta * verts[i*3+0];
+ (*J)(1,1) += dNi_deta * verts[i*3+1];
}
- (*J)(2,2) = 1.0; /* to make sure the Jacobian determinant is non-zero */
(*J) *= 0.25;
+ (*J)(2,2) = 1.0; /* to make sure the Jacobian determinant is non-zero */
return MB_SUCCESS;
}// LinearQuad::jacobian()
@@ -62,7 +62,7 @@ namespace moab
for(unsigned int j2 = 0; j2 < LinearQuad::gauss_count; ++j2) {
x[1] = LinearQuad::gauss[j2][1];
double w2 = LinearQuad::gauss[j2][0];
- rval = evalFcn(x.array(),field, ndim, num_tuples, NULL, tmp_result);
+ rval = evalFcn(x.array(), field, ndim, num_tuples, NULL, tmp_result);
if (MB_SUCCESS != rval) return rval;
rval = jacobianFcn(x.array(), verts, nverts, ndim, work, J[0]);
if (MB_SUCCESS != rval) return rval;
@@ -75,9 +75,11 @@ namespace moab
ErrorCode LinearQuad::reverseEvalFcn(EvalFcn eval, JacobianFcn jacob, InsideFcn ins,
const double *posn, const double *verts, const int nverts, const int ndim,
- const double tol, double *work, double *params, bool *is_inside)
+ const double iter_tol, const double inside_tol, double *work,
+ double *params, bool *is_inside)
{
- return EvalSet::evaluate_reverse(eval, jacob, ins, posn, verts, nverts, ndim, tol, work, params, is_inside);
+ return EvalSet::evaluate_reverse(eval, jacob, ins, posn, verts, nverts, ndim, iter_tol, inside_tol, work,
+ params, is_inside);
}
bool LinearQuad::insideFcn(const double *params, const int ndim, const double tol)
diff --git a/src/LocalDiscretization/LinearTet.cpp b/src/LocalDiscretization/LinearTet.cpp
index 705dbee..4965cfc 100644
--- a/src/LocalDiscretization/LinearTet.cpp
+++ b/src/LocalDiscretization/LinearTet.cpp
@@ -74,27 +74,29 @@ namespace moab
ErrorCode LinearTet::reverseEvalFcn(EvalFcn eval, JacobianFcn jacob, InsideFcn ins,
const double *posn, const double *verts, const int nverts, const int ndim,
- const double tol, double *work, double *params, bool *is_inside)
+ const double iter_tol, const double inside_tol, double *work,
+ double *params, bool *is_inside)
{
assert(posn && verts);
- return evaluate_reverse(eval, jacob, ins, posn, verts, nverts, ndim, tol, work, params, is_inside);
+ return evaluate_reverse(eval, jacob, ins, posn, verts, nverts, ndim, iter_tol, inside_tol,
+ work, params, is_inside);
}
bool LinearTet::insideFcn(const double *params, const int , const double tol)
{
return (params[0] >= -1.0-tol && params[1] >= -1.0-tol && params[2] >= -1.0-tol &&
- params[0] + params[1] + params[2] <= -1.0+tol);
+ params[0] + params[1] + params[2] <= 1.0+tol);
}
ErrorCode LinearTet::evaluate_reverse(EvalFcn eval, JacobianFcn jacob, InsideFcn inside_f,
const double *posn, const double *verts, const int nverts,
- const int ndim, const double tol, double *work,
- double *params, bool *inside) {
+ const int ndim, const double iter_tol, const double inside_tol,
+ double *work, double *params, bool *inside) {
// TODO: should differentiate between epsilons used for
// Newton Raphson iteration, and epsilons used for curved boundary geometry errors
// right now, fix the tolerance used for NR
- const double error_tol_sqr = tol*tol;
+ const double error_tol_sqr = iter_tol*iter_tol;
CartVect *cvparams = reinterpret_cast<CartVect*>(params);
const CartVect *cvposn = reinterpret_cast<const CartVect*>(posn);
@@ -138,7 +140,7 @@ namespace moab
}
if (inside)
- *inside = (*inside_f)(params, ndim, tol);
+ *inside = (*inside_f)(params, ndim, inside_tol);
return MB_SUCCESS;
}// Map::evaluate_reverse()
diff --git a/src/LocalDiscretization/LinearTri.cpp b/src/LocalDiscretization/LinearTri.cpp
new file mode 100644
index 0000000..8fe9d23
--- /dev/null
+++ b/src/LocalDiscretization/LinearTri.cpp
@@ -0,0 +1,145 @@
+#include "moab/LinearTri.hpp"
+#include "moab/Forward.hpp"
+#include <algorithm>
+
+namespace moab
+{
+
+ const double LinearTri::corner[3][2] = { {0,0},
+ {1,0},
+ {0,1}};
+
+ ErrorCode LinearTri::initFcn(const double *verts, const int /*nverts*/, double *&work) {
+ // allocate work array as:
+ // work[0..8] = T
+ // work[9..17] = Tinv
+ // work[18] = detT
+ // work[19] = detTinv
+ assert(!work && verts);
+ work = new double[20];
+ Matrix3 *T = reinterpret_cast<Matrix3*>(work),
+ *Tinv = reinterpret_cast<Matrix3*>(work+9);
+ double *detT = work+18, *detTinv = work+19;
+
+ *T = Matrix3(verts[1*3+0]-verts[0*3+0],verts[2*3+0]-verts[0*3+0],0.0,
+ verts[1*3+1]-verts[0*3+1],verts[2*3+1]-verts[0*3+1],0.0,
+ verts[1*3+2]-verts[0*3+2],verts[2*3+2]-verts[0*3+2],1.0);
+ *T *= 0.5;
+ (*T)(2,2) = 1.0;
+
+ *Tinv = T->inverse();
+ *detT = T->determinant();
+ *detTinv = (0.0 == *detT ? HUGE : 1.0 / *detT);
+
+ return MB_SUCCESS;
+ }
+
+ ErrorCode LinearTri::evalFcn(const double *params, const double *field, const int /*ndim*/, const int num_tuples,
+ double */*work*/, double *result) {
+ assert(params && field && num_tuples > 0);
+ // convert to [0,1]
+ double p1 = 0.5 * (1.0 + params[0]),
+ p2 = 0.5 * (1.0 + params[1]),
+ p0 = 1.0 - p1 - p2;
+
+ for (int j = 0; j < num_tuples; j++)
+ result[j] = p0 * field[0*num_tuples+j] + p1 * field[1*num_tuples+j] + p2 * field[2*num_tuples+j];
+
+ return MB_SUCCESS;
+ }
+
+ ErrorCode LinearTri::integrateFcn(const double *field, const double */*verts*/, const int /*nverts*/, const int /*ndim*/, const int num_tuples,
+ double *work, double *result)
+ {
+ assert(field && num_tuples > 0);
+ double tmp = work[18];
+
+ for (int i = 0; i < num_tuples; i++)
+ result[i] = tmp * (field[num_tuples+i] + field[2*num_tuples+i]);
+
+ return MB_SUCCESS;
+ }
+
+ ErrorCode LinearTri::jacobianFcn(const double *, const double *, const int, const int ,
+ double *work, double *result)
+ {
+ // jacobian is cached in work array
+ assert(work);
+ std::copy(work, work+9, result);
+ return MB_SUCCESS;
+ }
+
+ ErrorCode LinearTri::reverseEvalFcn(EvalFcn eval, JacobianFcn jacob, InsideFcn ins,
+ const double *posn, const double *verts, const int nverts, const int ndim,
+ const double iter_tol, const double inside_tol, double *work,
+ double *params, bool *is_inside)
+ {
+ assert(posn && verts);
+ return evaluate_reverse(eval, jacob, ins, posn, verts, nverts, ndim, iter_tol, inside_tol, work,
+ params, is_inside);
+ }
+
+ bool LinearTri::insideFcn(const double *params, const int , const double tol)
+ {
+ return (params[0] >= -1.0-tol && params[1] >= -1.0-tol &&
+ params[0] + params[1] <= 1.0+tol);
+
+ }
+
+ ErrorCode LinearTri::evaluate_reverse(EvalFcn eval, JacobianFcn jacob, InsideFcn inside_f,
+ const double *posn, const double *verts, const int nverts,
+ const int ndim, const double iter_tol, const double inside_tol,
+ double *work, double *params, bool *inside) {
+ // TODO: should differentiate between epsilons used for
+ // Newton Raphson iteration, and epsilons used for curved boundary geometry errors
+ // right now, fix the tolerance used for NR
+ const double error_tol_sqr = iter_tol*iter_tol;
+ CartVect *cvparams = reinterpret_cast<CartVect*>(params);
+ const CartVect *cvposn = reinterpret_cast<const CartVect*>(posn);
+
+ // find best initial guess to improve convergence
+ CartVect tmp_params[] = {CartVect(-1,-1,-1), CartVect(1,-1,-1), CartVect(-1,1,-1)};
+ double resl = HUGE;
+ CartVect new_pos, tmp_pos;
+ ErrorCode rval;
+ for (unsigned int i = 0; i < 3; i++) {
+ rval = (*eval)(tmp_params[i].array(), verts, ndim, 3, work, tmp_pos.array());
+ if (MB_SUCCESS != rval) return rval;
+ double tmp_resl = (tmp_pos-*cvposn).length_squared();
+ if (tmp_resl < resl) {
+ *cvparams = tmp_params[i];
+ new_pos = tmp_pos;
+ resl = tmp_resl;
+ }
+ }
+
+ // residual is diff between old and new pos; need to minimize that
+ CartVect res = new_pos - *cvposn;
+ Matrix3 J;
+ rval = (*jacob)(cvparams->array(), verts, nverts, ndim, work, J[0]);
+ double det = J.determinant();
+ assert(det > std::numeric_limits<double>::epsilon());
+ Matrix3 Ji = J.inverse(1.0/det);
+
+ int iters=0;
+ // while |res| larger than tol
+ while (res % res > error_tol_sqr) {
+ if(++iters>25)
+ return MB_FAILURE;
+
+ // new params tries to eliminate residual
+ *cvparams -= Ji * res;
+
+ // get the new forward-evaluated position, and its difference from the target pt
+ rval = (*eval)(params, verts, ndim, 3, work, new_pos.array());
+ if (MB_SUCCESS != rval) return rval;
+ res = new_pos - *cvposn;
+ }
+
+ if (inside)
+ *inside = (*inside_f)(params, ndim, inside_tol);
+
+ return MB_SUCCESS;
+ }// Map::evaluate_reverse()
+
+} // namespace moab
diff --git a/src/LocalDiscretization/Makefile.am b/src/LocalDiscretization/Makefile.am
index 7309aef..cb34071 100644
--- a/src/LocalDiscretization/Makefile.am
+++ b/src/LocalDiscretization/Makefile.am
@@ -17,6 +17,7 @@ MOAB_LOCALDISCR_SRCS = \
LinearHex.cpp \
LinearQuad.cpp \
LinearTet.cpp \
+ LinearTri.cpp \
QuadraticHex.cpp
MOAB_LOCALDISCR_HDRS = \
@@ -24,6 +25,7 @@ MOAB_LOCALDISCR_HDRS = \
moab/LinearHex.hpp \
moab/LinearQuad.hpp \
moab/LinearTet.hpp \
+ moab/LinearTri.hpp \
moab/QuadraticHex.hpp
# The list of source files, and any header files that do not need to be installed
diff --git a/src/LocalDiscretization/QuadraticHex.cpp b/src/LocalDiscretization/QuadraticHex.cpp
index ca4758b..5c621f6 100644
--- a/src/LocalDiscretization/QuadraticHex.cpp
+++ b/src/LocalDiscretization/QuadraticHex.cpp
@@ -107,10 +107,12 @@ namespace moab
ErrorCode QuadraticHex::reverseEvalFcn(EvalFcn eval, JacobianFcn jacob, InsideFcn ins,
const double *posn, const double *verts, const int nverts, const int ndim,
- const double tol, double *work, double *params, bool *is_inside)
+ const double iter_tol, const double inside_tol, double *work,
+ double *params, bool *is_inside)
{
assert(posn && verts);
- return EvalSet::evaluate_reverse(eval, jacob, ins, posn, verts, nverts, ndim, tol, work, params, is_inside);
+ return EvalSet::evaluate_reverse(eval, jacob, ins, posn, verts, nverts, ndim, iter_tol, inside_tol,
+ work, params, is_inside);
}
bool QuadraticHex::insideFcn(const double *params, const int ndim, const double tol)
diff --git a/src/LocalDiscretization/SpectralHex.cpp b/src/LocalDiscretization/SpectralHex.cpp
index 0e40eba..31f5391 100644
--- a/src/LocalDiscretization/SpectralHex.cpp
+++ b/src/LocalDiscretization/SpectralHex.cpp
@@ -73,7 +73,8 @@ CartVect SpectralHex::evaluate( const CartVect& params ) const
return result;
}
// replicate the functionality of hex_findpt
-bool SpectralHex::evaluate_reverse(CartVect const & xyz, CartVect ¶ms, double tol, const CartVect &init) const
+ bool SpectralHex::evaluate_reverse(CartVect const & xyz, CartVect ¶ms, double iter_tol, const double inside_tol,
+ const CartVect &init) const
{
params = init;
@@ -92,7 +93,7 @@ bool SpectralHex::evaluate_reverse(CartVect const & xyz, CartVect ¶ms, doubl
//copy parametric coords back
params = r;
- return is_inside(params, tol);
+ return is_inside(params, inside_tol);
}
Matrix3 SpectralHex::jacobian(const CartVect& params) const
{
diff --git a/src/LocalDiscretization/SpectralQuad.cpp b/src/LocalDiscretization/SpectralQuad.cpp
index cbd1251..5c11d13 100644
--- a/src/LocalDiscretization/SpectralQuad.cpp
+++ b/src/LocalDiscretization/SpectralQuad.cpp
@@ -95,7 +95,8 @@ CartVect SpectralQuad::evalFcn(const double *params, const double *field, const
}
// replicate the functionality of hex_findpt
bool SpectralQuad::reverseEvalFcn(const double *posn, const double *verts, const int nverts, const int ndim,
- const double tol, double *work, double *params, bool *is_inside)
+ const double iter_tol, const double inside_tol, double *work,
+ double *params, bool *is_inside)
{
params = init;
@@ -114,7 +115,7 @@ bool SpectralQuad::reverseEvalFcn(const double *posn, const double *verts, const
//copy parametric coords back
params = r;
- return insideFcn(params, 2, tol);
+ return insideFcn(params, 2, inside_tol);
}
diff --git a/src/LocalDiscretization/moab/ElemEvaluator.hpp b/src/LocalDiscretization/moab/ElemEvaluator.hpp
index 48d34ac..29d1051 100644
--- a/src/LocalDiscretization/moab/ElemEvaluator.hpp
+++ b/src/LocalDiscretization/moab/ElemEvaluator.hpp
@@ -25,7 +25,8 @@ namespace moab {
typedef ErrorCode (*ReverseEvalFcn)(EvalFcn eval, JacobianFcn jacob, InsideFcn ins,
const double *posn, const double *verts, const int nverts, const int ndim,
- const double tol, double *work, double *params, bool *is_inside);
+ const double iter_tol, const double inside_tol,
+ double *work, double *params, bool *is_inside);
class EvalSet
{
@@ -76,8 +77,8 @@ namespace moab {
/** \brief Common function to do reverse evaluation based on evaluation and jacobian functions */
static ErrorCode evaluate_reverse(EvalFcn eval, JacobianFcn jacob, InsideFcn inside_f,
const double *posn, const double *verts, const int nverts,
- const int ndim, const double tol, double *work, double *params,
- bool *inside);
+ const int ndim, const double iter_tol, const double inside_tol,
+ double *work, double *params, bool *inside);
/** \brief Common function that returns true if params is in [-1,1]^ndims */
static bool inside_function(const double *params, const int ndims, const double tol);
};
@@ -119,9 +120,9 @@ namespace moab {
* \param impl MOAB instance
* \param ent Entity handle to cache on the evaluator
* \param tag Tag to cache on the evaluator
- * \param tag_dim Tag dimension to cache on the evaluator
+ * \param tagged_ent_dim Dimension of entities to be tagged to cache on the evaluator
*/
- ElemEvaluator(Interface *impl, EntityHandle ent = 0, Tag tag = 0, int tag_dim = -1);
+ ElemEvaluator(Interface *impl, EntityHandle ent = 0, Tag tag = 0, int tagged_ent_dim = -1);
/** \brief Evaluate cached tag at a given parametric location within the cached entity
* If evaluating coordinates, call set_tag(0, 0), which indicates coords instead of a tag.
@@ -133,12 +134,14 @@ namespace moab {
/** \brief Reverse-evaluate the cached entity at a given physical position
* \param posn Position at which to evaluate parameters
- * \param tol Tolerance of reverse evaluation, usually 10^-6 or so
+ * \param iter_tol Tolerance of reverse evaluation non-linear iteration, usually 10^-10 or so
+ * \param inside_tol Tolerance of is_inside evaluation, usually 10^-6 or so
* \param params Result of evaluation
* \param is_inside If non-NULL, returns true of resulting parameters place the point inside the element
* (in most cases, within [-1]*(dim)
*/
- ErrorCode reverse_eval(const double *posn, double tol, double *params, bool *is_inside = NULL) const;
+ ErrorCode reverse_eval(const double *posn, double iter_tol, double inside_tol, double *params,
+ bool *is_inside = NULL) const;
/** \brief Evaluate the jacobian of the cached entity at a given parametric location
* \param params Parameters at which to evaluate jacobian
@@ -166,14 +169,16 @@ namespace moab {
* object is changed.
* \param entities Entities tested
* \param point Point tested, must have 3 dimensions, even for edge and face entities
- * \param tol Tolerance for is_inside test
+ * \param iter_tol Tolerance for non-linear reverse evaluation
+ * \param inside_tol Tolerance for is_inside test
* \param containing_ent Entity containing the point, returned 0 if no entity
* \param params Parameters of point in containing entity, unchanged if no containing entity
* \param num_evals If non-NULL, incremented each time reverse_eval is called
* \return Returns non-success only if evaulation failed for some reason (point not in element is NOT a
* reason for failure)
*/
- ErrorCode find_containing_entity(Range &entities, const double *point, double tol,
+ ErrorCode find_containing_entity(Range &entities, const double *point,
+ const double iter_tol, const double inside_tol,
EntityHandle &containing_ent, double *params,
unsigned int *num_evals = NULL);
@@ -186,14 +191,16 @@ namespace moab {
* object is changed.
* \param ent_set Entity set containing the entities to be tested
* \param point Point tested, must have 3 dimensions, even for edge and face entities
- * \param tol Tolerance for is_inside test
+ * \param iter_tol Tolerance for non-linear reverse evaluation
+ * \param inside_tol Tolerance for is_inside test
* \param containing_ent Entity containing the point, returned 0 if no entity
* \param params Parameters of point in containing entity, unchanged if no containing entity
* \param num_evals If non-NULL, incremented each time reverse_eval is called
* \return Returns non-success only if evaulation failed for some reason (point not in element is NOT a
* reason for failure)
*/
- ErrorCode find_containing_entity(EntityHandle ent_set, const double *point, double tol,
+ ErrorCode find_containing_entity(EntityHandle ent_set, const double *point,
+ const double iter_tol, const double inside_tol,
EntityHandle &containing_ent, double *params,
unsigned int *num_evals = NULL);
@@ -218,8 +225,12 @@ namespace moab {
inline ErrorCode set_ent_handle(EntityHandle ent);
/** \brief Get entity handle for this ElemEval */
- inline EntityHandle get_ent_handle() const {return entHandle;};
+ inline EntityHandle get_ent_handle() const {return entHandle;}
+ /* \brief Get vertex positions cached on this EE
+ */
+ inline double *get_vert_pos() {return vertPos[0].array();}
+
/* \brief Get the vertex handles cached here */
inline const EntityHandle *get_vert_handles() const {return vertHandles;}
@@ -233,24 +244,29 @@ namespace moab {
* To designate that vertex coordinates are the desired tag, pass in a tag handle of 0
* and a tag dimension of 0.
* \param tag Tag handle to cache, or 0 to cache vertex positions
- * \param tag_dim Dimension of entities tagged with this tag
+ * \param tagged_ent_dim Dimension of entities tagged with this tag
*/
- inline ErrorCode set_tag_handle(Tag tag, int tag_dim = -1);
+ inline ErrorCode set_tag_handle(Tag tag, int tagged_ent_dim = -1);
/* \brief Set the name of the tag to cache on this evaluator
* To designate that vertex coordinates are the desired tag, pass in "COORDS" as the name
* and a tag dimension of 0.
- * \param tag Tag handle to cache, or 0 to cache vertex positions
- * \param tag_dim Dimension of entities tagged with this tag
+ * \param tag_name Tag handle to cache, or 0 to cache vertex positions
+ * \param tagged_ent_dim Dimension of entities tagged with this tag
*/
- inline ErrorCode set_tag(const char *tag_name, int ent_dim = -1);
+ inline ErrorCode set_tag(const char *tag_name, int tagged_ent_dim = -1);
/* \brief Get the dimension of the entities on which tag is set */
- inline int get_tag_dim() const {return tagDim;};
+ inline int get_tagged_ent_dim() const {return taggedEntDim;};
/* \brief Set the dimension of entities having the tag */
- inline ErrorCode set_tag_dim(int dim);
+ inline ErrorCode set_tagged_ent_dim(int dim);
+ /* \brief Get work space, sometimes this is useful for evaluating data you don't want to set as tags on entities
+ * Can't be const because most of the functions (evaluate, integrate, etc.) take non-const work space *'s
+ */
+ inline double *get_work_space() {return workSpace;}
+
/* \brief MOAB interface cached on this evaluator */
inline Interface *get_moab() {return mbImpl;}
@@ -287,7 +303,7 @@ namespace moab {
int numTuples;
/** \brief Dimension of entities from which to grab tag */
- int tagDim;
+ int taggedEntDim;
/** \brief Tag space */
std::vector<unsigned char> tagSpace;
@@ -300,13 +316,13 @@ namespace moab {
}; // class ElemEvaluator
- inline ElemEvaluator::ElemEvaluator(Interface *impl, EntityHandle ent, Tag tag, int tag_dim)
+ inline ElemEvaluator::ElemEvaluator(Interface *impl, EntityHandle ent, Tag tag, int tagged_ent_dim)
: mbImpl(impl), entHandle(0), entType(MBMAXTYPE), entDim(-1), numVerts(0),
vertHandles(NULL), tagHandle(0), tagCoords(false), numTuples(0),
- tagDim(0), workSpace(NULL)
+ taggedEntDim(0), workSpace(NULL)
{
if (ent) set_ent_handle(ent);
- if (tag) set_tag_handle(tag, tag_dim);
+ if (tag) set_tag_handle(tag, tagged_ent_dim);
}
inline ErrorCode ElemEvaluator::set_ent_handle(EntityHandle ent)
@@ -329,18 +345,17 @@ namespace moab {
rval = set_tag_handle(tagHandle);
if (MB_SUCCESS != rval) return rval;
}
-
if (evalSets[entType].initFcn) return (*evalSets[entType].initFcn)(vertPos[0].array(), numVerts, workSpace);
return MB_SUCCESS;
}
- inline ErrorCode ElemEvaluator::set_tag_handle(Tag tag, int tag_dim)
+ inline ErrorCode ElemEvaluator::set_tag_handle(Tag tag, int tagged_ent_dim)
{
ErrorCode rval = MB_SUCCESS;
- if (!tag && !tag_dim) {
+ if (!tag && !tagged_ent_dim) {
tagCoords = true;
numTuples = 3;
- tagDim = 0;
+ taggedEntDim = 0;
tagHandle = 0;
return rval;
}
@@ -355,14 +370,14 @@ namespace moab {
tagCoords = false;
}
- tagDim = (-1 == tag_dim ? 0 : tag_dim);
+ taggedEntDim = (-1 == tagged_ent_dim ? 0 : tagged_ent_dim);
if (entHandle) {
- if (0 == tagDim) {
+ if (0 == taggedEntDim) {
rval = mbImpl->tag_get_data(tagHandle, vertHandles, numVerts, &tagSpace[0]);
if (MB_SUCCESS != rval) return rval;
}
- else if (tagDim == entDim) {
+ else if (taggedEntDim == entDim) {
rval = mbImpl->tag_get_data(tagHandle, &entHandle, 1, &tagSpace[0]);
if (MB_SUCCESS != rval) return rval;
}
@@ -371,14 +386,14 @@ namespace moab {
return rval;
}
- inline ErrorCode ElemEvaluator::set_tag(const char *tag_name, int tag_dim)
+ inline ErrorCode ElemEvaluator::set_tag(const char *tag_name, int tagged_ent_dim)
{
ErrorCode rval = MB_SUCCESS;
if (!tag_name || strlen(tag_name) == 0) return MB_FAILURE;
Tag tag;
if (!strcmp(tag_name, "COORDS")) {
tagCoords = true;
- tagDim = 0;
+ taggedEntDim = 0;
numTuples = 3;
tagHandle = 0;
// can return here, because vertex coords already cached when entity handle set
@@ -399,15 +414,15 @@ namespace moab {
tagCoords = false;
}
- tagDim = (-1 == tag_dim ? entDim : tag_dim);
+ taggedEntDim = (-1 == tagged_ent_dim ? entDim : tagged_ent_dim);
}
if (entHandle) {
- if (0 == tagDim) {
+ if (0 == taggedEntDim) {
rval = mbImpl->tag_get_data(tagHandle, vertHandles, numVerts, &tagSpace[0]);
if (MB_SUCCESS != rval) return rval;
}
- else if (tagDim == entDim) {
+ else if (taggedEntDim == entDim) {
rval = mbImpl->tag_get_data(tagHandle, &entHandle, 1, &tagSpace[0]);
if (MB_SUCCESS != rval) return rval;
}
@@ -431,16 +446,17 @@ namespace moab {
assert(entHandle && MBMAXTYPE != entType);
return (*evalSets[entType].evalFcn)(params,
(tagCoords ? (const double*) vertPos[0].array(): (const double*)&tagSpace[0]),
- entDim, (-1 == num_tuples ? numTuples : num_tuples),
- workSpace, result);
+ entDim,
+ (-1 == num_tuples ? numTuples : num_tuples), workSpace, result);
}
- inline ErrorCode ElemEvaluator::reverse_eval(const double *posn, const double tol, double *params, bool *ins) const
+ inline ErrorCode ElemEvaluator::reverse_eval(const double *posn, const double iter_tol, const double inside_tol,
+ double *params, bool *ins) const
{
assert(entHandle && MBMAXTYPE != entType);
return (*evalSets[entType].reverseEvalFcn)(evalSets[entType].evalFcn, evalSets[entType].jacobianFcn, evalSets[entType].insideFcn,
- posn, vertPos[0].array(), numVerts, entDim, tol, workSpace,
- params, ins);
+ posn, vertPos[0].array(), numVerts,
+ entDim, iter_tol, inside_tol, workSpace, params, ins);
}
/** \brief Evaluate the jacobian of the cached entity at a given parametric location */
@@ -456,7 +472,7 @@ namespace moab {
assert(entHandle && MBMAXTYPE != entType && (tagCoords || tagHandle));
ErrorCode rval = MB_SUCCESS;
if (!tagCoords) {
- if (0 == tagDim) rval = mbImpl->tag_get_data(tagHandle, vertHandles, numVerts, (void*)&tagSpace[0]);
+ if (0 == taggedEntDim) rval = mbImpl->tag_get_data(tagHandle, vertHandles, numVerts, (void*)&tagSpace[0]);
else rval = mbImpl->tag_get_data(tagHandle, &entHandle, 1, (void*)&tagSpace[0]);
if (MB_SUCCESS != rval) return rval;
}
@@ -465,7 +481,8 @@ namespace moab {
workSpace, result);
}
- inline ErrorCode ElemEvaluator::find_containing_entity(EntityHandle ent_set, const double *point, double tol,
+ inline ErrorCode ElemEvaluator::find_containing_entity(EntityHandle ent_set, const double *point,
+ const double iter_tol, const double inside_tol,
EntityHandle &containing_ent, double *params,
unsigned int *num_evals)
{
@@ -473,7 +490,7 @@ namespace moab {
Range entities;
ErrorCode rval = mbImpl->get_entities_by_handle(ent_set, entities);
if (MB_SUCCESS != rval) return rval;
- else return find_containing_entity(entities, point, tol, containing_ent, params, num_evals);
+ else return find_containing_entity(entities, point, iter_tol, inside_tol, containing_ent, params, num_evals);
}
inline bool ElemEvaluator::inside(const double *params, const double tol) const
diff --git a/src/LocalDiscretization/moab/LinearHex.hpp b/src/LocalDiscretization/moab/LinearHex.hpp
index e7ca3f4..03ed33f 100644
--- a/src/LocalDiscretization/moab/LinearHex.hpp
+++ b/src/LocalDiscretization/moab/LinearHex.hpp
@@ -17,7 +17,8 @@ public:
/** \brief Reverse-evaluation of parametric coordinates at physical space position */
static ErrorCode reverseEvalFcn(EvalFcn eval, JacobianFcn jacob, InsideFcn ins,
const double *posn, const double *verts, const int nverts, const int ndim,
- const double tol, double *work, double *params, bool *is_inside);
+ const double iter_tol, const double inside_tol, double *work,
+ double *params, bool *is_inside);
/** \brief Evaluate the jacobian at a specified parametric position */
static ErrorCode jacobianFcn(const double *params, const double *verts, const int nverts, const int ndim,
This diff is so big that we needed to truncate the remainder.
https://bitbucket.org/fathomteam/moab/commits/f4d6d80729a8/
Changeset: f4d6d80729a8
Branch: None
User: pshriwise
Date: 2014-01-07 21:24:32
Summary: Added a ReadCGM test for a simple cube file.
Affected #: 2 files
diff --git a/MeshFiles/unittest/io/cube.sat b/MeshFiles/unittest/io/cube.sat
new file mode 100644
index 0000000..58fd508
--- /dev/null
+++ b/MeshFiles/unittest/io/cube.sat
@@ -0,0 +1,117 @@
+1900 0 1 0
+10 Cubit 12.2 17 ACIS 19.0.2 Linux 24 Mon Jan 6 14:16:03 2014
+1 9.9999999999999995e-07 1e-10
+body $1 -1 -1 $-1 $2 $-1 $-1 T -5 -5 -5 5 5 5 #
+integer_attrib-name_attrib-gen-attrib $-1 -1 $-1 $-1 $0 2 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 @8 CUBIT_ID 1 #
+lump $3 -1 -1 $-1 $-1 $4 $0 T -5 -5 -5 5 5 5 #
+integer_attrib-name_attrib-gen-attrib $-1 -1 $-1 $-1 $2 2 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 @8 CUBIT_ID 1 #
+shell $-1 -1 -1 $-1 $-1 $-1 $5 $-1 $2 T -5 -5 -5 5 5 5 #
+face $6 -1 -1 $-1 $7 $8 $4 $-1 $9 forward single T -5 -5 5 5 5 5 F #
+integer_attrib-name_attrib-gen-attrib $-1 -1 $-1 $-1 $5 2 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 @8 CUBIT_ID 1 #
+face $10 -1 -1 $-1 $11 $12 $4 $-1 $13 reversed single T -5 -5 -5 5 5 -5 F #
+loop $-1 -1 -1 $-1 $-1 $14 $5 T -5 -5 5 5 5 5 unknown #
+plane-surface $-1 -1 -1 $-1 0 0 5 0 0 1 1 0 0 forward_v I I I I #
+integer_attrib-name_attrib-gen-attrib $-1 -1 $-1 $-1 $7 2 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 @8 CUBIT_ID 2 #
+face $15 -1 -1 $-1 $16 $17 $4 $-1 $18 reversed single T -5 -5 -5 5 -5 5 F #
+loop $-1 -1 -1 $-1 $-1 $19 $7 T -5 -5 -5 5 5 -5 unknown #
+plane-surface $-1 -1 -1 $-1 0 0 -5 0 0 1 1 0 0 forward_v I I I I #
+coedge $-1 -1 -1 $-1 $20 $21 $22 $23 forward $8 $-1 #
+integer_attrib-name_attrib-gen-attrib $-1 -1 $-1 $-1 $11 2 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 @8 CUBIT_ID 3 #
+face $24 -1 -1 $-1 $25 $26 $4 $-1 $27 reversed single T -5 -5 -5 -5 5 5 F #
+loop $-1 -1 -1 $-1 $-1 $28 $11 T -5 -5 -5 5 -5 5 unknown #
+plane-surface $-1 -1 -1 $-1 0 -5 0 0 1 -0 -0 0 1 forward_v I I I I #
+coedge $-1 -1 -1 $-1 $29 $30 $31 $32 forward $12 $-1 #
+coedge $-1 -1 -1 $-1 $33 $14 $34 $35 forward $8 $-1 #
+coedge $-1 -1 -1 $-1 $14 $33 $36 $37 forward $8 $-1 #
+coedge $-1 -1 -1 $-1 $38 $39 $14 $23 reversed $40 $-1 #
+edge $41 -1 -1 $-1 $42 -5 $43 5 $22 $44 forward @7 unknown T 5 -5 5 5 5 5 #
+integer_attrib-name_attrib-gen-attrib $-1 -1 $-1 $-1 $16 2 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 @8 CUBIT_ID 4 #
+face $45 -1 -1 $-1 $46 $47 $4 $-1 $48 reversed single T -5 5 -5 5 5 5 F #
+loop $-1 -1 -1 $-1 $-1 $49 $16 T -5 -5 -5 -5 5 5 unknown #
+plane-surface $-1 -1 -1 $-1 -5 0 0 1 0 0 0 0 -1 forward_v I I I I #
+coedge $-1 -1 -1 $-1 $50 $36 $51 $52 forward $17 $-1 #
+coedge $-1 -1 -1 $-1 $53 $19 $50 $54 forward $12 $-1 #
+coedge $-1 -1 -1 $-1 $19 $53 $55 $56 forward $12 $-1 #
+coedge $-1 -1 -1 $-1 $39 $38 $19 $32 reversed $40 $-1 #
+edge $57 -1 -1 $-1 $58 -5 $59 5 $31 $60 forward @7 unknown T 5 -5 -5 5 5 -5 #
+coedge $-1 -1 -1 $-1 $21 $20 $61 $62 forward $8 $-1 #
+coedge $-1 -1 -1 $-1 $63 $64 $20 $35 reversed $47 $-1 #
+edge $65 -1 -1 $-1 $43 -5 $66 5 $34 $67 forward @7 unknown T -5 5 5 5 5 5 #
+coedge $-1 -1 -1 $-1 $28 $68 $21 $37 reversed $17 $-1 #
+edge $69 -1 -1 $-1 $70 -5 $42 5 $36 $71 forward @7 unknown T -5 -5 5 5 -5 5 #
+coedge $-1 -1 -1 $-1 $31 $22 $68 $72 forward $40 $-1 #
+coedge $-1 -1 -1 $-1 $22 $31 $63 $73 reversed $40 $-1 #
+loop $-1 -1 -1 $-1 $-1 $38 $46 T 5 -5 -5 5 5 5 unknown #
+integer_attrib-name_attrib-gen-attrib $-1 -1 $-1 $-1 $23 2 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 @8 CUBIT_ID 1 #
+vertex $74 -1 -1 $-1 $23 $75 #
+vertex $76 -1 -1 $-1 $23 $77 #
+straight-curve $-1 -1 -1 $-1 5 0 5 0 1 0 I I #
+integer_attrib-name_attrib-gen-attrib $-1 -1 $-1 $-1 $25 2 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 @8 CUBIT_ID 5 #
+face $78 -1 -1 $-1 $-1 $40 $4 $-1 $79 reversed single T 5 -5 -5 5 5 5 F #
+loop $-1 -1 -1 $-1 $-1 $63 $25 T -5 5 -5 5 5 5 unknown #
+plane-surface $-1 -1 -1 $-1 0 5 0 0 -1 0 0 0 -1 forward_v I I I I #
+coedge $-1 -1 -1 $-1 $80 $61 $64 $81 forward $26 $-1 #
+coedge $-1 -1 -1 $-1 $68 $28 $29 $54 reversed $17 $-1 #
+coedge $-1 -1 -1 $-1 $61 $80 $28 $52 reversed $26 $-1 #
+edge $82 -1 -1 $-1 $70 -5 $83 5 $51 $84 forward @7 unknown T -5 -5 -5 -5 -5 5 #
+coedge $-1 -1 -1 $-1 $30 $29 $80 $85 forward $12 $-1 #
+edge $86 -1 -1 $-1 $59 -5 $83 5 $50 $87 forward @7 unknown T -5 -5 -5 5 -5 -5 #
+coedge $-1 -1 -1 $-1 $64 $63 $30 $56 reversed $47 $-1 #
+edge $88 -1 -1 $-1 $89 -5 $58 5 $55 $90 forward @7 unknown T -5 5 -5 5 5 -5 #
+integer_attrib-name_attrib-gen-attrib $-1 -1 $-1 $-1 $32 2 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 @8 CUBIT_ID 5 #
+vertex $91 -1 -1 $-1 $32 $92 #
+vertex $93 -1 -1 $-1 $72 $94 #
+straight-curve $-1 -1 -1 $-1 5 0 -5 0 -1 0 I I #
+coedge $-1 -1 -1 $-1 $49 $51 $33 $62 reversed $26 $-1 #
+edge $95 -1 -1 $-1 $66 -5 $70 5 $61 $96 forward @7 unknown T -5 -5 5 -5 5 5 #
+coedge $-1 -1 -1 $-1 $55 $34 $39 $73 forward $47 $-1 #
+coedge $-1 -1 -1 $-1 $34 $55 $49 $81 reversed $47 $-1 #
+integer_attrib-name_attrib-gen-attrib $-1 -1 $-1 $-1 $35 2 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 @8 CUBIT_ID 2 #
+vertex $97 -1 -1 $-1 $35 $98 #
+straight-curve $-1 -1 -1 $-1 0 5 5 -1 0 0 I I #
+coedge $-1 -1 -1 $-1 $36 $50 $38 $72 reversed $17 $-1 #
+integer_attrib-name_attrib-gen-attrib $-1 -1 $-1 $-1 $37 2 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 @8 CUBIT_ID 4 #
+vertex $99 -1 -1 $-1 $62 $100 #
+straight-curve $-1 -1 -1 $-1 0 -5 5 1 0 0 I I #
+edge $101 -1 -1 $-1 $42 -5 $59 5 $38 $102 forward @7 unknown T 5 -5 -5 5 -5 5 #
+edge $103 -1 -1 $-1 $43 -5 $58 5 $39 $104 forward @7 unknown T 5 5 -5 5 5 5 #
+integer_attrib-name_attrib-gen-attrib $-1 -1 $-1 $-1 $42 2 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 @8 CUBIT_ID 1 #
+point $-1 -1 -1 $-1 5 -5 5 #
+integer_attrib-name_attrib-gen-attrib $-1 -1 $-1 $-1 $43 2 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 @8 CUBIT_ID 2 #
+point $-1 -1 -1 $-1 5 5 5 #
+integer_attrib-name_attrib-gen-attrib $-1 -1 $-1 $-1 $46 2 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 @8 CUBIT_ID 6 #
+plane-surface $-1 -1 -1 $-1 5 0 0 -1 0 0 0 -0 1 forward_v I I I I #
+coedge $-1 -1 -1 $-1 $51 $49 $53 $85 reversed $26 $-1 #
+edge $105 -1 -1 $-1 $66 -5 $89 5 $64 $106 forward @7 unknown T -5 5 -5 -5 5 5 #
+integer_attrib-name_attrib-gen-attrib $-1 -1 $-1 $-1 $52 2 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 @8 CUBIT_ID 9 #
+vertex $107 -1 -1 $-1 $85 $108 #
+straight-curve $-1 -1 -1 $-1 -5 -5 0 0 0 -1 I I #
+edge $109 -1 -1 $-1 $83 -5 $89 5 $80 $110 forward @7 unknown T -5 -5 -5 -5 5 -5 #
+integer_attrib-name_attrib-gen-attrib $-1 -1 $-1 $-1 $54 2 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 @8 CUBIT_ID 6 #
+straight-curve $-1 -1 -1 $-1 0 -5 -5 -1 0 0 I I #
+integer_attrib-name_attrib-gen-attrib $-1 -1 $-1 $-1 $56 2 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 @8 CUBIT_ID 8 #
+vertex $111 -1 -1 $-1 $56 $112 #
+straight-curve $-1 -1 -1 $-1 0 5 -5 1 0 0 I I #
+integer_attrib-name_attrib-gen-attrib $-1 -1 $-1 $-1 $58 2 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 @8 CUBIT_ID 5 #
+point $-1 -1 -1 $-1 5 5 -5 #
+integer_attrib-name_attrib-gen-attrib $-1 -1 $-1 $-1 $59 2 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 @8 CUBIT_ID 6 #
+point $-1 -1 -1 $-1 5 -5 -5 #
+integer_attrib-name_attrib-gen-attrib $-1 -1 $-1 $-1 $62 2 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 @8 CUBIT_ID 3 #
+straight-curve $-1 -1 -1 $-1 -5 0 5 0 -1 0 I I #
+integer_attrib-name_attrib-gen-attrib $-1 -1 $-1 $-1 $66 2 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 @8 CUBIT_ID 3 #
+point $-1 -1 -1 $-1 -5 5 5 #
+integer_attrib-name_attrib-gen-attrib $-1 -1 $-1 $-1 $70 2 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 @8 CUBIT_ID 4 #
+point $-1 -1 -1 $-1 -5 -5 5 #
+integer_attrib-name_attrib-gen-attrib $-1 -1 $-1 $-1 $72 2 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 @8 CUBIT_ID 10 #
+straight-curve $-1 -1 -1 $-1 5 -5 0 0 0 -1 I I #
+integer_attrib-name_attrib-gen-attrib $-1 -1 $-1 $-1 $73 2 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 @8 CUBIT_ID 12 #
+straight-curve $-1 -1 -1 $-1 5 5 0 0 0 -1 I I #
+integer_attrib-name_attrib-gen-attrib $-1 -1 $-1 $-1 $81 2 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 @8 CUBIT_ID 11 #
+straight-curve $-1 -1 -1 $-1 -5 5 0 0 0 -1 I I #
+integer_attrib-name_attrib-gen-attrib $-1 -1 $-1 $-1 $83 2 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 @8 CUBIT_ID 7 #
+point $-1 -1 -1 $-1 -5 -5 -5 #
+integer_attrib-name_attrib-gen-attrib $-1 -1 $-1 $-1 $85 2 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 @8 CUBIT_ID 7 #
+straight-curve $-1 -1 -1 $-1 -5 0 -5 0 1 0 I I #
+integer_attrib-name_attrib-gen-attrib $-1 -1 $-1 $-1 $89 2 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 @8 CUBIT_ID 8 #
+point $-1 -1 -1 $-1 -5 5 -5 #
+End-of-ACIS-data
\ No newline at end of file
diff --git a/test/io/read_cgm_test.cpp b/test/io/read_cgm_test.cpp
index 5377220..8891f8c 100644
--- a/test/io/read_cgm_test.cpp
+++ b/test/io/read_cgm_test.cpp
@@ -7,6 +7,7 @@
#include "TestUtil.hpp"
#include "Internals.hpp"
#include "moab/Core.hpp"
+#include "MBTagConventions.hpp"
using namespace moab;
@@ -29,6 +30,7 @@ static const char input_file[] = "dum.sat";
#endif
#endif
+static const char input_cube[] = STRINGIFY(MESHDIR) "/io/cube.sat";
void read_multiple_test()
{
Core mb;
@@ -38,12 +40,37 @@ void read_multiple_test()
// second load
rval = mb.load_file(input_file);
CHECK_ERR(rval);
+}
+
+void read_cube_test()
+{
+ Core mb;
+
+ ErrorCode rval = mb.load_file(input_cube); CHECK_ERR(rval);
+
+ Tag geom_tag;
+
+ rval = mb.tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER,
+ geom_tag, MB_TAG_SPARSE|MB_TAG_CREAT);
+ CHECK_ERR(rval);
+
+ int number_of_tris;
+
+ rval = mb.get_number_entities_by_type(0, MBTRI , number_of_tris);
+ std::cout << "Number of Triangles = " << number_of_tris << std::endl;
+ CHECK_ERR(rval);
+
+ if( number_of_tris != 12) rval = MB_FAILURE; CHECK_ERR(rval);
+
+
}
int main(int /* argc */, char** /* argv */)
{
- int result = RUN_TEST( read_multiple_test );
+ int result = RUN_TEST( read_cube_test );
+
+ result = RUN_TEST( read_multiple_test ) ;
return result;
}
https://bitbucket.org/fathomteam/moab/commits/aadd4afe04d5/
Changeset: aadd4afe04d5
Branch: None
User: pshriwise
Date: 2014-02-09 17:14:15
Summary: Added a new test file for ReadCGM: read_cgm_basic_test.cpp
Affected #: 2 files
diff --git a/test/io/read_cgm_basic_test.cpp b/test/io/read_cgm_basic_test.cpp
new file mode 100644
index 0000000..3998d01
--- /dev/null
+++ b/test/io/read_cgm_basic_test.cpp
@@ -0,0 +1,69 @@
+
+#include <iostream>
+#include "moab/Interface.hpp"
+#ifndef IS_BUILDING_MB
+#define IS_BUILDING_MB
+#endif
+#include "TestUtil.hpp"
+#include "Internals.hpp"
+#include "moab/Core.hpp"
+#include "MBTagConventions.hpp"
+
+using namespace moab;
+
+#define CHKERR(A) do { if (MB_SUCCESS != (A)) { \
+ std::cerr << "Failure (error code " << (A) << ") at " __FILE__ ":" \
+ << __LINE__ << std::endl; \
+ return A; } } while(false)
+
+#ifdef MESHDIR
+#ifdef HAVE_OCC_STEP
+static const char input_file[] = STRINGIFY(MESHDIR) "/io/dum.stp";
+#else
+static const char input_file[] = STRINGIFY(MESHDIR) "/io/dum.sat";
+#endif
+#else
+#ifdef HAVE_OCC_STEP
+static const char input_file[] = "dum.stp";
+#else
+static const char input_file[] = "dum.sat";
+#endif
+#endif
+
+#ifdef MESHDIR
+static const char input_cube[] = STRINGIFY(MESHDIR) "/io/cube.sat";
+#else
+static const char input_cube[] = "/io/cube.sat";
+#endif
+
+
+void read_cube_test()
+{
+ Core mb;
+
+ ErrorCode rval = mb.load_file(input_cube); CHECK_ERR(rval);
+
+ Tag geom_tag;
+
+ rval = mb.tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER,
+ geom_tag, MB_TAG_SPARSE|MB_TAG_CREAT);
+ CHECK_ERR(rval);
+
+ int number_of_tris;
+
+ rval = mb.get_number_entities_by_type(0, MBTRI , number_of_tris);
+ std::cout << "Number of Triangles = " << number_of_tris << std::endl;
+ CHECK_ERR(rval);
+
+ if( number_of_tris != 12) rval = MB_FAILURE; CHECK_ERR(rval);
+
+}
+
+int main(int /* argc */, char** /* argv */)
+{
+ int result = 0;
+
+ result += RUN_TEST( read_cube_test );
+
+ return result;
+}
diff --git a/test/io/read_cgm_test.cpp b/test/io/read_cgm_test.cpp
index 8891f8c..5377220 100644
--- a/test/io/read_cgm_test.cpp
+++ b/test/io/read_cgm_test.cpp
@@ -7,7 +7,6 @@
#include "TestUtil.hpp"
#include "Internals.hpp"
#include "moab/Core.hpp"
-#include "MBTagConventions.hpp"
using namespace moab;
@@ -30,7 +29,6 @@ static const char input_file[] = "dum.sat";
#endif
#endif
-static const char input_cube[] = STRINGIFY(MESHDIR) "/io/cube.sat";
void read_multiple_test()
{
Core mb;
@@ -40,37 +38,12 @@ void read_multiple_test()
// second load
rval = mb.load_file(input_file);
CHECK_ERR(rval);
-}
-
-void read_cube_test()
-{
- Core mb;
-
- ErrorCode rval = mb.load_file(input_cube); CHECK_ERR(rval);
-
- Tag geom_tag;
-
- rval = mb.tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER,
- geom_tag, MB_TAG_SPARSE|MB_TAG_CREAT);
- CHECK_ERR(rval);
-
- int number_of_tris;
-
- rval = mb.get_number_entities_by_type(0, MBTRI , number_of_tris);
- std::cout << "Number of Triangles = " << number_of_tris << std::endl;
- CHECK_ERR(rval);
-
- if( number_of_tris != 12) rval = MB_FAILURE; CHECK_ERR(rval);
-
-
}
int main(int /* argc */, char** /* argv */)
{
- int result = RUN_TEST( read_cube_test );
-
- result = RUN_TEST( read_multiple_test ) ;
+ int result = RUN_TEST( read_multiple_test );
return result;
}
https://bitbucket.org/fathomteam/moab/commits/a7026c1cd865/
Changeset: a7026c1cd865
Branch: None
User: pshriwise
Date: 2014-02-09 17:15:44
Summary: Renamed old ReadCGM test file to reduce confusion. Updated Makefile to reflect this change. Added new test file name to the Makefile.
Affected #: 3 files
diff --git a/test/io/Makefile.am b/test/io/Makefile.am
index c776727..e467e78 100644
--- a/test/io/Makefile.am
+++ b/test/io/Makefile.am
@@ -25,7 +25,7 @@ else
endif
if HAVE_CGM
- CGM_TEST = read_cgm_test
+ CGM_TEST = read_cgm_load_test read_cgm_basic_test
else
CGM_TEST =
endif
diff --git a/test/io/read_cgm_load_test.cpp b/test/io/read_cgm_load_test.cpp
new file mode 100644
index 0000000..5377220
--- /dev/null
+++ b/test/io/read_cgm_load_test.cpp
@@ -0,0 +1,49 @@
+
+#include <iostream>
+#include "moab/Interface.hpp"
+#ifndef IS_BUILDING_MB
+#define IS_BUILDING_MB
+#endif
+#include "TestUtil.hpp"
+#include "Internals.hpp"
+#include "moab/Core.hpp"
+
+using namespace moab;
+
+#define CHKERR(A) do { if (MB_SUCCESS != (A)) { \
+ std::cerr << "Failure (error code " << (A) << ") at " __FILE__ ":" \
+ << __LINE__ << std::endl; \
+ return A; } } while(false)
+
+#ifdef MESHDIR
+#ifdef HAVE_OCC_STEP
+static const char input_file[] = STRINGIFY(MESHDIR) "/io/dum.stp";
+#else
+static const char input_file[] = STRINGIFY(MESHDIR) "/io/dum.sat";
+#endif
+#else
+#ifdef HAVE_OCC_STEP
+static const char input_file[] = "dum.stp";
+#else
+static const char input_file[] = "dum.sat";
+#endif
+#endif
+
+void read_multiple_test()
+{
+ Core mb;
+
+ ErrorCode rval = mb.load_file(input_file);
+ CHECK_ERR(rval);
+ // second load
+ rval = mb.load_file(input_file);
+ CHECK_ERR(rval);
+
+}
+
+int main(int /* argc */, char** /* argv */)
+{
+ int result = RUN_TEST( read_multiple_test );
+
+ return result;
+}
diff --git a/test/io/read_cgm_test.cpp b/test/io/read_cgm_test.cpp
deleted file mode 100644
index 5377220..0000000
--- a/test/io/read_cgm_test.cpp
+++ /dev/null
@@ -1,49 +0,0 @@
-
-#include <iostream>
-#include "moab/Interface.hpp"
-#ifndef IS_BUILDING_MB
-#define IS_BUILDING_MB
-#endif
-#include "TestUtil.hpp"
-#include "Internals.hpp"
-#include "moab/Core.hpp"
-
-using namespace moab;
-
-#define CHKERR(A) do { if (MB_SUCCESS != (A)) { \
- std::cerr << "Failure (error code " << (A) << ") at " __FILE__ ":" \
- << __LINE__ << std::endl; \
- return A; } } while(false)
-
-#ifdef MESHDIR
-#ifdef HAVE_OCC_STEP
-static const char input_file[] = STRINGIFY(MESHDIR) "/io/dum.stp";
-#else
-static const char input_file[] = STRINGIFY(MESHDIR) "/io/dum.sat";
-#endif
-#else
-#ifdef HAVE_OCC_STEP
-static const char input_file[] = "dum.stp";
-#else
-static const char input_file[] = "dum.sat";
-#endif
-#endif
-
-void read_multiple_test()
-{
- Core mb;
-
- ErrorCode rval = mb.load_file(input_file);
- CHECK_ERR(rval);
- // second load
- rval = mb.load_file(input_file);
- CHECK_ERR(rval);
-
-}
-
-int main(int /* argc */, char** /* argv */)
-{
- int result = RUN_TEST( read_multiple_test );
-
- return result;
-}
https://bitbucket.org/fathomteam/moab/commits/c9529c4759d3/
Changeset: c9529c4759d3
Branch: None
User: tautges
Date: 2014-02-09 17:16:34
Summary: Moving TupleList class from parallel to non-parallel location, since it's now used in SpatialLocator serial version too.
Mostly passes make check serial/parallel (though a few other problems appear there...)
Affected #: 6 files
diff --git a/src/Makefile.am b/src/Makefile.am
index 0d474b0..0a05543 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -124,6 +124,7 @@ libMOAB_la_SOURCES = \
TagInfo.cpp \
TagInfo.hpp \
Tree.cpp \
+ TupleList.cpp \
Types.cpp \
TypeSequenceManager.cpp \
TypeSequenceManager.hpp \
@@ -195,6 +196,7 @@ nobase_libMOAB_la_include_HEADERS = \
moab/SpectralMeshTool.hpp \
moab/Tree.hpp \
moab/TreeStats.hpp \
+ moab/TupleList.hpp \
moab/Types.hpp \
moab/UnknownInterface.hpp \
moab/Util.hpp \
diff --git a/src/TupleList.cpp b/src/TupleList.cpp
new file mode 100644
index 0000000..bf08946
--- /dev/null
+++ b/src/TupleList.cpp
@@ -0,0 +1,860 @@
+#include <string.h>
+#include <limits.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <stdarg.h>
+#include <iostream>
+
+#include "moab/TupleList.hpp"
+
+namespace moab {
+
+void fail(const char *fmt, ...)
+{
+ va_list ap;
+ va_start(ap, fmt);
+ vfprintf(stderr, fmt, ap);
+ va_end(ap);
+ exit(1);
+}
+
+TupleList::buffer::buffer(size_t sz)
+{
+ ptr = NULL;
+ buffSize = 0;
+ this->buffer_init_(sz, __FILE__);
+}
+
+TupleList::buffer::buffer()
+{
+ buffSize = 0;
+ ptr = NULL;
+}
+
+void TupleList::buffer::buffer_init_(size_t sizeIn, const char *file)
+{
+ this->buffSize = sizeIn;
+ void *res = malloc(this->buffSize);
+ if (!res && buffSize > 0)
+ fail("%s: allocation of %d bytes failed\n", file, (int) buffSize);
+ ptr = (char*) res;
+}
+
+void TupleList::buffer::buffer_reserve_(size_t min, const char *file)
+{
+ if (this->buffSize < min)
+ {
+ size_t newSize = this->buffSize;
+ newSize += newSize / 2 + 1;
+ if (newSize < min)
+ newSize = min;
+ void *res = realloc(ptr, newSize);
+ if (!res && newSize > 0)
+ fail("%s: allocation of %d bytes failed\n", file, (int) this->buffSize);
+ ptr = (char*) res;
+ this->buffSize = newSize;
+ }
+}
+
+void TupleList::buffer::reset()
+{
+ free(ptr);
+ ptr = NULL;
+ buffSize = 0;
+}
+
+TupleList::TupleList(uint p_mi, uint p_ml, uint p_mul, uint p_mr, uint p_max)
+{
+ vi = NULL;
+ vl = NULL;
+ vul = NULL;
+ vr = NULL;
+ initialize(p_mi, p_ml, p_mul, p_mr, p_max);
+}
+
+TupleList::TupleList()
+{
+ vi = NULL;
+ vl = NULL;
+ vul = NULL;
+ vr = NULL;
+ disableWriteAccess();
+}
+
+// Allocates space for the tuple list in memory according to parameters
+void TupleList::initialize(uint p_mi, uint p_ml, uint p_mul, uint p_mr,
+ uint p_max)
+{
+ this->n = 0;
+ this->max = p_max;
+ this->mi = p_mi;
+ this->ml = p_ml;
+ this->mul = p_mul;
+ this->mr = p_mr;
+ size_t sz;
+
+ if (max * mi > 0)
+ {
+ sz = max * mi * sizeof(sint);
+ void *resi = malloc(sz);
+ if (!resi && max * mi > 0)
+ fail("%s: allocation of %d bytes failed\n", __FILE__, (int) sz);
+ vi = (sint*) resi;
+ }
+ else
+ vi = NULL;
+ if (max * ml > 0)
+ {
+ sz = max * ml * sizeof(slong);
+ void *resl = malloc(sz);
+ if (!resl && max * ml > 0)
+ fail("%s: allocation of %d bytes failed\n", __FILE__, (int) sz);
+ vl = (slong*) resl;
+ }
+ else
+ vl = NULL;
+ if (max * mul > 0)
+ {
+ sz = max * mul * sizeof(ulong);
+ void *resu = malloc(sz);
+ if (!resu && max * mul > 0)
+ fail("%s: allocation of %d bytes failed\n", __FILE__, (int) sz);
+ vul = (ulong*) resu;
+ }
+ else
+ vul = NULL;
+ if (max * mr > 0)
+ {
+ sz = max * mr * sizeof(realType);
+ void *resr = malloc(sz);
+ if (!resr && max * ml > 0)
+ fail("%s: allocation of %d bytes failed\n", __FILE__, (int) sz);
+ vr = (realType*) resr;
+ }
+ else
+ vr = NULL;
+
+ //Begin with write access disabled
+ this->disableWriteAccess();
+
+ //Set read variables
+ vi_rd = vi;
+ vl_rd = vl;
+ vul_rd = vul;
+ vr_rd = vr;
+}
+
+// Resizes a tuplelist to the given uint max
+ErrorCode TupleList::resize(uint maxIn)
+{
+ this->max = maxIn;
+ size_t sz;
+
+ if (vi || (max * mi > 0))
+ {
+ sz = max * mi * sizeof(sint);
+ void *resi = realloc(vi, sz);
+ if (!resi && max * mi > 0)
+ {
+ fail("%s: allocation of %d bytes failed\n", __FILE__, (int) sz);
+ return moab::MB_MEMORY_ALLOCATION_FAILED;
+ }
+ vi = (sint*) resi;
+ }
+ if (vl || (max * ml > 0))
+ {
+ sz = max * ml * sizeof(slong);
+ void *resl = realloc(vl, sz);
+ if (!resl && max * ml > 0)
+ {
+ fail("%s: allocation of %d bytes failed\n", __FILE__, (int) sz);
+ return moab::MB_MEMORY_ALLOCATION_FAILED;
+ }
+ vl = (slong*) resl;
+ }
+ if (vul || (max * mul > 0))
+ {
+ sz = max * mul * sizeof(ulong);
+ void *resu = realloc(vul, sz);
+ if (!resu && max * mul > 0)
+ {
+ fail("%s: allocation of %d bytes failed\n", __FILE__, (int) sz);
+ return moab::MB_MEMORY_ALLOCATION_FAILED;
+ }
+ vul = (ulong*) resu;
+ }
+ if (vr || (max * mr > 0))
+ {
+ sz = max * mr * sizeof(realType);
+ void *resr = realloc(vr, sz);
+ if (!resr && max * mr > 0)
+ {
+ fail("%s: allocation of %d bytes failed\n", __FILE__, (int) sz);
+ return moab::MB_MEMORY_ALLOCATION_FAILED;
+ }
+ vr = (realType*) resr;
+ }
+
+ //Set read variables
+ vi_rd = vi;
+ vl_rd = vl;
+ vul_rd = vul;
+ vr_rd = vr;
+
+ //Set the write variables if necessary
+ if (writeEnabled)
+ {
+ vi_wr = vi;
+ vl_wr = vl;
+ vul_wr = vul;
+ vr_wr = vr;
+ }
+ return moab::MB_SUCCESS;
+}
+
+// Frees the memory used by the tuplelist
+void TupleList::reset()
+{
+ //free up the pointers
+ free(vi);
+ free(vl);
+ free(vul);
+ free(vr);
+ //Set them all to null
+ vr = NULL;
+ vi = NULL;
+ vul = NULL;
+ vl = NULL;
+ //Set the read and write pointers to null
+ disableWriteAccess();
+ vi_rd = NULL;
+ vl_rd = NULL;
+ vul_rd = NULL;
+ vr_rd = NULL;
+}
+
+// Increments n; if n>max, increase the size of the tuplelist
+void TupleList::reserve()
+{
+ n++;
+ while (n > max)
+ resize((max ? max + max / 2 + 1 : 2));
+ last_sorted = -1;
+}
+
+// Given the value and the position in the field, finds the index of the tuple
+// to which the value belongs
+int TupleList::find(unsigned int key_num, sint value)
+{
+ ulong uvalue = (ulong) value;
+ if (!(key_num > mi))
+ {
+ // Binary search: only if the tuple_list is sorted
+ if (last_sorted == (int) key_num)
+ {
+ int lb = 0, ub = n, index; // lb=lower bound, ub=upper bound, index=mid
+ for (; lb <= ub;)
+ {
+ index = (lb + ub) / 2;
+ if (vi[index * mi + key_num] == (long) uvalue)
+ return index;
+ else if (vi[index * mi + key_num] > (long) uvalue)
+ ub = index - 1;
+ else if (vi[index * mi + key_num] < (long) uvalue)
+ lb = index + 1;
+ }
+ }
+ else
+ {
+ // Sequential search: if tuple_list is not sorted
+ for (long index = 0; index < n; index++)
+ {
+ if (vi[index * mi + key_num] == (long) uvalue)
+ return index;
+ }
+ }
+ }
+ return -1; // If the value wasn't present or an invalid key was given
+}
+
+int TupleList::find(unsigned int key_num, slong value)
+{
+ ulong uvalue = (ulong) value;
+ if (!(key_num > ml))
+ {
+ if (last_sorted - mi == key_num)
+ {
+ int lb = 0, ub = n, index; // lb=lower bound, ub=upper bound, index=mid
+ for (; lb <= ub;)
+ {
+ index = (lb + ub) / 2;
+ if (vl[index * ml + key_num] == (long) uvalue)
+ return index;
+ else if (vl[index * ml + key_num] > (long) uvalue)
+ ub = index - 1;
+ else if (vl[index * ml + key_num] < (long) uvalue)
+ lb = index + 1;
+ }
+ }
+ else
+ {
+ // Sequential search: if tuple_list is not sorted
+ for (uint index = 0; index < n; index++)
+ {
+ if (vl[index * ml + key_num] == (long) uvalue)
+ return index;
+ }
+ }
+ }
+ return -1; // If the value wasn't present or an invalid key was given
+}
+
+int TupleList::find(unsigned int key_num, ulong value)
+{
+ if (!(key_num > mul))
+ {
+ if (last_sorted - mi - ml == key_num)
+ {
+ int lb = 0, ub = n - 1, index; // lb=lower bound, ub=upper bound, index=mid
+ for (; lb <= ub;)
+ {
+ index = (lb + ub) / 2;
+ if (vul[index * mul + key_num] == value)
+ return index;
+ else if (vul[index * mul + key_num] > value)
+ ub = index - 1;
+ else if (vul[index * mul + key_num] < value)
+ lb = index + 1;
+ }
+ }
+ else
+ {
+ // Sequential search: if tuple_list is not sorted
+ for (uint index = 0; index < n; index++)
+ {
+ if (vul[index * mul + key_num] == value)
+ return index;
+ }
+ }
+ }
+ return -1; // If the value wasn't present or an invalid key was given
+}
+
+int TupleList::find(unsigned int key_num, realType value)
+{
+ if (!key_num > mr)
+ {
+ // Sequential search: TupleList cannot be sorted by reals
+ for (uint index = 0; index < n; index++)
+ {
+ if (vr[index * mr + key_num] == value)
+ return index;
+ }
+ }
+ return -1; // If the value wasn't present or an invalid key was given
+}
+
+sint TupleList::get_sint(unsigned int index, unsigned int m)
+{
+ if (mi > m && n > index)
+ return vi[index * mi + m];
+ return 0;
+}
+
+slong TupleList::get_int(unsigned int index, unsigned int m)
+{
+ if (ml > m && n > index)
+ return vl[index * ml + m];
+ return 0;
+}
+
+ulong TupleList::get_ulong(unsigned int index, unsigned int m)
+{
+ if (mul > m && n > index)
+ return vul[index * mul + m];
+ return 0;
+}
+
+realType TupleList::get_double(unsigned int index, unsigned int m)
+{
+ if (mr > m && n > index)
+ return vr[index * mr + m];
+ return 0;
+}
+
+ErrorCode TupleList::get(unsigned int index, const sint *&sp, const slong *&ip,
+ const ulong *&lp, const realType *&dp)
+{
+ if (index <= n)
+ {
+ if (mi)
+ *&sp = &vi[index * mi];
+ else
+ *&sp = NULL;
+ if (ml)
+ *&ip = &vl[index * ml];
+ else
+ *&ip = NULL;
+ if (mul)
+ *&lp = &vul[index * mul];
+ else
+ *&lp = NULL;
+ if (mr)
+ *&dp = &vr[index * mr];
+ else
+ *&dp = NULL;
+
+ return MB_SUCCESS;
+ }
+ return MB_FAILURE;
+}
+
+unsigned int TupleList::push_back(sint *sp, slong *ip, ulong *lp, realType *dp)
+{
+ reserve();
+ if (mi)
+ memcpy(&vi[mi * (n - 1)], sp, mi * sizeof(sint));
+ if (ml)
+ memcpy(&vl[ml * (n - 1)], ip, ml * sizeof(long));
+ if (mul)
+ memcpy(&vul[mul * (n - 1)], lp, mul * sizeof(ulong));
+ if (mr)
+ memcpy(&vr[mr * (n - 1)], dp, mr * sizeof(realType));
+
+ last_sorted = -1;
+ return n - 1;
+}
+
+void TupleList::enableWriteAccess()
+{
+ writeEnabled = true;
+ last_sorted = -1;
+ vi_wr = vi;
+ vl_wr = vl;
+ vul_wr = vul;
+ vr_wr = vr;
+}
+
+void TupleList::disableWriteAccess()
+{
+ writeEnabled = false;
+ vi_wr = NULL;
+ vl_wr = NULL;
+ vul_wr = NULL;
+ vr_wr = NULL;
+}
+
+void TupleList::getTupleSize(uint &mi_out, uint &ml_out, uint &mul_out,
+ uint &mr_out) const
+{
+ mi_out = mi;
+ ml_out = ml;
+ mul_out = mul;
+ mr_out = mr;
+}
+
+uint TupleList::inc_n()
+{
+ //Check for direct write access
+ if (!writeEnabled)
+ {
+ enableWriteAccess();
+ }
+ n++;
+ return n;
+}
+
+void TupleList::set_n(uint n_in)
+{
+ //Check for direct write access;
+ if (!writeEnabled)
+ {
+ enableWriteAccess();
+ }
+ n = n_in;
+}
+
+void TupleList::print(const char *name) const
+{
+ std::cout << "Printing Tuple " << name << "===================" << std::endl;
+ unsigned long i = 0, l = 0, ul = 0, r = 0;
+ for (uint k = 0; k < n; k++)
+ {
+ for (uint j = 0; j < mi; j++)
+ {
+ std::cout << vi[i++] << " | ";
+ }
+ for (uint j = 0; j < ml; j++)
+ {
+ std::cout << vl[l++] << " | ";
+ }
+ for (uint j = 0; j < mul; j++)
+ {
+ std::cout << vul[ul++] << " | ";
+ }
+ for (uint j = 0; j < mr; j++)
+ {
+ std::cout << vr[r++] << " | ";
+ }
+ std::cout << std::endl;
+ }
+ std::cout << "=======================================" << std::endl
+ << std::endl;
+}
+
+void TupleList::permute(uint *perm, void *work)
+{
+ const unsigned int_size = mi * sizeof(sint), long_size = ml * sizeof(slong),
+ ulong_size = mul * sizeof(ulong), real_size = mr * sizeof(realType);
+ if (mi)
+ {
+ uint *p = perm, *pe = p + n;
+ char *sorted = (char *) work;
+ while (p != pe)
+ memcpy((void *) sorted, &vi[mi * (*p++)], int_size), sorted += int_size;
+ memcpy(vi, work, int_size * n);
+ }
+ if (ml)
+ {
+ uint *p = perm, *pe = p + n;
+ char *sorted = (char *) work;
+ while (p != pe)
+ memcpy((void *) sorted, &vl[ml * (*p++)], long_size), sorted += long_size;
+ memcpy(vl, work, long_size * n);
+ }
+ if (mul)
+ {
+ uint *p = perm, *pe = p + n;
+ char *sorted = (char *) work;
+ while (p != pe)
+ memcpy((void *) sorted, &vul[mul * (*p++)], ulong_size), sorted +=
+ ulong_size;
+ memcpy(vul, work, ulong_size * n);
+ }
+ if (mr)
+ {
+ uint *p = perm, *pe = p + n;
+ char *sorted = (char *) work;
+ while (p != pe)
+ memcpy((void *) sorted, &vr[mr * (*p++)], real_size), sorted += real_size;
+ memcpy(vr, work, real_size * n);
+ }
+}
+
+# define umax_2(a, b) (((a)>(b)) ? (a):(b))
+
+ErrorCode TupleList::sort(uint key, TupleList::buffer *buf)
+{
+ const unsigned int_size = mi * sizeof(sint);
+ const unsigned long_size = ml * sizeof(slong);
+ const unsigned ulong_size = mul * sizeof(ulong);
+ const unsigned real_size = mr * sizeof(realType);
+ const unsigned width = umax_2(umax_2(int_size,long_size),
+ umax_2(ulong_size,real_size));
+ const unsigned data_size =
+ key >= mi ? sizeof(SortData<long> ) : sizeof(SortData<uint> );
+
+ uint work_min = n * umax_2(2*data_size,sizeof(sint)+width);
+ uint *work;
+ buf->buffer_reserve(work_min);
+ work = (uint *) buf->ptr;
+ if (key < mi)
+ index_sort((uint *) &vi[key], n, mi, work, (SortData<uint>*) work);
+ else if (key < mi + ml)
+ index_sort((long*) &vl[key - mi], n, ml, work, (SortData<long>*) work);
+ else if (key < mi + ml + mul)
+ index_sort((ulong*) &vul[key - mi - ml], n, mul, work,
+ (SortData<ulong>*) work);
+ else
+ return MB_NOT_IMPLEMENTED;
+
+ permute(work, work + n);
+
+ if (!writeEnabled)
+ last_sorted = key;
+ return MB_SUCCESS;
+}
+
+#undef umax_2
+
+#define DIGIT_BITS 8
+#define DIGIT_VALUES (1<<DIGIT_BITS)
+#define DIGIT_MASK ((Value)(DIGIT_VALUES-1))
+#define CEILDIV(a,b) (((a)+(b)-1)/(b))
+#define DIGITS CEILDIV(CHAR_BIT*sizeof(Value),DIGIT_BITS)
+#define VALUE_BITS (DIGIT_BITS*DIGITS)
+#define COUNT_SIZE (DIGITS*DIGIT_VALUES)
+
+/* used to unroll a tiny loop: */
+#define COUNT_DIGIT_01(n,i) \
+ if(n>i) count[i][val&DIGIT_MASK]++, val>>=DIGIT_BITS
+#define COUNT_DIGIT_02(n,i) COUNT_DIGIT_01(n,i); COUNT_DIGIT_01(n,i+ 1)
+#define COUNT_DIGIT_04(n,i) COUNT_DIGIT_02(n,i); COUNT_DIGIT_02(n,i+ 2)
+#define COUNT_DIGIT_08(n,i) COUNT_DIGIT_04(n,i); COUNT_DIGIT_04(n,i+ 4)
+#define COUNT_DIGIT_16(n,i) COUNT_DIGIT_08(n,i); COUNT_DIGIT_08(n,i+ 8)
+#define COUNT_DIGIT_32(n,i) COUNT_DIGIT_16(n,i); COUNT_DIGIT_16(n,i+16)
+#define COUNT_DIGIT_64(n,i) COUNT_DIGIT_32(n,i); COUNT_DIGIT_32(n,i+32)
+
+template<class Value>
+Value TupleList::radix_count(const Value *A, const Value *end, Index stride,
+ Index count[DIGITS][DIGIT_VALUES])
+{
+ Value bitorkey = 0;
+ memset(count, 0, COUNT_SIZE * sizeof(Index));
+ do
+ {
+ Value val = *A;
+ bitorkey |= val;
+ COUNT_DIGIT_64(DIGITS, 0);
+ // above macro expands to:
+ //if(DIGITS> 0) count[ 0][val&DIGIT_MASK]++, val>>=DIGIT_BITS;
+ //if(DIGITS> 1) count[ 1][val&DIGIT_MASK]++, val>>=DIGIT_BITS;
+ // ...
+ //if(DIGITS>63) count[63][val&DIGIT_MASK]++, val>>=DIGIT_BITS;
+
+ } while (A += stride, A != end);
+ return bitorkey;
+}
+
+#undef COUNT_DIGIT_01
+#undef COUNT_DIGIT_02
+#undef COUNT_DIGIT_04
+#undef COUNT_DIGIT_08
+#undef COUNT_DIGIT_16
+#undef COUNT_DIGIT_32
+#undef COUNT_DIGIT_64
+
+void TupleList::radix_offsets(Index *c)
+{
+ Index sum = 0, t, *ce = c + DIGIT_VALUES;
+ do
+ t = *c, *c++ = sum, sum += t;
+ while (c != ce);
+}
+
+template<class Value>
+unsigned TupleList::radix_zeros(Value bitorkey,
+ Index count[DIGITS][DIGIT_VALUES], unsigned *shift, Index **offsets)
+{
+ unsigned digits = 0, sh = 0;
+ Index *c = &count[0][0];
+ do
+ {
+ if (bitorkey & DIGIT_MASK)
+ *shift++ = sh, *offsets++ = c, ++digits, radix_offsets(c);
+ } while (bitorkey >>= DIGIT_BITS,sh += DIGIT_BITS,c += DIGIT_VALUES,sh
+ != VALUE_BITS);
+ return digits;
+}
+
+template<class Value>
+void TupleList::radix_index_pass_b(const Value *A, Index n, Index stride,
+ unsigned sh, Index *off, SortData<Value> *out)
+{
+ Index i = 0;
+ do
+ {
+ Value v = *A;
+ SortData<Value> *d = &out[off[(v >> sh) & DIGIT_MASK]++];
+ d->v = v, d->i = i++;
+ } while (A += stride, i != n);
+}
+
+template<class Value>
+void TupleList::radix_index_pass_m(const SortData<Value> *src,
+ const SortData<Value> *end, unsigned sh, Index *off, SortData<Value> *out)
+{
+ do
+ {
+ SortData<Value> *d = &out[off[(src->v >> sh) & DIGIT_MASK]++];
+ d->v = src->v, d->i = src->i;
+ } while (++src != end);
+}
+
+template<class Value>
+void TupleList::radix_index_pass_e(const SortData<Value> *src,
+ const SortData<Value> *end, unsigned sh, Index *off, Index *out)
+{
+ do
+ out[off[(src->v >> sh) & DIGIT_MASK]++] = src->i;
+ while (++src != end);
+}
+
+template<class Value>
+void TupleList::radix_index_pass_be(const Value *A, Index n, Index stride,
+ unsigned sh, Index *off, Index *out)
+{
+ Index i = 0;
+ do
+ out[off[(*A >> sh) & DIGIT_MASK]++] = i++;
+ while (A += stride, i != n);
+}
+
+template<class Value>
+void TupleList::radix_index_sort(const Value *A, Index n, Index stride,
+ Index *idx, SortData<Value> *work)
+{
+ Index count[DIGITS][DIGIT_VALUES];
+ Value bitorkey = radix_count(A, A + n * stride, stride, count);
+ unsigned shift[DIGITS];
+ Index *offsets[DIGITS];
+ unsigned digits = radix_zeros(bitorkey, count, shift, offsets);
+ if (digits == 0)
+ {
+ Index i = 0;
+ do
+ *idx++ = i++;
+ while (i != n);
+ }
+ else if (digits == 1)
+ {
+ radix_index_pass_be(A, n, stride, shift[0], offsets[0], idx);
+ }
+ else
+ {
+ SortData<Value> *src, *dst;
+ unsigned d;
+ if ((digits & 1) == 0)
+ dst = work, src = dst + n;
+ else
+ src = work, dst = src + n;
+ radix_index_pass_b(A, n, stride, shift[0], offsets[0], src);
+ for (d = 1; d != digits - 1; ++d)
+ {
+ SortData<Value> *t;
+ radix_index_pass_m(src, src + n, shift[d], offsets[d], dst);
+ t = src, src = dst, dst = t;
+ }
+ radix_index_pass_e(src, src + n, shift[d], offsets[d], idx);
+ }
+}
+
+template<class Value>
+void TupleList::merge_index_sort(const Value *A, const Index An, Index stride,
+ Index *idx, SortData<Value> *work)
+{
+ SortData<Value> * const buf[2] = { work + An, work };
+ Index n = An, base = -n, odd = 0, c = 0, b = 1;
+ Index i = 0;
+ for (;;)
+ {
+ SortData<Value> *p;
+ if ((c & 1) == 0)
+ {
+ base += n, n += (odd & 1), c |= 1, b ^= 1;
+ while (n > 3)
+ odd <<= 1, odd |= (n & 1), n >>= 1, c <<= 1, b ^= 1;
+ }
+ else
+ base -= n - (odd & 1), n <<= 1, n -= (odd & 1), odd >>= 1, c >>= 1;
+ if (c == 0)
+ break;
+ p = buf[b] + base;
+ if (n == 2)
+ {
+ Value v[2];
+ v[0] = *A, A += stride, v[1] = *A, A += stride;
+ if (v[1] < v[0])
+ p[0].v = v[1], p[0].i = i + 1, p[1].v = v[0], p[1].i = i;
+ else
+ p[0].v = v[0], p[0].i = i, p[1].v = v[1], p[1].i = i + 1;
+ i += 2;
+ }
+ else if (n == 3)
+ {
+ Value v[3];
+ v[0] = *A, A += stride, v[1] = *A, A += stride, v[2] = *A, A += stride;
+ if (v[1] < v[0])
+ {
+ if (v[2] < v[1])
+ p[0].v = v[2], p[1].v = v[1], p[2].v = v[0], p[0].i = i + 2, p[1].i =
+ i + 1, p[2].i = i;
+ else
+ {
+ if (v[2] < v[0])
+ p[0].v = v[1], p[1].v = v[2], p[2].v = v[0], p[0].i = i + 1, p[1].i =
+ i + 2, p[2].i = i;
+ else
+ p[0].v = v[1], p[1].v = v[0], p[2].v = v[2], p[0].i = i + 1, p[1].i =
+ i, p[2].i = i + 2;
+ }
+ }
+ else
+ {
+ if (v[2] < v[0])
+ p[0].v = v[2], p[1].v = v[0], p[2].v = v[1], p[0].i = i + 2, p[1].i =
+ i, p[2].i = i + 1;
+ else
+ {
+ if (v[2] < v[1])
+ p[0].v = v[0], p[1].v = v[2], p[2].v = v[1], p[0].i = i, p[1].i = i
+ + 2, p[2].i = i + 1;
+ else
+ p[0].v = v[0], p[1].v = v[1], p[2].v = v[2], p[0].i = i, p[1].i = i
+ + 1, p[2].i = i + 2;
+ }
+ }
+ i += 3;
+ }
+ else
+ {
+ const Index na = n >> 1, nb = (n + 1) >> 1;
+ const SortData<Value> *ap = buf[b ^ 1] + base, *ae = ap + na;
+ SortData<Value> *bp = p + na, *be = bp + nb;
+ for (;;)
+ {
+ if (bp->v < ap->v)
+ {
+ *p++ = *bp++;
+ if (bp != be)
+ continue;
+ do
+ *p++ = *ap++;
+ while (ap != ae);
+ break;
+ }
+ else
+ {
+ *p++ = *ap++;
+ if (ap == ae)
+ break;
+ }
+ }
+ }
+ }
+ {
+ const SortData<Value> *p = buf[0], *pe = p + An;
+ do
+ *idx++ = (p++)->i;
+ while (p != pe);
+ }
+}
+
+template<class Value>
+void TupleList::index_sort(const Value *A, Index n, Index stride, Index *idx,
+ SortData<Value> *work)
+{
+ if (n < DIGIT_VALUES)
+ {
+ if (n == 0)
+ return;
+ if (n == 1)
+ *idx = 0;
+ else
+ merge_index_sort(A, n, stride, idx, work);
+ }
+ else
+ radix_index_sort(A, n, stride, idx, work);
+}
+
+#undef DIGIT_BITS
+#undef DIGIT_VALUES
+#undef DIGIT_MASK
+#undef CEILDIV
+#undef DIGITS
+#undef VALUE_BITS
+#undef COUNT_SIZE
+#undef sort_data_long
+
+} //namespace
+
diff --git a/src/moab/TupleList.hpp b/src/moab/TupleList.hpp
new file mode 100644
index 0000000..8cbd388
--- /dev/null
+++ b/src/moab/TupleList.hpp
@@ -0,0 +1,413 @@
+/*-----------------------------------------------------------------------------
+
+ Tuple list definition and utilities
+
+ Conceptually, a tuple list is a list of n records or tuples,
+ each with mi integers, ml longs, mul ulongs, and mr reals
+ (these types are defined in "types.h" as sint, slong, sulong, real;
+ it may be that sint==slong)
+
+ There are four arrays, one for each type (vi,vl,vul,vr),
+ with records layed out contiguously within each array
+
+ -----------------------------------------------------------------------------*/
+
+#ifndef TUPLE_LIST_HPP
+#define TUPLE_LIST_HPP
+
+#include <limits.h>
+#include <stdlib.h>
+
+#include "moab/Types.hpp"
+#include <string>
+
+/* Integral types defined here to ensure variable type sizes are consistent */
+/* integer type to use for everything */
+#if defined(USE_LONG)
+# define INTEGER long
+#elif defined(USE_LONG_LONG)
+# define INTEGER long long
+#elif defined(USE_SHORT)
+# define INTEGER short
+#else
+# define INTEGER int
+#endif
+
+/* when defined, use the given type for global indices instead of INTEGER */
+#if defined(USE_GLOBAL_LONG_LONG)
+# define GLOBAL_INT long long
+#elif defined(USE_GLOBAL_LONG)
+# define GLOBAL_INT long
+#else
+# define GLOBAL_INT long
+#endif
+
+/* floating point type to use for everything */
+#if defined(USE_FLOAT)
+typedef float realType;
+# define floorr floorf
+# define ceilr ceilf
+# define sqrtr sqrtf
+# define fabsr fabsf
+# define cosr cosf
+# define sinr sinf
+# define EPS (128*FLT_EPSILON)
+# define PI 3.1415926535897932384626433832795028841971693993751058209749445923F
+#elif defined(USE_LONG_DOUBLE)
+typedef long double realType;
+# define floorr floorl
+# define ceilr ceill
+# define sqrtr sqrtl
+# define fabsr fabsl
+# define cosr cosl
+# define sinr sinl
+# define EPS (128*LDBL_EPSILON)
+# define PI 3.1415926535897932384626433832795028841971693993751058209749445923L
+#else
+typedef double realType;
+# define floorr floor
+# define ceilr ceil
+# define sqrtr sqrt
+# define fabsr fabs
+# define cosr cos
+# define sinr sin
+# define EPS (128*DBL_EPSILON)
+# define PI 3.1415926535897932384626433832795028841971693993751058209749445923
+#endif
+
+/* apparently uint and ulong can be defined already in standard headers */
+#define uint uint_
+#define ulong ulong_
+#define sint sint_
+#define slong slong_
+
+typedef signed INTEGER sint;
+typedef unsigned INTEGER uint;
+#undef INTEGER
+
+#ifdef GLOBAL_INT
+typedef signed GLOBAL_INT slong;
+typedef unsigned GLOBAL_INT ulong;
+#else
+typedef sint slong;
+typedef uint ulong;
+#endif
+
+
+
+namespace moab
+{
+ void fail(const char *fmt, ...);
+
+ class TupleList
+ {
+ public:
+ /*---------------------------------------------------------------------------
+
+ buffer: a simple class which can be used to store data
+
+ The ptr points to the chunk of memory allocated for the buffer's use. The
+ size denotes the size of the allocated memory; the user must take care to
+ note how much of the buffer they are using.
+
+ ---------------------------------------------------------------------------*/
+ class buffer
+ {
+ public:
+ size_t buffSize;
+ char *ptr;
+
+ /**Constructor which sets an initial capacity of the buffer
+ */
+ buffer(size_t sz);
+
+ /**Default constructor (Note: buffer must be initialized before use!)
+ */
+ buffer();
+
+ ~buffer () { this->reset(); };
+
+ /**Initializes the buffer to have a capacity of size
+ */
+ void buffer_init_(size_t sz, const char *file);
+
+ /**Ensures that the buffer has at least a capacity of min
+ */
+ void buffer_reserve_(size_t min, const char *file);
+
+ /**Frees any allocated memory used by the buffer
+ */
+ void reset ();
+
+ //Aliases for using the buffer methods
+#define buffer_init(sz) buffer_init_(sz,__FILE__)
+#define buffer_reserve(min) buffer_reserve_(min,__FILE__)
+
+ };
+
+ public:
+
+ /**Constructor that takes all parameters and initializes the TupleList
+ */
+ TupleList(uint mi, uint ml,
+ uint mul, uint mr, uint max);
+
+ /**Default constructor (Note: TupleList must be initialized before use!)
+ */
+ TupleList();
+
+ ~TupleList() { reset(); };
+
+ /**Initializes the starting memory to be used by the TupleList
+ * Note: TupleLists must be initialized before they can be used
+ *
+ * param mi number of signed ints in each tuple
+ * param ml number of long ints in each tuple
+ * param mul number of unsigned long ints in each tuple
+ * param mr number of reals in each tuple
+ * param max starting capacity of max tuples in the TupleList
+ */
+ void initialize(uint mi, uint ml,
+ uint mul, uint mr, uint max);
+
+ /**Resizes the TupleList to a new max
+ *
+ * param max the new max capacity of the TupleList
+ */
+ ErrorCode resize(uint max);
+
+ /**Sorts the TupleList by 'key'
+ * if key<mi: key represents the key'th index in vi
+ * if mi<key<ml: key represents the (key-mi)'th index in vl
+ * if ml<key<mul: key represents the (key-mi-ml)'th index in vul
+ *
+ * param key index to be sorted by
+ * param *buf buffer space used for sorting
+ */
+ /*------------------------------------------------------------------------------
+
+ Hybrid Stable Sort
+
+ low-overhead merge sort when n is small,
+ otherwise asymptotically superior radix sort
+
+ result = O(n) sort with good performance for all n
+
+ A, n, stride : specifices the input
+
+ sort:
+ uint out[n] : the sorted values (output)
+ uint work[n]: scratch area
+
+ index_sort:
+ uint idx[n] : the sorted indices (output)
+ sort_data work[2*n]: scratch area
+
+ ----------------------------------------------------------------------------*/
+ ErrorCode sort(uint key, TupleList::buffer *buf);
+
+ /**Frees all allocated memory in use by the TupleList
+ */
+ void reset();
+
+ /**Adds one to the total number of in-use tuples and resizes if necessary
+ */
+ void reserve();
+
+ /**Finds index of the tuple containing 'value' at the key_numth index of
+ * said tuple; return -1 if key_num is out of bounds or if 'value' not found
+ * Uses binary search if TupleList is sorted by the key_numth field, seqential
+ * otherwise (very slow for large TupleLists; please sort before you search)
+ *
+ * param key_num index of the tuple where to search for value
+ * param value value to search for at the given key_num
+ * return the index of the tuple that contains value
+ */
+ int find(unsigned int key_num, sint value);
+ int find(unsigned int key_num, slong value);
+ int find(unsigned int key_num, ulong value);
+ int find(unsigned int key_num, realType value);
+
+ /**get the mth number of return type in the index'th tuple
+ * returns 0 if m or index is out of bounds
+ *
+ * param index index of the tuple within the TupleList
+ * param m index of the value within the tuple
+ * return the value at the given position
+ */
+ sint get_sint(unsigned int index, unsigned int m);
+ slong get_int(unsigned int index, unsigned int m);
+ ulong get_ulong(unsigned int index, unsigned int m);
+ realType get_double(unsigned int index, unsigned int m);
+
+ /**get pointers to the data for the index'th tuple; ptr is
+ * NULL if that type is not part of this tuple
+ *
+ * param index index of the tuple needed
+ * param *&sp, *&ip, *&lp, *&dp pointers to each piece of the tuple
+ */
+ ErrorCode get(unsigned int index, const sint *&sp,
+ const slong *&ip, const ulong *&lp, const realType *&dp);
+
+ /**push back a new tuple on the TupleList;
+ *
+ * param *sp pointer to a list of signed ints
+ * param *ip pointer to a list of signed longs
+ * param *lp pointer to a list of unsigned longs
+ * param *dp pointer to a list of reals
+ * return index of that tuple
+ */
+ unsigned int push_back(sint *sp, slong *ip,
+ ulong *lp, realType *dp);
+
+ /*Enable or disable direct write access to arrays
+ This is important so that we know whether or not
+ the user of this class can directly modify data
+ which can affect operations such as
+ whether or not we know the tuple list is sorted
+ (for a binary search)*/
+ void enableWriteAccess();
+ void disableWriteAccess();
+
+ /*Get information on the Tuple Sizes
+ * param &mi_out Count of uints in a tuple
+ * param &ml_out Count of uints in a tuple
+ * param &mul_out Count of uints in a tuple
+ * param &mr_out Count of uints in a tuple
+ */
+ void getTupleSize(uint &mi_out, uint &ml_out,
+ uint &mul_out, uint &mr_out) const;
+
+ /*Set the count of Tuples in the Tuple List
+ * Warning, automatically calls enableWriteAccess()
+ * param n_in New count of Tuples
+ */
+ void set_n(uint n_in);
+
+ /* Get the count of Tuples in the Tuple List */
+ uint get_n() const;
+
+ /*Get the maximum number of Tuples currently allocated for*/
+ uint get_max() const;
+
+ bool get_writeEnabled() const;
+
+ /*Increment n by 1
+ * Warning, automatically calls enableWriteAccess()
+ * returns current TupleList.n after the increment */
+ uint inc_n();
+
+ void print(const char *) const;
+
+ //Variables to allow for direct write access
+ sint *vi_wr; slong *vl_wr; ulong *vul_wr; realType *vr_wr;
+
+ //Variables to allow for direct read access
+ const sint *vi_rd; slong *vl_rd; ulong *vul_rd; realType *vr_rd;
+
+ private:
+ /* storage layed out as: vi[max][mi], vl[max][ml], vul[max][mul],
+ * vr[max][mr] where "tuple" i is given by
+ * (vi[i][0:mi-1],vl[i][0:ml-1],vul[i][0:mul-1],vr[i][0:mr-1]).
+ * only the first n tuples are in use */
+ uint mi,ml,mul,mr;
+ uint n, max;
+ sint *vi; slong *vl; ulong *vul; realType *vr;
+
+ // Used by sort: see .cpp for more details
+ //void sort_bits(uint *work, uint key);
+ void permute(uint *perm, void *work);
+
+ /* last_sorted = the last sorted position in the tuple (if the
+ * TupleList has not been sorted, or has become unsorted--i.e.
+ * by adding a tuple--last_sorted = -1) */
+ int last_sorted;
+ //Whether or not the object is currently allowing direct
+ //write access to the arrays
+ bool writeEnabled;
+
+ typedef uint Index;
+
+ template <typename Value>
+ struct SortData {Value v; Index i;};
+
+#define DIGIT_BITS 8
+#define DIGIT_VALUES (1<<DIGIT_BITS)
+#define DIGIT_MASK ((Value)(DIGIT_VALUES-1))
+#define CEILDIV(a,b) (((a)+(b)-1)/(b))
+#define DIGITS CEILDIV(CHAR_BIT*sizeof(Value),DIGIT_BITS)
+#define VALUE_BITS (DIGIT_BITS*DIGITS)
+#define COUNT_SIZE (DIGITS*DIGIT_VALUES)
+
+ template<class Value>
+ static Value radix_count(const Value *A, const Value *end, Index stride,
+ Index count[DIGITS][DIGIT_VALUES]);
+
+ static void radix_offsets(Index *c);
+
+ template<class Value>
+ static unsigned radix_zeros(Value bitorkey, Index count[DIGITS][DIGIT_VALUES],
+ unsigned *shift, Index **offsets);
+
+ template<class Value>
+ static void radix_index_pass_b(const Value *A, Index n, Index stride,
+ unsigned sh, Index *off, SortData<Value> *out);
+
+ template<class Value>
+ static void radix_index_pass_m(const SortData<Value> *src, const SortData<Value> *end,
+ unsigned sh, Index *off, SortData<Value> *out);
+
+ template<class Value>
+ static void radix_index_pass_e(const SortData<Value> *src, const SortData<Value> *end,
+ unsigned sh, Index *off, Index *out);
+
+ template<class Value>
+ static void radix_index_pass_be(const Value *A, Index n, Index stride,
+ unsigned sh, Index *off, Index *out);
+
+
+ /*------------------------------------------------------------------------------
+
+
+ Radix Sort
+
+ stable; O(n) time
+
+ ----------------------------------------------------------------------------*/
+ template<class Value>
+ static void radix_index_sort(const Value *A, Index n, Index stride,
+ Index *idx, SortData<Value> *work);
+
+ /*------------------------------------------------------------------------------
+
+
+ Merge Sort
+
+ stable; O(n log n) time
+
+ ----------------------------------------------------------------------------*/
+ template<class Value>
+ static void merge_index_sort(const Value *A, const Index An, Index stride,
+ Index *idx, SortData<Value> *work);
+
+ template<class Value>
+ static void index_sort(const Value *A, Index n, Index stride,
+ Index *idx, SortData<Value> *work);
+
+
+#undef DIGIT_BITS
+#undef DIGIT_VALUES
+#undef DIGIT_MASK
+#undef CEILDIV
+#undef DIGITS
+#undef VALUE_BITS
+#undef COUNT_SIZE
+ };
+
+ inline uint TupleList::get_max() const{ return max; }
+ inline uint TupleList::get_n() const{ return n; }
+ inline bool TupleList::get_writeEnabled() const{ return writeEnabled; }
+
+} //namespace
+#endif
+#include <stdlib.h>
diff --git a/src/parallel/Makefile.am b/src/parallel/Makefile.am
index 2718263..b80b85e 100644
--- a/src/parallel/Makefile.am
+++ b/src/parallel/Makefile.am
@@ -34,8 +34,7 @@ if PARALLEL
ReadParallel.hpp \
SharedSetData.cpp \
SharedSetData.hpp \
- gs.cpp \
- TupleList.cpp
+ gs.cpp
# The list of header files which are to be installed
@@ -47,8 +46,7 @@ if PARALLEL
moab/ParallelMergeMesh.hpp \
moab/ProcConfig.hpp \
moab/ParallelData.hpp \
- MBParallelConventions.h \
- moab/TupleList.hpp
+ MBParallelConventions.h
if PARALLEL_HDF5
# libMOABpar_la_LIBADD = $(top_builddir)/mhdf/libmhdf.la
This diff is so big that we needed to truncate the remainder.
https://bitbucket.org/fathomteam/moab/commits/a665d9b23fa9/
Changeset: a665d9b23fa9
Branch: None
User: tautges
Date: 2014-02-09 17:16:34
Summary: Removing reference to install include location, since that may refer to stale files.
Should never have to refer there, since all includes in that location come from within source anyway.
Affected #: 1 file
diff --git a/test/perf/point_location/Makefile.am b/test/perf/point_location/Makefile.am
index 308aec4..8aa6258 100644
--- a/test/perf/point_location/Makefile.am
+++ b/test/perf/point_location/Makefile.am
@@ -2,7 +2,6 @@
AM_CPPFLAGS += -I$(top_srcdir)/src \
-I$(top_srcdir)/src/io/mhdf/include \
-I$(top_builddir)/ \
- -I$(top_builddir)/include \
-I$(top_srcdir)/src/LocalDiscretization/ \
-I$(top_srcdir)/tools/mbcoupler \
-I$(top_srcdir)/src/moab/point_locator/lotte \
https://bitbucket.org/fathomteam/moab/commits/05bb9dabf142/
Changeset: 05bb9dabf142
Branch: None
User: tautges
Date: 2014-02-09 17:16:34
Summary: Adding installation of Version.h.
Technically, I wonder whether this should go in lib/ instead, since this is typically where compiled things go; however, there's a precedent for putting .in-generated files in include, so I'll stick with that.
Affected #: 1 file
diff --git a/src/Makefile.am b/src/Makefile.am
index 0a05543..5772d12 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -200,6 +200,7 @@ nobase_libMOAB_la_include_HEADERS = \
moab/Types.hpp \
moab/UnknownInterface.hpp \
moab/Util.hpp \
+ moab/Version.h \
moab/WriteUtilIface.hpp \
moab/WriterIface.hpp \
MBEntityType.h \
https://bitbucket.org/fathomteam/moab/commits/67c25bc45bc5/
Changeset: 67c25bc45bc5
Branch: None
User: pshriwise
Date: 2014-02-09 17:16:34
Summary: Added number of vertices check to the simple cube test for ReadCGM
Affected #: 1 file
diff --git a/test/io/read_cgm_basic_test.cpp b/test/io/read_cgm_basic_test.cpp
index 3998d01..88f5d0d 100644
--- a/test/io/read_cgm_basic_test.cpp
+++ b/test/io/read_cgm_basic_test.cpp
@@ -55,8 +55,15 @@ void read_cube_test()
std::cout << "Number of Triangles = " << number_of_tris << std::endl;
CHECK_ERR(rval);
+ int number_of_vertices;
+ rval = mb.get_number_entities_by_type(0, MBVERTEX, number_of_vertices);
+ CHECK_ERR(rval);
+
+
if( number_of_tris != 12) rval = MB_FAILURE; CHECK_ERR(rval);
+ if( number_of_vertices !=8) rval = MB_FAILURE; CHECK_ERR(rval);
+
}
int main(int /* argc */, char** /* argv */)
https://bitbucket.org/fathomteam/moab/commits/9d21a558ba6b/
Changeset: 9d21a558ba6b
Branch: None
User: pshriwise
Date: 2014-02-09 17:16:34
Summary: Removed unused lines for the old dum.sat file from read_cgm_basic_tests.cpp
Affected #: 1 file
diff --git a/test/io/read_cgm_basic_test.cpp b/test/io/read_cgm_basic_test.cpp
index 88f5d0d..156e5bd 100644
--- a/test/io/read_cgm_basic_test.cpp
+++ b/test/io/read_cgm_basic_test.cpp
@@ -16,19 +16,6 @@ using namespace moab;
<< __LINE__ << std::endl; \
return A; } } while(false)
-#ifdef MESHDIR
-#ifdef HAVE_OCC_STEP
-static const char input_file[] = STRINGIFY(MESHDIR) "/io/dum.stp";
-#else
-static const char input_file[] = STRINGIFY(MESHDIR) "/io/dum.sat";
-#endif
-#else
-#ifdef HAVE_OCC_STEP
-static const char input_file[] = "dum.stp";
-#else
-static const char input_file[] = "dum.sat";
-#endif
-#endif
#ifdef MESHDIR
static const char input_cube[] = STRINGIFY(MESHDIR) "/io/cube.sat";
https://bitbucket.org/fathomteam/moab/commits/44c59b0c4a8d/
Changeset: 44c59b0c4a8d
Branch: None
User: pshriwise
Date: 2014-02-09 17:35:57
Summary: Updated io test Makefile sources to include new files.
Affected #: 1 file
diff --git a/test/io/Makefile.am b/test/io/Makefile.am
index e467e78..a0b3a97 100644
--- a/test/io/Makefile.am
+++ b/test/io/Makefile.am
@@ -70,8 +70,10 @@ ideas_test_SOURCES = $(srcdir)/../TestUtil.hpp ideas_test.cc
stl_test_SOURCES = $(srcdir)/../TestUtil.hpp stl_test.cc
smf_test_SOURCES = $(srcdir)/../TestUtil.hpp smf_test.cc
ccmio_test_SOURCES = $(srcdir)/../TestUtil.hpp ccmio_test.cpp
-read_cgm_test_SOURCES = $(srcdir)/../TestUtil.hpp read_cgm_test.cpp
-read_cgm_test_CXXFLAGS = $(CGM_CPPFLAGS) $(CXXFLAGS)
+read_cgm__load_test_SOURCES = $(srcdir)/../TestUtil.hpp read_cgm_load_test.cpp
+read_cgm_load_test_CXXFLAGS = $(CGM_CPPFLAGS) $(CXXFLAGS)
+read_cgm_basic_test_SOURCES = $(srcdir)/../TestUtil.hpp read_cgm_basic_test.cpp
+read_cgm_basic_test_CXXFLAGS = $(CGM_CPPFLAGS) $(CXXFLAGS)
readutil_test_SOURCES = $(srcdir)/../TestUtil.hpp readutil_test.cpp
cgns_test_SOURCES=$(srcdir)/../TestUtil.hpp cgns_test.cpp
https://bitbucket.org/fathomteam/moab/commits/e2f51e5ec299/
Changeset: e2f51e5ec299
Branch: None
User: pshriwise
Date: 2014-02-10 00:45:28
Summary: Correction to Makefile for read_cgm_load_test sources.
Affected #: 1 file
diff --git a/test/io/Makefile.am b/test/io/Makefile.am
index a0b3a97..6cf835a 100644
--- a/test/io/Makefile.am
+++ b/test/io/Makefile.am
@@ -70,7 +70,7 @@ ideas_test_SOURCES = $(srcdir)/../TestUtil.hpp ideas_test.cc
stl_test_SOURCES = $(srcdir)/../TestUtil.hpp stl_test.cc
smf_test_SOURCES = $(srcdir)/../TestUtil.hpp smf_test.cc
ccmio_test_SOURCES = $(srcdir)/../TestUtil.hpp ccmio_test.cpp
-read_cgm__load_test_SOURCES = $(srcdir)/../TestUtil.hpp read_cgm_load_test.cpp
+read_cgm_load_test_SOURCES = $(srcdir)/../TestUtil.hpp read_cgm_load_test.cpp
read_cgm_load_test_CXXFLAGS = $(CGM_CPPFLAGS) $(CXXFLAGS)
read_cgm_basic_test_SOURCES = $(srcdir)/../TestUtil.hpp read_cgm_basic_test.cpp
read_cgm_basic_test_CXXFLAGS = $(CGM_CPPFLAGS) $(CXXFLAGS)
https://bitbucket.org/fathomteam/moab/commits/183ef203f60f/
Changeset: 183ef203f60f
Branch: master
User: gonuke
Date: 2014-02-10 21:06:31
Summary: Merge remote-tracking branch 'pshriwise/add_cgm_tests_2'
This PR tests successfully.
Affected #: 5 files
diff --git a/MeshFiles/unittest/io/cube.sat b/MeshFiles/unittest/io/cube.sat
new file mode 100644
index 0000000..58fd508
--- /dev/null
+++ b/MeshFiles/unittest/io/cube.sat
@@ -0,0 +1,117 @@
+1900 0 1 0
+10 Cubit 12.2 17 ACIS 19.0.2 Linux 24 Mon Jan 6 14:16:03 2014
+1 9.9999999999999995e-07 1e-10
+body $1 -1 -1 $-1 $2 $-1 $-1 T -5 -5 -5 5 5 5 #
+integer_attrib-name_attrib-gen-attrib $-1 -1 $-1 $-1 $0 2 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 @8 CUBIT_ID 1 #
+lump $3 -1 -1 $-1 $-1 $4 $0 T -5 -5 -5 5 5 5 #
+integer_attrib-name_attrib-gen-attrib $-1 -1 $-1 $-1 $2 2 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 @8 CUBIT_ID 1 #
+shell $-1 -1 -1 $-1 $-1 $-1 $5 $-1 $2 T -5 -5 -5 5 5 5 #
+face $6 -1 -1 $-1 $7 $8 $4 $-1 $9 forward single T -5 -5 5 5 5 5 F #
+integer_attrib-name_attrib-gen-attrib $-1 -1 $-1 $-1 $5 2 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 @8 CUBIT_ID 1 #
+face $10 -1 -1 $-1 $11 $12 $4 $-1 $13 reversed single T -5 -5 -5 5 5 -5 F #
+loop $-1 -1 -1 $-1 $-1 $14 $5 T -5 -5 5 5 5 5 unknown #
+plane-surface $-1 -1 -1 $-1 0 0 5 0 0 1 1 0 0 forward_v I I I I #
+integer_attrib-name_attrib-gen-attrib $-1 -1 $-1 $-1 $7 2 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 @8 CUBIT_ID 2 #
+face $15 -1 -1 $-1 $16 $17 $4 $-1 $18 reversed single T -5 -5 -5 5 -5 5 F #
+loop $-1 -1 -1 $-1 $-1 $19 $7 T -5 -5 -5 5 5 -5 unknown #
+plane-surface $-1 -1 -1 $-1 0 0 -5 0 0 1 1 0 0 forward_v I I I I #
+coedge $-1 -1 -1 $-1 $20 $21 $22 $23 forward $8 $-1 #
+integer_attrib-name_attrib-gen-attrib $-1 -1 $-1 $-1 $11 2 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 @8 CUBIT_ID 3 #
+face $24 -1 -1 $-1 $25 $26 $4 $-1 $27 reversed single T -5 -5 -5 -5 5 5 F #
+loop $-1 -1 -1 $-1 $-1 $28 $11 T -5 -5 -5 5 -5 5 unknown #
+plane-surface $-1 -1 -1 $-1 0 -5 0 0 1 -0 -0 0 1 forward_v I I I I #
+coedge $-1 -1 -1 $-1 $29 $30 $31 $32 forward $12 $-1 #
+coedge $-1 -1 -1 $-1 $33 $14 $34 $35 forward $8 $-1 #
+coedge $-1 -1 -1 $-1 $14 $33 $36 $37 forward $8 $-1 #
+coedge $-1 -1 -1 $-1 $38 $39 $14 $23 reversed $40 $-1 #
+edge $41 -1 -1 $-1 $42 -5 $43 5 $22 $44 forward @7 unknown T 5 -5 5 5 5 5 #
+integer_attrib-name_attrib-gen-attrib $-1 -1 $-1 $-1 $16 2 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 @8 CUBIT_ID 4 #
+face $45 -1 -1 $-1 $46 $47 $4 $-1 $48 reversed single T -5 5 -5 5 5 5 F #
+loop $-1 -1 -1 $-1 $-1 $49 $16 T -5 -5 -5 -5 5 5 unknown #
+plane-surface $-1 -1 -1 $-1 -5 0 0 1 0 0 0 0 -1 forward_v I I I I #
+coedge $-1 -1 -1 $-1 $50 $36 $51 $52 forward $17 $-1 #
+coedge $-1 -1 -1 $-1 $53 $19 $50 $54 forward $12 $-1 #
+coedge $-1 -1 -1 $-1 $19 $53 $55 $56 forward $12 $-1 #
+coedge $-1 -1 -1 $-1 $39 $38 $19 $32 reversed $40 $-1 #
+edge $57 -1 -1 $-1 $58 -5 $59 5 $31 $60 forward @7 unknown T 5 -5 -5 5 5 -5 #
+coedge $-1 -1 -1 $-1 $21 $20 $61 $62 forward $8 $-1 #
+coedge $-1 -1 -1 $-1 $63 $64 $20 $35 reversed $47 $-1 #
+edge $65 -1 -1 $-1 $43 -5 $66 5 $34 $67 forward @7 unknown T -5 5 5 5 5 5 #
+coedge $-1 -1 -1 $-1 $28 $68 $21 $37 reversed $17 $-1 #
+edge $69 -1 -1 $-1 $70 -5 $42 5 $36 $71 forward @7 unknown T -5 -5 5 5 -5 5 #
+coedge $-1 -1 -1 $-1 $31 $22 $68 $72 forward $40 $-1 #
+coedge $-1 -1 -1 $-1 $22 $31 $63 $73 reversed $40 $-1 #
+loop $-1 -1 -1 $-1 $-1 $38 $46 T 5 -5 -5 5 5 5 unknown #
+integer_attrib-name_attrib-gen-attrib $-1 -1 $-1 $-1 $23 2 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 @8 CUBIT_ID 1 #
+vertex $74 -1 -1 $-1 $23 $75 #
+vertex $76 -1 -1 $-1 $23 $77 #
+straight-curve $-1 -1 -1 $-1 5 0 5 0 1 0 I I #
+integer_attrib-name_attrib-gen-attrib $-1 -1 $-1 $-1 $25 2 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 @8 CUBIT_ID 5 #
+face $78 -1 -1 $-1 $-1 $40 $4 $-1 $79 reversed single T 5 -5 -5 5 5 5 F #
+loop $-1 -1 -1 $-1 $-1 $63 $25 T -5 5 -5 5 5 5 unknown #
+plane-surface $-1 -1 -1 $-1 0 5 0 0 -1 0 0 0 -1 forward_v I I I I #
+coedge $-1 -1 -1 $-1 $80 $61 $64 $81 forward $26 $-1 #
+coedge $-1 -1 -1 $-1 $68 $28 $29 $54 reversed $17 $-1 #
+coedge $-1 -1 -1 $-1 $61 $80 $28 $52 reversed $26 $-1 #
+edge $82 -1 -1 $-1 $70 -5 $83 5 $51 $84 forward @7 unknown T -5 -5 -5 -5 -5 5 #
+coedge $-1 -1 -1 $-1 $30 $29 $80 $85 forward $12 $-1 #
+edge $86 -1 -1 $-1 $59 -5 $83 5 $50 $87 forward @7 unknown T -5 -5 -5 5 -5 -5 #
+coedge $-1 -1 -1 $-1 $64 $63 $30 $56 reversed $47 $-1 #
+edge $88 -1 -1 $-1 $89 -5 $58 5 $55 $90 forward @7 unknown T -5 5 -5 5 5 -5 #
+integer_attrib-name_attrib-gen-attrib $-1 -1 $-1 $-1 $32 2 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 @8 CUBIT_ID 5 #
+vertex $91 -1 -1 $-1 $32 $92 #
+vertex $93 -1 -1 $-1 $72 $94 #
+straight-curve $-1 -1 -1 $-1 5 0 -5 0 -1 0 I I #
+coedge $-1 -1 -1 $-1 $49 $51 $33 $62 reversed $26 $-1 #
+edge $95 -1 -1 $-1 $66 -5 $70 5 $61 $96 forward @7 unknown T -5 -5 5 -5 5 5 #
+coedge $-1 -1 -1 $-1 $55 $34 $39 $73 forward $47 $-1 #
+coedge $-1 -1 -1 $-1 $34 $55 $49 $81 reversed $47 $-1 #
+integer_attrib-name_attrib-gen-attrib $-1 -1 $-1 $-1 $35 2 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 @8 CUBIT_ID 2 #
+vertex $97 -1 -1 $-1 $35 $98 #
+straight-curve $-1 -1 -1 $-1 0 5 5 -1 0 0 I I #
+coedge $-1 -1 -1 $-1 $36 $50 $38 $72 reversed $17 $-1 #
+integer_attrib-name_attrib-gen-attrib $-1 -1 $-1 $-1 $37 2 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 @8 CUBIT_ID 4 #
+vertex $99 -1 -1 $-1 $62 $100 #
+straight-curve $-1 -1 -1 $-1 0 -5 5 1 0 0 I I #
+edge $101 -1 -1 $-1 $42 -5 $59 5 $38 $102 forward @7 unknown T 5 -5 -5 5 -5 5 #
+edge $103 -1 -1 $-1 $43 -5 $58 5 $39 $104 forward @7 unknown T 5 5 -5 5 5 5 #
+integer_attrib-name_attrib-gen-attrib $-1 -1 $-1 $-1 $42 2 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 @8 CUBIT_ID 1 #
+point $-1 -1 -1 $-1 5 -5 5 #
+integer_attrib-name_attrib-gen-attrib $-1 -1 $-1 $-1 $43 2 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 @8 CUBIT_ID 2 #
+point $-1 -1 -1 $-1 5 5 5 #
+integer_attrib-name_attrib-gen-attrib $-1 -1 $-1 $-1 $46 2 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 @8 CUBIT_ID 6 #
+plane-surface $-1 -1 -1 $-1 5 0 0 -1 0 0 0 -0 1 forward_v I I I I #
+coedge $-1 -1 -1 $-1 $51 $49 $53 $85 reversed $26 $-1 #
+edge $105 -1 -1 $-1 $66 -5 $89 5 $64 $106 forward @7 unknown T -5 5 -5 -5 5 5 #
+integer_attrib-name_attrib-gen-attrib $-1 -1 $-1 $-1 $52 2 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 @8 CUBIT_ID 9 #
+vertex $107 -1 -1 $-1 $85 $108 #
+straight-curve $-1 -1 -1 $-1 -5 -5 0 0 0 -1 I I #
+edge $109 -1 -1 $-1 $83 -5 $89 5 $80 $110 forward @7 unknown T -5 -5 -5 -5 5 -5 #
+integer_attrib-name_attrib-gen-attrib $-1 -1 $-1 $-1 $54 2 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 @8 CUBIT_ID 6 #
+straight-curve $-1 -1 -1 $-1 0 -5 -5 -1 0 0 I I #
+integer_attrib-name_attrib-gen-attrib $-1 -1 $-1 $-1 $56 2 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 @8 CUBIT_ID 8 #
+vertex $111 -1 -1 $-1 $56 $112 #
+straight-curve $-1 -1 -1 $-1 0 5 -5 1 0 0 I I #
+integer_attrib-name_attrib-gen-attrib $-1 -1 $-1 $-1 $58 2 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 @8 CUBIT_ID 5 #
+point $-1 -1 -1 $-1 5 5 -5 #
+integer_attrib-name_attrib-gen-attrib $-1 -1 $-1 $-1 $59 2 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 @8 CUBIT_ID 6 #
+point $-1 -1 -1 $-1 5 -5 -5 #
+integer_attrib-name_attrib-gen-attrib $-1 -1 $-1 $-1 $62 2 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 @8 CUBIT_ID 3 #
+straight-curve $-1 -1 -1 $-1 -5 0 5 0 -1 0 I I #
+integer_attrib-name_attrib-gen-attrib $-1 -1 $-1 $-1 $66 2 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 @8 CUBIT_ID 3 #
+point $-1 -1 -1 $-1 -5 5 5 #
+integer_attrib-name_attrib-gen-attrib $-1 -1 $-1 $-1 $70 2 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 @8 CUBIT_ID 4 #
+point $-1 -1 -1 $-1 -5 -5 5 #
+integer_attrib-name_attrib-gen-attrib $-1 -1 $-1 $-1 $72 2 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 @8 CUBIT_ID 10 #
+straight-curve $-1 -1 -1 $-1 5 -5 0 0 0 -1 I I #
+integer_attrib-name_attrib-gen-attrib $-1 -1 $-1 $-1 $73 2 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 @8 CUBIT_ID 12 #
+straight-curve $-1 -1 -1 $-1 5 5 0 0 0 -1 I I #
+integer_attrib-name_attrib-gen-attrib $-1 -1 $-1 $-1 $81 2 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 @8 CUBIT_ID 11 #
+straight-curve $-1 -1 -1 $-1 -5 5 0 0 0 -1 I I #
+integer_attrib-name_attrib-gen-attrib $-1 -1 $-1 $-1 $83 2 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 @8 CUBIT_ID 7 #
+point $-1 -1 -1 $-1 -5 -5 -5 #
+integer_attrib-name_attrib-gen-attrib $-1 -1 $-1 $-1 $85 2 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 @8 CUBIT_ID 7 #
+straight-curve $-1 -1 -1 $-1 -5 0 -5 0 1 0 I I #
+integer_attrib-name_attrib-gen-attrib $-1 -1 $-1 $-1 $89 2 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 @8 CUBIT_ID 8 #
+point $-1 -1 -1 $-1 -5 5 -5 #
+End-of-ACIS-data
\ No newline at end of file
diff --git a/test/io/Makefile.am b/test/io/Makefile.am
index c776727..6cf835a 100644
--- a/test/io/Makefile.am
+++ b/test/io/Makefile.am
@@ -25,7 +25,7 @@ else
endif
if HAVE_CGM
- CGM_TEST = read_cgm_test
+ CGM_TEST = read_cgm_load_test read_cgm_basic_test
else
CGM_TEST =
endif
@@ -70,8 +70,10 @@ ideas_test_SOURCES = $(srcdir)/../TestUtil.hpp ideas_test.cc
stl_test_SOURCES = $(srcdir)/../TestUtil.hpp stl_test.cc
smf_test_SOURCES = $(srcdir)/../TestUtil.hpp smf_test.cc
ccmio_test_SOURCES = $(srcdir)/../TestUtil.hpp ccmio_test.cpp
-read_cgm_test_SOURCES = $(srcdir)/../TestUtil.hpp read_cgm_test.cpp
-read_cgm_test_CXXFLAGS = $(CGM_CPPFLAGS) $(CXXFLAGS)
+read_cgm_load_test_SOURCES = $(srcdir)/../TestUtil.hpp read_cgm_load_test.cpp
+read_cgm_load_test_CXXFLAGS = $(CGM_CPPFLAGS) $(CXXFLAGS)
+read_cgm_basic_test_SOURCES = $(srcdir)/../TestUtil.hpp read_cgm_basic_test.cpp
+read_cgm_basic_test_CXXFLAGS = $(CGM_CPPFLAGS) $(CXXFLAGS)
readutil_test_SOURCES = $(srcdir)/../TestUtil.hpp readutil_test.cpp
cgns_test_SOURCES=$(srcdir)/../TestUtil.hpp cgns_test.cpp
diff --git a/test/io/read_cgm_basic_test.cpp b/test/io/read_cgm_basic_test.cpp
new file mode 100644
index 0000000..156e5bd
--- /dev/null
+++ b/test/io/read_cgm_basic_test.cpp
@@ -0,0 +1,63 @@
+
+#include <iostream>
+#include "moab/Interface.hpp"
+#ifndef IS_BUILDING_MB
+#define IS_BUILDING_MB
+#endif
+#include "TestUtil.hpp"
+#include "Internals.hpp"
+#include "moab/Core.hpp"
+#include "MBTagConventions.hpp"
+
+using namespace moab;
+
+#define CHKERR(A) do { if (MB_SUCCESS != (A)) { \
+ std::cerr << "Failure (error code " << (A) << ") at " __FILE__ ":" \
+ << __LINE__ << std::endl; \
+ return A; } } while(false)
+
+
+#ifdef MESHDIR
+static const char input_cube[] = STRINGIFY(MESHDIR) "/io/cube.sat";
+#else
+static const char input_cube[] = "/io/cube.sat";
+#endif
+
+
+void read_cube_test()
+{
+ Core mb;
+
+ ErrorCode rval = mb.load_file(input_cube); CHECK_ERR(rval);
+
+ Tag geom_tag;
+
+ rval = mb.tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER,
+ geom_tag, MB_TAG_SPARSE|MB_TAG_CREAT);
+ CHECK_ERR(rval);
+
+ int number_of_tris;
+
+ rval = mb.get_number_entities_by_type(0, MBTRI , number_of_tris);
+ std::cout << "Number of Triangles = " << number_of_tris << std::endl;
+ CHECK_ERR(rval);
+
+ int number_of_vertices;
+ rval = mb.get_number_entities_by_type(0, MBVERTEX, number_of_vertices);
+ CHECK_ERR(rval);
+
+
+ if( number_of_tris != 12) rval = MB_FAILURE; CHECK_ERR(rval);
+
+ if( number_of_vertices !=8) rval = MB_FAILURE; CHECK_ERR(rval);
+
+}
+
+int main(int /* argc */, char** /* argv */)
+{
+ int result = 0;
+
+ result += RUN_TEST( read_cube_test );
+
+ return result;
+}
diff --git a/test/io/read_cgm_load_test.cpp b/test/io/read_cgm_load_test.cpp
new file mode 100644
index 0000000..5377220
--- /dev/null
+++ b/test/io/read_cgm_load_test.cpp
@@ -0,0 +1,49 @@
+
+#include <iostream>
+#include "moab/Interface.hpp"
+#ifndef IS_BUILDING_MB
+#define IS_BUILDING_MB
+#endif
+#include "TestUtil.hpp"
+#include "Internals.hpp"
+#include "moab/Core.hpp"
+
+using namespace moab;
+
+#define CHKERR(A) do { if (MB_SUCCESS != (A)) { \
+ std::cerr << "Failure (error code " << (A) << ") at " __FILE__ ":" \
+ << __LINE__ << std::endl; \
+ return A; } } while(false)
+
+#ifdef MESHDIR
+#ifdef HAVE_OCC_STEP
+static const char input_file[] = STRINGIFY(MESHDIR) "/io/dum.stp";
+#else
+static const char input_file[] = STRINGIFY(MESHDIR) "/io/dum.sat";
+#endif
+#else
+#ifdef HAVE_OCC_STEP
+static const char input_file[] = "dum.stp";
+#else
+static const char input_file[] = "dum.sat";
+#endif
+#endif
+
+void read_multiple_test()
+{
+ Core mb;
+
+ ErrorCode rval = mb.load_file(input_file);
+ CHECK_ERR(rval);
+ // second load
+ rval = mb.load_file(input_file);
+ CHECK_ERR(rval);
+
+}
+
+int main(int /* argc */, char** /* argv */)
+{
+ int result = RUN_TEST( read_multiple_test );
+
+ return result;
+}
diff --git a/test/io/read_cgm_test.cpp b/test/io/read_cgm_test.cpp
deleted file mode 100644
index 5377220..0000000
--- a/test/io/read_cgm_test.cpp
+++ /dev/null
@@ -1,49 +0,0 @@
-
-#include <iostream>
-#include "moab/Interface.hpp"
-#ifndef IS_BUILDING_MB
-#define IS_BUILDING_MB
-#endif
-#include "TestUtil.hpp"
-#include "Internals.hpp"
-#include "moab/Core.hpp"
-
-using namespace moab;
-
-#define CHKERR(A) do { if (MB_SUCCESS != (A)) { \
- std::cerr << "Failure (error code " << (A) << ") at " __FILE__ ":" \
- << __LINE__ << std::endl; \
- return A; } } while(false)
-
-#ifdef MESHDIR
-#ifdef HAVE_OCC_STEP
-static const char input_file[] = STRINGIFY(MESHDIR) "/io/dum.stp";
-#else
-static const char input_file[] = STRINGIFY(MESHDIR) "/io/dum.sat";
-#endif
-#else
-#ifdef HAVE_OCC_STEP
-static const char input_file[] = "dum.stp";
-#else
-static const char input_file[] = "dum.sat";
-#endif
-#endif
-
-void read_multiple_test()
-{
- Core mb;
-
- ErrorCode rval = mb.load_file(input_file);
- CHECK_ERR(rval);
- // second load
- rval = mb.load_file(input_file);
- CHECK_ERR(rval);
-
-}
-
-int main(int /* argc */, char** /* argv */)
-{
- int result = RUN_TEST( read_multiple_test );
-
- return result;
-}
Repository URL: https://bitbucket.org/fathomteam/moab/
--
This is a commit notification from bitbucket.org. You are receiving
this because you have the service enabled, addressing the recipient of
this email.
More information about the moab-dev
mailing list