moab
|
Use direct access to MOAB data to avoid calling through API
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.
To compile:
make DirectAccessWithHoles MOAB_DIR=<installdir>
To run: ./DirectAccess [-nquads <# quads>] [-holes <# holes>]
#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; }