MOAB: Mesh Oriented datABase  (version 5.5.0)
adj_moab_test.cpp
Go to the documentation of this file.
1 /*This function tests the AHF datastructures on CST meshes*/
2 #include <iostream>
3 #include <vector>
4 #include <algorithm>
5 #include "moab/Core.hpp"
6 #include "moab/Range.hpp"
7 #include "moab/MeshTopoUtil.hpp"
8 #include "moab/HalfFacetRep.hpp"
9 #include "TestUtil.hpp"
10 
11 #ifdef MOAB_HAVE_MPI
12 #include "moab/ParallelComm.hpp"
13 #include "MBParallelConventions.h"
14 #include "moab/FileOptions.hpp"
15 #include "MBTagConventions.hpp"
16 #include "moab_mpi.h"
17 #endif
18 
19 using namespace moab;
20 
21 #ifdef MOAB_HAVE_MPI
22 std::string read_options;
23 #endif
24 
25 int number_tests_successful = 0;
26 int number_tests_failed = 0;
27 
28 void handle_error_code( ErrorCode rv, int& number_failed, int& number_successful )
29 {
30  if( rv == MB_SUCCESS )
31  {
32 #ifdef MOAB_HAVE_MPI
33  int rank = 0;
34  MPI_Comm_rank( MPI_COMM_WORLD, &rank );
35  if( rank == 0 ) std::cout << "Success";
36 #else
37  std::cout << "Success";
38 #endif
39  number_successful++;
40  }
41  else
42  {
43  std::cout << "Failure";
44  number_failed++;
45  }
46 }
47 
48 ErrorCode ahf_test( const char* filename )
49 {
50 
51  Core moab;
52  Interface* mbImpl = &moab;
53  ParallelComm* pc = NULL;
54  MeshTopoUtil mtu( mbImpl );
55  ErrorCode error;
56  EntityHandle fileset;
57  error = mbImpl->create_meshset( moab::MESHSET_SET, fileset );CHECK_ERR( error );
58 
59 #ifdef MOAB_HAVE_MPI
60  int procs = 1;
61  MPI_Comm_size( MPI_COMM_WORLD, &procs );
62 
63  if( procs > 1 )
64  {
65  read_options = "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARALLEL_RESOLVE_SHARED_ENTS;";
66 
67  error = mbImpl->load_file( filename, &fileset, read_options.c_str() );CHECK_ERR( error );
68  }
69  else if( procs == 1 )
70  {
71 #endif
72  error = mbImpl->load_file( filename, &fileset );CHECK_ERR( error );
73 #ifdef MOAB_HAVE_MPI
74  }
75 #endif
76 
77  /*Create ranges for handles of explicit elements of the mixed mesh*/
78  Range verts, edges, faces, cells;
79  error = mbImpl->get_entities_by_dimension( fileset, 0, verts );CHECK_ERR( error );
80  error = mbImpl->get_entities_by_dimension( fileset, 1, edges );CHECK_ERR( error );
81  error = mbImpl->get_entities_by_dimension( fileset, 2, faces );CHECK_ERR( error );
82  error = mbImpl->get_entities_by_dimension( fileset, 3, cells );CHECK_ERR( error );
83 
84  // Create an ahf instance
85 #ifdef MOAB_HAVE_MPI
86  pc = ParallelComm::get_pcomm( mbImpl, 0 );
87  if( !pc ) pc = new moab::ParallelComm( &moab, MPI_COMM_WORLD );
88 #endif
89 
90  HalfFacetRep ahf( &moab, pc, fileset );
91 
92  // Call the initialize function which creates the maps for each dimension
93  error = ahf.initialize();CHECK_ERR( error );
94 
95  std::cout << "Finished AHF initialization" << std::endl;
96 
97  // Perform queries
98  std::vector< EntityHandle > adjents;
99  Range mbents, ahfents;
100 
101  // 1D Queries //
102  // IQ1: For every vertex, obtain incident edges
103  if( edges.size() )
104  {
105  Range everts;
106  error = mbImpl->get_connectivity( edges, everts, true );CHECK_ERR( error );
107  for( Range::iterator i = everts.begin(); i != everts.end(); ++i )
108  {
109  adjents.clear();
110  error = ahf.get_up_adjacencies( *i, 1, adjents );CHECK_ERR( error );
111  mbents.clear();
112  error = mbImpl->get_adjacencies( &*i, 1, 1, false, mbents );CHECK_ERR( error );
113 
114  CHECK_EQUAL( mbents.size(), adjents.size() );
115  std::sort( adjents.begin(), adjents.end() );
116  std::copy( adjents.begin(), adjents.end(), range_inserter( ahfents ) );
117  mbents = subtract( mbents, ahfents );
118  CHECK( !mbents.size() );
119  }
120  }
121 
122  // NQ1: For every edge, obtain neighbor edges
123  if( edges.size() )
124  {
125  for( Range::iterator i = edges.begin(); i != edges.end(); ++i )
126  {
127  adjents.clear();
128  error = ahf.get_neighbor_adjacencies( *i, adjents );CHECK_ERR( error );
129  mbents.clear();
130  error = mtu.get_bridge_adjacencies( *i, 0, 1, mbents );CHECK_ERR( error );
131 
132  CHECK_EQUAL( mbents.size(), adjents.size() );
133  std::sort( adjents.begin(), adjents.end() );
134  std::copy( adjents.begin(), adjents.end(), range_inserter( ahfents ) );
135  mbents = subtract( mbents, ahfents );
136  CHECK( !mbents.size() );
137  }
138  }
139 
140  std::cout << "Finished 1D queries" << std::endl;
141 
142  // 2D Queries
143  // IQ21: For every vertex, obtain incident faces
144  if( faces.size() )
145  {
146  Range fverts;
147  error = mbImpl->get_connectivity( faces, fverts, true );CHECK_ERR( error );
148  for( Range::iterator i = fverts.begin(); i != fverts.end(); ++i )
149  {
150  adjents.clear();
151  error = ahf.get_up_adjacencies( *i, 2, adjents );CHECK_ERR( error );
152  mbents.clear();
153  error = mbImpl->get_adjacencies( &*i, 1, 2, false, mbents );CHECK_ERR( error );
154 
155  CHECK_EQUAL( mbents.size(), adjents.size() );
156  std::sort( adjents.begin(), adjents.end() );
157  std::copy( adjents.begin(), adjents.end(), range_inserter( ahfents ) );
158  mbents = subtract( mbents, ahfents );
159  CHECK( !mbents.size() );
160  }
161  }
162 
163  // IQ22: For every edge, obtain incident faces
164  if( edges.size() && faces.size() )
165  {
166  for( Range::iterator i = edges.begin(); i != edges.end(); ++i )
167  {
168  adjents.clear();
169  error = ahf.get_up_adjacencies( *i, 2, adjents );CHECK_ERR( error );
170  mbents.clear();
171  error = mbImpl->get_adjacencies( &*i, 1, 2, false, mbents );CHECK_ERR( error );
172 
173  CHECK_EQUAL( mbents.size(), adjents.size() );
174  std::sort( adjents.begin(), adjents.end() );
175  std::copy( adjents.begin(), adjents.end(), range_inserter( ahfents ) );
176  mbents = subtract( mbents, ahfents );
177  CHECK( !mbents.size() );
178  }
179  }
180 
181  // NQ2: For every face, obtain neighbor faces
182  if( faces.size() )
183  {
184  for( Range::iterator i = faces.begin(); i != faces.end(); ++i )
185  {
186  adjents.clear();
187  error = ahf.get_neighbor_adjacencies( *i, adjents );CHECK_ERR( error );
188  mbents.clear();
189  error = mtu.get_bridge_adjacencies( *i, 1, 2, mbents );CHECK_ERR( error );
190 
191  CHECK_EQUAL( mbents.size(), adjents.size() );
192  std::sort( adjents.begin(), adjents.end() );
193  std::copy( adjents.begin(), adjents.end(), range_inserter( ahfents ) );
194  mbents = subtract( mbents, ahfents );
195  CHECK( !mbents.size() );
196  }
197  }
198 
199  // DQ 21: For every face, obtain its edges
200  if( edges.size() && faces.size() )
201  {
202  for( Range::iterator i = faces.begin(); i != faces.end(); ++i )
203  {
204  adjents.clear();
205  error = ahf.get_down_adjacencies( *i, 1, adjents );CHECK_ERR( error );
206  mbents.clear();
207  error = mbImpl->get_adjacencies( &*i, 1, 1, false, mbents );CHECK_ERR( error );
208 
209  CHECK_EQUAL( mbents.size(), adjents.size() );
210  std::sort( adjents.begin(), adjents.end() );
211  std::copy( adjents.begin(), adjents.end(), range_inserter( ahfents ) );
212  mbents = subtract( mbents, ahfents );
213  CHECK( !mbents.size() );
214  }
215  }
216 
217  std::cout << "Finished 2D queries: " << std::endl;
218 
219  // 3D Queries
220  // IQ 31: For every vertex, obtain incident cells
221  if( cells.size() )
222  {
223  Range cverts;
224  error = mbImpl->get_connectivity( cells, cverts, true );CHECK_ERR( error );
225  for( Range::iterator i = cverts.begin(); i != cverts.end(); ++i )
226  {
227  adjents.clear();
228  error = ahf.get_up_adjacencies( *i, 3, adjents );CHECK_ERR( error );
229  mbents.clear();
230  error = mbImpl->get_adjacencies( &*i, 1, 3, false, mbents );CHECK_ERR( error );
231 
232  CHECK_EQUAL( mbents.size(), adjents.size() );
233  std::sort( adjents.begin(), adjents.end() );
234  std::copy( adjents.begin(), adjents.end(), range_inserter( ahfents ) );
235  mbents = subtract( mbents, ahfents );
236  CHECK( !mbents.size() );
237  }
238  }
239 
240  // IQ 32: For every edge, obtain incident cells
241  if( edges.size() && cells.size() )
242  {
243  for( Range::iterator i = edges.begin(); i != edges.end(); ++i )
244  {
245  adjents.clear();
246  error = ahf.get_up_adjacencies( *i, 3, adjents );CHECK_ERR( error );
247  mbents.clear();
248  error = mbImpl->get_adjacencies( &*i, 1, 3, false, mbents );CHECK_ERR( error );
249 
250  if( adjents.size() != mbents.size() )
251  {
252  // std::cout<<"ahf results = "<<std::endl;
253  // ahfents.print();
254  // std::cout<<"native results = "<<std::endl;
255  // mbents.print();
256  }
257 
258  CHECK_EQUAL( mbents.size(), adjents.size() );
259  std::sort( adjents.begin(), adjents.end() );
260  std::copy( adjents.begin(), adjents.end(), range_inserter( ahfents ) );
261  mbents = subtract( mbents, ahfents );
262  CHECK( !mbents.size() );
263  }
264  }
265 
266  // IQ33: For every face, obtain incident cells
267  if( faces.size() && cells.size() )
268  {
269  for( Range::iterator i = faces.begin(); i != faces.end(); ++i )
270  {
271  adjents.clear();
272  error = ahf.get_up_adjacencies( *i, 3, adjents );CHECK_ERR( error );
273  mbents.clear();
274  error = mbImpl->get_adjacencies( &*i, 1, 3, false, mbents );CHECK_ERR( error );
275 
276  if( adjents.size() != mbents.size() )
277  {
278  // std::cout<<"ahf results = "<<std::endl;
279  // ahfents.print();
280  // std::cout<<"native results = "<<std::endl;
281  // mbents.print();
282  }
283 
284  CHECK_EQUAL( mbents.size(), adjents.size() );
285  std::sort( adjents.begin(), adjents.end() );
286  std::copy( adjents.begin(), adjents.end(), range_inserter( ahfents ) );
287  mbents = subtract( mbents, ahfents );
288  CHECK( !mbents.size() );
289  }
290  }
291 
292  // NQ3: For every cell, obtain neighbor cells
293  if( cells.size() )
294  {
295  for( Range::iterator i = cells.begin(); i != cells.end(); ++i )
296  {
297  adjents.clear();
298  error = ahf.get_neighbor_adjacencies( *i, adjents );CHECK_ERR( error );
299  mbents.clear();
300  error = mtu.get_bridge_adjacencies( *i, 2, 3, mbents );CHECK_ERR( error );
301 
302  CHECK_EQUAL( mbents.size(), adjents.size() );
303  std::sort( adjents.begin(), adjents.end() );
304  std::copy( adjents.begin(), adjents.end(), range_inserter( ahfents ) );
305  mbents = subtract( mbents, ahfents );
306  CHECK( !mbents.size() );
307  }
308  }
309 
310  // DQ 31: For every cell, obtain its edges
311  if( edges.size() && cells.size() )
312  {
313  for( Range::iterator i = cells.begin(); i != cells.end(); ++i )
314  {
315  adjents.clear();
316  error = ahf.get_down_adjacencies( *i, 1, adjents );CHECK_ERR( error );
317  mbents.clear();
318  error = mbImpl->get_adjacencies( &*i, 1, 1, false, mbents );CHECK_ERR( error );
319 
320  CHECK_EQUAL( mbents.size(), adjents.size() );
321  std::sort( adjents.begin(), adjents.end() );
322  std::copy( adjents.begin(), adjents.end(), range_inserter( ahfents ) );
323  mbents = subtract( mbents, ahfents );
324  CHECK( !mbents.size() );
325  }
326  }
327 
328  // DQ 32: For every cell, obtain its faces
329  if( faces.size() && cells.size() )
330  {
331  for( Range::iterator i = cells.begin(); i != cells.end(); ++i )
332  {
333  adjents.clear();
334  error = ahf.get_down_adjacencies( *i, 2, adjents );CHECK_ERR( error );
335  mbents.clear();
336  error = mbImpl->get_adjacencies( &*i, 1, 2, false, mbents );CHECK_ERR( error );
337 
338  CHECK_EQUAL( mbents.size(), adjents.size() );
339  std::sort( adjents.begin(), adjents.end() );
340  std::copy( adjents.begin(), adjents.end(), range_inserter( ahfents ) );
341  mbents = subtract( mbents, ahfents );
342  CHECK( !mbents.size() );
343  }
344  }
345 
346  std::cout << "Finished 3D queries" << std::endl;
347 
348  error = ahf.deinitialize();CHECK_ERR( error );
349 
350  return MB_SUCCESS;
351 }
352 
353 int main( int argc, char* argv[] )
354 {
355 
356 #ifdef MOAB_HAVE_MPI
357  MPI_Init( &argc, &argv );
358 
359  int nprocs, rank;
360  MPI_Comm_size( MPI_COMM_WORLD, &nprocs );
361  MPI_Comm_rank( MPI_COMM_WORLD, &rank );
362 #endif
363 
364  std::string filename;
365 #ifdef MOAB_HAVE_HDF5
366 #ifdef MOAB_HAVE_AHF
367  filename = TestDir + "unittest/spectral.h5m";
368 #else
369  filename = TestDir + "unittest/32hex_ef.h5m";
370 #endif
371 #else
372  filename = TestDir + "unittest/hexes_mixed.vtk";
373 #endif
374 
375  if( argc == 1 )
376  {
377 #ifdef MOAB_HAVE_MPI
378  if( rank == 0 ) std::cout << "Using default input file:" << filename << std::endl;
379 #else
380  std::cout << "Using default input file:" << filename << std::endl;
381 #endif
382  }
383 
384  else if( argc == 2 )
385  filename = std::string( argv[1] );
386  else
387  {
388  std::cerr << "Usage: " << argv[0] << " [filename]" << std::endl;
389  return 1;
390  }
391 
392  ErrorCode result;
393 
394 #ifdef MOAB_HAVE_MPI
395  if( rank == 0 ) std::cout << " para_ahf_test: ";
396 #else
397  std::cout << "ahf_test:";
398 #endif
399 
400  result = ahf_test( filename.c_str() );
401  handle_error_code( result, number_tests_failed, number_tests_successful );
402  std::cout << "\n";
403 
404 #ifdef MOAB_HAVE_MPI
405  MPI_Finalize();
406 #endif
407 
408  return number_tests_failed;
409 }