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 "WriteNCDF.hpp"
25 : :
26 : : #include "netcdf.h"
27 : : #include <utility>
28 : : #include <algorithm>
29 : : #include <time.h>
30 : : #include <string>
31 : : #include <vector>
32 : : #include <stdio.h>
33 : : #include <string.h>
34 : : #include <assert.h>
35 : :
36 : : #include "moab/Interface.hpp"
37 : : #include "moab/Range.hpp"
38 : : #include "moab/CN.hpp"
39 : : #include "moab/FileOptions.hpp"
40 : : #include "MBTagConventions.hpp"
41 : : #include "Internals.hpp"
42 : : #include "ExoIIUtil.hpp"
43 : : #include "moab/WriteUtilIface.hpp"
44 : : #include "exodus_order.h"
45 : :
46 : : #ifndef MOAB_HAVE_NETCDF
47 : : #error Attempt to compile WriteNCDF with NetCDF support disabled
48 : : #endif
49 : :
50 : : namespace moab
51 : : {
52 : :
53 : : const int TIME_STR_LEN = 11;
54 : :
55 : : #define INS_ID( stringvar, prefix, id ) sprintf( stringvar, prefix, id )
56 : :
57 : : #define GET_DIM( ncdim, name, val ) \
58 : : { \
59 : : int gdfail = nc_inq_dimid( ncFile, name, &ncdim ); \
60 : : if( NC_NOERR == gdfail ) \
61 : : { \
62 : : size_t tmp_val; \
63 : : gdfail = nc_inq_dimlen( ncFile, ncdim, &tmp_val ); \
64 : : if( NC_NOERR != gdfail ) { MB_SET_ERR( MB_FAILURE, "WriteNCDF:: couldn't get dimension length" ); } \
65 : : else \
66 : : val = tmp_val; \
67 : : } \
68 : : else \
69 : : val = 0; \
70 : : }
71 : :
72 : : #define GET_DIMB( ncdim, name, varname, id, val ) \
73 : : INS_ID( name, varname, id ); \
74 : : GET_DIM( ncdim, name, val );
75 : :
76 : : #define GET_VAR( name, id, dims ) \
77 : : { \
78 : : id = -1; \
79 : : int gvfail = nc_inq_varid( ncFile, name, &id ); \
80 : : if( NC_NOERR == gvfail ) \
81 : : { \
82 : : int ndims; \
83 : : gvfail = nc_inq_varndims( ncFile, id, &ndims ); \
84 : : if( NC_NOERR == gvfail ) \
85 : : { \
86 : : dims.resize( ndims ); \
87 : : gvfail = nc_inq_vardimid( ncFile, id, &dims[0] ); \
88 : : } \
89 : : } \
90 : : }
91 : :
92 : 9 : WriterIface* WriteNCDF::factory( Interface* iface )
93 : : {
94 [ + - ]: 9 : return new WriteNCDF( iface );
95 : : }
96 : :
97 : 9 : WriteNCDF::WriteNCDF( Interface* impl )
98 [ + - ]: 9 : : mdbImpl( impl ), ncFile( 0 ), mCurrentMeshHandle( 0 ), mGeomDimensionTag( 0 ), repeat_face_blocks( 0 )
99 : : {
100 [ - + ]: 9 : assert( impl != NULL );
101 : :
102 [ + - ]: 9 : impl->query_interface( mWriteIface );
103 : :
104 : : // Initialize in case tag_get_handle fails below
105 : : //! Get and cache predefined tag handles
106 : 9 : int negone = -1;
107 : : impl->tag_get_handle( MATERIAL_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mMaterialSetTag, MB_TAG_SPARSE | MB_TAG_CREAT,
108 [ + - ]: 9 : &negone );
109 : :
110 : : impl->tag_get_handle( DIRICHLET_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mDirichletSetTag, MB_TAG_SPARSE | MB_TAG_CREAT,
111 [ + - ]: 9 : &negone );
112 : :
113 : : impl->tag_get_handle( NEUMANN_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mNeumannSetTag, MB_TAG_SPARSE | MB_TAG_CREAT,
114 [ + - ]: 9 : &negone );
115 : :
116 [ + - ]: 9 : mGlobalIdTag = impl->globalId_tag();
117 : :
118 : 9 : int dum_val_array[] = { -1, -1, -1, -1 };
119 : : impl->tag_get_handle( HAS_MID_NODES_TAG_NAME, 4, MB_TYPE_INTEGER, mHasMidNodesTag, MB_TAG_SPARSE | MB_TAG_CREAT,
120 [ + - ]: 9 : dum_val_array );
121 : :
122 : : impl->tag_get_handle( "distFactor", 0, MB_TYPE_DOUBLE, mDistFactorTag,
123 [ + - ]: 9 : MB_TAG_SPARSE | MB_TAG_VARLEN | MB_TAG_CREAT );
124 : :
125 [ + - ]: 9 : impl->tag_get_handle( "qaRecord", 0, MB_TYPE_OPAQUE, mQaRecordTag, MB_TAG_SPARSE | MB_TAG_VARLEN | MB_TAG_CREAT );
126 : :
127 [ + - ]: 9 : impl->tag_get_handle( "WriteNCDF element mark", 1, MB_TYPE_BIT, mEntityMark, MB_TAG_CREAT );
128 : 9 : }
129 : :
130 : 27 : WriteNCDF::~WriteNCDF()
131 : : {
132 : 9 : mdbImpl->release_interface( mWriteIface );
133 : :
134 : 9 : mdbImpl->tag_delete( mEntityMark );
135 : :
136 [ + - ]: 9 : if( 0 != ncFile ) ncFile = 0;
137 [ - + ]: 18 : }
138 : :
139 : 0 : void WriteNCDF::reset_block( std::vector< MaterialSetData >& block_info )
140 : : {
141 : 0 : std::vector< MaterialSetData >::iterator iter;
142 : :
143 [ # # ][ # # ]: 0 : for( iter = block_info.begin(); iter != block_info.end(); ++iter )
[ # # ]
144 : : {
145 [ # # ][ # # ]: 0 : iter->elements.clear();
146 : : }
147 : 0 : }
148 : :
149 : 18 : void WriteNCDF::time_and_date( char* time_string, char* date_string )
150 : : {
151 : : struct tm* local_time;
152 : : time_t calendar_time;
153 : :
154 : 18 : calendar_time = time( NULL );
155 : 18 : local_time = localtime( &calendar_time );
156 : :
157 [ + - ][ - + ]: 18 : assert( NULL != time_string && NULL != date_string );
158 : :
159 : 18 : strftime( time_string, TIME_STR_LEN, "%H:%M:%S", local_time );
160 : 18 : strftime( date_string, TIME_STR_LEN, "%m/%d/%Y", local_time );
161 : :
162 : : // Terminate with NULL character
163 : 18 : time_string[10] = (char)NULL;
164 : 18 : date_string[10] = (char)NULL;
165 : 18 : }
166 : :
167 : 9 : ErrorCode WriteNCDF::write_file( const char* exodus_file_name, const bool overwrite, const FileOptions& opts,
168 : : const EntityHandle* ent_handles, const int num_sets,
169 : : const std::vector< std::string >& qa_records, const Tag*, int, int user_dimension )
170 : : {
171 [ + - ][ + - ]: 9 : assert( 0 != mMaterialSetTag && 0 != mNeumannSetTag && 0 != mDirichletSetTag );
[ - + ]
172 : :
173 [ - + ][ # # ]: 9 : if( user_dimension == 0 ) mdbImpl->get_dimension( user_dimension );
174 : :
175 [ + - ][ - + ]: 9 : if( opts.get_null_option( "REPEAT_FACE_BLOCKS" ) == MB_SUCCESS ) repeat_face_blocks = 1;
176 : :
177 [ + - ][ + - ]: 18 : std::vector< EntityHandle > blocks, nodesets, sidesets, entities;
[ + - ][ + - ]
178 : :
179 : : // Separate into blocks, nodesets, sidesets
180 : :
181 [ + + ]: 9 : if( num_sets == 0 )
182 : : {
183 : : // Default to all defined block, nodeset and sideset-type sets
184 [ + - ]: 7 : Range this_range;
185 [ + - ]: 7 : mdbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, &mMaterialSetTag, NULL, 1, this_range );
186 [ + - ][ + - ]: 7 : std::copy( this_range.begin(), this_range.end(), std::back_inserter( blocks ) );
[ + - ][ + - ]
187 [ + - ]: 7 : this_range.clear();
188 [ + - ]: 7 : mdbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, &mDirichletSetTag, NULL, 1, this_range );
189 [ + - ][ + - ]: 7 : std::copy( this_range.begin(), this_range.end(), std::back_inserter( nodesets ) );
[ + - ][ + - ]
190 [ + - ]: 7 : this_range.clear();
191 [ + - ]: 7 : mdbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, &mNeumannSetTag, NULL, 1, this_range );
192 [ + - ][ + - ]: 7 : std::copy( this_range.begin(), this_range.end(), std::back_inserter( sidesets ) );
[ + - ][ + - ]
193 : :
194 : : // If there is nothing to write, write everything as one block.
195 [ - + ][ # # ]: 7 : if( blocks.empty() && nodesets.empty() && sidesets.empty() )
[ # # ][ - + ]
196 : : {
197 [ # # ]: 0 : this_range.clear();
198 [ # # ][ # # ]: 0 : for( int d = user_dimension; d > 0 && this_range.empty(); --d )
[ # # ][ # # ]
199 [ # # ]: 0 : mdbImpl->get_entities_by_dimension( 0, d, this_range, false );
200 [ # # ][ # # ]: 0 : if( this_range.empty() ) return MB_FILE_WRITE_ERROR;
201 : :
202 : : EntityHandle block_handle;
203 : 0 : int block_id = 1;
204 [ # # ]: 0 : mdbImpl->create_meshset( MESHSET_SET, block_handle );
205 [ # # ]: 0 : mdbImpl->tag_set_data( mMaterialSetTag, &block_handle, 1, &block_id );
206 [ # # ]: 0 : mdbImpl->add_entities( block_handle, this_range );
207 [ # # ][ + - ]: 7 : blocks.push_back( block_handle );
208 : 7 : }
209 : : }
210 : : else
211 : : {
212 : : int dummy;
213 [ + + ]: 4 : for( const EntityHandle* iter = ent_handles; iter < ent_handles + num_sets; ++iter )
214 : : {
215 [ + - ][ + - ]: 2 : if( MB_SUCCESS == mdbImpl->tag_get_data( mMaterialSetTag, &( *iter ), 1, &dummy ) && -1 != dummy )
[ + - ][ + - ]
216 [ + - ]: 2 : blocks.push_back( *iter );
217 [ # # ][ # # ]: 0 : else if( MB_SUCCESS == mdbImpl->tag_get_data( mDirichletSetTag, &( *iter ), 1, &dummy ) && -1 != dummy )
[ # # ][ # # ]
218 [ # # ]: 0 : nodesets.push_back( *iter );
219 [ # # ][ # # ]: 0 : else if( MB_SUCCESS == mdbImpl->tag_get_data( mNeumannSetTag, &( *iter ), 1, &dummy ) && -1 != dummy )
[ # # ][ # # ]
220 [ # # ]: 0 : sidesets.push_back( *iter );
221 : : }
222 : : }
223 : :
224 : : // If there is nothing to write just return.
225 [ - + ][ # # ]: 9 : if( blocks.empty() && nodesets.empty() && sidesets.empty() ) return MB_FILE_WRITE_ERROR;
[ # # ][ - + ]
226 : :
227 : : // Try to get mesh information
228 [ + - ]: 18 : ExodusMeshInfo mesh_info;
229 : :
230 [ + - ]: 18 : std::vector< MaterialSetData > block_info;
231 [ + - ]: 18 : std::vector< NeumannSetData > sideset_info;
232 [ + - ]: 18 : std::vector< DirichletSetData > nodeset_info;
233 : :
234 : 9 : mesh_info.num_dim = user_dimension;
235 : :
236 [ + - ]: 9 : if( qa_records.empty() )
237 : : {
238 : : // qa records are empty - initialize some MB-standard ones
239 [ + - ][ + - ]: 9 : mesh_info.qaRecords.push_back( "MB" );
240 [ + - ][ + - ]: 9 : mesh_info.qaRecords.push_back( "0.99" );
241 : : char string1[80], string2[80];
242 [ + - ]: 9 : time_and_date( string2, string1 );
243 [ + - ][ + - ]: 9 : mesh_info.qaRecords.push_back( string2 );
244 [ + - ][ + - ]: 9 : mesh_info.qaRecords.push_back( string1 );
245 : : }
246 : : else
247 : : {
248 : : // Constrained to multiples of 4 qa records
249 [ # # ]: 0 : assert( qa_records.size() % 4 == 0 );
250 : :
251 [ # # ][ # # ]: 0 : std::copy( qa_records.begin(), qa_records.end(), std::back_inserter( mesh_info.qaRecords ) );
252 : : }
253 : :
254 : 9 : block_info.clear();
255 [ + - ][ - + ]: 9 : if( gather_mesh_information( mesh_info, block_info, sideset_info, nodeset_info, blocks, sidesets, nodesets ) !=
256 : : MB_SUCCESS )
257 : : {
258 [ # # ]: 0 : reset_block( block_info );
259 : 0 : return MB_FAILURE;
260 : : }
261 : :
262 : : // Try to open the file after gather mesh info succeeds
263 [ + - ][ + - ]: 9 : int fail = nc_create( exodus_file_name, overwrite ? NC_CLOBBER : NC_NOCLOBBER, &ncFile );
264 [ - + ]: 9 : if( NC_NOERR != fail )
265 : : {
266 [ # # ]: 0 : reset_block( block_info );
267 : 0 : return MB_FAILURE;
268 : : }
269 : :
270 [ + - ][ - + ]: 9 : if( write_header( mesh_info, block_info, sideset_info, nodeset_info, exodus_file_name ) != MB_SUCCESS )
271 : : {
272 [ # # ]: 0 : reset_block( block_info );
273 : 0 : return MB_FAILURE;
274 : : }
275 : :
276 : : {
277 : : // write dummy time_whole
278 : 9 : double timev = 0.0; // dummy, to make paraview happy
279 : 9 : size_t start = 0, count = 1;
280 : : int nc_var;
281 [ + - ]: 9 : std::vector< int > dims;
282 [ + - ][ + - ]: 9 : GET_VAR( "time_whole", nc_var, dims );
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
283 [ + - ]: 9 : fail = nc_put_vara_double( ncFile, nc_var, &start, &count, &timev );
284 [ - + ][ # # ]: 9 : if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "Failed writing dist factor array" ); }
[ # # ][ # # ]
[ # # ][ # # ]
[ + - ]
285 : : }
286 : :
287 [ + - ][ - + ]: 9 : if( write_nodes( mesh_info.num_nodes, mesh_info.nodes, mesh_info.num_dim ) != MB_SUCCESS )
288 : : {
289 [ # # ]: 0 : reset_block( block_info );
290 : 0 : return MB_FAILURE;
291 : : }
292 : :
293 [ + - ][ - + ]: 9 : if( !mesh_info.polyhedronFaces.empty() )
294 : : {
295 [ # # ][ # # ]: 0 : if( write_poly_faces( mesh_info ) != MB_SUCCESS )
296 : : {
297 [ # # ]: 0 : reset_block( block_info );
298 : 0 : return MB_FAILURE;
299 : : }
300 : : }
301 : :
302 [ + - ][ - + ]: 9 : if( write_elementblocks( mesh_info, block_info ) )
303 : : {
304 [ # # ]: 0 : reset_block( block_info );
305 : 0 : return MB_FAILURE;
306 : : }
307 : :
308 : : // Write the three maps
309 [ + - ][ - + ]: 9 : if( write_global_node_order_map( mesh_info.num_nodes, mesh_info.nodes ) != MB_SUCCESS )
310 : : {
311 [ # # ]: 0 : reset_block( block_info );
312 : 0 : return MB_FAILURE;
313 : : }
314 : :
315 [ + - ][ - + ]: 9 : if( write_global_element_order_map( mesh_info.num_elements ) != MB_SUCCESS )
316 : : {
317 [ # # ]: 0 : reset_block( block_info );
318 : 0 : return MB_FAILURE;
319 : : }
320 : :
321 [ + - ][ - + ]: 9 : if( write_element_order_map( mesh_info.num_elements ) != MB_SUCCESS )
322 : : {
323 [ # # ]: 0 : reset_block( block_info );
324 : 0 : return MB_FAILURE;
325 : : }
326 : :
327 : : /*
328 : : if (write_elementmap(mesh_info) != MB_SUCCESS)
329 : : return MB_FAILURE;
330 : : */
331 : :
332 [ + - ][ - + ]: 9 : if( write_BCs( sideset_info, nodeset_info ) != MB_SUCCESS )
333 : : {
334 [ # # ]: 0 : reset_block( block_info );
335 : 0 : return MB_FAILURE;
336 : : }
337 : :
338 [ + - ][ - + ]: 9 : if( write_qa_records( mesh_info.qaRecords ) != MB_SUCCESS ) return MB_FAILURE;
339 : :
340 : : // Copy the qa records into the argument
341 : : // mesh_info.qaRecords.swap(qa_records);
342 : : // Close the file
343 [ + - ]: 9 : fail = nc_close( ncFile );
344 [ - + ][ # # ]: 9 : if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "Trouble closing file" ); }
[ # # ][ # # ]
[ # # ][ # # ]
345 : :
346 : 18 : return MB_SUCCESS;
347 : : }
348 : :
349 : 9 : ErrorCode WriteNCDF::gather_mesh_information( ExodusMeshInfo& mesh_info, std::vector< MaterialSetData >& block_info,
350 : : std::vector< NeumannSetData >& sideset_info,
351 : : std::vector< DirichletSetData >& nodeset_info,
352 : : std::vector< EntityHandle >& blocks,
353 : : std::vector< EntityHandle >& sidesets,
354 : : std::vector< EntityHandle >& nodesets )
355 : : {
356 : : ErrorCode rval;
357 : 9 : std::vector< EntityHandle >::iterator vector_iter, end_vector_iter;
358 : :
359 : 9 : mesh_info.num_nodes = 0;
360 : 9 : mesh_info.num_elements = 0;
361 : 9 : mesh_info.num_elementblocks = 0;
362 : 9 : mesh_info.num_polyhedra_blocks = 0;
363 : :
364 : 9 : int id = 0;
365 : :
366 : 9 : vector_iter = blocks.begin();
367 : 9 : end_vector_iter = blocks.end();
368 : :
369 [ + - ]: 9 : std::vector< EntityHandle > parent_meshsets;
370 : :
371 : : // Clean out the bits for the element mark
372 [ + - ]: 9 : rval = mdbImpl->tag_delete( mEntityMark );
373 [ - + ]: 9 : if( MB_SUCCESS != rval ) return rval;
374 [ + - ]: 9 : rval = mdbImpl->tag_get_handle( "WriteNCDF element mark", 1, MB_TYPE_BIT, mEntityMark, MB_TAG_CREAT );
375 [ - + ]: 9 : if( MB_SUCCESS != rval ) return rval;
376 : :
377 : 9 : int highest_dimension_of_element_blocks = 0;
378 : :
379 [ + - ][ + - ]: 34 : for( vector_iter = blocks.begin(); vector_iter != blocks.end(); ++vector_iter )
[ + + ]
380 : : {
381 [ + - ]: 25 : MaterialSetData block_data;
382 : :
383 : : // For the purpose of qa records, get the parents of these blocks
384 [ + - ][ + - ]: 25 : if( mdbImpl->get_parent_meshsets( *vector_iter, parent_meshsets ) != MB_SUCCESS ) return MB_FAILURE;
[ - + ]
385 : :
386 : : // Get all Entity Handles in the mesh set
387 [ + - ]: 50 : Range dummy_range;
[ + - - ]
388 [ + - ][ + - ]: 25 : rval = mdbImpl->get_entities_by_handle( *vector_iter, dummy_range, true );
389 [ - + ]: 25 : if( MB_SUCCESS != rval ) return rval;
390 : :
391 : : // Skip empty blocks
392 [ + - ][ - + ]: 25 : if( dummy_range.empty() ) continue;
393 : :
394 : : // Get the block's id
395 [ + - ][ + - ]: 25 : if( mdbImpl->tag_get_data( mMaterialSetTag, &( *vector_iter ), 1, &id ) != MB_SUCCESS )
[ - + ]
396 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "Couldn't get block id from a tag for an element block" ); }
[ # # ][ # # ]
[ # # ]
397 : :
398 : 25 : block_data.id = id;
399 : 25 : block_data.number_attributes = 0;
400 : :
401 : : // Wait a minute, we are doing some filtering here that doesn't make sense at this level CJS
402 : :
403 : : // Find the dimension of the last entity in this range
404 [ + - ][ + - ]: 25 : int this_dim = CN::Dimension( TYPE_FROM_HANDLE( dummy_range.back() ) );
[ + - ]
405 [ - + ][ # # ]: 25 : if( this_dim > 3 ) { MB_SET_ERR( MB_TYPE_OUT_OF_RANGE, "Block " << id << " contains entity sets" ); }
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
406 [ + - ][ + - ]: 25 : block_data.elements = dummy_range.subset_by_dimension( this_dim );
407 : :
408 : : // End of -- wait a minute, we are doing some filtering here that doesn't make sense at this
409 : : // level CJS
410 : :
411 : : // Get the entity type for this block, verifying that it's the same for all elements
412 [ + - ][ + - ]: 25 : EntityType entity_type = TYPE_FROM_HANDLE( block_data.elements.front() );
413 [ + - ][ - + ]: 25 : if( !block_data.elements.all_of_type( entity_type ) )
414 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "Entities in block " << id << " not of common type" ); }
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
415 : :
416 : 25 : int dimension = -1;
417 [ + + ][ + + ]: 25 : if( entity_type == MBQUAD || entity_type == MBTRI )
418 : 7 : dimension = 2; // Output shells by default
419 [ + + ]: 18 : else if( entity_type == MBEDGE )
420 : 3 : dimension = 1;
421 : : else
422 [ + - ]: 15 : dimension = CN::Dimension( entity_type );
423 : :
424 [ + + ]: 25 : if( dimension > highest_dimension_of_element_blocks ) highest_dimension_of_element_blocks = dimension;
425 : :
426 [ + - ]: 50 : std::vector< EntityHandle > tmp_conn;
[ + - - ]
427 [ + - ][ + - ]: 25 : rval = mdbImpl->get_connectivity( &( block_data.elements.front() ), 1, tmp_conn );
428 [ - + ]: 25 : if( MB_SUCCESS != rval ) return rval;
429 [ + - ]: 25 : block_data.element_type = ExoIIUtil::get_element_type_from_num_verts( tmp_conn.size(), entity_type, dimension );
430 : :
431 [ - + ]: 25 : if( block_data.element_type == EXOII_MAX_ELEM_TYPE )
432 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "Element type in block " << id << " didn't get set correctly" ); }
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
433 : :
434 [ - + ]: 25 : if( block_data.element_type == EXOII_POLYGON )
435 : : {
436 : : // get all poly connectivity
437 : 0 : int numconn = 0;
438 [ # # ][ # # ]: 0 : for( Range::iterator eit = block_data.elements.begin(); eit != block_data.elements.end(); eit++ )
[ # # ][ # # ]
[ # # ]
439 : : {
440 [ # # ]: 0 : EntityHandle polg = *eit;
441 : 0 : int nnodes = 0;
442 : 0 : const EntityHandle* conn = NULL;
443 [ # # ][ # # ]: 0 : rval = mdbImpl->get_connectivity( polg, conn, nnodes );MB_CHK_ERR( rval );
[ # # ][ # # ]
444 : 0 : numconn += nnodes;
445 : : }
446 : 0 : block_data.number_nodes_per_element = numconn;
447 : : }
448 : : else
449 : 25 : block_data.number_nodes_per_element = ExoIIUtil::VerticesPerElement[block_data.element_type];
450 : :
451 : : // Number of nodes for this block
452 [ + - ]: 25 : block_data.number_elements = block_data.elements.size();
453 : :
454 : : // Total number of elements
455 : 25 : mesh_info.num_elements += block_data.number_elements;
456 : :
457 : : // Get the nodes for the elements
458 [ + - ]: 25 : rval = mWriteIface->gather_nodes_from_elements( block_data.elements, mEntityMark, mesh_info.nodes );
459 [ - + ]: 25 : if( MB_SUCCESS != rval ) return rval;
460 : :
461 : : // if polyhedra block
462 [ - + ]: 25 : if( EXOII_POLYHEDRON == block_data.element_type )
463 : : {
464 [ # # ][ # # ]: 0 : rval = mdbImpl->get_connectivity( block_data.elements, mesh_info.polyhedronFaces );MB_CHK_ERR( rval );
[ # # ][ # # ]
465 : 0 : mesh_info.num_polyhedra_blocks++;
466 : : }
467 : :
468 [ + + ]: 25 : if( !sidesets.empty() )
469 : : {
470 : : // If there are sidesets, keep track of which elements are being written out
471 [ + - ][ + - ]: 147 : for( Range::iterator iter = block_data.elements.begin(); iter != block_data.elements.end(); ++iter )
[ + - ][ + - ]
[ + + ]
472 : : {
473 : 132 : unsigned char bit = 0x1;
474 [ + - ][ + - ]: 132 : rval = mdbImpl->tag_set_data( mEntityMark, &( *iter ), 1, &bit );
475 [ - + ]: 132 : if( MB_SUCCESS != rval ) return rval;
476 : : }
477 : : }
478 : :
479 [ + - ]: 25 : block_info.push_back( block_data );
480 : :
481 : 25 : const void* data = NULL;
482 : 25 : int size = 0;
483 [ + - ][ + - ]: 25 : if( MB_SUCCESS == mdbImpl->tag_get_by_ptr( mQaRecordTag, &( *vector_iter ), 1, &data, &size ) && NULL != data )
[ - + ][ # # ]
[ - + ]
484 : : {
485 : : // There are qa records on this block - copy over to mesh qa records
486 : 0 : const char* qa_rec = static_cast< const char* >( data );
487 : 0 : int start = 0;
488 : 0 : int count = 0;
489 [ # # ]: 0 : for( int i = 0; i < size; i++ )
490 : : {
491 [ # # ]: 0 : if( qa_rec[i] == '\0' )
492 : : {
493 [ # # ]: 0 : std::string qa_string( &qa_rec[start], i - start );
494 [ # # ]: 0 : mesh_info.qaRecords.push_back( qa_string );
495 : 0 : start = i + 1;
496 : 0 : count++;
497 : : }
498 : : }
499 : :
500 : : // Constrained to multiples of 4 qa records
501 [ # # ][ # # ]: 25 : if( count > 0 ) assert( count % 4 == 0 );
[ + - ]
502 : : }
503 : 25 : }
504 : :
505 : 9 : mesh_info.num_elementblocks = block_info.size();
506 : :
507 [ + - ][ + - ]: 34 : for( std::vector< MaterialSetData >::iterator blit = block_info.begin(); blit != block_info.end(); blit++ )
[ + + ]
508 : : {
509 [ + - ]: 25 : MaterialSetData& block = *blit;
510 [ + - ]: 25 : if( block.element_type != EXOII_POLYHEDRON )
511 [ + - ][ + - ]: 25 : { mesh_info.polyhedronFaces = subtract( mesh_info.polyhedronFaces, block.elements ); }
512 : : }
513 : :
514 : : // If user hasn't entered dimension, we figure it out
515 [ - + ]: 9 : if( mesh_info.num_dim == 0 )
516 : : {
517 : : // Never want 1 or zero dimensions
518 [ # # ]: 0 : if( highest_dimension_of_element_blocks < 2 )
519 : 0 : mesh_info.num_dim = 3;
520 : : else
521 : 0 : mesh_info.num_dim = highest_dimension_of_element_blocks;
522 : : }
523 : :
524 [ + - ][ + - ]: 9 : Range::iterator range_iter, end_range_iter;
525 [ + - ]: 9 : range_iter = mesh_info.nodes.begin();
526 [ + - ]: 9 : end_range_iter = mesh_info.nodes.end();
527 : :
528 [ + - ]: 9 : mesh_info.num_nodes = mesh_info.nodes.size();
529 : :
530 : : //------nodesets--------
531 : :
532 : 9 : vector_iter = nodesets.begin();
533 : 9 : end_vector_iter = nodesets.end();
534 : :
535 [ + - ][ + - ]: 18 : for( ; vector_iter != end_vector_iter; ++vector_iter )
[ + + ]
536 : : {
537 [ + - ]: 9 : DirichletSetData nodeset_data;
538 : 9 : nodeset_data.id = 0;
539 : 9 : nodeset_data.number_nodes = 0;
540 : :
541 : : // Get the nodeset's id
542 [ + - ][ + - ]: 9 : if( mdbImpl->tag_get_data( mDirichletSetTag, &( *vector_iter ), 1, &id ) != MB_SUCCESS )
[ - + ]
543 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "Couldn't get id tag for nodeset " << id ); }
[ # # ][ # # ]
[ # # ][ # # ]
544 : :
545 : 9 : nodeset_data.id = id;
546 : :
547 [ + - ][ + - ]: 18 : std::vector< EntityHandle > node_vector;
548 : : // Get the nodes of the nodeset that are in mesh_info.nodes
549 [ + - ][ + - ]: 9 : if( mdbImpl->get_entities_by_handle( *vector_iter, node_vector, true ) != MB_SUCCESS )
[ - + ]
550 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "Couldn't get nodes in nodeset " << id ); }
[ # # ][ # # ]
[ # # ][ # # ]
551 : :
552 : : // Get the tag for distribution factors
553 : : const double* dist_factor_vector;
554 : : int dist_factor_size;
555 : 9 : const void* ptr = 0;
556 : :
557 : 9 : int has_dist_factors = 0;
558 [ + - ][ + - ]: 9 : if( mdbImpl->tag_get_by_ptr( mDistFactorTag, &( *vector_iter ), 1, &ptr, &dist_factor_size ) == MB_SUCCESS &&
[ - + ][ # # ]
[ - + ]
559 : : dist_factor_size )
560 : 0 : has_dist_factors = 1;
561 : 9 : dist_factor_size /= sizeof( double );
562 : 9 : dist_factor_vector = reinterpret_cast< const double* >( ptr );
563 : 9 : std::vector< EntityHandle >::iterator iter, end_iter;
564 : 9 : iter = node_vector.begin();
565 : 9 : end_iter = node_vector.end();
566 : :
567 : 9 : int j = 0;
568 : 9 : unsigned char node_marked = 0;
569 : : ErrorCode result;
570 [ + - ][ + - ]: 111 : for( ; iter != end_iter; ++iter )
[ + + ]
571 : : {
572 [ + - ][ + - ]: 102 : if( TYPE_FROM_HANDLE( *iter ) != MBVERTEX ) continue;
[ - + ]
573 [ + - ][ + - ]: 102 : result = mdbImpl->tag_get_data( mEntityMark, &( *iter ), 1, &node_marked );MB_CHK_SET_ERR( result, "Couldn't get mark data" );
[ - + ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
574 : :
575 [ + - ]: 102 : if( 0x1 == node_marked )
576 : : {
577 [ + - ][ + - ]: 102 : nodeset_data.nodes.push_back( *iter );
578 [ - + ]: 102 : if( 0 != has_dist_factors )
579 [ # # ]: 0 : nodeset_data.node_dist_factors.push_back( dist_factor_vector[j] );
580 : : else
581 [ + - ]: 102 : nodeset_data.node_dist_factors.push_back( 1.0 );
582 : : }
583 : 102 : j++;
584 : : }
585 : :
586 : 9 : nodeset_data.number_nodes = nodeset_data.nodes.size();
587 [ + - ][ + - ]: 9 : nodeset_info.push_back( nodeset_data );
588 : 9 : }
589 : :
590 : : //------sidesets--------
591 : 9 : vector_iter = sidesets.begin();
592 : 9 : end_vector_iter = sidesets.end();
593 : :
594 [ + - ][ + - ]: 21 : for( ; vector_iter != end_vector_iter; ++vector_iter )
[ + + ]
595 : : {
596 [ + - ]: 12 : NeumannSetData sideset_data;
597 : :
598 : : // Get the sideset's id
599 [ + - ][ + - ]: 12 : if( mdbImpl->tag_get_data( mNeumannSetTag, &( *vector_iter ), 1, &id ) != MB_SUCCESS ) return MB_FAILURE;
[ - + ]
600 : :
601 : 12 : sideset_data.id = id;
602 [ + - ]: 12 : sideset_data.mesh_set_handle = *vector_iter;
603 : :
604 : : // Get the sides in two lists, one forward the other reverse; starts with forward sense
605 : : // by convention
606 [ + - ][ + - ]: 24 : Range forward_elems, reverse_elems;
[ + - ][ + - ]
607 [ + - ][ + - ]: 12 : if( get_sideset_elems( *vector_iter, 0, forward_elems, reverse_elems ) == MB_FAILURE ) return MB_FAILURE;
[ - + ]
608 : :
609 [ + - ][ - + ]: 12 : ErrorCode result = get_valid_sides( forward_elems, mesh_info, 1, sideset_data );MB_CHK_SET_ERR( result, "Couldn't get valid sides data" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
610 [ + - ][ - + ]: 12 : result = get_valid_sides( reverse_elems, mesh_info, -1, sideset_data );MB_CHK_SET_ERR( result, "Couldn't get valid sides data" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
611 : :
612 : 12 : sideset_data.number_elements = sideset_data.elements.size();
613 [ + - ][ + - ]: 12 : sideset_info.push_back( sideset_data );
614 : 12 : }
615 : :
616 : 9 : return MB_SUCCESS;
617 : : }
618 : :
619 : 24 : ErrorCode WriteNCDF::get_valid_sides( Range& elems, ExodusMeshInfo& /*mesh_info*/, const int sense,
620 : : NeumannSetData& sideset_data )
621 : : {
622 : : // This is where we see if underlying element of side set element is included in output
623 : :
624 : : // Get the sideset-based info for distribution factors
625 : 24 : const double* dist_factor_vector = 0;
626 : 24 : int dist_factor_size = 0;
627 : :
628 : : // Initialize dist_fac_iter to get rid of compiler warning
629 : 24 : const double* dist_fac_iter = 0;
630 : 24 : const void* ptr = 0;
631 : 24 : bool has_dist_factors = false;
632 [ + - ][ + - ]: 48 : if( mdbImpl->tag_get_by_ptr( mDistFactorTag, &( sideset_data.mesh_set_handle ), 1, &ptr, &dist_factor_size ) ==
633 [ + - ][ + - ]: 24 : MB_SUCCESS &&
634 : : dist_factor_size )
635 : : {
636 : 24 : has_dist_factors = true;
637 : 24 : dist_factor_vector = reinterpret_cast< const double* >( ptr );
638 : 24 : dist_fac_iter = dist_factor_vector;
639 : 24 : dist_factor_size /= sizeof( double );
640 : : }
641 : :
642 : 24 : unsigned char element_marked = 0;
643 : : ErrorCode result;
644 [ + - ][ + - ]: 192 : for( Range::iterator iter = elems.begin(); iter != elems.end(); ++iter )
[ + - ][ + - ]
[ + + ]
645 : : {
646 : : // Should insert here if "side" is a quad/tri on a quad/tri mesh
647 [ + - ][ + - ]: 168 : result = mdbImpl->tag_get_data( mEntityMark, &( *iter ), 1, &element_marked );MB_CHK_SET_ERR( result, "Couldn't get mark data" );
[ - + ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
648 : :
649 [ + + ]: 168 : if( 0x1 == element_marked )
650 : : {
651 [ + - ][ + - ]: 24 : sideset_data.elements.push_back( *iter );
652 : :
653 : : // TJT TODO: the sense should really be # edges + 1or2
654 [ + + ][ + - ]: 24 : sideset_data.side_numbers.push_back( ( sense == 1 ? 1 : 2 ) );
655 : : }
656 : : else
657 : : { // Then "side" is probably a quad/tri on a hex/tet mesh
658 [ + - ]: 144 : std::vector< EntityHandle > parents;
659 [ + - ][ + - ]: 144 : int dimension = CN::Dimension( TYPE_FROM_HANDLE( *iter ) );
[ + - ]
660 : :
661 : : // Get the adjacent parent element of "side"
662 [ + - ][ + - ]: 144 : if( mdbImpl->get_adjacencies( &( *iter ), 1, dimension + 1, false, parents ) != MB_SUCCESS )
663 : : {
664 : : #if 0
665 : : // This is not treated as an error, print warning messages for
666 : : // debugging only
667 : : fprintf(stderr, "[Warning]: Couldn't get adjacencies for sideset.\n");
668 : : #endif
669 : : }
670 : :
671 [ + - ]: 144 : if( !parents.empty() )
672 : : {
673 : : // Make sure the adjacent parent element will be output
674 [ + + ][ + - ]: 294 : for( unsigned int k = 0; k < parents.size(); k++ )
675 : : {
676 [ + - ][ + - ]: 150 : result = mdbImpl->tag_get_data( mEntityMark, &( parents[k] ), 1, &element_marked );MB_CHK_SET_ERR( result, "Couldn't get mark data" );
[ - + ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
677 : :
678 : : int side_no, this_sense, this_offset;
679 [ + - ][ + + ]: 294 : if( 0x1 == element_marked &&
680 [ + + ][ + - ]: 294 : mdbImpl->side_number( parents[k], *iter, side_no, this_sense, this_offset ) == MB_SUCCESS &&
[ + - ][ + - ]
[ + + ]
681 : 144 : this_sense == sense )
682 : : {
683 [ + - ][ + - ]: 72 : sideset_data.elements.push_back( parents[k] );
684 [ + - ]: 72 : sideset_data.side_numbers.push_back( side_no + 1 );
685 : 72 : break;
686 : : }
687 : : }
688 : : }
689 : : else
690 : : {
691 : : #if 0
692 : : // This is not treated as an error, print warning messages for
693 : : // debugging only
694 : : fprintf(stderr, "[Warning]: No parent element exists for element in sideset %i\n", sideset_data.id);
695 : : #endif
696 : 144 : }
697 : : }
698 : :
699 [ + - ]: 168 : if( sideset_data.elements.size() != 0 )
700 : : {
701 : : // Distribution factors
702 [ + - ][ + - ]: 168 : int num_nodes = CN::VerticesPerEntity( TYPE_FROM_HANDLE( *iter ) );
[ + - ]
703 : : // put some dummy dist factors for polygons; why do we need them?
704 [ + - ][ + - ]: 168 : if( TYPE_FROM_HANDLE( *iter ) == MBPOLYGON ) num_nodes = 1; // dummy
[ - + ]
705 [ + - ]: 168 : if( has_dist_factors )
706 : : {
707 : : std::copy( dist_fac_iter, dist_fac_iter + num_nodes,
708 [ + - ][ + - ]: 168 : std::back_inserter( sideset_data.ss_dist_factors ) );
709 : 168 : dist_fac_iter += num_nodes;
710 : : }
711 : : else
712 : : {
713 [ # # ]: 168 : for( int j = 0; j < num_nodes; j++ )
714 [ # # ]: 0 : sideset_data.ss_dist_factors.push_back( 1.0 );
715 : : }
716 : : }
717 : : }
718 : :
719 : 24 : return MB_SUCCESS;
720 : : }
721 : :
722 : 9 : ErrorCode WriteNCDF::write_qa_records( std::vector< std::string >& qa_record_list )
723 : : {
724 : 9 : int i = 0;
725 : :
726 [ + - ][ + + ]: 18 : for( std::vector< std::string >::iterator string_it = qa_record_list.begin(); string_it != qa_record_list.end(); )
727 : : {
728 [ + + ]: 45 : for( int j = 0; j < 4; j++ )
729 [ + - ][ + - ]: 36 : write_qa_string( ( *string_it++ ).c_str(), i, j );
[ + - ]
730 : 9 : i++;
731 : : }
732 : :
733 : 9 : return MB_SUCCESS;
734 : : }
735 : :
736 : 36 : ErrorCode WriteNCDF::write_qa_string( const char* string, int record_number, int record_position )
737 : : {
738 : : // Get the variable id in the exodus file
739 : :
740 [ + - ]: 36 : std::vector< int > dims;
741 : 36 : int temp_var = -1;
742 [ + - ][ + - ]: 36 : GET_VAR( "qa_records", temp_var, dims );
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
743 [ - + ][ # # ]: 36 : if( -1 == temp_var ) { MB_SET_ERR( MB_FAILURE, "WriteNCDF:: Problem getting qa record variable" ); }
[ # # ][ # # ]
[ # # ][ # # ]
744 : : size_t count[3], start[3];
745 : :
746 : : // Write out the record
747 : 36 : start[0] = record_number;
748 : 36 : start[1] = record_position;
749 : 36 : start[2] = 0;
750 : :
751 : 36 : count[0] = 1;
752 : 36 : count[1] = 1;
753 : 36 : count[2] = (long)strlen( string ) + 1;
754 [ + - ]: 36 : int fail = nc_put_vara_text( ncFile, temp_var, start, count, string );
755 [ - + ][ # # ]: 36 : if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "Failed to position qa string variable" ); }
[ # # ][ # # ]
[ # # ][ # # ]
756 : :
757 : 36 : return MB_SUCCESS;
758 : : }
759 : :
760 : 9 : ErrorCode WriteNCDF::write_nodes( int num_nodes, Range& nodes, int dimension )
761 : : {
762 : : // Write coordinates names
763 : 9 : int nc_var = -1;
764 [ + - ]: 9 : std::vector< int > dims;
765 [ + - ][ + - ]: 9 : GET_VAR( "coor_names", nc_var, dims );
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
766 [ - + ][ # # ]: 9 : if( -1 == nc_var ) { MB_SET_ERR( MB_FAILURE, "Trouble getting coordinate name variable" ); }
[ # # ][ # # ]
[ # # ][ # # ]
767 : :
768 : 9 : size_t start[2] = { 0, 0 }, count[2] = { 1, ExoIIInterface::MAX_STR_LENGTH };
769 : : char dum_str[ExoIIInterface::MAX_STR_LENGTH];
770 : 9 : strcpy( dum_str, "x" );
771 [ + - ]: 9 : int fail = nc_put_vara_text( ncFile, nc_var, start, count, dum_str );
772 [ - + ]: 9 : if( NC_NOERR != fail )
773 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "Trouble adding x coordinate name; netcdf message: " << nc_strerror( fail ) ); }
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
774 : :
775 : 9 : start[0] = 1;
776 : 9 : strcpy( dum_str, "y" );
777 [ + - ]: 9 : fail = nc_put_vara_text( ncFile, nc_var, start, count, dum_str );
778 [ - + ]: 9 : if( NC_NOERR != fail )
779 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "Trouble adding y coordinate name; netcdf message: " << nc_strerror( fail ) ); }
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
780 : :
781 : 9 : start[0] = 2;
782 : 9 : strcpy( dum_str, "z" );
783 [ + - ]: 9 : fail = nc_put_vara_text( ncFile, nc_var, start, count, dum_str );
784 [ - + ]: 9 : if( NC_NOERR != fail )
785 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "Trouble adding z coordinate name; netcdf message: " << nc_strerror( fail ) ); }
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
786 : :
787 : : // See if should transform coordinates
788 : : ErrorCode result;
789 : : Tag trans_tag;
790 [ + - ]: 9 : result = mdbImpl->tag_get_handle( MESH_TRANSFORM_TAG_NAME, 16, MB_TYPE_DOUBLE, trans_tag );
791 : 9 : bool transform_needed = true;
792 [ + - ]: 9 : if( result == MB_TAG_NOT_FOUND ) transform_needed = false;
793 : :
794 [ - + ]: 9 : int num_coords_to_fill = transform_needed ? 3 : dimension;
795 : :
796 [ + - ]: 18 : std::vector< double* > coord_arrays( 3 );
797 [ + - ][ + - ]: 9 : coord_arrays[0] = new double[num_nodes];
[ + - ]
798 [ + - ][ + - ]: 9 : coord_arrays[1] = new double[num_nodes];
[ + - ]
799 [ + - ]: 9 : coord_arrays[2] = NULL;
800 : :
801 [ + - ][ + - ]: 9 : if( num_coords_to_fill == 3 ) coord_arrays[2] = new double[num_nodes];
[ + - ][ + - ]
802 : :
803 [ + - ]: 9 : result = mWriteIface->get_node_coords( dimension, num_nodes, nodes, mGlobalIdTag, 1, coord_arrays );
804 [ - + ]: 9 : if( result != MB_SUCCESS )
805 : : {
806 [ # # ][ # # ]: 0 : delete[] coord_arrays[0];
807 [ # # ][ # # ]: 0 : delete[] coord_arrays[1];
808 [ # # ][ # # ]: 0 : if( coord_arrays[2] ) delete[] coord_arrays[2];
[ # # ][ # # ]
809 : 0 : return result;
810 : : }
811 : :
812 [ - + ]: 9 : if( transform_needed )
813 : : {
814 : : double trans_matrix[16];
815 : 0 : const EntityHandle mesh = 0;
816 [ # # ][ # # ]: 0 : result = mdbImpl->tag_get_data( trans_tag, &mesh, 0, trans_matrix );MB_CHK_SET_ERR( result, "Couldn't get transform data" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
817 : :
818 [ # # ]: 0 : for( int i = 0; i < num_nodes; i++ )
819 : : {
820 : : double vec1[3];
821 : : double vec2[3];
822 : :
823 [ # # ]: 0 : vec2[0] = coord_arrays[0][i];
824 [ # # ]: 0 : vec2[1] = coord_arrays[1][i];
825 [ # # ]: 0 : vec2[2] = coord_arrays[2][i];
826 : :
827 [ # # ]: 0 : for( int row = 0; row < 3; row++ )
828 : : {
829 : 0 : vec1[row] = 0.0;
830 [ # # ]: 0 : for( int col = 0; col < 3; col++ )
831 : 0 : vec1[row] += ( trans_matrix[( row * 4 ) + col] * vec2[col] );
832 : : }
833 : :
834 [ # # ]: 0 : coord_arrays[0][i] = vec1[0];
835 [ # # ]: 0 : coord_arrays[1][i] = vec1[1];
836 [ # # ]: 0 : coord_arrays[2][i] = vec1[2];
837 : : }
838 : : }
839 : :
840 : : // Write the nodes
841 : 9 : nc_var = -1;
842 [ + - ][ + - ]: 9 : GET_VAR( "coord", nc_var, dims );
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
843 [ - + ][ # # ]: 9 : if( -1 == nc_var ) { MB_SET_ERR( MB_FAILURE, "Trouble getting coordinate variable" ); }
[ # # ][ # # ]
[ # # ][ # # ]
844 : 9 : start[0] = 0;
845 : 9 : count[1] = num_nodes;
846 [ + - ][ + - ]: 9 : fail = nc_put_vara_double( ncFile, nc_var, start, count, &( coord_arrays[0][0] ) );
847 [ - + ][ # # ]: 9 : if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "Trouble writing x coordinate" ); }
[ # # ][ # # ]
[ # # ][ # # ]
848 : :
849 : 9 : start[0] = 1;
850 [ + - ][ + - ]: 9 : fail = nc_put_vara_double( ncFile, nc_var, start, count, &( coord_arrays[1][0] ) );
851 [ - + ][ # # ]: 9 : if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "Trouble writing y coordinate" ); }
[ # # ][ # # ]
[ # # ][ # # ]
852 : :
853 : 9 : start[0] = 2;
854 [ + - ][ + - ]: 9 : fail = nc_put_vara_double( ncFile, nc_var, start, count, &( coord_arrays[2][0] ) );
855 [ - + ][ # # ]: 9 : if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "Trouble writing z coordinate" ); }
[ # # ][ # # ]
[ # # ][ # # ]
856 : :
857 [ + - ][ + - ]: 9 : delete[] coord_arrays[0];
858 [ + - ][ + - ]: 9 : delete[] coord_arrays[1];
859 [ + - ][ + - ]: 9 : if( coord_arrays[2] ) delete[] coord_arrays[2];
[ + - ][ + - ]
860 : :
861 : 18 : return MB_SUCCESS;
862 : : }
863 : :
864 : 0 : ErrorCode WriteNCDF::write_poly_faces( ExodusMeshInfo& mesh_info )
865 : : {
866 : : // write all polygons that are not in another element block;
867 : : // usually they are nowhere else, but be sure, write in this block only ones that are not in the
868 : : // other blocks
869 [ # # ]: 0 : Range pfaces = mesh_info.polyhedronFaces;
870 : :
871 : : /*
872 : : * int fbconn1(num_nod_per_fa1) ;
873 : : fbconn1:elem_type = "nsided" ;
874 : : int fbepecnt1(num_fa_in_blk1) ;
875 : : fbepecnt1:entity_type1 = "NODE" ;
876 : : fbepecnt1:entity_type2 = "FACE" ;
877 : : */
878 [ # # ][ # # ]: 0 : if( pfaces.empty() ) return MB_SUCCESS;
879 : : char wname[80];
880 : 0 : int nc_var = -1;
881 [ # # ]: 0 : std::vector< int > dims;
882 : :
883 : : // write one for each element block, to make paraview and visit happy
884 [ # # ]: 0 : int num_faces_in_block = (int)pfaces.size();
885 [ # # ]: 0 : for( unsigned int bl = 0; bl < mesh_info.num_polyhedra_blocks; bl++ )
886 : : {
887 : 0 : INS_ID( wname, "fbconn%u", bl + 1 ); // it is the first block
888 [ # # ][ # # ]: 0 : GET_VAR( wname, nc_var, dims ); // fbconn# variable, 1 dimensional
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
889 : :
890 : 0 : INS_ID( wname, "num_nod_per_fa%u", bl + 1 );
891 : : int ncdim, num_nod_per_face;
892 [ # # ][ # # ]: 0 : GET_DIM( ncdim, wname, num_nod_per_face );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
893 [ # # ][ # # ]: 0 : int* connectivity = new int[num_nod_per_face];
894 : 0 : int ixcon = 0, j = 0;
895 [ # # ]: 0 : std::vector< int > fbepe( num_faces_in_block ); // fbepecnt1
896 [ # # ][ # # ]: 0 : for( Range::iterator eit = pfaces.begin(); eit != pfaces.end(); eit++ )
[ # # ][ # # ]
[ # # ]
897 : : {
898 [ # # ]: 0 : EntityHandle polyg = *eit;
899 : 0 : int nnodes = 0;
900 : 0 : const EntityHandle* conn = NULL;
901 [ # # ][ # # ]: 0 : ErrorCode rval = mdbImpl->get_connectivity( polyg, conn, nnodes );MB_CHK_ERR( rval );
[ # # ][ # # ]
902 [ # # ]: 0 : for( int k = 0; k < nnodes; k++ )
903 : 0 : connectivity[ixcon++] = conn[k];
904 [ # # ]: 0 : fbepe[j++] = nnodes;
905 : : }
906 : 0 : size_t start[1] = { 0 }, count[1] = { 0 };
907 : 0 : count[0] = ixcon;
908 [ # # ]: 0 : int fail = nc_put_vara_int( ncFile, nc_var, start, count, connectivity );
909 [ # # ]: 0 : if( NC_NOERR != fail )
910 : : {
911 [ # # ]: 0 : delete[] connectivity;
912 [ # # ][ # # ]: 0 : MB_SET_ERR( MB_FAILURE, "Couldn't write fbconn variable" );
[ # # ][ # # ]
[ # # ]
913 : : }
914 : :
915 : 0 : INS_ID( wname, "fbepecnt%u", bl + 1 );
916 [ # # ][ # # ]: 0 : GET_VAR( wname, nc_var, dims ); // fbconn# variable, 1 dimensional
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
917 : 0 : count[0] = num_faces_in_block;
918 : :
919 [ # # ][ # # ]: 0 : fail = nc_put_vara_int( ncFile, nc_var, start, count, &fbepe[0] );
920 [ # # ]: 0 : if( NC_NOERR != fail )
921 : : {
922 [ # # ]: 0 : delete[] connectivity;
923 [ # # ][ # # ]: 0 : MB_SET_ERR( MB_FAILURE, "Couldn't write fbepecnt variable" );
[ # # ][ # # ]
[ # # ]
924 : : }
925 : :
926 : 0 : int id = bl + 1;
927 [ # # ][ # # ]: 0 : if( write_exodus_integer_variable( "fa_prop1", &id, bl, 1 ) != MB_SUCCESS )
928 [ # # ][ # # ]: 0 : { MB_SET_ERR_CONT( "Problem writing element block id " << id ); }
[ # # ][ # # ]
[ # # ][ # # ]
929 : :
930 : 0 : int status = 1;
931 [ # # ][ # # ]: 0 : if( write_exodus_integer_variable( "fa_status", &status, bl, 1 ) != MB_SUCCESS )
932 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "Problem writing face block status" ); }
[ # # ][ # # ]
[ # # ]
933 : :
934 [ # # ]: 0 : delete[] connectivity;
935 [ # # ]: 0 : if( 0 == repeat_face_blocks ) break; // do not repeat face blocks
[ # # # ]
936 : 0 : }
937 : :
938 : 0 : return MB_SUCCESS;
939 : : }
940 : 9 : ErrorCode WriteNCDF::write_header( ExodusMeshInfo& mesh_info, std::vector< MaterialSetData >& block_info,
941 : : std::vector< NeumannSetData >& sideset_info,
942 : : std::vector< DirichletSetData >& nodeset_info, const char* filename )
943 : : {
944 : : // Get the date and time
945 : : char time[TIME_STR_LEN];
946 : : char date[TIME_STR_LEN];
947 [ + - ]: 9 : time_and_date( time, date );
948 : :
949 [ + - ]: 9 : std::string title_string = "MOAB";
950 [ + - ]: 9 : title_string.append( "(" );
951 [ + - ]: 9 : title_string.append( filename );
952 [ + - ]: 9 : title_string.append( "): " );
953 [ + - ]: 9 : title_string.append( date );
954 [ + - ]: 9 : title_string.append( ": " );
955 [ + - ]: 9 : title_string.append( "time " );
956 : :
957 [ - + ]: 9 : if( title_string.length() > ExoIIInterface::MAX_LINE_LENGTH )
958 [ # # ]: 0 : title_string.resize( ExoIIInterface::MAX_LINE_LENGTH );
959 : :
960 : : // Initialize the exodus file
961 : :
962 [ + - ]: 9 : int result = initialize_exodus_file( mesh_info, block_info, sideset_info, nodeset_info, title_string.c_str() );
963 : :
964 [ - + ]: 9 : if( result == MB_FAILURE ) return MB_FAILURE;
965 : :
966 : 9 : return MB_SUCCESS;
967 : : }
968 : :
969 : 9 : ErrorCode WriteNCDF::write_elementblocks( ExodusMeshInfo& mesh_info, std::vector< MaterialSetData >& block_data )
970 : : {
971 : : unsigned int i;
972 : 9 : int block_index = 0; // Index into block list, may != 1 if there are inactive blocks
973 : 9 : int exodus_id = 1;
974 : :
975 [ + + ]: 34 : for( i = 0; i < block_data.size(); i++ )
976 : : {
977 [ + - ]: 25 : MaterialSetData& block = block_data[i];
978 : :
979 : 25 : unsigned int num_nodes_per_elem = block.number_nodes_per_element;
980 : :
981 : : // Write out the id for the block
982 : :
983 : 25 : int id = block.id;
984 : 25 : int num_values = 1;
985 : :
986 [ + - ][ - + ]: 25 : if( write_exodus_integer_variable( "eb_prop1", &id, block_index, num_values ) != MB_SUCCESS )
987 [ # # ][ # # ]: 0 : { MB_SET_ERR_CONT( "Problem writing element block id " << id ); }
[ # # ][ # # ]
[ # # ][ # # ]
988 : :
989 : : // Write out the block status
990 : :
991 : 25 : int status = 1;
992 [ - + ][ # # ]: 25 : if( 0 == block.number_elements ) { MB_SET_ERR( MB_FAILURE, "No elements in block " << id ); }
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
993 : :
994 [ + - ][ - + ]: 25 : if( write_exodus_integer_variable( "eb_status", &status, block_index, num_values ) != MB_SUCCESS )
995 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "Problem writing element block status" ); }
[ # # ][ # # ]
[ # # ]
996 : :
997 : : //
998 : : // Map the connectivity to the new nodes
999 : 25 : const unsigned int num_elem = block.number_elements;
1000 : 25 : unsigned int num_nodes = num_nodes_per_elem * num_elem;
1001 [ + - ][ - + ]: 25 : if( EXOII_POLYGON == block.element_type || EXOII_POLYHEDRON == block.element_type )
1002 : 0 : { num_nodes = num_nodes_per_elem; }
1003 [ + - ][ + - ]: 25 : int* connectivity = new int[num_nodes];
1004 : :
1005 : 25 : ErrorCode result = MB_SUCCESS;
1006 [ + - ]: 25 : if( block.element_type != EXOII_POLYHEDRON )
1007 : : mWriteIface->get_element_connect( num_elem, num_nodes_per_elem, mGlobalIdTag, block.elements, mGlobalIdTag,
1008 [ + - ]: 25 : exodus_id, connectivity );
1009 [ - + ]: 25 : if( result != MB_SUCCESS )
1010 : : {
1011 [ # # ]: 0 : delete[] connectivity;
1012 [ # # ][ # # ]: 0 : MB_SET_ERR( result, "Couldn't get element array to write from" );
[ # # ][ # # ]
[ # # ]
1013 : : }
1014 : :
1015 : : // If necessary, convert from EXODUS to CN node order
1016 : 25 : const EntityType elem_type = ExoIIUtil::ExoIIElementMBEntity[block.element_type];
1017 [ + - ][ - + ]: 25 : assert( block.elements.all_of_type( elem_type ) );
1018 : 25 : const int* reorder = 0;
1019 [ + - ][ + - ]: 25 : if( block.element_type != EXOII_POLYHEDRON && block.element_type != EXOII_POLYGON )
1020 : 25 : reorder = exodus_elem_order_map[elem_type][block.number_nodes_per_element];
1021 [ + + ]: 25 : if( reorder )
1022 [ + - ]: 4 : WriteUtilIface::reorder( reorder, connectivity, block.number_elements, block.number_nodes_per_element );
1023 : :
1024 : : char wname[80];
1025 : 25 : int nc_var = -1;
1026 [ + - ]: 25 : std::vector< int > dims;
1027 [ + - ]: 25 : if( block.element_type != EXOII_POLYHEDRON )
1028 : : {
1029 : 25 : exodus_id += num_elem;
1030 : 25 : INS_ID( wname, "connect%u", i + 1 );
1031 : :
1032 [ + - ][ + - ]: 25 : GET_VAR( wname, nc_var, dims );
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
1033 [ - + ]: 25 : if( -1 == nc_var )
1034 : : {
1035 [ # # ]: 0 : delete[] connectivity;
1036 [ # # ][ # # ]: 25 : MB_SET_ERR( MB_FAILURE, "Couldn't get connectivity variable" );
[ # # ][ # # ]
[ # # ]
1037 : : }
1038 : : }
1039 : :
1040 [ - + ]: 25 : if( EXOII_POLYGON == block.element_type )
1041 : : {
1042 : 0 : size_t start[1] = { 0 }, count[1] = { num_nodes_per_elem };
1043 [ # # ]: 0 : int fail = nc_put_vara_int( ncFile, nc_var, start, count, connectivity );
1044 [ # # ]: 0 : if( NC_NOERR != fail )
1045 : : {
1046 [ # # ]: 0 : delete[] connectivity;
1047 [ # # ][ # # ]: 0 : MB_SET_ERR( MB_FAILURE, "Couldn't write connectivity variable" );
[ # # ][ # # ]
[ # # ]
1048 : : }
1049 : : // now put also number ebepecnt1
1050 : 0 : INS_ID( wname, "ebepecnt%u", i + 1 );
1051 [ # # ][ # # ]: 0 : GET_VAR( wname, nc_var, dims );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1052 : 0 : count[0] = block.number_elements;
1053 : 0 : start[0] = 0;
1054 : : // reuse connectivity array, to not allocate another one
1055 : 0 : int j = 0;
1056 [ # # ][ # # ]: 0 : for( Range::iterator eit = block.elements.begin(); eit != block.elements.end(); j++, eit++ )
[ # # ][ # # ]
[ # # ]
1057 : : {
1058 [ # # ]: 0 : EntityHandle polg = *eit;
1059 : 0 : int nnodes = 0;
1060 : 0 : const EntityHandle* conn = NULL;
1061 [ # # ][ # # ]: 0 : ErrorCode rval = mdbImpl->get_connectivity( polg, conn, nnodes );MB_CHK_ERR( rval );
[ # # ][ # # ]
1062 : 0 : connectivity[j] = nnodes;
1063 : : }
1064 [ # # ]: 0 : fail = nc_put_vara_int( ncFile, nc_var, start, count, connectivity );
1065 [ # # ]: 0 : if( NC_NOERR != fail )
1066 : : {
1067 [ # # ]: 0 : delete[] connectivity;
1068 [ # # ][ # # ]: 0 : MB_SET_ERR( MB_FAILURE, "Couldn't write ebepecnt variable" );
[ # # ][ # # ]
[ # # ]
1069 : : }
1070 : : }
1071 [ + - ]: 25 : else if( block.element_type != EXOII_POLYHEDRON )
1072 : : {
1073 : 25 : size_t start[2] = { 0, 0 }, count[2] = { num_elem, num_nodes_per_elem };
1074 [ + - ]: 25 : int fail = nc_put_vara_int( ncFile, nc_var, start, count, connectivity );
1075 [ - + ]: 25 : if( NC_NOERR != fail )
1076 : : {
1077 [ # # ]: 0 : delete[] connectivity;
1078 [ # # ][ # # ]: 25 : MB_SET_ERR( MB_FAILURE, "Couldn't write connectivity variable" );
[ # # ][ # # ]
[ # # ]
1079 : : }
1080 : : }
1081 : : else // if (block.element_type == EXOII_POLYHEDRON)
1082 : : {
1083 : : /* write a lot of stuff // faconn
1084 : : num_fa_in_blk1 = 15 ;
1085 : : num_nod_per_fa1 = 58 ;
1086 : : num_el_in_blk1 = 3 ;
1087 : : num_fac_per_el1 = 17 ;
1088 : : int fbconn1(num_nod_per_fa1) ;
1089 : : fbconn1:elem_type = "nsided" ;
1090 : : int fbepecnt1(num_fa_in_blk1) ;
1091 : : fbepecnt1:entity_type1 = "NODE" ;
1092 : : fbepecnt1:entity_type2 = "FACE" ;
1093 : : int facconn1(num_fac_per_el1) ;
1094 : : facconn1:elem_type = "NFACED" ;
1095 : : int ebepecnt1(num_el_in_blk1) ;
1096 : : ebepecnt1:entity_type1 = "FACE" ;
1097 : : ebepecnt1:entity_type2 = "ELEM" ;
1098 : :
1099 : : */
1100 : 0 : Range& block_faces = mesh_info.polyhedronFaces;
1101 : : // ErrorCode rval = mdbImpl->get_connectivity(block.elements, block_faces);
1102 : : // MB_CHK_ERR(rval);
1103 : :
1104 : : // reuse now connectivity for facconn1
1105 : 0 : INS_ID( wname, "facconn%u", i + 1 );
1106 [ # # ][ # # ]: 0 : GET_VAR( wname, nc_var, dims ); // fbconn# variable, 1 dimensional
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1107 : :
1108 [ # # ][ # # ]: 0 : std::vector< int > ebepe( block.elements.size() ); // ebepecnt1
1109 : 0 : int ixcon = 0, j = 0;
1110 : 0 : size_t start[1] = { 0 }, count[1] = { 0 };
1111 : :
1112 [ # # ][ # # ]: 0 : for( Range::iterator eit = block.elements.begin(); eit != block.elements.end(); eit++ )
[ # # ][ # # ]
[ # # ]
1113 : : {
1114 [ # # ]: 0 : EntityHandle polyh = *eit;
1115 : 0 : int nfaces = 0;
1116 : 0 : const EntityHandle* conn = NULL;
1117 [ # # ][ # # ]: 0 : ErrorCode rval = mdbImpl->get_connectivity( polyh, conn, nfaces );MB_CHK_ERR( rval );
[ # # ][ # # ]
1118 [ # # ]: 0 : for( int k = 0; k < nfaces; k++ )
1119 : : {
1120 [ # # ]: 0 : int index = block_faces.index( conn[k] );
1121 [ # # ][ # # ]: 0 : if( index == -1 ) MB_SET_ERR( MB_FAILURE, "Couldn't find face in polyhedron" );
[ # # ][ # # ]
[ # # ][ # # ]
1122 : 0 : connectivity[ixcon++] = index + 1;
1123 : : }
1124 [ # # ]: 0 : ebepe[j++] = nfaces;
1125 : : // num_faces+=nfaces;
1126 : : }
1127 : 0 : count[0] = ixcon; // facconn1
1128 [ # # ]: 0 : int fail = nc_put_vara_int( ncFile, nc_var, start, count, connectivity );
1129 [ # # ]: 0 : if( NC_NOERR != fail )
1130 : : {
1131 [ # # ]: 0 : delete[] connectivity;
1132 [ # # ][ # # ]: 0 : MB_SET_ERR( MB_FAILURE, "Couldn't write fbconn variable" );
[ # # ][ # # ]
[ # # ]
1133 : : }
1134 : :
1135 : 0 : INS_ID( wname, "ebepecnt%u", i + 1 );
1136 [ # # ][ # # ]: 0 : GET_VAR( wname, nc_var, dims ); // ebepecnt# variable, 1 dimensional
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1137 [ # # ]: 0 : count[0] = block.elements.size();
1138 : :
1139 [ # # ][ # # ]: 0 : fail = nc_put_vara_int( ncFile, nc_var, start, count, &ebepe[0] );
1140 [ # # ]: 0 : if( NC_NOERR != fail )
1141 : : {
1142 [ # # ]: 0 : delete[] connectivity;
1143 [ # # ][ # # ]: 0 : MB_SET_ERR( MB_FAILURE, "Couldn't write fbepecnt variable" );
[ # # ][ # # ]
[ # # ][ # # ]
1144 : 0 : }
1145 : : }
1146 : 25 : block_index++;
1147 [ + - ][ + - ]: 25 : delete[] connectivity;
1148 : 25 : }
1149 : :
1150 : 9 : return MB_SUCCESS;
1151 : : }
1152 : :
1153 : 9 : ErrorCode WriteNCDF::write_global_node_order_map( int num_nodes, Range& nodes )
1154 : : {
1155 : : // Note: this routine bypasses the standard exodusII interface for efficiency!
1156 : :
1157 : : // Node order map
1158 [ + - ][ + - ]: 9 : int* map = new int[num_nodes];
1159 : :
1160 : : // For now, output a dummy map!
1161 : :
1162 [ + - ][ + - ]: 9 : Range::iterator range_iter, end_iter;
1163 [ + - ]: 9 : range_iter = nodes.begin();
1164 [ + - ]: 9 : end_iter = nodes.end();
1165 : :
1166 : 9 : int i = 0;
1167 : :
1168 [ + - ][ + - ]: 2345551 : for( ; range_iter != end_iter; ++range_iter )
[ + + ]
1169 : : {
1170 : : // TODO -- do we really want to cast this to an int?
1171 [ + - ][ + - ]: 2345542 : map[i++] = (int)ID_FROM_HANDLE( *range_iter );
1172 : : }
1173 : :
1174 : : // Output array and cleanup
1175 : :
1176 [ + - ]: 9 : int error = write_exodus_integer_variable( "node_num_map", map, 0, num_nodes );
1177 : :
1178 [ + - ][ + - ]: 9 : if( map ) delete[] map;
1179 : :
1180 [ - + ][ # # ]: 9 : if( error < 0 ) { MB_SET_ERR( MB_FAILURE, "Failed writing global node order map" ); }
[ # # ][ # # ]
[ # # ][ # # ]
1181 : :
1182 : 9 : return MB_SUCCESS;
1183 : : }
1184 : :
1185 : 9 : ErrorCode WriteNCDF::write_global_element_order_map( int num_elements )
1186 : : {
1187 : : // Allocate map array
1188 [ + - ]: 9 : int* map = new int[num_elements];
1189 : :
1190 : : // Many Sandia codes assume this map is unique, and CUBIT does not currently
1191 : : // have unique ids for all elements. Therefore, to make sure nothing crashes,
1192 : : // insert a dummy map...
1193 : :
1194 [ + + ]: 512249 : for( int i = 0; i < num_elements; i++ )
1195 : 512240 : map[i] = i + 1;
1196 : :
1197 : : // Output array and cleanup
1198 : :
1199 : 9 : int error = write_exodus_integer_variable( "elem_num_map", map, 0, num_elements );
1200 : :
1201 [ + - ][ + - ]: 9 : if( map ) delete[] map;
1202 : :
1203 [ - + ][ # # ]: 9 : if( error < 0 ) { MB_SET_ERR( MB_FAILURE, "Failed writing global element order map" ); }
[ # # ][ # # ]
[ # # ][ # # ]
1204 : :
1205 : 9 : return MB_SUCCESS;
1206 : : }
1207 : :
1208 : 9 : ErrorCode WriteNCDF::write_element_order_map( int num_elements )
1209 : : {
1210 : : // Note: this routine bypasses the standard exodusII interface for efficiency!
1211 : :
1212 : : // Element order map
1213 [ + - ]: 9 : int* map = new int[num_elements];
1214 : :
1215 : : // For now, output a dummy map!
1216 : :
1217 [ + + ]: 512249 : for( int i = 0; i < num_elements; i++ )
1218 : : {
1219 : 512240 : map[i] = i + 1;
1220 : : }
1221 : :
1222 : : // Output array and cleanup
1223 : :
1224 : 9 : int error = write_exodus_integer_variable( "elem_map", map, 0, num_elements );
1225 : :
1226 [ + - ][ + - ]: 9 : if( map ) delete[] map;
1227 : :
1228 [ - + ][ # # ]: 9 : if( error < 0 ) { MB_SET_ERR( MB_FAILURE, "Failed writing element map" ); }
[ # # ][ # # ]
[ # # ][ # # ]
1229 : :
1230 : 9 : return MB_SUCCESS;
1231 : : }
1232 : :
1233 : 119 : ErrorCode WriteNCDF::write_exodus_integer_variable( const char* variable_name, int* variable_array, int start_position,
1234 : : int number_values )
1235 : : {
1236 : : // Note: this routine bypasses the standard exodusII interface for efficiency!
1237 : :
1238 : : // Write directly to netcdf interface for efficiency
1239 : :
1240 : : // Get the variable id of the element map
1241 : 119 : int nc_var = -1;
1242 [ + - ]: 119 : std::vector< int > dims;
1243 [ + - ][ + - ]: 119 : GET_VAR( variable_name, nc_var, dims );
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
1244 [ - + ]: 119 : if( -1 == nc_var )
1245 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to locate variable " << variable_name << " in file" ); }
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1246 : : // This contortion is necessary because netCDF is expecting nclongs;
1247 : : // fortunately it's necessary only when ints and nclongs aren't the same size
1248 : :
1249 : : size_t start[1], count[1];
1250 : 119 : start[0] = start_position;
1251 : 119 : count[0] = number_values;
1252 : :
1253 : 119 : int fail = NC_NOERR;
1254 : : if( sizeof( int ) == sizeof( long ) ) { fail = nc_put_vara_int( ncFile, nc_var, start, count, variable_array ); }
1255 : : else
1256 : : {
1257 [ + - ][ + - ]: 119 : long* lptr = new long[number_values];
1258 [ + + ]: 3370233 : for( int jj = 0; jj < number_values; jj++ )
1259 : 3370114 : lptr[jj] = variable_array[jj];
1260 [ + - ]: 119 : fail = nc_put_vara_long( ncFile, nc_var, start, count, lptr );
1261 [ + - ]: 119 : delete[] lptr;
1262 : : }
1263 : :
1264 [ - + ][ # # ]: 119 : if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "Failed to store variable " << variable_name ); }
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1265 : :
1266 : 119 : return MB_SUCCESS;
1267 : : }
1268 : :
1269 : 9 : ErrorCode WriteNCDF::write_BCs( std::vector< NeumannSetData >& sidesets, std::vector< DirichletSetData >& nodesets )
1270 : : {
1271 : : unsigned int i, j;
1272 : : int id;
1273 : 9 : int ns_index = -1;
1274 : :
1275 [ + - ][ + - ]: 18 : for( std::vector< DirichletSetData >::iterator ns_it = nodesets.begin(); ns_it != nodesets.end(); ++ns_it )
[ + + ]
1276 : : {
1277 : : // Get number of nodes in set
1278 [ + - ]: 9 : int number_nodes = ( *ns_it ).number_nodes;
1279 [ - + ]: 9 : if( 0 == number_nodes ) continue;
1280 : :
1281 : : // If we're here, we have a non-empty nodeset; increment the index
1282 : 9 : ns_index++;
1283 : :
1284 : : // Get the node set id
1285 [ + - ]: 9 : id = ( *ns_it ).id;
1286 : :
1287 : : // Build new array to old exodus ids
1288 [ + - ][ + - ]: 9 : int* exodus_id_array = new int[number_nodes];
1289 [ + - ][ + - ]: 9 : double* dist_factor_array = new double[number_nodes];
1290 : :
1291 : 9 : std::vector< EntityHandle >::iterator begin_iter, end_iter;
1292 : 9 : std::vector< double >::iterator other_iter;
1293 [ + - ]: 9 : begin_iter = ( *ns_it ).nodes.begin();
1294 [ + - ]: 9 : end_iter = ( *ns_it ).nodes.end();
1295 [ + - ]: 9 : other_iter = ( *ns_it ).node_dist_factors.begin();
1296 : :
1297 : 9 : j = 0;
1298 : : int exodus_id;
1299 : : ErrorCode result;
1300 : : // Fill up node array and dist. factor array at the same time
1301 [ + - ][ + - ]: 111 : for( ; begin_iter != end_iter; ++begin_iter )
[ + + ]
1302 : : {
1303 [ + - ][ + - ]: 102 : result = mdbImpl->tag_get_data( mGlobalIdTag, &( *begin_iter ), 1, &exodus_id );MB_CHK_SET_ERR( result, "Problem getting id tag data" );
[ - + ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
1304 : :
1305 : 102 : exodus_id_array[j] = exodus_id;
1306 [ + - ]: 102 : dist_factor_array[j] = *( other_iter );
1307 [ + - ]: 102 : ++other_iter;
1308 : 102 : j++;
1309 : : }
1310 : :
1311 : : // Write out the id for the nodeset
1312 : :
1313 : 9 : int num_values = 1;
1314 : :
1315 [ + - ][ - + ]: 9 : result = write_exodus_integer_variable( "ns_prop1", &id, ns_index, num_values );MB_CHK_SET_ERR_RET_VAL( result, "Problem writing node set id " << id, MB_FAILURE );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
1316 : :
1317 : : // Write out the nodeset status
1318 : :
1319 : 9 : int status = 1;
1320 [ - + ]: 9 : if( !number_nodes ) status = 0;
1321 : :
1322 [ + - ][ - + ]: 9 : result = write_exodus_integer_variable( "ns_status", &status, ns_index, num_values );MB_CHK_SET_ERR_RET_VAL( result, "Problem writing node set status", MB_FAILURE );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1323 : :
1324 : : // Write it out
1325 : : char wname[80];
1326 : 9 : int nc_var = -1;
1327 [ + - ]: 9 : std::vector< int > dims;
1328 : 9 : INS_ID( wname, "node_ns%d", ns_index + 1 );
1329 [ + - ][ + - ]: 9 : GET_VAR( wname, nc_var, dims );
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
1330 [ - + ][ # # ]: 9 : if( -1 == nc_var ) { MB_SET_ERR( MB_FAILURE, "Failed to get node_ns variable" ); }
[ # # ][ # # ]
[ # # ][ # # ]
1331 : :
1332 : 9 : size_t start = 0, count = number_nodes;
1333 [ + - ]: 9 : int fail = nc_put_vara_int( ncFile, nc_var, &start, &count, exodus_id_array );
1334 [ - + ][ # # ]: 9 : if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "Failed writing exodus id array" ); }
[ # # ][ # # ]
[ # # ][ # # ]
1335 : :
1336 : : // Write out nodeset distribution factors
1337 : 9 : INS_ID( wname, "dist_fact_ns%d", ns_index + 1 );
1338 : 9 : nc_var = -1;
1339 [ + - ][ + - ]: 9 : GET_VAR( wname, nc_var, dims );
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
1340 [ - + ][ # # ]: 9 : if( -1 == nc_var ) { MB_SET_ERR( MB_FAILURE, "Failed to get dist_fact variable" ); }
[ # # ][ # # ]
[ # # ][ # # ]
1341 [ + - ]: 9 : fail = nc_put_vara_double( ncFile, nc_var, &start, &count, dist_factor_array );
1342 [ - + ][ # # ]: 9 : if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "Failed writing dist factor array" ); }
[ # # ][ # # ]
[ # # ][ # # ]
1343 : :
1344 [ + - ]: 9 : delete[] dist_factor_array;
1345 [ + - ][ + - ]: 9 : delete[] exodus_id_array;
1346 : 9 : }
1347 : :
1348 : : // Now do sidesets
1349 : 9 : int ss_index = 0; // Index of sideset - not the same as 'i' because
1350 : : // only writing non-empty side sets
1351 [ + + ]: 21 : for( i = 0; i < sidesets.size(); i++ )
1352 : : {
1353 [ + - ][ + - ]: 12 : NeumannSetData sideset_data = sidesets[i];
1354 : :
1355 : : // Get the side set id
1356 : 12 : int side_set_id = sideset_data.id;
1357 : :
1358 : : // Get number of elements in set
1359 : 12 : int number_elements = sideset_data.number_elements;
1360 [ - + ]: 12 : if( 0 == number_elements ) continue;
1361 : :
1362 : : // Build new array to old exodus ids
1363 [ + - ][ + - ]: 12 : int* output_element_ids = new int[number_elements];
1364 [ + - ][ + - ]: 12 : int* output_element_side_numbers = new int[number_elements];
1365 : :
1366 : 12 : std::vector< EntityHandle >::iterator begin_iter, end_iter;
1367 : 12 : begin_iter = sideset_data.elements.begin();
1368 : 12 : end_iter = sideset_data.elements.end();
1369 : 12 : std::vector< int >::iterator side_iter = sideset_data.side_numbers.begin();
1370 : :
1371 : : // Get the tag handle
1372 : 12 : j = 0;
1373 : : int exodus_id;
1374 : :
1375 : : // For each "side"
1376 [ + - ][ + - ]: 108 : for( ; begin_iter != end_iter; ++begin_iter, ++side_iter )
[ + - ][ + + ]
1377 : : {
1378 [ + - ][ + - ]: 96 : ErrorCode result = mdbImpl->tag_get_data( mGlobalIdTag, &( *begin_iter ), 1, &exodus_id );MB_CHK_SET_ERR( result, "Problem getting exodus id for sideset element "
[ - + ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1379 : : << (long unsigned int)ID_FROM_HANDLE( *begin_iter ) );
1380 : :
1381 : 96 : output_element_ids[j] = exodus_id;
1382 [ + - ]: 96 : output_element_side_numbers[j++] = *side_iter;
1383 : : }
1384 : :
1385 [ + - ]: 12 : if( 0 != number_elements )
1386 : : {
1387 : : // Write out the id for the nodeset
1388 : :
1389 : 12 : int num_values = 1;
1390 : :
1391 : : // ss_prop1[ss_index] = side_set_id
1392 [ + - ][ - + ]: 12 : ErrorCode result = write_exodus_integer_variable( "ss_prop1", &side_set_id, ss_index, num_values );MB_CHK_SET_ERR_RET_VAL( result, "Problem writing node set id " << id, MB_FAILURE );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
1393 : :
1394 : : // FIXME : Something seems wrong here. The we are within a block
1395 : : // started with if (0 != number_elements), so this condition is always
1396 : : // false. This code seems to indicate that we want to write all
1397 : : // sidesets, not just empty ones. But the code both here and in
1398 : : // initialize_exodus_file() skip empty side sets.
1399 : : // - j.k. 2007-03-09
1400 : 12 : int status = 1;
1401 [ - + ]: 12 : if( 0 == number_elements ) status = 0;
1402 : :
1403 : : // ss_status[ss_index] = status
1404 [ + - ][ - + ]: 12 : result = write_exodus_integer_variable( "ss_status", &status, ss_index, num_values );MB_CHK_SET_ERR_RET_VAL( result, "Problem writing side set status", MB_FAILURE );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1405 : :
1406 : : // Increment ss_index now because we want a) we need to
1407 : : // increment it somewhere within the if (0 != number_elements)
1408 : : // block and b) the above calls need a zero-based index
1409 : : // while the following use a one-based index.
1410 : 12 : ++ss_index;
1411 : :
1412 : : char wname[80];
1413 : : int nc_var;
1414 [ + - ]: 12 : std::vector< int > dims;
1415 : 12 : INS_ID( wname, "elem_ss%d", ss_index );
1416 [ + - ][ + - ]: 12 : GET_VAR( wname, nc_var, dims );
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
1417 [ - + ][ # # ]: 12 : if( -1 == nc_var ) { MB_SET_ERR( MB_FAILURE, "Failed to get elem_ss variable" ); }
[ # # ][ # # ]
[ # # ][ # # ]
1418 : 12 : size_t start = 0, count = number_elements;
1419 [ + - ]: 12 : int fail = nc_put_vara_int( ncFile, nc_var, &start, &count, output_element_ids );
1420 [ - + ][ # # ]: 12 : if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "Failed writing sideset element array" ); }
[ # # ][ # # ]
[ # # ][ # # ]
1421 : :
1422 : 12 : INS_ID( wname, "side_ss%d", ss_index );
1423 : 12 : nc_var = -1;
1424 [ + - ][ + - ]: 12 : GET_VAR( wname, nc_var, dims );
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
1425 [ - + ][ # # ]: 12 : if( -1 == nc_var ) { MB_SET_ERR( MB_FAILURE, "Failed to get side_ss variable" ); }
[ # # ][ # # ]
[ # # ][ # # ]
1426 [ + - ]: 12 : fail = nc_put_vara_int( ncFile, nc_var, &start, &count, output_element_side_numbers );
1427 [ - + ][ # # ]: 12 : if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "Failed writing sideset side array" ); }
[ # # ][ # # ]
[ # # ][ # # ]
1428 : :
1429 : 12 : INS_ID( wname, "dist_fact_ss%d", ss_index );
1430 : 12 : nc_var = -1;
1431 [ + - ][ + - ]: 12 : GET_VAR( wname, nc_var, dims );
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
1432 [ - + ][ # # ]: 12 : if( -1 == nc_var ) { MB_SET_ERR( MB_FAILURE, "Failed to get sideset dist factors variable" ); }
[ # # ][ # # ]
[ # # ][ # # ]
1433 : 12 : count = sideset_data.ss_dist_factors.size();
1434 [ + - ][ + - ]: 12 : fail = nc_put_vara_double( ncFile, nc_var, &start, &count, &( sideset_data.ss_dist_factors[0] ) );
1435 [ - + ][ # # ]: 12 : if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "Failed writing sideset dist factors array" ); }
[ # # ][ # # ]
[ # # ][ # # ]
[ + - ]
1436 : : }
1437 : :
1438 [ + - ]: 12 : delete[] output_element_ids;
1439 [ + - ]: 12 : delete[] output_element_side_numbers;
[ + - - ]
1440 : 12 : }
1441 : :
1442 : 9 : return MB_SUCCESS;
1443 : : }
1444 : :
1445 : 9 : ErrorCode WriteNCDF::initialize_exodus_file( ExodusMeshInfo& mesh_info, std::vector< MaterialSetData >& block_data,
1446 : : std::vector< NeumannSetData >& sideset_data,
1447 : : std::vector< DirichletSetData >& nodeset_data, const char* title_string,
1448 : : bool write_maps, bool /* write_sideset_distribution_factors */ )
1449 : : {
1450 : : // This routine takes the place of the exodusii routine ex_put_init,
1451 : : // and additionally pre-defines variables such as qa, element blocks,
1452 : : // sidesets and nodesets in a single pass.
1453 : : //
1454 : : // This is necessary because of the way exodusII works. Each time the
1455 : : // netcdf routine endef is called to take the file out of define mode,
1456 : : // the entire file is copied before the new information is added.
1457 : : //
1458 : : // With very large files, this is simply not workable. This routine takes
1459 : : // the definition portions of all applicable exodus routines and puts them
1460 : : // in a single definition, preventing repeated copying of the file.
1461 : : //
1462 : : // Most of the code is copied directly from the applicable exodus routine,
1463 : : // and thus this routine may not seem very consistent in usage or variable
1464 : : // naming, etc.
1465 : :
1466 : : // Perform the initializations
1467 : :
1468 : : int element_block_index;
1469 : :
1470 : : // Inquire on defined string dimension and general dimension for qa
1471 : :
1472 : : int dim_str, dim_four, dim_line, dim_time;
1473 [ + - ][ - + ]: 9 : if( nc_def_dim( ncFile, "len_string", ExoIIInterface::MAX_STR_LENGTH, &dim_str ) != NC_NOERR )
1474 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to get string length in file" ); }
[ # # ][ # # ]
[ # # ]
1475 : :
1476 [ + - ][ - + ]: 9 : if( nc_def_dim( ncFile, "len_line", ExoIIInterface::MAX_STR_LENGTH, &dim_line ) != NC_NOERR )
1477 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to get line length in file" ); }
[ # # ][ # # ]
[ # # ]
1478 : :
1479 [ + - ][ - + ]: 9 : if( nc_def_dim( ncFile, "four", 4, &dim_four ) != NC_NOERR )
1480 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to locate four in file" ); }
[ # # ][ # # ]
[ # # ]
1481 : :
1482 [ + - ][ - + ]: 9 : if( nc_def_dim( ncFile, "time_step", 1, &dim_time ) != NC_NOERR )
1483 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to locate time step in file" ); }
[ # # ][ # # ]
[ # # ]
1484 : : // some whole_time dummy :(
1485 : : int dtime;
1486 [ + - ][ - + ]: 9 : if( NC_NOERR != nc_def_var( ncFile, "time_whole", NC_DOUBLE, 1, &dim_time, &dtime ) )
1487 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define time whole array" ); }
[ # # ][ # # ]
[ # # ]
1488 : :
1489 : : /* Put file into define mode */
1490 : :
1491 : : // It is possible that an NT filename using backslashes is in the title string
1492 : : // this can give fits to unix codes where the backslash is assumed to be an escape
1493 : : // sequence. For the exodus file, just flip the backslash to a slash to prevent
1494 : : // this problem
1495 : :
1496 : : // Get a working copy of the title_string;
1497 : :
1498 : : char working_title[80];
1499 : 9 : strncpy( working_title, title_string, 79 );
1500 : :
1501 : 9 : int length = strlen( working_title );
1502 [ + + ]: 422 : for( int pos = 0; pos < length; pos++ )
1503 : : {
1504 [ - + ]: 413 : if( working_title[pos] == '\\' ) strncpy( &working_title[pos], "/", 1 );
1505 : : }
1506 : :
1507 [ + - ][ - + ]: 9 : if( NC_NOERR != nc_put_att_text( ncFile, NC_GLOBAL, "title", length, working_title ) )
1508 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define title attribute" ); }
[ # # ][ # # ]
[ # # ]
1509 : :
1510 : : // Add other attributes while we're at it
1511 : 9 : float dum_vers = 6.28F;
1512 [ + - ][ - + ]: 9 : if( NC_NOERR != nc_put_att_float( ncFile, NC_GLOBAL, "api_version", NC_FLOAT, 1, &dum_vers ) )
1513 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define api_version attribute" ); }
[ # # ][ # # ]
[ # # ]
1514 : 9 : dum_vers = 6.28F;
1515 [ + - ][ - + ]: 9 : if( NC_NOERR != nc_put_att_float( ncFile, NC_GLOBAL, "version", NC_FLOAT, 1, &dum_vers ) )
1516 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define version attribute" ); }
[ # # ][ # # ]
[ # # ]
1517 : 9 : int dum_siz = sizeof( double );
1518 [ + - ][ - + ]: 9 : if( NC_NOERR != nc_put_att_int( ncFile, NC_GLOBAL, "floating_point_word_size", NC_INT, 1, &dum_siz ) )
1519 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define floating pt word size attribute" ); }
[ # # ][ # # ]
[ # # ]
1520 : :
1521 : : // Set up number of dimensions
1522 : :
1523 : : int num_el_blk, num_elem, num_nodes, num_dim, num_fa_blk, num_faces;
1524 [ + - ][ - + ]: 9 : if( nc_def_dim( ncFile, "num_dim", (size_t)mesh_info.num_dim, &num_dim ) != NC_NOERR )
1525 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of dimensions" ); }
[ # # ][ # # ]
[ # # ]
1526 : :
1527 [ + - ][ - + ]: 9 : if( nc_def_dim( ncFile, "num_nodes", mesh_info.num_nodes, &num_nodes ) != NC_NOERR )
1528 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of nodes" ); }
[ # # ][ # # ]
[ # # ]
1529 : :
1530 : : int num_nod_per_fa; // it is needed for polyhedron only; need to compute it (connectivity of
1531 : : // faces of polyhedra)
1532 [ + - ][ - + ]: 9 : if( mesh_info.polyhedronFaces.size() > 0 )
1533 [ # # ][ # # ]: 0 : if( nc_def_dim( ncFile, "num_faces", (int)mesh_info.polyhedronFaces.size(), &num_faces ) != NC_NOERR )
[ # # ]
1534 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of nodes" ); }
[ # # ][ # # ]
[ # # ]
1535 : :
1536 [ + - ][ - + ]: 9 : if( nc_def_dim( ncFile, "num_elem", mesh_info.num_elements, &num_elem ) != NC_NOERR )
1537 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of elements" ); }
[ # # ][ # # ]
[ # # ]
1538 : :
1539 [ + - ][ - + ]: 9 : if( nc_def_dim( ncFile, "num_el_blk", mesh_info.num_elementblocks, &num_el_blk ) != NC_NOERR )
1540 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of element blocks" ); }
[ # # ][ # # ]
[ # # ]
1541 : :
1542 : : /* ...and some variables */
1543 : :
1544 : : /* Element block id status array */
1545 : 9 : int idstat = -1;
1546 [ + - ][ - + ]: 9 : if( NC_NOERR != nc_def_var( ncFile, "eb_status", NC_LONG, 1, &num_el_blk, &idstat ) )
1547 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define element block status array" ); }
[ # # ][ # # ]
[ # # ]
1548 : :
1549 : : /* Element block id array */
1550 : :
1551 : 9 : int idarr = -1;
1552 [ + - ][ - + ]: 9 : if( NC_NOERR != nc_def_var( ncFile, "eb_prop1", NC_LONG, 1, &num_el_blk, &idarr ) )
1553 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define element block id array" ); }
[ # # ][ # # ]
[ # # ]
1554 : :
1555 : : /* store property name as attribute of property array variable */
1556 [ + - ][ - + ]: 9 : if( NC_NOERR != nc_put_att_text( ncFile, idarr, "name", strlen( "ID" ), "ID" ) )
1557 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to store element block property name ID" ); }
[ # # ][ # # ]
[ # # ]
1558 : :
1559 : : // count how many are polyhedron blocks
1560 : 9 : int num_fa_blocks = 0, num_polyh_blocks = 0;
1561 [ + + ]: 34 : for( unsigned int i = 0; i < block_data.size(); i++ )
1562 : : {
1563 [ + - ]: 25 : MaterialSetData& block = block_data[i];
1564 [ - + ]: 25 : if( EXOII_POLYHEDRON == block.element_type )
1565 : : {
1566 : 0 : num_fa_blocks++;
1567 : 0 : num_polyh_blocks++;
1568 : : }
1569 : : }
1570 [ + - ][ - + ]: 9 : if( 0 == this->repeat_face_blocks && num_fa_blocks > 1 ) num_fa_blocks = 1;
1571 : : char wname[80];
1572 : :
1573 [ - + ]: 9 : if( num_fa_blocks > 0 )
1574 : : {
1575 : : /* face block id status array */
1576 [ # # ][ # # ]: 0 : if( nc_def_dim( ncFile, "num_fa_blk", num_fa_blocks, &num_fa_blk ) != NC_NOERR )
1577 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of face blocks" ); }
[ # # ][ # # ]
[ # # ]
1578 : :
1579 : 0 : int idstatf = -1;
1580 [ # # ][ # # ]: 0 : if( NC_NOERR != nc_def_var( ncFile, "fa_status", NC_LONG, 1, &num_fa_blk, &idstatf ) )
1581 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define face block status array" ); }
[ # # ][ # # ]
[ # # ]
1582 : :
1583 : : /* Element block id array */
1584 : :
1585 : 0 : int idarrf = -1;
1586 [ # # ][ # # ]: 0 : if( NC_NOERR != nc_def_var( ncFile, "fa_prop1", NC_LONG, 1, &num_fa_blk, &idarrf ) )
1587 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define face block id array" ); }
[ # # ][ # # ]
[ # # ]
1588 : :
1589 : : /* store property name as attribute of property array variable */
1590 [ # # ][ # # ]: 0 : if( NC_NOERR != nc_put_att_text( ncFile, idarrf, "name", strlen( "ID" ), "ID" ) )
1591 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to store face block property name ID" ); }
[ # # ][ # # ]
[ # # ]
1592 : : // determine the number of num_nod_per_face
1593 : : /*
1594 : : num_fa_in_blk1 = 15 ;
1595 : : num_nod_per_fa1 = 58 ;
1596 : :
1597 : : int fbconn1(num_nod_per_fa1) ;
1598 : : fbconn1:elem_type = "nsided" ;
1599 : : int fbepecnt1(num_fa_in_blk1) ;
1600 : : fbepecnt1:entity_type1 = "NODE" ;
1601 : : fbepecnt1:entity_type2 = "FACE" ;
1602 : : */
1603 : :
1604 : 0 : int num_nodes_per_face = 0;
1605 : :
1606 : : int dims[1]; // maybe 1 is enough here
1607 [ # # ][ # # ]: 0 : for( Range::iterator eit = mesh_info.polyhedronFaces.begin(); eit != mesh_info.polyhedronFaces.end(); eit++ )
[ # # ][ # # ]
[ # # ]
1608 : : {
1609 [ # # ]: 0 : EntityHandle polyg = *eit;
1610 : 0 : int nnodes = 0;
1611 : 0 : const EntityHandle* conn = NULL;
1612 [ # # ][ # # ]: 0 : ErrorCode rval = mdbImpl->get_connectivity( polyg, conn, nnodes );MB_CHK_ERR( rval );
[ # # ][ # # ]
1613 : 0 : num_nodes_per_face += nnodes;
1614 : : }
1615 : :
1616 : : // duplicate if needed; default is not duplicate
1617 [ # # ]: 0 : for( int j = 1; j <= num_fa_blocks; j++ )
1618 : : {
1619 : 0 : INS_ID( wname, "num_nod_per_fa%d", j );
1620 [ # # ][ # # ]: 0 : if( nc_def_dim( ncFile, wname, (size_t)num_nodes_per_face, &num_nod_per_fa ) != NC_NOERR )
1621 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of nodes for face block " ); }
[ # # ][ # # ]
[ # # ]
1622 : 0 : dims[0] = num_nod_per_fa;
1623 : 0 : INS_ID( wname, "fbconn%d", j ); // first one, or more
1624 : : int fbconn;
1625 [ # # ][ # # ]: 0 : if( NC_NOERR != nc_def_var( ncFile, wname, NC_LONG, 1, dims, &fbconn ) )
1626 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to create connectivity array for face block " << 1 ); }
[ # # ][ # # ]
[ # # ][ # # ]
1627 [ # # ]: 0 : std::string element_type_string( "nsided" );
1628 [ # # ]: 0 : if( NC_NOERR != nc_put_att_text( ncFile, fbconn, "elem_type", element_type_string.length(),
1629 [ # # ]: 0 : element_type_string.c_str() ) )
1630 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to store element type nsided " ); }
[ # # ][ # # ]
[ # # ]
1631 : :
1632 : 0 : INS_ID( wname, "num_fa_in_blk%d", j ); // first one, or more
1633 : : int num_fa_in_blk;
1634 [ # # ][ # # ]: 0 : if( nc_def_dim( ncFile, wname, (size_t)mesh_info.polyhedronFaces.size(), &num_fa_in_blk ) != NC_NOERR )
[ # # ]
1635 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of nodes for face block " ); }
[ # # ][ # # ]
[ # # ]
1636 : :
1637 : : // fbepecnt
1638 : 0 : INS_ID( wname, "fbepecnt%d", j ); // first one, or more
1639 : : int fbepecnt;
1640 : 0 : dims[0] = num_fa_in_blk;
1641 [ # # ][ # # ]: 0 : if( NC_NOERR != nc_def_var( ncFile, wname, NC_LONG, 1, dims, &fbepecnt ) )
1642 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to create fbepecnt array for block " << 1 ); }
[ # # ][ # # ]
[ # # ][ # # ]
1643 [ # # ][ # # ]: 0 : std::string enttype1( "NODE" );
1644 [ # # ][ # # ]: 0 : if( NC_NOERR != nc_put_att_text( ncFile, fbepecnt, "entity_type1", enttype1.length(), enttype1.c_str() ) )
1645 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to store entity type 1 " ); }
[ # # ][ # # ]
[ # # ]
1646 [ # # ][ # # ]: 0 : std::string enttype2( "FACE" );
1647 [ # # ][ # # ]: 0 : if( NC_NOERR != nc_put_att_text( ncFile, fbepecnt, "entity_type2", enttype2.length(), enttype2.c_str() ) )
1648 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to store entity type 2 " ); }
[ # # ][ # # ]
[ # # ][ # # ]
1649 : 0 : }
1650 : : }
1651 : :
1652 : : // Define element blocks
1653 : :
1654 [ + + ]: 34 : for( unsigned int i = 0; i < block_data.size(); i++ )
1655 : : {
1656 [ + - ]: 25 : MaterialSetData& block = block_data[i];
1657 : :
1658 : 25 : element_block_index = i + 1;
1659 : 25 : int num_el_in_blk = -1, num_att_in_blk = -1;
1660 : : int blk_attrib, connect;
1661 : :
1662 : : /* Define number of elements in this block */
1663 : :
1664 : 25 : INS_ID( wname, "num_el_in_blk%d", element_block_index );
1665 [ + - ][ - + ]: 25 : if( nc_def_dim( ncFile, wname, (size_t)block.number_elements, &num_el_in_blk ) != NC_NOERR )
1666 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of elements/block for block " << i + 1 ); }
[ # # ][ # # ]
[ # # ][ # # ]
1667 : :
1668 : : /* Define number of nodes per element for this block */
1669 : 25 : INS_ID( wname, "num_nod_per_el%d", element_block_index );
1670 : 25 : int num_nod_per_el = -1;
1671 [ + - ]: 25 : if( EXOII_POLYHEDRON != block.element_type )
1672 [ + - ][ - + ]: 25 : if( nc_def_dim( ncFile, wname, (size_t)block.number_nodes_per_element, &num_nod_per_el ) != NC_NOERR )
1673 : : {
1674 [ # # ][ # # ]: 25 : MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of nodes/element for block " << block.id );
[ # # ][ # # ]
[ # # ][ # # ]
1675 : : }
1676 : :
1677 : : /* Define element attribute array for this block */
1678 : : int dims[3];
1679 [ - + ]: 25 : if( block.number_attributes > 0 )
1680 : : {
1681 : 0 : INS_ID( wname, "num_att_in_blk%d", element_block_index );
1682 [ # # ][ # # ]: 0 : if( nc_def_dim( ncFile, wname, (size_t)block.number_attributes, &num_att_in_blk ) != NC_NOERR )
1683 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of attributes in block " << block.id ); }
[ # # ][ # # ]
[ # # ][ # # ]
1684 : :
1685 : 0 : INS_ID( wname, "attrib%d", element_block_index );
1686 : 0 : dims[0] = num_el_in_blk;
1687 : 0 : dims[1] = num_att_in_blk;
1688 [ # # ][ # # ]: 0 : if( NC_NOERR != nc_def_var( ncFile, wname, NC_DOUBLE, 2, dims, &blk_attrib ) )
1689 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define attributes for element block " << block.id ); }
[ # # ][ # # ]
[ # # ][ # # ]
1690 : : }
1691 : :
1692 : : /* Define element connectivity array for this block */
1693 : :
1694 [ + - ][ + - ]: 25 : if( EXOII_POLYGON != block.element_type && EXOII_POLYHEDRON != block.element_type )
1695 : : {
1696 : 25 : INS_ID( wname, "connect%d", element_block_index );
1697 : 25 : dims[0] = num_el_in_blk;
1698 : 25 : dims[1] = num_nod_per_el;
1699 [ + - ][ - + ]: 25 : if( NC_NOERR != nc_def_var( ncFile, wname, NC_LONG, 2, dims, &connect ) )
1700 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to create connectivity array for block " << i + 1 ); }
[ # # ][ # # ]
[ # # ][ # # ]
1701 : :
1702 : : /* Store element type as attribute of connectivity variable */
1703 : :
1704 [ + - ]: 25 : std::string element_type_string( ExoIIUtil::ElementTypeNames[block.element_type] );
1705 [ - + ]: 25 : if( NC_NOERR != nc_put_att_text( ncFile, connect, "elem_type", element_type_string.length(),
1706 [ + - ]: 25 : element_type_string.c_str() ) )
1707 [ # # ][ # # ]: 25 : { MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to store element type name " << (int)block.element_type ); }
[ # # ][ # # ]
[ # # ][ # # ]
[ + - ]
1708 : : }
1709 [ # # ]: 0 : else if( EXOII_POLYGON == block.element_type )
1710 : : {
1711 : 0 : INS_ID( wname, "connect%d", element_block_index );
1712 : : // need to define num_nod_per_el as total number of nodes
1713 : : // ebepecnt1 as number of nodes per polygon
1714 : : /*
1715 : : * int connect1(num_nod_per_el1) ;
1716 : : connect1:elem_type = "nsided" ;
1717 : : int ebepecnt1(num_el_in_blk1) ;
1718 : : ebepecnt1:entity_type1 = "NODE" ;
1719 : : ebepecnt1:entity_type2 = "ELEM" ;*/
1720 : 0 : dims[0] = num_nod_per_el;
1721 [ # # ][ # # ]: 0 : if( NC_NOERR != nc_def_var( ncFile, wname, NC_LONG, 1, dims, &connect ) )
1722 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to create connectivity array for block " << i + 1 ); }
[ # # ][ # # ]
[ # # ][ # # ]
1723 [ # # ]: 0 : std::string element_type_string( "nsided" );
1724 [ # # ]: 0 : if( NC_NOERR != nc_put_att_text( ncFile, connect, "elem_type", element_type_string.length(),
1725 [ # # ]: 0 : element_type_string.c_str() ) )
1726 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to store element type name " << (int)block.element_type ); }
[ # # ][ # # ]
[ # # ][ # # ]
1727 : 0 : INS_ID( wname, "ebepecnt%d", element_block_index );
1728 : : int ebepecnt;
1729 : 0 : dims[0] = num_el_in_blk;
1730 [ # # ][ # # ]: 0 : if( NC_NOERR != nc_def_var( ncFile, wname, NC_LONG, 1, dims, &ebepecnt ) )
1731 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to create ebepecnt array for block " << i + 1 ); }
[ # # ][ # # ]
[ # # ][ # # ]
1732 [ # # ][ # # ]: 0 : std::string etype1( "NODE" );
1733 [ # # ][ # # ]: 0 : if( NC_NOERR != nc_put_att_text( ncFile, ebepecnt, "entity_type1", etype1.length(), etype1.c_str() ) )
1734 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to store entity type1 " << (int)block.element_type ); }
[ # # ][ # # ]
[ # # ][ # # ]
1735 [ # # ][ # # ]: 0 : std::string etype2( "ELEM" );
1736 [ # # ][ # # ]: 0 : if( NC_NOERR != nc_put_att_text( ncFile, ebepecnt, "entity_type2", etype2.length(), etype2.c_str() ) )
1737 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to store entity type2 " << (int)block.element_type ); }
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1738 : : }
1739 [ # # ]: 0 : else if( EXOII_POLYHEDRON == block.element_type )
1740 : : {
1741 : : // INS_ID(wname, "connect%d", element_block_index);
1742 : : /*
1743 : : testn face example: 3 polyh, 15 total faces, 2 shared; 15+2 = 17
1744 : : num_elem = 3 ;
1745 : : num_face = 15 ; // not needed?
1746 : : num_el_blk = 1 ;
1747 : :
1748 : : num_el_in_blk1 = 3 ;
1749 : : num_fac_per_el1 = 17 ;
1750 : :
1751 : : * num faces will be total face conn
1752 : : * num_faces_in_block will be number of faces (non repeated)
1753 : : * num_nodes_per_face will have total face connectivity
1754 : : */
1755 : 0 : int num_faces2 = 0;
1756 : :
1757 [ # # ][ # # ]: 0 : for( Range::iterator eit = block.elements.begin(); eit != block.elements.end(); eit++ )
[ # # ][ # # ]
[ # # ]
1758 : : {
1759 [ # # ]: 0 : EntityHandle polyh = *eit;
1760 : 0 : int nfaces = 0;
1761 : 0 : const EntityHandle* conn = NULL;
1762 [ # # ][ # # ]: 0 : ErrorCode rval = mdbImpl->get_connectivity( polyh, conn, nfaces );MB_CHK_ERR( rval );
[ # # ][ # # ]
1763 : 0 : num_faces2 += nfaces;
1764 : : }
1765 : :
1766 : : int num_fac_per_el;
1767 : :
1768 : 0 : INS_ID( wname, "num_fac_per_el%d", element_block_index );
1769 [ # # ][ # # ]: 0 : if( nc_def_dim( ncFile, wname, (size_t)num_faces2, &num_fac_per_el ) != NC_NOERR )
1770 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of faces per block " << block.id ); }
[ # # ][ # # ]
[ # # ][ # # ]
1771 : :
1772 : : /*
1773 : : int facconn1(num_fac_per_el1) ;
1774 : : facconn1:elem_type = "NFACED" ;
1775 : : int ebepecnt1(num_el_in_blk1) ;
1776 : : ebepecnt1:entity_type1 = "FACE" ;
1777 : : ebepecnt1:entity_type2 = "ELEM" ;
1778 : : */
1779 : :
1780 : : // facconn
1781 : 0 : INS_ID( wname, "facconn%d", element_block_index );
1782 : : int facconn;
1783 : 0 : dims[0] = num_fac_per_el;
1784 [ # # ][ # # ]: 0 : if( NC_NOERR != nc_def_var( ncFile, wname, NC_LONG, 1, dims, &facconn ) )
1785 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to create facconn array for block " << i + 1 ); }
[ # # ][ # # ]
[ # # ][ # # ]
1786 [ # # ]: 0 : std::string etype( "NFACED" );
1787 [ # # ][ # # ]: 0 : if( NC_NOERR != nc_put_att_text( ncFile, facconn, "elem_type", etype.length(), etype.c_str() ) )
1788 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to store elem type " << (int)block.element_type ); }
[ # # ][ # # ]
[ # # ][ # # ]
1789 : :
1790 : : // ebepecnt
1791 : 0 : INS_ID( wname, "ebepecnt%d", element_block_index );
1792 : : int ebepecnt;
1793 : 0 : dims[0] = num_el_in_blk;
1794 [ # # ][ # # ]: 0 : if( NC_NOERR != nc_def_var( ncFile, wname, NC_LONG, 1, dims, &ebepecnt ) )
1795 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to create ebepecnt array for block " << i + 1 ); }
[ # # ][ # # ]
[ # # ][ # # ]
1796 [ # # ][ # # ]: 0 : std::string etype1( "FACE" );
1797 [ # # ][ # # ]: 0 : if( NC_NOERR != nc_put_att_text( ncFile, ebepecnt, "entity_type1", etype1.length(), etype1.c_str() ) )
1798 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to store entity type1 " << (int)block.element_type ); }
[ # # ][ # # ]
[ # # ][ # # ]
1799 [ # # ][ # # ]: 0 : std::string etype2( "ELEM" );
1800 [ # # ][ # # ]: 0 : if( NC_NOERR != nc_put_att_text( ncFile, ebepecnt, "entity_type2", etype2.length(), etype2.c_str() ) )
1801 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to store entity type2 " << (int)block.element_type ); }
[ # # ][ # # ]
[ # # ][ # # ]
1802 : :
1803 [ # # ]: 0 : block.number_nodes_per_element = num_faces2; // connectivity for all polyhedra in block
1804 : : }
1805 : : }
1806 : :
1807 : : /* Node set id array: */
1808 : :
1809 : 9 : int non_empty_nss = 0;
1810 : : // Need to go through nodesets to compute # nodesets, some might be empty
1811 : 9 : std::vector< DirichletSetData >::iterator ns_it;
1812 [ + - ][ + - ]: 18 : for( ns_it = nodeset_data.begin(); ns_it != nodeset_data.end(); ++ns_it )
[ + + ]
1813 : : {
1814 [ + - ][ + - ]: 9 : if( 0 != ( *ns_it ).number_nodes ) non_empty_nss++;
1815 : : }
1816 : :
1817 : 9 : int num_ns = -1;
1818 : 9 : int ns_idstat = -1, ns_idarr = -1;
1819 [ + + ]: 9 : if( non_empty_nss > 0 )
1820 : : {
1821 [ + - ][ - + ]: 3 : if( nc_def_dim( ncFile, "num_node_sets", ( size_t )( non_empty_nss ), &num_ns ) != NC_NOERR )
1822 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of node sets" ); }
[ # # ][ # # ]
[ # # ]
1823 : :
1824 : : /* Node set id status array: */
1825 : :
1826 [ + - ][ - + ]: 3 : if( NC_NOERR != nc_def_var( ncFile, "ns_status", NC_LONG, 1, &num_ns, &ns_idstat ) )
1827 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to create node sets status array" ); }
[ # # ][ # # ]
[ # # ]
1828 : :
1829 : : /* Node set id array: */
1830 [ + - ][ - + ]: 3 : if( NC_NOERR != nc_def_var( ncFile, "ns_prop1", NC_LONG, 1, &num_ns, &ns_idarr ) )
1831 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to create node sets property array" ); }
[ # # ][ # # ]
[ # # ]
1832 : :
1833 : : /* Store property name as attribute of property array variable */
1834 [ + - ][ - + ]: 3 : if( NC_NOERR != nc_put_att_text( ncFile, NC_GLOBAL, "name", strlen( "ID" ), "ID" ) )
1835 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to store node set property name ID" ); }
[ # # ][ # # ]
[ # # ]
1836 : :
1837 : : // Now, define the arrays needed for each node set
1838 : :
1839 : 3 : int index = 0;
1840 : :
1841 [ + + ]: 12 : for( unsigned i = 0; i < nodeset_data.size(); i++ )
1842 : : {
1843 [ + - ][ + - ]: 9 : DirichletSetData node_set = nodeset_data[i];
1844 : :
1845 [ - + ]: 9 : if( 0 == node_set.number_nodes )
1846 : : {
1847 [ # # ][ # # ]: 0 : MB_SET_ERR_CONT( "WriteNCDF: empty nodeset " << node_set.id );
[ # # ][ # # ]
[ # # ][ # # ]
1848 : 0 : continue;
1849 : : }
1850 : 9 : index++;
1851 : :
1852 : 9 : int num_nod_ns = -1;
1853 : 9 : INS_ID( wname, "num_nod_ns%d", index );
1854 [ + - ][ - + ]: 9 : if( nc_def_dim( ncFile, wname, (size_t)node_set.number_nodes, &num_nod_ns ) != NC_NOERR )
1855 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of nodes for set " << node_set.id ); }
[ # # ][ # # ]
[ # # ][ # # ]
1856 : :
1857 : : /* Create variable array in which to store the node set node list */
1858 : 9 : int node_ns = -1;
1859 : 9 : INS_ID( wname, "node_ns%d", index );
1860 [ + - ][ - + ]: 9 : if( NC_NOERR != nc_def_var( ncFile, wname, NC_LONG, 1, &num_nod_ns, &node_ns ) )
1861 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to create node set " << node_set.id << " node list" ); }
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1862 : :
1863 : : // Create distribution factor array
1864 : 9 : int fact_ns = -1;
1865 : 9 : INS_ID( wname, "dist_fact_ns%d", index );
1866 [ + - ][ - + ]: 9 : if( NC_NOERR != nc_def_var( ncFile, wname, NC_DOUBLE, 1, &num_nod_ns, &fact_ns ) )
1867 : : {
1868 [ # # ][ # # ]: 9 : MB_SET_ERR( MB_FAILURE,
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
[ + - - ]
1869 : : "WriteNCDF: failed to create node set " << node_set.id << " distribution factor list" );
1870 : : }
1871 : 9 : }
1872 : : }
1873 : :
1874 : : /* Side set id array: */
1875 : :
1876 : 9 : long non_empty_ss = 0;
1877 : : // Need to go through nodesets to compute # nodesets, some might be empty
1878 : 9 : std::vector< NeumannSetData >::iterator ss_it;
1879 [ + - ][ + - ]: 21 : for( ss_it = sideset_data.begin(); ss_it != sideset_data.end(); ++ss_it )
[ + + ]
1880 : : {
1881 [ + - ][ + - ]: 12 : if( 0 != ( *ss_it ).number_elements ) non_empty_ss++;
1882 : : }
1883 : :
1884 [ + + ]: 9 : if( non_empty_ss > 0 )
1885 : : {
1886 : 3 : int num_ss = -1;
1887 [ + - ][ - + ]: 3 : if( nc_def_dim( ncFile, "num_side_sets", non_empty_ss, &num_ss ) != NC_NOERR )
1888 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of side sets" ); }
[ # # ][ # # ]
[ # # ]
1889 : :
1890 : : /* Side set id status array: */
1891 : 3 : int ss_idstat = -1, ss_idarr = -1;
1892 [ + - ][ - + ]: 3 : if( NC_NOERR != nc_def_var( ncFile, "ss_status", NC_LONG, 1, &num_ss, &ss_idstat ) )
1893 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define side set status" ); }
[ # # ][ # # ]
[ # # ]
1894 : :
1895 : : /* Side set id array: */
1896 [ + - ][ - + ]: 3 : if( NC_NOERR != nc_def_var( ncFile, "ss_prop1", NC_LONG, 1, &num_ss, &ss_idarr ) )
1897 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define side set property" ); }
[ # # ][ # # ]
[ # # ]
1898 : :
1899 : : /* Store property name as attribute of property array variable */
1900 [ + - ][ - + ]: 3 : if( NC_NOERR != nc_put_att_text( ncFile, ss_idarr, "name", strlen( "ID" ), "ID" ) )
1901 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to store side set property name ID" ); }
[ # # ][ # # ]
[ # # ]
1902 : :
1903 : : // Now, define the arrays needed for each side set
1904 : :
1905 : 3 : int index = 0;
1906 [ + + ]: 15 : for( unsigned int i = 0; i < sideset_data.size(); i++ )
1907 : : {
1908 [ + - ][ + - ]: 12 : NeumannSetData side_set = sideset_data[i];
1909 : :
1910 : : // Don't define an empty set
1911 [ - + ]: 12 : if( 0 == side_set.number_elements ) continue;
1912 : :
1913 : 12 : index++;
1914 : :
1915 : 12 : int num_side_ss = -1;
1916 : 12 : int elem_ss = -1, side_ss = -1;
1917 : 12 : INS_ID( wname, "num_side_ss%d", index );
1918 [ + - ][ - + ]: 12 : if( nc_def_dim( ncFile, wname, (size_t)side_set.number_elements, &num_side_ss ) != NC_NOERR )
1919 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of sides in side set " << side_set.id ); }
[ # # ][ # # ]
[ # # ][ # # ]
1920 : :
1921 : 12 : INS_ID( wname, "elem_ss%d", index );
1922 [ + - ][ - + ]: 12 : if( NC_NOERR != nc_def_var( ncFile, wname, NC_LONG, 1, &num_side_ss, &elem_ss ) )
1923 : : {
1924 [ # # ][ # # ]: 0 : MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to create element list for side set "
[ # # ][ # # ]
[ # # ][ # # ]
1925 : : << side_set.id ); /* Exit define mode and return */
1926 : : }
1927 : 12 : INS_ID( wname, "side_ss%d", index );
1928 [ + - ][ - + ]: 12 : if( NC_NOERR != nc_def_var( ncFile, wname, NC_LONG, 1, &num_side_ss, &side_ss ) )
1929 : : {
1930 [ # # ][ # # ]: 0 : MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to create side list for side set "
[ # # ][ # # ]
[ # # ][ # # ]
1931 : : << side_set.id ); /* Exit define mode and return */
1932 : : }
1933 : :
1934 : : // sideset distribution factors
1935 : 12 : int num_df_ss = -1;
1936 : 12 : INS_ID( wname, "num_df_ss%d", index );
1937 [ + - ][ - + ]: 12 : if( nc_def_dim( ncFile, wname, (size_t)side_set.ss_dist_factors.size(), &num_df_ss ) != NC_NOERR )
1938 : : {
1939 [ # # ][ # # ]: 0 : MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of dist factors in side set "
[ # # ][ # # ]
[ # # ][ # # ]
1940 : : << side_set.id ); /* Exit define mode and return */
1941 : : }
1942 : :
1943 : : /* Create variable array in which to store the side set distribution factors */
1944 : :
1945 : 12 : int fact_ss = -1;
1946 : 12 : INS_ID( wname, "dist_fact_ss%d", index );
1947 [ + - ][ - + ]: 12 : if( NC_NOERR != nc_def_var( ncFile, wname, NC_LONG, 1, &num_df_ss, &fact_ss ) )
1948 : : {
1949 [ # # ][ # # ]: 12 : MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to create dist factors list for side set "
[ # # ][ # # ]
[ # # ][ # # ]
[ + - - ]
1950 : : << side_set.id ); /* Exit define mode and return */
1951 : : }
1952 : 12 : }
1953 : : }
1954 : :
1955 : : /* Node coordinate arrays: */
1956 : :
1957 : : int coord, name_coord, dims[3];
1958 : 9 : dims[0] = num_dim;
1959 : 9 : dims[1] = num_nodes;
1960 [ + - ][ - + ]: 9 : if( NC_NOERR != nc_def_var( ncFile, "coord", NC_DOUBLE, 2, dims, &coord ) )
1961 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define node coordinate array" ); }
[ # # ][ # # ]
[ # # ]
1962 : :
1963 : : /* Coordinate names array */
1964 : :
1965 : 9 : dims[0] = num_dim;
1966 : 9 : dims[1] = dim_str;
1967 [ + - ][ - + ]: 9 : if( NC_NOERR != nc_def_var( ncFile, "coor_names", NC_CHAR, 2, dims, &name_coord ) )
1968 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define coordinate name array" ); }
[ # # ][ # # ]
[ # # ]
1969 : :
1970 : : // Define genesis maps if required
1971 : :
1972 [ + - ]: 9 : if( write_maps )
1973 : : {
1974 : : // Element map
1975 : 9 : int elem_map = -1, elem_map2 = -1, node_map = -1;
1976 [ + - ][ - + ]: 9 : if( NC_NOERR != nc_def_var( ncFile, "elem_map", NC_LONG, 1, &num_elem, &elem_map ) )
1977 : : {
1978 [ # # ][ # # ]: 0 : MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to create element map array" ); /* Exit define mode and return */
[ # # ][ # # ]
[ # # ]
1979 : : }
1980 : :
1981 : : // Create the element number map
1982 [ + - ][ - + ]: 9 : if( NC_NOERR != nc_def_var( ncFile, "elem_num_map", NC_LONG, 1, &num_elem, &elem_map2 ) )
1983 : : {
1984 [ # # ][ # # ]: 0 : MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to create element numbering map" ); /* Exit define mode
[ # # ][ # # ]
[ # # ]
1985 : : and return */
1986 : : }
1987 : :
1988 : : // Create node number map
1989 [ + - ][ - + ]: 9 : if( NC_NOERR != nc_def_var( ncFile, "node_num_map", NC_LONG, 1, &num_nodes, &node_map ) )
1990 : : {
1991 [ # # ][ # # ]: 9 : MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to create node numbering map array" ); /* Exit define mode and
[ # # ][ # # ]
[ # # ]
1992 : : return */
1993 : : }
1994 : : }
1995 : :
1996 : : // Define qa records to be used
1997 : :
1998 : 9 : int num_qa_rec = mesh_info.qaRecords.size() / 4;
1999 : 9 : int num_qa = -1;
2000 : :
2001 [ + - ][ - + ]: 9 : if( nc_def_dim( ncFile, "num_qa_rec", (long)num_qa_rec, &num_qa ) != NC_NOERR )
2002 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define qa record array size" ); }
[ # # ][ # # ]
[ # # ]
2003 : :
2004 : : // Define qa array
2005 : : int qa_title;
2006 : 9 : dims[0] = num_qa;
2007 : 9 : dims[1] = dim_four;
2008 : 9 : dims[2] = dim_str;
2009 [ + - ][ - + ]: 9 : if( NC_NOERR != nc_def_var( ncFile, "qa_records", NC_CHAR, 3, dims, &qa_title ) )
2010 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define qa record array" ); }
[ # # ][ # # ]
[ # # ]
2011 : :
2012 : : // Take it out of define mode
2013 [ + - ][ - + ]: 9 : if( NC_NOERR != nc_enddef( ncFile ) ) { MB_SET_ERR( MB_FAILURE, "WriteNCDF: Trouble leaving define mode" ); }
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
2014 : :
2015 : 9 : return MB_SUCCESS;
2016 : : }
2017 : :
2018 : 0 : ErrorCode WriteNCDF::open_file( const char* filename )
2019 : : {
2020 : : // Not a valid filename
2021 [ # # ][ # # ]: 0 : if( strlen( (const char*)filename ) == 0 ) { MB_SET_ERR( MB_FAILURE, "Output Exodus filename not specified" ); }
[ # # ][ # # ]
[ # # ][ # # ]
2022 : :
2023 : 0 : int fail = nc_create( filename, NC_CLOBBER, &ncFile );
2024 : :
2025 : : // File couldn't be opened
2026 [ # # ][ # # ]: 0 : if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "Cannot open " << filename ); }
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
2027 : :
2028 : 0 : return MB_SUCCESS;
2029 : : }
2030 : :
2031 : 12 : ErrorCode WriteNCDF::get_sideset_elems( EntityHandle sideset, int current_sense, Range& forward_elems,
2032 : : Range& reverse_elems )
2033 : : {
2034 [ + - ][ + - ]: 24 : Range ss_elems, ss_meshsets;
2035 : :
2036 : : // Get the sense tag; don't need to check return, might be an error if the tag
2037 : : // hasn't been created yet
2038 : 12 : Tag sense_tag = 0;
2039 [ + - ]: 12 : mdbImpl->tag_get_handle( "SENSE", 1, MB_TYPE_INTEGER, sense_tag );
2040 : :
2041 : : // Get the entities in this set
2042 [ + - ]: 12 : ErrorCode result = mdbImpl->get_entities_by_handle( sideset, ss_elems, true );
2043 [ - + ]: 12 : if( MB_FAILURE == result ) return result;
2044 : :
2045 : : // Now remove the meshsets into the ss_meshsets; first find the first meshset,
2046 [ + - ]: 12 : Range::iterator range_iter = ss_elems.begin();
2047 [ + - ][ + - ]: 96 : while( TYPE_FROM_HANDLE( *range_iter ) != MBENTITYSET && range_iter != ss_elems.end() )
[ + - ][ + - ]
[ + - ][ + + ]
[ + - ]
[ + + # # ]
2048 [ + - ]: 84 : ++range_iter;
2049 : :
2050 : : // Then, if there are some, copy them into ss_meshsets and erase from ss_elems
2051 [ + - ][ + - ]: 12 : if( range_iter != ss_elems.end() )
[ - + ]
2052 : : {
2053 [ # # ][ # # ]: 0 : std::copy( range_iter, ss_elems.end(), range_inserter( ss_meshsets ) );
[ # # ]
2054 [ # # ][ # # ]: 0 : ss_elems.erase( range_iter, ss_elems.end() );
2055 : : }
2056 : :
2057 : : // OK, for the elements, check the sense of this set and copy into the right range
2058 : : // (if the sense is 0, copy into both ranges)
2059 : :
2060 : : // Need to step forward on list until we reach the right dimension
2061 [ + - ]: 12 : Range::iterator dum_it = ss_elems.end();
2062 [ + - ]: 12 : --dum_it;
2063 [ + - ][ + - ]: 12 : int target_dim = CN::Dimension( TYPE_FROM_HANDLE( *dum_it ) );
[ + - ]
2064 [ + - ]: 12 : dum_it = ss_elems.begin();
2065 [ + - ][ + - ]: 12 : while( target_dim != CN::Dimension( TYPE_FROM_HANDLE( *dum_it ) ) && dum_it != ss_elems.end() )
[ + - ][ - + ]
[ # # ][ # # ]
[ # # ][ - + ]
[ - + # # ]
2066 [ # # ]: 0 : ++dum_it;
2067 : :
2068 [ + - ][ + - ]: 12 : if( current_sense == 1 || current_sense == 0 ) std::copy( dum_it, ss_elems.end(), range_inserter( forward_elems ) );
[ + - ][ + - ]
[ + - ]
2069 [ + - ][ + - ]: 12 : if( current_sense == -1 || current_sense == 0 )
2070 [ + - ][ + - ]: 12 : std::copy( dum_it, ss_elems.end(), range_inserter( reverse_elems ) );
[ + - ]
2071 : :
2072 : : // Now loop over the contained meshsets, getting the sense of those and calling this
2073 : : // function recursively
2074 [ + - ][ # # ]: 12 : for( range_iter = ss_meshsets.begin(); range_iter != ss_meshsets.end(); ++range_iter )
[ + - ][ + - ]
[ - + ]
2075 : : {
2076 : : // First get the sense; if it's not there, by convention it's forward
2077 : : int this_sense;
2078 [ # # ][ # # ]: 0 : if( 0 == sense_tag || MB_FAILURE == mdbImpl->tag_get_data( sense_tag, &( *range_iter ), 1, &this_sense ) )
[ # # ][ # # ]
[ # # ]
2079 : 0 : this_sense = 1;
2080 : :
2081 : : // Now get all the entities on this meshset, with the proper (possibly reversed) sense
2082 [ # # ][ # # ]: 0 : get_sideset_elems( *range_iter, this_sense * current_sense, forward_elems, reverse_elems );
2083 : : }
2084 : :
2085 : 24 : return result;
2086 : : }
2087 : :
2088 [ + - ][ + - ]: 228 : } // namespace moab
|
