DirectAccessWithHoles.cpp

Go to the documentation of this file.

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