[MOAB-dev] commit/MOAB: rajeeja: Important changes needed in Version4.6
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Fri May 17 11:57:54 CDT 2013
1 new commit in MOAB:
https://bitbucket.org/fathomteam/moab/commits/c743ac29fc66/
Changeset: c743ac29fc66
Branch: Version4.6
User: rajeeja
Date: 2013-05-17 18:56:10
Summary: Important changes needed in Version4.6
Merge branch 'master' into Version4.6
Affected #: 129 files
diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..6d29b5e
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,187 @@
+MOABConfig.cmake
+moab.files
+config/test-driver
+moab.config
+moab.includes
+moab.creator*
+*.ccm
+*.cub
+aclocal.m4
+autom4te.cache/
+config.h
+config.h.in
+config.log
+config.lt
+config.status
+config/config.guess
+config/config.sub
+config/depcomp
+config/install-sh
+config/libtool.m4
+config/ltmain.sh
+config/ltoptions.m4
+config/ltsugar.m4
+config/ltversion.m4
+config/lt~obsolete.m4
+config/missing
+configure
+doc/config.tex
+doc/dev.dox
+doc/user.dox
+doc/user/*
+examples/examples.make
+itaps/iBase_f.h
+itaps/igeom/FBiGeom-Defs.inc
+itaps/imesh/iMesh-Defs.inc
+itaps/imesh/iMeshP_extensions_protos.h
+itaps/imesh/iMeshP_protos.h
+itaps/imesh/iMesh_extensions_protos.h
+itaps/imesh/iMesh_protos.h
+libtool
+moab.make
+src/FCDefs.h
+src/MBCN_protos.h
+src/MOAB_FCDefs.h
+src/moab/EntityHandle.hpp
+src/moab/Version.h
+src/moab/stamp-h2
+src/moab/stamp-h3
+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/coords_connect_iterate
+test/cropvol_test
+test/dual/dual_test
+test/file_options_test
+test/geom_util_test
+test/gttool_test
+test/h5file/dump_sets
+test/h5file/h5legacy
+test/h5file/h5partial
+test/h5file/h5portable
+test/h5file/h5regression
+test/h5file/h5sets_test
+test/h5file/h5test
+test/h5file/h5varlen
+test/*.g
+test/homxform_test
+test/io/cub_file_test
+test/io/exodus_test
+test/io/gmsh_test
+test/io/ideas_test
+test/io/*.g
+test/io/nastran_test
+test/io/read_cgm_test
+test/io/read_nc
+test/io/read_ucd_nc
+test/io/readutil_test
+test/io/smf_test
+test/io/stl_test
+test/io/tqdcfr
+test/io/vtk_test
+test/kd_tree_test
+test/kd_tree_time
+test/kd_tree_tool
+test/mbcn_test
+test/mbfacet_test
+test/mbground_test
+test/mesh_set_test
+test/moab_test
+test/obb/obb_test
+test/obb/obb_time
+test/obb/obb_tree_tool
+test/obb_test
+test/oldinc/test_oldinc
+test/parallel/*.h5m
+test/parallel/*.vtk
+test/parallel/mbparallelcomm_test
+test/parallel/mhdf_parallel
+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/parmerge
+test/parallel/partcheck
+test/parallel/pcomm_serial
+test/parallel/pcomm_unit
+test/parallel/read_nc_par
+test/parallel/scdpart
+test/parallel/scdtest
+test/parallel/structured3
+test/parallel/uber_parallel_test
+test/parallel/ucdtrvpart
+test/perf/adj_time
+test/perf/perf
+test/perf/perftool
+test/perf/seqperf
+test/perf/tstt_perf_binding
+test/range_test
+test/reorder_test
+test/scdseq_test
+test/seq_man_test
+test/tag_test
+test/test_adj
+test/test_prog_opt
+test/var_len_test
+test/var_len_test_no_template
+test/xform_test
+tools/dagmc/pt_vol_test
+tools/dagmc/ray_fire_test
+tools/dagmc/test_geom
+tools/dagmc/update_coords
+tools/mbcoupler/*.h5m
+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/lagr.h5m
+tools/mbcslam/spec_visu_test
+tools/mbcslam/spectral.vtk
+tools/mbcslam/spherical_area_test
+.project
+.cproject
+examples/*.h5m
+examples/ReduceExchangeTags
+
diff --git a/Makefile.am b/Makefile.am
index 33a99cf..1c91622 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -5,7 +5,7 @@ AUTOMAKE_OPTIONS = foreign
ACLOCAL_AMFLAGS = -I config
# Subdirectories to build
-SUBDIRS = src MeshFiles $(itaps_dir) tools test doc examples
+SUBDIRS = src MeshFiles $(itaps_dir) tools test doc
if ENABLE_igeom
itaps_dir_igeom = itaps
diff --git a/MeshFiles/unittest/disk.h5m b/MeshFiles/unittest/disk.h5m
new file mode 100644
index 0000000..ef19883
Binary files /dev/null and b/MeshFiles/unittest/disk.h5m differ
diff --git a/RELEASE_NOTES b/RELEASE_NOTES
index ca26bef..d10c95e 100644
--- a/RELEASE_NOTES
+++ b/RELEASE_NOTES
@@ -1,3 +1,5 @@
+Version 4.7pre:
+
Version 4.6:
* Removed deprecated functions from the Interface: (some overloaded variants of) query_interface, release_interface,
tag_get_data, tag_set_data, tag_get_size.
diff --git a/config/compiler.m4 b/config/compiler.m4
index ce00d10..18df411 100644
--- a/config/compiler.m4
+++ b/config/compiler.m4
@@ -126,6 +126,8 @@ fi
if test "xno" != "x$CHECK_FC"; then
AC_PROG_FC
AC_PROG_F77
+ AC_F77_LIBRARY_LDFLAGS
+ AC_FC_LIBRARY_LDFLAGS
fi
]) # FATHOM_CHECK_COMPILERS
diff --git a/configure.ac b/configure.ac
index 7494691..1524b06 100644
--- a/configure.ac
+++ b/configure.ac
@@ -1,7 +1,7 @@
################################################################################
# Standard Stuff
################################################################################
-AC_INIT(MOAB, 4.6.0)
+AC_INIT(MOAB, 4.7.0pre)
AC_CONFIG_SRCDIR([moab.make.in])
AC_CONFIG_SRCDIR([MOABConfig.cmake.in])
AC_CONFIG_MACRO_DIR([config])
@@ -67,6 +67,14 @@ if test "xyes" = "x$ENABLE_FORTRAN" && test "x" != "x$FC"; then
AC_FC_WRAPPERS
fi
+if test "xyes" = "x$ENABLE_FORTRAN" && test "x" != "x$FLIBS"; then
+ LIBS="$LIBS $FLIBS"
+fi
+
+if test "xyes" = "x$ENABLE_FORTRAN" && test "x" != "x$FCLIBS"; then
+ LIBS="$LIBS $FCLIBS"
+fi
+
################################################################################
# Check for need for extra flags to support cray pointers
################################################################################
@@ -1163,6 +1171,8 @@ AC_SUBST([EXPORT_LTFLAGS])
AC_SUBST([EXPORT_LDFLAGS])
AC_SUBST([AM_CXXFLAGS])
AC_SUBST([AM_CFLAGS])
+ AC_SUBST([AM_FCFLAGS])
+ AC_SUBST([AM_FFLAGS])
AC_SUBST([DISTCHECK_CONFIGURE_FLAGS])
AC_ARG_VAR([FC], [FORTRAN compiler command])
@@ -1208,10 +1218,6 @@ AC_CONFIG_FILES([Makefile
doc/user.dox
doc/dev.dox
doc/config.tex
- examples/Makefile
- examples/examples.make
- examples/simple/makefile
- examples/itaps/Makefile
MeshFiles/Makefile
MeshFiles/unittest/Makefile
MeshFiles/unittest/io/Makefile
diff --git a/doc/user.dox.in b/doc/user.dox.in
index 703ce6b..4fed19e 100644
--- a/doc/user.dox.in
+++ b/doc/user.dox.in
@@ -52,7 +52,7 @@ EXTRACT_ALL = YES
# If the EXTRACT_PRIVATE tag is set to YES all private members of a class
# will be included in the documentation.
-EXTRACT_PRIVATE = NO
+EXTRACT_PRIVATE = YES
# If the EXTRACT_STATIC tag is set to YES all static members of a file
# will be included in the documentation.
@@ -112,7 +112,7 @@ STRIP_FROM_PATH =
# to NO (the default) then the documentation will be excluded.
# Set it to YES to include the internal documentation.
-INTERNAL_DOCS = NO
+INTERNAL_DOCS = YES
# If the CLASS_DIAGRAMS tag is set to YES (the default) Doxygen will
# generate a class diagram (in Html and LaTeX) for classes with base or
@@ -123,12 +123,12 @@ CLASS_DIAGRAMS = YES
# If the SOURCE_BROWSER tag is set to YES then a list of source files will
# be generated. Documented entities will be cross-referenced with these sources.
-SOURCE_BROWSER = NO
+SOURCE_BROWSER = YES
# Setting the INLINE_SOURCES tag to YES will include the body
# of functions and classes directly in the documentation.
-INLINE_SOURCES = NO
+INLINE_SOURCES = YES
# Setting the STRIP_CODE_COMMENTS tag to YES (the default) will instruct
# doxygen to hide any special comment blocks from generated source code
@@ -199,12 +199,12 @@ SORT_MEMBER_DOCS = YES
# member in the group (if any) for the other members of the group. By default
# all members of a group must be documented explicitly.
-DISTRIBUTE_GROUP_DOC = NO
+DISTRIBUTE_GROUP_DOC = YES
# The TAB_SIZE tag can be used to set the number of spaces in a tab.
# Doxygen uses this value to replace tabs by spaces in code fragments.
-TAB_SIZE = 8
+TAB_SIZE = 4
# The GENERATE_TODOLIST tag can be used to enable (YES) or
# disable (NO) the todo list. This list is created by putting \todo
@@ -304,14 +304,14 @@ WARN_LOGFILE =
# directories like "/usr/src/myproject". Separate the files or directories
# with spaces.
-INPUT = @top_srcdir@/src/moab \
+INPUT = @top_srcdir@/src @top_srcdir@/src/moab \
@top_srcdir@/src/parallel/moab \
@top_srcdir@/src/MBTagConventions.hpp \
@top_srcdir@/src/MBCN.h \
@top_srcdir@/src/MBEntityType.h \
@top_srcdir@/src/parallel/MBParallelConventions.h \
- @top_srcdir@/tools/dagmc/DagMC.hpp \
- @top_srcdir@/itaps @top_srcdir@/itaps/imesh
+ @top_srcdir@/tools/ @top_srcdir@/tools/mbcoupler @top_srcdir@/tools/dagmc/DagMC.hpp \
+ @top_srcdir@/itaps @top_srcdir@/itaps/imesh @top_srcdir@/examples
# If the value of the INPUT tag contains directories, you can use the
# FILE_PATTERNS tag to specify one or more wildcard pattern (like *.cpp
@@ -339,11 +339,11 @@ EXCLUDE =
EXCLUDE_PATTERNS =
-# The EXAMPLE_PATH tag can be used to specify one or more files or
+# The EXAMPLE_PATH tag can be used to specify one or more files or
# directories that contain example code fragments that are included (see
# the \include command).
-EXAMPLE_PATH =
+EXAMPLE_PATH = @top_srcdir@/examples/
# If the value of the EXAMPLE_PATH tag contains directories, you can use the
# EXAMPLE_PATTERNS tag to specify one or more wildcard pattern (like *.cpp
@@ -458,7 +458,7 @@ BINARY_TOC = NO
# The TOC_EXPAND flag can be set to YES to add extra items for group members
# to the contents of the Html help documentation and to the tree view.
-TOC_EXPAND = NO
+TOC_EXPAND = YES
# The DISABLE_INDEX tag can be used to turn on/off the condensed index at
# top of each HTML page. The value NO (the default) enables the index and
@@ -647,7 +647,7 @@ MACRO_EXPANSION = NO
# then the macro expansion is limited to the macros specified with the
# PREDEFINED and EXPAND_AS_PREDEFINED tags.
-EXPAND_ONLY_PREDEF = NO
+EXPAND_ONLY_PREDEF = YES
# If the SEARCH_INCLUDES tag is set to YES (the default) the includes files
# in the INCLUDE_PATH (see below) will be search if a #include is found.
@@ -808,7 +808,7 @@ DOT_CLEANUP = YES
# The SEARCHENGINE tag specifies whether or not a search engine should be
# used. If set to NO the values of all tags below this one will be ignored.
-SEARCHENGINE = NO
+SEARCHENGINE = YES
# The CGI_NAME tag should be the name of the CGI script that
# starts the search engine (doxysearch) with the correct parameters.
diff --git a/examples/A.1.ele b/examples/A.1.ele
deleted file mode 100644
index d742678..0000000
--- a/examples/A.1.ele
+++ /dev/null
@@ -1,31 +0,0 @@
-29 3 0
- 1 29 2 1
- 2 2 29 23
- 3 25 24 23
- 4 23 22 2
- 5 25 23 29
- 6 2 22 3
- 7 3 21 4
- 8 21 3 22
- 9 4 21 20
- 10 5 4 26
- 11 19 26 4
- 12 26 19 18
- 13 19 4 20
- 14 5 26 28
- 15 12 14 13
- 16 14 12 11
- 17 11 10 9
- 18 8 14 9
- 19 8 15 14
- 20 9 14 11
- 21 6 27 7
- 22 26 18 27
- 23 5 28 6
- 24 27 18 7
- 25 28 27 6
- 26 15 7 16
- 27 7 15 8
- 28 17 7 18
- 29 7 17 16
-# Generated by triangle A.poly
diff --git a/examples/A.1.node b/examples/A.1.node
deleted file mode 100644
index 31d09f1..0000000
--- a/examples/A.1.node
+++ /dev/null
@@ -1,31 +0,0 @@
-29 2 1 1
- 1 0.20000000000000001 -0.77639999999999998 -0.56999999999999995 1
- 2 0.22 -0.7732 -0.55000000000000004 1
- 3 0.24560000000000001 -0.75639999999999996 -0.51000000000000001 1
- 4 0.27760000000000001 -0.70199999999999996 -0.53000000000000003 1
- 5 0.48880000000000001 -0.20760000000000001 0.28000000000000003 1
- 6 0.50480000000000003 -0.20760000000000001 0.29999999999999999 1
- 7 0.74080000000000001 -0.73960000000000004 0 1
- 8 0.75600000000000001 -0.76119999999999999 -0.01 1
- 9 0.77439999999999998 -0.77239999999999998 0 1
- 10 0.80000000000000004 -0.77639999999999998 0.02 1
- 11 0.80000000000000004 -0.79239999999999999 0.01 1
- 12 0.57920000000000005 -0.79239999999999999 -0.20999999999999999 1
- 13 0.57920000000000005 -0.77639999999999998 -0.20000000000000001 1
- 14 0.62160000000000004 -0.77159999999999995 -0.14999999999999999 1
- 15 0.63360000000000005 -0.76280000000000003 -0.13 1
- 16 0.63919999999999999 -0.74439999999999995 -0.10000000000000001 1
- 17 0.62080000000000002 -0.68440000000000001 -0.059999999999999998 1
- 18 0.58720000000000006 -0.60440000000000005 -0.01 1
- 19 0.36080000000000001 -0.60440000000000005 -0.23999999999999999 1
- 20 0.31919999999999998 -0.70679999999999998 -0.39000000000000001 1
- 21 0.312 -0.73960000000000004 -0.42999999999999999 1
- 22 0.31840000000000002 -0.76119999999999999 -0.44 1
- 23 0.33439999999999998 -0.77159999999999995 -0.44 1
- 24 0.37119999999999997 -0.77639999999999998 -0.40999999999999998 1
- 25 0.37119999999999997 -0.79239999999999999 -0.41999999999999998 1
- 26 0.37440000000000001 -0.56999999999999995 -0.20000000000000001 1
- 27 0.57440000000000002 -0.56999999999999995 0 1
- 28 0.47360000000000002 -0.33079999999999998 0.14000000000000001 1
- 29 0.20000000000000001 -0.79239999999999999 -0.58999999999999997 1
-# Generated by triangle A.poly
diff --git a/examples/DirectAccessNoHoles.cpp b/examples/DirectAccessNoHoles.cpp
new file mode 100644
index 0000000..340efae
--- /dev/null
+++ b/examples/DirectAccessNoHoles.cpp
@@ -0,0 +1,186 @@
+/** @example DirectAccessNoHoles.cpp \n
+ * \brief Use direct access to MOAB data to avoid calling through API \n
+ *
+ * This example creates a 1d row of quad elements, such that all quad and vertex handles
+ * are contiguous in the handle space and in the database. Then it shows how to get access
+ * to pointers to MOAB-native data for vertex coordinates, quad connectivity, tag storage,
+ * and vertex to quad adjacency lists. This allows applications to access this data directly
+ * without going through MOAB's API. In cases where the mesh is not changing (or only mesh
+ * vertices are moving), this can save significant execution time in applications.
+ *
+ * ----------------------
+ * | | | |
+ * | | | | ...
+ * | | | |
+ * ----------------------
+ *
+ * -# Initialize MOAB \n
+ * -# Create a quad mesh, as depicted above
+ * -# Create 2 dense tags (tag1, tag2) for avg position to assign to quads, and # verts per quad (tag3)
+ * -# Get connectivity, coordinate, tag1 iterators
+ * -# Iterate through quads, computing midpoint based on vertex positions, set on quad-based tag1
+ * -# Iterate through vertices, summing positions into tag2 on connected quads and incrementing vertex count
+ * -# Iterate through quads, normalizing tag2 by vertex count and comparing values of tag1 and tag2
+ *
+ * <b>To compile</b>: \n
+ * make DirectAccessNoHoles MOAB_DIR=<installdir> \n
+ * <b>To run</b>: ./DirectAccessNoHoles [-nquads <# quads>]\n
+ *
+ */
+
+#include "moab/Core.hpp"
+#include "moab/ProgOptions.hpp"
+#include "moab/ReadUtilIface.hpp"
+#include <map>
+#include <iostream>
+#include <assert.h>
+
+// Error routines for use with MOAB API
+#define CHKERR(CODE, MSG) do {if (MB_SUCCESS != (CODE)) {std::string errstr; mbImpl->get_last_error(errstr); \
+ std::cerr << errstr << std::endl; std::cerr << MSG << std::endl; return CODE;}} while(false)
+
+using namespace moab;
+
+ErrorCode create_mesh_no_holes(Interface *mbImpl, int nquads);
+
+int main(int argc, char **argv)
+{
+ // get MOAB instance
+ Interface *mbImpl = new Core;
+ if (NULL == mbImpl) return 1;
+
+ int nquads = 1000;
+
+ // parse options
+ ProgOptions opts;
+ opts.addOpt<int>(std::string("nquads,n"), std::string("Number of quads in the mesh (default = 1000"), &nquads);
+ opts.parseCommandLine(argc, argv);
+
+ // create simple structured mesh with hole, but using unstructured representation
+ ErrorCode rval = create_mesh_no_holes(mbImpl, nquads); CHKERR(rval, "Trouble creating mesh.");
+
+ // get all vertices and non-vertex entities
+ Range verts, quads;
+ rval = mbImpl->get_entities_by_handle(0, quads); CHKERR(rval, "Getting all entities.");
+ verts = quads.subset_by_dimension(0);
+ quads -= verts;
+
+ // create tag1 (element-based avg), tag2 (vertex-based avg), tag3 (# connected verts)
+ Tag tag1, tag2, tag3;
+ rval = mbImpl->tag_get_handle("tag1", 3, MB_TYPE_DOUBLE, tag1, MB_TAG_DENSE | MB_TAG_CREAT); CHKERR(rval, "Creating tag1.");
+ double def_val[3] = {0.0, 0.0, 0.0}; // need a default value for tag2 because we sum into it
+ rval = mbImpl->tag_get_handle("tag2", 3, MB_TYPE_DOUBLE, tag2, MB_TAG_DENSE | MB_TAG_CREAT, def_val); CHKERR(rval, "Creating tag2.");
+ int def_val_int = 0; // need a default value for tag3 because we increment it
+ rval = mbImpl->tag_get_handle("tag3", 1, MB_TYPE_INTEGER, tag3, MB_TAG_DENSE | MB_TAG_CREAT, &def_val_int); CHKERR(rval, "Creating tag3.");
+
+ // Get pointers to connectivity, coordinate, tag, and adjacency arrays; each of these returns a count,
+ // which should be compared to the # entities you expect to verify there's only one chunk (no holes)
+ int count, vpere;
+ EntityHandle *conn_ptr;
+ rval = mbImpl->connect_iterate(quads.begin(), quads.end(), conn_ptr, vpere, count); CHKERR(rval, "connect_iterate.");
+ assert(count == (int) quads.size()); // should end up with just one contiguous chunk of quads
+
+ double *x_ptr, *y_ptr, *z_ptr;
+ rval = mbImpl->coords_iterate(verts.begin(), verts.end(), x_ptr, y_ptr, z_ptr, count); CHKERR(rval, "coords_iterate.");
+ assert(count == (int) verts.size()); // should end up with just one contiguous chunk of vertices
+
+ double *tag1_ptr, *tag2_ptr;
+ int *tag3_ptr;
+ rval = mbImpl->tag_iterate(tag1, quads.begin(), quads.end(), count, (void*&)tag1_ptr); CHKERR(rval, "tag1_iterate.");
+ assert(count == (int) quads.size()); // should end up with just one contiguous chunk of quads
+ rval = mbImpl->tag_iterate(tag2, quads.begin(), quads.end(), count, (void*&)tag2_ptr); CHKERR(rval, "tag2_iterate.");
+ assert(count == (int) quads.size()); // should end up with just one contiguous chunk of quads
+ rval = mbImpl->tag_iterate(tag3, quads.begin(), quads.end(), count, (void*&)tag3_ptr); CHKERR(rval, "tag3_iterate.");
+ assert(count == (int) quads.size()); // should end up with just one contiguous chunk of quads
+
+ const std::vector<EntityHandle> **adjs_ptr;
+ rval = mbImpl->adjacencies_iterate(verts.begin(), verts.end(), adjs_ptr, count); CHKERR(rval, "adjacencies_iterate.");
+ assert(count == (int) verts.size()); // should end up with just one contiguous chunk of vertices
+
+ // Start_ handles used to compute indices into vertex/quad arrays; can use direct subtraction because we know
+ // there aren't any holes in the handle spaces for verts or quads
+ EntityHandle start_vert = *verts.begin(), start_quad = *quads.begin();
+
+ // iterate over elements, computing tag1 from coords positions
+ for (int i = 0; i < nquads; i++) {
+ tag1_ptr[3*i+0] = tag1_ptr[3*i+1] = tag1_ptr[3*i+2] = 0.0; // initialize position for this element
+ for (int j = 0; j < vpere; j++) { // loop over vertices in this element
+ int v_index = conn_ptr[vpere*i+j] - start_vert; // vert index is just the offset from start vertex
+ tag1_ptr[3*i+0] += x_ptr[v_index];
+ tag1_ptr[3*i+1] += y_ptr[v_index]; // sum vertex positions into tag1...
+ tag1_ptr[3*i+2] += z_ptr[v_index];
+ }
+ tag1_ptr[3*i+0] /= vpere;
+ tag1_ptr[3*i+1] /= vpere; // then normalize
+ tag1_ptr[3*i+2] /= vpere;
+ } // loop over elements in chunk
+
+ // Iterate through vertices, summing positions into tag2 on connected elements and incrementing vertex count
+ for (int v = 0; v < count; v++) {
+ const std::vector<EntityHandle> *avec = *(adjs_ptr+v);
+ for (std::vector<EntityHandle>::const_iterator ait = avec->begin(); ait != avec->end(); ait++) {
+ // *ait is the quad handle, its index is computed by subtracting the start quad handle
+ int a_ind = *ait - start_quad;
+ tag2_ptr[3*a_ind+0] += x_ptr[v]; // tag on each element is 3 doubles, x/y/z
+ tag2_ptr[3*a_ind+1] += y_ptr[v];
+ tag2_ptr[3*a_ind+2] += z_ptr[v];
+ tag3_ptr[a_ind]++; // increment the vertex count
+ }
+ }
+
+ // Normalize tag2 by vertex count (tag3); loop over elements using same approach as before
+ // At the same time, compare values of tag1 and tag2
+ int n_dis = 0;
+ for (Range::iterator q_it = quads.begin(); q_it != quads.end(); q_it++) {
+ int i = *q_it - start_quad;
+ for (int j = 0; j < 3; j++) tag2_ptr[3*i+j] /= (double)tag3_ptr[i]; // normalize by # verts
+ if (tag1_ptr[3*i] != tag2_ptr[3*i] || tag1_ptr[3*i+1] != tag2_ptr[3*i+1] || tag1_ptr[3*i+2] != tag2_ptr[3*i+2]) {
+ std::cout << "Tag1, tag2 disagree for element " << *q_it + i << std::endl;
+ n_dis++;
+ }
+ }
+ if (!n_dis) std::cout << "All tags agree, success!" << std::endl;
+
+ // Ok, we're done, shut down MOAB
+ delete mbImpl;
+
+ return 0;
+}
+
+ErrorCode create_mesh_no_holes(Interface *mbImpl, int nquads)
+{
+ // first make the mesh, a 1d array of quads with left hand side x = elem_num; vertices are numbered in layers
+ ReadUtilIface *read_iface;
+ ErrorCode rval = mbImpl->query_interface(read_iface); CHKERR(rval, "query_interface");
+ std::vector<double *> coords;
+ EntityHandle start_vert, start_elem, *connect;
+ // create verts, num is 2(nquads+1) because they're in a 1d row; will initialize coords in loop over quads later
+ rval = read_iface->get_node_coords (3, 2*(nquads+1), 0, start_vert, coords); CHKERR(rval, "get_node_arrays");
+ // create quads
+ rval = read_iface->get_element_connect(nquads, 4, MBQUAD, 0, start_elem, connect); CHKERR(rval, "get_element_connect");
+ for (int i = 0; i < nquads; i++) {
+ coords[0][2*i] = coords[0][2*i+1] = (double) i; // x values are all i
+ coords[1][2*i] = 0.0; coords[1][2*i+1] = 1.0; // y coords
+ coords[2][2*i] = coords[2][2*i+1] = (double) 0.0; // z values, all zero (2d mesh)
+ EntityHandle quad_v = start_vert + 2*i;
+ connect[4*i+0] = quad_v;
+ connect[4*i+1] = quad_v+2;
+ connect[4*i+2] = quad_v+3;
+ connect[4*i+3] = quad_v+1;
+ }
+
+ // last two vertices
+ coords[0][2*nquads] = coords[0][2*nquads+1] = (double) nquads;
+ coords[1][2*nquads] = 0.0; coords[1][2*nquads+1] = 1.0; // y coords
+ coords[2][2*nquads] = coords[2][2*nquads+1] = (double) 0.0; // z values, all zero (2d mesh)
+
+ // call a vertex-quad adjacencies function to generate vertex-element adjacencies in MOAB
+ Range dum_range;
+ rval = mbImpl->get_adjacencies(&start_vert, 1, 2, false, dum_range); CHKERR(rval, "get_adjacencies");
+ assert(!dum_range.empty());
+
+ return MB_SUCCESS;
+}
+
+
+
diff --git a/examples/DirectAccessNoHolesF90.F90 b/examples/DirectAccessNoHolesF90.F90
new file mode 100644
index 0000000..3cf5fa7
--- /dev/null
+++ b/examples/DirectAccessNoHolesF90.F90
@@ -0,0 +1,231 @@
+! @example DirectAccessNoHolesF90.F90 \n
+! \brief Use direct access to MOAB data to avoid calling through API, in Fortran90 \n
+!
+! This example creates a 1d row of quad elements, such that all quad and vertex handles
+! are contiguous in the handle space and in the database. Then it shows how to get access
+! to pointers to MOAB-native data for vertex coordinates, quad connectivity, and tag storage
+! (vertex to quad adjacency lists aren't accessible from Fortran because they are std::vector's).
+! This allows applications to access this data directly
+! without going through MOAB's API. In cases where the mesh is not changing (or only mesh
+! vertices are moving), this can save significant execution time in applications.
+!
+! Using direct access (or MOAB in general) from Fortran is complicated in two specific ways:
+! 1) There is no way to use arrays with specified dimension starting/ending values with ISO_C_BINDING;
+! therefore, all arrays passed back from MOAB/iMesh must use 1-based indices; this makes it a bit
+! more difficult to traverse the direct arrays. In this example, I explicitly use indices that
+! look like my_array(1+v_ind...) to remind users of this.
+! 2) Arithmetic on handles is complicated by the fact that Fortran has no unsigned integer type. I get
+! around this by assigning to a C-typed variable first, before the handle arithmetic. This seems to
+! work fine. C-typed variables are part of the Fortran95 standard.
+!
+! ----------------------
+! | | | |
+! | | | | ...
+! | | | |
+! ----------------------
+!
+! -# Initialize MOAB \n
+! -# Create a quad mesh, as depicted above
+! -# Create 2 dense tags (tag1, tag2) for avg position to assign to quads, and # verts per quad (tag3)
+! -# Get connectivity, coordinate, tag1 iterators
+! -# Iterate through quads, computing midpoint based on vertex positions, set on quad-based tag1
+! -# Iterate through vertices, summing positions into tag2 on connected quads and incrementing vertex count
+! -# Iterate through quads, normalizing tag2 by vertex count and comparing values of tag1 and tag2
+!
+! <b>To compile</b>: \n
+! make DirectAccessNoHolesF90 MOAB_DIR=<installdir> \n
+! <b>To run</b>: ./DirectAccessNoHolesF90 [-nquads <# quads>]\n
+!
+!
+
+#define CHECK(a) \
+ if (a .ne. iBase_SUCCESS) call exit(a)
+
+module freem
+ interface
+ subroutine c_free(ptr) bind(C, name='free')
+ use ISO_C_BINDING
+ type(C_ptr), value, intent(in) :: ptr
+ end subroutine c_free
+ end interface
+end module freem
+
+program DirectAccessNoHolesF90
+ use ISO_C_BINDING
+ use freem
+ implicit none
+#include "iMesh_f.h"
+
+ ! declarations
+ iMesh_Instance imesh
+ iBase_EntitySetHandle root_set
+ integer ierr, nquads, nquads_tmp, nverts
+ iBase_TagHandle tag1h, tag2h, tag3h
+ TYPE(C_PTR) verts_ptr, quads_ptr, conn_ptr, x_ptr, y_ptr, z_ptr, tag1_ptr, tag2_ptr, tag3_ptr, &
+ tmp_quads_ptr, offsets_ptr ! pointers that are equivalence'd to arrays using c_f_pointer
+ real(C_DOUBLE), pointer :: x(:), y(:), z(:), tag1(:), tag2(:) ! arrays equivalence'd to pointers
+ integer, pointer :: tag3(:), offsets(:) ! arrays equivalence'd to pointers
+ iBase_EntityHandle, pointer :: quads(:), verts(:), conn(:), tmp_quads(:) ! arrays equivalence'd to pointers
+ iBase_EntityArrIterator viter, qiter
+ integer(C_SIZE_T) :: startv, startq, v_ind, e_ind ! needed to do handle arithmetic
+ integer count, vpere, e, i, j, v, offsets_size, tmp_quads_size, n_dis
+ character*120 opt
+
+ ! initialize interface and get root set
+ call iMesh_newMesh("", imesh, ierr); CHECK(ierr)
+ call iMesh_getRootSet(%VAL(imesh), root_set, ierr); CHECK(ierr)
+
+ ! create mesh
+ nquads_tmp = 1000
+ call create_mesh_no_holes(imesh, nquads_tmp, ierr); CHECK(ierr)
+
+ ! get all vertices and quads as arrays
+ nverts = 0
+ call iMesh_getEntities(%VAL(imesh), %VAL(root_set), %VAL(0), %VAL(iBase_VERTEX), &
+ verts_ptr, nverts, nverts, ierr); CHECK(ierr)
+ call c_f_pointer(verts_ptr, verts, [nverts])
+ nquads = 0
+ call iMesh_getEntities(%VAL(imesh), %VAL(root_set), %VAL(2), %VAL(iMesh_QUADRILATERAL), &
+ quads_ptr, nquads, nquads, ierr); CHECK(ierr)
+ call c_f_pointer(quads_ptr, quads, [nquads])
+
+ ! First vertex/quad is at start of vertex/quad list, and is offset for vertex/quad index calculation
+ startv = verts(1); startq = quads(1)
+ call c_free(quads_ptr) ! free memory we allocated in MOAB
+
+ ! create tag1 (element-based avg), tag2 (vertex-based avg)
+ opt = 'moab:TAG_STORAGE_TYPE=DENSE moab:TAG_DEFAULT_VALUE=0'
+ call iMesh_createTagWithOptions(%VAL(imesh), 'tag1', opt, %VAL(3), %VAL(iBase_DOUBLE), &
+ tag1h, ierr, %VAL(4), %VAL(56)); CHECK(ierr)
+ call iMesh_createTagWithOptions(%VAL(imesh), 'tag2', opt, %VAL(3), %VAL(iBase_DOUBLE), &
+ tag2h, ierr, %VAL(4), %VAL(56)); CHECK(ierr)
+
+ ! Get iterator to verts, then pointer to coordinates
+ call iMesh_initEntArrIter(%VAL(imesh), %VAL(root_set), %VAL(iBase_VERTEX), %VAL(iMesh_ALL_TOPOLOGIES), &
+ %VAL(nverts), %VAL(0), viter, ierr); CHECK(ierr)
+ call iMesh_coordsIterate(%VAL(imesh), %VAL(viter), x_ptr, y_ptr, z_ptr, count, ierr); CHECK(ierr)
+ CHECK(count-nverts) ! check that all verts are in one contiguous chunk
+ call c_f_pointer(x_ptr, x, [nverts]); call c_f_pointer(y_ptr, y, [nverts]); call c_f_pointer(z_ptr, z, [nverts])
+
+ ! Get iterator to quads, then pointers to connectivity and tags
+ call iMesh_initEntArrIter(%VAL(imesh), %VAL(root_set), %VAL(iBase_FACE), %VAL(iMesh_ALL_TOPOLOGIES), &
+ %VAL(nquads), %VAL(0), qiter, ierr); CHECK(ierr)
+ call iMesh_connectIterate(%VAL(imesh), %VAL(qiter), conn_ptr, vpere, count, ierr); CHECK(ierr)
+ call c_f_pointer(conn_ptr, conn, [vpere*nquads])
+ call iMesh_tagIterate(%VAL(imesh), %VAL(tag1h), %VAL(qiter), tag1_ptr, count, ierr); CHECK(ierr)
+ call c_f_pointer(tag1_ptr, tag1, [3*nquads])
+ call iMesh_tagIterate(%VAL(imesh), %VAL(tag2h), %VAL(qiter), tag2_ptr, count, ierr); CHECK(ierr)
+ call c_f_pointer(tag2_ptr, tag2, [3*nquads])
+
+ ! iterate over elements, computing tag1 from coords positions; use (1+... in array indices for 1-based indexing
+ do i = 0, nquads-1
+ tag1(1+3*i+0) = 0.0; tag1(1+3*i+1) = 0.0; tag1(1+3*i+2) = 0.0 ! initialize position for this element
+ do j = 0, vpere-1 ! loop over vertices in this quad
+ v_ind = conn(1+vpere*i+j) ! assign to v_ind to allow handle arithmetic
+ v_ind = v_ind - startv ! vert index is just the offset from start vertex
+ tag1(1+3*i+0) = tag1(1+3*i+0) + x(1+v_ind)
+ tag1(1+3*i+1) = tag1(1+3*i+1) + y(1+v_ind) ! sum vertex positions into tag1...
+ tag1(1+3*i+2) = tag1(1+3*i+2) + z(1+v_ind)
+ end do
+ tag1(1+3*i+0) = tag1(1+3*i+0) / vpere;
+ tag1(1+3*i+1) = tag1(1+3*i+1) / vpere; ! then normalize
+ tag1(1+3*i+2) = tag1(1+3*i+2) / vpere;
+ end do ! loop over elements in chunk
+
+ ! Iterate through vertices, summing positions into tag2 on connected elements; get adjacencies all at once,
+ ! could also get them per vertex (time/memory tradeoff)
+ tmp_quads_size = 0
+ offsets_size = 0
+ call iMesh_getEntArrAdj(%VAL(imesh), verts, %VAL(nverts), %VAL(iMesh_QUADRILATERAL), &
+ tmp_quads_ptr, tmp_quads_size, tmp_quads_size, &
+ offsets_ptr, offsets_size, offsets_size, ierr); CHECK(ierr)
+ call c_f_pointer(tmp_quads_ptr, tmp_quads, [tmp_quads_size])
+ call c_f_pointer(offsets_ptr, offsets, [offsets_size])
+ call c_free(verts_ptr) ! free memory we allocated in MOAB
+ do v = 0, nverts-1
+ do e = 1+offsets(1+v), 1+offsets(1+v+1)-1 ! 1+ because offsets are in terms of 0-based arrays
+ e_ind = tmp_quads(e) ! assign to e_ind to allow handle arithmetic
+ e_ind = e_ind - startq ! e_ind is the quad handle minus the starting quad handle
+ tag2(1+3*e_ind+0) = tag2(1+3*e_ind+0) + x(1+v)
+ tag2(1+3*e_ind+1) = tag2(1+3*e_ind+1) + y(1+v) ! tag on each element is 3 doubles, x/y/z
+ tag2(1+3*e_ind+2) = tag2(1+3*e_ind+2) + z(1+v)
+ end do
+ end do
+ call c_free(tmp_quads_ptr) ! free memory we allocated in MOAB
+ call c_free(offsets_ptr) ! free memory we allocated in MOAB
+
+ ! Normalize tag2 by vertex count (vpere); loop over elements using same approach as before
+ ! At the same time, compare values of tag1 and tag2
+ n_dis = 0
+ do e = 0, nquads-1
+ do j = 0, 2
+ tag2(1+3*e+j) = tag2(1+3*e+j) / vpere; ! normalize by # verts
+ end do
+ if (tag1(1+3*e) .ne. tag2(1+3*e) .or. tag1(1+3*e+1) .ne. tag2(1+3*e+1) .or. tag1(1+3*e+2) .ne. tag2(1+3*i+2)) then
+ print *, "Tag1, tag2 disagree for element ", e
+ print *, "tag1: ", tag1(1+3*e), tag1(1+3*e+1), tag1(1+3*e+2)
+ print *, "tag2: ", tag2(1+3*e), tag2(1+3*e+1), tag2(1+3*e+2)
+ n_dis = n_dis + 1
+ end if
+ end do
+ if (n_dis .eq. 0) print *, 'All tags agree, success!'
+
+ ! Ok, we are done, shut down MOAB
+ call iMesh_dtor(%VAL(imesh), ierr)
+ return
+end program DirectAccessNoHolesF90
+
+subroutine create_mesh_no_holes(imesh, nquads, ierr)
+ use ISO_C_BINDING
+ use freem
+ implicit none
+#include "iMesh_f.h"
+
+ iMesh_Instance imesh
+ integer nquads, ierr
+
+ real(C_DOUBLE), pointer :: coords(:,:)
+ TYPE(C_PTR) connect_ptr, tmp_ents_ptr, stat_ptr
+ iBase_EntityHandle, pointer :: connect(:), tmp_ents(:)
+ integer, pointer :: stat(:)
+ integer nverts, tmp_size, stat_size, i
+
+ ! first make the mesh, a 1d array of quads with left hand side x = elem_num; vertices are numbered in layers
+ ! create verts, num is 2(nquads+1) because they're in a 1d row
+ nverts = 2*(nquads+1)
+ allocate(coords(0:2, 0:nverts-1))
+ do i = 0, nquads-1
+ coords(0,2*i) = i; coords(0,2*i+1) = i ! x values are all i
+ coords(1,2*i) = 0.0; coords(1,2*i+1) = 1.0 ! y coords
+ coords(2,2*i) = 0.0; coords(2,2*i+1) = 0.0 ! z values, all zero (2d mesh)
+ end do
+ ! last two vertices
+ coords(0, nverts-2) = nquads; coords(0, nverts-1) = nquads
+ coords(1, nverts-2) = 0.0; coords(1, nverts-1) = 1.0 ! y coords
+ coords(2, nverts-2) = 0.0; coords(2, nverts-1) = 0.0 ! z values, all zero (2d mesh)
+ tmp_size = 0
+ call iMesh_createVtxArr(%VAL(imesh), %VAL(nverts), %VAL(iBase_INTERLEAVED), &
+ coords, %VAL(3*nverts), tmp_ents_ptr, tmp_size, tmp_size, ierr); CHECK(ierr)
+ call c_f_pointer(tmp_ents_ptr, tmp_ents, [nverts])
+ deallocate(coords)
+ allocate(connect(0:4*nquads-1))
+ do i = 0, nquads-1
+ connect(4*i+0) = tmp_ents(1+2*i)
+ connect(4*i+1) = tmp_ents(1+2*i+2)
+ connect(4*i+2) = tmp_ents(1+2*i+3)
+ connect(4*i+3) = tmp_ents(1+2*i+1)
+ end do
+ stat_size = 0
+ stat_ptr = C_NULL_PTR
+ ! re-use tmp_ents here; we know it'll always be larger than nquads so it's ok
+ call iMesh_createEntArr(%VAL(imesh), %VAL(iMesh_QUADRILATERAL), connect, %VAL(4*nquads), &
+ tmp_ents_ptr, tmp_size, tmp_size, stat_ptr, stat_size, stat_size, ierr); CHECK(ierr)
+
+ ierr = iBase_SUCCESS
+
+ call c_free(tmp_ents_ptr)
+ call c_free(stat_ptr)
+ deallocate(connect)
+
+ return
+end subroutine create_mesh_no_holes
diff --git a/examples/DirectAccessWithHoles.cpp b/examples/DirectAccessWithHoles.cpp
new file mode 100644
index 0000000..524c30f
--- /dev/null
+++ b/examples/DirectAccessWithHoles.cpp
@@ -0,0 +1,239 @@
+/** @example DirectAccess.cpp \n
+ * \brief Use direct access to MOAB data to avoid calling through API \n
+ *
+ * This example creates a 1d row of quad elements, with a user-specified number of "holes" (missing quads) in the row:
+ *
+ * ---------------------- ---------------------- --------
+ * | | | | | | | | | |
+ * | | | |(hole)| | | |(hole)| | ...
+ * | | | | | | | | | |
+ * ---------------------- ---------------------- --------
+ *
+ * This makes (nholes+1) contiguous runs of quad handles in the handle space
+ * This example shows how to use the xxx_iterate functions in MOAB (xxx = coords, connect, tag, adjacencies) to get
+ * direct pointer access to MOAB internal storage, which can be used without calling through the MOAB API.
+ *
+ * -# Initialize MOAB \n
+ * -# Create a quad mesh with holes, as depicted above
+ * -# Create 2 dense tags (tag1, tag2) for avg position to assign to quads, and # verts per quad (tag3)
+ * -# Get connectivity, coordinate, tag1 iterators
+ * -# Iterate through quads, computing midpoint based on vertex positions, set on quad-based tag1
+ * -# Set up map from starting quad handle for a chunk to struct of (tag1_ptr, tag2_ptr, tag3_ptr), pointers to
+ * the dense tag storage for those tags for the chunk
+ * -# Iterate through vertices, summing positions into tag2 on connected quads and incrementing vertex count
+ * -# Iterate through quads, normalizing tag2 by vertex count and comparing values of tag1 and tag2
+ *
+ * <b>To compile</b>: \n
+ * make DirectAccessWithHoles MOAB_DIR=<installdir> \n
+ * <b>To run</b>: ./DirectAccess [-nquads <# quads>] [-holes <# holes>]\n
+ *
+ */
+
+#include "moab/Core.hpp"
+#include "moab/ProgOptions.hpp"
+#include "moab/ReadUtilIface.hpp"
+#include <map>
+#include <iostream>
+#include <assert.h>
+
+// Error routines for use with MOAB API
+#define CHKERR(CODE, MSG) do {if (MB_SUCCESS != (CODE)) {std::string errstr; mbImpl->get_last_error(errstr); \
+ std::cerr << errstr << std::endl; std::cerr << MSG << std::endl; return CODE;}} while(false)
+
+using namespace moab;
+
+ErrorCode create_mesh_with_holes(Interface *mbImpl, int nquads, int nholes);
+
+struct tag_struct {double *avg_ptr; int *nv_ptr;};
+
+int main(int argc, char **argv)
+{
+ // get MOAB instance
+ Interface *mbImpl = new Core;
+ if (NULL == mbImpl) return 1;
+
+ int nquads = 1000, nholes = 1;
+
+ // parse options
+ ProgOptions opts;
+ opts.addOpt<int>(std::string("nquads,n"), std::string("Number of quads in the mesh (default = 1000"), &nquads);
+ opts.addOpt<int>(std::string("holes,H"), std::string("Number of holes in the element handle space (default = 1"), &nholes);
+ opts.parseCommandLine(argc, argv);
+ if (nholes >= nquads) {
+ std::cerr << "Number of holes needs to be < number of elements." << std::endl;
+ return 1;
+ }
+
+ // create simple structured mesh with hole, but using unstructured representation
+ ErrorCode rval = create_mesh_with_holes(mbImpl, nquads, nholes); CHKERR(rval, "Trouble creating mesh.");
+
+ // get all vertices and non-vertex entities
+ Range verts, elems;
+ rval = mbImpl->get_entities_by_handle(0, elems); CHKERR(rval, "Getting all entities.");
+ verts = elems.subset_by_dimension(0);
+ elems -= verts;
+
+ // create tag1 (element-based avg), tag2 (vertex-based avg), tag3 (# connected verts)
+ Tag tag1, tag2, tag3;
+ rval = mbImpl->tag_get_handle("tag1", 3, MB_TYPE_DOUBLE, tag1, MB_TAG_DENSE | MB_TAG_CREAT); CHKERR(rval, "Creating tag1.");
+ double def_val[3] = {0.0, 0.0, 0.0}; // need a default value for tag2 because we sum into it
+ rval = mbImpl->tag_get_handle("tag2", 3, MB_TYPE_DOUBLE, tag2, MB_TAG_DENSE | MB_TAG_CREAT, def_val); CHKERR(rval, "Creating tag2.");
+ int def_val_int = 0; // need a default value for tag3 because we increment it
+ rval = mbImpl->tag_get_handle("tag3", 1, MB_TYPE_INTEGER, tag3, MB_TAG_DENSE | MB_TAG_CREAT, &def_val_int); CHKERR(rval, "Creating tag3.");
+
+ // Get connectivity, coordinate, tag, and adjacency iterators
+ EntityHandle *conn_ptr;
+ double *x_ptr, *y_ptr, *z_ptr, *tag1_ptr, *tag2_ptr;
+ int *tag3_ptr;
+
+ // First vertex is at start of range (ranges are sorted), and is offset for vertex index calculation
+ EntityHandle first_vert = *verts.begin();
+
+ // When iterating over elements, each chunk can have a different # vertices; also, count tells you how many
+ // elements are in the current chunk
+ int vpere, count;
+
+ // Get coordinates iterator, just need this once because we know verts handle space doesn't have holes
+ rval = mbImpl->coords_iterate(verts.begin(), verts.end(), x_ptr, y_ptr, z_ptr, count); CHKERR(rval, "coords_iterate.");
+ assert(count == (int) verts.size()); // should end up with just one contiguous chunk of vertices
+
+ // Iterate through elements, computing midpoint based on vertex positions, set on element-based tag1
+ // Control loop by iterator over elem range
+ Range::iterator e_it = elems.begin();
+
+ while (e_it != elems.end()) {
+ // get conn_ptr, tag1_ptr for next contiguous chunk of element handles, and coords pointers for all verts
+ rval = mbImpl->connect_iterate(e_it, elems.end(), conn_ptr, vpere, count); CHKERR(rval, "connect_iterate.");
+ rval = mbImpl->tag_iterate(tag1, e_it, elems.end(), count, (void*&)tag1_ptr); CHKERR(rval, "tag1_iterate.");
+
+ // iterate over elements in this chunk
+ for (int i = 0; i < count; i++) {
+ tag1_ptr[0] = tag1_ptr[1] = tag1_ptr[2] = 0.0; // initialize position for this element
+ for (int j = 0; j < vpere; j++) { // loop over vertices in this element
+ int v_index = conn_ptr[j] - first_vert; // vert index is just the offset from first vertex
+ tag1_ptr[0] += x_ptr[v_index];
+ tag1_ptr[1] += y_ptr[v_index]; // sum vertex positions into tag1...
+ tag1_ptr[2] += z_ptr[v_index];
+ }
+ tag1_ptr[0] /= vpere;
+ tag1_ptr[1] /= vpere; // then normalize
+ tag1_ptr[2] /= vpere;
+
+ // done with this element; advance connect_ptr and tag1_ptr to next element
+ conn_ptr += vpere;
+ tag1_ptr += 3;
+ } // loop over elements in chunk
+
+ // done with chunk; advance range iterator by count; will skip over gaps in range
+ e_it += count;
+ } // while loop over all elements
+
+ // Iterate through vertices, summing positions into tag2 on connected elements and incrementing vertex count
+ // Iterate over chunks the same as elements, even though we know we have only one chunk here, just to show
+ // how it's done
+
+ // Create a std::map from EntityHandle (first entity handle in chunk) to
+ // tag_struct (ptrs to start of avg/#verts tags for that chunk); then for a given entity handle, we can quickly
+ // find the chunk it's in using map::lower_bound; could have set up this map in earlier loop over elements, but do
+ // it here for clarity
+
+ std::map< EntityHandle, tag_struct> elem_map;
+ e_it = elems.begin();
+ while (e_it != elems.end()) {
+ tag_struct ts = {NULL, NULL};
+ rval = mbImpl->tag_iterate(tag2, e_it, elems.end(), count, (void*&)ts.avg_ptr); CHKERR(rval, "tag2_iterate.");
+ rval = mbImpl->tag_iterate(tag3, e_it, elems.end(), count, (void*&)ts.nv_ptr); CHKERR(rval, "tag3_iterate.");
+ elem_map[*e_it] = ts;
+ e_it += count;
+ }
+
+ // call a vertex-quad adjacencies function to generate vertex-element adjacencies in MOAB
+ Range::iterator v_it = verts.begin();
+ Range dum_range;
+ rval = mbImpl->get_adjacencies(&(*v_it), 1, 2, false, dum_range); CHKERR(rval, "get_adjacencies");
+ const std::vector<EntityHandle> **adjs_ptr;
+ while (v_it != verts.end()) {
+ // get coords ptrs, adjs_ptr; can't set tag2_ptr by direct access, because of hole in element handle space
+ rval = mbImpl->coords_iterate(v_it, verts.end(), x_ptr, y_ptr, z_ptr, count); CHKERR(rval, "coords_iterate.");
+ rval = mbImpl->adjacencies_iterate(v_it, verts.end(), adjs_ptr, count); CHKERR(rval, "adjacencies_iterate.");
+
+ for (int v = 0; v < count; v++) {
+ const std::vector<EntityHandle> *avec = *(adjs_ptr+v);
+ for (std::vector<EntityHandle>::const_iterator ait = avec->begin(); ait != avec->end(); ait++) {
+ // get chunk that this element resides in; upper_bound points to the first element strictly > key, so get that and decrement
+ // (would work the same as lower_bound with an if-test and decrement)
+ std::map<EntityHandle, tag_struct>::iterator mit = elem_map.upper_bound(*ait); mit--;
+ // index of *ait in that chunk
+ int a_ind = *ait - (*mit).first;
+ tag_struct ts = (*mit).second;
+ ts.avg_ptr[3*a_ind] += x_ptr[v]; // tag on each element is 3 doubles, x/y/z
+ ts.avg_ptr[3*a_ind+1] += y_ptr[v];
+ ts.avg_ptr[3*a_ind+2] += z_ptr[v];
+ ts.nv_ptr[a_ind]++; // increment the vertex count
+ }
+ }
+
+ v_it += count;
+ }
+
+ // Normalize tag2 by vertex count; loop over elements using same approach as before
+ // At the same time, compare values of tag1 and tag2
+ e_it = elems.begin();
+ while (e_it != elems.end()) {
+ // get conn_ptr, tag1_ptr for next contiguous chunk of element handles, and coords pointers for all verts
+ rval = mbImpl->tag_iterate(tag1, e_it, elems.end(), count, (void*&)tag1_ptr); CHKERR(rval, "tag1_iterate.");
+ rval = mbImpl->tag_iterate(tag2, e_it, elems.end(), count, (void*&)tag2_ptr); CHKERR(rval, "tag2_iterate.");
+ rval = mbImpl->tag_iterate(tag3, e_it, elems.end(), count, (void*&)tag3_ptr); CHKERR(rval, "tag3_iterate.");
+
+ // iterate over elements in this chunk
+ for (int i = 0; i < count; i++) {
+ for (int j = 0; j < 3; j++) tag2_ptr[3*i+j] /= (double)tag3_ptr[i]; // normalize by # verts
+ if (tag1_ptr[3*i] != tag2_ptr[3*i] || tag1_ptr[3*i+1] != tag2_ptr[3*i+1] || tag1_ptr[3*i+2] != tag2_ptr[3*i+2])
+ std::cout << "Tag1, tag2 disagree for element " << *e_it + i << std::endl;
+ }
+
+ e_it += count;
+ }
+
+ // Ok, we're done, shut down MOAB
+ delete mbImpl;
+
+ return 0;
+}
+
+ErrorCode create_mesh_with_holes(Interface *mbImpl, int nquads, int nholes)
+{
+ // first make the mesh, a 1d array of quads with left hand side x = elem_num; vertices are numbered in layers
+ ReadUtilIface *read_iface;
+ ErrorCode rval = mbImpl->query_interface(read_iface); CHKERR(rval, "query_interface");
+ std::vector<double *> coords;
+ EntityHandle start_vert, start_elem, *connect;
+ // create verts, num is 4(nquads+1) because they're in a 1d row; will initialize coords in loop over elems later
+ rval = read_iface->get_node_coords (3, 2*(nquads+1), 0, start_vert, coords); CHKERR(rval, "get_node_arrays");
+ // create elems
+ rval = read_iface->get_element_connect(nquads, 4, MBQUAD, 0, start_elem, connect); CHKERR(rval, "get_element_connect");
+ for (int i = 0; i < nquads; i++) {
+ coords[0][2*i] = coords[0][2*i+1] = (double) i; // x values are all i
+ coords[1][2*i] = 0.0; coords[1][2*i+1] = 1.0; // y coords
+ coords[2][2*i] = coords[2][2*i+1] = (double) 0.0; // z values, all zero (2d mesh)
+ EntityHandle quad_v = start_vert + 2*i;
+ for (int j = 0; j < 4; j++) connect[4*i+j] = quad_v+j; // connectivity of each quad is a sequence starting from quad_v
+ }
+ // last two vertices
+ coords[0][2*nquads] = coords[0][2*nquads+1] = (double) nquads;
+ coords[1][2*nquads] = 0.0; coords[1][2*nquads+1] = 1.0; // y coords
+ coords[2][2*nquads] = coords[2][2*nquads+1] = (double) 0.0; // z values, all zero (2d mesh)
+
+ // now delete nholes elements, spaced approximately equally through mesh, so contiguous size is about nquads/(nholes+1)
+ // reinterpret start_elem as the next element to be deleted
+ int de = nquads / (nholes + 1);
+ for (int i = 0; i < nholes; i++) {
+ start_elem += de;
+ rval = mbImpl->delete_entities(&start_elem, 1); CHKERR(rval, "delete_entities");
+ }
+
+ return MB_SUCCESS;
+}
+
+
+
diff --git a/examples/FileRead.cpp b/examples/FileRead.cpp
deleted file mode 100644
index bff9e28..0000000
--- a/examples/FileRead.cpp
+++ /dev/null
@@ -1,164 +0,0 @@
-#include <iostream>
-#include <fstream>
-#include <vector>
-#include <string>
-#include <sstream>
-
-#include "moab/Core.hpp"
-#include "moab/ReadUtilIface.hpp"
-
-using namespace std;
-using namespace moab;
-
-int comment(string & line)
-{
- // if a line starts with '#' is a comment
- // eat white space characters
- size_t found=line.find_first_not_of(" \t");
- if (found==string::npos)
- return 1; // empty line
- if ('#'==line[found])
- return 1; // a comment indeed
- return 0; // a line with some data in it, then
-
-}
-ErrorCode ReadTriangleOutput( Interface *mb, string fileBase ) {
-
- //
- // get the read interface from moab
- ReadUtilIface *iface;
- ErrorCode rval = mb->query_interface(iface);
- //
- if (MB_SUCCESS != rval)
- {
- cout<<"Can't get interface.\n";
- return MB_FAILURE;
- }
- // Triangle default <name>.node
- string nodeFileName = fileBase+".node";
- ifstream nodeFile (nodeFileName.c_str());
- if (!nodeFile.is_open())
- {
- cout<<"can't open node file .\n";
- return MB_FILE_DOES_NOT_EXIST;
- }
- cout << "reading nodes from file " << nodeFileName.c_str() << endl;
-
- string eleFileName = fileBase+".ele";
- ifstream eleFile (eleFileName.c_str());
- if (!eleFile.is_open())
- {
- cout<<"can't open element file .\n";
- return MB_FILE_DOES_NOT_EXIST;
- }
- cout << "reading elements from file " << eleFileName.c_str() << endl;
-
- string line;
-
- // ignore comment lines that start with #
-
- int num_nodes=0, num_triangles=0;
- while(num_nodes==0)
- {
- getline(nodeFile, line);
- if (comment(line))
- continue;
- stringstream tks(line);
- tks >> num_nodes; // ignore the rest of the first line
- // maybe will read attributes some other time
- cout << "num nodes:" << num_nodes << endl;
- }
-
- // allocate a block of vertex handles and read xyz’s into them
- // we know the size of the node arrays, and this call will allocate
- // needed arrays, coordinate arrays
- // also, it will return a starting handle for the node sequence
- vector<double*> arrays;
- EntityHandle startv;
- rval = iface->get_node_coords(2, num_nodes, 0, startv, arrays);
- for (int i = 0; i < num_nodes; i++)
- {
- getline(nodeFile, line);
- if (comment(line))
- {
- i--;// read one more line
- continue;
- }
- stringstream tokens(line);
- int nodeId;
- tokens >> nodeId >> arrays[0][i] >> arrays[1][i] ;
- }
-
- // now read the element data from a different file
- // first, find out how many elements are out there
- // first line with data should have it
- while(num_triangles==0)
- {
- getline(eleFile, line);
- if (comment(line))
- continue;
- stringstream tks(line);
- tks >> num_triangles; // ignore the rest of the line
- cout << "num triangles:" << num_triangles << endl;
- }
-
- EntityHandle starte;
- EntityHandle *starth; // the connectivity array that will get populated
- // with triangle data
- // allocate block of triangle handles and read connectivity into them
- rval = iface->get_element_connect(num_triangles, 3, MBTRI, 0, starte, starth);
-
- for (int j = 0; j < num_triangles; j++)
- {
- getline(eleFile, line);
- if (comment(line))
- {
- j--;// read one more line
- continue;
- }
- stringstream tokens(line);
- int eleId;
- unsigned int node;
- tokens >> eleId;
- for (int k=0; k<3; k++)
- {
- tokens >> node;
- // vertex handles start at startv
- starth[3*j+k] = startv + node - 1;
- }
- }
-
- mb->release_interface(iface);
- //
- return MB_SUCCESS;
-}
-
-
-// .
-// Read Triangle output files
-// Assume that the format is <filename>.node and <filename>.ele
-// see http://www.cs.cmu.edu/~quake/triangle.html for details
-//
-int main(int argc, char **argv) {
- if (3!=argc) {
- cout << "Usage: " << argv[0] << " <filename><outFile> " << endl;
- cout << " <filename> is the base file name; *.ele and *.node file are read; outFile is a file with an extension recognized by MOAB " << endl;
- return 0;
- }
-
- string filename = argv[1];
- char * outfile = argv[2];
-
-
- // get MOAB instance and read the file
- Core *mb = new Core();
-
- ErrorCode rval = ReadTriangleOutput(mb, filename);
-
- if (rval==MB_SUCCESS)
- {
- cout << "Writing output file " << outfile << endl;
- mb->write_file(outfile);
- }
- return 0;
-}
diff --git a/examples/GeomSetHierarchy.cpp b/examples/GeomSetHierarchy.cpp
deleted file mode 100644
index a49dafe..0000000
--- a/examples/GeomSetHierarchy.cpp
+++ /dev/null
@@ -1,83 +0,0 @@
-#include "moab/Core.hpp"
-#include "moab/Range.hpp"
-#include "MBCN.hpp"
-#include "MBTagConventions.hpp"
-#include "moab/GeomTopoTool.hpp"
-#include <iostream>
-
-const char *ent_names[] = {"Vertex", "Edge", "Face", "Region"};
-
-int main(int argc, char **argv) {
- if (1 == argc) {
- std::cout << "Usage: " << argv[0] << " <filename>" << std::endl;
- return 0;
- }
-
- // instantiate & load a file
- moab::Interface *mb = new moab::Core();
- moab::ErrorCode rval = mb->load_file(argv[1]);
-
- // get the geometric topology tag handle
- moab::Tag geom_tag, gid_tag;
- rval = mb->tag_get_handle(GEOM_DIMENSION_TAG_NAME, 1, moab::MB_TYPE_INTEGER, geom_tag);
- rval = mb->tag_get_handle(GLOBAL_ID_TAG_NAME, 1, moab::MB_TYPE_INTEGER, gid_tag);
-
- // traverse the model, from dimension 3 downward
- moab::Range psets, chsets;
- std::vector<moab::EntityHandle> sense_ents;
- std::vector<int> senses, pgids;
- int dim, pgid, chgid;
- void *dim_ptr = &dim;
- int sense;
-
- moab::GeomTopoTool gt(mb, true);
-
- for (dim = 3; dim >= 0; dim--) {
- // get parents at this dimension
- chsets.clear();
- rval = mb->get_entities_by_type_and_tag(0, MBENTITYSET,
- &geom_tag, &dim_ptr, 1,
- chsets, 1, false);
-
- // for each child, get parents and do something with them
- moab::Range::iterator ch_it, p_it;
- for (ch_it = chsets.begin(); ch_it != chsets.end(); ch_it++) {
- // get the children and put in child set list
- psets.clear();
- rval = mb ->get_parent_meshsets(*ch_it, psets);
-
- rval = mb->tag_get_data(gid_tag, &(*ch_it), 1, &chgid);
-
- // print # parents
- std::cout << ent_names[dim] << " " << chgid << " has " << psets.size()
- << " parents." << std::endl;
-
- if (2 == dim) {
- for (p_it = psets.begin(); p_it != psets.end(); p_it++) {
- rval = mb->tag_get_data(gid_tag, &(*p_it), 1, &pgid);
- rval = gt.get_sense(*ch_it, *p_it, sense);
- if (moab::MB_SUCCESS != rval) continue;
- std::cout << ent_names[dim+1] << " " << pgid << ", "
- << ent_names[dim] << " " << chgid << " sense is: ";
- if (1==sense) std::cout << "FORWARD" << std::endl;
- else std::cout << "REVERSE" << std::endl;
- }
- }
- else if (1 == dim) {
- sense_ents.clear();
- senses.clear();
- rval = gt.get_senses(*ch_it, sense_ents, senses);
- if (moab::MB_SUCCESS != rval) continue;
- for (unsigned int i = 0; i < sense_ents.size(); i++) {
- rval = mb->tag_get_data(gid_tag, &sense_ents[i], 1, &pgid);
- std::cout << ent_names[dim+1] << " " << pgid << ", "
- << ent_names[dim] << " " << chgid << " sense is: ";
- if (-1 == senses[i]) std::cout << "REVERSED" << std::endl;
- else if (0 == senses[i]) std::cout << "BOTH" << std::endl;
- else if (1 == senses[i]) std::cout << "FORWARD" << std::endl;
- else std::cout << "(invalid)" << std::endl;
- }
- }
- }
- }
-}
diff --git a/examples/HelloMOAB.cpp b/examples/HelloMOAB.cpp
new file mode 100644
index 0000000..ebab971
--- /dev/null
+++ b/examples/HelloMOAB.cpp
@@ -0,0 +1,64 @@
+/** @example HelloMOAB.cpp
+ * Description: read a mesh, get the entities.\n
+ * HelloMOAB is a simple test file which is used to read meshes from VTK file and test how many entities there are.\n
+ *
+ * To run: ./HelloMOAB <meshfile>\n
+ * (default values can run if users don't specify a mesh file)
+ */
+
+
+#include "moab/Core.hpp"
+#include <iostream>
+#include <assert.h>
+
+using namespace moab;
+using namespace std;
+
+string test_file_name = string(MESH_DIR) + string("/3k-tri-sphere.vtk");
+
+int main( int argc, char** argv )
+{
+ Interface *iface = new Core;
+
+ // need option handling here for input filename
+ if (argc > 1){
+ //user has input a mesh file
+ test_file_name = argv[1];
+ }
+ //load the mesh from vtk file
+ ErrorCode rval = iface->load_mesh( test_file_name.c_str() );
+ assert( rval == MB_SUCCESS);
+
+ //get verts entities
+ Range verts;
+ rval = iface->get_entities_by_type(0, MBVERTEX, verts);
+ assert( rval == MB_SUCCESS);
+ //get edge entities
+ Range edges;
+ rval = iface->get_entities_by_type(0, MBEDGE, edges);
+ assert(rval == MB_SUCCESS);
+
+ //get triangular entities
+ Range tri;
+ rval = iface->get_entities_by_type(0, MBTRI, tri);
+ assert( rval == MB_SUCCESS);
+
+ //get quad entities
+ Range quads;
+ rval = iface->get_entities_by_type(0, MBQUAD, quads);
+ assert(rval == MB_SUCCESS);
+
+ //get hex entities
+ Range hex;
+ rval = iface->get_entities_by_type(0, MBHEX, hex);
+ assert(rval == MB_SUCCESS);
+
+ //output the number of entities
+ cout << "Number of vertices is " << verts.size() << endl;
+ cout << "Number of edges is " << edges.size() << endl;
+ cout << "Number of triangular faces is " << tri.size() << endl;
+ cout << "Number of quad faces is " << quads.size() << endl;
+ cout << "Number of hex is " << hex.size() << endl;
+
+ return 0;
+}
diff --git a/examples/HelloMoabPar.cpp b/examples/HelloMoabPar.cpp
new file mode 100644
index 0000000..a5a71aa
--- /dev/null
+++ b/examples/HelloMoabPar.cpp
@@ -0,0 +1,154 @@
+/** @example HelloMoabPar.cpp \n
+ * \brief Read mesh into MOAB in parallel \n
+ * This example shows the simplest way of telling MOAB to read in parallel.
+ *
+ * -# Initialize MPI and get the rank and number of processors \n
+ * -# Process arguments (file name and options for parallel read) \n
+ * -# Initialize MOAB \n
+ * -# Load a partitioned file in parallel; \n
+ * -# retrieve shared entities on each processor \n
+ * -# Filter owned entities among shared ones on each processor \n
+ * -# Exchange ghost layers, and repeat the reports \n
+ *
+ * <b>To compile</b>: \n
+ * make HelloMoabPar MOAB_DIR=<installdir> \n
+ * <b>To run</b>: mpiexec -np 4 HelloMoabPar \n
+ * (depending on your configuration, LD_LIBRARY_PATH may need to contain <hdf5>/lib folder)
+ *
+ */
+
+#include "moab/ParallelComm.hpp"
+#include "MBParallelConventions.h"
+#include "moab/Core.hpp"
+#include <iostream>
+
+using namespace moab;
+
+int main(int argc, char **argv)
+{
+ MPI_Init(&argc, &argv);
+
+ int nprocs, rank;
+ MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
+ MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+
+ std::string filename;
+ std::string options;
+ if (3 != argc)
+ {
+ if (rank == 0)
+ {
+ std::cout << "Usage: " << argv[0] << " <filename><options (separated by;)>\n ";
+ }
+ /* this file has a partition with 4 parts */
+ filename = "../MeshFiles/unittest/disk.h5m";
+ /* Options for reading
+ * - read in parallel
+ * - use PARALLEL_PARTITION tag
+ * - resolve shared entities after reading
+ */
+ options = "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARALLEL_RESOLVE_SHARED_ENTS";
+ }
+ else
+ {
+ filename = argv[1];
+ options = argv[2];
+ }
+ if (rank == 0)
+ std::cout << "reading file " << filename << "\n with options:" << options <<
+ "\n on " << nprocs << " processors\n";
+
+ // get MOAB instance and read the file with the specified options
+ Interface *mbImpl = new Core;
+ if (NULL == mbImpl) return 1;
+ ErrorCode rval = mbImpl->load_file(filename.c_str(), 0, options.c_str());
+ if (rval != MB_SUCCESS) return 1;
+
+ // get the ParallelComm instance
+ ParallelComm* pcomm = ParallelComm::get_pcomm(mbImpl, 0);
+ MPI_Comm comm = pcomm->comm();
+ if (0 == pcomm) return 1;
+
+ // get all shared entities with other processors
+ Range shared_ents;
+ rval = pcomm->get_shared_entities(-1, // -1 means all other processors Range &shared_ents,
+ shared_ents);
+ if (rval != MB_SUCCESS) return 1;
+ /* Among shared entities, get those owned by the current processor
+ * For this, use a filter operation;
+ * Each shared entity is owned by exactly one processor;
+ * An entity could be simply-shared (with exactly one other processor) or
+ * multi-shared.
+ */
+ Range owned_entities;
+ rval = pcomm->filter_pstatus(shared_ents, // pass entities that we want to filter
+ PSTATUS_NOT_OWNED, // status we are looking for
+ PSTATUS_NOT, // operation applied ; so it will return owned entities (!not_owned = owned)
+ -1, // this means all processors
+ &owned_entities);
+ if (rval != MB_SUCCESS) return 1;
+ unsigned int nums[4]={0}; // to store the owned entities per dimension
+ for (int i=0; i<3; i++)
+ {
+ nums[i]=(int)owned_entities.num_of_dimension(i);
+ }
+ int * rbuf;
+ if (rank==0)
+ rbuf = (int *)malloc(nprocs*4*sizeof(int));
+ MPI_Gather( nums, 4, MPI_INT, rbuf, 4, MPI_INT, 0, comm);
+ // print the stats gathered:
+ if (rank == 0)
+ {
+ for (int i=0; i<nprocs; i++)
+ {
+ std::cout << " shared, owned entities on proc " << i << " :" << rbuf[4*i] << " verts, " <<
+ rbuf[4*i+1] << " edges, " << rbuf[4*i+2] << " faces\n";
+ }
+
+ }
+
+ /*
+ * Now exchange 1 layer of ghost elements, using vertices as bridge
+ * we could have done this as part of reading process, by passing an extra read option
+ * ";PARALLEL_GHOSTS=2.0.1.0"
+ */
+ rval = pcomm->exchange_ghost_cells(2, // int ghost_dim,
+ 0, // int bridge_dim,
+ 1, //int num_layers,
+ 0, //int addl_ents,
+ true); // bool store_remote_handles);
+ if (rval != MB_SUCCESS) return 1;
+
+ // repeat the reports, after ghost exchange
+ shared_ents.clear();
+ owned_entities.clear();
+ rval = pcomm->get_shared_entities(-1, // -1 means all other processors Range &shared_ents,
+ shared_ents);
+ if (rval != MB_SUCCESS) return 1;
+ rval = pcomm->filter_pstatus(shared_ents,
+ PSTATUS_NOT_OWNED,
+ PSTATUS_NOT,
+ -1,
+ &owned_entities);
+ if (rval != MB_SUCCESS) return 1;
+
+ // find out how many shared entities of each dimension are owned on this processor
+ for (int i=0; i<3; i++)
+ nums[i]=(int)owned_entities.num_of_dimension(i);
+
+ // gather the statistics on processor 0
+ MPI_Gather( nums, 4, MPI_INT, rbuf, 4, MPI_INT, 0, comm);
+ if (rank == 0)
+ {
+ std::cout << " \n\n After exchanging one ghost layer: \n";
+ for (int i=0; i<nprocs; i++)
+ {
+ std::cout << " shared, owned entities on proc " << i << " :" << rbuf[4*i] << " verts, " <<
+ rbuf[4*i+1] << " edges, " << rbuf[4*i+2] << " faces\n";
+ }
+ free(rbuf);
+ }
+ MPI_Finalize();
+
+ return 0;
+}
diff --git a/examples/KDTree.cpp b/examples/KDTree.cpp
deleted file mode 100644
index f402ca1..0000000
--- a/examples/KDTree.cpp
+++ /dev/null
@@ -1,154 +0,0 @@
-/* Simple example of use of moab::AdaptiveKDTree class.
-
- Given a hexahedral mesh, find the hexahedron containing each
- input position.
- */
-
-
-#include "moab/Core.hpp"
-#include "moab/AdaptiveKDTree.hpp"
-#include "moab/Range.hpp"
-#include "moab/GeomUtil.hpp"
-
-#include <iostream>
-#include <string>
-
-const double EPSILON = 1e-6; // tolerance to use in intersection checks
-
-// Help with error handling. Given ErrorCode, print
-// corresponding string and any available message.
-void print_error( moab::Interface& mb, moab::ErrorCode err )
-{
- std::string message;
- std::string code;
- if (moab::MB_SUCCESS != mb.get_last_error( message ))
- message.clear();
- code = mb.get_error_string(err);
- std::cerr << "Error: " << code << std::endl;
- if (!message.empty())
- std::cerr << " " << message << std::endl;
-}
-
-// Print diagnostic info for unexpected failures.
-#define CHKERR(err) do { if (moab::MB_SUCCESS != (err)) { \
- print_error( mb, (err) ); \
- std::cerr << "Unexpected failure at: " << __FILE__ << ":" << __LINE__ << std::endl; \
- return 2; \
- } } while (false)
-
-// Given an entity set and a point, find the hex contained in the
-// entity set which in turn contains the specified point. Returns
-// 0 if point is not in any hexahedron.
-moab::EntityHandle hex_containing_point( moab::Interface& mb,
- moab::EntityHandle set,
- const double point[3] );
-
-// Print hex containing point.
-void print_hex( moab::Interface& mb, moab::EntityHandle hex );
-
-
-int main( )
-{
- // Ask user for file to read
- std::string filename;
- std::cout << "Hex mesh file name: ";
- std::cin >> filename;
-
- // Read file into MOAB instance
- moab::ErrorCode rval;
- moab::Core moab;
- moab::Interface& mb = moab;
- rval = mb.load_file( filename.c_str() );
- if (moab::MB_SUCCESS != rval) {
- print_error(mb,rval);
- std::cerr << filename << ": file load failed" << std::endl;
- return 1;
- }
-
- // Get all hex elemeents
- moab::Range elems;
- rval = mb.get_entities_by_type( 0, moab::MBHEX, elems ); CHKERR(rval);
- if (elems.empty()) {
- std::cerr << filename << ": file containd no hexahedra" << std::endl;
- return 1;
- }
-
- // Build a kD-tree from hex elements
- moab::EntityHandle tree_root;
- moab::AdaptiveKDTree tool( &mb );
- rval = tool.build_tree( elems, tree_root ); CHKERR(rval);
-
- // Loop forever (or until EOF), asking user for a point
- // to query and printing the hex element containing that
- // point.
- for (;;) {
- double point[3];
- std::cout << "Point coordinates: ";
- if (!(std::cin >> point[0] >> point[1] >> point[2]))
- break;
-
- moab::EntityHandle leaf;
- rval = tool.leaf_containing_point( tree_root, point, leaf ); CHKERR(rval);
- moab::EntityHandle hex = hex_containing_point( mb, leaf, point );
- if (0 == hex)
- std::cout << "Point is not contained in any hexahedron." << std::endl;
- else
- print_hex( mb, hex );
- }
-
- return 0;
-}
-
-moab::EntityHandle hex_containing_point( moab::Interface& mb,
- moab::EntityHandle set,
- const double point[3] )
-{
- moab::ErrorCode rval;
- moab::CartVect pt(point); // input location
- moab::CartVect coords[8]; // coordinates of corners of hexahedron
- const moab::EntityHandle* conn; // hex connectivity
- int conn_len;
-
- // Get hexes in leaf
- std::vector<moab::EntityHandle> hexes;
- rval = mb.get_entities_by_type( set, moab::MBHEX, hexes ); CHKERR(rval);
-
- // Check which hex the point is in
- std::vector<moab::EntityHandle>::const_iterator i;
- for (i = hexes.begin(); i != hexes.end(); ++i) {
- rval = mb.get_connectivity( *i, conn, conn_len ); CHKERR(rval);
- rval = mb.get_coords( conn, 8, &coords[0][0] ); CHKERR(rval);
- if (moab::GeomUtil::point_in_trilinear_hex( coords, pt, EPSILON ))
- return *i;
- }
-
- // Return 0 if no hex contains point.
- return 0;
-}
-
-void print_hex( moab::Interface& mb, moab::EntityHandle hex )
-{
- // Get MOAB's internal ID for hex element
- int id = mb.id_from_handle(hex);
-
- // Get vertex handles for hex corners
- const moab::EntityHandle* conn; // hex connectivity
- int conn_len;
- mb.get_connectivity( hex, conn, conn_len );
-
- // Get coordinates of vertices
- double coords[3*8];
- mb.get_coords( conn, 8, coords );
-
- // Print
- std::cout << " Point is in hex " << id << " with corners: " << std::endl;
- for (int i = 0; i < 8; ++i) {
- std::cout << " (" << coords[3*i]
- << ", " << coords[3*i+1]
- << ", " << coords[3*i+2]
- << ")" << std::endl;
- }
-}
-
-
-
diff --git a/examples/Makefile.am b/examples/Makefile.am
deleted file mode 100644
index d0a33b9..0000000
--- a/examples/Makefile.am
+++ /dev/null
@@ -1,68 +0,0 @@
-# NOTE: The 'makefile' in the *installed* examples directory is generated
-# using the install-data-hook target below. This logic assumes that:
-# 1) All examples are listed in check_PROGRAMS
-# 2) All examples have a single source file with the same name (case
-# sensitive) as the example executable with a .cpp suffix.
-
-check_PROGRAMS = FileRead \
- GeomSetHierarchy \
- GetEntities \
- KDTree \
- ObbTree \
- SetsNTags \
- SkinMesh \
- SurfArea
-
-SRCDIR = $(top_srcdir)/examples
-
-LDADD = ../src/libMOAB.la
-AM_CPPFLAGS += -I../src \
- -I$(srcdir)/../src \
- -I$(srcdir)/../src/parallel \
- -I$(srcdir)/../src/oldinc \
- -DSRCDIR=$(SRCDIR)
-
-FileRead_SOURCES = FileRead.cpp
-GeomSetHierarchy_SOURCES = GeomSetHierarchy.cpp
-GetEntities_SOURCES = simple/GetEntities.cpp
-SetsNTags_SOURCES = SetsNTags.cpp
-SkinMesh_SOURCES = SkinMesh.cpp
-SurfArea_SOURCES = SurfArea.cpp
-KDTree_SOURCES = KDTree.cpp
-ObbTree_SOURCES = ObbTree.cpp
-
-exampledir = ${docdir}/examples
-nobase_example_DATA = \
- examples.make \
- FileRead.cpp \
- GeomSetHierarchy.cpp \
- simple/GetEntities.cpp \
- simple/makefile \
- KDTree.cpp \
- ObbTree.cpp \
- SetsNTags.cpp \
- SkinMesh.cpp \
- SurfArea.cpp
-
-if ENABLE_imesh
- imesh_DIR = itaps
-else
- imesh_DIR =
-endif
-
-SUBDIRS = $(imesh_DIR)
-
-ex_make = $(DESTDIR)$(exampledir)/makefile
-install-data-hook:
- $(AM_V_at)rm -f $(ex_make)
- $(AM_V_at)mv $(DESTDIR)$(exampledir)/examples.make $(ex_make)
- $(AM_V_at)echo "all: $(check_PROGRAMS)" >>$(ex_make)
- $(AM_V_at)rule=' $$(CXX) -o $$@ $$< $$(CXXFLAGS) $$(MOAB_INCLUDES) $$(MOAB_LIBS_LINK)'; \
- for example in $(check_PROGRAMS); do \
- echo >>$(ex_make); \
- echo "$${example}: $${example}.cpp" >>$(ex_make); \
- echo "$$rule" >>$(ex_make); \
- done
-
-uninstall-hook:
- $(AM_V_at)rm -f $(DESTDIR)$(exampledir)/makefile
diff --git a/examples/ObbTree.cpp b/examples/ObbTree.cpp
deleted file mode 100644
index 220681f..0000000
--- a/examples/ObbTree.cpp
+++ /dev/null
@@ -1,80 +0,0 @@
-// simple example construct obb tree and ray-tracing the tree
-// it reads triangle mesh, construct obb tree and get intersection distances
-
-#include "moab/Core.hpp"
-#include "moab/Range.hpp"
-#include "moab/OrientedBoxTreeTool.hpp"
-#include <iostream>
-#include <math.h>
-
-int main(int argc, char **argv) {
- if (1 == argc) {
- std::cout << "Usage: " << argv[0] << " <filename>" << std::endl;
- return 0;
- }
-
- // instantiate & load a mesh from a file
- moab::Core *mb = new moab::Core();
- moab::ErrorCode rval = mb->load_mesh(argv[1]);
-
- // get all triangles
- moab::EntityHandle tree_root;
- moab::Range tris;
- //moab::OrientedBoxTreeTool::Settings settings;
-
- rval = mb->get_entities_by_type( 0, moab::MBTRI, tris );
- if (rval != moab::MB_SUCCESS) {
- std::cerr << "Couldn't get triangles." << std::endl;
- return 1;
- }
-
- // build OBB trees for all triangles
- moab::OrientedBoxTreeTool tool(mb);
- //rval = tool.build(tris, tree_root, &settings);
- rval = tool.build(tris, tree_root);
- if (rval != moab::MB_SUCCESS) {
- std::cerr << "Could'nt build tree." << std::endl;
- return 1;
- }
-
- // build box
- double box_center[3], box_axis1[3], box_axis2[3], box_axis3[3], pnt_start[3], ray_length;
- rval = tool.box(tree_root, box_center, box_axis1, box_axis2, box_axis3);
- if (rval != moab::MB_SUCCESS) {
- std::cerr << "Couldn't get box for tree root set.";
- return 1;
- }
-
- ray_length = 2.*sqrt(box_axis1[0]*box_axis1[0] + box_axis1[1]*box_axis1[1] +
- box_axis1[2]*box_axis1[2]);
-
- // do ray-tracing from box center side to x direction
- std::vector<double> intersections;
- std::vector<moab::EntityHandle> intersection_facets;
-
- for (int i=0; i<3; i++)
- pnt_start[i] = box_center[i]-box_axis1[i];
-
- if (ray_length>0) // normalize ray direction
- for (int j=0; j<3; j++)
- box_axis1[j]=2*box_axis1[j]/ray_length;
- rval = tool.ray_intersect_triangles(intersections, intersection_facets,
- tree_root, 10e-12, pnt_start, box_axis1,
- &ray_length);
- if (rval != moab::MB_SUCCESS) {
- std::cerr << "Couldn't ray tracing.";
- return 1;
- }
-
- std::cout << "ray start point: " << pnt_start[0] << " "
- << pnt_start[1] << " " << pnt_start[2] << std::endl;
- std::cout << " ray direction: " << box_axis1[0] << " " << box_axis1[1] << " " << box_axis1[2] << "\n";
- std::cout << "# of intersections : " << intersections.size() << std::endl;
- std::cout << "intersection distances are on";
- for (unsigned int i = 0; i < intersections.size(); i++) {
- std::cout << " " << intersections[i];
- }
- std::cout << " of ray length " << ray_length << std::endl;
-
- return 0;
-}
This diff is so big that we needed to truncate the remainder.
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