Mesh Oriented datABase  (version 5.4.1)
Array-based unstructured mesh datastructure
DirectAccessWithHoles.cpp
Go to the documentation of this file.
00001 /** @example DirectAccessWithHoles.cpp \n
00002  * \brief Use direct access to MOAB data to avoid calling through API \n
00003  *
00004  * This example creates a 1d row of quad elements, with a user-specified number of "holes" (missing
00005  * quads) in the row: \verbatim
00006  *  ----------------------      ----------------------      --------
00007  *  |      |      |      |      |      |      |      |      |      |
00008  *  |      |      |      |(hole)|      |      |      |(hole)|      | ...
00009  *  |      |      |      |      |      |      |      |      |      |
00010  *  ----------------------      ----------------------      --------
00011  * \endverbatim
00012  * This makes (nholes+1) contiguous runs of quad handles in the handle space
00013  * This example shows how to use the xxx_iterate functions in MOAB (xxx = coords, connect, tag,
00014  * adjacencies) to get direct pointer access to MOAB internal storage, which can be used without
00015  * calling through the MOAB API.
00016  *
00017  *    -#  Initialize MOAB \n
00018  *    -#  Create a quad mesh with holes, as depicted above
00019  *    -#  Create 2 dense tags (tag1, tag2) for avg position to assign to quads, and # verts per quad
00020  * (tag3)
00021  *    -#  Get connectivity, coordinate, tag1 iterators
00022  *    -#  Iterate through quads, computing midpoint based on vertex positions, set on quad-based
00023  * tag1
00024  *    -#  Set up map from starting quad handle for a chunk to struct of (tag1_ptr, tag2_ptr,
00025  * tag3_ptr), pointers to the dense tag storage for those tags for the chunk
00026  *    -#  Iterate through vertices, summing positions into tag2 on connected quads and incrementing
00027  * vertex count
00028  *    -#  Iterate through quads, normalizing tag2 by vertex count and comparing values of tag1 and
00029  * tag2
00030  *
00031  * <b>To compile</b>: \n
00032  *    make DirectAccessWithHoles \n
00033  * <b>To run</b>: ./DirectAccess [-nquads <# quads>] [-holes <# holes>]\n
00034  *
00035  */
00036 
00037 #include "moab/Core.hpp"
00038 #include "moab/ProgOptions.hpp"
00039 #include "moab/ReadUtilIface.hpp"
00040 #include <map>
00041 #include <iostream>
00042 #include <cassert>
00043 
00044 using namespace moab;
00045 using namespace std;
00046 
00047 ErrorCode create_mesh_with_holes( Interface* mbImpl, int nquads, int nholes );
00048 
00049 struct tag_struct
00050 {
00051     double* avg_ptr;
00052     int* nv_ptr;
00053 };
00054 
00055 int main( int argc, char** argv )
00056 {
00057     // Get MOAB instance
00058     Interface* mbImpl = new( std::nothrow ) Core;
00059     if( NULL == mbImpl ) return 1;
00060 
00061     int nquads = 1000, nholes = 1;
00062 
00063     // Parse options
00064     ProgOptions opts;
00065     opts.addOpt< int >( string( "nquads,n" ), string( "Number of quads in the mesh (default = 1000" ), &nquads );
00066     opts.addOpt< int >( string( "holes,H" ), string( "Number of holes in the element handle space (default = 1" ),
00067                         &nholes );
00068     opts.parseCommandLine( argc, argv );
00069     if( nholes >= nquads )
00070     {
00071         cerr << "Number of holes needs to be < number of elements." << endl;
00072         return 1;
00073     }
00074 
00075     // Create simple structured mesh with hole, but using unstructured representation
00076     ErrorCode rval = create_mesh_with_holes( mbImpl, nquads, nholes );MB_CHK_SET_ERR( rval, "Trouble creating mesh" );
00077 
00078     // Get all vertices and non-vertex entities
00079     Range verts, elems;
00080     rval = mbImpl->get_entities_by_handle( 0, elems );MB_CHK_SET_ERR( rval, "Trouble getting all entities" );
00081     verts = elems.subset_by_dimension( 0 );
00082     elems -= verts;
00083 
00084     // Create tag1 (element-based avg), tag2 (vertex-based avg), tag3 (# connected verts)
00085     Tag tag1, tag2, tag3;
00086     rval = mbImpl->tag_get_handle( "tag1", 3, MB_TYPE_DOUBLE, tag1, MB_TAG_DENSE | MB_TAG_CREAT );MB_CHK_SET_ERR( rval, "Trouble creating tag1" );
00087     double def_val[3] = { 0.0, 0.0, 0.0 };  // Need a default value for tag2 because we sum into it
00088     rval              = mbImpl->tag_get_handle( "tag2", 3, MB_TYPE_DOUBLE, tag2, MB_TAG_DENSE | MB_TAG_CREAT, def_val );MB_CHK_SET_ERR( rval, "Trouble creating tag2" );
00089     int def_val_int = 0;  // Need a default value for tag3 because we increment it
00090     rval = mbImpl->tag_get_handle( "tag3", 1, MB_TYPE_INTEGER, tag3, MB_TAG_DENSE | MB_TAG_CREAT, &def_val_int );MB_CHK_SET_ERR( rval, "Trouble creating tag3" );
00091 
00092     // Get connectivity, coordinate, tag, and adjacency iterators
00093     EntityHandle* conn_ptr;
00094     double *x_ptr, *y_ptr, *z_ptr, *tag1_ptr, *tag2_ptr;
00095     int* tag3_ptr;
00096 
00097     // First vertex is at start of range (ranges are sorted), and is offset for vertex index
00098     // calculation
00099     EntityHandle first_vert = *verts.begin();
00100 
00101     // When iterating over elements, each chunk can have a different # vertices; also, count tells
00102     // you how many elements are in the current chunk
00103     int vpere, count;
00104 
00105     // Get coordinates iterator, just need this once because we know verts handle space doesn't have
00106     // holes
00107     rval = mbImpl->coords_iterate( verts.begin(), verts.end(), x_ptr, y_ptr, z_ptr, count );MB_CHK_SET_ERR( rval, "Error in coords_iterate" );
00108     assert( count == (int)verts.size() );  // Should end up with just one contiguous chunk of vertices
00109 
00110     // Iterate through elements, computing midpoint based on vertex positions, set on element-based
00111     // tag1 Control loop by iterator over elem range
00112     Range::iterator e_it = elems.begin();
00113 
00114     while( e_it != elems.end() )
00115     {
00116         // Get conn_ptr, tag1_ptr for next contiguous chunk of element handles, and coords pointers
00117         // for all verts
00118         rval = mbImpl->connect_iterate( e_it, elems.end(), conn_ptr, vpere, count );MB_CHK_SET_ERR( rval, "Error in connect_iterate" );
00119         rval = mbImpl->tag_iterate( tag1, e_it, elems.end(), count, (void*&)tag1_ptr );MB_CHK_SET_ERR( rval, "Error in tag1_iterate" );
00120 
00121         // Iterate over elements in this chunk
00122         for( int i = 0; i < count; i++ )
00123         {
00124             tag1_ptr[0] = tag1_ptr[1] = tag1_ptr[2] = 0.0;  // Initialize position for this element
00125             for( int j = 0; j < vpere; j++ )
00126             {                                            // Loop over vertices in this element
00127                 int v_index = conn_ptr[j] - first_vert;  // vert index is just the offset from first vertex
00128                 tag1_ptr[0] += x_ptr[v_index];
00129                 tag1_ptr[1] += y_ptr[v_index];  // Sum vertex positions into tag1...
00130                 tag1_ptr[2] += z_ptr[v_index];
00131             }
00132             tag1_ptr[0] /= vpere;
00133             tag1_ptr[1] /= vpere;  // Then normalize
00134             tag1_ptr[2] /= vpere;
00135 
00136             // Done with this element; advance connect_ptr and tag1_ptr to next element
00137             conn_ptr += vpere;
00138             tag1_ptr += 3;
00139         }  // Loop over elements in chunk
00140 
00141         // Done with chunk; advance range iterator by count; will skip over gaps in range
00142         e_it += count;
00143     }  // While loop over all elements
00144 
00145     // Iterate through vertices, summing positions into tag2 on connected elements and incrementing
00146     // vertex count Iterate over chunks the same as elements, even though we know we have only one
00147     // chunk here, just to show how it's done
00148 
00149     // Create a std::map from EntityHandle (first entity handle in chunk) to
00150     // tag_struct (ptrs to start of avg/#verts tags for that chunk); then for a given entity handle,
00151     // we can quickly find the chunk it's in using map::lower_bound; could have set up this map in
00152     // earlier loop over elements, but do it here for clarity
00153 
00154     map< EntityHandle, tag_struct > elem_map;
00155     e_it = elems.begin();
00156     while( e_it != elems.end() )
00157     {
00158         tag_struct ts = { NULL, NULL };
00159         rval          = mbImpl->tag_iterate( tag2, e_it, elems.end(), count, (void*&)ts.avg_ptr );MB_CHK_SET_ERR( rval, "Error in tag2_iterate" );
00160         rval = mbImpl->tag_iterate( tag3, e_it, elems.end(), count, (void*&)ts.nv_ptr );MB_CHK_SET_ERR( rval, "Error in tag3_iterate" );
00161         elem_map[*e_it] = ts;
00162         e_it += count;
00163     }
00164 
00165     // Call a vertex-quad adjacencies function to generate vertex-element adjacencies in MOAB
00166     Range::iterator v_it = verts.begin();
00167     Range dum_range;
00168     rval = mbImpl->get_adjacencies( &( *v_it ), 1, 2, false, dum_range );MB_CHK_SET_ERR( rval, "Error in get_adjacencies" );
00169     const vector< EntityHandle >** adjs_ptr;
00170     while( v_it != verts.end() )
00171     {
00172         // Get coords ptrs, adjs_ptr; can't set tag2_ptr by direct access, because of hole in
00173         // element handle space
00174         rval = mbImpl->coords_iterate( v_it, verts.end(), x_ptr, y_ptr, z_ptr, count );MB_CHK_SET_ERR( rval, "Error in coords_iterate" );
00175         rval = mbImpl->adjacencies_iterate( v_it, verts.end(), adjs_ptr, count );MB_CHK_SET_ERR( rval, "Error in adjacencies_iterate" );
00176 
00177         for( int v = 0; v < count; v++ )
00178         {
00179             const vector< EntityHandle >* avec = *( adjs_ptr + v );
00180             for( vector< EntityHandle >::const_iterator ait = avec->begin(); ait != avec->end(); ++ait )
00181             {
00182                 // Get chunk that this element resides in; upper_bound points to the first element
00183                 // strictly > key, so get that and decrement (would work the same as lower_bound
00184                 // with an if-test and decrement)
00185                 map< EntityHandle, tag_struct >::iterator mit = elem_map.upper_bound( *ait );
00186                 --mit;
00187                 // Index of *ait in that chunk
00188                 int a_ind     = *ait - ( *mit ).first;
00189                 tag_struct ts = ( *mit ).second;
00190                 ts.avg_ptr[3 * a_ind + 0] += x_ptr[v];  // Tag on each element is 3 doubles, x/y/z
00191                 ts.avg_ptr[3 * a_ind + 1] += y_ptr[v];
00192                 ts.avg_ptr[3 * a_ind + 2] += z_ptr[v];
00193                 ts.nv_ptr[a_ind]++;  // Increment the vertex count
00194             }
00195         }
00196 
00197         v_it += count;
00198     }
00199 
00200     // Normalize tag2 by vertex count; loop over elements using same approach as before
00201     // At the same time, compare values of tag1 and tag2
00202     e_it = elems.begin();
00203     while( e_it != elems.end() )
00204     {
00205         // Get conn_ptr, tag1_ptr for next contiguous chunk of element handles, and coords pointers
00206         // for all verts
00207         rval = mbImpl->tag_iterate( tag1, e_it, elems.end(), count, (void*&)tag1_ptr );MB_CHK_SET_ERR( rval, "Error in tag1_iterate" );
00208         rval = mbImpl->tag_iterate( tag2, e_it, elems.end(), count, (void*&)tag2_ptr );MB_CHK_SET_ERR( rval, "Error in tag2_iterate" );
00209         rval = mbImpl->tag_iterate( tag3, e_it, elems.end(), count, (void*&)tag3_ptr );MB_CHK_SET_ERR( rval, "Error in tag3_iterate" );
00210 
00211         // Iterate over elements in this chunk
00212         for( int i = 0; i < count; i++ )
00213         {
00214             for( int j = 0; j < 3; j++ )
00215                 tag2_ptr[3 * i + j] /= (double)tag3_ptr[i];  // Normalize by # verts
00216             if( tag1_ptr[3 * i] != tag2_ptr[3 * i] || tag1_ptr[3 * i + 1] != tag2_ptr[3 * i + 1] ||
00217                 tag1_ptr[3 * i + 2] != tag2_ptr[3 * i + 2] )
00218                 cout << "Tag1, tag2 disagree for element " << *e_it + i << endl;
00219         }
00220 
00221         e_it += count;
00222     }
00223 
00224     // Ok, we're done, shut down MOAB
00225     delete mbImpl;
00226 
00227     return 0;
00228 }
00229 
00230 ErrorCode create_mesh_with_holes( Interface* mbImpl, int nquads, int nholes )
00231 {
00232     // First make the mesh, a 1d array of quads with left hand side x = elem_num; vertices are
00233     // numbered in layers
00234     ReadUtilIface* read_iface;
00235     ErrorCode rval = mbImpl->query_interface( read_iface );MB_CHK_SET_ERR( rval, "Error in query_interface" );
00236     vector< double* > coords;
00237     EntityHandle start_vert, start_elem, *connect;
00238     // Create verts, num is 4(nquads+1) because they're in a 1d row; will initialize coords in loop
00239     // over elems later
00240     rval = read_iface->get_node_coords( 3, 2 * ( nquads + 1 ), 0, start_vert, coords );MB_CHK_SET_ERR( rval, "Error in get_node_arrays" );
00241     // Create elems
00242     rval = read_iface->get_element_connect( nquads, 4, MBQUAD, 0, start_elem, connect );MB_CHK_SET_ERR( rval, "Error in get_element_connect" );
00243     for( int i = 0; i < nquads; i++ )
00244     {
00245         coords[0][2 * i] = coords[0][2 * i + 1] = (double)i;  // x values are all i
00246         coords[1][2 * i]                        = 0.0;
00247         coords[1][2 * i + 1]                    = 1.0;          // y coords
00248         coords[2][2 * i] = coords[2][2 * i + 1] = (double)0.0;  // z values, all zero (2d mesh)
00249         EntityHandle quad_v                     = start_vert + 2 * i;
00250         for( int j = 0; j < 4; j++ )
00251             connect[4 * i + j] = quad_v + j;  // Connectivity of each quad is a sequence starting from quad_v
00252     }
00253     // Last two vertices
00254     // Cppcheck warning (false positive): variable coords is assigned a value that is never used
00255     coords[0][2 * nquads] = coords[0][2 * nquads + 1] = (double)nquads;
00256     coords[1][2 * nquads]                             = 0.0;
00257     coords[1][2 * nquads + 1]                         = 1.0;          // y coords
00258     coords[2][2 * nquads] = coords[2][2 * nquads + 1] = (double)0.0;  // z values, all zero (2d mesh)
00259 
00260     // Now delete nholes elements, spaced approximately equally through mesh, so contiguous size is
00261     // about nquads/(nholes + 1) reinterpret start_elem as the next element to be deleted
00262     int de = nquads / ( nholes + 1 );
00263     for( int i = 0; i < nholes; i++ )
00264     {
00265         start_elem += de;
00266         rval = mbImpl->delete_entities( &start_elem, 1 );MB_CHK_SET_ERR( rval, "Error in delete_entities" );
00267     }
00268 
00269     return MB_SUCCESS;
00270 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines