Mesh Oriented datABase  (version 5.6.0)
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 
68 WriterIface* WriteSLAC::factory( Interface* iface )
69 {
70  return new WriteSLAC( iface );
71 }
72 
73 WriteSLAC::WriteSLAC( Interface* impl ) : mbImpl( impl ), ncFile( 0 )
74 {
75  assert( impl != NULL );
76 
77  impl->query_interface( mWriteIface );
78 
79  // Initialize in case tag_get_handle fails below
80  //! get and cache predefined tag handles
81  int negone = -1;
82  impl->tag_get_handle( MATERIAL_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mMaterialSetTag, MB_TAG_SPARSE | MB_TAG_CREAT,
83  &negone );
84 
85  impl->tag_get_handle( DIRICHLET_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mDirichletSetTag, MB_TAG_SPARSE | MB_TAG_CREAT,
86  &negone );
87 
88  impl->tag_get_handle( NEUMANN_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mNeumannSetTag, MB_TAG_SPARSE | MB_TAG_CREAT,
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 
99 WriteSLAC::~WriteSLAC()
100 {
101  mbImpl->release_interface( mWriteIface );
102  mbImpl->tag_delete( mEntityMark );
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 
210 ErrorCode WriteSLAC::gather_mesh_information( MeshInfo& mesh_info,
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
234  mbImpl->tag_delete( mEntityMark );
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 
314  matset_data.number_nodes_per_element = ExoIIUtil::VerticesPerElement[matset_data.element_type];
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 );
390  MB_CHK_SET_ERR( result, "Couldn't get mark data" );
391 
392  if( 0x1 == node_marked ) dirset_data.nodes.push_back( *iter );
393  }
394 
395  dirset_data.number_nodes = dirset_data.nodes.size();
396  dirset_info.push_back( dirset_data );
397  }
398 
399  //------neusets--------
400  vector_iter = neusets.begin();
401  end_vector_iter = neusets.end();
402 
403  for( ; vector_iter != end_vector_iter; ++vector_iter )
404  {
405  WriteSLAC::NeumannSetData neuset_data;
406 
407  // Get the neuset's id
408  if( mbImpl->tag_get_data( mNeumannSetTag, &( *vector_iter ), 1, &id ) != MB_SUCCESS ) return MB_FAILURE;
409 
410  neuset_data.id = id;
411  neuset_data.mesh_set_handle = *vector_iter;
412 
413  // Get the sides in two lists, one forward the other reverse; starts with forward sense
414  // by convention
415  Range forward_elems, reverse_elems;
416  if( get_neuset_elems( *vector_iter, 0, forward_elems, reverse_elems ) == MB_FAILURE ) return MB_FAILURE;
417 
418  ErrorCode result = get_valid_sides( forward_elems, 1, neuset_data );
419  MB_CHK_SET_ERR( result, "Couldn't get valid sides data" );
420  result = get_valid_sides( reverse_elems, -1, neuset_data );
421  MB_CHK_SET_ERR( result, "Couldn't get valid sides data" );
422 
423  neuset_data.number_elements = neuset_data.elements.size();
424  neuset_info.push_back( neuset_data );
425  }
426 
427  // Get information about interior/exterior tets/hexes, and mark matset ids
428  return gather_interior_exterior( mesh_info, matset_info, neuset_info );
429 }
430 
431 ErrorCode WriteSLAC::get_valid_sides( Range& elems, const int sense, WriteSLAC::NeumannSetData& neuset_data )
432 {
433  // This is where we see if underlying element of side set element is included in output
434 
435  unsigned char element_marked = 0;
436  ErrorCode result;
437  for( Range::iterator iter = elems.begin(); iter != elems.end(); ++iter )
438  {
439  // Should insert here if "side" is a quad/tri on a quad/tri mesh
440  result = mbImpl->tag_get_data( mEntityMark, &( *iter ), 1, &element_marked );
441  MB_CHK_SET_ERR( result, "Couldn't get mark data" );
442 
443  if( 0x1 == element_marked )
444  {
445  neuset_data.elements.push_back( *iter );
446 
447  // TJT TODO: the sense should really be # edges + 1or2
448  neuset_data.side_numbers.push_back( ( sense == 1 ? 1 : 2 ) );
449  }
450  else
451  { // Then "side" is probably a quad/tri on a hex/tet mesh
452  std::vector< EntityHandle > parents;
453  int dimension = CN::Dimension( TYPE_FROM_HANDLE( *iter ) );
454 
455  // Get the adjacent parent element of "side"
456  if( mbImpl->get_adjacencies( &( *iter ), 1, dimension + 1, false, parents ) != MB_SUCCESS )
457  {
458  MB_SET_ERR( MB_FAILURE, "Couldn't get adjacencies for neuset" );
459  }
460 
461  if( !parents.empty() )
462  {
463  // Make sure the adjacent parent element will be output
464  for( unsigned int k = 0; k < parents.size(); k++ )
465  {
466  result = mbImpl->tag_get_data( mEntityMark, &( parents[k] ), 1, &element_marked );
467  MB_CHK_SET_ERR( result, "Couldn't get mark data" );
468 
469  int side_no, this_sense, this_offset;
470  if( 0x1 == element_marked &&
471  mbImpl->side_number( parents[k], *iter, side_no, this_sense, this_offset ) == MB_SUCCESS &&
472  this_sense == sense )
473  {
474  neuset_data.elements.push_back( parents[k] );
475  neuset_data.side_numbers.push_back( side_no + 1 );
476  break;
477  }
478  }
479  }
480  else
481  {
482  MB_SET_ERR( MB_FAILURE, "No parent element exists for element in neuset " << neuset_data.id );
483  }
484  }
485  }
486 
487  return MB_SUCCESS;
488 }
489 
490 ErrorCode WriteSLAC::write_nodes( const int num_nodes, const Range& nodes, const int dimension )
491 {
492  // See if should transform coordinates
493  ErrorCode result;
494  Tag trans_tag;
495  result = mbImpl->tag_get_handle( MESH_TRANSFORM_TAG_NAME, 16, MB_TYPE_DOUBLE, trans_tag );
496  bool transform_needed = true;
497  if( result == MB_TAG_NOT_FOUND ) transform_needed = false;
498 
499  int num_coords_to_fill = transform_needed ? 3 : dimension;
500 
501  std::vector< double* > coord_arrays( 3 );
502  coord_arrays[0] = new double[num_nodes];
503  coord_arrays[1] = new double[num_nodes];
504  coord_arrays[2] = NULL;
505 
506  if( num_coords_to_fill == 3 ) coord_arrays[2] = new double[num_nodes];
507 
508  result = mWriteIface->get_node_coords( dimension, num_nodes, nodes, mGlobalIdTag, 0, coord_arrays );
509  if( result != MB_SUCCESS )
510  {
511  delete[] coord_arrays[0];
512  delete[] coord_arrays[1];
513  if( coord_arrays[2] ) delete[] coord_arrays[2];
514  return result;
515  }
516 
517  if( transform_needed )
518  {
519  double trans_matrix[16];
520  const EntityHandle mesh = 0;
521  result = mbImpl->tag_get_data( trans_tag, &mesh, 1, trans_matrix );
522  MB_CHK_SET_ERR( result, "Couldn't get transform data" );
523 
524  for( int i = 0; i < num_nodes; i++ )
525  {
526  double vec1[3];
527  double vec2[3];
528 
529  vec2[0] = coord_arrays[0][i];
530  vec2[1] = coord_arrays[1][i];
531  vec2[2] = coord_arrays[2][i];
532 
533  for( int row = 0; row < 3; row++ )
534  {
535  vec1[row] = 0.0;
536  for( int col = 0; col < 3; col++ )
537  vec1[row] += ( trans_matrix[( row * 4 ) + col] * vec2[col] );
538  }
539 
540  coord_arrays[0][i] = vec1[0];
541  coord_arrays[1][i] = vec1[1];
542  coord_arrays[2][i] = vec1[2];
543  }
544  }
545 
546  // Write the nodes
547  int nc_var = -1;
548  std::vector< int > dims;
549  GET_VAR( "coords", nc_var, dims );
550  if( -1 == nc_var ) return MB_FAILURE;
551  size_t start[2] = { 0, 0 }, count[2] = { static_cast< size_t >( num_nodes ), 1 };
552  int fail = nc_put_vara_double( ncFile, nc_var, start, count, coord_arrays[0] );
553  if( NC_NOERR != fail ) return MB_FAILURE;
554  start[1] = 1;
555  fail = nc_put_vara_double( ncFile, nc_var, start, count, coord_arrays[1] );
556  if( NC_NOERR != fail ) return MB_FAILURE;
557  start[1] = 2;
558  fail = nc_put_vara_double( ncFile, nc_var, start, count, coord_arrays[2] );
559  if( NC_NOERR != fail ) return MB_FAILURE;
560 
561  delete[] coord_arrays[0];
562  delete[] coord_arrays[1];
563  if( coord_arrays[2] ) delete[] coord_arrays[2];
564 
565  return MB_SUCCESS;
566 }
567 
568 ErrorCode WriteSLAC::gather_interior_exterior( MeshInfo& mesh_info,
569  std::vector< WriteSLAC::MaterialSetData >& matset_data,
570  std::vector< WriteSLAC::NeumannSetData >& neuset_data )
571 {
572  // Need to assign a tag with the matset id
573  Tag matset_id_tag;
574  unsigned int i;
575  int dum = -1;
576  ErrorCode result =
577  mbImpl->tag_get_handle( "__matset_id", 4, MB_TYPE_INTEGER, matset_id_tag, MB_TAG_DENSE | MB_TAG_CREAT, &dum );
578  if( MB_SUCCESS != result ) return result;
579 
580  Range::iterator rit;
581  mesh_info.num_int_hexes = mesh_info.num_int_tets = 0;
582 
583  for( i = 0; i < matset_data.size(); i++ )
584  {
585  WriteSLAC::MaterialSetData matset = matset_data[i];
586  if( matset.moab_type == MBHEX )
587  mesh_info.num_int_hexes += matset.elements->size();
588  else if( matset.moab_type == MBTET )
589  mesh_info.num_int_tets += matset.elements->size();
590  else
591  {
592  std::cout << "WriteSLAC doesn't support elements of type " << CN::EntityTypeName( matset.moab_type )
593  << std::endl;
594  continue;
595  }
596 
597  for( rit = matset.elements->begin(); rit != matset.elements->end(); ++rit )
598  {
599  result = mbImpl->tag_set_data( mMatSetIdTag, &( *rit ), 1, &( matset.id ) );
600  if( MB_SUCCESS != result ) return result;
601  }
602  }
603 
604  // Now go through the neumann sets, pulling out the hexes with faces on the
605  // boundary
606  std::vector< EntityHandle >::iterator vit;
607  for( i = 0; i < neuset_data.size(); i++ )
608  {
609  WriteSLAC::NeumannSetData neuset = neuset_data[i];
610  for( vit = neuset.elements.begin(); vit != neuset.elements.end(); ++vit )
611  {
612  if( TYPE_FROM_HANDLE( *vit ) == MBHEX )
613  mesh_info.bdy_hexes.insert( *vit );
614  else if( TYPE_FROM_HANDLE( *vit ) == MBTET )
615  mesh_info.bdy_tets.insert( *vit );
616  }
617  }
618 
619  // Now we have the number of bdy hexes and tets, we know how many interior ones
620  // there are too
621  mesh_info.num_int_hexes -= mesh_info.bdy_hexes.size();
622  mesh_info.num_int_tets -= mesh_info.bdy_tets.size();
623 
624  return MB_SUCCESS;
625 }
626 
627 ErrorCode WriteSLAC::write_matsets( MeshInfo& mesh_info,
628  std::vector< WriteSLAC::MaterialSetData >& matset_data,
629  std::vector< WriteSLAC::NeumannSetData >& neuset_data )
630 {
631  unsigned int i;
632  std::vector< int > connect;
633  const EntityHandle* connecth;
634  int num_connecth;
635  ErrorCode result;
636 
637  // First write the interior hexes
638  int hex_conn = -1;
639  std::vector< int > dims;
640  if( mesh_info.bdy_hexes.size() != 0 || mesh_info.num_int_hexes != 0 )
641  {
642  GET_VAR( "hexahedron_interior", hex_conn, dims );
643  if( -1 == hex_conn ) return MB_FAILURE;
644  }
645  connect.reserve( 13 );
646  Range::iterator rit;
647 
648  int elem_num = 0;
649  WriteSLAC::MaterialSetData matset;
650  size_t start[2] = { 0, 0 }, count[2] = { 1, 1 };
651  int fail;
652  for( i = 0; i < matset_data.size(); i++ )
653  {
654  matset = matset_data[i];
655  if( matset.moab_type != MBHEX ) continue;
656 
657  int id = matset.id;
658  connect[0] = id;
659 
660  for( rit = matset.elements->begin(); rit != matset.elements->end(); ++rit )
661  {
662  // Skip if it's on the bdy
663  if( mesh_info.bdy_hexes.find( *rit ) != mesh_info.bdy_hexes.end() ) continue;
664 
665  // Get the connectivity of this element
666  result = mbImpl->get_connectivity( *rit, connecth, num_connecth );
667  if( MB_SUCCESS != result ) return result;
668 
669  // Get the vertex ids
670  result = mbImpl->tag_get_data( mGlobalIdTag, connecth, num_connecth, &connect[1] );
671  if( MB_SUCCESS != result ) return result;
672 
673  // Put the variable at the right position
674  start[0] = elem_num++;
675  count[1] = 9;
676 
677  // Write the data
678  fail = nc_put_vara_int( ncFile, hex_conn, start, count, &connect[0] );
679  if( NC_NOERR != fail ) return MB_FAILURE;
680  }
681  }
682 
683  int tet_conn = -1;
684  if( mesh_info.bdy_tets.size() != 0 || mesh_info.num_int_tets != 0 )
685  {
686  GET_VAR( "tetrahedron_interior", tet_conn, dims );
687  if( -1 == tet_conn ) return MB_FAILURE;
688  }
689 
690  // Now the interior tets
691  elem_num = 0;
692  for( i = 0; i < matset_data.size(); i++ )
693  {
694  matset = matset_data[i];
695  if( matset.moab_type != MBTET ) continue;
696 
697  int id = matset.id;
698  connect[0] = id;
699  elem_num = 0;
700  for( rit = matset.elements->begin(); rit != matset.elements->end(); ++rit )
701  {
702  // Skip if it's on the bdy
703  if( mesh_info.bdy_tets.find( *rit ) != mesh_info.bdy_tets.end() ) continue;
704 
705  // Get the connectivity of this element
706  result = mbImpl->get_connectivity( *rit, connecth, num_connecth );
707  if( MB_SUCCESS != result ) return result;
708 
709  // Get the vertex ids
710  result = mbImpl->tag_get_data( mGlobalIdTag, connecth, num_connecth, &connect[1] );
711  if( MB_SUCCESS != result ) return result;
712 
713  // Put the variable at the right position
714  start[0] = elem_num++;
715  count[1] = 5;
716  fail = nc_put_vara_int( ncFile, tet_conn, start, count, &connect[0] );
717  // Write the data
718  if( NC_NOERR != fail ) return MB_FAILURE;
719  }
720  }
721 
722  // Now the exterior hexes
723  if( mesh_info.bdy_hexes.size() != 0 )
724  {
725  hex_conn = -1;
726  GET_VAR( "hexahedron_exterior", hex_conn, dims );
727  if( -1 == hex_conn ) return MB_FAILURE;
728 
729  connect.reserve( 15 );
730  elem_num = 0;
731 
732  // Write the elements
733  for( rit = mesh_info.bdy_hexes.begin(); rit != mesh_info.bdy_hexes.end(); ++rit )
734  {
735  // Get the material set for this hex
736  result = mbImpl->tag_get_data( mMatSetIdTag, &( *rit ), 1, &connect[0] );
737  if( MB_SUCCESS != result ) return result;
738 
739  // Get the connectivity of this element
740  result = mbImpl->get_connectivity( *rit, connecth, num_connecth );
741  if( MB_SUCCESS != result ) return result;
742 
743  // Get the vertex ids
744  result = mbImpl->tag_get_data( mGlobalIdTag, connecth, num_connecth, &connect[1] );
745  if( MB_SUCCESS != result ) return result;
746 
747  // Preset side numbers
748  for( i = 9; i < 15; i++ )
749  connect[i] = -1;
750 
751  // Now write the side numbers
752  for( i = 0; i < neuset_data.size(); i++ )
753  {
754  std::vector< EntityHandle >::iterator vit =
755  std::find( neuset_data[i].elements.begin(), neuset_data[i].elements.end(), *rit );
756  while( vit != neuset_data[i].elements.end() )
757  {
758  // Have a side - get the side # and put in connect array
759  int side_no = neuset_data[i].side_numbers[vit - neuset_data[i].elements.begin()];
760  connect[9 + side_no] = neuset_data[i].id;
761  ++vit;
762  vit = std::find( vit, neuset_data[i].elements.end(), *rit );
763  }
764  }
765 
766  // Put the variable at the right position
767  start[0] = elem_num++;
768  count[1] = 15;
769  fail = nc_put_vara_int( ncFile, hex_conn, start, count, &connect[0] );
770  // Write the data
771  if( NC_NOERR != fail ) return MB_FAILURE;
772  }
773  }
774 
775  // Now the exterior tets
776  if( mesh_info.bdy_tets.size() != 0 )
777  {
778  tet_conn = -1;
779  GET_VAR( "tetrahedron_exterior", tet_conn, dims );
780  if( -1 == tet_conn ) return MB_FAILURE;
781 
782  connect.reserve( 9 );
783  elem_num = 0;
784 
785  // Write the elements
786  for( rit = mesh_info.bdy_tets.begin(); rit != mesh_info.bdy_tets.end(); ++rit )
787  {
788  // Get the material set for this tet
789  result = mbImpl->tag_get_data( mMatSetIdTag, &( *rit ), 1, &connect[0] );
790  if( MB_SUCCESS != result ) return result;
791 
792  // Get the connectivity of this element
793  result = mbImpl->get_connectivity( *rit, connecth, num_connecth );
794  if( MB_SUCCESS != result ) return result;
795 
796  // Get the vertex ids
797  result = mbImpl->tag_get_data( mGlobalIdTag, connecth, num_connecth, &connect[1] );
798  if( MB_SUCCESS != result ) return result;
799 
800  // Preset side numbers
801  for( i = 5; i < 9; i++ )
802  connect[i] = -1;
803 
804  // Now write the side numbers
805  for( i = 0; i < neuset_data.size(); i++ )
806  {
807  std::vector< EntityHandle >::iterator vit =
808  std::find( neuset_data[i].elements.begin(), neuset_data[i].elements.end(), *rit );
809  while( vit != neuset_data[i].elements.end() )
810  {
811  // Have a side - get the side # and put in connect array
812  int side_no = neuset_data[i].side_numbers[vit - neuset_data[i].elements.begin()];
813  connect[5 + side_no] = neuset_data[i].id;
814  ++vit;
815  vit = std::find( vit, neuset_data[i].elements.end(), *rit );
816  }
817  }
818 
819  // Put the variable at the right position
820  start[0] = elem_num++;
821  count[1] = 9;
822  fail = nc_put_vara_int( ncFile, tet_conn, start, count, &connect[0] );
823  // Write the data
824  if( NC_NOERR != fail ) return MB_FAILURE;
825  }
826  }
827 
828  return MB_SUCCESS;
829 }
830 
831 ErrorCode WriteSLAC::initialize_file( MeshInfo& mesh_info )
832 {
833  // Perform the initializations
834 
835  int coord_size = -1, ncoords = -1;
836  // Initialization to avoid warnings on Linux
837  int hexinterior = -1, hexinteriorsize, hexexterior = -1, hexexteriorsize = -1;
838  int tetinterior = -1, tetinteriorsize, tetexterior = -1, tetexteriorsize = -1;
839 
840  if( nc_def_dim( ncFile, "coord_size", (size_t)mesh_info.num_dim, &coord_size ) != NC_NOERR )
841  {
842  MB_SET_ERR( MB_FAILURE, "WriteSLAC: failed to define number of dimensions" );
843  }
844 
845  if( nc_def_dim( ncFile, "ncoords", (size_t)mesh_info.num_nodes, &ncoords ) != NC_NOERR )
846  {
847  MB_SET_ERR( MB_FAILURE, "WriteSLAC: failed to define number of nodes" );
848  }
849 
850  if( 0 != mesh_info.num_int_hexes &&
851  nc_def_dim( ncFile, "hexinterior", (size_t)mesh_info.num_int_hexes, &hexinterior ) != NC_NOERR )
852  {
853  MB_SET_ERR( MB_FAILURE, "WriteSLAC: failed to define number of interior hex elements" );
854  }
855 
856  if( nc_def_dim( ncFile, "hexinteriorsize", (size_t)9, &hexinteriorsize ) != NC_NOERR )
857  {
858  MB_SET_ERR( MB_FAILURE, "WriteSLAC: failed to define interior hex element size" );
859  }
860 
861  if( 0 != mesh_info.bdy_hexes.size() &&
862  nc_def_dim( ncFile, "hexexterior", (size_t)mesh_info.bdy_hexes.size(), &hexexterior ) != NC_NOERR )
863  {
864  MB_SET_ERR( MB_FAILURE, "WriteSLAC: failed to define number of exterior hex elements" );
865  }
866 
867  if( nc_def_dim( ncFile, "hexexteriorsize", (size_t)15, &hexexteriorsize ) != NC_NOERR )
868  {
869  MB_SET_ERR( MB_FAILURE, "WriteSLAC: failed to define exterior hex element size" );
870  }
871 
872  if( 0 != mesh_info.num_int_tets &&
873  nc_def_dim( ncFile, "tetinterior", (size_t)mesh_info.num_int_tets, &tetinterior ) != NC_NOERR )
874  {
875  MB_SET_ERR( MB_FAILURE, "WriteSLAC: failed to define number of interior tet elements" );
876  }
877 
878  if( nc_def_dim( ncFile, "tetinteriorsize", (size_t)5, &tetinteriorsize ) != NC_NOERR )
879  {
880  MB_SET_ERR( MB_FAILURE, "WriteSLAC: failed to define interior tet element size" );
881  }
882 
883  if( 0 != mesh_info.bdy_tets.size() &&
884  nc_def_dim( ncFile, "tetexterior", (size_t)mesh_info.bdy_tets.size(), &tetexterior ) != NC_NOERR )
885  {
886  MB_SET_ERR( MB_FAILURE, "WriteSLAC: failed to define number of exterior tet elements" );
887  }
888 
889  if( nc_def_dim( ncFile, "tetexteriorsize", (size_t)9, &tetexteriorsize ) != NC_NOERR )
890  {
891  MB_SET_ERR( MB_FAILURE, "WriteSLAC: failed to define exterior tet element size" );
892  }
893 
894  /* ...and some variables */
895 
896  int dims[2];
897  dims[0] = hexinterior;
898  dims[1] = hexinteriorsize;
899  int dum_var;
900  if( 0 != mesh_info.num_int_hexes &&
901  NC_NOERR != nc_def_var( ncFile, "hexahedron_interior", NC_LONG, 2, dims, &dum_var ) )
902  {
903  MB_SET_ERR( MB_FAILURE, "WriteSLAC: failed to create connectivity array for interior hexes" );
904  }
905 
906  dims[0] = hexexterior;
907  dims[1] = hexexteriorsize;
908  if( 0 != mesh_info.bdy_hexes.size() &&
909  NC_NOERR != nc_def_var( ncFile, "hexahedron_exterior", NC_LONG, 2, dims, &dum_var ) )
910  {
911  MB_SET_ERR( MB_FAILURE, "WriteSLAC: failed to create connectivity array for exterior hexes" );
912  }
913 
914  dims[0] = tetinterior;
915  dims[1] = tetinteriorsize;
916  if( 0 != mesh_info.num_int_tets &&
917  NC_NOERR != nc_def_var( ncFile, "tetrahedron_exterior", NC_LONG, 2, dims, &dum_var ) )
918  {
919  MB_SET_ERR( MB_FAILURE, "WriteSLAC: failed to create connectivity array for interior tets" );
920  }
921 
922  dims[0] = tetexterior;
923  dims[1] = tetexteriorsize;
924  if( 0 != mesh_info.bdy_tets.size() &&
925  NC_NOERR != nc_def_var( ncFile, "tetrahedron_exterior", NC_LONG, 2, dims, &dum_var ) )
926  {
927  MB_SET_ERR( MB_FAILURE, "WriteSLAC: failed to create connectivity array for exterior tets" );
928  }
929 
930  /* Node coordinate arrays: */
931 
932  dims[0] = ncoords;
933  dims[1] = coord_size;
934  if( NC_NOERR != nc_def_var( ncFile, "coords", NC_DOUBLE, 2, dims, &dum_var ) )
935  {
936  MB_SET_ERR( MB_FAILURE, "WriteSLAC: failed to define node coordinate array" );
937  }
938 
939  return MB_SUCCESS;
940 }
941 
942 ErrorCode WriteSLAC::open_file( const char* filename )
943 {
944  // Not a valid filname
945  if( strlen( (const char*)filename ) == 0 )
946  {
947  MB_SET_ERR( MB_FAILURE, "Output filename not specified" );
948  }
949 
950  int fail = nc_create( filename, NC_CLOBBER, &ncFile );
951  // File couldn't be opened
952  if( NC_NOERR != fail )
953  {
954  MB_SET_ERR( MB_FAILURE, "Cannot open " << filename );
955  }
956 
957  return MB_SUCCESS;
958 }
959 
960 ErrorCode WriteSLAC::get_neuset_elems( EntityHandle neuset,
961  int current_sense,
962  Range& forward_elems,
963  Range& reverse_elems )
964 {
965  Range ss_elems, ss_meshsets;
966 
967  // Get the sense tag; don't need to check return, might be an error if the tag
968  // hasn't been created yet
969  Tag sense_tag = 0;
970  mbImpl->tag_get_handle( "SENSE", 1, MB_TYPE_INTEGER, sense_tag );
971 
972  // Get the entities in this set
973  ErrorCode result = mbImpl->get_entities_by_handle( neuset, ss_elems, true );
974  if( MB_FAILURE == result ) return result;
975 
976  // Now remove the meshsets into the ss_meshsets; first find the first meshset,
977  Range::iterator range_iter = ss_elems.begin();
978  while( TYPE_FROM_HANDLE( *range_iter ) != MBENTITYSET && range_iter != ss_elems.end() )
979  ++range_iter;
980 
981  // Then, if there are some, copy them into ss_meshsets and erase from ss_elems
982  if( range_iter != ss_elems.end() )
983  {
984  std::copy( range_iter, ss_elems.end(), range_inserter( ss_meshsets ) );
985  ss_elems.erase( range_iter, ss_elems.end() );
986  }
987 
988  // OK, for the elements, check the sense of this set and copy into the right range
989  // (if the sense is 0, copy into both ranges)
990 
991  // Need to step forward on list until we reach the right dimension
992  Range::iterator dum_it = ss_elems.end();
993  --dum_it;
994  int target_dim = CN::Dimension( TYPE_FROM_HANDLE( *dum_it ) );
995  dum_it = ss_elems.begin();
996  while( target_dim != CN::Dimension( TYPE_FROM_HANDLE( *dum_it ) ) && dum_it != ss_elems.end() )
997  ++dum_it;
998 
999  if( current_sense == 1 || current_sense == 0 ) std::copy( dum_it, ss_elems.end(), range_inserter( forward_elems ) );
1000  if( current_sense == -1 || current_sense == 0 )
1001  std::copy( dum_it, ss_elems.end(), range_inserter( reverse_elems ) );
1002 
1003  // Now loop over the contained meshsets, getting the sense of those and calling this
1004  // function recursively
1005  for( range_iter = ss_meshsets.begin(); range_iter != ss_meshsets.end(); ++range_iter )
1006  {
1007  // First get the sense; if it's not there, by convention it's forward
1008  int this_sense;
1009  if( 0 == sense_tag || MB_FAILURE == mbImpl->tag_get_data( sense_tag, &( *range_iter ), 1, &this_sense ) )
1010  this_sense = 1;
1011 
1012  // Now get all the entities on this meshset, with the proper (possibly reversed) sense
1013  get_neuset_elems( *range_iter, this_sense * current_sense, forward_elems, reverse_elems );
1014  }
1015 
1016  return result;
1017 }
1018 
1019 } // namespace moab