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