|
moab
|
#include "moab/Core.hpp"#include "moab/ProgOptions.hpp"#include "moab/ReadUtilIface.hpp"#include <map>#include <iostream>#include <assert.h>Go to the source code of this file.
Classes | |
| struct | tag_struct |
Defines | |
| #define | CHKERR(CODE, MSG) |
Functions | |
| ErrorCode | create_mesh_with_holes (Interface *mbImpl, int nquads, int nholes) |
| int | main (int argc, char **argv) |
| #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)
Definition at line 40 of file DirectAccessWithHoles.cpp.
| ErrorCode create_mesh_with_holes | ( | Interface * | mbImpl, |
| int | nquads, | ||
| int | nholes | ||
| ) |
Definition at line 204 of file DirectAccessWithHoles.cpp.
{
// 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;
}
| int main | ( | int | argc, |
| char ** | argv | ||
| ) |
Definition at line 49 of file DirectAccessWithHoles.cpp.
{
// 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;
}
1.7.6.1