Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
WriteSLAC.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 #ifdef _DEBUG
18 // turn off warnings that say they debugging identifier has been truncated
19 // this warning comes up when using some STL containers
20 #pragma warning( disable : 4786 )
21 #endif
22 #endif
23 
24 #include "WriteSLAC.hpp"
25 
26 #include <utility>
27 #include <algorithm>
28 #include <ctime>
29 #include <string>
30 #include <vector>
31 #include <cstdio>
32 #include <cstring>
33 #include <iostream>
34 #include <cassert>
35 
36 #include "netcdf.h"
37 #include "moab/Interface.hpp"
38 #include "moab/Range.hpp"
39 #include "moab/CN.hpp"
40 #include "Internals.hpp"
41 #include "ExoIIUtil.hpp"
42 #include "MBTagConventions.hpp"
43 #include "moab/WriteUtilIface.hpp"
44 
45 #ifndef MOAB_HAVE_NETCDF
46 #error Attempt to compile WriteSLAC with NetCDF disabled.
47 #endif
48 
49 namespace moab
50 {
51 
52 #define GET_VAR( name, id, dims ) \
53  { \
54  ( id ) = -1; \
55  int gvfail = nc_inq_varid( ncFile, name, &( id ) ); \
56  if( NC_NOERR == gvfail ) \
57  { \
58  int ndims; \
59  gvfail = nc_inq_varndims( ncFile, id, &ndims ); \
60  if( NC_NOERR == gvfail ) \
61  { \
62  ( dims ).resize( ndims ); \
63  gvfail = nc_inq_vardimid( ncFile, id, &( dims )[0] ); \
64  } \
65  } \
66  }
67 
69 {
70  return new WriteSLAC( iface );
71 }
72 
73 WriteSLAC::WriteSLAC( Interface* impl ) : mbImpl( impl ), ncFile( 0 )
74 {
75  assert( impl != NULL );
76 
78 
79  // Initialize in case tag_get_handle fails below
80  //! get and cache predefined tag handles
81  int negone = -1;
83  &negone );
84 
86  &negone );
87 
89  &negone );
90 
91  mGlobalIdTag = impl->globalId_tag();
92 
93  int dum_val = -1;
94  impl->tag_get_handle( "__matSetIdTag", 1, MB_TYPE_INTEGER, mMatSetIdTag, MB_TAG_DENSE | MB_TAG_CREAT, &dum_val );
95 
96  impl->tag_get_handle( "WriteSLAC element mark", 1, MB_TYPE_BIT, mEntityMark, MB_TAG_CREAT );
97 }
98 
100 {
103 }
104 
105 void WriteSLAC::reset_matset( std::vector< WriteSLAC::MaterialSetData >& matset_info )
106 {
107  std::vector< WriteSLAC::MaterialSetData >::iterator iter;
108 
109  for( iter = matset_info.begin(); iter != matset_info.end(); ++iter )
110  delete( *iter ).elements;
111 }
112 
113 ErrorCode WriteSLAC::write_file( const char* file_name,
114  const bool overwrite,
115  const FileOptions&,
116  const EntityHandle* ent_handles,
117  const int num_sets,
118  const std::vector< std::string >& /* qa_list */,
119  const Tag* /* tag_list */,
120  int /* num_tags */,
121  int /* export_dimension */ )
122 {
123  assert( 0 != mMaterialSetTag && 0 != mNeumannSetTag && 0 != mDirichletSetTag );
124 
125  // Check the file name
126  if( NULL == strstr( file_name, ".ncdf" ) ) return MB_FAILURE;
127 
128  std::vector< EntityHandle > matsets, dirsets, neusets, entities;
129 
130  fileName = file_name;
131 
132  // Separate into material sets, dirichlet sets, neumann sets
133 
134  if( num_sets == 0 )
135  {
136  // Default to all defined sets
137  Range this_range;
138  mbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, &mMaterialSetTag, NULL, 1, this_range );
139  std::copy( this_range.begin(), this_range.end(), std::back_inserter( matsets ) );
140  this_range.clear();
141  mbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, &mDirichletSetTag, NULL, 1, this_range );
142  std::copy( this_range.begin(), this_range.end(), std::back_inserter( dirsets ) );
143  this_range.clear();
144  mbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, &mNeumannSetTag, NULL, 1, this_range );
145  std::copy( this_range.begin(), this_range.end(), std::back_inserter( neusets ) );
146  }
147  else
148  {
149  int dummy;
150  for( const EntityHandle* iter = ent_handles; iter < ent_handles + num_sets; ++iter )
151  {
152  if( MB_SUCCESS == mbImpl->tag_get_data( mMaterialSetTag, &( *iter ), 1, &dummy ) )
153  matsets.push_back( *iter );
154  else if( MB_SUCCESS == mbImpl->tag_get_data( mDirichletSetTag, &( *iter ), 1, &dummy ) )
155  dirsets.push_back( *iter );
156  else if( MB_SUCCESS == mbImpl->tag_get_data( mNeumannSetTag, &( *iter ), 1, &dummy ) )
157  neusets.push_back( *iter );
158  }
159  }
160 
161  // If there is nothing to write just return.
162  if( matsets.empty() && dirsets.empty() && neusets.empty() ) return MB_FILE_WRITE_ERROR;
163 
164  std::vector< WriteSLAC::MaterialSetData > matset_info;
165  std::vector< WriteSLAC::DirichletSetData > dirset_info;
166  std::vector< WriteSLAC::NeumannSetData > neuset_info;
167 
168  MeshInfo mesh_info;
169 
170  matset_info.clear();
171  if( gather_mesh_information( mesh_info, matset_info, neuset_info, dirset_info, matsets, neusets, dirsets ) !=
172  MB_SUCCESS )
173  {
174  reset_matset( matset_info );
175  return MB_FAILURE;
176  }
177 
178  // Try to open the file after gather mesh info succeeds
179  int fail = nc_create( file_name, overwrite ? NC_CLOBBER : NC_NOCLOBBER, &ncFile );
180  if( NC_NOERR != fail )
181  {
182  reset_matset( matset_info );
183  return MB_FAILURE;
184  }
185 
186  if( initialize_file( mesh_info ) != MB_SUCCESS )
187  {
188  reset_matset( matset_info );
189  return MB_FAILURE;
190  }
191 
192  if( write_nodes( mesh_info.num_nodes, mesh_info.nodes, mesh_info.num_dim ) != MB_SUCCESS )
193  {
194  reset_matset( matset_info );
195  return MB_FAILURE;
196  }
197 
198  if( write_matsets( mesh_info, matset_info, neuset_info ) )
199  {
200  reset_matset( matset_info );
201  return MB_FAILURE;
202  }
203 
204  fail = nc_close( ncFile );
205  if( NC_NOERR != fail ) return MB_FAILURE;
206 
207  return MB_SUCCESS;
208 }
209 
211  std::vector< WriteSLAC::MaterialSetData >& matset_info,
212  std::vector< WriteSLAC::NeumannSetData >& neuset_info,
213  std::vector< WriteSLAC::DirichletSetData >& dirset_info,
214  std::vector< EntityHandle >& matsets,
215  std::vector< EntityHandle >& neusets,
216  std::vector< EntityHandle >& dirsets )
217 {
218  std::vector< EntityHandle >::iterator vector_iter, end_vector_iter;
219 
220  mesh_info.num_nodes = 0;
221  mesh_info.num_elements = 0;
222  mesh_info.num_matsets = 0;
223 
224  int id = 0;
225 
226  vector_iter = matsets.begin();
227  end_vector_iter = matsets.end();
228 
229  mesh_info.num_matsets = matsets.size();
230 
231  std::vector< EntityHandle > parent_meshsets;
232 
233  // Clean out the bits for the element mark
235  mbImpl->tag_get_handle( "WriteSLAC element mark", 1, MB_TYPE_BIT, mEntityMark, MB_TAG_CREAT );
236 
237  int highest_dimension_of_element_matsets = 0;
238 
239  for( vector_iter = matsets.begin(); vector_iter != matsets.end(); ++vector_iter )
240  {
241  WriteSLAC::MaterialSetData matset_data;
242  matset_data.elements = new Range;
243 
244  // For the purpose of qa records, get the parents of these matsets
245  if( mbImpl->get_parent_meshsets( *vector_iter, parent_meshsets ) != MB_SUCCESS ) return MB_FAILURE;
246 
247  // Get all Entity Handles in the mesh set
248  Range dummy_range;
249  mbImpl->get_entities_by_handle( *vector_iter, dummy_range, true );
250 
251  // Wait a minute, we are doing some filtering here that doesn't make sense at this level CJS
252 
253  // Find the dimension of the last entity in this range
254  Range::iterator entity_iter = dummy_range.end();
255  entity_iter = dummy_range.end();
256  --entity_iter;
257  int this_dim = CN::Dimension( TYPE_FROM_HANDLE( *entity_iter ) );
258  entity_iter = dummy_range.begin();
259  while( entity_iter != dummy_range.end() && CN::Dimension( TYPE_FROM_HANDLE( *entity_iter ) ) != this_dim )
260  ++entity_iter;
261 
262  if( entity_iter != dummy_range.end() )
263  std::copy( entity_iter, dummy_range.end(), range_inserter( *( matset_data.elements ) ) );
264 
265  assert( matset_data.elements->begin() == matset_data.elements->end() ||
266  CN::Dimension( TYPE_FROM_HANDLE( *( matset_data.elements->begin() ) ) ) == this_dim );
267 
268  // Get the matset's id
269  if( mbImpl->tag_get_data( mMaterialSetTag, &( *vector_iter ), 1, &id ) != MB_SUCCESS )
270  {
271  MB_SET_ERR( MB_FAILURE, "Couldn't get matset id from a tag for an element matset" );
272  }
273 
274  matset_data.id = id;
275  matset_data.number_attributes = 0;
276 
277  // Iterate through all the elements in the meshset
278  Range::iterator elem_range_iter, end_elem_range_iter;
279  elem_range_iter = matset_data.elements->begin();
280  end_elem_range_iter = matset_data.elements->end();
281 
282  // Get the entity type for this matset, verifying that it's the same for all elements
283  // THIS ASSUMES HANDLES SORT BY TYPE!!!
284  EntityType entity_type = TYPE_FROM_HANDLE( *elem_range_iter );
285  --end_elem_range_iter;
286  if( entity_type != TYPE_FROM_HANDLE( *( end_elem_range_iter++ ) ) )
287  {
288  MB_SET_ERR( MB_FAILURE, "Entities in matset " << id << " not of common type" );
289  }
290 
291  int dimension = -1;
292  if( entity_type == MBQUAD || entity_type == MBTRI )
293  dimension = 3; // Output shells by default
294  else if( entity_type == MBEDGE )
295  dimension = 2;
296  else
297  dimension = CN::Dimension( entity_type );
298 
299  if( dimension > highest_dimension_of_element_matsets ) highest_dimension_of_element_matsets = dimension;
300 
301  matset_data.moab_type = mbImpl->type_from_handle( *( matset_data.elements->begin() ) );
302  if( MBMAXTYPE == matset_data.moab_type ) return MB_FAILURE;
303 
304  std::vector< EntityHandle > tmp_conn;
305  mbImpl->get_connectivity( &( *( matset_data.elements->begin() ) ), 1, tmp_conn );
306  matset_data.element_type =
307  ExoIIUtil::get_element_type_from_num_verts( tmp_conn.size(), entity_type, dimension );
308 
309  if( matset_data.element_type == EXOII_MAX_ELEM_TYPE )
310  {
311  MB_SET_ERR( MB_FAILURE, "Element type in matset " << id << " didn't get set correctly" );
312  }
313 
315 
316  // Number of nodes for this matset
317  matset_data.number_elements = matset_data.elements->size();
318 
319  // Total number of elements
320  mesh_info.num_elements += matset_data.number_elements;
321 
322  // Get the nodes for the elements
323  mWriteIface->gather_nodes_from_elements( *matset_data.elements, mEntityMark, mesh_info.nodes );
324 
325  if( !neusets.empty() )
326  {
327  // If there are neusets, keep track of which elements are being written out
328  for( Range::iterator iter = matset_data.elements->begin(); iter != matset_data.elements->end(); ++iter )
329  {
330  unsigned char bit = 0x1;
331  mbImpl->tag_set_data( mEntityMark, &( *iter ), 1, &bit );
332  }
333  }
334 
335  matset_info.push_back( matset_data );
336  }
337 
338  // If user hasn't entered dimension, we figure it out
339  if( mesh_info.num_dim == 0 )
340  {
341  // Never want 1 or zero dimensions
342  if( highest_dimension_of_element_matsets < 2 )
343  mesh_info.num_dim = 3;
344  else
345  mesh_info.num_dim = highest_dimension_of_element_matsets;
346  }
347 
348  Range::iterator range_iter, end_range_iter;
349  range_iter = mesh_info.nodes.begin();
350  end_range_iter = mesh_info.nodes.end();
351 
352  mesh_info.num_nodes = mesh_info.nodes.size();
353 
354  //------dirsets--------
355 
356  vector_iter = dirsets.begin();
357  end_vector_iter = dirsets.end();
358 
359  for( ; vector_iter != end_vector_iter; ++vector_iter )
360  {
361  WriteSLAC::DirichletSetData dirset_data;
362  dirset_data.id = 0;
363  dirset_data.number_nodes = 0;
364 
365  // Get the dirset's id
366  if( mbImpl->tag_get_data( mDirichletSetTag, &( *vector_iter ), 1, &id ) != MB_SUCCESS )
367  {
368  MB_SET_ERR( MB_FAILURE, "Couldn't get id tag for dirset " << id );
369  }
370 
371  dirset_data.id = id;
372 
373  std::vector< EntityHandle > node_vector;
374  // Get the nodes of the dirset that are in mesh_info.nodes
375  if( mbImpl->get_entities_by_handle( *vector_iter, node_vector, true ) != MB_SUCCESS )
376  {
377  MB_SET_ERR( MB_FAILURE, "Couldn't get nodes in dirset " << id );
378  }
379 
380  std::vector< EntityHandle >::iterator iter, end_iter;
381  iter = node_vector.begin();
382  end_iter = node_vector.end();
383 
384  unsigned char node_marked = 0;
385  ErrorCode result;
386  for( ; iter != end_iter; ++iter )
387  {
388  if( TYPE_FROM_HANDLE( *iter ) != MBVERTEX ) continue;
389  result = mbImpl->tag_get_data( mEntityMark, &( *iter ), 1, &node_marked );MB_CHK_SET_ERR( result, "Couldn't get mark data" );
390 
391  if( 0x1 == node_marked ) dirset_data.nodes.push_back( *iter );
392  }
393 
394  dirset_data.number_nodes = dirset_data.nodes.size();
395  dirset_info.push_back( dirset_data );
396  }
397 
398  //------neusets--------
399  vector_iter = neusets.begin();
400  end_vector_iter = neusets.end();
401 
402  for( ; vector_iter != end_vector_iter; ++vector_iter )
403  {
404  WriteSLAC::NeumannSetData neuset_data;
405 
406  // Get the neuset's id
407  if( mbImpl->tag_get_data( mNeumannSetTag, &( *vector_iter ), 1, &id ) != MB_SUCCESS ) return MB_FAILURE;
408 
409  neuset_data.id = id;
410  neuset_data.mesh_set_handle = *vector_iter;
411 
412  // Get the sides in two lists, one forward the other reverse; starts with forward sense
413  // by convention
414  Range forward_elems, reverse_elems;
415  if( get_neuset_elems( *vector_iter, 0, forward_elems, reverse_elems ) == MB_FAILURE ) return MB_FAILURE;
416 
417  ErrorCode result = get_valid_sides( forward_elems, 1, neuset_data );MB_CHK_SET_ERR( result, "Couldn't get valid sides data" );
418  result = get_valid_sides( reverse_elems, -1, neuset_data );MB_CHK_SET_ERR( result, "Couldn't get valid sides data" );
419 
420  neuset_data.number_elements = neuset_data.elements.size();
421  neuset_info.push_back( neuset_data );
422  }
423 
424  // Get information about interior/exterior tets/hexes, and mark matset ids
425  return gather_interior_exterior( mesh_info, matset_info, neuset_info );
426 }
427 
428 ErrorCode WriteSLAC::get_valid_sides( Range& elems, const int sense, WriteSLAC::NeumannSetData& neuset_data )
429 {
430  // This is where we see if underlying element of side set element is included in output
431 
432  unsigned char element_marked = 0;
433  ErrorCode result;
434  for( Range::iterator iter = elems.begin(); iter != elems.end(); ++iter )
435  {
436  // Should insert here if "side" is a quad/tri on a quad/tri mesh
437  result = mbImpl->tag_get_data( mEntityMark, &( *iter ), 1, &element_marked );MB_CHK_SET_ERR( result, "Couldn't get mark data" );
438 
439  if( 0x1 == element_marked )
440  {
441  neuset_data.elements.push_back( *iter );
442 
443  // TJT TODO: the sense should really be # edges + 1or2
444  neuset_data.side_numbers.push_back( ( sense == 1 ? 1 : 2 ) );
445  }
446  else
447  { // Then "side" is probably a quad/tri on a hex/tet mesh
448  std::vector< EntityHandle > parents;
449  int dimension = CN::Dimension( TYPE_FROM_HANDLE( *iter ) );
450 
451  // Get the adjacent parent element of "side"
452  if( mbImpl->get_adjacencies( &( *iter ), 1, dimension + 1, false, parents ) != MB_SUCCESS )
453  {
454  MB_SET_ERR( MB_FAILURE, "Couldn't get adjacencies for neuset" );
455  }
456 
457  if( !parents.empty() )
458  {
459  // Make sure the adjacent parent element will be output
460  for( unsigned int k = 0; k < parents.size(); k++ )
461  {
462  result = mbImpl->tag_get_data( mEntityMark, &( parents[k] ), 1, &element_marked );MB_CHK_SET_ERR( result, "Couldn't get mark data" );
463 
464  int side_no, this_sense, this_offset;
465  if( 0x1 == element_marked &&
466  mbImpl->side_number( parents[k], *iter, side_no, this_sense, this_offset ) == MB_SUCCESS &&
467  this_sense == sense )
468  {
469  neuset_data.elements.push_back( parents[k] );
470  neuset_data.side_numbers.push_back( side_no + 1 );
471  break;
472  }
473  }
474  }
475  else
476  {
477  MB_SET_ERR( MB_FAILURE, "No parent element exists for element in neuset " << neuset_data.id );
478  }
479  }
480  }
481 
482  return MB_SUCCESS;
483 }
484 
485 ErrorCode WriteSLAC::write_nodes( const int num_nodes, const Range& nodes, const int dimension )
486 {
487  // See if should transform coordinates
488  ErrorCode result;
489  Tag trans_tag;
490  result = mbImpl->tag_get_handle( MESH_TRANSFORM_TAG_NAME, 16, MB_TYPE_DOUBLE, trans_tag );
491  bool transform_needed = true;
492  if( result == MB_TAG_NOT_FOUND ) transform_needed = false;
493 
494  int num_coords_to_fill = transform_needed ? 3 : dimension;
495 
496  std::vector< double* > coord_arrays( 3 );
497  coord_arrays[0] = new double[num_nodes];
498  coord_arrays[1] = new double[num_nodes];
499  coord_arrays[2] = NULL;
500 
501  if( num_coords_to_fill == 3 ) coord_arrays[2] = new double[num_nodes];
502 
503  result = mWriteIface->get_node_coords( dimension, num_nodes, nodes, mGlobalIdTag, 0, coord_arrays );
504  if( result != MB_SUCCESS )
505  {
506  delete[] coord_arrays[0];
507  delete[] coord_arrays[1];
508  if( coord_arrays[2] ) delete[] coord_arrays[2];
509  return result;
510  }
511 
512  if( transform_needed )
513  {
514  double trans_matrix[16];
515  const EntityHandle mesh = 0;
516  result = mbImpl->tag_get_data( trans_tag, &mesh, 1, trans_matrix );MB_CHK_SET_ERR( result, "Couldn't get transform data" );
517 
518  for( int i = 0; i < num_nodes; i++ )
519  {
520  double vec1[3];
521  double vec2[3];
522 
523  vec2[0] = coord_arrays[0][i];
524  vec2[1] = coord_arrays[1][i];
525  vec2[2] = coord_arrays[2][i];
526 
527  for( int row = 0; row < 3; row++ )
528  {
529  vec1[row] = 0.0;
530  for( int col = 0; col < 3; col++ )
531  vec1[row] += ( trans_matrix[( row * 4 ) + col] * vec2[col] );
532  }
533 
534  coord_arrays[0][i] = vec1[0];
535  coord_arrays[1][i] = vec1[1];
536  coord_arrays[2][i] = vec1[2];
537  }
538  }
539 
540  // Write the nodes
541  int nc_var = -1;
542  std::vector< int > dims;
543  GET_VAR( "coords", nc_var, dims );
544  if( -1 == nc_var ) return MB_FAILURE;
545  size_t start[2] = { 0, 0 }, count[2] = { static_cast< size_t >( num_nodes ), 1 };
546  int fail = nc_put_vara_double( ncFile, nc_var, start, count, coord_arrays[0] );
547  if( NC_NOERR != fail ) return MB_FAILURE;
548  start[1] = 1;
549  fail = nc_put_vara_double( ncFile, nc_var, start, count, coord_arrays[1] );
550  if( NC_NOERR != fail ) return MB_FAILURE;
551  start[1] = 2;
552  fail = nc_put_vara_double( ncFile, nc_var, start, count, coord_arrays[2] );
553  if( NC_NOERR != fail ) return MB_FAILURE;
554 
555  delete[] coord_arrays[0];
556  delete[] coord_arrays[1];
557  if( coord_arrays[2] ) delete[] coord_arrays[2];
558 
559  return MB_SUCCESS;
560 }
561 
563  std::vector< WriteSLAC::MaterialSetData >& matset_data,
564  std::vector< WriteSLAC::NeumannSetData >& neuset_data )
565 {
566  // Need to assign a tag with the matset id
567  Tag matset_id_tag;
568  unsigned int i;
569  int dum = -1;
570  ErrorCode result =
571  mbImpl->tag_get_handle( "__matset_id", 4, MB_TYPE_INTEGER, matset_id_tag, MB_TAG_DENSE | MB_TAG_CREAT, &dum );
572  if( MB_SUCCESS != result ) return result;
573 
574  Range::iterator rit;
575  mesh_info.num_int_hexes = mesh_info.num_int_tets = 0;
576 
577  for( i = 0; i < matset_data.size(); i++ )
578  {
579  WriteSLAC::MaterialSetData matset = matset_data[i];
580  if( matset.moab_type == MBHEX )
581  mesh_info.num_int_hexes += matset.elements->size();
582  else if( matset.moab_type == MBTET )
583  mesh_info.num_int_tets += matset.elements->size();
584  else
585  {
586  std::cout << "WriteSLAC doesn't support elements of type " << CN::EntityTypeName( matset.moab_type )
587  << std::endl;
588  continue;
589  }
590 
591  for( rit = matset.elements->begin(); rit != matset.elements->end(); ++rit )
592  {
593  result = mbImpl->tag_set_data( mMatSetIdTag, &( *rit ), 1, &( matset.id ) );
594  if( MB_SUCCESS != result ) return result;
595  }
596  }
597 
598  // Now go through the neumann sets, pulling out the hexes with faces on the
599  // boundary
600  std::vector< EntityHandle >::iterator vit;
601  for( i = 0; i < neuset_data.size(); i++ )
602  {
603  WriteSLAC::NeumannSetData neuset = neuset_data[i];
604  for( vit = neuset.elements.begin(); vit != neuset.elements.end(); ++vit )
605  {
606  if( TYPE_FROM_HANDLE( *vit ) == MBHEX )
607  mesh_info.bdy_hexes.insert( *vit );
608  else if( TYPE_FROM_HANDLE( *vit ) == MBTET )
609  mesh_info.bdy_tets.insert( *vit );
610  }
611  }
612 
613  // Now we have the number of bdy hexes and tets, we know how many interior ones
614  // there are too
615  mesh_info.num_int_hexes -= mesh_info.bdy_hexes.size();
616  mesh_info.num_int_tets -= mesh_info.bdy_tets.size();
617 
618  return MB_SUCCESS;
619 }
620 
622  std::vector< WriteSLAC::MaterialSetData >& matset_data,
623  std::vector< WriteSLAC::NeumannSetData >& neuset_data )
624 {
625  unsigned int i;
626  std::vector< int > connect;
627  const EntityHandle* connecth;
628  int num_connecth;
629  ErrorCode result;
630 
631  // First write the interior hexes
632  int hex_conn = -1;
633  std::vector< int > dims;
634  if( mesh_info.bdy_hexes.size() != 0 || mesh_info.num_int_hexes != 0 )
635  {
636  GET_VAR( "hexahedron_interior", hex_conn, dims );
637  if( -1 == hex_conn ) return MB_FAILURE;
638  }
639  connect.reserve( 13 );
640  Range::iterator rit;
641 
642  int elem_num = 0;
644  size_t start[2] = { 0, 0 }, count[2] = { 1, 1 };
645  int fail;
646  for( i = 0; i < matset_data.size(); i++ )
647  {
648  matset = matset_data[i];
649  if( matset.moab_type != MBHEX ) continue;
650 
651  int id = matset.id;
652  connect[0] = id;
653 
654  for( rit = matset.elements->begin(); rit != matset.elements->end(); ++rit )
655  {
656  // Skip if it's on the bdy
657  if( mesh_info.bdy_hexes.find( *rit ) != mesh_info.bdy_hexes.end() ) continue;
658 
659  // Get the connectivity of this element
660  result = mbImpl->get_connectivity( *rit, connecth, num_connecth );
661  if( MB_SUCCESS != result ) return result;
662 
663  // Get the vertex ids
664  result = mbImpl->tag_get_data( mGlobalIdTag, connecth, num_connecth, &connect[1] );
665  if( MB_SUCCESS != result ) return result;
666 
667  // Put the variable at the right position
668  start[0] = elem_num++;
669  count[1] = 9;
670 
671  // Write the data
672  fail = nc_put_vara_int( ncFile, hex_conn, start, count, &connect[0] );
673  if( NC_NOERR != fail ) return MB_FAILURE;
674  }
675  }
676 
677  int tet_conn = -1;
678  if( mesh_info.bdy_tets.size() != 0 || mesh_info.num_int_tets != 0 )
679  {
680  GET_VAR( "tetrahedron_interior", tet_conn, dims );
681  if( -1 == tet_conn ) return MB_FAILURE;
682  }
683 
684  // Now the interior tets
685  elem_num = 0;
686  for( i = 0; i < matset_data.size(); i++ )
687  {
688  matset = matset_data[i];
689  if( matset.moab_type != MBTET ) continue;
690 
691  int id = matset.id;
692  connect[0] = id;
693  elem_num = 0;
694  for( rit = matset.elements->begin(); rit != matset.elements->end(); ++rit )
695  {
696  // Skip if it's on the bdy
697  if( mesh_info.bdy_tets.find( *rit ) != mesh_info.bdy_tets.end() ) continue;
698 
699  // Get the connectivity of this element
700  result = mbImpl->get_connectivity( *rit, connecth, num_connecth );
701  if( MB_SUCCESS != result ) return result;
702 
703  // Get the vertex ids
704  result = mbImpl->tag_get_data( mGlobalIdTag, connecth, num_connecth, &connect[1] );
705  if( MB_SUCCESS != result ) return result;
706 
707  // Put the variable at the right position
708  start[0] = elem_num++;
709  count[1] = 5;
710  fail = nc_put_vara_int( ncFile, tet_conn, start, count, &connect[0] );
711  // Write the data
712  if( NC_NOERR != fail ) return MB_FAILURE;
713  }
714  }
715 
716  // Now the exterior hexes
717  if( mesh_info.bdy_hexes.size() != 0 )
718  {
719  hex_conn = -1;
720  GET_VAR( "hexahedron_exterior", hex_conn, dims );
721  if( -1 == hex_conn ) return MB_FAILURE;
722 
723  connect.reserve( 15 );
724  elem_num = 0;
725 
726  // Write the elements
727  for( rit = mesh_info.bdy_hexes.begin(); rit != mesh_info.bdy_hexes.end(); ++rit )
728  {
729  // Get the material set for this hex
730  result = mbImpl->tag_get_data( mMatSetIdTag, &( *rit ), 1, &connect[0] );
731  if( MB_SUCCESS != result ) return result;
732 
733  // Get the connectivity of this element
734  result = mbImpl->get_connectivity( *rit, connecth, num_connecth );
735  if( MB_SUCCESS != result ) return result;
736 
737  // Get the vertex ids
738  result = mbImpl->tag_get_data( mGlobalIdTag, connecth, num_connecth, &connect[1] );
739  if( MB_SUCCESS != result ) return result;
740 
741  // Preset side numbers
742  for( i = 9; i < 15; i++ )
743  connect[i] = -1;
744 
745  // Now write the side numbers
746  for( i = 0; i < neuset_data.size(); i++ )
747  {
748  std::vector< EntityHandle >::iterator vit =
749  std::find( neuset_data[i].elements.begin(), neuset_data[i].elements.end(), *rit );
750  while( vit != neuset_data[i].elements.end() )
751  {
752  // Have a side - get the side # and put in connect array
753  int side_no = neuset_data[i].side_numbers[vit - neuset_data[i].elements.begin()];
754  connect[9 + side_no] = neuset_data[i].id;
755  ++vit;
756  vit = std::find( vit, neuset_data[i].elements.end(), *rit );
757  }
758  }
759 
760  // Put the variable at the right position
761  start[0] = elem_num++;
762  count[1] = 15;
763  fail = nc_put_vara_int( ncFile, hex_conn, start, count, &connect[0] );
764  // Write the data
765  if( NC_NOERR != fail ) return MB_FAILURE;
766  }
767  }
768 
769  // Now the exterior tets
770  if( mesh_info.bdy_tets.size() != 0 )
771  {
772  tet_conn = -1;
773  GET_VAR( "tetrahedron_exterior", tet_conn, dims );
774  if( -1 == tet_conn ) return MB_FAILURE;
775 
776  connect.reserve( 9 );
777  elem_num = 0;
778 
779  // Write the elements
780  for( rit = mesh_info.bdy_tets.begin(); rit != mesh_info.bdy_tets.end(); ++rit )
781  {
782  // Get the material set for this tet
783  result = mbImpl->tag_get_data( mMatSetIdTag, &( *rit ), 1, &connect[0] );
784  if( MB_SUCCESS != result ) return result;
785 
786  // Get the connectivity of this element
787  result = mbImpl->get_connectivity( *rit, connecth, num_connecth );
788  if( MB_SUCCESS != result ) return result;
789 
790  // Get the vertex ids
791  result = mbImpl->tag_get_data( mGlobalIdTag, connecth, num_connecth, &connect[1] );
792  if( MB_SUCCESS != result ) return result;
793 
794  // Preset side numbers
795  for( i = 5; i < 9; i++ )
796  connect[i] = -1;
797 
798  // Now write the side numbers
799  for( i = 0; i < neuset_data.size(); i++ )
800  {
801  std::vector< EntityHandle >::iterator vit =
802  std::find( neuset_data[i].elements.begin(), neuset_data[i].elements.end(), *rit );
803  while( vit != neuset_data[i].elements.end() )
804  {
805  // Have a side - get the side # and put in connect array
806  int side_no = neuset_data[i].side_numbers[vit - neuset_data[i].elements.begin()];
807  connect[5 + side_no] = neuset_data[i].id;
808  ++vit;
809  vit = std::find( vit, neuset_data[i].elements.end(), *rit );
810  }
811  }
812 
813  // Put the variable at the right position
814  start[0] = elem_num++;
815  count[1] = 9;
816  fail = nc_put_vara_int( ncFile, tet_conn, start, count, &connect[0] );
817  // Write the data
818  if( NC_NOERR != fail ) return MB_FAILURE;
819  }
820  }
821 
822  return MB_SUCCESS;
823 }
824 
826 {
827  // Perform the initializations
828 
829  int coord_size = -1, ncoords = -1;
830  // Initialization to avoid warnings on Linux
831  int hexinterior = -1, hexinteriorsize, hexexterior = -1, hexexteriorsize = -1;
832  int tetinterior = -1, tetinteriorsize, tetexterior = -1, tetexteriorsize = -1;
833 
834  if( nc_def_dim( ncFile, "coord_size", (size_t)mesh_info.num_dim, &coord_size ) != NC_NOERR )
835  {
836  MB_SET_ERR( MB_FAILURE, "WriteSLAC: failed to define number of dimensions" );
837  }
838 
839  if( nc_def_dim( ncFile, "ncoords", (size_t)mesh_info.num_nodes, &ncoords ) != NC_NOERR )
840  {
841  MB_SET_ERR( MB_FAILURE, "WriteSLAC: failed to define number of nodes" );
842  }
843 
844  if( 0 != mesh_info.num_int_hexes &&
845  nc_def_dim( ncFile, "hexinterior", (size_t)mesh_info.num_int_hexes, &hexinterior ) != NC_NOERR )
846  {
847  MB_SET_ERR( MB_FAILURE, "WriteSLAC: failed to define number of interior hex elements" );
848  }
849 
850  if( nc_def_dim( ncFile, "hexinteriorsize", (size_t)9, &hexinteriorsize ) != NC_NOERR )
851  {
852  MB_SET_ERR( MB_FAILURE, "WriteSLAC: failed to define interior hex element size" );
853  }
854 
855  if( 0 != mesh_info.bdy_hexes.size() &&
856  nc_def_dim( ncFile, "hexexterior", (size_t)mesh_info.bdy_hexes.size(), &hexexterior ) != NC_NOERR )
857  {
858  MB_SET_ERR( MB_FAILURE, "WriteSLAC: failed to define number of exterior hex elements" );
859  }
860 
861  if( nc_def_dim( ncFile, "hexexteriorsize", (size_t)15, &hexexteriorsize ) != NC_NOERR )
862  {
863  MB_SET_ERR( MB_FAILURE, "WriteSLAC: failed to define exterior hex element size" );
864  }
865 
866  if( 0 != mesh_info.num_int_tets &&
867  nc_def_dim( ncFile, "tetinterior", (size_t)mesh_info.num_int_tets, &tetinterior ) != NC_NOERR )
868  {
869  MB_SET_ERR( MB_FAILURE, "WriteSLAC: failed to define number of interior tet elements" );
870  }
871 
872  if( nc_def_dim( ncFile, "tetinteriorsize", (size_t)5, &tetinteriorsize ) != NC_NOERR )
873  {
874  MB_SET_ERR( MB_FAILURE, "WriteSLAC: failed to define interior tet element size" );
875  }
876 
877  if( 0 != mesh_info.bdy_tets.size() &&
878  nc_def_dim( ncFile, "tetexterior", (size_t)mesh_info.bdy_tets.size(), &tetexterior ) != NC_NOERR )
879  {
880  MB_SET_ERR( MB_FAILURE, "WriteSLAC: failed to define number of exterior tet elements" );
881  }
882 
883  if( nc_def_dim( ncFile, "tetexteriorsize", (size_t)9, &tetexteriorsize ) != NC_NOERR )
884  {
885  MB_SET_ERR( MB_FAILURE, "WriteSLAC: failed to define exterior tet element size" );
886  }
887 
888  /* ...and some variables */
889 
890  int dims[2];
891  dims[0] = hexinterior;
892  dims[1] = hexinteriorsize;
893  int dum_var;
894  if( 0 != mesh_info.num_int_hexes &&
895  NC_NOERR != nc_def_var( ncFile, "hexahedron_interior", NC_LONG, 2, dims, &dum_var ) )
896  {
897  MB_SET_ERR( MB_FAILURE, "WriteSLAC: failed to create connectivity array for interior hexes" );
898  }
899 
900  dims[0] = hexexterior;
901  dims[1] = hexexteriorsize;
902  if( 0 != mesh_info.bdy_hexes.size() &&
903  NC_NOERR != nc_def_var( ncFile, "hexahedron_exterior", NC_LONG, 2, dims, &dum_var ) )
904  {
905  MB_SET_ERR( MB_FAILURE, "WriteSLAC: failed to create connectivity array for exterior hexes" );
906  }
907 
908  dims[0] = tetinterior;
909  dims[1] = tetinteriorsize;
910  if( 0 != mesh_info.num_int_tets &&
911  NC_NOERR != nc_def_var( ncFile, "tetrahedron_exterior", NC_LONG, 2, dims, &dum_var ) )
912  {
913  MB_SET_ERR( MB_FAILURE, "WriteSLAC: failed to create connectivity array for interior tets" );
914  }
915 
916  dims[0] = tetexterior;
917  dims[1] = tetexteriorsize;
918  if( 0 != mesh_info.bdy_tets.size() &&
919  NC_NOERR != nc_def_var( ncFile, "tetrahedron_exterior", NC_LONG, 2, dims, &dum_var ) )
920  {
921  MB_SET_ERR( MB_FAILURE, "WriteSLAC: failed to create connectivity array for exterior tets" );
922  }
923 
924  /* Node coordinate arrays: */
925 
926  dims[0] = ncoords;
927  dims[1] = coord_size;
928  if( NC_NOERR != nc_def_var( ncFile, "coords", NC_DOUBLE, 2, dims, &dum_var ) )
929  {
930  MB_SET_ERR( MB_FAILURE, "WriteSLAC: failed to define node coordinate array" );
931  }
932 
933  return MB_SUCCESS;
934 }
935 
936 ErrorCode WriteSLAC::open_file( const char* filename )
937 {
938  // Not a valid filname
939  if( strlen( (const char*)filename ) == 0 )
940  {
941  MB_SET_ERR( MB_FAILURE, "Output filename not specified" );
942  }
943 
944  int fail = nc_create( filename, NC_CLOBBER, &ncFile );
945  // File couldn't be opened
946  if( NC_NOERR != fail )
947  {
948  MB_SET_ERR( MB_FAILURE, "Cannot open " << filename );
949  }
950 
951  return MB_SUCCESS;
952 }
953 
955  int current_sense,
956  Range& forward_elems,
957  Range& reverse_elems )
958 {
959  Range ss_elems, ss_meshsets;
960 
961  // Get the sense tag; don't need to check return, might be an error if the tag
962  // hasn't been created yet
963  Tag sense_tag = 0;
964  mbImpl->tag_get_handle( "SENSE", 1, MB_TYPE_INTEGER, sense_tag );
965 
966  // Get the entities in this set
967  ErrorCode result = mbImpl->get_entities_by_handle( neuset, ss_elems, true );
968  if( MB_FAILURE == result ) return result;
969 
970  // Now remove the meshsets into the ss_meshsets; first find the first meshset,
971  Range::iterator range_iter = ss_elems.begin();
972  while( TYPE_FROM_HANDLE( *range_iter ) != MBENTITYSET && range_iter != ss_elems.end() )
973  ++range_iter;
974 
975  // Then, if there are some, copy them into ss_meshsets and erase from ss_elems
976  if( range_iter != ss_elems.end() )
977  {
978  std::copy( range_iter, ss_elems.end(), range_inserter( ss_meshsets ) );
979  ss_elems.erase( range_iter, ss_elems.end() );
980  }
981 
982  // OK, for the elements, check the sense of this set and copy into the right range
983  // (if the sense is 0, copy into both ranges)
984 
985  // Need to step forward on list until we reach the right dimension
986  Range::iterator dum_it = ss_elems.end();
987  --dum_it;
988  int target_dim = CN::Dimension( TYPE_FROM_HANDLE( *dum_it ) );
989  dum_it = ss_elems.begin();
990  while( target_dim != CN::Dimension( TYPE_FROM_HANDLE( *dum_it ) ) && dum_it != ss_elems.end() )
991  ++dum_it;
992 
993  if( current_sense == 1 || current_sense == 0 ) std::copy( dum_it, ss_elems.end(), range_inserter( forward_elems ) );
994  if( current_sense == -1 || current_sense == 0 )
995  std::copy( dum_it, ss_elems.end(), range_inserter( reverse_elems ) );
996 
997  // Now loop over the contained meshsets, getting the sense of those and calling this
998  // function recursively
999  for( range_iter = ss_meshsets.begin(); range_iter != ss_meshsets.end(); ++range_iter )
1000  {
1001  // First get the sense; if it's not there, by convention it's forward
1002  int this_sense;
1003  if( 0 == sense_tag || MB_FAILURE == mbImpl->tag_get_data( sense_tag, &( *range_iter ), 1, &this_sense ) )
1004  this_sense = 1;
1005 
1006  // Now get all the entities on this meshset, with the proper (possibly reversed) sense
1007  get_neuset_elems( *range_iter, this_sense * current_sense, forward_elems, reverse_elems );
1008  }
1009 
1010  return result;
1011 }
1012 
1013 } // namespace moab