Mesh Oriented datABase  (version 5.4.1)
Array-based unstructured mesh datastructure
DirectAccessNoHoles.cpp
Go to the documentation of this file.
00001 /** @example DirectAccessNoHoles.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, such that all quad and vertex handles
00005  * are contiguous in the handle space and in the database.  Then it shows how to get access
00006  * to pointers to MOAB-native data for vertex coordinates, quad connectivity, tag storage,
00007  * and vertex to quad adjacency lists.  This allows applications to access this data directly
00008  * without going through MOAB's API.  In cases where the mesh is not changing (or only mesh
00009  * vertices are moving), this can save significant execution time in applications.
00010  * \verbatim
00011  *  ----------------------
00012  *  |      |      |      |
00013  *  |      |      |      | ...
00014  *  |      |      |      |
00015  *  ----------------------
00016  * \endverbatim
00017  *    -#  Initialize MOAB \n
00018  *    -#  Create a quad mesh, 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  *    -#  Iterate through vertices, summing positions into tag2 on connected quads and incrementing
00025  * vertex count
00026  *    -#  Iterate through quads, normalizing tag2 by vertex count and comparing values of tag1 and
00027  * tag2
00028  *
00029  * <b>To compile</b>: \n
00030  *    make DirectAccessNoHoles \n
00031  * <b>To run</b>: ./DirectAccessNoHoles [-nquads <# quads>]\n
00032  *
00033  */
00034 
00035 #include "moab/Core.hpp"
00036 #include "moab/ProgOptions.hpp"
00037 #include "moab/ReadUtilIface.hpp"
00038 #include <map>
00039 #include <iostream>
00040 #include <cassert>
00041 
00042 using namespace moab;
00043 using namespace std;
00044 
00045 ErrorCode create_mesh_no_holes( Interface* mbImpl, int nquads );
00046 
00047 int main( int argc, char** argv )
00048 {
00049     // Get MOAB instance
00050     Interface* mbImpl = new( std::nothrow ) Core;
00051     if( NULL == mbImpl ) return 1;
00052 
00053     int nquads = 1000;
00054 
00055     // Parse options
00056     ProgOptions opts;
00057     opts.addOpt< int >( string( "nquads,n" ), string( "Number of quads in the mesh (default = 1000" ), &nquads );
00058     opts.parseCommandLine( argc, argv );
00059 
00060     // Create simple structured mesh with hole, but using unstructured representation
00061     ErrorCode rval = create_mesh_no_holes( mbImpl, nquads );MB_CHK_SET_ERR( rval, "Trouble creating mesh" );
00062 
00063     // Get all vertices and non-vertex entities
00064     Range verts, quads;
00065     rval = mbImpl->get_entities_by_handle( 0, quads );MB_CHK_SET_ERR( rval, "Trouble getting all entities" );
00066     verts = quads.subset_by_dimension( 0 );
00067     quads -= verts;
00068 
00069     // Create tag1 (element-based avg), tag2 (vertex-based avg), tag3 (# connected verts)
00070     Tag tag1, tag2, tag3;
00071     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" );
00072     double def_val[3] = { 0.0, 0.0, 0.0 };  // need a default value for tag2 because we sum into it
00073     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" );
00074     int def_val_int = 0;  // need a default value for tag3 because we increment it
00075     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" );
00076 
00077     // Get pointers to connectivity, coordinate, tag, and adjacency arrays; each of these returns a
00078     // count, which should be compared to the # entities you expect to verify there's only one chunk
00079     // (no holes)
00080     int count, vpere;
00081     EntityHandle* conn_ptr;
00082     rval = mbImpl->connect_iterate( quads.begin(), quads.end(), conn_ptr, vpere, count );MB_CHK_SET_ERR( rval, "Error in connect_iterate" );
00083     assert( count == (int)quads.size() );  // Should end up with just one contiguous chunk of quads
00084 
00085     double *x_ptr, *y_ptr, *z_ptr;
00086     rval = mbImpl->coords_iterate( verts.begin(), verts.end(), x_ptr, y_ptr, z_ptr, count );MB_CHK_SET_ERR( rval, "Error in coords_iterate" );
00087     assert( count == (int)verts.size() );  // Should end up with just one contiguous chunk of vertices
00088 
00089     double *tag1_ptr, *tag2_ptr;
00090     int* tag3_ptr;
00091     rval = mbImpl->tag_iterate( tag1, quads.begin(), quads.end(), count, (void*&)tag1_ptr );MB_CHK_SET_ERR( rval, "Error in tag1_iterate" );
00092     assert( count == (int)quads.size() );  // Should end up with just one contiguous chunk of quads
00093     rval = mbImpl->tag_iterate( tag2, quads.begin(), quads.end(), count, (void*&)tag2_ptr );MB_CHK_SET_ERR( rval, "Error in tag2_iterate" );
00094     assert( count == (int)quads.size() );  // Should end up with just one contiguous chunk of quads
00095     rval = mbImpl->tag_iterate( tag3, quads.begin(), quads.end(), count, (void*&)tag3_ptr );MB_CHK_SET_ERR( rval, "Error in tag3_iterate" );
00096     assert( count == (int)quads.size() );  // Should end up with just one contiguous chunk of quads
00097 
00098     const vector< EntityHandle >** adjs_ptr;
00099     rval = mbImpl->adjacencies_iterate( verts.begin(), verts.end(), adjs_ptr, count );MB_CHK_SET_ERR( rval, "Error in adjacencies_iterate" );
00100     assert( count == (int)verts.size() );  // Should end up with just one contiguous chunk of vertices
00101     // Start_ handles used to compute indices into vertex/quad arrays; can use direct subtraction
00102     // because we know there aren't any holes in the handle spaces for verts or quads
00103     EntityHandle start_vert = *verts.begin(), start_quad = *quads.begin();
00104 
00105     // Iterate over elements, computing tag1 from coords positions
00106     for( int i = 0; i < nquads; i++ )
00107     {
00108         tag1_ptr[3 * i + 0] = tag1_ptr[3 * i + 1] = tag1_ptr[3 * i + 2] = 0.0;  // Initialize position for this element
00109         for( int j = 0; j < vpere; j++ )
00110         {                                                        // Loop over vertices in this element
00111             int v_index = conn_ptr[vpere * i + j] - start_vert;  // vert index is just the offset from start vertex
00112             tag1_ptr[3 * i + 0] += x_ptr[v_index];
00113             tag1_ptr[3 * i + 1] += y_ptr[v_index];  // Sum vertex positions into tag1...
00114             tag1_ptr[3 * i + 2] += z_ptr[v_index];
00115         }
00116         tag1_ptr[3 * i + 0] /= vpere;
00117         tag1_ptr[3 * i + 1] /= vpere;  // Then normalize
00118         tag1_ptr[3 * i + 2] /= vpere;
00119     }  // Loop over elements in chunk
00120 
00121     // Iterate through vertices, summing positions into tag2 on connected elements and incrementing
00122     // vertex count
00123     for( int v = 0; v < count; v++ )
00124     {
00125         const vector< EntityHandle >* avec = *( adjs_ptr + v );
00126         for( vector< EntityHandle >::const_iterator ait = avec->begin(); ait != avec->end(); ++ait )
00127         {
00128             // *ait is the quad handle, its index is computed by subtracting the start quad handle
00129             int a_ind = *ait - start_quad;
00130             tag2_ptr[3 * a_ind + 0] += x_ptr[v];  // Tag on each element is 3 doubles, x/y/z
00131             tag2_ptr[3 * a_ind + 1] += y_ptr[v];
00132             tag2_ptr[3 * a_ind + 2] += z_ptr[v];
00133             tag3_ptr[a_ind]++;  // Increment the vertex count
00134         }
00135     }
00136 
00137     // Normalize tag2 by vertex count (tag3); loop over elements using same approach as before
00138     // At the same time, compare values of tag1 and tag2
00139     int n_dis = 0;
00140     for( Range::iterator q_it = quads.begin(); q_it != quads.end(); ++q_it )
00141     {
00142         int i = *q_it - start_quad;
00143         for( int j = 0; j < 3; j++ )
00144             tag2_ptr[3 * i + j] /= (double)tag3_ptr[i];  // Normalize by # verts
00145         if( tag1_ptr[3 * i] != tag2_ptr[3 * i] || tag1_ptr[3 * i + 1] != tag2_ptr[3 * i + 1] ||
00146             tag1_ptr[3 * i + 2] != tag2_ptr[3 * i + 2] )
00147         {
00148             cout << "Tag1, tag2 disagree for element " << *q_it + i << endl;
00149             n_dis++;
00150         }
00151     }
00152     if( !n_dis ) cout << "All tags agree, success!" << endl;
00153 
00154     // Ok, we're done, shut down MOAB
00155     delete mbImpl;
00156 
00157     return 0;
00158 }
00159 
00160 ErrorCode create_mesh_no_holes( Interface* mbImpl, int nquads )
00161 {
00162     // First make the mesh, a 1d array of quads with left hand side x = elem_num; vertices are
00163     // numbered in layers
00164     ReadUtilIface* read_iface;
00165     ErrorCode rval = mbImpl->query_interface( read_iface );MB_CHK_SET_ERR( rval, "Error in query_interface" );
00166     vector< double* > coords;
00167     EntityHandle start_vert, start_elem, *connect;
00168     // Create verts, num is 2(nquads+1) because they're in a 1d row; will initialize coords in loop
00169     // over quads later
00170     rval = read_iface->get_node_coords( 3, 2 * ( nquads + 1 ), 0, start_vert, coords );MB_CHK_SET_ERR( rval, "Error in get_node_arrays" );
00171     // Create quads
00172     rval = read_iface->get_element_connect( nquads, 4, MBQUAD, 0, start_elem, connect );MB_CHK_SET_ERR( rval, "Error in get_element_connect" );
00173     for( int i = 0; i < nquads; i++ )
00174     {
00175         coords[0][2 * i] = coords[0][2 * i + 1] = (double)i;  // x values are all i
00176         coords[1][2 * i]                        = 0.0;
00177         coords[1][2 * i + 1]                    = 1.0;          // y coords
00178         coords[2][2 * i] = coords[2][2 * i + 1] = (double)0.0;  // z values, all zero (2d mesh)
00179         EntityHandle quad_v                     = start_vert + 2 * i;
00180         connect[4 * i + 0]                      = quad_v;
00181         connect[4 * i + 1]                      = quad_v + 2;
00182         connect[4 * i + 2]                      = quad_v + 3;
00183         connect[4 * i + 3]                      = quad_v + 1;
00184     }
00185 
00186     // Last two vertices
00187     // Cppcheck warning (false positive): variable coords is assigned a value that is never used
00188     coords[0][2 * nquads] = coords[0][2 * nquads + 1] = (double)nquads;
00189     coords[1][2 * nquads]                             = 0.0;
00190     coords[1][2 * nquads + 1]                         = 1.0;          // y coords
00191     coords[2][2 * nquads] = coords[2][2 * nquads + 1] = (double)0.0;  // z values, all zero (2d mesh)
00192 
00193     // Call a vertex-quad adjacencies function to generate vertex-element adjacencies in MOAB
00194     Range dum_range;
00195     rval = mbImpl->get_adjacencies( &start_vert, 1, 2, false, dum_range );MB_CHK_SET_ERR( rval, "Error in get_adjacencies" );
00196     assert( !dum_range.empty() );
00197 
00198     return MB_SUCCESS;
00199 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines