Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
ReadNCDF.cpp
Go to the documentation of this file.
1 /**
2  * MOAB, a Mesh-Oriented datABase, is a software component for creating,
3  * storing and accessing finite element mesh data.
4  *
5  * Copyright 2004 Sandia Corporation. Under the terms of Contract
6  * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
7  * retains certain rights in this software.
8  *
9  * This library is free software; you can redistribute it and/or
10  * modify it under the terms of the GNU Lesser General Public
11  * License as published by the Free Software Foundation; either
12  * version 2.1 of the License, or (at your option) any later version.
13  *
14  */
15 
16 #ifdef WIN32
17 #pragma warning( disable : 4786 )
18 #endif
19 
20 #include "ReadNCDF.hpp"
21 #include "netcdf.h"
22 
23 #include <algorithm>
24 #include <string>
25 #include <cassert>
26 #include <cstdio>
27 #include <cstring>
28 #include <cmath>
29 #include <sstream>
30 #include <iostream>
31 #include <map>
32 
33 #include "moab/CN.hpp"
34 #include "moab/Range.hpp"
35 #include "moab/Interface.hpp"
36 #include "ExoIIUtil.hpp"
37 #include "MBTagConventions.hpp"
38 #include "Internals.hpp"
39 #include "moab/ReadUtilIface.hpp"
40 #include "exodus_order.h"
41 #include "moab/FileOptions.hpp"
42 #include "moab/AdaptiveKDTree.hpp"
43 #include "moab/CartVect.hpp"
44 
45 namespace moab
46 {
47 
48 #define INS_ID( stringvar, prefix, id, length ) snprintf( stringvar, length, prefix, id )
49 
50 #define GET_DIM( ncdim, name, val ) \
51  { \
52  int gdfail = nc_inq_dimid( ncFile, name, &( ncdim ) ); \
53  if( NC_NOERR == gdfail ) \
54  { \
55  size_t tmp_val; \
56  gdfail = nc_inq_dimlen( ncFile, ncdim, &tmp_val ); \
57  if( NC_NOERR != gdfail ) \
58  { \
59  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Couldn't get dimension length" ); \
60  } \
61  else \
62  ( val ) = tmp_val; \
63  } \
64  else \
65  ( val ) = 0; \
66  }
67 
68 #define GET_DIMB( ncdim, name, varname, id, val ) \
69  INS_ID( name, varname, id, max_str_length + 1 ); \
70  GET_DIM( ncdim, name, val );
71 
72 #define GET_VAR( name, id, dims ) \
73  { \
74  ( id ) = -1; \
75  int gvfail = nc_inq_varid( ncFile, name, &( id ) ); \
76  if( NC_NOERR == gvfail ) \
77  { \
78  int ndims; \
79  gvfail = nc_inq_varndims( ncFile, id, &ndims ); \
80  if( NC_NOERR == gvfail ) \
81  { \
82  ( dims ).resize( ndims ); \
83  gvfail = nc_inq_vardimid( ncFile, id, &( dims )[0] ); \
84  if( NC_NOERR != gvfail ) \
85  { \
86  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Couldn't get variable dimension IDs" ); \
87  } \
88  } \
89  } \
90  }
91 
92 #define GET_1D_INT_VAR( name, id, vals ) \
93  { \
94  GET_VAR( name, id, vals ); \
95  if( -1 != ( id ) ) \
96  { \
97  size_t ntmp; \
98  int ivfail = nc_inq_dimlen( ncFile, ( vals )[0], &ntmp ); \
99  if( NC_NOERR != ivfail ) \
100  { \
101  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Couldn't get dimension length" ); \
102  } \
103  ( vals ).resize( ntmp ); \
104  size_t ntmp1 = 0; \
105  ivfail = nc_get_vara_int( ncFile, id, &ntmp1, &ntmp, &( vals )[0] ); \
106  if( NC_NOERR != ivfail ) \
107  { \
108  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting variable " << ( name ) ); \
109  } \
110  } \
111  }
112 
113 #define GET_1D_DBL_VAR( name, id, vals ) \
114  { \
115  std::vector< int > dum_dims; \
116  GET_VAR( name, id, dum_dims ); \
117  if( -1 != ( id ) ) \
118  { \
119  size_t ntmp; \
120  int dvfail = nc_inq_dimlen( ncFile, dum_dims[0], &ntmp ); \
121  if( NC_NOERR != dvfail ) \
122  { \
123  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Couldn't get dimension length" ); \
124  } \
125  ( vals ).resize( ntmp ); \
126  size_t ntmp1 = 0; \
127  dvfail = nc_get_vara_double( ncFile, id, &ntmp1, &ntmp, &( vals )[0] ); \
128  if( NC_NOERR != dvfail ) \
129  { \
130  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting variable " << ( name ) ); \
131  } \
132  } \
133  }
134 
135 ReaderIface* ReadNCDF::factory( Interface* iface )
136 {
137  return new ReadNCDF( iface );
138 }
139 
140 ReadNCDF::ReadNCDF( Interface* impl ) : mdbImpl( impl ), max_line_length( -1 ), max_str_length( -1 )
141 {
142  assert( impl != NULL );
143  reset();
144 
145  impl->query_interface( readMeshIface );
146 
147  // Initialize in case tag_get_handle fails below
148  mMaterialSetTag = 0;
149  mDirichletSetTag = 0;
150  mNeumannSetTag = 0;
151  mHasMidNodesTag = 0;
152  mDistFactorTag = 0;
153  mQaRecordTag = 0;
154  mGlobalIdTag = 0;
155 
156  //! Get and cache predefined tag handles
157  ErrorCode result;
158  const int negone = -1;
159  result = impl->tag_get_handle( MATERIAL_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mMaterialSetTag,
160  MB_TAG_SPARSE | MB_TAG_CREAT, &negone );
161  assert( MB_SUCCESS == result );
162  result = impl->tag_get_handle( DIRICHLET_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mDirichletSetTag,
163  MB_TAG_SPARSE | MB_TAG_CREAT, &negone );
164  assert( MB_SUCCESS == result );
165  result = impl->tag_get_handle( NEUMANN_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mNeumannSetTag,
166  MB_TAG_SPARSE | MB_TAG_CREAT, &negone );
167  assert( MB_SUCCESS == result );
168  const int mids[] = { -1, -1, -1, -1 };
169  result = impl->tag_get_handle( HAS_MID_NODES_TAG_NAME, 4, MB_TYPE_INTEGER, mHasMidNodesTag,
170  MB_TAG_SPARSE | MB_TAG_CREAT, mids );
171  assert( MB_SUCCESS == result );
172  result = impl->tag_get_handle( "distFactor", 0, MB_TYPE_DOUBLE, mDistFactorTag,
173  MB_TAG_SPARSE | MB_TAG_VARLEN | MB_TAG_CREAT );
174  assert( MB_SUCCESS == result );
175  result = impl->tag_get_handle( "qaRecord", 0, MB_TYPE_OPAQUE, mQaRecordTag,
176  MB_TAG_SPARSE | MB_TAG_VARLEN | MB_TAG_CREAT );
177  assert( MB_SUCCESS == result );
178  mGlobalIdTag = impl->globalId_tag();
179 
180  assert( MB_SUCCESS == result );
181 #ifdef NDEBUG
182  if( MB_SUCCESS == result )
183  {
184  }; // Line to avoid compiler warning about unused variable
185 #endif
186  ncFile = 0;
187 }
188 
189 void ReadNCDF::reset()
190 {
191  mCurrentMeshHandle = 0;
192  vertexOffset = 0;
193 
194  numberNodes_loading = 0;
195  numberElements_loading = 0;
196  numberElementBlocks_loading = 0;
197  numberFaceBlocks_loading = 0;
198  numberNodeSets_loading = 0;
199  numberSideSets_loading = 0;
200  numberDimensions_loading = 0;
201 
202  if( !blocksLoading.empty() ) blocksLoading.clear();
203 
204  if( !nodesInLoadedBlocks.empty() ) nodesInLoadedBlocks.clear();
205 
206  if( !faceBlocksLoading.empty() ) faceBlocksLoading.clear();
207 }
208 
209 ReadNCDF::~ReadNCDF()
210 {
211  mdbImpl->release_interface( readMeshIface );
212 }
213 
214 ErrorCode ReadNCDF::read_tag_values( const char* file_name,
215  const char* tag_name,
216  const FileOptions& /* opts */,
217  std::vector< int >& id_array,
218  const SubsetList* subset_list )
219 {
220  if( subset_list )
221  {
222  MB_SET_ERR( MB_UNSUPPORTED_OPERATION, "ExodusII reader supports subset read only by material ID" );
223  }
224 
225  // Open netcdf/exodus file
226  int fail = nc_open( file_name, 0, &ncFile );
227  if( NC_NOWRITE != fail )
228  {
229  MB_SET_ERR( MB_FILE_DOES_NOT_EXIST, "ReadNCDF:: problem opening Netcdf/Exodus II file " << file_name );
230  }
231 
232  // 1. Read the header
233  ErrorCode rval = read_exodus_header();
234  if( MB_FAILURE == rval ) return rval;
235 
236  int count = 0;
237  const char* prop;
238  const char* blocks = "eb_prop1";
239  const char* nodesets = "ns_prop1";
240  const char* sidesets = "ss_prop1";
241 
242  if( !strcmp( tag_name, MATERIAL_SET_TAG_NAME ) )
243  {
244  count = numberElementBlocks_loading;
245  prop = blocks;
246  }
247  else if( !strcmp( tag_name, DIRICHLET_SET_TAG_NAME ) )
248  {
249  count = numberNodeSets_loading;
250  prop = nodesets;
251  }
252  else if( !strcmp( tag_name, NEUMANN_SET_TAG_NAME ) )
253  {
254  count = numberSideSets_loading;
255  prop = sidesets;
256  }
257  else
258  {
259  ncFile = 0;
260  return MB_TAG_NOT_FOUND;
261  }
262 
263  if( count )
264  {
265  int nc_var = -1;
266  GET_1D_INT_VAR( prop, nc_var, id_array );
267  if( !nc_var )
268  {
269  MB_SET_ERR( MB_FAILURE, "Problem getting prop variable" );
270  }
271  }
272 
273  // Close the file
274  fail = nc_close( ncFile );
275  if( NC_NOERR != fail )
276  {
277  MB_SET_ERR( MB_FAILURE, "Trouble closing file" );
278  }
279 
280  ncFile = 0;
281  return MB_SUCCESS;
282 }
283 
284 ErrorCode ReadNCDF::load_file( const char* exodus_file_name,
285  const EntityHandle* file_set,
286  const FileOptions& opts,
287  const ReaderIface::SubsetList* subset_list,
288  const Tag* file_id_tag )
289 {
290  ErrorCode status;
291  int fail;
292 
293  int num_blocks = 0;
294  const int* blocks_to_load = 0;
295  if( subset_list )
296  {
297  if( subset_list->tag_list_length > 1 || !strcmp( subset_list->tag_list[0].tag_name, MATERIAL_SET_TAG_NAME ) )
298  {
299  MB_SET_ERR( MB_UNSUPPORTED_OPERATION, "ExodusII reader supports subset read only by material ID" );
300  }
301  if( subset_list->num_parts )
302  {
303  MB_SET_ERR( MB_UNSUPPORTED_OPERATION, "ExodusII reader does not support mesh partitioning" );
304  }
305  blocks_to_load = subset_list->tag_list[0].tag_values;
306  num_blocks = subset_list->tag_list[0].num_tag_values;
307  }
308 
309  // This function directs the reading of an exoii file, but doesn't do any of
310  // the actual work
311 
312  // See if opts has tdata.
313  ErrorCode rval;
314  std::string s;
315  rval = opts.get_str_option( "tdata", s );
316  if( MB_SUCCESS == rval && !s.empty() )
317  return update( exodus_file_name, opts, num_blocks, blocks_to_load, *file_set );
318 
319  reset();
320 
321  // 0. Open the file.
322 
323  // open netcdf/exodus file
324  fail = nc_open( exodus_file_name, 0, &ncFile );
325  if( NC_NOERR != fail )
326  {
327  MB_SET_ERR( MB_FILE_DOES_NOT_EXIST, "ReadNCDF:: problem opening Netcdf/Exodus II file " << exodus_file_name );
328  }
329 
330  // 1. Read the header
331  status = read_exodus_header();
332  if( MB_FAILURE == status ) return status;
333 
334  status = mdbImpl->get_entities_by_handle( 0, initRange );
335  if( MB_FAILURE == status ) return status;
336 
337  // 2. Read the nodes unless they've already been read before
338  status = read_nodes( file_id_tag );
339  if( MB_FAILURE == status ) return status;
340 
341  // 3.
342  // extra for polyhedra blocks
343  if( numberFaceBlocks_loading > 0 )
344  {
345  status = read_face_blocks_headers();
346  if( MB_FAILURE == status ) return status;
347  }
348 
349  status = read_block_headers( blocks_to_load, num_blocks );
350  if( MB_FAILURE == status ) return status;
351 
352  // 4. Read elements (might not read them, depending on active blocks)
353  if( numberFaceBlocks_loading > 0 )
354  {
355  status = read_polyhedra_faces();
356  if( MB_FAILURE == status ) return status;
357  }
358  status = read_elements( file_id_tag );
359  if( MB_FAILURE == status ) return status;
360 
361  // 5. Read global ids
362  status = read_global_ids();
363  if( MB_FAILURE == status ) return status;
364 
365  // 6. Read nodesets
366  status = read_nodesets();
367  if( MB_FAILURE == status ) return status;
368 
369  // 7. Read sidesets
370  status = read_sidesets();
371  if( MB_FAILURE == status ) return status;
372 
373  // 8. Read qa records
374  if( file_set )
375  {
376  status = read_qa_records( *file_set );
377  if( MB_FAILURE == status ) return status;
378  }
379 
380  // What about properties???
381 
382  // Close the file
383  fail = nc_close( ncFile );
384  if( NC_NOERR != fail )
385  {
386  MB_SET_ERR( MB_FAILURE, "Trouble closing file" );
387  }
388 
389  ncFile = 0;
390  return MB_SUCCESS;
391 }
392 
393 ErrorCode ReadNCDF::read_exodus_header()
394 {
395  CPU_WORD_SIZE = sizeof( double ); // With ExodusII version 2, all floats
396  IO_WORD_SIZE = sizeof( double ); // should be changed to doubles
397 
398  // NetCDF doesn't check its own limits on file read, so check
399  // them here so it doesn't corrupt memory any more than absolutely
400  // necessary.
401  int ndims;
402  int fail = nc_inq_ndims( ncFile, &ndims );
403  if( NC_NOERR != fail || ndims > NC_MAX_DIMS )
404  {
405  MB_SET_ERR( MB_FAILURE, "ReadNCDF: File contains " << ndims << " dims but NetCDF library supports only "
406  << (int)NC_MAX_DIMS );
407  }
408  int nvars;
409  fail = nc_inq_nvars( ncFile, &nvars );
410  if( nvars > NC_MAX_VARS )
411  {
412  MB_SET_ERR( MB_FAILURE, "ReadNCDF: File contains " << nvars << " vars but NetCDF library supports only "
413  << (int)NC_MAX_VARS );
414  }
415 
416  // Get the attributes
417 
418  // Get the word size, scalar value
419  nc_type att_type;
420  size_t att_len;
421  fail = nc_inq_att( ncFile, NC_GLOBAL, "floating_point_word_size", &att_type, &att_len );
422  if( NC_NOERR != fail )
423  {
424  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting floating_point_word_size attribute" );
425  }
426  if( att_type != NC_INT || att_len != 1 )
427  {
428  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Word size didn't have type int or size 1" );
429  }
430  fail = nc_get_att_int( ncFile, NC_GLOBAL, "floating_point_word_size", &IO_WORD_SIZE );
431  if( NC_NOERR != fail )
432  {
433  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Trouble getting word size" );
434  }
435 
436  // Exodus version
437  fail = nc_inq_att( ncFile, NC_GLOBAL, "version", &att_type, &att_len );
438  if( NC_NOERR != fail )
439  {
440  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting version attribute" );
441  }
442  if( att_type != NC_FLOAT || att_len != 1 )
443  {
444  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Version didn't have type float or size 1" );
445  }
446  // float version = temp_att->as_float(0);
447 
448  // Read in initial variables
449  int temp_dim;
450  GET_DIM( temp_dim, "num_dim", numberDimensions_loading );
451  GET_DIM( temp_dim, "num_nodes", numberNodes_loading );
452  GET_DIM( temp_dim, "num_elem", numberElements_loading );
453  GET_DIM( temp_dim, "num_el_blk", numberElementBlocks_loading );
454  GET_DIM( temp_dim, "num_fa_blk", numberFaceBlocks_loading ); // for polyhedra blocks
455  GET_DIM( temp_dim, "num_elem", numberElements_loading );
456  GET_DIM( temp_dim, "num_node_sets", numberNodeSets_loading );
457  GET_DIM( temp_dim, "num_side_sets", numberSideSets_loading );
458  GET_DIM( temp_dim, "len_string", max_str_length );
459  GET_DIM( temp_dim, "len_line", max_line_length );
460 
461  // Title; why are we even bothering if we're not going to keep it???
462  fail = nc_inq_att( ncFile, NC_GLOBAL, "title", &att_type, &att_len );
463  if( NC_NOERR != fail )
464  {
465  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting title attribute" );
466  }
467  if( att_type != NC_CHAR )
468  {
469  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: title didn't have type char" );
470  }
471  char* title = new char[att_len + 1];
472  fail = nc_get_att_text( ncFile, NC_GLOBAL, "title", title );
473  if( NC_NOERR != fail )
474  {
475  delete[] title;
476  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: trouble getting title" );
477  }
478  delete[] title;
479 
480  return MB_SUCCESS;
481 }
482 
483 ErrorCode ReadNCDF::read_nodes( const Tag* file_id_tag )
484 {
485  // Read the nodes into memory
486 
487  assert( 0 != ncFile );
488 
489  // Create a sequence to hold the node coordinates
490  // Get the current number of entities and start at the next slot
491 
492  EntityHandle node_handle = 0;
493  std::vector< double* > arrays;
494  readMeshIface->get_node_coords( 3, numberNodes_loading, MB_START_ID, node_handle, arrays );
495 
496  vertexOffset = ID_FROM_HANDLE( node_handle ) - MB_START_ID;
497 
498  // Read in the coordinates
499  int fail;
500  int coord = 0;
501  nc_inq_varid( ncFile, "coord", &coord );
502 
503  // Single var for all coords
504  size_t start[2] = { 0, 0 }, count[2] = { 1, static_cast< size_t >( numberNodes_loading ) };
505  if( coord )
506  {
507  for( int d = 0; d < numberDimensions_loading; ++d )
508  {
509  start[0] = d;
510  fail = nc_get_vara_double( ncFile, coord, start, count, arrays[d] );
511  if( NC_NOERR != fail )
512  {
513  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting " << (char)( 'x' + d ) << " coord array" );
514  }
515  }
516  }
517  // Var for each coord
518  else
519  {
520  char varname[] = "coord ";
521  for( int d = 0; d < numberDimensions_loading; ++d )
522  {
523  varname[5] = 'x' + (char)d;
524  fail = nc_inq_varid( ncFile, varname, &coord );
525  if( NC_NOERR != fail )
526  {
527  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting " << (char)( 'x' + d ) << " coord variable" );
528  }
529 
530  fail = nc_get_vara_double( ncFile, coord, start, &count[1], arrays[d] );
531  if( NC_NOERR != fail )
532  {
533  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting " << (char)( 'x' + d ) << " coord array" );
534  }
535  }
536  }
537 
538  // Zero out any coord values that are in database but not in file
539  // (e.g. if MOAB has 3D coords but file is 2D then set Z coords to zero.)
540  for( unsigned d = numberDimensions_loading; d < arrays.size(); ++d )
541  std::fill( arrays[d], arrays[d] + numberNodes_loading, 0.0 );
542 
543  if( file_id_tag )
544  {
545  Range nodes;
546  nodes.insert( node_handle, node_handle + numberNodes_loading - 1 );
547  readMeshIface->assign_ids( *file_id_tag, nodes, vertexOffset );
548  }
549 
550  return MB_SUCCESS;
551 }
552 
553 ErrorCode ReadNCDF::read_block_headers( const int* blocks_to_load, const int num_blocks )
554 {
555  // Get the element block ids; keep this in a separate list,
556  // which is not offset by blockIdOffset; this list used later for
557  // reading block connectivity
558 
559  // Get the ids of all the blocks of this file we're reading in
560  std::vector< int > block_ids( numberElementBlocks_loading );
561  int nc_block_ids = -1;
562  GET_1D_INT_VAR( "eb_prop1", nc_block_ids, block_ids );
563  if( -1 == nc_block_ids )
564  {
565  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting eb_prop1 variable" );
566  }
567 
568  int exodus_id = 1;
569 
570  // If the active_block_id_list is NULL all blocks are active.
571  if( NULL == blocks_to_load || 0 == num_blocks )
572  {
573  blocks_to_load = &block_ids[0];
574  }
575 
576  std::vector< int > new_blocks( blocks_to_load, blocks_to_load + numberElementBlocks_loading );
577 
578  std::vector< int >::iterator iter, end_iter;
579  iter = block_ids.begin();
580  end_iter = block_ids.end();
581 
582  // Read header information and initialize header-type block information
583  int temp_dim;
584  std::vector< char > temp_string_storage( max_str_length + 1 );
585  char* temp_string = &temp_string_storage[0];
586  int block_seq_id = 1;
587 
588  for( ; iter != end_iter; ++iter, block_seq_id++ )
589  {
590  int num_elements;
591 
592  GET_DIMB( temp_dim, temp_string, "num_el_in_blk%d", block_seq_id, num_elements );
593 
594  // Don't read element type string for now, since it's an attrib
595  // on the connectivity
596  // Get the entity type corresponding to this ExoII element type
597  // ExoIIElementType elem_type =
598  // ExoIIUtil::static_element_name_to_type(element_type_string);
599 
600  // Tag each element block(mesh set) with enum for ElementType (SHELL, QUAD4, ....etc)
601  ReadBlockData block_data;
602  block_data.elemType = EXOII_MAX_ELEM_TYPE;
603  block_data.blockId = *iter;
604  block_data.startExoId = exodus_id;
605  block_data.numElements = num_elements;
606 
607  // If block is in 'blocks_to_load'----load it!
608  if( std::find( new_blocks.begin(), new_blocks.end(), *iter ) != new_blocks.end() )
609  block_data.reading_in = true;
610  else
611  block_data.reading_in = false;
612 
613  blocksLoading.push_back( block_data );
614  exodus_id += num_elements;
615  }
616 
617  return MB_SUCCESS;
618 }
619 
620 ErrorCode ReadNCDF::read_face_blocks_headers()
621 {
622  // Get the ids of all the blocks of this file we're reading in
623  std::vector< int > fblock_ids( numberFaceBlocks_loading );
624  int nc_fblock_ids = -1;
625  GET_1D_INT_VAR( "fa_prop1", nc_fblock_ids, fblock_ids );
626  if( -1 == nc_fblock_ids )
627  {
628  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting fa_prop1 variable" );
629  }
630 
631  int temp_dim;
632  std::vector< char > temp_string_storage( max_str_length + 1 );
633  char* temp_string = &temp_string_storage[0];
634  // int block_seq_id = 1;
635  int exodus_id = 1;
636 
637  for( int i = 1; i <= numberFaceBlocks_loading; i++ )
638  {
639  int num_elements;
640 
641  GET_DIMB( temp_dim, temp_string, "num_fa_in_blk%d", i, num_elements );
642 
643  // Don't read element type string for now, since it's an attrib
644  // on the connectivity
645  // Get the entity type corresponding to this ExoII element type
646  // ExoIIElementType elem_type =
647  // ExoIIUtil::static_element_name_to_type(element_type_string);
648 
649  // Tag each element block(mesh set) with enum for ElementType (SHELL, QUAD4, ....etc)
650  ReadFaceBlockData fblock_data;
651  fblock_data.faceBlockId = i;
652  fblock_data.startExoId = exodus_id;
653  fblock_data.numElements = num_elements;
654 
655  faceBlocksLoading.push_back( fblock_data );
656  exodus_id += num_elements;
657  }
658  return MB_SUCCESS;
659 }
660 ErrorCode ReadNCDF::read_polyhedra_faces()
661 {
662  int temp_dim;
663  std::vector< char > temp_string_storage( max_str_length + 1 );
664  char* temp_string = &temp_string_storage[0];
665  nodesInLoadedBlocks.resize( numberNodes_loading + 1 );
666  size_t start[1] = { 0 }, count[1]; // one dim arrays only here!
667  std::vector< ReadFaceBlockData >::iterator this_it;
668  this_it = faceBlocksLoading.begin();
669 
670  int fblock_seq_id = 1;
671  std::vector< int > dims;
672  int nc_var, fail;
673 
674  for( ; this_it != faceBlocksLoading.end(); ++this_it, fblock_seq_id++ )
675  {
676  int num_fa_in_blk;
677  GET_DIMB( temp_dim, temp_string, "num_fa_in_blk%d", fblock_seq_id, num_fa_in_blk );
678  int num_nod_per_fa;
679  GET_DIMB( temp_dim, temp_string, "num_nod_per_fa%d", fblock_seq_id, num_nod_per_fa );
680  // Get the ncdf connect variable and the element type
681  INS_ID( temp_string, "fbconn%d", fblock_seq_id, max_str_length + 1 );
682  GET_VAR( temp_string, nc_var, dims );
683  std::vector< int > fbconn;
684  fbconn.resize( num_nod_per_fa );
685  count[0] = num_nod_per_fa;
686 
687  fail = nc_get_vara_int( ncFile, nc_var, start, count, &fbconn[0] );
688  if( NC_NOERR != fail )
689  {
690  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting face polyhedra connectivity " );
691  }
692  std::vector< int > fbepecnt;
693  fbepecnt.resize( num_fa_in_blk );
694  INS_ID( temp_string, "fbepecnt%d", fblock_seq_id, max_str_length + 1 );
695  GET_VAR( temp_string, nc_var, dims );
696  count[0] = num_fa_in_blk;
697  fail = nc_get_vara_int( ncFile, nc_var, start, count, &fbepecnt[0] );
698  if( NC_NOERR != fail )
699  {
700  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting face polyhedra connectivity " );
701  }
702  // now we can create some polygons
703  std::vector< EntityHandle > polyConn;
704  int ix = 0;
705  for( int i = 0; i < num_fa_in_blk; i++ )
706  {
707  polyConn.resize( fbepecnt[i] );
708  for( int j = 0; j < fbepecnt[i]; j++ )
709  {
710  nodesInLoadedBlocks[fbconn[ix]] = 1;
711  polyConn[j] = vertexOffset + fbconn[ix++];
712  }
713  EntityHandle newp;
714  /*
715  * ErrorCode create_element(const EntityType type,
716  const EntityHandle *connectivity,
717  const int num_vertices,
718  EntityHandle &element_handle)
719  */
720  if( mdbImpl->create_element( MBPOLYGON, &polyConn[0], fbepecnt[i], newp ) != MB_SUCCESS ) return MB_FAILURE;
721 
722  // add the element in order
723  polyfaces.push_back( newp );
724  }
725  }
726  return MB_SUCCESS;
727 }
728 
729 ErrorCode ReadNCDF::read_elements( const Tag* file_id_tag )
730 {
731  // Read in elements
732 
733  int result = 0;
734 
735  // Initialize the nodeInLoadedBlocks vector
736  if( nodesInLoadedBlocks.size() < (size_t)( numberNodes_loading + 1 ) )
737  nodesInLoadedBlocks.resize( numberNodes_loading + 1 ); // this could be repeated?
738  memset( &nodesInLoadedBlocks[0], 0, ( numberNodes_loading + 1 ) * sizeof( char ) );
739 
740  std::vector< ReadBlockData >::iterator this_it;
741  this_it = blocksLoading.begin();
742 
743  int temp_dim;
744  std::vector< char > temp_string_storage( max_str_length + 1 );
745  char* temp_string = &temp_string_storage[0];
746 
747  int nc_var;
748  int block_seq_id = 1;
749  std::vector< int > dims;
750  size_t start[2] = { 0, 0 }, count[2];
751 
752  for( ; this_it != blocksLoading.end(); ++this_it, block_seq_id++ )
753  {
754  // If this block isn't to be read in --- continue
755  if( !( *this_it ).reading_in ) continue;
756 
757  // Get some information about this block
758  int block_id = ( *this_it ).blockId;
759  EntityHandle* conn = 0;
760 
761  // Get the ncdf connect variable and the element type
762  INS_ID( temp_string, "connect%d", block_seq_id, max_str_length + 1 );
763  GET_VAR( temp_string, nc_var, dims );
764  if( -1 == nc_var || 0 == nc_var )
765  { // try other var, for polyhedra blocks
766  // it could be polyhedra block, defined by fbconn and NFACED attribute
767  INS_ID( temp_string, "facconn%d", block_seq_id, max_str_length + 1 );
768  GET_VAR( temp_string, nc_var, dims );
769 
770  if( -1 == nc_var || 0 == nc_var )
771  {
772  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting connec or faccon variable" );
773  }
774  }
775  nc_type att_type;
776  size_t att_len;
777  int fail = nc_inq_att( ncFile, nc_var, "elem_type", &att_type, &att_len );
778  if( NC_NOERR != fail )
779  {
780  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting elem type attribute" );
781  }
782  std::vector< char > dum_str( att_len + 1 );
783  fail = nc_get_att_text( ncFile, nc_var, "elem_type", &dum_str[0] );
784  if( NC_NOERR != fail )
785  {
786  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting elem type" );
787  }
788  ExoIIElementType elem_type = ExoIIUtil::static_element_name_to_type( &dum_str[0] );
789  ( *this_it ).elemType = elem_type;
790 
791  int verts_per_element = ExoIIUtil::VerticesPerElement[( *this_it ).elemType];
792  int number_nodes = ( *this_it ).numElements * verts_per_element;
793  const EntityType mb_type = ExoIIUtil::ExoIIElementMBEntity[( *this_it ).elemType];
794 
795  if( mb_type == MBPOLYGON )
796  {
797  // need to read another variable, num_nod_per_el, and
798  int num_nodes_poly_conn;
799  GET_DIMB( temp_dim, temp_string, "num_nod_per_el%d", block_seq_id, num_nodes_poly_conn );
800  // read the connec1 array, which is of dimensions num_nodes_poly_conn (only one )
801  std::vector< int > connec;
802  connec.resize( num_nodes_poly_conn );
803  count[0] = num_nodes_poly_conn;
804 
805  fail = nc_get_vara_int( ncFile, nc_var, start, count, &connec[0] );
806  if( NC_NOERR != fail )
807  {
808  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting polygon connectivity " );
809  }
810  // ebepecnt1
811  std::vector< int > ebec;
812  ebec.resize( this_it->numElements );
813  // Get the ncdf connect variable and the element type
814  INS_ID( temp_string, "ebepecnt%d", block_seq_id, max_str_length + 1 );
815  GET_VAR( temp_string, nc_var, dims );
816  count[0] = this_it->numElements;
817  fail = nc_get_vara_int( ncFile, nc_var, start, count, &ebec[0] );
818  if( NC_NOERR != fail )
819  {
820  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting polygon nodes per elem " );
821  }
822  EntityHandle ms_handle;
823  if( mdbImpl->create_meshset( MESHSET_SET | MESHSET_TRACK_OWNER, ms_handle ) != MB_SUCCESS )
824  return MB_FAILURE;
825  // create polygons one by one, and put them in the list
826  // also need to read the number of nodes for each polygon, then create
827  std::vector< EntityHandle > polyConn;
828  int ix = 0;
829  for( int i = 0; i < this_it->numElements; i++ )
830  {
831  polyConn.resize( ebec[i] );
832  for( int j = 0; j < ebec[i]; j++ )
833  {
834  nodesInLoadedBlocks[connec[ix]] = 1;
835  polyConn[j] = vertexOffset + connec[ix++];
836  }
837  EntityHandle newp;
838  /*
839  * ErrorCode create_element(const EntityType type,
840  const EntityHandle *connectivity,
841  const int num_vertices,
842  EntityHandle &element_handle)
843  */
844  if( mdbImpl->create_element( MBPOLYGON, &polyConn[0], ebec[i], newp ) != MB_SUCCESS ) return MB_FAILURE;
845 
846  // add the element in order
847  this_it->polys.push_back( newp );
848  if( mdbImpl->add_entities( ms_handle, &newp, 1 ) != MB_SUCCESS ) return MB_FAILURE;
849  }
850  // Set the block id with an offset
851  if( mdbImpl->tag_set_data( mMaterialSetTag, &ms_handle, 1, &block_id ) != MB_SUCCESS ) return MB_FAILURE;
852  if( mdbImpl->tag_set_data( mGlobalIdTag, &ms_handle, 1, &block_id ) != MB_SUCCESS ) return MB_FAILURE;
853  }
854  else if( mb_type == MBPOLYHEDRON )
855  {
856  // need to read another variable, num_fac_per_el
857  int num_fac_per_el;
858  GET_DIMB( temp_dim, temp_string, "num_fac_per_el%d", block_seq_id, num_fac_per_el );
859  // read the fbconn array, which is of dimensions num_nod_per_fa (only one dimension)
860  std::vector< int > facconn;
861  facconn.resize( num_fac_per_el );
862  count[0] = num_fac_per_el;
863 
864  fail = nc_get_vara_int( ncFile, nc_var, start, count, &facconn[0] );
865  if( NC_NOERR != fail )
866  {
867  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting polygon connectivity " );
868  }
869  // ebepecnt1
870  std::vector< int > ebepecnt;
871  ebepecnt.resize( this_it->numElements );
872  // Get the ncdf connect variable and the element type
873  INS_ID( temp_string, "ebepecnt%d", block_seq_id, max_str_length + 1 );
874  GET_VAR( temp_string, nc_var, dims );
875  count[0] = this_it->numElements;
876  fail = nc_get_vara_int( ncFile, nc_var, start, count, &ebepecnt[0] );
877  if( NC_NOERR != fail )
878  {
879  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting polygon nodes per elem " );
880  }
881  EntityHandle ms_handle;
882  if( mdbImpl->create_meshset( MESHSET_SET | MESHSET_TRACK_OWNER, ms_handle ) != MB_SUCCESS )
883  return MB_FAILURE;
884  // create polygons one by one, and put them in the list
885  // also need to read the number of nodes for each polygon, then create
886  std::vector< EntityHandle > polyConn;
887  int ix = 0;
888  for( int i = 0; i < this_it->numElements; i++ )
889  {
890  polyConn.resize( ebepecnt[i] );
891  for( int j = 0; j < ebepecnt[i]; j++ )
892  {
893  polyConn[j] = polyfaces[facconn[ix++] - 1];
894  }
895  EntityHandle newp;
896  /*
897  * ErrorCode create_element(const EntityType type,
898  const EntityHandle *connectivity,
899  const int num_vertices,
900  EntityHandle &element_handle)
901  */
902  if( mdbImpl->create_element( MBPOLYHEDRON, &polyConn[0], ebepecnt[i], newp ) != MB_SUCCESS )
903  return MB_FAILURE;
904 
905  // add the element in order
906  this_it->polys.push_back( newp );
907  if( mdbImpl->add_entities( ms_handle, &newp, 1 ) != MB_SUCCESS ) return MB_FAILURE;
908  }
909  // Set the block id with an offset
910  if( mdbImpl->tag_set_data( mMaterialSetTag, &ms_handle, 1, &block_id ) != MB_SUCCESS ) return MB_FAILURE;
911  if( mdbImpl->tag_set_data( mGlobalIdTag, &ms_handle, 1, &block_id ) != MB_SUCCESS ) return MB_FAILURE;
912  }
913  else // this is regular
914  {
915  // Allocate an array to read in connectivity data
916  readMeshIface->get_element_connect( this_it->numElements, verts_per_element, mb_type, this_it->startExoId,
917  this_it->startMBId, conn );
918 
919  // Create a range for this sequence of elements
920  EntityHandle start_range, end_range;
921  start_range = ( *this_it ).startMBId;
922  end_range = start_range + ( *this_it ).numElements - 1;
923 
924  Range new_range( start_range, end_range );
925  // Range<EntityHandle> new_range((*this_it).startMBId,
926  // (*this_it).startMBId + (*this_it).numElements - 1);
927 
928  // Create a MBSet for this block and set the material tag
929 
930  EntityHandle ms_handle;
931  if( mdbImpl->create_meshset( MESHSET_SET | MESHSET_TRACK_OWNER, ms_handle ) != MB_SUCCESS )
932  return MB_FAILURE;
933 
934  if( mdbImpl->add_entities( ms_handle, new_range ) != MB_SUCCESS ) return MB_FAILURE;
935 
936  int mid_nodes[4];
937  CN::HasMidNodes( mb_type, verts_per_element, mid_nodes );
938  if( mdbImpl->tag_set_data( mHasMidNodesTag, &ms_handle, 1, mid_nodes ) != MB_SUCCESS ) return MB_FAILURE;
939 
940  // Just a check because the following code won't work if this case fails
941  assert( sizeof( EntityHandle ) >= sizeof( int ) );
942 
943  // tmp_ptr is of type int* and points at the same place as conn
944  int* tmp_ptr = reinterpret_cast< int* >( conn );
945 
946  // Read the connectivity into that memory, this will take only part of the array
947  // 1/2 if sizeof(EntityHandle) == 64 bits.
948  count[0] = this_it->numElements;
949  count[1] = verts_per_element;
950  fail = nc_get_vara_int( ncFile, nc_var, start, count, tmp_ptr );
951  if( NC_NOERR != fail )
952  {
953  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting connectivity" );
954  }
955 
956  // Convert from exodus indices to vertex handles.
957  // Iterate backwards in case handles are larger than ints.
958  for( int i = number_nodes - 1; i >= 0; --i )
959  {
960  if( (unsigned)tmp_ptr[i] >= nodesInLoadedBlocks.size() )
961  {
962  MB_SET_ERR( MB_FAILURE, "Invalid node ID in block connectivity" );
963  }
964  nodesInLoadedBlocks[tmp_ptr[i]] = 1;
965  conn[i] = static_cast< EntityHandle >( tmp_ptr[i] ) + vertexOffset;
966  }
967 
968  // Adjust connectivity order if necessary
969  const int* reorder = exodus_elem_order_map[mb_type][verts_per_element];
970  if( reorder ) ReadUtilIface::reorder( reorder, conn, this_it->numElements, verts_per_element );
971 
972  readMeshIface->update_adjacencies( ( *this_it ).startMBId, ( *this_it ).numElements,
973  ExoIIUtil::VerticesPerElement[( *this_it ).elemType], conn );
974 
975  if( result == -1 )
976  {
977  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: error getting element connectivity for block " << block_id );
978  }
979 
980  // Set the block id with an offset
981  if( mdbImpl->tag_set_data( mMaterialSetTag, &ms_handle, 1, &block_id ) != MB_SUCCESS ) return MB_FAILURE;
982  if( mdbImpl->tag_set_data( mGlobalIdTag, &ms_handle, 1, &block_id ) != MB_SUCCESS ) return MB_FAILURE;
983 
984  if( file_id_tag )
985  {
986  Range range;
987  range.insert( this_it->startMBId, this_it->startMBId + this_it->numElements - 1 );
988  readMeshIface->assign_ids( *file_id_tag, range, this_it->startExoId );
989  }
990  } // end regular block (not polygon)
991  }
992 
993  return MB_SUCCESS;
994 }
995 
996 ErrorCode ReadNCDF::read_global_ids()
997 {
998  // Read in the map from the exodus file
999  std::vector< int > ids( std::max( numberElements_loading, numberNodes_loading ) );
1000 
1001  int varid = -1;
1002  GET_1D_INT_VAR( "elem_map", varid, ids );
1003  if( -1 != varid )
1004  {
1005  std::vector< ReadBlockData >::iterator iter;
1006  int id_pos = 0;
1007  for( iter = blocksLoading.begin(); iter != blocksLoading.end(); ++iter )
1008  {
1009  if( iter->reading_in )
1010  {
1011  if( iter->elemType != EXOII_POLYGON && iter->elemType != EXOII_POLYHEDRON )
1012  {
1013  if( iter->startMBId != 0 )
1014  {
1015  Range range( iter->startMBId, iter->startMBId + iter->numElements - 1 );
1016  ErrorCode error = mdbImpl->tag_set_data( mGlobalIdTag, range, &ids[id_pos] );
1017  if( error != MB_SUCCESS ) return error;
1018  id_pos += iter->numElements;
1019  }
1020  else
1021  return MB_FAILURE;
1022  }
1023  else // polygons or polyhedrons; use id from elements
1024  {
1025  for( std::vector< EntityHandle >::iterator eit = iter->polys.begin(); eit != iter->polys.end();
1026  eit++, id_pos++ )
1027  {
1028  EntityHandle peh = *eit;
1029  if( mdbImpl->tag_set_data( mGlobalIdTag, &peh, 1, &ids[id_pos] ) != MB_SUCCESS )
1030  return MB_FAILURE;
1031  }
1032  }
1033  }
1034  }
1035  }
1036 
1037  // Read in node map next
1038  varid = -1;
1039  GET_1D_INT_VAR( "node_num_map", varid, ids );
1040  if( -1 != varid )
1041  {
1042  Range range( MB_START_ID + vertexOffset, MB_START_ID + vertexOffset + numberNodes_loading - 1 );
1043  ErrorCode error = mdbImpl->tag_set_data( mGlobalIdTag, range, &ids[0] );
1044  if( MB_SUCCESS != error ) MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem setting node global ids" );
1045  }
1046 
1047  return MB_SUCCESS;
1048 }
1049 
1050 ErrorCode ReadNCDF::read_nodesets()
1051 {
1052  // Read in the nodesets for the model
1053 
1054  if( 0 == numberNodeSets_loading ) return MB_SUCCESS;
1055  std::vector< int > id_array( numberNodeSets_loading );
1056 
1057  // Read in the nodeset ids
1058  int nc_var;
1059  GET_1D_INT_VAR( "ns_prop1", nc_var, id_array );
1060  if( -1 == nc_var )
1061  {
1062  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting ns_prop1 variable" );
1063  }
1064 
1065  // Use a vector of ints to read node handles
1066  std::vector< int > node_handles;
1067 
1068  int i;
1069  std::vector< char > temp_string_storage( max_str_length + 1 );
1070  char* temp_string = &temp_string_storage[0];
1071  int temp_dim;
1072  for( i = 0; i < numberNodeSets_loading; i++ )
1073  {
1074  // Get nodeset parameters
1075  int number_nodes_in_set;
1076  int number_dist_factors_in_set;
1077 
1078  GET_DIMB( temp_dim, temp_string, "num_nod_ns%d", i + 1, number_nodes_in_set );
1079  GET_DIMB( temp_dim, temp_string, "num_df_ns%d", i + 1, number_dist_factors_in_set );
1080 
1081  // Need to new a vector to store dist. factors
1082  // This vector * gets stored as a tag on the sideset meshset
1083  std::vector< double > temp_dist_factor_vector( number_nodes_in_set );
1084  if( number_dist_factors_in_set != 0 )
1085  {
1086  INS_ID( temp_string, "dist_fact_ns%d", i + 1, max_str_length + 1 );
1087  GET_1D_DBL_VAR( temp_string, temp_dim, temp_dist_factor_vector );
1088  if( -1 == temp_dim )
1089  {
1090  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting dist fact variable" );
1091  }
1092  }
1093 
1094  // Size new arrays and get ids and distribution factors
1095  if( node_handles.size() < (unsigned int)number_nodes_in_set )
1096  {
1097  node_handles.reserve( number_nodes_in_set );
1098  node_handles.resize( number_nodes_in_set );
1099  }
1100 
1101  INS_ID( temp_string, "node_ns%d", i + 1, max_str_length + 1 );
1102  int temp_var = -1;
1103  GET_1D_INT_VAR( temp_string, temp_var, node_handles );
1104  if( -1 == temp_var )
1105  {
1106  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting nodeset node variable" );
1107  }
1108 
1109  // Maybe there is already a nodesets meshset here we can append to
1110  Range child_meshsets;
1111  if( mdbImpl->get_entities_by_handle( 0, child_meshsets ) != MB_SUCCESS ) return MB_FAILURE;
1112 
1113  child_meshsets = subtract( child_meshsets, initRange );
1114 
1115  Range::iterator iter, end_iter;
1116  iter = child_meshsets.begin();
1117  end_iter = child_meshsets.end();
1118 
1119  EntityHandle ns_handle = 0;
1120  for( ; iter != end_iter; ++iter )
1121  {
1122  int nodeset_id;
1123  if( mdbImpl->tag_get_data( mDirichletSetTag, &( *iter ), 1, &nodeset_id ) != MB_SUCCESS ) continue;
1124 
1125  if( id_array[i] == nodeset_id )
1126  {
1127  // Found the meshset
1128  ns_handle = *iter;
1129  break;
1130  }
1131  }
1132 
1133  std::vector< EntityHandle > nodes_of_nodeset;
1134  if( ns_handle )
1135  if( mdbImpl->get_entities_by_handle( ns_handle, nodes_of_nodeset, true ) != MB_SUCCESS ) return MB_FAILURE;
1136 
1137  // Make these into entity handles
1138  // TODO: could we have read it into EntityHandle sized array in the first place?
1139  int j, temp;
1140  std::vector< EntityHandle > nodes;
1141  std::vector< double > dist_factor_vector;
1142  for( j = 0; j < number_nodes_in_set; j++ )
1143  {
1144  // See if this node is one we're currently reading in
1145  if( nodesInLoadedBlocks[node_handles[j]] == 1 )
1146  {
1147  // Make sure that it already isn't in a nodeset
1148  unsigned int node_id = CREATE_HANDLE( MBVERTEX, node_handles[j] + vertexOffset, temp );
1149  if( !ns_handle ||
1150  std::find( nodes_of_nodeset.begin(), nodes_of_nodeset.end(), node_id ) == nodes_of_nodeset.end() )
1151  {
1152  nodes.push_back( node_id );
1153 
1154  if( number_dist_factors_in_set != 0 ) dist_factor_vector.push_back( temp_dist_factor_vector[j] );
1155  }
1156  }
1157  }
1158 
1159  // No nodes to add
1160  if( nodes.empty() ) continue;
1161 
1162  // If there was no meshset found --> create one
1163  if( ns_handle == 0 )
1164  {
1165  if( mdbImpl->create_meshset( MESHSET_ORDERED | MESHSET_TRACK_OWNER, ns_handle ) != MB_SUCCESS )
1166  return MB_FAILURE;
1167 
1168  // Set a tag signifying dirichlet bc
1169  // TODO: create this tag another way
1170 
1171  int nodeset_id = id_array[i];
1172  if( mdbImpl->tag_set_data( mDirichletSetTag, &ns_handle, 1, &nodeset_id ) != MB_SUCCESS ) return MB_FAILURE;
1173  if( mdbImpl->tag_set_data( mGlobalIdTag, &ns_handle, 1, &nodeset_id ) != MB_SUCCESS ) return MB_FAILURE;
1174 
1175  if( !dist_factor_vector.empty() )
1176  {
1177  int size = dist_factor_vector.size();
1178  const void* data = &dist_factor_vector[0];
1179  if( mdbImpl->tag_set_by_ptr( mDistFactorTag, &ns_handle, 1, &data, &size ) != MB_SUCCESS )
1180  return MB_FAILURE;
1181  }
1182  }
1183  else if( !dist_factor_vector.empty() )
1184  {
1185  // Append dist factors to vector
1186  const void* ptr = 0;
1187  int size = 0;
1188  if( mdbImpl->tag_get_by_ptr( mDistFactorTag, &ns_handle, 1, &ptr, &size ) != MB_SUCCESS ) return MB_FAILURE;
1189 
1190  const double* data = reinterpret_cast< const double* >( ptr );
1191  dist_factor_vector.reserve( dist_factor_vector.size() + size );
1192  std::copy( data, data + size, std::back_inserter( dist_factor_vector ) );
1193  size = dist_factor_vector.size();
1194  ptr = &dist_factor_vector[0];
1195  if( mdbImpl->tag_set_by_ptr( mDistFactorTag, &ns_handle, 1, &ptr, &size ) != MB_SUCCESS ) return MB_FAILURE;
1196  }
1197 
1198  // Add the nodes to the meshset
1199  if( mdbImpl->add_entities( ns_handle, &nodes[0], nodes.size() ) != MB_SUCCESS ) return MB_FAILURE;
1200  }
1201 
1202  return MB_SUCCESS;
1203 }
1204 
1205 ErrorCode ReadNCDF::read_sidesets()
1206 {
1207  // Uhhh if you read the same file (previously_read==true) then blow away the
1208  // sidesets pertaining to this file? How do you do that? If you read a
1209  // new file make sure all the offsets are correct.
1210 
1211  // If not loading any sidesets -- exit
1212  if( 0 == numberSideSets_loading ) return MB_SUCCESS;
1213 
1214  // Read in the sideset ids
1215  std::vector< int > id_array( numberSideSets_loading );
1216  int temp_var = -1;
1217  GET_1D_INT_VAR( "ss_prop1", temp_var, id_array );
1218  if( -1 == temp_var )
1219  {
1220  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting ss_prop1 variable" );
1221  }
1222 
1223  // Create a side set for each one
1224  int number_sides_in_set;
1225  int number_dist_factors_in_set;
1226 
1227  // Maybe there is already a sidesets meshset here we can append to
1228  Range child_meshsets;
1229  if( mdbImpl->get_entities_by_type( 0, MBENTITYSET, child_meshsets ) != MB_SUCCESS ) return MB_FAILURE;
1230 
1231  child_meshsets = subtract( child_meshsets, initRange );
1232 
1233  Range::iterator iter, end_iter;
1234 
1235  int i;
1236  std::vector< char > temp_string_storage( max_str_length + 1 );
1237  char* temp_string = &temp_string_storage[0];
1238  int temp_dim;
1239  for( i = 0; i < numberSideSets_loading; i++ )
1240  {
1241  // Get sideset parameters
1242  GET_DIMB( temp_dim, temp_string, "num_side_ss%d", i + 1, number_sides_in_set );
1243  GET_DIMB( temp_dim, temp_string, "num_df_ss%d", i + 1, number_dist_factors_in_set );
1244 
1245  // Size new arrays and get element and side lists
1246  std::vector< int > side_list( number_sides_in_set );
1247  std::vector< int > element_list( number_sides_in_set );
1248  INS_ID( temp_string, "side_ss%d", i + 1, max_str_length + 1 );
1249  temp_var = -1;
1250  GET_1D_INT_VAR( temp_string, temp_var, side_list );
1251  if( -1 == temp_var )
1252  {
1253  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting sideset side variable" );
1254  }
1255 
1256  INS_ID( temp_string, "elem_ss%d", i + 1, max_str_length + 1 );
1257  temp_var = -1;
1258  GET_1D_INT_VAR( temp_string, temp_var, element_list );
1259  if( -1 == temp_var )
1260  {
1261  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting sideset elem variable" );
1262  }
1263 
1264  std::vector< double > temp_dist_factor_vector;
1265  std::vector< EntityHandle > entities_to_add, reverse_entities;
1266  // Create the sideset entities
1267  if( create_ss_elements( &element_list[0], &side_list[0], number_sides_in_set, number_dist_factors_in_set,
1268  entities_to_add, reverse_entities, temp_dist_factor_vector, i + 1 ) != MB_SUCCESS )
1269  return MB_FAILURE;
1270 
1271  // If there are elements to add
1272  if( !entities_to_add.empty() || !reverse_entities.empty() )
1273  {
1274  iter = child_meshsets.begin();
1275  end_iter = child_meshsets.end();
1276 
1277  EntityHandle ss_handle = 0;
1278  for( ; iter != end_iter; ++iter )
1279  {
1280  int sideset_id;
1281  if( mdbImpl->tag_get_data( mNeumannSetTag, &( *iter ), 1, &sideset_id ) != MB_SUCCESS ) continue;
1282 
1283  if( id_array[i] == sideset_id )
1284  {
1285  // Found the meshset
1286  ss_handle = *iter;
1287  break;
1288  }
1289  }
1290 
1291  // If we didn't find a sideset already
1292  if( ss_handle == 0 )
1293  {
1294  if( mdbImpl->create_meshset( MESHSET_ORDERED | MESHSET_TRACK_OWNER, ss_handle ) != MB_SUCCESS )
1295  return MB_FAILURE;
1296 
1297  if( ss_handle == 0 ) return MB_FAILURE;
1298 
1299  int sideset_id = id_array[i];
1300  if( mdbImpl->tag_set_data( mNeumannSetTag, &ss_handle, 1, &sideset_id ) != MB_SUCCESS )
1301  return MB_FAILURE;
1302  if( mdbImpl->tag_set_data( mGlobalIdTag, &ss_handle, 1, &sideset_id ) != MB_SUCCESS ) return MB_FAILURE;
1303 
1304  if( !reverse_entities.empty() )
1305  {
1306  // Also make a reverse set to put in this set
1307  EntityHandle reverse_set;
1308  if( mdbImpl->create_meshset( MESHSET_SET | MESHSET_TRACK_OWNER, reverse_set ) != MB_SUCCESS )
1309  return MB_FAILURE;
1310 
1311  // Add the reverse set to the sideset set and the entities to the reverse set
1312  ErrorCode result = mdbImpl->add_entities( ss_handle, &reverse_set, 1 );
1313  if( MB_SUCCESS != result ) return result;
1314 
1315  result = mdbImpl->add_entities( reverse_set, &reverse_entities[0], reverse_entities.size() );
1316  if( MB_SUCCESS != result ) return result;
1317 
1318  // Set the reverse tag
1319  Tag sense_tag;
1320  int dum_sense = 0;
1321  result = mdbImpl->tag_get_handle( "NEUSET_SENSE", 1, MB_TYPE_INTEGER, sense_tag,
1322  MB_TAG_SPARSE | MB_TAG_CREAT, &dum_sense );
1323  if( result != MB_SUCCESS ) return result;
1324  dum_sense = -1;
1325  result = mdbImpl->tag_set_data( sense_tag, &reverse_set, 1, &dum_sense );
1326  if( result != MB_SUCCESS ) return result;
1327  }
1328  }
1329 
1330  if( mdbImpl->add_entities( ss_handle, ( entities_to_add.empty() ) ? NULL : &entities_to_add[0],
1331  entities_to_add.size() ) != MB_SUCCESS )
1332  return MB_FAILURE;
1333 
1334  // Distribution factor stuff
1335  if( number_dist_factors_in_set )
1336  {
1337  // If this sideset does not already has a distribution factor array...set one
1338  const void* ptr = 0;
1339  int size = 0;
1340  if( MB_SUCCESS == mdbImpl->tag_get_by_ptr( mDistFactorTag, &ss_handle, 1, &ptr, &size ) )
1341  {
1342  const double* data = reinterpret_cast< const double* >( ptr );
1343  std::copy( data, data + size, std::back_inserter( temp_dist_factor_vector ) );
1344  }
1345 
1346  ptr = &temp_dist_factor_vector[0];
1347  size = temp_dist_factor_vector.size();
1348  if( mdbImpl->tag_set_by_ptr( mDistFactorTag, &ss_handle, 1, &ptr, &size ) != MB_SUCCESS )
1349  return MB_FAILURE;
1350  }
1351  }
1352  }
1353 
1354  return MB_SUCCESS;
1355 }
1356 
1357 ErrorCode ReadNCDF::create_ss_elements( int* element_ids,
1358  int* side_list,
1359  int num_sides,
1360  int num_dist_factors,
1361  std::vector< EntityHandle >& entities_to_add,
1362  std::vector< EntityHandle >& reverse_entities,
1363  std::vector< double >& dist_factor_vector,
1364  int ss_seq_id )
1365 {
1366  // Determine entity type from element id
1367 
1368  // If there are dist. factors, create a vector to hold the array
1369  // and place this array as a tag onto the sideset meshset
1370 
1371  std::vector< double > temp_dist_factor_vector( num_dist_factors );
1372  std::vector< char > temp_string_storage( max_str_length + 1 );
1373  char* temp_string = &temp_string_storage[0];
1374  int temp_var;
1375  if( num_dist_factors )
1376  {
1377  INS_ID( temp_string, "dist_fact_ss%d", ss_seq_id, max_str_length + 1 );
1378  temp_var = -1;
1379  GET_1D_DBL_VAR( temp_string, temp_var, temp_dist_factor_vector );
1380  if( -1 == temp_var )
1381  {
1382  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting dist fact variable" );
1383  }
1384  }
1385 
1386  EntityType subtype;
1387  int num_side_nodes, num_elem_nodes;
1388  const EntityHandle* nodes;
1389  std::vector< EntityHandle > connectivity;
1390  int side_node_idx[32];
1391 
1392  int df_index = 0;
1393  int sense = 0;
1394  for( int i = 0; i < num_sides; i++ )
1395  {
1396  ExoIIElementType exoii_type;
1397  ReadBlockData block_data;
1398  block_data.elemType = EXOII_MAX_ELEM_TYPE;
1399 
1400  if( find_side_element_type( element_ids[i], exoii_type, block_data, df_index, side_list[i] ) != MB_SUCCESS )
1401  continue; // Isn't being read in this time
1402 
1403  EntityType type = ExoIIUtil::ExoIIElementMBEntity[exoii_type];
1404 
1405  EntityHandle ent_handle = element_ids[i] - block_data.startExoId + block_data.startMBId;
1406 
1407  const int side_num = side_list[i] - 1;
1408  if( type == MBHEX )
1409  {
1410  // Get the nodes of the element
1411  if( mdbImpl->get_connectivity( ent_handle, nodes, num_elem_nodes ) != MB_SUCCESS ) return MB_FAILURE;
1412 
1413  CN::SubEntityNodeIndices( type, num_elem_nodes, 2, side_num, subtype, num_side_nodes, side_node_idx );
1414  if( num_side_nodes <= 0 ) return MB_FAILURE;
1415 
1416  connectivity.resize( num_side_nodes );
1417  for( int k = 0; k < num_side_nodes; ++k )
1418  connectivity[k] = nodes[side_node_idx[k]];
1419 
1420  if( MB_SUCCESS != create_sideset_element( connectivity, subtype, ent_handle, sense ) ) return MB_FAILURE;
1421  if( 1 == sense )
1422  entities_to_add.push_back( ent_handle );
1423  else
1424  reverse_entities.push_back( ent_handle );
1425 
1426  // Read in distribution factor array
1427  if( num_dist_factors )
1428  {
1429  for( int k = 0; k < 4; k++ )
1430  dist_factor_vector.push_back( temp_dist_factor_vector[df_index++] );
1431  }
1432  }
1433 
1434  // If it is a Tet
1435  else if( type == MBTET )
1436  {
1437  // Get the nodes of the element
1438  if( mdbImpl->get_connectivity( ent_handle, nodes, num_elem_nodes ) != MB_SUCCESS ) return MB_FAILURE;
1439 
1440  CN::SubEntityNodeIndices( type, num_elem_nodes, 2, side_num, subtype, num_side_nodes, side_node_idx );
1441  if( num_side_nodes <= 0 ) return MB_FAILURE;
1442 
1443  connectivity.resize( num_side_nodes );
1444  for( int k = 0; k < num_side_nodes; ++k )
1445  connectivity[k] = nodes[side_node_idx[k]];
1446 
1447  if( MB_SUCCESS != create_sideset_element( connectivity, subtype, ent_handle, sense ) ) return MB_FAILURE;
1448  if( 1 == sense )
1449  entities_to_add.push_back( ent_handle );
1450  else
1451  reverse_entities.push_back( ent_handle );
1452 
1453  // Read in distribution factor array
1454  if( num_dist_factors )
1455  {
1456  for( int k = 0; k < 3; k++ )
1457  dist_factor_vector.push_back( temp_dist_factor_vector[df_index++] );
1458  }
1459  }
1460  else if( type == MBQUAD && exoii_type >= EXOII_SHELL && exoii_type <= EXOII_SHELL9 )
1461  {
1462  // ent_handle = CREATE_HANDLE(MBQUAD, base_id, error);
1463 
1464  // Just use this quad
1465  if( side_list[i] == 1 )
1466  {
1467  if( 1 == sense )
1468  entities_to_add.push_back( ent_handle );
1469  else
1470  reverse_entities.push_back( ent_handle );
1471 
1472  if( num_dist_factors )
1473  {
1474  for( int k = 0; k < 4; k++ )
1475  dist_factor_vector.push_back( temp_dist_factor_vector[df_index++] );
1476  }
1477 
1478  continue;
1479  }
1480  else if( side_list[i] == 2 )
1481  {
1482  reverse_entities.push_back( ent_handle );
1483 
1484  if( num_dist_factors )
1485  {
1486  for( int k = 0; k < 4; k++ )
1487  dist_factor_vector.push_back( temp_dist_factor_vector[df_index++] );
1488  }
1489 
1490  continue;
1491  }
1492  else
1493  {
1494  // Get the nodes of the element
1495  if( mdbImpl->get_connectivity( ent_handle, nodes, num_elem_nodes ) != MB_SUCCESS ) return MB_FAILURE;
1496 
1497  CN::SubEntityNodeIndices( type, num_elem_nodes, 1, side_num - 2, subtype, num_side_nodes,
1498  side_node_idx );
1499  if( num_side_nodes <= 0 ) return MB_FAILURE;
1500 
1501  connectivity.resize( num_side_nodes );
1502  for( int k = 0; k < num_side_nodes; ++k )
1503  connectivity[k] = nodes[side_node_idx[k]];
1504 
1505  if( MB_SUCCESS != create_sideset_element( connectivity, subtype, ent_handle, sense ) )
1506  return MB_FAILURE;
1507  if( 1 == sense )
1508  entities_to_add.push_back( ent_handle );
1509  else
1510  reverse_entities.push_back( ent_handle );
1511 
1512  if( num_dist_factors )
1513  {
1514  for( int k = 0; k < 2; k++ )
1515  dist_factor_vector.push_back( temp_dist_factor_vector[df_index++] );
1516  }
1517  }
1518  }
1519  // If it is a Quad
1520  else if( type == MBQUAD )
1521  {
1522  // Get the nodes of the element
1523  if( mdbImpl->get_connectivity( ent_handle, nodes, num_elem_nodes ) != MB_SUCCESS ) return MB_FAILURE;
1524 
1525  CN::SubEntityNodeIndices( type, num_elem_nodes, 1, side_num, subtype, num_side_nodes, side_node_idx );
1526  if( num_side_nodes <= 0 ) return MB_FAILURE;
1527 
1528  connectivity.resize( num_side_nodes );
1529  for( int k = 0; k < num_side_nodes; ++k )
1530  connectivity[k] = nodes[side_node_idx[k]];
1531 
1532  if( MB_SUCCESS != create_sideset_element( connectivity, subtype, ent_handle, sense ) ) return MB_FAILURE;
1533  if( 1 == sense )
1534  entities_to_add.push_back( ent_handle );
1535  else
1536  reverse_entities.push_back( ent_handle );
1537 
1538  // Read in distribution factor array
1539  if( num_dist_factors )
1540  {
1541  for( int k = 0; k < 2; k++ )
1542  dist_factor_vector.push_back( temp_dist_factor_vector[df_index++] );
1543  }
1544  }
1545  else if( type == MBTRI )
1546  {
1547  int side_offset = 0;
1548  if( number_dimensions() == 3 && side_list[i] <= 2 )
1549  {
1550  if( 1 == sense )
1551  entities_to_add.push_back( ent_handle );
1552  else
1553  reverse_entities.push_back( ent_handle );
1554  if( num_dist_factors )
1555  {
1556  for( int k = 0; k < 3; k++ )
1557  dist_factor_vector.push_back( temp_dist_factor_vector[df_index++] );
1558  }
1559  }
1560  else
1561  {
1562  if( number_dimensions() == 3 )
1563  {
1564  if( side_list[i] > 2 ) side_offset = 2;
1565  }
1566 
1567  // Get the nodes of the element
1568  if( mdbImpl->get_connectivity( ent_handle, nodes, num_elem_nodes ) != MB_SUCCESS ) return MB_FAILURE;
1569 
1570  CN::SubEntityNodeIndices( type, num_elem_nodes, 1, side_num - side_offset, subtype, num_side_nodes,
1571  side_node_idx );
1572  if( num_side_nodes <= 0 ) return MB_FAILURE;
1573 
1574  connectivity.resize( num_side_nodes );
1575  for( int k = 0; k < num_side_nodes; ++k )
1576  connectivity[k] = nodes[side_node_idx[k]];
1577 
1578  if( MB_SUCCESS != create_sideset_element( connectivity, subtype, ent_handle, sense ) )
1579  return MB_FAILURE;
1580  if( 1 == sense )
1581  entities_to_add.push_back( ent_handle );
1582  else
1583  reverse_entities.push_back( ent_handle );
1584 
1585  if( num_dist_factors )
1586  {
1587  for( int k = 0; k < 2; k++ )
1588  dist_factor_vector.push_back( temp_dist_factor_vector[df_index++] );
1589  }
1590  }
1591  }
1592  }
1593 
1594  return MB_SUCCESS;
1595 }
1596 
1597 ErrorCode ReadNCDF::create_sideset_element( const std::vector< EntityHandle >& connectivity,
1598  EntityType type,
1599  EntityHandle& handle,
1600  int& sense )
1601 {
1602  // Get adjacent entities
1603  ErrorCode error = MB_SUCCESS;
1604  int to_dim = CN::Dimension( type );
1605  // for matching purposes , consider only corner nodes
1606  // basically, assume everything is matching if the corner nodes are matching, and
1607  // the number of total nodes is the same
1608  int nb_corner_nodes = CN::VerticesPerEntity( type );
1609  std::vector< EntityHandle > adj_ent;
1610  mdbImpl->get_adjacencies( &( connectivity[0] ), 1, to_dim, false, adj_ent );
1611 
1612  // For each entity, see if we can find a match
1613  // If we find a match, return it
1614  bool match_found = false;
1615  std::vector< EntityHandle > match_conn;
1616  // By default, assume positive sense
1617  sense = 1;
1618  for( unsigned int i = 0; i < adj_ent.size() && match_found == false; i++ )
1619  {
1620  // Get the connectivity
1621  error = mdbImpl->get_connectivity( &( adj_ent[i] ), 1, match_conn );
1622  if( error != MB_SUCCESS ) continue;
1623 
1624  // Make sure they have the same number of vertices (higher order elements ?)
1625  if( match_conn.size() != connectivity.size() ) continue;
1626 
1627  // Find a matching node
1628  std::vector< EntityHandle >::iterator iter;
1629  std::vector< EntityHandle >::iterator end_corner = match_conn.begin() + nb_corner_nodes;
1630  iter = std::find( match_conn.begin(), end_corner, connectivity[0] );
1631  if( iter == end_corner ) continue;
1632 
1633  // Rotate to match connectivity
1634  // rotate only corner nodes, ignore the rest from now on
1635  std::rotate( match_conn.begin(), iter, end_corner );
1636 
1637  bool they_match = true;
1638  for( int j = 1; j < nb_corner_nodes; j++ )
1639  {
1640  if( connectivity[j] != match_conn[j] )
1641  {
1642  they_match = false;
1643  break;
1644  }
1645  }
1646 
1647  // If we didn't get a match
1648  if( !they_match )
1649  {
1650  // Try the opposite sense
1651  they_match = true;
1652 
1653  int k = nb_corner_nodes - 1;
1654  for( int j = 1; j < nb_corner_nodes; )
1655  {
1656  if( connectivity[j] != match_conn[k] )
1657  {
1658  they_match = false;
1659  break;
1660  }
1661  ++j;
1662  --k;
1663  }
1664  // If they matched here, sense is reversed
1665  if( they_match ) sense = -1;
1666  }
1667  match_found = they_match;
1668  if( match_found ) handle = adj_ent[i];
1669  }
1670 
1671  // If we didn't find a match, create an element
1672  if( !match_found ) error = mdbImpl->create_element( type, &connectivity[0], connectivity.size(), handle );
1673 
1674  return error;
1675 }
1676 
1677 ErrorCode ReadNCDF::find_side_element_type( const int exodus_id,
1678  ExoIIElementType& elem_type,
1679  ReadBlockData& block_data,
1680  int& df_index,
1681  int side_id )
1682 {
1683  std::vector< ReadBlockData >::iterator iter, end_iter;
1684  iter = blocksLoading.begin();
1685  end_iter = blocksLoading.end();
1686  elem_type = EXOII_MAX_ELEM_TYPE;
1687 
1688  for( ; iter != end_iter; ++iter )
1689  {
1690  if( exodus_id >= ( *iter ).startExoId && exodus_id < ( *iter ).startExoId + ( *iter ).numElements )
1691  {
1692  elem_type = ( *iter ).elemType;
1693 
1694  // If we're not reading this block in
1695  if( !( *iter ).reading_in )
1696  {
1697  // Offset df_indes according to type
1698  if( elem_type >= EXOII_HEX && elem_type <= EXOII_HEX27 )
1699  df_index += 4;
1700  else if( elem_type >= EXOII_TETRA && elem_type <= EXOII_TETRA14 )
1701  df_index += 3;
1702  else if( elem_type >= EXOII_QUAD && elem_type <= EXOII_QUAD9 )
1703  df_index += 2;
1704  else if( elem_type >= EXOII_SHELL && elem_type <= EXOII_SHELL9 )
1705  {
1706  if( side_id == 1 || side_id == 2 )
1707  df_index += 4;
1708  else
1709  df_index += 2;
1710  }
1711  else if( elem_type >= EXOII_TRI && elem_type <= EXOII_TRI7 )
1712  df_index += 3;
1713 
1714  return MB_FAILURE;
1715  }
1716 
1717  block_data = *iter;
1718 
1719  return MB_SUCCESS;
1720  }
1721  }
1722  return MB_FAILURE;
1723 }
1724 
1725 ErrorCode ReadNCDF::read_qa_records( EntityHandle file_set )
1726 {
1727  std::vector< std::string > qa_records;
1728  read_qa_information( qa_records );
1729 
1730  std::vector< char > tag_data;
1731  for( std::vector< std::string >::iterator i = qa_records.begin(); i != qa_records.end(); ++i )
1732  {
1733  std::copy( i->begin(), i->end(), std::back_inserter( tag_data ) );
1734  tag_data.push_back( '\0' );
1735  }
1736 
1737  // If there were qa_records...tag them to the mCurrentMeshHandle
1738  if( !tag_data.empty() )
1739  {
1740  const void* ptr = &tag_data[0];
1741  int size = tag_data.size();
1742  if( mdbImpl->tag_set_by_ptr( mQaRecordTag, &file_set, 1, &ptr, &size ) != MB_SUCCESS )
1743  {
1744  return MB_FAILURE;
1745  }
1746  }
1747 
1748  return MB_SUCCESS;
1749 }
1750 
1751 ErrorCode ReadNCDF::read_qa_information( std::vector< std::string >& qa_record_list )
1752 {
1753  // Inquire on the genesis file to find the number of qa records
1754 
1755  int number_records = 0;
1756 
1757  int temp_dim;
1758  GET_DIM( temp_dim, "num_qa_rec", number_records );
1759  std::vector< char > data( max_str_length + 1 );
1760 
1761  for( int i = 0; i < number_records; i++ )
1762  {
1763  for( int j = 0; j < 4; j++ )
1764  {
1765  data[max_str_length] = '\0';
1766  if( read_qa_string( &data[0], i, j ) != MB_SUCCESS ) return MB_FAILURE;
1767 
1768  qa_record_list.push_back( &data[0] );
1769  }
1770  }
1771 
1772  return MB_SUCCESS;
1773 }
1774 
1775 ErrorCode ReadNCDF::read_qa_string( char* temp_string, int record_number, int record_position )
1776 {
1777  std::vector< int > dims;
1778  int temp_var = -1;
1779  GET_VAR( "qa_records", temp_var, dims );
1780  if( -1 == temp_var )
1781  {
1782  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting qa record variable" );
1783  return MB_FAILURE;
1784  }
1785  size_t count[3], start[3];
1786  start[0] = record_number;
1787  start[1] = record_position;
1788  start[2] = 0;
1789  count[0] = 1;
1790  count[1] = 1;
1791  count[2] = max_str_length;
1792  int fail = nc_get_vara_text( ncFile, temp_var, start, count, temp_string );
1793  if( NC_NOERR != fail )
1794  {
1795  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem setting current record number variable position" );
1796  }
1797  // Get the variable id in the exodus file
1798 
1799  return MB_SUCCESS;
1800 }
1801 
1802 ErrorCode ReadNCDF::update( const char* exodus_file_name,
1803  const FileOptions& opts,
1804  const int num_blocks,
1805  const int* blocks_to_load,
1806  const EntityHandle cub_file_set )
1807 {
1808  // Function : updating current database from new exodus_file.
1809  // Creator: Jane Hu
1810  // opts is currently designed as following
1811  // tdata = <var_name>[, time][,op][,destination]
1812  // where var_name show the tag name to be updated, this version just takes
1813  // coord.
1814  // time is the optional, and it gives time step of each of the mesh
1815  // info in exodus file. It start from 1.
1816  // op is the operation that is going to be performed on the var_name info.
1817  // currently support 'set'
1818  // destination shows where to store the updated info, currently assume it is
1819  // stored in the same database by replacing the old info if there's no input
1820  // for destination, or the destination data is given in exodus format and just
1821  // need to update the coordinates.
1822  // Assumptions:
1823  // 1. Assume the num_el_blk's in both DB and update exodus file are the same.
1824  // 2. Assume num_el_in_blk1...num_el_in_blk(num_el_blk) numbers are matching, may in
1825  // different order. example: num_el_in_blk11 = num_el_in_blk22 && num_el_in_blk12 =
1826  // num_el_in_blk21.
1827  // 3. In exodus file, get node_num_map
1828  // 4. loop through the node_num_map, use it to find the node in the cub file.
1829  // 5. Replace coord[0][n] with coordx[m] + vals_nod_var1(time_step, m) for all directions for
1830  // matching nodes.
1831  // Test: try to get hexes
1832 
1833  // *******************************************************************
1834  // Move nodes to their deformed locations.
1835  // *******************************************************************
1836  ErrorCode rval;
1837  std::string s;
1838  rval = opts.get_str_option( "tdata", s );
1839  if( MB_SUCCESS != rval )
1840  {
1841  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem reading file options" );
1842  }
1843  std::vector< std::string > tokens;
1844  tokenize( s, tokens, "," );
1845 
1846  // 1. Check for time step to find the match time
1847  int time_step = 1;
1848  if( tokens.size() > 1 && !tokens[1].empty() )
1849  {
1850  const char* time_s = tokens[1].c_str();
1851  char* endptr;
1852  long int pval = strtol( time_s, &endptr, 0 );
1853  std::string end = endptr;
1854  if( !end.empty() ) // Syntax error
1855  return MB_TYPE_OUT_OF_RANGE;
1856 
1857  // Check for overflow (parsing long int, returning int)
1858  time_step = pval;
1859  if( pval != (long int)time_step ) return MB_TYPE_OUT_OF_RANGE;
1860  if( time_step <= 0 ) return MB_TYPE_OUT_OF_RANGE;
1861  }
1862 
1863  // 2. Check for the operations, currently support set.
1864  const char* op;
1865  if( tokens.size() < 3 || ( tokens[2] != "set" && tokens[2] != "add" ) )
1866  {
1867  MB_SET_ERR( MB_TYPE_OUT_OF_RANGE, "ReadNCDF: invalid operation specified for update" );
1868  }
1869  op = tokens[2].c_str();
1870 
1871  // Open netcdf/exodus file
1872  ncFile = 0;
1873  int fail = nc_open( exodus_file_name, 0, &ncFile );
1874  if( !ncFile )
1875  {
1876  MB_SET_ERR( MB_FILE_DOES_NOT_EXIST, "ReadNCDF:: problem opening Netcdf/Exodus II file " << exodus_file_name );
1877  }
1878 
1879  rval = read_exodus_header();MB_CHK_ERR( rval );
1880 
1881  // Check to make sure that the requested time step exists
1882  int ncdim = -1;
1883  int max_time_steps;
1884  GET_DIM( ncdim, "time_step", max_time_steps );
1885  if( -1 == ncdim )
1886  {
1887  std::cout << "ReadNCDF: could not get number of time steps" << std::endl;
1888  return MB_FAILURE;
1889  }
1890  std::cout << " Maximum time step=" << max_time_steps << std::endl;
1891  if( max_time_steps < time_step )
1892  {
1893  std::cout << "ReadNCDF: time step is greater than max_time_steps" << std::endl;
1894  return MB_FAILURE;
1895  }
1896 
1897  // Get the time
1898  std::vector< double > times( max_time_steps );
1899  int nc_var = -1;
1900  GET_1D_DBL_VAR( "time_whole", nc_var, times );
1901  if( -1 == nc_var )
1902  {
1903  std::cout << "ReadNCDF: unable to get time variable" << std::endl;
1904  }
1905  else
1906  {
1907  std::cout << " Step " << time_step << " is at " << times[time_step - 1] << " seconds" << std::endl;
1908  }
1909 
1910  // Read in the node_num_map.
1911  std::vector< int > ptr( numberNodes_loading );
1912 
1913  int varid = -1;
1914  GET_1D_INT_VAR( "node_num_map", varid, ptr );
1915  if( -1 == varid )
1916  {
1917  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting node number map data" );
1918  }
1919 
1920  // Read in the deformations.
1921  std::vector< std::vector< double > > deformed_arrays( 3 );
1922  std::vector< std::vector< double > > orig_coords( 3 );
1923  deformed_arrays[0].reserve( numberNodes_loading );
1924  deformed_arrays[1].reserve( numberNodes_loading );
1925  deformed_arrays[2].reserve( numberNodes_loading );
1926  orig_coords[0].reserve( numberNodes_loading );
1927  orig_coords[1].reserve( numberNodes_loading );
1928  orig_coords[2].reserve( numberNodes_loading );
1929  size_t start[2] = { static_cast< size_t >( time_step - 1 ), 0 },
1930  count[2] = { 1, static_cast< size_t >( numberNodes_loading ) };
1931  std::vector< int > dims;
1932  int coordx = -1, coordy = -1, coordz = -1;
1933  GET_VAR( "vals_nod_var1", coordx, dims );
1934  GET_VAR( "vals_nod_var2", coordy, dims );
1935  if( numberDimensions_loading == 3 ) GET_VAR( "vals_nod_var3", coordz, dims );
1936  if( -1 == coordx || -1 == coordy || ( numberDimensions_loading == 3 && -1 == coordz ) )
1937  {
1938  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting coords variable" );
1939  }
1940 
1941  fail = nc_get_vara_double( ncFile, coordx, start, count, &deformed_arrays[0][0] );
1942  if( NC_NOERR != fail )
1943  {
1944  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting x deformation array" );
1945  }
1946 
1947  fail = nc_get_vara_double( ncFile, coordy, start, count, &deformed_arrays[1][0] );
1948  if( NC_NOERR != fail )
1949  {
1950  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting y deformation array" );
1951  }
1952  if( numberDimensions_loading == 3 )
1953  {
1954  fail = nc_get_vara_double( ncFile, coordz, start, count, &deformed_arrays[2][0] );
1955  if( NC_NOERR != fail )
1956  {
1957  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting z deformation array" );
1958  }
1959  }
1960 
1961  int coord1 = -1, coord2 = -1, coord3 = -1;
1962  GET_1D_DBL_VAR( "coordx", coord1, orig_coords[0] );
1963  if( -1 == coord1 )
1964  {
1965  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting x coord array" );
1966  }
1967  GET_1D_DBL_VAR( "coordy", coord2, orig_coords[1] );
1968  if( -1 == coord2 )
1969  {
1970  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting y coord array" );
1971  }
1972  if( numberDimensions_loading == 3 )
1973  {
1974  GET_1D_DBL_VAR( "coordz", coord3, orig_coords[2] );
1975  if( -1 == coord3 )
1976  {
1977  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting z coord array" );
1978  }
1979  }
1980 
1981  // b. Deal with DB file : get node info. according to node_num_map.
1982  if( tokens[0] != "coord" && tokens[0] != "COORD" ) return MB_NOT_IMPLEMENTED;
1983 
1984  if( strcmp( op, "set" ) != 0 && strcmp( op, " set" ) != 0 ) return MB_NOT_IMPLEMENTED;
1985 
1986  // Two methods of matching nodes (id vs. proximity)
1987  const bool match_node_ids = true;
1988 
1989  // Get nodes in cubit file
1990  Range cub_verts;
1991  rval = mdbImpl->get_entities_by_type( cub_file_set, MBVERTEX, cub_verts );MB_CHK_ERR( rval );
1992  std::cout << " cub_file_set contains " << cub_verts.size() << " nodes." << std::endl;
1993 
1994  // Some accounting
1995  std::cout << " exodus file contains " << numberNodes_loading << " nodes." << std::endl;
1996  double max_magnitude = 0;
1997  double average_magnitude = 0;
1998  int found = 0;
1999  int lost = 0;
2000  std::map< int, EntityHandle > cub_verts_id_map;
2001  AdaptiveKDTree kdtree( mdbImpl );
2002  EntityHandle root;
2003 
2004  // Should not use cub verts unless they have been matched. Place in a map
2005  // for fast handle_by_id lookup.
2006  std::map< int, EntityHandle > matched_cub_vert_id_map;
2007 
2008  // Place cub verts in a map for searching by id
2009  if( match_node_ids )
2010  {
2011  std::vector< int > cub_ids( cub_verts.size() );
2012  rval = mdbImpl->tag_get_data( mGlobalIdTag, cub_verts, &cub_ids[0] );MB_CHK_ERR( rval );
2013  for( unsigned i = 0; i != cub_verts.size(); ++i )
2014  {
2015  cub_verts_id_map.insert( std::pair< int, EntityHandle >( cub_ids[i], cub_verts[i] ) );
2016  }
2017 
2018  // Place cub verts in a kdtree for searching by proximity
2019  }
2020  else
2021  {
2022  FileOptions tree_opts( "MAX_PER_LEAF=1;SPLITS_PER_DIR=1;CANDIDATE_PLANE_SET=0" );
2023  rval = kdtree.build_tree( cub_verts, &root, &tree_opts );MB_CHK_ERR( rval );
2024  AdaptiveKDTreeIter tree_iter;
2025  rval = kdtree.get_tree_iterator( root, tree_iter );MB_CHK_ERR( rval );
2026  }
2027 
2028  // For each exo vert, find the matching cub vert
2029  for( int i = 0; i < numberNodes_loading; ++i )
2030  {
2031  int exo_id = ptr[i];
2032  CartVect exo_coords( orig_coords[0][i], orig_coords[1][i], orig_coords[2][i] );
2033  EntityHandle cub_vert = 0;
2034  bool found_match = false;
2035 
2036  // By id
2037  if( match_node_ids )
2038  {
2039  std::map< int, EntityHandle >::iterator i_iter;
2040  i_iter = cub_verts_id_map.find( exo_id );
2041  if( i_iter != cub_verts_id_map.end() )
2042  {
2043  found_match = true;
2044  cub_vert = i_iter->second;
2045  }
2046 
2047  // By proximity
2048  }
2049  else
2050  {
2051  // The MAX_NODE_DIST is the farthest distance to search for a node.
2052  // For the 1/12th symmetry 85 pin model, the max node dist could not be less
2053  // than 1e-1 (March 26, 2010).
2054  const double MAX_NODE_DIST = 1e-1;
2055 
2056  std::vector< EntityHandle > leaves;
2057  double min_dist = MAX_NODE_DIST;
2058  rval = kdtree.distance_search( exo_coords.array(), MAX_NODE_DIST, leaves );MB_CHK_ERR( rval );
2059  for( std::vector< EntityHandle >::const_iterator j = leaves.begin(); j != leaves.end(); ++j )
2060  {
2061  std::vector< EntityHandle > leaf_verts;
2062  rval = mdbImpl->get_entities_by_type( *j, MBVERTEX, leaf_verts );MB_CHK_ERR( rval );
2063  for( std::vector< EntityHandle >::const_iterator k = leaf_verts.begin(); k != leaf_verts.end(); ++k )
2064  {
2065  CartVect orig_cub_coords, difference;
2066  rval = mdbImpl->get_coords( &( *k ), 1, orig_cub_coords.array() );MB_CHK_ERR( rval );
2067  difference = orig_cub_coords - exo_coords;
2068  double dist = difference.length();
2069  if( dist < min_dist )
2070  {
2071  min_dist = dist;
2072  cub_vert = *k;
2073  }
2074  }
2075  }
2076  if( 0 != cub_vert ) found_match = true;
2077  }
2078 
2079  // If a match is found, update it with the deformed coords from the exo file.
2080  if( found_match )
2081  {
2082  CartVect updated_exo_coords;
2083  matched_cub_vert_id_map.insert( std::pair< int, EntityHandle >( exo_id, cub_vert ) );
2084  updated_exo_coords[0] = orig_coords[0][i] + deformed_arrays[0][i];
2085  updated_exo_coords[1] = orig_coords[1][i] + deformed_arrays[1][i];
2086 
2087  if( numberDimensions_loading == 3 ) updated_exo_coords[2] = orig_coords[2][i] + deformed_arrays[2][i];
2088  rval = mdbImpl->set_coords( &cub_vert, 1, updated_exo_coords.array() );MB_CHK_ERR( rval );
2089  ++found;
2090 
2091  double magnitude =
2092  sqrt( deformed_arrays[0][i] * deformed_arrays[0][i] + deformed_arrays[1][i] * deformed_arrays[1][i] +
2093  deformed_arrays[2][i] * deformed_arrays[2][i] );
2094  if( magnitude > max_magnitude ) max_magnitude = magnitude;
2095  average_magnitude += magnitude;
2096  }
2097  else
2098  {
2099  ++lost;
2100  std::cout << "cannot match exo node " << exo_id << " " << exo_coords << std::endl;
2101  }
2102  }
2103 
2104  // Exo verts that could not be matched have already been printed. Now print the
2105  // cub verts that could not be matched.
2106  if( matched_cub_vert_id_map.size() < cub_verts.size() )
2107  {
2108  Range unmatched_cub_verts = cub_verts;
2109  for( std::map< int, EntityHandle >::const_iterator i = matched_cub_vert_id_map.begin();
2110  i != matched_cub_vert_id_map.end(); ++i )
2111  {
2112  unmatched_cub_verts.erase( i->second );
2113  }
2114 
2115  for( Range::const_iterator i = unmatched_cub_verts.begin(); i != unmatched_cub_verts.end(); ++i )
2116  {
2117  int cub_id;
2118  rval = mdbImpl->tag_get_data( mGlobalIdTag, &( *i ), 1, &cub_id );MB_CHK_ERR( rval );
2119 
2120  CartVect cub_coords;
2121  rval = mdbImpl->get_coords( &( *i ), 1, cub_coords.array() );MB_CHK_ERR( rval );
2122  std::cout << "cannot match cub node " << cub_id << " " << cub_coords << std::endl;
2123  }
2124 
2125  std::cout << " " << unmatched_cub_verts.size() << " nodes from the cub file could not be matched."
2126  << std::endl;
2127  }
2128 
2129  // Summarize statistics
2130  std::cout << " " << found << " nodes from the exodus file were matched in the cub_file_set ";
2131  if( match_node_ids )
2132  {
2133  std::cout << "by id." << std::endl;
2134  }
2135  else
2136  {
2137  std::cout << "by proximity." << std::endl;
2138  }
2139 
2140  // Fail if all of the nodes could not be matched.
2141  if( 0 != lost )
2142  {
2143  std::cout << "Error: " << lost << " nodes from the exodus file could not be matched." << std::endl;
2144  // return MB_FAILURE;
2145  }
2146  std::cout << " maximum node displacement magnitude: " << max_magnitude << " cm" << std::endl;
2147  std::cout << " average node displacement magnitude: " << average_magnitude / found << " cm" << std::endl;
2148 
2149  // *******************************************************************
2150  // Remove dead elements from the MOAB instance.
2151  // *******************************************************************
2152 
2153  // How many element variables are in the file?
2154  int n_elem_var;
2155  GET_DIM( ncdim, "num_elem_var", n_elem_var );
2156 
2157  // Get element variable names
2158  varid = -1;
2159  int cstatus = nc_inq_varid( ncFile, "name_elem_var", &varid );
2160  std::vector< char > names_memory( n_elem_var * max_str_length );
2161  std::vector< char* > names( n_elem_var );
2162  for( int i = 0; i < n_elem_var; ++i )
2163  names[i] = &names_memory[i * max_str_length];
2164  if( cstatus != NC_NOERR || varid == -1 )
2165  {
2166  std::cout << "ReadNCDF: name_elem_var does not exist" << std::endl;
2167  return MB_FAILURE;
2168  }
2169  int status = nc_get_var_text( ncFile, varid, &names_memory[0] );
2170  if( NC_NOERR != status )
2171  {
2172  MB_SET_ERR( MB_FAILURE, "ReadNCDF: Problem getting element variable names" );
2173  }
2174 
2175  // Is one of the element variable names "death_status"? If so, get its index
2176  // in the element variable array.
2177  int death_index;
2178  bool found_death_index = false;
2179  for( int i = 0; i < n_elem_var; ++i )
2180  {
2181  std::string temp( names[i] );
2182  std::string::size_type pos0 = temp.find( "death_status" );
2183  std::string::size_type pos1 = temp.find( "Death_Status" );
2184  std::string::size_type pos2 = temp.find( "DEATH_STATUS" );
2185  if( std::string::npos != pos0 || std::string::npos != pos1 || std::string::npos != pos2 )
2186  {
2187  found_death_index = true;
2188  death_index = i + 1; // NetCDF variables start with 1
2189  break;
2190  }
2191  }
2192  if( !found_death_index )
2193  {
2194  MB_SET_ERR( MB_FAILURE, "ReadNCDF: Problem getting index of death_status variable" );
2195  }
2196 
2197  // The exodus header has already been read. This contains the number of element
2198  // blocks.
2199 
2200  // Dead elements are listed by block. Read the block headers to determine how
2201  // many elements are in each block.
2202  rval = read_block_headers( blocks_to_load, num_blocks );
2203  if( MB_FAILURE == rval )
2204  {
2205  MB_SET_ERR( MB_FAILURE, "ReadNCDF: Problem reading block headers" );
2206  }
2207 
2208  // Dead elements from the Exodus file can be located in the cub_file_set by id
2209  // or by connectivity. Currently, finding elements by id requires careful book
2210  // keeping when constructing the model in Cubit. To avoid this, one can match
2211  // elements by connectivity instead.
2212  const bool match_elems_by_connectivity = true;
2213 
2214  // Get the element id map. The ids in the map are from the elements in the blocks.
2215  // elem_num_map(blk1 elem ids, blk2 elem ids, blk3 elem ids, ...)
2216  std::vector< int > elem_ids( numberNodes_loading );
2217  if( !match_elems_by_connectivity )
2218  {
2219  GET_1D_INT_VAR( "elem_num_map", varid, elem_ids );
2220  if( -1 == varid )
2221  {
2222  MB_SET_ERR( MB_FAILURE, "ReadNCDF: Problem getting element number map data" );
2223  }
2224  }
2225 
2226  // For each block
2227  int first_elem_id_in_block = 0;
2228  int block_count = 1; // NetCDF variables start with 1
2229  int total_elems = 0;
2230  int total_dead_elems = 0;
2231  for( std::vector< ReadBlockData >::iterator i = blocksLoading.begin(); i != blocksLoading.end(); ++i )
2232  {
2233  // Get the ncdf connect variable
2234  std::string temp_string( "connect" );
2235  std::stringstream temp_ss;
2236  temp_ss << block_count;
2237  temp_string += temp_ss.str();
2238  temp_string += "\0";
2239  std::vector< int > ldims;
2240  GET_VAR( temp_string.c_str(), nc_var, ldims );
2241  if( !nc_var )
2242  {
2243  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting connect variable" );
2244  }
2245  // The element type is an attribute of the connectivity variable
2246  nc_type att_type;
2247  size_t att_len;
2248  fail = nc_inq_att( ncFile, nc_var, "elem_type", &att_type, &att_len );
2249  if( NC_NOERR != fail )
2250  {
2251  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting elem type attribute" );
2252  }
2253  std::vector< char > dum_str( att_len + 1 );
2254  fail = nc_get_att_text( ncFile, nc_var, "elem_type", &dum_str[0] );
2255  if( NC_NOERR != fail )
2256  {
2257  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting elem type" );
2258  }
2259  ExoIIElementType elem_type = ExoIIUtil::static_element_name_to_type( &dum_str[0] );
2260  ( *i ).elemType = elem_type;
2261  const EntityType mb_type = ExoIIUtil::ExoIIElementMBEntity[( *i ).elemType];
2262 
2263  // Get the number of nodes per element
2264  unsigned int nodes_per_element = ExoIIUtil::VerticesPerElement[( *i ).elemType];
2265 
2266  // Read the connectivity into that memory.
2267  // int exo_conn[i->numElements][nodes_per_element];
2268  int* exo_conn = new int[i->numElements * nodes_per_element];
2269  // NcBool status = temp_var->get(&exo_conn[0][0], i->numElements, nodes_per_element);
2270  size_t lstart[2] = { 0, 0 }, lcount[2];
2271  lcount[0] = i->numElements;
2272  lcount[1] = nodes_per_element;
2273  fail = nc_get_vara_int( ncFile, nc_var, lstart, lcount, exo_conn );
2274  if( NC_NOERR != fail )
2275  {
2276  delete[] exo_conn;
2277  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting connectivity" );
2278  }
2279 
2280  // Get the death_status at the correct time step.
2281  std::vector< double > death_status( i->numElements ); // It seems wrong, but it uses doubles
2282  std::string array_name( "vals_elem_var" );
2283  temp_ss.str( "" ); // stringstream won't clear by temp.clear()
2284  temp_ss << death_index;
2285  array_name += temp_ss.str();
2286  array_name += "eb";
2287  temp_ss.str( "" ); // stringstream won't clear by temp.clear()
2288  temp_ss << block_count;
2289  array_name += temp_ss.str();
2290  array_name += "\0";
2291  GET_VAR( array_name.c_str(), nc_var, ldims );
2292  if( !nc_var )
2293  {
2294  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting death status variable" );
2295  }
2296  lstart[0] = time_step - 1;
2297  lstart[1] = 0;
2298  lcount[0] = 1;
2299  lcount[1] = i->numElements;
2300  status = nc_get_vara_double( ncFile, nc_var, lstart, lcount, &death_status[0] );
2301  if( NC_NOERR != status )
2302  {
2303  delete[] exo_conn;
2304  MB_SET_ERR( MB_FAILURE, "ReadNCDF: Problem setting time step for death_status" );
2305  }
2306 
2307  // Look for dead elements. If there is too many dead elements and this starts
2308  // to take too long, I should put the elems in a kd-tree for more efficient
2309  // searching. Alternatively I could get the exo connectivity and match nodes.
2310  int dead_elem_counter = 0, missing_elem_counter = 0;
2311  for( int j = 0; j < i->numElements; ++j )
2312  {
2313  if( 1 != death_status[j] )
2314  {
2315  Range cub_elem, cub_nodes;
2316  if( match_elems_by_connectivity )
2317  {
2318  // Get exodus nodes for the element
2319  std::vector< int > elem_conn( nodes_per_element );
2320  for( unsigned int k = 0; k < nodes_per_element; ++k )
2321  {
2322  // elem_conn[k] = exo_conn[j][k];
2323  elem_conn[k] = exo_conn[j * nodes_per_element + k];
2324  }
2325  // Get the ids of the nodes (assume we are matching by id)
2326  // Remember that the exodus array locations start with 1 (not 0).
2327  std::vector< int > elem_conn_node_ids( nodes_per_element );
2328  for( unsigned int k = 0; k < nodes_per_element; ++k )
2329  {
2330  elem_conn_node_ids[k] = ptr[elem_conn[k] - 1];
2331  }
2332  // Get the cub_file_set nodes by id
2333  // The map is a log search and takes almost no time.
2334  // MOAB's linear tag search takes 5-10 minutes.
2335  for( unsigned int k = 0; k < nodes_per_element; ++k )
2336  {
2337  std::map< int, EntityHandle >::iterator k_iter;
2338  k_iter = matched_cub_vert_id_map.find( elem_conn_node_ids[k] );
2339 
2340  if( k_iter == matched_cub_vert_id_map.end() )
2341  {
2342  std::cout << "ReadNCDF: Found no cub node with id=" << elem_conn_node_ids[k]
2343  << ", but expected to find only 1." << std::endl;
2344  break;
2345  }
2346  cub_nodes.insert( k_iter->second );
2347  }
2348 
2349  if( nodes_per_element != cub_nodes.size() )
2350  {
2351  std::cout << "ReadNCDF: nodes_per_elemenet != cub_nodes.size()" << std::endl;
2352  delete[] exo_conn;
2353  return MB_INVALID_SIZE;
2354  }
2355 
2356  // Get the cub_file_set element with the same nodes
2357  int to_dim = CN::Dimension( mb_type );
2358  rval = mdbImpl->get_adjacencies( cub_nodes, to_dim, false, cub_elem );
2359  if( MB_SUCCESS != rval )
2360  {
2361  std::cout << "ReadNCDF: could not get dead cub element" << std::endl;
2362  delete[] exo_conn;
2363  return rval;
2364  }
2365 
2366  // Pronto/Presto renumbers elements, so matching cub and exo elements by
2367  // id is not possible at this time.
2368  }
2369  else
2370  {
2371  // Get dead element's id
2372  int elem_id = elem_ids[first_elem_id_in_block + j];
2373  void* id[] = { &elem_id };
2374  // Get the element by id
2375  rval = mdbImpl->get_entities_by_type_and_tag( cub_file_set, mb_type, &mGlobalIdTag, id, 1, cub_elem,
2376  Interface::INTERSECT );
2377  if( MB_SUCCESS != rval )
2378  {
2379  delete[] exo_conn;
2380  return rval;
2381  }
2382  }
2383 
2384  if( 1 == cub_elem.size() )
2385  {
2386  // Delete the dead element from the cub file. It will be removed from sets
2387  // ONLY if they are tracking meshsets.
2388  rval = mdbImpl->remove_entities( cub_file_set, cub_elem );
2389  if( MB_SUCCESS != rval )
2390  {
2391  delete[] exo_conn;
2392  return rval;
2393  }
2394  rval = mdbImpl->delete_entities( cub_elem );
2395  if( MB_SUCCESS != rval )
2396  {
2397  delete[] exo_conn;
2398  return rval;
2399  }
2400  }
2401  else
2402  {
2403  std::cout << "ReadNCDF: Should have found 1 element with type=" << mb_type
2404  << " in cub_file_set, but instead found " << cub_elem.size() << std::endl;
2405  rval = mdbImpl->list_entities( cub_nodes );
2406  ++missing_elem_counter;
2407  delete[] exo_conn;
2408  return MB_FAILURE;
2409  }
2410  ++dead_elem_counter;
2411  }
2412  }
2413  // Print some statistics
2414  temp_ss.str( "" );
2415  temp_ss << i->blockId;
2416  total_dead_elems += dead_elem_counter;
2417  total_elems += i->numElements;
2418  std::cout << " Block " << temp_ss.str() << " has " << dead_elem_counter << "/" << i->numElements
2419  << " dead elements." << std::endl;
2420  if( 0 != missing_elem_counter )
2421  {
2422  std::cout << " " << missing_elem_counter
2423  << " dead elements in this block were not found in the cub_file_set." << std::endl;
2424  }
2425 
2426  // Advance the pointers into element ids and block_count. Memory cleanup.
2427  first_elem_id_in_block += i->numElements;
2428  ++block_count;
2429  delete[] exo_conn;
2430  }
2431 
2432  std::cout << " Total: " << total_dead_elems << "/" << total_elems << " dead elements." << std::endl;
2433 
2434  // Close the file
2435  fail = nc_close( ncFile );
2436  if( NC_NOERR != fail )
2437  {
2438  MB_SET_ERR( MB_FAILURE, "Trouble closing file" );
2439  }
2440 
2441  ncFile = 0;
2442  return MB_SUCCESS;
2443 }
2444 
2445 void ReadNCDF::tokenize( const std::string& str, std::vector< std::string >& tokens, const char* delimiters )
2446 {
2447  std::string::size_type last = str.find_first_not_of( delimiters, 0 );
2448  std::string::size_type pos = str.find_first_of( delimiters, last );
2449  while( std::string::npos != pos && std::string::npos != last )
2450  {
2451  tokens.push_back( str.substr( last, pos - last ) );
2452  last = str.find_first_not_of( delimiters, pos );
2453  pos = str.find_first_of( delimiters, last );
2454  if( std::string::npos == pos ) pos = str.size();
2455  }
2456 }
2457 
2458 } // namespace moab