DirectAccessWithHoles.cpp

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.

1. Initialize MOAB
   
2. Create a quad mesh with holes, as depicted above
3. Create 2 dense tags (tag1, tag2) for avg position to assign to quads, and # verts per quad (tag3)
4. Get connectivity, coordinate, tag1 iterators
5. Iterate through quads, computing midpoint based on vertex positions, set on quad-based tag1
6. 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
7. Iterate through vertices, summing positions into tag2 on connected quads and incrementing vertex count
8. Iterate through quads, normalizing tag2 by vertex count and comparing values of tag1 and tag2

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;
}