[MOAB-dev] commit/MOAB: 5 new changesets
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Tue May 7 11:01:14 CDT 2013
5 new commits in MOAB:
https://bitbucket.org/fathomteam/moab/commits/e38337a1e926/
Changeset: e38337a1e926
Branch: None
User: iulian07
Date: 2013-04-24 21:23:31
Summary: Merge branch 'master' of bitbucket.org:fathomteam/moab
Affected #: 5 files
diff --git a/examples/DirectAccessNoHoles.cpp b/examples/DirectAccessNoHoles.cpp
new file mode 100644
index 0000000..445c517
--- /dev/null
+++ b/examples/DirectAccessNoHoles.cpp
@@ -0,0 +1,179 @@
+/** @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_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
+ 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;
+ }
+
+ // 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 4(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;
+ 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)
+
+ // 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/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/makefile b/examples/makefile
index 4ec6154..fd6585e 100644
--- a/examples/makefile
+++ b/examples/makefile
@@ -21,6 +21,12 @@ HelloMoabPar: HelloMoabPar.o
TestExodusII: TestExodusII.o
${MOAB_CXX} -o $@ $< ${MOAB_LIBS_LINK}
+DirectAccessWithHoles: DirectAccessWithHoles.o
+ ${MOAB_CXX} -o $@ $< ${MOAB_LIBS_LINK}
+
+DirectAccessNoHoles: DirectAccessNoHoles.o
+ ${MOAB_CXX} -o $@ $< ${MOAB_LIBS_LINK}
+
.cpp.o :
${MOAB_CXX} ${MOAB_CXXFLAGS} ${MOAB_INCLUDES} -DMESH_DIR=\"${MESH_DIR}\" -c $<
diff --git a/src/io/ReadNC.cpp b/src/io/ReadNC.cpp
index d7bcb35..774f7a6 100644
--- a/src/io/ReadNC.cpp
+++ b/src/io/ReadNC.cpp
@@ -2321,6 +2321,10 @@ ErrorCode ReadNC::init_FVCDscd_vals(const FileOptions &opts, EntityHandle file_s
dbgOut.tprintf(2, "Tag created for variable %s\n", tag_name.c_str());
}
+ // hack: create dummy tags, if needed, for dimensions like nbnd
+ // with no corresponding variables
+ init_dims_with_no_cvars_info();
+
return MB_SUCCESS;
}
@@ -2699,9 +2703,39 @@ ErrorCode ReadNC::init_EulSpcscd_vals(const FileOptions &opts, EntityHandle file
dbgOut.tprintf(2, "Tag created for variable %s\n", tag_name.c_str());
}
+ // hack: create dummy tags, if needed, for variables like nbnd
+ // with no corresponding variables
+ init_dims_with_no_cvars_info();
+
return MB_SUCCESS;
}
+void ReadNC::init_dims_with_no_cvars_info() {
+ // hack: look at all dimensions, and see if we have one that does not appear in the list of varInfo names
+ // right now, candidates are ncol and nbnd
+ // for them, create dummy tags
+ for (unsigned int i=0; i<dimNames.size(); i++)
+ {
+ // if there is a var with this name, skip, we are fine; if not, create a varInfo...
+ if ( varInfo.find(dimNames[i])!=varInfo.end())
+ continue; // we already have a variable with this dimension name
+
+ int sizeTotalVar = varInfo.size();
+ std::string var_name(dimNames[i]);
+ VarData &data = varInfo[var_name];
+ data.varName = std::string(var_name);
+ data.varId =sizeTotalVar;
+ data.varTags.resize(1, 0);
+ data.varDataType = NC_DOUBLE; // could be int, actually, but we do not really need the type
+ data.varDims.resize(1);
+ data.varDims[0]= (int)i;
+ data.numAtts=0;
+ data.entLoc = ENTLOCSET;
+ dbgOut.tprintf(2, "Dummy varInfo created for dimension %s\n", dimNames[i].c_str());
+ dummyVarNames.insert(dimNames[i]);
+ }
+}
+
ErrorCode ReadNC::init_HOMMEucd_vals() {
ErrorCode rval;
unsigned int idx;
@@ -2812,29 +2846,10 @@ ErrorCode ReadNC::init_HOMMEucd_vals() {
// don't read coordinates of columns until we actually create the mesh
- // hack: look at all dimensions, and see if we have one that does not appear in the list of varInfo names
- // right now, candidates are ncol and nbnd
- // for them, create dummy tags
- for (unsigned int i=0; i<dimNames.size(); i++)
- {
- // if there is a var with this name, skip, we are fine; if not, create a varInfo...
- if ( varInfo.find(dimNames[i])!=varInfo.end())
- continue; // we already have a variable with this dimension name
+ // hack: create dummy tags, if needed, for variables like ncol and nbnd
+ // with no corresponding variables
+ init_dims_with_no_cvars_info();
- int sizeTotalVar = varInfo.size();
- std::string var_name(dimNames[i]);
- VarData &data = varInfo[var_name];
- data.varName = std::string(var_name);
- data.varId =sizeTotalVar;
- data.varTags.resize(1, 0);
- data.varDataType = NC_DOUBLE; // could be int, actually, but we do not really need the type
- data.varDims.resize(1);
- data.varDims[0]= (int)i;
- data.numAtts=0;
- data.entLoc = ENTLOCSET;
- dbgOut.tprintf(2, "Dummy varInfo created for dimension %s\n", dimNames[i].c_str());
- dummyVarNames.insert(dimNames[i]);
- }
return MB_SUCCESS;
}
diff --git a/src/io/ReadNC.hpp b/src/io/ReadNC.hpp
index d5a18eb..7bc07d5 100644
--- a/src/io/ReadNC.hpp
+++ b/src/io/ReadNC.hpp
@@ -236,6 +236,10 @@ private:
//! create COORDS tag for quads coordinate
ErrorCode create_quad_coordinate_tag(EntityHandle file_set);
+ //! Init info for dimensions that don't have corresponding
+ //! coordinate variables - this info is used for creating tags
+ void init_dims_with_no_cvars_info();
+
ErrorCode init_HOMMEucd_vals();
ErrorCode create_ucd_verts_quads(const FileOptions &opts, EntityHandle tmp_set, Range &quads);
https://bitbucket.org/fathomteam/moab/commits/dd971866630a/
Changeset: dd971866630a
Branch: None
User: iulian07
Date: 2013-04-25 19:59:20
Summary: Merge branch 'master' of bitbucket.org:iulian07/moab
Affected #: 0 files
https://bitbucket.org/fathomteam/moab/commits/12a0b8de1e11/
Changeset: 12a0b8de1e11
Branch: None
User: iulian07
Date: 2013-05-07 16:16:22
Summary: Merge branch 'master' of bitbucket.org:iulian07/moab
Affected #: 9 files
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 66fcf81..1524b06 100644
--- a/configure.ac
+++ b/configure.ac
@@ -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])
diff --git a/examples/DirectAccessNoHoles.cpp b/examples/DirectAccessNoHoles.cpp
index 445c517..340efae 100644
--- a/examples/DirectAccessNoHoles.cpp
+++ b/examples/DirectAccessNoHoles.cpp
@@ -1,31 +1,30 @@
-/** @example DirectAccess.cpp \n
+/** @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, with a user-specified number of "holes" (missing quads) in the row:
+ * 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.
*
- * ---------------------- ---------------------- --------
- * | | | | | | | | | |
- * | | | |(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 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
- * -# 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
+ * make DirectAccessNoHoles MOAB_DIR=<installdir> \n
+ * <b>To run</b>: ./DirectAccessNoHoles [-nquads <# quads>]\n
*
*/
@@ -131,12 +130,16 @@ int main(int argc, char **argv)
// 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])
+ 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;
@@ -151,7 +154,7 @@ ErrorCode create_mesh_no_holes(Interface *mbImpl, int nquads)
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 quads later
+ // 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");
@@ -160,8 +163,12 @@ ErrorCode create_mesh_no_holes(Interface *mbImpl, int nquads)
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
+ 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
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/makefile b/examples/makefile
index fd6585e..3402f96 100644
--- a/examples/makefile
+++ b/examples/makefile
@@ -1,7 +1,8 @@
# MOAB_DIR points to top-level install dir, below which MOAB's lib/ and include/ are located
include ${MOAB_DIR}/lib/moab.make
+include ${MOAB_DIR}/lib/iMesh-Defs.inc
-.SUFFIXES: .o .cpp
+.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="../MeshFiles/unittest"
@@ -27,6 +28,12 @@ DirectAccessWithHoles: DirectAccessWithHoles.o
DirectAccessNoHoles: DirectAccessNoHoles.o
${MOAB_CXX} -o $@ $< ${MOAB_LIBS_LINK}
+DirectAccessNoHolesF90: DirectAccessNoHolesF90.o
+ ${MOAB_CXX} -o $@ $< ${IMESH_LIBS}
+
.cpp.o :
${MOAB_CXX} ${MOAB_CXXFLAGS} ${MOAB_INCLUDES} -DMESH_DIR=\"${MESH_DIR}\" -c $<
+.F90.o :
+ ${IMESH_FC} ${IMESH_FCFLAGS} ${IMESH_INCLUDES} ${IMESH_FCDEFS} -DMESH_DIR=\"${MESH_DIR}\" -c $<
+
diff --git a/itaps/iBase_f.h.in b/itaps/iBase_f.h.in
index 4a70fd3..90fde77 100644
--- a/itaps/iBase_f.h.in
+++ b/itaps/iBase_f.h.in
@@ -50,7 +50,27 @@
parameter (iBase_REGION = 3)
parameter (iBase_ALL_TYPES = 4)
-
+ integer iBase_AdjacencyCost_MIN
+ integer iBase_UNAVAILABLE
+ integer iBase_ALL_ORDER_1
+ integer iBase_ALL_ORDER_LOGN
+ integer iBase_ALL_ORDER_N
+ integer iBase_SOME_ORDER_1
+ integer iBase_SOME_ORDER_LOGN
+ integer iBase_SOME_ORDER_N
+ integer iBase_AVAILABLE
+ integer iBase_AdjacencyCost_MAX
+
+ parameter (iBase_AdjacencyCost_MIN = 0)
+ parameter (iBase_UNAVAILABLE = 0)
+ parameter (iBase_ALL_ORDER_1 = 1)
+ parameter (iBase_ALL_ORDER_LOGN = 2)
+ parameter (iBase_ALL_ORDER_N = 3)
+ parameter (iBase_SOME_ORDER_1 = 4)
+ parameter (iBase_SOME_ORDER_LOGN = 5)
+ parameter (iBase_SOME_ORDER_N = 6)
+ parameter (iBase_AVAILABLE = 7)
+ parameter (iBase_AdjacencyCost_MAX = 7)
integer iBase_NEW
integer iBase_ALREADY_EXISTED
diff --git a/itaps/imesh/Makefile.am b/itaps/imesh/Makefile.am
index 647b8c1..97c76bb 100644
--- a/itaps/imesh/Makefile.am
+++ b/itaps/imesh/Makefile.am
@@ -45,6 +45,9 @@ if PARALLEL
# ftest_DEPENDENCIES = libiMesh.la $(top_builddir)/libMOAB.la
endif
+FCLINK = $(CXXLINK)
+F77LINK = $(CXXLINK)
+
TESTS = $(check_PROGRAMS)
LDADD = libiMesh.la $(top_builddir)/src/libMOAB.la ${MOAB_CXX_LINKFLAGS} ${MOAB_CXX_LIBS}
TESTDEPS = libiMesh.la $(top_builddir)/src/libMOAB.la
diff --git a/moab.make.in b/moab.make.in
index 1e294b3..e6dc42d 100644
--- a/moab.make.in
+++ b/moab.make.in
@@ -10,12 +10,16 @@ MOAB_INCLUDES = -I at abs_srcdir@/src \
MOAB_CXXFLAGS = @CXXFLAGS@ @AM_CXXFLAGS@
MOAB_CFLAGS = @CFLAGS@ @AM_CFLAGS@
+MOAB_FFLAGS = @FFLAGS@
+MOAB_FCFLAGS = @FCFLAGS@
MOAB_LDFLAGS = @EXPORT_LDFLAGS@
MOAB_LIBS_LINK = ${MOAB_LDFLAGS} -L${MOAB_LIBDIR} -lMOAB @LIBS@ @PNETCDF_LIBS@ @NETCDF_LIBS@ @HDF5_LIBS@ @CCMIO_LIBS@ @CGM_LIBS@
MOAB_CXX = @CXX@
MOAB_CC = @CC@
+MOAB_FC = @FC@
+MOAB_F77 = @F77@
# Override MOAB_LIBDIR and MOAB_INCLUDES from above with the correct
# values for the installed MOAB.
diff --git a/src/parallel/WriteHDF5Parallel.cpp b/src/parallel/WriteHDF5Parallel.cpp
index da16c7b..e96d933 100644
--- a/src/parallel/WriteHDF5Parallel.cpp
+++ b/src/parallel/WriteHDF5Parallel.cpp
@@ -1680,23 +1680,17 @@ static void convert_to_ranged_ids( const TYPE* buffer,
}
result.resize( len*2 );
- std::copy( buffer, buffer+len, result.begin()+len );
- std::sort( result.begin()+len, result.end() );
- std::vector<WriteHDF5::id_t>::iterator w = result.begin();
- std::vector<WriteHDF5::id_t>::const_iterator r = result.begin()+len;
- *w = *r; ++w;
- *w = *r; ++r;
- for (; r != result.end(); ++r) {
- if (*r == *w + 1) {
- ++*w;
- }
- else {
- ++w; *w = *r;
- ++w; *w = *r;
- }
+ Range tmp;
+ for (size_t i=0; i<len; i++)
+ tmp.insert( (EntityHandle)buffer[i]);
+ result.resize(tmp.psize()*2);
+ int j=0;
+ for (Range::const_pair_iterator pit=tmp.const_pair_begin();
+ pit!=tmp.const_pair_end(); pit++, j++)
+ {
+ result[2*j]=pit->first;
+ result[2*j+1]=pit->second-pit->first+1;
}
- ++w;
- result.erase( w, result.end() );
}
static void merge_ranged_ids( const unsigned long* range_list,
@@ -1710,22 +1704,21 @@ static void merge_ranged_ids( const unsigned long* range_list,
result.insert( result.end(), range_list, range_list+len );
size_t plen = result.size()/2;
- if (plen > 1) {
- std::pair<id_t,id_t> *plist = reinterpret_cast<std::pair<id_t,id_t>*>(&result[0]);
- std::sort( plist, plist+plen );
- std::pair<id_t,id_t> *wr = plist, *pend = plist+plen;
- std::sort( plist, pend );
- for (std::pair<id_t,id_t> *iter = plist+1; iter != pend; ++iter) {
- if (wr->second+1 < iter->first) {
- ++wr;
- *wr = *iter;
- }
- else {
- wr->second = std::max( wr->second, iter->second );
- }
- }
- ++wr;
- result.resize( 2*(wr - plist) );
+ Range tmp;
+ for (size_t i=0; i<plen; i++)
+ {
+ EntityHandle starth=(EntityHandle)result[2*i];
+ EntityHandle endh=(EntityHandle)result[2*i]+(id_t)result[2*i+1]-1; // id+count-1
+ tmp.insert( starth, endh);
+ }
+ // now convert back to std::vector<WriteHDF5::id_t>, compressed range format
+ result.resize(tmp.psize()*2);
+ int j=0;
+ for (Range::const_pair_iterator pit=tmp.const_pair_begin();
+ pit!=tmp.const_pair_end(); pit++, j++)
+ {
+ result[2*j]=pit->first;
+ result[2*j+1]=pit->second-pit->first+1;
}
}
@@ -1767,7 +1760,7 @@ ErrorCode WriteHDF5Parallel::unpack_set( EntityHandle set,
if (flags & mhdf_SET_RANGE_BIT) {
tmp = data->contentIds;
convert_to_ranged_ids( &tmp[0], tmp.size(), data->contentIds );
- data->setFlags &= mhdf_SET_RANGE_BIT;
+ data->setFlags |= mhdf_SET_RANGE_BIT;
}
else {
tmp.clear();
https://bitbucket.org/fathomteam/moab/commits/a486afb25464/
Changeset: a486afb25464
Branch: None
User: iulian07
Date: 2013-05-07 17:54:42
Summary: move the cslam_par_test driver to cslam tools folder
Affected #: 2 files
Diff not available.
https://bitbucket.org/fathomteam/moab/commits/c90ca7aadbc5/
Changeset: c90ca7aadbc5
Branch: master
User: iulian07
Date: 2013-05-07 18:00:11
Summary: complete the move correctly, for cslam_par_test
Affected #: 2 files
Diff not available.
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