 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
1 /**
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 );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
562 ErrorCode WriteSLAC::gather_interior_exterior( MeshInfo& mesh_info,
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
621 ErrorCode WriteSLAC::write_matsets( MeshInfo& mesh_info,
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;
643 WriteSLAC::MaterialSetData matset;
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
825 ErrorCode WriteSLAC::initialize_file( MeshInfo& mesh_info )
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
954 ErrorCode WriteSLAC::get_neuset_elems( EntityHandle neuset,
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
1.9.1.