 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
1 /** @example DirectAccessNoHoles.cpp \n
2 * \brief Use direct access to MOAB data to avoid calling through API \n
3 *
4 * This example creates a 1d row of quad elements, such that all quad and vertex handles
5 * are contiguous in the handle space and in the database. Then it shows how to get access
6 * to pointers to MOAB-native data for vertex coordinates, quad connectivity, tag storage,
7 * and vertex to quad adjacency lists. This allows applications to access this data directly
8 * without going through MOAB's API. In cases where the mesh is not changing (or only mesh
9 * vertices are moving), this can save significant execution time in applications.
10 * \verbatim
11 * ----------------------
12 * | | | |
13 * | | | | ...
14 * | | | |
15 * ----------------------
16 * \endverbatim
17 * -# Initialize MOAB \n
18 * -# Create a quad mesh, as depicted above
19 * -# Create 2 dense tags (tag1, tag2) for avg position to assign to quads, and # verts per quad
20 * (tag3)
21 * -# Get connectivity, coordinate, tag1 iterators
22 * -# Iterate through quads, computing midpoint based on vertex positions, set on quad-based
23 * tag1
24 * -# Iterate through vertices, summing positions into tag2 on connected quads and incrementing
25 * vertex count
26 * -# Iterate through quads, normalizing tag2 by vertex count and comparing values of tag1 and
27 * tag2
28 *
29 * <b>To compile</b>: \n
30 * make DirectAccessNoHoles \n
31 * <b>To run</b>: ./DirectAccessNoHoles [-nquads <# quads>]\n
32 *
33 */
34
35 #include "moab/Core.hpp"
36 #include "moab/ProgOptions.hpp"
37 #include "moab/ReadUtilIface.hpp"
38 #include <map>
39 #include <iostream>
40 #include <cassert>
41
42 using namespace moab;
43 using namespace std;
44
45 ErrorCode create_mesh_no_holes( Interface* mbImpl, int nquads );
46
47 int main( int argc, char** argv )
48 {
49 // Get MOAB instance
50 Interface* mbImpl = new( std::nothrow ) Core;
51 if( NULL == mbImpl ) return 1;
52
53 int nquads = 1000;
54
55 // Parse options
56 ProgOptions opts;
57 opts.addOpt< int >( string( "nquads,n" ), string( "Number of quads in the mesh (default = 1000" ), &nquads );
58 opts.parseCommandLine( argc, argv );
59
60 // Create simple structured mesh with hole, but using unstructured representation
61 ErrorCode rval = create_mesh_no_holes( mbImpl, nquads );MB_CHK_SET_ERR( rval, "Trouble creating mesh" );
62
63 // Get all vertices and non-vertex entities
64 Range verts, quads;
65 rval = mbImpl->get_entities_by_handle( 0, quads );MB_CHK_SET_ERR( rval, "Trouble getting all entities" );
66 verts = quads.subset_by_dimension( 0 );
67 quads -= verts;
68
69 // Create tag1 (element-based avg), tag2 (vertex-based avg), tag3 (# connected verts)
70 Tag tag1, tag2, tag3;
71 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" );
72 double def_val[3] = { 0.0, 0.0, 0.0 }; // need a default value for tag2 because we sum into it
73 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" );
74 int def_val_int = 0; // need a default value for tag3 because we increment it
75 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" );
76
77 // Get pointers to connectivity, coordinate, tag, and adjacency arrays; each of these returns a
78 // count, which should be compared to the # entities you expect to verify there's only one chunk
79 // (no holes)
80 int count, vpere;
81 EntityHandle* conn_ptr;
82 rval = mbImpl->connect_iterate( quads.begin(), quads.end(), conn_ptr, vpere, count );MB_CHK_SET_ERR( rval, "Error in connect_iterate" );
83 assert( count == (int)quads.size() ); // Should end up with just one contiguous chunk of quads
84
85 double *x_ptr, *y_ptr, *z_ptr;
86 rval = mbImpl->coords_iterate( verts.begin(), verts.end(), x_ptr, y_ptr, z_ptr, count );MB_CHK_SET_ERR( rval, "Error in coords_iterate" );
87 assert( count == (int)verts.size() ); // Should end up with just one contiguous chunk of vertices
88
89 double *tag1_ptr, *tag2_ptr;
90 int* tag3_ptr;
91 rval = mbImpl->tag_iterate( tag1, quads.begin(), quads.end(), count, (void*&)tag1_ptr );MB_CHK_SET_ERR( rval, "Error in tag1_iterate" );
92 assert( count == (int)quads.size() ); // Should end up with just one contiguous chunk of quads
93 rval = mbImpl->tag_iterate( tag2, quads.begin(), quads.end(), count, (void*&)tag2_ptr );MB_CHK_SET_ERR( rval, "Error in tag2_iterate" );
94 assert( count == (int)quads.size() ); // Should end up with just one contiguous chunk of quads
95 rval = mbImpl->tag_iterate( tag3, quads.begin(), quads.end(), count, (void*&)tag3_ptr );MB_CHK_SET_ERR( rval, "Error in tag3_iterate" );
96 assert( count == (int)quads.size() ); // Should end up with just one contiguous chunk of quads
97
98 const vector< EntityHandle >** adjs_ptr;
99 rval = mbImpl->adjacencies_iterate( verts.begin(), verts.end(), adjs_ptr, count );MB_CHK_SET_ERR( rval, "Error in adjacencies_iterate" );
100 assert( count == (int)verts.size() ); // Should end up with just one contiguous chunk of vertices
101 // Start_ handles used to compute indices into vertex/quad arrays; can use direct subtraction
102 // because we know there aren't any holes in the handle spaces for verts or quads
103 EntityHandle start_vert = *verts.begin(), start_quad = *quads.begin();
104
105 // Iterate over elements, computing tag1 from coords positions
106 for( int i = 0; i < nquads; i++ )
107 {
108 tag1_ptr[3 * i + 0] = tag1_ptr[3 * i + 1] = tag1_ptr[3 * i + 2] = 0.0; // Initialize position for this element
109 for( int j = 0; j < vpere; j++ )
110 { // Loop over vertices in this element
111 int v_index = conn_ptr[vpere * i + j] - start_vert; // vert index is just the offset from start vertex
112 tag1_ptr[3 * i + 0] += x_ptr[v_index];
113 tag1_ptr[3 * i + 1] += y_ptr[v_index]; // Sum vertex positions into tag1...
114 tag1_ptr[3 * i + 2] += z_ptr[v_index];
115 }
116 tag1_ptr[3 * i + 0] /= vpere;
117 tag1_ptr[3 * i + 1] /= vpere; // Then normalize
118 tag1_ptr[3 * i + 2] /= vpere;
119 } // Loop over elements in chunk
120
121 // Iterate through vertices, summing positions into tag2 on connected elements and incrementing
122 // vertex count
123 for( int v = 0; v < count; v++ )
124 {
125 const vector< EntityHandle >* avec = *( adjs_ptr + v );
126 for( vector< EntityHandle >::const_iterator ait = avec->begin(); ait != avec->end(); ++ait )
127 {
128 // *ait is the quad handle, its index is computed by subtracting the start quad handle
129 int a_ind = *ait - start_quad;
130 tag2_ptr[3 * a_ind + 0] += x_ptr[v]; // Tag on each element is 3 doubles, x/y/z
131 tag2_ptr[3 * a_ind + 1] += y_ptr[v];
132 tag2_ptr[3 * a_ind + 2] += z_ptr[v];
133 tag3_ptr[a_ind]++; // Increment the vertex count
134 }
135 }
136
137 // Normalize tag2 by vertex count (tag3); loop over elements using same approach as before
138 // At the same time, compare values of tag1 and tag2
139 int n_dis = 0;
140 for( Range::iterator q_it = quads.begin(); q_it != quads.end(); ++q_it )
141 {
142 int i = *q_it - start_quad;
143 for( int j = 0; j < 3; j++ )
144 tag2_ptr[3 * i + j] /= (double)tag3_ptr[i]; // Normalize by # verts
145 if( tag1_ptr[3 * i] != tag2_ptr[3 * i] || tag1_ptr[3 * i + 1] != tag2_ptr[3 * i + 1] ||
146 tag1_ptr[3 * i + 2] != tag2_ptr[3 * i + 2] )
147 {
148 cout << "Tag1, tag2 disagree for element " << *q_it + i << endl;
149 n_dis++;
150 }
151 }
152 if( !n_dis ) cout << "All tags agree, success!" << endl;
153
154 // Ok, we're done, shut down MOAB
155 delete mbImpl;
156
157 return 0;
158 }
159
160 ErrorCode create_mesh_no_holes( Interface* mbImpl, int nquads )
161 {
162 // First make the mesh, a 1d array of quads with left hand side x = elem_num; vertices are
163 // numbered in layers
164 ReadUtilIface* read_iface;
165 ErrorCode rval = mbImpl->query_interface( read_iface );MB_CHK_SET_ERR( rval, "Error in query_interface" );
166 vector< double* > coords;
167 EntityHandle start_vert, start_elem, *connect;
168 // Create verts, num is 2(nquads+1) because they're in a 1d row; will initialize coords in loop
169 // over quads later
170 rval = read_iface->get_node_coords( 3, 2 * ( nquads + 1 ), 0, start_vert, coords );MB_CHK_SET_ERR( rval, "Error in get_node_arrays" );
171 // Create quads
172 rval = read_iface->get_element_connect( nquads, 4, MBQUAD, 0, start_elem, connect );MB_CHK_SET_ERR( rval, "Error in get_element_connect" );
173 for( int i = 0; i < nquads; i++ )
174 {
175 coords[0][2 * i] = coords[0][2 * i + 1] = (double)i; // x values are all i
176 coords[1][2 * i] = 0.0;
177 coords[1][2 * i + 1] = 1.0; // y coords
178 coords[2][2 * i] = coords[2][2 * i + 1] = (double)0.0; // z values, all zero (2d mesh)
179 EntityHandle quad_v = start_vert + 2 * i;
180 connect[4 * i + 0] = quad_v;
181 connect[4 * i + 1] = quad_v + 2;
182 connect[4 * i + 2] = quad_v + 3;
183 connect[4 * i + 3] = quad_v + 1;
184 }
185
186 // Last two vertices
187 // Cppcheck warning (false positive): variable coords is assigned a value that is never used
188 coords[0][2 * nquads] = coords[0][2 * nquads + 1] = (double)nquads;
189 coords[1][2 * nquads] = 0.0;
190 coords[1][2 * nquads + 1] = 1.0; // y coords
191 coords[2][2 * nquads] = coords[2][2 * nquads + 1] = (double)0.0; // z values, all zero (2d mesh)
192
193 // Call a vertex-quad adjacencies function to generate vertex-element adjacencies in MOAB
194 Range dum_range;
195 rval = mbImpl->get_adjacencies( &start_vert, 1, 2, false, dum_range );MB_CHK_SET_ERR( rval, "Error in get_adjacencies" );
196 assert( !dum_range.empty() );
197
198 return MB_SUCCESS;
199 }
1.9.1.