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 : : #pragma warning( disable : 4786 )
18 : : #endif
19 : :
20 : : #include "ReadNCDF.hpp"
21 : : #include "netcdf.h"
22 : :
23 : : #include <algorithm>
24 : : #include <string>
25 : : #include <assert.h>
26 : : #include <stdio.h>
27 : : #include <string.h>
28 : : #include <cmath>
29 : : #include <sstream>
30 : : #include <iostream>
31 : : #include <map>
32 : :
33 : : #include "moab/CN.hpp"
34 : : #include "moab/Range.hpp"
35 : : #include "moab/Interface.hpp"
36 : : #include "ExoIIUtil.hpp"
37 : : #include "MBTagConventions.hpp"
38 : : #include "Internals.hpp"
39 : : #include "moab/ReadUtilIface.hpp"
40 : : #include "exodus_order.h"
41 : : #include "moab/FileOptions.hpp"
42 : : #include "moab/AdaptiveKDTree.hpp"
43 : : #include "moab/CartVect.hpp"
44 : :
45 : : namespace moab
46 : : {
47 : :
48 : : #define INS_ID( stringvar, prefix, id ) sprintf( stringvar, prefix, id )
49 : :
50 : : #define GET_DIM( ncdim, name, val ) \
51 : : { \
52 : : int gdfail = nc_inq_dimid( ncFile, name, &ncdim ); \
53 : : if( NC_NOERR == gdfail ) \
54 : : { \
55 : : size_t tmp_val; \
56 : : gdfail = nc_inq_dimlen( ncFile, ncdim, &tmp_val ); \
57 : : if( NC_NOERR != gdfail ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Couldn't get dimension length" ); } \
58 : : else \
59 : : val = tmp_val; \
60 : : } \
61 : : else \
62 : : val = 0; \
63 : : }
64 : :
65 : : #define GET_DIMB( ncdim, name, varname, id, val ) \
66 : : INS_ID( name, varname, id ); \
67 : : GET_DIM( ncdim, name, val );
68 : :
69 : : #define GET_VAR( name, id, dims ) \
70 : : { \
71 : : id = -1; \
72 : : int gvfail = nc_inq_varid( ncFile, name, &id ); \
73 : : if( NC_NOERR == gvfail ) \
74 : : { \
75 : : int ndims; \
76 : : gvfail = nc_inq_varndims( ncFile, id, &ndims ); \
77 : : if( NC_NOERR == gvfail ) \
78 : : { \
79 : : dims.resize( ndims ); \
80 : : gvfail = nc_inq_vardimid( ncFile, id, &dims[0] ); \
81 : : if( NC_NOERR != gvfail ) \
82 : : { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Couldn't get variable dimension IDs" ); } \
83 : : } \
84 : : } \
85 : : }
86 : :
87 : : #define GET_1D_INT_VAR( name, id, vals ) \
88 : : { \
89 : : GET_VAR( name, id, vals ); \
90 : : if( -1 != id ) \
91 : : { \
92 : : size_t ntmp; \
93 : : int ivfail = nc_inq_dimlen( ncFile, vals[0], &ntmp ); \
94 : : if( NC_NOERR != ivfail ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Couldn't get dimension length" ); } \
95 : : vals.resize( ntmp ); \
96 : : size_t ntmp1 = 0; \
97 : : ivfail = nc_get_vara_int( ncFile, id, &ntmp1, &ntmp, &vals[0] ); \
98 : : if( NC_NOERR != ivfail ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting variable " << name ); } \
99 : : } \
100 : : }
101 : :
102 : : #define GET_1D_DBL_VAR( name, id, vals ) \
103 : : { \
104 : : std::vector< int > dum_dims; \
105 : : GET_VAR( name, id, dum_dims ); \
106 : : if( -1 != id ) \
107 : : { \
108 : : size_t ntmp; \
109 : : int dvfail = nc_inq_dimlen( ncFile, dum_dims[0], &ntmp ); \
110 : : if( NC_NOERR != dvfail ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Couldn't get dimension length" ); } \
111 : : vals.resize( ntmp ); \
112 : : size_t ntmp1 = 0; \
113 : : dvfail = nc_get_vara_double( ncFile, id, &ntmp1, &ntmp, &vals[0] ); \
114 : : if( NC_NOERR != dvfail ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting variable " << name ); } \
115 : : } \
116 : : }
117 : :
118 : 26 : ReaderIface* ReadNCDF::factory( Interface* iface )
119 : : {
120 [ + - ]: 26 : return new ReadNCDF( iface );
121 : : }
122 : :
123 [ + - ][ + - ]: 52 : ReadNCDF::ReadNCDF( Interface* impl ) : mdbImpl( impl ), max_line_length( -1 ), max_str_length( -1 )
[ + - ][ + - ]
[ + - ][ + - ]
124 : : {
125 [ - + ]: 26 : assert( impl != NULL );
126 [ + - ]: 26 : reset();
127 : :
128 [ + - ]: 26 : impl->query_interface( readMeshIface );
129 : :
130 : : // Initialize in case tag_get_handle fails below
131 : 26 : mMaterialSetTag = 0;
132 : 26 : mDirichletSetTag = 0;
133 : 26 : mNeumannSetTag = 0;
134 : 26 : mHasMidNodesTag = 0;
135 : 26 : mDistFactorTag = 0;
136 : 26 : mQaRecordTag = 0;
137 : 26 : mGlobalIdTag = 0;
138 : :
139 : : //! Get and cache predefined tag handles
140 : : ErrorCode result;
141 : 26 : const int negone = -1;
142 : : result = impl->tag_get_handle( MATERIAL_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mMaterialSetTag,
143 [ + - ]: 26 : MB_TAG_SPARSE | MB_TAG_CREAT, &negone );
144 [ - + ]: 26 : assert( MB_SUCCESS == result );
145 : : result = impl->tag_get_handle( DIRICHLET_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mDirichletSetTag,
146 [ + - ]: 26 : MB_TAG_SPARSE | MB_TAG_CREAT, &negone );
147 [ - + ]: 26 : assert( MB_SUCCESS == result );
148 : : result = impl->tag_get_handle( NEUMANN_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mNeumannSetTag,
149 [ + - ]: 26 : MB_TAG_SPARSE | MB_TAG_CREAT, &negone );
150 [ - + ]: 26 : assert( MB_SUCCESS == result );
151 : 26 : const int mids[] = { -1, -1, -1, -1 };
152 : : result = impl->tag_get_handle( HAS_MID_NODES_TAG_NAME, 4, MB_TYPE_INTEGER, mHasMidNodesTag,
153 [ + - ]: 26 : MB_TAG_SPARSE | MB_TAG_CREAT, mids );
154 [ - + ]: 26 : assert( MB_SUCCESS == result );
155 : : result = impl->tag_get_handle( "distFactor", 0, MB_TYPE_DOUBLE, mDistFactorTag,
156 [ + - ]: 26 : MB_TAG_SPARSE | MB_TAG_VARLEN | MB_TAG_CREAT );
157 [ - + ]: 26 : assert( MB_SUCCESS == result );
158 : : result = impl->tag_get_handle( "qaRecord", 0, MB_TYPE_OPAQUE, mQaRecordTag,
159 [ + - ]: 26 : MB_TAG_SPARSE | MB_TAG_VARLEN | MB_TAG_CREAT );
160 [ - + ]: 26 : assert( MB_SUCCESS == result );
161 [ + - ]: 26 : mGlobalIdTag = impl->globalId_tag();
162 : :
163 [ - + ]: 26 : assert( MB_SUCCESS == result );
164 : : #ifdef NDEBUG
165 : : if( MB_SUCCESS == result ) {}; // Line to avoid compiler warning about unused variable
166 : : #endif
167 : 26 : ncFile = 0;
168 : 26 : }
169 : :
170 : 52 : void ReadNCDF::reset()
171 : : {
172 : 52 : mCurrentMeshHandle = 0;
173 : 52 : vertexOffset = 0;
174 : :
175 : 52 : numberNodes_loading = 0;
176 : 52 : numberElements_loading = 0;
177 : 52 : numberElementBlocks_loading = 0;
178 : 52 : numberFaceBlocks_loading = 0;
179 : 52 : numberNodeSets_loading = 0;
180 : 52 : numberSideSets_loading = 0;
181 : 52 : numberDimensions_loading = 0;
182 : :
183 [ - + ]: 52 : if( !blocksLoading.empty() ) blocksLoading.clear();
184 : :
185 [ - + ]: 52 : if( !nodesInLoadedBlocks.empty() ) nodesInLoadedBlocks.clear();
186 : :
187 [ - + ]: 52 : if( !faceBlocksLoading.empty() ) faceBlocksLoading.clear();
188 : 52 : }
189 : :
190 : 78 : ReadNCDF::~ReadNCDF()
191 : : {
192 : 26 : mdbImpl->release_interface( readMeshIface );
193 [ - + ]: 52 : }
194 : :
195 : 0 : ErrorCode ReadNCDF::read_tag_values( const char* file_name, const char* tag_name, const FileOptions& /* opts */,
196 : : std::vector< int >& id_array, const SubsetList* subset_list )
197 : : {
198 [ # # ]: 0 : if( subset_list )
199 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_UNSUPPORTED_OPERATION, "ExodusII reader supports subset read only by material ID" ); }
[ # # ][ # # ]
[ # # ]
200 : :
201 : : // Open netcdf/exodus file
202 : 0 : int fail = nc_open( file_name, 0, &ncFile );
203 [ # # ]: 0 : if( NC_NOWRITE != fail )
204 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FILE_DOES_NOT_EXIST, "ReadNCDF:: problem opening Netcdf/Exodus II file " << file_name ); }
[ # # ][ # # ]
[ # # ][ # # ]
205 : :
206 : : // 1. Read the header
207 : 0 : ErrorCode rval = read_exodus_header();
208 [ # # ]: 0 : if( MB_FAILURE == rval ) return rval;
209 : :
210 : 0 : int count = 0;
211 : : const char* prop;
212 : 0 : const char* blocks = "eb_prop1";
213 : 0 : const char* nodesets = "ns_prop1";
214 : 0 : const char* sidesets = "ss_prop1";
215 : :
216 [ # # ]: 0 : if( !strcmp( tag_name, MATERIAL_SET_TAG_NAME ) )
217 : : {
218 : 0 : count = numberElementBlocks_loading;
219 : 0 : prop = blocks;
220 : : }
221 [ # # ]: 0 : else if( !strcmp( tag_name, DIRICHLET_SET_TAG_NAME ) )
222 : : {
223 : 0 : count = numberNodeSets_loading;
224 : 0 : prop = nodesets;
225 : : }
226 [ # # ]: 0 : else if( !strcmp( tag_name, NEUMANN_SET_TAG_NAME ) )
227 : : {
228 : 0 : count = numberSideSets_loading;
229 : 0 : prop = sidesets;
230 : : }
231 : : else
232 : : {
233 : 0 : ncFile = 0;
234 : 0 : return MB_TAG_NOT_FOUND;
235 : : }
236 : :
237 [ # # ]: 0 : if( count )
238 : : {
239 : 0 : int nc_var = -1;
240 [ # # ][ # # ]: 0 : GET_1D_INT_VAR( prop, nc_var, id_array );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
241 [ # # ][ # # ]: 0 : if( !nc_var ) { MB_SET_ERR( MB_FAILURE, "Problem getting prop variable" ); }
[ # # ][ # # ]
[ # # ][ # # ]
242 : : }
243 : :
244 : : // Close the file
245 : 0 : fail = nc_close( ncFile );
246 [ # # ][ # # ]: 0 : if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "Trouble closing file" ); }
[ # # ][ # # ]
[ # # ][ # # ]
247 : :
248 : 0 : ncFile = 0;
249 : 0 : return MB_SUCCESS;
250 : : }
251 : :
252 : 26 : ErrorCode ReadNCDF::load_file( const char* exodus_file_name, const EntityHandle* file_set, const FileOptions& opts,
253 : : const ReaderIface::SubsetList* subset_list, const Tag* file_id_tag )
254 : : {
255 : : ErrorCode status;
256 : : int fail;
257 : :
258 : 26 : int num_blocks = 0;
259 : 26 : const int* blocks_to_load = 0;
260 [ - + ]: 26 : if( subset_list )
261 : : {
262 [ # # ][ # # ]: 0 : if( subset_list->tag_list_length > 1 || !strcmp( subset_list->tag_list[0].tag_name, MATERIAL_SET_TAG_NAME ) )
263 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_UNSUPPORTED_OPERATION, "ExodusII reader supports subset read only by material ID" ); }
[ # # ][ # # ]
[ # # ]
264 [ # # ]: 0 : if( subset_list->num_parts )
265 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_UNSUPPORTED_OPERATION, "ExodusII reader does not support mesh partitioning" ); }
[ # # ][ # # ]
[ # # ]
266 : 0 : blocks_to_load = subset_list->tag_list[0].tag_values;
267 : 0 : num_blocks = subset_list->tag_list[0].num_tag_values;
268 : : }
269 : :
270 : : // This function directs the reading of an exoii file, but doesn't do any of
271 : : // the actual work
272 : :
273 : : // See if opts has tdata.
274 : : ErrorCode rval;
275 [ + - ]: 26 : std::string s;
276 [ + - ]: 26 : rval = opts.get_str_option( "tdata", s );
277 [ - + ][ # # ]: 26 : if( MB_SUCCESS == rval && !s.empty() )
[ - + ]
278 [ # # ]: 0 : return update( exodus_file_name, opts, num_blocks, blocks_to_load, *file_set );
279 : :
280 [ + - ]: 26 : reset();
281 : :
282 : : // 0. Open the file.
283 : :
284 : : // open netcdf/exodus file
285 [ + - ]: 26 : fail = nc_open( exodus_file_name, 0, &ncFile );
286 [ + + ]: 26 : if( NC_NOERR != fail )
287 [ + - ][ + - ]: 1 : { MB_SET_ERR( MB_FILE_DOES_NOT_EXIST, "ReadNCDF:: problem opening Netcdf/Exodus II file " << exodus_file_name ); }
[ + - ][ + - ]
[ - + ][ + - ]
288 : :
289 : : // 1. Read the header
290 [ + - ]: 25 : status = read_exodus_header();
291 [ - + ]: 25 : if( MB_FAILURE == status ) return status;
292 : :
293 [ + - ]: 25 : status = mdbImpl->get_entities_by_handle( 0, initRange );
294 [ - + ]: 25 : if( MB_FAILURE == status ) return status;
295 : :
296 : : // 2. Read the nodes unless they've already been read before
297 [ + - ]: 25 : status = read_nodes( file_id_tag );
298 [ - + ]: 25 : if( MB_FAILURE == status ) return status;
299 : :
300 : : // 3.
301 : : // extra for polyhedra blocks
302 [ - + ]: 25 : if( numberFaceBlocks_loading > 0 )
303 : : {
304 [ # # ]: 0 : status = read_face_blocks_headers();
305 [ # # ]: 0 : if( MB_FAILURE == status ) return status;
306 : : }
307 : :
308 [ + - ]: 25 : status = read_block_headers( blocks_to_load, num_blocks );
309 [ - + ]: 25 : if( MB_FAILURE == status ) return status;
310 : :
311 : : // 4. Read elements (might not read them, depending on active blocks)
312 [ - + ]: 25 : if( numberFaceBlocks_loading > 0 )
313 : : {
314 [ # # ]: 0 : status = read_polyhedra_faces();
315 [ # # ]: 0 : if( MB_FAILURE == status ) return status;
316 : : }
317 [ + - ]: 25 : status = read_elements( file_id_tag );
318 [ - + ]: 25 : if( MB_FAILURE == status ) return status;
319 : :
320 : : // 5. Read global ids
321 [ + - ]: 25 : status = read_global_ids();
322 [ - + ]: 25 : if( MB_FAILURE == status ) return status;
323 : :
324 : : // 6. Read nodesets
325 [ + - ]: 25 : status = read_nodesets();
326 [ - + ]: 25 : if( MB_FAILURE == status ) return status;
327 : :
328 : : // 7. Read sidesets
329 [ + - ]: 25 : status = read_sidesets();
330 [ - + ]: 25 : if( MB_FAILURE == status ) return status;
331 : :
332 : : // 8. Read qa records
333 [ + + ]: 25 : if( file_set )
334 : : {
335 [ + - ]: 2 : status = read_qa_records( *file_set );
336 [ - + ]: 2 : if( MB_FAILURE == status ) return status;
337 : : }
338 : :
339 : : // What about properties???
340 : :
341 : : // Close the file
342 [ + - ]: 25 : fail = nc_close( ncFile );
343 [ - + ][ # # ]: 25 : if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "Trouble closing file" ); }
[ # # ][ # # ]
[ # # ][ # # ]
344 : :
345 : 25 : ncFile = 0;
346 : 26 : return MB_SUCCESS;
347 : : }
348 : :
349 : 25 : ErrorCode ReadNCDF::read_exodus_header()
350 : : {
351 : 25 : CPU_WORD_SIZE = sizeof( double ); // With ExodusII version 2, all floats
352 : 25 : IO_WORD_SIZE = sizeof( double ); // should be changed to doubles
353 : :
354 : : // NetCDF doesn't check its own limits on file read, so check
355 : : // them here so it doesn't corrupt memory any more than absolutely
356 : : // necessary.
357 : : int ndims;
358 [ + - ]: 25 : int fail = nc_inq_ndims( ncFile, &ndims );
359 [ + - ][ - + ]: 25 : if( NC_NOERR != fail || ndims > NC_MAX_DIMS )
360 : : {
361 [ # # ][ # # ]: 0 : MB_SET_ERR( MB_FAILURE, "ReadNCDF: File contains " << ndims << " dims but NetCDF library supports only "
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
362 : : << (int)NC_MAX_DIMS );
363 : : }
364 : : int nvars;
365 [ + - ]: 25 : fail = nc_inq_nvars( ncFile, &nvars );
366 [ - + ]: 25 : if( nvars > NC_MAX_VARS )
367 : : {
368 [ # # ][ # # ]: 0 : MB_SET_ERR( MB_FAILURE, "ReadNCDF: File contains " << nvars << " vars but NetCDF library supports only "
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
369 : : << (int)NC_MAX_VARS );
370 : : }
371 : :
372 : : // Get the attributes
373 : :
374 : : // Get the word size, scalar value
375 : : nc_type att_type;
376 : : size_t att_len;
377 [ + - ]: 25 : fail = nc_inq_att( ncFile, NC_GLOBAL, "floating_point_word_size", &att_type, &att_len );
378 [ - + ]: 25 : if( NC_NOERR != fail )
379 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting floating_point_word_size attribute" ); }
[ # # ][ # # ]
[ # # ]
380 [ + - ][ - + ]: 25 : if( att_type != NC_INT || att_len != 1 )
381 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Word size didn't have type int or size 1" ); }
[ # # ][ # # ]
[ # # ]
382 [ + - ]: 25 : fail = nc_get_att_int( ncFile, NC_GLOBAL, "floating_point_word_size", &IO_WORD_SIZE );
383 [ - + ][ # # ]: 25 : if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Trouble getting word size" ); }
[ # # ][ # # ]
[ # # ][ # # ]
384 : :
385 : : // Exodus version
386 [ + - ]: 25 : fail = nc_inq_att( ncFile, NC_GLOBAL, "version", &att_type, &att_len );
387 [ - + ][ # # ]: 25 : if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting version attribute" ); }
[ # # ][ # # ]
[ # # ][ # # ]
388 [ + - ][ - + ]: 25 : if( att_type != NC_FLOAT || att_len != 1 )
389 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Version didn't have type float or size 1" ); }
[ # # ][ # # ]
[ # # ]
390 : : // float version = temp_att->as_float(0);
391 : :
392 : : // Read in initial variables
393 : : int temp_dim;
394 [ + - ][ + - ]: 25 : GET_DIM( temp_dim, "num_dim", numberDimensions_loading );
[ + - ][ - + ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
395 [ + - ][ + - ]: 25 : GET_DIM( temp_dim, "num_nodes", numberNodes_loading );
[ + - ][ - + ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
396 [ + - ][ + - ]: 25 : GET_DIM( temp_dim, "num_elem", numberElements_loading );
[ + - ][ - + ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
397 [ + - ][ + - ]: 25 : GET_DIM( temp_dim, "num_el_blk", numberElementBlocks_loading );
[ + - ][ - + ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
398 [ + - ][ - + ]: 25 : GET_DIM( temp_dim, "num_fa_blk", numberFaceBlocks_loading ); // for polyhedra blocks
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
399 [ + - ][ + - ]: 25 : GET_DIM( temp_dim, "num_elem", numberElements_loading );
[ + - ][ - + ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
400 [ + - ][ + + ]: 25 : GET_DIM( temp_dim, "num_node_sets", numberNodeSets_loading );
[ + - ][ - + ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
401 [ + - ][ + + ]: 25 : GET_DIM( temp_dim, "num_side_sets", numberSideSets_loading );
[ + - ][ - + ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
402 [ + - ][ + - ]: 25 : GET_DIM( temp_dim, "len_string", max_str_length );
[ + - ][ - + ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
403 [ + - ][ + - ]: 25 : GET_DIM( temp_dim, "len_line", max_line_length );
[ + - ][ - + ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
404 : :
405 : : // Title; why are we even bothering if we're not going to keep it???
406 [ + - ]: 25 : fail = nc_inq_att( ncFile, NC_GLOBAL, "title", &att_type, &att_len );
407 [ - + ][ # # ]: 25 : if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting title attribute" ); }
[ # # ][ # # ]
[ # # ][ # # ]
408 [ - + ][ # # ]: 25 : if( att_type != NC_CHAR ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: title didn't have type char" ); }
[ # # ][ # # ]
[ # # ][ # # ]
409 [ + - ]: 25 : char* title = new char[att_len + 1];
410 [ + - ]: 25 : fail = nc_get_att_text( ncFile, NC_GLOBAL, "title", title );
411 [ - + ]: 25 : if( NC_NOERR != fail )
412 : : {
413 [ # # ]: 0 : delete[] title;
414 [ # # ][ # # ]: 0 : MB_SET_ERR( MB_FAILURE, "ReadNCDF:: trouble getting title" );
[ # # ][ # # ]
[ # # ]
415 : : }
416 [ + - ]: 25 : delete[] title;
417 : :
418 : 25 : return MB_SUCCESS;
419 : : }
420 : :
421 : 25 : ErrorCode ReadNCDF::read_nodes( const Tag* file_id_tag )
422 : : {
423 : : // Read the nodes into memory
424 : :
425 [ - + ]: 25 : assert( 0 != ncFile );
426 : :
427 : : // Create a sequence to hold the node coordinates
428 : : // Get the current number of entities and start at the next slot
429 : :
430 : 25 : EntityHandle node_handle = 0;
431 [ + - ]: 25 : std::vector< double* > arrays;
432 [ + - ]: 25 : readMeshIface->get_node_coords( 3, numberNodes_loading, MB_START_ID, node_handle, arrays );
433 : :
434 [ + - ]: 25 : vertexOffset = ID_FROM_HANDLE( node_handle ) - MB_START_ID;
435 : :
436 : : // Read in the coordinates
437 : : int fail;
438 : 25 : int coord = 0;
439 [ + - ]: 25 : nc_inq_varid( ncFile, "coord", &coord );
440 : :
441 : : // Single var for all coords
442 : 25 : size_t start[2] = { 0, 0 }, count[2] = { 1, static_cast< size_t >( numberNodes_loading ) };
443 [ + - ]: 25 : if( coord )
444 : : {
445 [ + + ]: 100 : for( int d = 0; d < numberDimensions_loading; ++d )
446 : : {
447 : 75 : start[0] = d;
448 [ + - ][ + - ]: 75 : fail = nc_get_vara_double( ncFile, coord, start, count, arrays[d] );
449 [ - + ]: 75 : if( NC_NOERR != fail )
450 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting " << (char)( 'x' + d ) << " coord array" ); }
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
451 : : }
452 : : }
453 : : // Var for each coord
454 : : else
455 : : {
456 : 0 : char varname[] = "coord ";
457 [ # # ]: 0 : for( int d = 0; d < numberDimensions_loading; ++d )
458 : : {
459 : 0 : varname[5] = 'x' + (char)d;
460 [ # # ]: 0 : fail = nc_inq_varid( ncFile, varname, &coord );
461 [ # # ]: 0 : if( NC_NOERR != fail )
462 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting " << (char)( 'x' + d ) << " coord variable" ); }
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
463 : :
464 [ # # ][ # # ]: 0 : fail = nc_get_vara_double( ncFile, coord, start, &count[1], arrays[d] );
465 [ # # ]: 0 : if( NC_NOERR != fail )
466 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting " << (char)( 'x' + d ) << " coord array" ); }
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
467 : : }
468 : : }
469 : :
470 : : // Zero out any coord values that are in database but not in file
471 : : // (e.g. if MOAB has 3D coords but file is 2D then set Z coords to zero.)
472 [ - + ]: 25 : for( unsigned d = numberDimensions_loading; d < arrays.size(); ++d )
473 [ # # ][ # # ]: 0 : std::fill( arrays[d], arrays[d] + numberNodes_loading, 0.0 );
[ # # ]
474 : :
475 [ - + ]: 25 : if( file_id_tag )
476 : : {
477 [ # # ]: 0 : Range nodes;
478 [ # # ]: 0 : nodes.insert( node_handle, node_handle + numberNodes_loading - 1 );
479 [ # # ]: 0 : readMeshIface->assign_ids( *file_id_tag, nodes, vertexOffset );
480 : : }
481 : :
482 : 25 : return MB_SUCCESS;
483 : : }
484 : :
485 : 25 : ErrorCode ReadNCDF::read_block_headers( const int* blocks_to_load, const int num_blocks )
486 : : {
487 : : // Get the element block ids; keep this in a separate list,
488 : : // which is not offset by blockIdOffset; this list used later for
489 : : // reading block connectivity
490 : :
491 : : // Get the ids of all the blocks of this file we're reading in
492 [ + - ]: 25 : std::vector< int > block_ids( numberElementBlocks_loading );
493 : 25 : int nc_block_ids = -1;
494 [ + - ][ + - ]: 25 : GET_1D_INT_VAR( "eb_prop1", nc_block_ids, block_ids );
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ - + ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ + - ]
[ + - ][ + - ]
[ - + ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ + - ][ + - ]
[ + - ][ - + ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
495 [ - + ][ # # ]: 25 : if( -1 == nc_block_ids ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting eb_prop1 variable" ); }
[ # # ][ # # ]
[ # # ][ # # ]
496 : :
497 : 25 : int exodus_id = 1;
498 : :
499 : : // If the active_block_id_list is NULL all blocks are active.
500 [ - + ][ # # ]: 25 : if( NULL == blocks_to_load || 0 == num_blocks ) { blocks_to_load = &block_ids[0]; }
[ + - ]
501 : :
502 [ + - ]: 50 : std::vector< int > new_blocks( blocks_to_load, blocks_to_load + numberElementBlocks_loading );
503 : :
504 : 25 : std::vector< int >::iterator iter, end_iter;
505 : 25 : iter = block_ids.begin();
506 : 25 : end_iter = block_ids.end();
507 : :
508 : : // Read header information and initialize header-type block information
509 : : int temp_dim;
510 [ + - ]: 50 : std::vector< char > temp_string_storage( max_str_length + 1 );
511 [ + - ]: 25 : char* temp_string = &temp_string_storage[0];
512 : 25 : int block_seq_id = 1;
513 : :
514 [ + - ][ + - ]: 135 : for( ; iter != end_iter; ++iter, block_seq_id++ )
[ + + ]
515 : : {
516 : : int num_elements;
517 : :
518 [ + - ][ + - ]: 110 : GET_DIMB( temp_dim, temp_string, "num_el_in_blk%d", block_seq_id, num_elements );
[ + - ][ - + ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
519 : :
520 : : // Don't read element type string for now, since it's an attrib
521 : : // on the connectivity
522 : : // Get the entity type corresponding to this ExoII element type
523 : : // ExoIIElementType elem_type =
524 : : // ExoIIUtil::static_element_name_to_type(element_type_string);
525 : :
526 : : // Tag each element block(mesh set) with enum for ElementType (SHELL, QUAD4, ....etc)
527 [ + - ]: 110 : ReadBlockData block_data;
528 : 110 : block_data.elemType = EXOII_MAX_ELEM_TYPE;
529 [ + - ]: 110 : block_data.blockId = *iter;
530 : 110 : block_data.startExoId = exodus_id;
531 : 110 : block_data.numElements = num_elements;
532 : :
533 : : // If block is in 'blocks_to_load'----load it!
534 [ + - ][ + - ]: 110 : if( std::find( new_blocks.begin(), new_blocks.end(), *iter ) != new_blocks.end() )
[ + - ][ + - ]
535 : 110 : block_data.reading_in = true;
536 : : else
537 : 0 : block_data.reading_in = false;
538 : :
539 [ + - ]: 110 : blocksLoading.push_back( block_data );
540 : 110 : exodus_id += num_elements;
541 : 110 : }
542 : :
543 : 50 : return MB_SUCCESS;
544 : : }
545 : :
546 : 0 : ErrorCode ReadNCDF::read_face_blocks_headers()
547 : : {
548 : : // Get the ids of all the blocks of this file we're reading in
549 [ # # ]: 0 : std::vector< int > fblock_ids( numberFaceBlocks_loading );
550 : 0 : int nc_fblock_ids = -1;
551 [ # # ][ # # ]: 0 : GET_1D_INT_VAR( "fa_prop1", nc_fblock_ids, fblock_ids );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
552 [ # # ][ # # ]: 0 : if( -1 == nc_fblock_ids ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting fa_prop1 variable" ); }
[ # # ][ # # ]
[ # # ][ # # ]
553 : :
554 : : int temp_dim;
555 [ # # ]: 0 : std::vector< char > temp_string_storage( max_str_length + 1 );
556 [ # # ]: 0 : char* temp_string = &temp_string_storage[0];
557 : : // int block_seq_id = 1;
558 : 0 : int exodus_id = 1;
559 : :
560 [ # # ]: 0 : for( int i = 1; i <= numberFaceBlocks_loading; i++ )
561 : : {
562 : : int num_elements;
563 : :
564 [ # # ][ # # ]: 0 : GET_DIMB( temp_dim, temp_string, "num_fa_in_blk%d", i, num_elements );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
565 : :
566 : : // Don't read element type string for now, since it's an attrib
567 : : // on the connectivity
568 : : // Get the entity type corresponding to this ExoII element type
569 : : // ExoIIElementType elem_type =
570 : : // ExoIIUtil::static_element_name_to_type(element_type_string);
571 : :
572 : : // Tag each element block(mesh set) with enum for ElementType (SHELL, QUAD4, ....etc)
573 : : ReadFaceBlockData fblock_data;
574 : 0 : fblock_data.faceBlockId = i;
575 : 0 : fblock_data.startExoId = exodus_id;
576 : 0 : fblock_data.numElements = num_elements;
577 : :
578 [ # # ]: 0 : faceBlocksLoading.push_back( fblock_data );
579 : 0 : exodus_id += num_elements;
580 : : }
581 : 0 : return MB_SUCCESS;
582 : : }
583 : 0 : ErrorCode ReadNCDF::read_polyhedra_faces()
584 : : {
585 : : int temp_dim;
586 [ # # ]: 0 : std::vector< char > temp_string_storage( max_str_length + 1 );
587 [ # # ]: 0 : char* temp_string = &temp_string_storage[0];
588 [ # # ]: 0 : nodesInLoadedBlocks.resize( numberNodes_loading + 1 );
589 : 0 : size_t start[1] = { 0 }, count[1]; // one dim arrays only here!
590 : 0 : std::vector< ReadFaceBlockData >::iterator this_it;
591 : 0 : this_it = faceBlocksLoading.begin();
592 : :
593 : 0 : int fblock_seq_id = 1;
594 [ # # ]: 0 : std::vector< int > dims;
595 : : int nc_var, fail;
596 : :
597 [ # # ][ # # ]: 0 : for( ; this_it != faceBlocksLoading.end(); ++this_it, fblock_seq_id++ )
[ # # ]
598 : : {
599 : : int num_fa_in_blk;
600 [ # # ][ # # ]: 0 : GET_DIMB( temp_dim, temp_string, "num_fa_in_blk%d", fblock_seq_id, num_fa_in_blk );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
601 : : int num_nod_per_fa;
602 [ # # ][ # # ]: 0 : GET_DIMB( temp_dim, temp_string, "num_nod_per_fa%d", fblock_seq_id, num_nod_per_fa );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
603 : : // Get the ncdf connect variable and the element type
604 : 0 : INS_ID( temp_string, "fbconn%d", fblock_seq_id );
605 [ # # ][ # # ]: 0 : GET_VAR( temp_string, nc_var, dims );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
606 [ # # ]: 0 : std::vector< int > fbconn;
607 [ # # ]: 0 : fbconn.resize( num_nod_per_fa );
608 : 0 : count[0] = num_nod_per_fa;
609 : :
610 [ # # ][ # # ]: 0 : fail = nc_get_vara_int( ncFile, nc_var, start, count, &fbconn[0] );
611 [ # # ][ # # ]: 0 : if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting face polyhedra connectivity " ); }
[ # # ][ # # ]
[ # # ][ # # ]
612 [ # # ][ # # ]: 0 : std::vector< int > fbepecnt;
613 [ # # ]: 0 : fbepecnt.resize( num_fa_in_blk );
614 : 0 : INS_ID( temp_string, "fbepecnt%d", fblock_seq_id );
615 [ # # ][ # # ]: 0 : GET_VAR( temp_string, nc_var, dims );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
616 : 0 : count[0] = num_fa_in_blk;
617 [ # # ][ # # ]: 0 : fail = nc_get_vara_int( ncFile, nc_var, start, count, &fbepecnt[0] );
618 [ # # ][ # # ]: 0 : if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting face polyhedra connectivity " ); }
[ # # ][ # # ]
[ # # ][ # # ]
619 : : // now we can create some polygons
620 [ # # ][ # # ]: 0 : std::vector< EntityHandle > polyConn;
621 : 0 : int ix = 0;
622 [ # # ][ # # ]: 0 : for( int i = 0; i < num_fa_in_blk; i++ )
623 : : {
624 [ # # ][ # # ]: 0 : polyConn.resize( fbepecnt[i] );
625 [ # # ][ # # ]: 0 : for( int j = 0; j < fbepecnt[i]; j++ )
626 : : {
627 [ # # ][ # # ]: 0 : nodesInLoadedBlocks[fbconn[ix]] = 1;
628 [ # # ][ # # ]: 0 : polyConn[j] = vertexOffset + fbconn[ix++];
629 : : }
630 : : EntityHandle newp;
631 : : /*
632 : : * ErrorCode create_element(const EntityType type,
633 : : const EntityHandle *connectivity,
634 : : const int num_vertices,
635 : : EntityHandle &element_handle)
636 : : */
637 [ # # ][ # # ]: 0 : if( mdbImpl->create_element( MBPOLYGON, &polyConn[0], fbepecnt[i], newp ) != MB_SUCCESS ) return MB_FAILURE;
[ # # ][ # # ]
638 : :
639 : : // add the element in order
640 [ # # ]: 0 : polyfaces.push_back( newp );
641 : : }
642 : 0 : }
643 : 0 : return MB_SUCCESS;
644 : : }
645 : :
646 : 25 : ErrorCode ReadNCDF::read_elements( const Tag* file_id_tag )
647 : : {
648 : : // Read in elements
649 : :
650 : 25 : int result = 0;
651 : :
652 : : // Initialize the nodeInLoadedBlocks vector
653 [ + - ]: 25 : if( nodesInLoadedBlocks.size() < ( size_t )( numberNodes_loading + 1 ) )
654 [ + - ]: 25 : nodesInLoadedBlocks.resize( numberNodes_loading + 1 ); // this could be repeated?
655 [ + - ]: 25 : memset( &nodesInLoadedBlocks[0], 0, ( numberNodes_loading + 1 ) * sizeof( char ) );
656 : :
657 : 25 : std::vector< ReadBlockData >::iterator this_it;
658 : 25 : this_it = blocksLoading.begin();
659 : :
660 : : int temp_dim;
661 [ + - ]: 25 : std::vector< char > temp_string_storage( max_str_length + 1 );
662 [ + - ]: 25 : char* temp_string = &temp_string_storage[0];
663 : :
664 : : int nc_var;
665 : 25 : int block_seq_id = 1;
666 [ + - ]: 50 : std::vector< int > dims;
667 : 25 : size_t start[2] = { 0, 0 }, count[2];
668 : :
669 [ + - ][ + - ]: 135 : for( ; this_it != blocksLoading.end(); ++this_it, block_seq_id++ )
[ + + ]
670 : : {
671 : : // If this block isn't to be read in --- continue
672 [ + - ][ - + ]: 110 : if( !( *this_it ).reading_in ) continue;
673 : :
674 : : // Get some information about this block
675 [ + - ]: 110 : int block_id = ( *this_it ).blockId;
676 : 110 : EntityHandle* conn = 0;
677 : :
678 : : // Get the ncdf connect variable and the element type
679 : 110 : INS_ID( temp_string, "connect%d", block_seq_id );
680 [ + - ][ + - ]: 110 : GET_VAR( temp_string, nc_var, dims );
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ - + ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
681 [ + - ][ - + ]: 110 : if( -1 == nc_var || 0 == nc_var )
682 : : { // try other var, for polyhedra blocks
683 : : // it could be polyhedra block, defined by fbconn and NFACED attribute
684 : 0 : INS_ID( temp_string, "facconn%d", block_seq_id );
685 [ # # ][ # # ]: 0 : GET_VAR( temp_string, nc_var, dims );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
686 : :
687 [ # # ][ # # ]: 0 : if( -1 == nc_var || 0 == nc_var )
688 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting connec or faccon variable" ); }
[ # # ][ # # ]
[ # # ]
689 : : }
690 : : nc_type att_type;
691 : : size_t att_len;
692 [ + - ]: 110 : int fail = nc_inq_att( ncFile, nc_var, "elem_type", &att_type, &att_len );
693 [ - + ][ # # ]: 110 : if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting elem type attribute" ); }
[ # # ][ # # ]
[ # # ][ # # ]
694 [ + - ]: 110 : std::vector< char > dum_str( att_len + 1 );
695 [ + - ][ + - ]: 110 : fail = nc_get_att_text( ncFile, nc_var, "elem_type", &dum_str[0] );
696 [ - + ][ # # ]: 110 : if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting elem type" ); }
[ # # ][ # # ]
[ # # ][ # # ]
697 [ + - ][ + - ]: 110 : ExoIIElementType elem_type = ExoIIUtil::static_element_name_to_type( &dum_str[0] );
698 [ + - ]: 110 : ( *this_it ).elemType = elem_type;
699 : :
700 [ + - ]: 110 : int verts_per_element = ExoIIUtil::VerticesPerElement[( *this_it ).elemType];
701 [ + - ]: 110 : int number_nodes = ( *this_it ).numElements * verts_per_element;
702 [ + - ]: 110 : const EntityType mb_type = ExoIIUtil::ExoIIElementMBEntity[( *this_it ).elemType];
703 : :
704 [ - + ]: 110 : if( mb_type == MBPOLYGON )
705 : : {
706 : : // need to read another variable, num_nod_per_el, and
707 : : int num_nodes_poly_conn;
708 [ # # ][ # # ]: 0 : GET_DIMB( temp_dim, temp_string, "num_nod_per_el%d", block_seq_id, num_nodes_poly_conn );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
709 : : // read the connec1 array, which is of dimensions num_nodes_poly_conn (only one )
710 [ # # ]: 0 : std::vector< int > connec;
711 [ # # ]: 0 : connec.resize( num_nodes_poly_conn );
712 : 0 : count[0] = num_nodes_poly_conn;
713 : :
714 [ # # ][ # # ]: 0 : fail = nc_get_vara_int( ncFile, nc_var, start, count, &connec[0] );
715 [ # # ][ # # ]: 0 : if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting polygon connectivity " ); }
[ # # ][ # # ]
[ # # ][ # # ]
716 : : // ebepecnt1
717 [ # # ][ # # ]: 0 : std::vector< int > ebec;
718 [ # # ][ # # ]: 0 : ebec.resize( this_it->numElements );
719 : : // Get the ncdf connect variable and the element type
720 : 0 : INS_ID( temp_string, "ebepecnt%d", block_seq_id );
721 [ # # ][ # # ]: 0 : GET_VAR( temp_string, nc_var, dims );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
722 [ # # ]: 0 : count[0] = this_it->numElements;
723 [ # # ][ # # ]: 0 : fail = nc_get_vara_int( ncFile, nc_var, start, count, &ebec[0] );
724 [ # # ][ # # ]: 0 : if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting polygon nodes per elem " ); }
[ # # ][ # # ]
[ # # ][ # # ]
725 : : EntityHandle ms_handle;
726 [ # # ][ # # ]: 0 : if( mdbImpl->create_meshset( MESHSET_SET | MESHSET_TRACK_OWNER, ms_handle ) != MB_SUCCESS )
727 : 0 : return MB_FAILURE;
728 : : // create polygons one by one, and put them in the list
729 : : // also need to read the number of nodes for each polygon, then create
730 [ # # ][ # # ]: 0 : std::vector< EntityHandle > polyConn;
731 : 0 : int ix = 0;
732 [ # # ][ # # ]: 0 : for( int i = 0; i < this_it->numElements; i++ )
733 : : {
734 [ # # ][ # # ]: 0 : polyConn.resize( ebec[i] );
735 [ # # ][ # # ]: 0 : for( int j = 0; j < ebec[i]; j++ )
736 : : {
737 [ # # ][ # # ]: 0 : nodesInLoadedBlocks[connec[ix]] = 1;
738 [ # # ][ # # ]: 0 : polyConn[j] = vertexOffset + connec[ix++];
739 : : }
740 : : EntityHandle newp;
741 : : /*
742 : : * ErrorCode create_element(const EntityType type,
743 : : const EntityHandle *connectivity,
744 : : const int num_vertices,
745 : : EntityHandle &element_handle)
746 : : */
747 [ # # ][ # # ]: 0 : if( mdbImpl->create_element( MBPOLYGON, &polyConn[0], ebec[i], newp ) != MB_SUCCESS ) return MB_FAILURE;
[ # # ][ # # ]
748 : :
749 : : // add the element in order
750 [ # # ][ # # ]: 0 : this_it->polys.push_back( newp );
751 [ # # ][ # # ]: 0 : if( mdbImpl->add_entities( ms_handle, &newp, 1 ) != MB_SUCCESS ) return MB_FAILURE;
752 : : }
753 : : // Set the block id with an offset
754 [ # # ][ # # ]: 0 : if( mdbImpl->tag_set_data( mMaterialSetTag, &ms_handle, 1, &block_id ) != MB_SUCCESS ) return MB_FAILURE;
755 [ # # ][ # # ]: 0 : if( mdbImpl->tag_set_data( mGlobalIdTag, &ms_handle, 1, &block_id ) != MB_SUCCESS ) return MB_FAILURE;
[ # # ]
756 : : }
757 [ - + ]: 110 : else if( mb_type == MBPOLYHEDRON )
758 : : {
759 : : // need to read another variable, num_fac_per_el
760 : : int num_fac_per_el;
761 [ # # ][ # # ]: 0 : GET_DIMB( temp_dim, temp_string, "num_fac_per_el%d", block_seq_id, num_fac_per_el );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
762 : : // read the fbconn array, which is of dimensions num_nod_per_fa (only one dimension)
763 [ # # ]: 0 : std::vector< int > facconn;
764 [ # # ]: 0 : facconn.resize( num_fac_per_el );
765 : 0 : count[0] = num_fac_per_el;
766 : :
767 [ # # ][ # # ]: 0 : fail = nc_get_vara_int( ncFile, nc_var, start, count, &facconn[0] );
768 [ # # ][ # # ]: 0 : if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting polygon connectivity " ); }
[ # # ][ # # ]
[ # # ][ # # ]
769 : : // ebepecnt1
770 [ # # ][ # # ]: 0 : std::vector< int > ebepecnt;
771 [ # # ][ # # ]: 0 : ebepecnt.resize( this_it->numElements );
772 : : // Get the ncdf connect variable and the element type
773 : 0 : INS_ID( temp_string, "ebepecnt%d", block_seq_id );
774 [ # # ][ # # ]: 0 : GET_VAR( temp_string, nc_var, dims );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
775 [ # # ]: 0 : count[0] = this_it->numElements;
776 [ # # ][ # # ]: 0 : fail = nc_get_vara_int( ncFile, nc_var, start, count, &ebepecnt[0] );
777 [ # # ][ # # ]: 0 : if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting polygon nodes per elem " ); }
[ # # ][ # # ]
[ # # ][ # # ]
778 : : EntityHandle ms_handle;
779 [ # # ][ # # ]: 0 : if( mdbImpl->create_meshset( MESHSET_SET | MESHSET_TRACK_OWNER, ms_handle ) != MB_SUCCESS )
780 : 0 : return MB_FAILURE;
781 : : // create polygons one by one, and put them in the list
782 : : // also need to read the number of nodes for each polygon, then create
783 [ # # ][ # # ]: 0 : std::vector< EntityHandle > polyConn;
784 : 0 : int ix = 0;
785 [ # # ][ # # ]: 0 : for( int i = 0; i < this_it->numElements; i++ )
786 : : {
787 [ # # ][ # # ]: 0 : polyConn.resize( ebepecnt[i] );
788 [ # # ][ # # ]: 0 : for( int j = 0; j < ebepecnt[i]; j++ )
789 : : {
790 [ # # ][ # # ]: 0 : polyConn[j] = polyfaces[facconn[ix++] - 1];
[ # # ]
791 : : }
792 : : EntityHandle newp;
793 : : /*
794 : : * ErrorCode create_element(const EntityType type,
795 : : const EntityHandle *connectivity,
796 : : const int num_vertices,
797 : : EntityHandle &element_handle)
798 : : */
799 [ # # ][ # # ]: 0 : if( mdbImpl->create_element( MBPOLYHEDRON, &polyConn[0], ebepecnt[i], newp ) != MB_SUCCESS )
[ # # ][ # # ]
800 : 0 : return MB_FAILURE;
801 : :
802 : : // add the element in order
803 [ # # ][ # # ]: 0 : this_it->polys.push_back( newp );
804 [ # # ][ # # ]: 0 : if( mdbImpl->add_entities( ms_handle, &newp, 1 ) != MB_SUCCESS ) return MB_FAILURE;
805 : : }
806 : : // Set the block id with an offset
807 [ # # ][ # # ]: 0 : if( mdbImpl->tag_set_data( mMaterialSetTag, &ms_handle, 1, &block_id ) != MB_SUCCESS ) return MB_FAILURE;
808 [ # # ][ # # ]: 0 : if( mdbImpl->tag_set_data( mGlobalIdTag, &ms_handle, 1, &block_id ) != MB_SUCCESS ) return MB_FAILURE;
[ # # ]
809 : : }
810 : : else // this is regular
811 : : {
812 : : // Allocate an array to read in connectivity data
813 [ + - ][ + - ]: 110 : readMeshIface->get_element_connect( this_it->numElements, verts_per_element, mb_type, this_it->startExoId,
814 [ + - ][ + - ]: 220 : this_it->startMBId, conn );
815 : :
816 : : // Create a range for this sequence of elements
817 : : EntityHandle start_range, end_range;
818 [ + - ]: 110 : start_range = ( *this_it ).startMBId;
819 [ + - ]: 110 : end_range = start_range + ( *this_it ).numElements - 1;
820 : :
821 [ + - ]: 110 : Range new_range( start_range, end_range );
822 : : // Range<EntityHandle> new_range((*this_it).startMBId,
823 : : // (*this_it).startMBId + (*this_it).numElements - 1);
824 : :
825 : : // Create a MBSet for this block and set the material tag
826 : :
827 : : EntityHandle ms_handle;
828 [ + - ][ - + ]: 110 : if( mdbImpl->create_meshset( MESHSET_SET | MESHSET_TRACK_OWNER, ms_handle ) != MB_SUCCESS )
829 : 0 : return MB_FAILURE;
830 : :
831 [ + - ][ - + ]: 110 : if( mdbImpl->add_entities( ms_handle, new_range ) != MB_SUCCESS ) return MB_FAILURE;
832 : :
833 : : int mid_nodes[4];
834 [ + - ]: 110 : CN::HasMidNodes( mb_type, verts_per_element, mid_nodes );
835 [ + - ][ - + ]: 110 : if( mdbImpl->tag_set_data( mHasMidNodesTag, &ms_handle, 1, mid_nodes ) != MB_SUCCESS ) return MB_FAILURE;
836 : :
837 : : // Just a check because the following code won't work if this case fails
838 : : assert( sizeof( EntityHandle ) >= sizeof( int ) );
839 : :
840 : : // tmp_ptr is of type int* and points at the same place as conn
841 : 110 : int* tmp_ptr = reinterpret_cast< int* >( conn );
842 : :
843 : : // Read the connectivity into that memory, this will take only part of the array
844 : : // 1/2 if sizeof(EntityHandle) == 64 bits.
845 [ + - ]: 110 : count[0] = this_it->numElements;
846 : 110 : count[1] = verts_per_element;
847 [ + - ]: 110 : fail = nc_get_vara_int( ncFile, nc_var, start, count, tmp_ptr );
848 [ - + ][ # # ]: 110 : if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting connectivity" ); }
[ # # ][ # # ]
[ # # ][ # # ]
849 : :
850 : : // Convert from exodus indices to vertex handles.
851 : : // Iterate backwards in case handles are larger than ints.
852 [ + + ]: 2063806 : for( int i = number_nodes - 1; i >= 0; --i )
853 : : {
854 [ - + ]: 2063696 : if( (unsigned)tmp_ptr[i] >= nodesInLoadedBlocks.size() )
855 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "Invalid node ID in block connectivity" ); }
[ # # ][ # # ]
[ # # ]
856 [ + - ]: 2063696 : nodesInLoadedBlocks[tmp_ptr[i]] = 1;
857 : 2063696 : conn[i] = static_cast< EntityHandle >( tmp_ptr[i] ) + vertexOffset;
858 : : }
859 : :
860 : : // Adjust connectivity order if necessary
861 : 110 : const int* reorder = exodus_elem_order_map[mb_type][verts_per_element];
862 [ + + ][ + - ]: 110 : if( reorder ) ReadUtilIface::reorder( reorder, conn, this_it->numElements, verts_per_element );
[ + - ]
863 : :
864 [ + - ][ + - ]: 110 : readMeshIface->update_adjacencies( ( *this_it ).startMBId, ( *this_it ).numElements,
865 [ + - ][ + - ]: 220 : ExoIIUtil::VerticesPerElement[( *this_it ).elemType], conn );
866 : :
867 [ - + ]: 110 : if( result == -1 )
868 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: error getting element connectivity for block " << block_id ); }
[ # # ][ # # ]
[ # # ][ # # ]
869 : :
870 : : // Set the block id with an offset
871 [ + - ][ - + ]: 110 : if( mdbImpl->tag_set_data( mMaterialSetTag, &ms_handle, 1, &block_id ) != MB_SUCCESS ) return MB_FAILURE;
872 [ + - ][ - + ]: 110 : if( mdbImpl->tag_set_data( mGlobalIdTag, &ms_handle, 1, &block_id ) != MB_SUCCESS ) return MB_FAILURE;
873 : :
874 [ - + ]: 110 : if( file_id_tag )
875 : : {
876 [ # # ]: 0 : Range range;
877 [ # # ][ # # ]: 0 : range.insert( this_it->startMBId, this_it->startMBId + this_it->numElements - 1 );
[ # # ][ # # ]
878 [ # # ][ # # ]: 110 : readMeshIface->assign_ids( *file_id_tag, range, this_it->startExoId );
[ + - ]
879 [ + - ]: 110 : }
880 : : } // end regular block (not polygon)
881 : 110 : }
882 : :
883 : 50 : return MB_SUCCESS;
884 : : }
885 : :
886 : 25 : ErrorCode ReadNCDF::read_global_ids()
887 : : {
888 : : // Read in the map from the exodus file
889 [ + - ][ + - ]: 25 : std::vector< int > ids( std::max( numberElements_loading, numberNodes_loading ) );
890 : :
891 : 25 : int varid = -1;
892 [ + - ][ + - ]: 25 : GET_1D_INT_VAR( "elem_map", varid, ids );
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ - + ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ + - ]
[ + - ][ + - ]
[ - + ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ + - ][ + - ]
[ + - ][ - + ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
893 [ + - ]: 25 : if( -1 != varid )
894 : : {
895 : 25 : std::vector< ReadBlockData >::iterator iter;
896 : 25 : int id_pos = 0;
897 [ + - ][ + - ]: 135 : for( iter = blocksLoading.begin(); iter != blocksLoading.end(); ++iter )
[ + + ]
898 : : {
899 [ + - ][ + - ]: 110 : if( iter->reading_in )
900 : : {
901 [ + - ][ + - ]: 110 : if( iter->elemType != EXOII_POLYGON && iter->elemType != EXOII_POLYHEDRON )
[ + - ][ + - ]
[ + - ]
902 : : {
903 [ + - ][ + - ]: 110 : if( iter->startMBId != 0 )
904 : : {
905 [ + - ][ + - ]: 110 : Range range( iter->startMBId, iter->startMBId + iter->numElements - 1 );
[ + - ][ + - ]
906 [ + - ][ + - ]: 110 : ErrorCode error = mdbImpl->tag_set_data( mGlobalIdTag, range, &ids[id_pos] );
907 [ - + ]: 110 : if( error != MB_SUCCESS ) return error;
908 [ + - ][ + - ]: 110 : id_pos += iter->numElements;
909 : : }
910 : : else
911 : 110 : return MB_FAILURE;
912 : : }
913 : : else // polygons or polyhedrons; use id from elements
914 : : {
915 [ # # ][ # # ]: 110 : for( std::vector< EntityHandle >::iterator eit = iter->polys.begin(); eit != iter->polys.end();
[ # # ][ # # ]
[ # # ]
916 : : eit++, id_pos++ )
917 : : {
918 [ # # ]: 0 : EntityHandle peh = *eit;
919 [ # # ][ # # ]: 0 : if( mdbImpl->tag_set_data( mGlobalIdTag, &peh, 1, &ids[id_pos] ) != MB_SUCCESS )
[ # # ]
920 : 0 : return MB_FAILURE;
921 : : }
922 : : }
923 : : }
924 : : }
925 : : }
926 : :
927 : : // Read in node map next
928 : 25 : varid = -1;
929 [ + - ][ + + ]: 25 : GET_1D_INT_VAR( "node_num_map", varid, ids );
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ - + ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ + + ]
[ + - ][ + - ]
[ - + ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ + - ][ + - ]
[ + - ][ - + ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
930 [ + + ]: 25 : if( -1 != varid )
931 : : {
932 [ + - ]: 22 : Range range( MB_START_ID + vertexOffset, MB_START_ID + vertexOffset + numberNodes_loading - 1 );
933 [ + - ][ + - ]: 22 : ErrorCode error = mdbImpl->tag_set_data( mGlobalIdTag, range, &ids[0] );
934 [ - + ][ # # ]: 22 : if( MB_SUCCESS != error ) MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem setting node global ids" );
[ # # ][ # # ]
[ # # ][ # # ]
[ + - ]
935 : : }
936 : :
937 : 25 : return MB_SUCCESS;
938 : : }
939 : :
940 : 25 : ErrorCode ReadNCDF::read_nodesets()
941 : : {
942 : : // Read in the nodesets for the model
943 : :
944 [ + + ]: 25 : if( 0 == numberNodeSets_loading ) return MB_SUCCESS;
945 [ + - ]: 13 : std::vector< int > id_array( numberNodeSets_loading );
946 : :
947 : : // Read in the nodeset ids
948 : : int nc_var;
949 [ + - ][ + - ]: 13 : GET_1D_INT_VAR( "ns_prop1", nc_var, id_array );
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ - + ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ + - ]
[ + - ][ + - ]
[ - + ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ + - ][ + - ]
[ + - ][ - + ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
950 [ - + ][ # # ]: 13 : if( -1 == nc_var ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting ns_prop1 variable" ); }
[ # # ][ # # ]
[ # # ][ # # ]
951 : :
952 : : // Use a vector of ints to read node handles
953 [ + - ]: 26 : std::vector< int > node_handles;
954 : :
955 : : int i;
956 [ + - ]: 26 : std::vector< char > temp_string_storage( max_str_length + 1 );
957 [ + - ]: 13 : char* temp_string = &temp_string_storage[0];
958 : : int temp_dim;
959 [ + + ]: 50 : for( i = 0; i < numberNodeSets_loading; i++ )
960 : : {
961 : : // Get nodeset parameters
962 : : int number_nodes_in_set;
963 : : int number_dist_factors_in_set;
964 : :
965 [ + - ][ + - ]: 37 : GET_DIMB( temp_dim, temp_string, "num_nod_ns%d", i + 1, number_nodes_in_set );
[ + - ][ - + ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
966 [ + - ][ - + ]: 37 : GET_DIMB( temp_dim, temp_string, "num_df_ns%d", i + 1, number_dist_factors_in_set );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
967 : :
968 : : // Need to new a vector to store dist. factors
969 : : // This vector * gets stored as a tag on the sideset meshset
970 [ + - ]: 37 : std::vector< double > temp_dist_factor_vector( number_nodes_in_set );
971 [ - + ]: 37 : if( number_dist_factors_in_set != 0 )
972 : : {
973 : 0 : INS_ID( temp_string, "dist_fact_ns%d", i + 1 );
974 [ # # ][ # # ]: 0 : GET_1D_DBL_VAR( temp_string, temp_dim, temp_dist_factor_vector );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
975 [ # # ][ # # ]: 0 : if( -1 == temp_dim ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting dist fact variable" ); }
[ # # ][ # # ]
[ # # ][ # # ]
976 : : }
977 : :
978 : : // Size new arrays and get ids and distribution factors
979 [ + + ]: 37 : if( node_handles.size() < (unsigned int)number_nodes_in_set )
980 : : {
981 [ + - ]: 13 : node_handles.reserve( number_nodes_in_set );
982 [ + - ]: 13 : node_handles.resize( number_nodes_in_set );
983 : : }
984 : :
985 : 37 : INS_ID( temp_string, "node_ns%d", i + 1 );
986 : 37 : int temp_var = -1;
987 [ + - ][ + - ]: 37 : GET_1D_INT_VAR( temp_string, temp_var, node_handles );
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ - + ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ + - ]
[ + - ][ + - ]
[ - + ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ + - ][ + - ]
[ + - ][ - + ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
988 [ - + ][ # # ]: 37 : if( -1 == temp_var ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting nodeset node variable" ); }
[ # # ][ # # ]
[ # # ][ # # ]
989 : :
990 : : // Maybe there is already a nodesets meshset here we can append to
991 [ + - ]: 74 : Range child_meshsets;
[ + - - ]
992 [ + - ][ - + ]: 37 : if( mdbImpl->get_entities_by_handle( 0, child_meshsets ) != MB_SUCCESS ) return MB_FAILURE;
993 : :
994 [ + - ][ + - ]: 37 : child_meshsets = subtract( child_meshsets, initRange );
995 : :
996 [ + - ][ + - ]: 37 : Range::iterator iter, end_iter;
997 [ + - ]: 37 : iter = child_meshsets.begin();
998 [ + - ]: 37 : end_iter = child_meshsets.end();
999 : :
1000 : 37 : EntityHandle ns_handle = 0;
1001 [ + - ][ + - ]: 4597 : for( ; iter != end_iter; ++iter )
[ + + ]
1002 : : {
1003 : : int nodeset_id;
1004 [ + - ][ + - ]: 4560 : if( mdbImpl->tag_get_data( mDirichletSetTag, &( *iter ), 1, &nodeset_id ) != MB_SUCCESS ) continue;
[ - + ]
1005 : :
1006 [ + - ][ - + ]: 4560 : if( id_array[i] == nodeset_id )
1007 : : {
1008 : : // Found the meshset
1009 [ # # ]: 0 : ns_handle = *iter;
1010 : 4560 : break;
1011 : : }
1012 : : }
1013 : :
1014 [ + - ]: 74 : std::vector< EntityHandle > nodes_of_nodeset;
[ + - - ]
1015 [ - + ]: 37 : if( ns_handle )
1016 [ # # ][ # # ]: 0 : if( mdbImpl->get_entities_by_handle( ns_handle, nodes_of_nodeset, true ) != MB_SUCCESS ) return MB_FAILURE;
1017 : :
1018 : : // Make these into entity handles
1019 : : // TODO: could we have read it into EntityHandle sized array in the first place?
1020 : : int j, temp;
1021 [ + - ]: 74 : std::vector< EntityHandle > nodes;
[ + - - ]
1022 [ + - ]: 74 : std::vector< double > dist_factor_vector;
[ + - - ]
1023 [ + + ]: 481 : for( j = 0; j < number_nodes_in_set; j++ )
1024 : : {
1025 : : // See if this node is one we're currently reading in
1026 [ + - ][ + - ]: 444 : if( nodesInLoadedBlocks[node_handles[j]] == 1 )
[ + - ]
1027 : : {
1028 : : // Make sure that it already isn't in a nodeset
1029 [ + - ][ + - ]: 444 : unsigned int node_id = CREATE_HANDLE( MBVERTEX, node_handles[j] + vertexOffset, temp );
1030 [ - + ][ # # ]: 888 : if( !ns_handle ||
[ + - ]
1031 [ # # ][ # # ]: 444 : std::find( nodes_of_nodeset.begin(), nodes_of_nodeset.end(), node_id ) == nodes_of_nodeset.end() )
[ - + ][ + - ]
[ # # # # ]
1032 : : {
1033 [ + - ]: 444 : nodes.push_back( node_id );
1034 : :
1035 [ - + ][ # # ]: 444 : if( number_dist_factors_in_set != 0 ) dist_factor_vector.push_back( temp_dist_factor_vector[j] );
[ # # ]
1036 : : }
1037 : : }
1038 : : }
1039 : :
1040 : : // No nodes to add
1041 [ - + ]: 37 : if( nodes.empty() ) continue;
1042 : :
1043 : : // If there was no meshset found --> create one
1044 [ + - ]: 37 : if( ns_handle == 0 )
1045 : : {
1046 [ + - ][ - + ]: 37 : if( mdbImpl->create_meshset( MESHSET_ORDERED | MESHSET_TRACK_OWNER, ns_handle ) != MB_SUCCESS )
1047 : 0 : return MB_FAILURE;
1048 : :
1049 : : // Set a tag signifying dirichlet bc
1050 : : // TODO: create this tag another way
1051 : :
1052 [ + - ]: 37 : int nodeset_id = id_array[i];
1053 [ + - ][ - + ]: 37 : if( mdbImpl->tag_set_data( mDirichletSetTag, &ns_handle, 1, &nodeset_id ) != MB_SUCCESS ) return MB_FAILURE;
1054 [ + - ][ - + ]: 37 : if( mdbImpl->tag_set_data( mGlobalIdTag, &ns_handle, 1, &nodeset_id ) != MB_SUCCESS ) return MB_FAILURE;
1055 : :
1056 [ - + ]: 37 : if( !dist_factor_vector.empty() )
1057 : : {
1058 : 0 : int size = dist_factor_vector.size();
1059 [ # # ]: 0 : const void* data = &dist_factor_vector[0];
1060 [ # # ][ # # ]: 0 : if( mdbImpl->tag_set_by_ptr( mDistFactorTag, &ns_handle, 1, &data, &size ) != MB_SUCCESS )
1061 : 37 : return MB_FAILURE;
1062 : : }
1063 : : }
1064 [ # # ]: 0 : else if( !dist_factor_vector.empty() )
1065 : : {
1066 : : // Append dist factors to vector
1067 : 0 : const void* ptr = 0;
1068 : 0 : int size = 0;
1069 [ # # ][ # # ]: 0 : if( mdbImpl->tag_get_by_ptr( mDistFactorTag, &ns_handle, 1, &ptr, &size ) != MB_SUCCESS ) return MB_FAILURE;
1070 : :
1071 : 0 : const double* data = reinterpret_cast< const double* >( ptr );
1072 [ # # ]: 0 : dist_factor_vector.reserve( dist_factor_vector.size() + size );
1073 [ # # ][ # # ]: 0 : std::copy( data, data + size, std::back_inserter( dist_factor_vector ) );
1074 : 0 : size = dist_factor_vector.size();
1075 [ # # ]: 0 : ptr = &dist_factor_vector[0];
1076 [ # # ][ # # ]: 0 : if( mdbImpl->tag_set_by_ptr( mDistFactorTag, &ns_handle, 1, &ptr, &size ) != MB_SUCCESS ) return MB_FAILURE;
1077 : : }
1078 : :
1079 : : // Add the nodes to the meshset
1080 [ + - ][ + - ]: 37 : if( mdbImpl->add_entities( ns_handle, &nodes[0], nodes.size() ) != MB_SUCCESS ) return MB_FAILURE;
[ - + ]
[ + - - ]
1081 : 37 : }
1082 : :
1083 : 38 : return MB_SUCCESS;
1084 : : }
1085 : :
1086 : 25 : ErrorCode ReadNCDF::read_sidesets()
1087 : : {
1088 : : // Uhhh if you read the same file (previously_read==true) then blow away the
1089 : : // sidesets pertaining to this file? How do you do that? If you read a
1090 : : // new file make sure all the offsets are correct.
1091 : :
1092 : : // If not loading any sidesets -- exit
1093 [ + + ]: 25 : if( 0 == numberSideSets_loading ) return MB_SUCCESS;
1094 : :
1095 : : // Read in the sideset ids
1096 [ + - ]: 13 : std::vector< int > id_array( numberSideSets_loading );
1097 : 13 : int temp_var = -1;
1098 [ + - ][ + - ]: 13 : GET_1D_INT_VAR( "ss_prop1", temp_var, id_array );
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ - + ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ + - ]
[ + - ][ + - ]
[ - + ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ + - ][ + - ]
[ + - ][ - + ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
1099 [ - + ][ # # ]: 13 : if( -1 == temp_var ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting ss_prop1 variable" ); }
[ # # ][ # # ]
[ # # ][ # # ]
1100 : :
1101 : : // Create a side set for each one
1102 : : int number_sides_in_set;
1103 : : int number_dist_factors_in_set;
1104 : :
1105 : : // Maybe there is already a sidesets meshset here we can append to
1106 [ + - ]: 26 : Range child_meshsets;
1107 [ + - ][ - + ]: 13 : if( mdbImpl->get_entities_by_type( 0, MBENTITYSET, child_meshsets ) != MB_SUCCESS ) return MB_FAILURE;
1108 : :
1109 [ + - ][ + - ]: 13 : child_meshsets = subtract( child_meshsets, initRange );
1110 : :
1111 [ + - ][ + - ]: 13 : Range::iterator iter, end_iter;
1112 : :
1113 : : int i;
1114 [ + - ]: 26 : std::vector< char > temp_string_storage( max_str_length + 1 );
1115 [ + - ]: 13 : char* temp_string = &temp_string_storage[0];
1116 : : int temp_dim;
1117 [ + + ]: 62 : for( i = 0; i < numberSideSets_loading; i++ )
1118 : : {
1119 : : // Get sideset parameters
1120 [ + - ][ + - ]: 49 : GET_DIMB( temp_dim, temp_string, "num_side_ss%d", i + 1, number_sides_in_set );
[ + - ][ - + ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1121 [ + - ][ + - ]: 49 : GET_DIMB( temp_dim, temp_string, "num_df_ss%d", i + 1, number_dist_factors_in_set );
[ + - ][ - + ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1122 : :
1123 : : // Size new arrays and get element and side lists
1124 [ + - ]: 49 : std::vector< int > side_list( number_sides_in_set );
1125 [ + - ][ + - ]: 98 : std::vector< int > element_list( number_sides_in_set );
1126 : 49 : INS_ID( temp_string, "side_ss%d", i + 1 );
1127 : 49 : temp_var = -1;
1128 [ + - ][ + - ]: 49 : GET_1D_INT_VAR( temp_string, temp_var, side_list );
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ - + ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ + - ]
[ + - ][ + - ]
[ - + ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ + - ][ + - ]
[ + - ][ - + ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
1129 [ - + ][ # # ]: 49 : if( -1 == temp_var ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting sideset side variable" ); }
[ # # ][ # # ]
[ # # ][ # # ]
1130 : :
1131 : 49 : INS_ID( temp_string, "elem_ss%d", i + 1 );
1132 : 49 : temp_var = -1;
1133 [ + - ][ + - ]: 49 : GET_1D_INT_VAR( temp_string, temp_var, element_list );
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ - + ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ + - ]
[ + - ][ + - ]
[ - + ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ + - ][ + - ]
[ + - ][ - + ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
1134 [ - + ][ # # ]: 49 : if( -1 == temp_var ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting sideset elem variable" ); }
[ # # ][ # # ]
[ # # ][ # # ]
1135 : :
1136 [ + - ][ + - ]: 98 : std::vector< double > temp_dist_factor_vector;
1137 [ + - ][ + - ]: 98 : std::vector< EntityHandle > entities_to_add, reverse_entities;
[ + - ][ + - ]
1138 : : // Create the sideset entities
1139 [ + - ][ + - ]: 98 : if( create_ss_elements( &element_list[0], &side_list[0], number_sides_in_set, number_dist_factors_in_set,
[ - + ]
1140 [ + - ]: 49 : entities_to_add, reverse_entities, temp_dist_factor_vector, i + 1 ) != MB_SUCCESS )
1141 : 0 : return MB_FAILURE;
1142 : :
1143 : : // If there are elements to add
1144 [ - + ][ # # ]: 49 : if( !entities_to_add.empty() || !reverse_entities.empty() )
[ + - ]
1145 : : {
1146 [ + - ]: 49 : iter = child_meshsets.begin();
1147 [ + - ]: 49 : end_iter = child_meshsets.end();
1148 : :
1149 : 49 : EntityHandle ss_handle = 0;
1150 [ + - ][ + - ]: 435 : for( ; iter != end_iter; ++iter )
[ + + ]
1151 : : {
1152 : : int sideset_id;
1153 [ + - ][ + - ]: 386 : if( mdbImpl->tag_get_data( mNeumannSetTag, &( *iter ), 1, &sideset_id ) != MB_SUCCESS ) continue;
[ - + ]
1154 : :
1155 [ + - ][ - + ]: 386 : if( id_array[i] == sideset_id )
1156 : : {
1157 : : // Found the meshset
1158 [ # # ]: 0 : ss_handle = *iter;
1159 : 386 : break;
1160 : : }
1161 : : }
1162 : :
1163 : : // If we didn't find a sideset already
1164 [ + - ]: 49 : if( ss_handle == 0 )
1165 : : {
1166 [ + - ][ - + ]: 49 : if( mdbImpl->create_meshset( MESHSET_ORDERED | MESHSET_TRACK_OWNER, ss_handle ) != MB_SUCCESS )
1167 : 0 : return MB_FAILURE;
1168 : :
1169 [ - + ]: 49 : if( ss_handle == 0 ) return MB_FAILURE;
1170 : :
1171 [ + - ]: 49 : int sideset_id = id_array[i];
1172 [ + - ][ - + ]: 49 : if( mdbImpl->tag_set_data( mNeumannSetTag, &ss_handle, 1, &sideset_id ) != MB_SUCCESS )
1173 : 0 : return MB_FAILURE;
1174 [ + - ][ - + ]: 49 : if( mdbImpl->tag_set_data( mGlobalIdTag, &ss_handle, 1, &sideset_id ) != MB_SUCCESS ) return MB_FAILURE;
1175 : :
1176 [ - + ]: 49 : if( !reverse_entities.empty() )
1177 : : {
1178 : : // Also make a reverse set to put in this set
1179 : : EntityHandle reverse_set;
1180 [ # # ][ # # ]: 0 : if( mdbImpl->create_meshset( MESHSET_SET | MESHSET_TRACK_OWNER, reverse_set ) != MB_SUCCESS )
1181 : 49 : return MB_FAILURE;
1182 : :
1183 : : // Add the reverse set to the sideset set and the entities to the reverse set
1184 [ # # ]: 0 : ErrorCode result = mdbImpl->add_entities( ss_handle, &reverse_set, 1 );
1185 [ # # ]: 0 : if( MB_SUCCESS != result ) return result;
1186 : :
1187 [ # # ][ # # ]: 0 : result = mdbImpl->add_entities( reverse_set, &reverse_entities[0], reverse_entities.size() );
1188 [ # # ]: 0 : if( MB_SUCCESS != result ) return result;
1189 : :
1190 : : // Set the reverse tag
1191 : : Tag sense_tag;
1192 : 0 : int dum_sense = 0;
1193 : : result = mdbImpl->tag_get_handle( "NEUSET_SENSE", 1, MB_TYPE_INTEGER, sense_tag,
1194 [ # # ]: 0 : MB_TAG_SPARSE | MB_TAG_CREAT, &dum_sense );
1195 [ # # ]: 0 : if( result != MB_SUCCESS ) return result;
1196 : 0 : dum_sense = -1;
1197 [ # # ]: 0 : result = mdbImpl->tag_set_data( sense_tag, &reverse_set, 1, &dum_sense );
1198 [ # # ]: 0 : if( result != MB_SUCCESS ) return result;
1199 : : }
1200 : : }
1201 : :
1202 [ - + ]: 98 : if( mdbImpl->add_entities( ss_handle, ( entities_to_add.empty() ) ? NULL : &entities_to_add[0],
1203 [ - + ][ + - ]: 98 : entities_to_add.size() ) != MB_SUCCESS )
[ + - ]
1204 : 0 : return MB_FAILURE;
1205 : :
1206 : : // Distribution factor stuff
1207 [ + - ]: 49 : if( number_dist_factors_in_set )
1208 : : {
1209 : : // If this sideset does not already has a distribution factor array...set one
1210 : 49 : const void* ptr = 0;
1211 : 49 : int size = 0;
1212 [ + - ][ - + ]: 49 : if( MB_SUCCESS == mdbImpl->tag_get_by_ptr( mDistFactorTag, &ss_handle, 1, &ptr, &size ) )
1213 : : {
1214 : 0 : const double* data = reinterpret_cast< const double* >( ptr );
1215 [ # # ][ # # ]: 0 : std::copy( data, data + size, std::back_inserter( temp_dist_factor_vector ) );
1216 : : }
1217 : :
1218 [ + - ]: 49 : ptr = &temp_dist_factor_vector[0];
1219 : 49 : size = temp_dist_factor_vector.size();
1220 [ + - ][ - + ]: 49 : if( mdbImpl->tag_set_by_ptr( mDistFactorTag, &ss_handle, 1, &ptr, &size ) != MB_SUCCESS )
1221 [ + - ]: 49 : return MB_FAILURE;
1222 : : }
1223 : : }
1224 : 49 : }
1225 : :
1226 : 38 : return MB_SUCCESS;
1227 : : }
1228 : :
1229 : 49 : ErrorCode ReadNCDF::create_ss_elements( int* element_ids, int* side_list, int num_sides, int num_dist_factors,
1230 : : std::vector< EntityHandle >& entities_to_add,
1231 : : std::vector< EntityHandle >& reverse_entities,
1232 : : std::vector< double >& dist_factor_vector, int ss_seq_id )
1233 : : {
1234 : : // Determine entity type from element id
1235 : :
1236 : : // If there are dist. factors, create a vector to hold the array
1237 : : // and place this array as a tag onto the sideset meshset
1238 : :
1239 [ + - ]: 49 : std::vector< double > temp_dist_factor_vector( num_dist_factors );
1240 [ + - ]: 98 : std::vector< char > temp_string_storage( max_str_length + 1 );
1241 [ + - ]: 49 : char* temp_string = &temp_string_storage[0];
1242 : : int temp_var;
1243 [ + - ]: 49 : if( num_dist_factors )
1244 : : {
1245 : 49 : INS_ID( temp_string, "dist_fact_ss%d", ss_seq_id );
1246 : 49 : temp_var = -1;
1247 [ + - ][ + - ]: 49 : GET_1D_DBL_VAR( temp_string, temp_var, temp_dist_factor_vector );
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ - + ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ + - ][ + - ]
[ + - ][ - + ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ + - ]
[ + - ][ + - ]
[ - + ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ + - ]
1248 [ - + ][ # # ]: 49 : if( -1 == temp_var ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting dist fact variable" ); }
[ # # ][ # # ]
[ # # ][ # # ]
1249 : : }
1250 : :
1251 : : EntityType subtype;
1252 : : int num_side_nodes, num_elem_nodes;
1253 : : const EntityHandle* nodes;
1254 [ + - ]: 98 : std::vector< EntityHandle > connectivity;
1255 : : int side_node_idx[32];
1256 : :
1257 : 49 : int df_index = 0;
1258 : 49 : int sense = 0;
1259 [ + + ]: 422 : for( int i = 0; i < num_sides; i++ )
1260 : : {
1261 : : ExoIIElementType exoii_type;
1262 [ + - ]: 373 : ReadBlockData block_data;
1263 : 373 : block_data.elemType = EXOII_MAX_ELEM_TYPE;
1264 : :
1265 [ + - ][ - + ]: 373 : if( find_side_element_type( element_ids[i], exoii_type, block_data, df_index, side_list[i] ) != MB_SUCCESS )
1266 : 0 : continue; // Isn't being read in this time
1267 : :
1268 : 373 : EntityType type = ExoIIUtil::ExoIIElementMBEntity[exoii_type];
1269 : :
1270 : 373 : EntityHandle ent_handle = element_ids[i] - block_data.startExoId + block_data.startMBId;
1271 : :
1272 : 373 : const int side_num = side_list[i] - 1;
1273 [ + + ]: 373 : if( type == MBHEX )
1274 : : {
1275 : : // Get the nodes of the element
1276 [ + - ][ - + ]: 109 : if( mdbImpl->get_connectivity( ent_handle, nodes, num_elem_nodes ) != MB_SUCCESS ) return MB_FAILURE;
1277 : :
1278 [ + - ]: 109 : CN::SubEntityNodeIndices( type, num_elem_nodes, 2, side_num, subtype, num_side_nodes, side_node_idx );
1279 [ - + ]: 109 : if( num_side_nodes <= 0 ) return MB_FAILURE;
1280 : :
1281 [ + - ]: 109 : connectivity.resize( num_side_nodes );
1282 [ + + ]: 545 : for( int k = 0; k < num_side_nodes; ++k )
1283 [ + - ]: 436 : connectivity[k] = nodes[side_node_idx[k]];
1284 : :
1285 [ + - ][ - + ]: 109 : if( MB_SUCCESS != create_sideset_element( connectivity, subtype, ent_handle, sense ) ) return MB_FAILURE;
1286 [ + - ]: 109 : if( 1 == sense )
1287 [ + - ]: 109 : entities_to_add.push_back( ent_handle );
1288 : : else
1289 [ # # ]: 0 : reverse_entities.push_back( ent_handle );
1290 : :
1291 : : // Read in distribution factor array
1292 [ + - ]: 109 : if( num_dist_factors )
1293 : : {
1294 [ + + ]: 545 : for( int k = 0; k < 4; k++ )
1295 [ + - ][ + - ]: 436 : dist_factor_vector.push_back( temp_dist_factor_vector[df_index++] );
1296 : : }
1297 : : }
1298 : :
1299 : : // If it is a Tet
1300 [ + + ]: 264 : else if( type == MBTET )
1301 : : {
1302 : : // Get the nodes of the element
1303 [ + - ][ - + ]: 144 : if( mdbImpl->get_connectivity( ent_handle, nodes, num_elem_nodes ) != MB_SUCCESS ) return MB_FAILURE;
1304 : :
1305 [ + - ]: 144 : CN::SubEntityNodeIndices( type, num_elem_nodes, 2, side_num, subtype, num_side_nodes, side_node_idx );
1306 [ - + ]: 144 : if( num_side_nodes <= 0 ) return MB_FAILURE;
1307 : :
1308 [ + - ]: 144 : connectivity.resize( num_side_nodes );
1309 [ + + ]: 672 : for( int k = 0; k < num_side_nodes; ++k )
1310 [ + - ]: 528 : connectivity[k] = nodes[side_node_idx[k]];
1311 : :
1312 [ + - ][ - + ]: 144 : if( MB_SUCCESS != create_sideset_element( connectivity, subtype, ent_handle, sense ) ) return MB_FAILURE;
1313 [ + - ]: 144 : if( 1 == sense )
1314 [ + - ]: 144 : entities_to_add.push_back( ent_handle );
1315 : : else
1316 [ # # ]: 0 : reverse_entities.push_back( ent_handle );
1317 : :
1318 : : // Read in distribution factor array
1319 [ + - ]: 144 : if( num_dist_factors )
1320 : : {
1321 [ + + ]: 576 : for( int k = 0; k < 3; k++ )
1322 [ + - ][ + - ]: 432 : dist_factor_vector.push_back( temp_dist_factor_vector[df_index++] );
1323 : : }
1324 : : }
1325 [ + - ][ + + ]: 120 : else if( type == MBQUAD && exoii_type >= EXOII_SHELL && exoii_type <= EXOII_SHELL9 )
[ + - ]
1326 : : {
1327 : : // ent_handle = CREATE_HANDLE(MBQUAD, base_id, error);
1328 : :
1329 : : // Just use this quad
1330 [ - + ]: 72 : if( side_list[i] == 1 )
1331 : : {
1332 [ # # ]: 0 : if( 1 == sense )
1333 [ # # ]: 0 : entities_to_add.push_back( ent_handle );
1334 : : else
1335 [ # # ]: 0 : reverse_entities.push_back( ent_handle );
1336 : :
1337 [ # # ]: 0 : if( num_dist_factors )
1338 : : {
1339 [ # # ]: 0 : for( int k = 0; k < 4; k++ )
1340 [ # # ][ # # ]: 0 : dist_factor_vector.push_back( temp_dist_factor_vector[df_index++] );
1341 : : }
1342 : :
1343 : 0 : continue;
1344 : : }
1345 [ - + ]: 72 : else if( side_list[i] == 2 )
1346 : : {
1347 [ # # ]: 0 : reverse_entities.push_back( ent_handle );
1348 : :
1349 [ # # ]: 0 : if( num_dist_factors )
1350 : : {
1351 [ # # ]: 0 : for( int k = 0; k < 4; k++ )
1352 [ # # ][ # # ]: 0 : dist_factor_vector.push_back( temp_dist_factor_vector[df_index++] );
1353 : : }
1354 : :
1355 : 0 : continue;
1356 : : }
1357 : : else
1358 : : {
1359 : : // Get the nodes of the element
1360 [ + - ][ - + ]: 72 : if( mdbImpl->get_connectivity( ent_handle, nodes, num_elem_nodes ) != MB_SUCCESS ) return MB_FAILURE;
1361 : :
1362 : : CN::SubEntityNodeIndices( type, num_elem_nodes, 1, side_num - 2, subtype, num_side_nodes,
1363 [ + - ]: 72 : side_node_idx );
1364 [ - + ]: 72 : if( num_side_nodes <= 0 ) return MB_FAILURE;
1365 : :
1366 [ + - ]: 72 : connectivity.resize( num_side_nodes );
1367 [ + + ]: 216 : for( int k = 0; k < num_side_nodes; ++k )
1368 [ + - ]: 144 : connectivity[k] = nodes[side_node_idx[k]];
1369 : :
1370 [ + - ][ - + ]: 72 : if( MB_SUCCESS != create_sideset_element( connectivity, subtype, ent_handle, sense ) )
1371 : 0 : return MB_FAILURE;
1372 [ + - ]: 72 : if( 1 == sense )
1373 [ + - ]: 72 : entities_to_add.push_back( ent_handle );
1374 : : else
1375 [ # # ]: 0 : reverse_entities.push_back( ent_handle );
1376 : :
1377 [ + - ]: 72 : if( num_dist_factors )
1378 : : {
1379 [ + + ]: 216 : for( int k = 0; k < 2; k++ )
1380 [ + - ][ + - ]: 144 : dist_factor_vector.push_back( temp_dist_factor_vector[df_index++] );
1381 : : }
1382 : 72 : }
1383 : : }
1384 : : // If it is a Quad
1385 [ + - ]: 48 : else if( type == MBQUAD )
1386 : : {
1387 : : // Get the nodes of the element
1388 [ + - ][ - + ]: 48 : if( mdbImpl->get_connectivity( ent_handle, nodes, num_elem_nodes ) != MB_SUCCESS ) return MB_FAILURE;
1389 : :
1390 [ + - ]: 48 : CN::SubEntityNodeIndices( type, num_elem_nodes, 1, side_num, subtype, num_side_nodes, side_node_idx );
1391 [ - + ]: 48 : if( num_side_nodes <= 0 ) return MB_FAILURE;
1392 : :
1393 [ + - ]: 48 : connectivity.resize( num_side_nodes );
1394 [ + + ]: 144 : for( int k = 0; k < num_side_nodes; ++k )
1395 [ + - ]: 96 : connectivity[k] = nodes[side_node_idx[k]];
1396 : :
1397 [ + - ][ - + ]: 48 : if( MB_SUCCESS != create_sideset_element( connectivity, subtype, ent_handle, sense ) ) return MB_FAILURE;
1398 [ + - ]: 48 : if( 1 == sense )
1399 [ + - ]: 48 : entities_to_add.push_back( ent_handle );
1400 : : else
1401 [ # # ]: 0 : reverse_entities.push_back( ent_handle );
1402 : :
1403 : : // Read in distribution factor array
1404 [ + - ]: 48 : if( num_dist_factors )
1405 : : {
1406 [ + + ]: 144 : for( int k = 0; k < 2; k++ )
1407 [ + - ][ + - ]: 96 : dist_factor_vector.push_back( temp_dist_factor_vector[df_index++] );
1408 : : }
1409 : : }
1410 [ # # ]: 0 : else if( type == MBTRI )
1411 : : {
1412 : 0 : int side_offset = 0;
1413 [ # # ][ # # ]: 0 : if( number_dimensions() == 3 && side_list[i] <= 2 )
[ # # ][ # # ]
1414 : : {
1415 [ # # ]: 0 : if( 1 == sense )
1416 [ # # ]: 0 : entities_to_add.push_back( ent_handle );
1417 : : else
1418 [ # # ]: 0 : reverse_entities.push_back( ent_handle );
1419 [ # # ]: 0 : if( num_dist_factors )
1420 : : {
1421 [ # # ]: 0 : for( int k = 0; k < 3; k++ )
1422 [ # # ][ # # ]: 0 : dist_factor_vector.push_back( temp_dist_factor_vector[df_index++] );
1423 : : }
1424 : : }
1425 : : else
1426 : : {
1427 [ # # ][ # # ]: 0 : if( number_dimensions() == 3 )
1428 : : {
1429 [ # # ]: 0 : if( side_list[i] > 2 ) side_offset = 2;
1430 : : }
1431 : :
1432 : : // Get the nodes of the element
1433 [ # # ][ # # ]: 0 : if( mdbImpl->get_connectivity( ent_handle, nodes, num_elem_nodes ) != MB_SUCCESS ) return MB_FAILURE;
1434 : :
1435 : : CN::SubEntityNodeIndices( type, num_elem_nodes, 1, side_num - side_offset, subtype, num_side_nodes,
1436 [ # # ]: 0 : side_node_idx );
1437 [ # # ]: 0 : if( num_side_nodes <= 0 ) return MB_FAILURE;
1438 : :
1439 [ # # ]: 0 : connectivity.resize( num_side_nodes );
1440 [ # # ]: 0 : for( int k = 0; k < num_side_nodes; ++k )
1441 [ # # ]: 0 : connectivity[k] = nodes[side_node_idx[k]];
1442 : :
1443 [ # # ][ # # ]: 0 : if( MB_SUCCESS != create_sideset_element( connectivity, subtype, ent_handle, sense ) )
1444 : 0 : return MB_FAILURE;
1445 [ # # ]: 0 : if( 1 == sense )
1446 [ # # ]: 0 : entities_to_add.push_back( ent_handle );
1447 : : else
1448 [ # # ]: 0 : reverse_entities.push_back( ent_handle );
1449 : :
1450 [ # # ]: 0 : if( num_dist_factors )
1451 : : {
1452 [ # # ]: 373 : for( int k = 0; k < 2; k++ )
[ + - - ]
1453 [ # # ][ # # ]: 0 : dist_factor_vector.push_back( temp_dist_factor_vector[df_index++] );
1454 : : }
1455 : : }
1456 : : }
1457 : 373 : }
1458 : :
1459 : 98 : return MB_SUCCESS;
1460 : : }
1461 : :
1462 : 373 : ErrorCode ReadNCDF::create_sideset_element( const std::vector< EntityHandle >& connectivity, EntityType type,
1463 : : EntityHandle& handle, int& sense )
1464 : : {
1465 : : // Get adjacent entities
1466 : 373 : ErrorCode error = MB_SUCCESS;
1467 [ + - ]: 373 : int to_dim = CN::Dimension( type );
1468 : : // for matching purposes , consider only corner nodes
1469 : : // basically, assume everything is matching if the corner nodes are matching, and
1470 : : // the number of total nodes is the same
1471 [ + - ]: 373 : int nb_corner_nodes = CN::VerticesPerEntity( type );
1472 [ + - ]: 373 : std::vector< EntityHandle > adj_ent;
1473 [ + - ][ + - ]: 373 : mdbImpl->get_adjacencies( &( connectivity[0] ), 1, to_dim, false, adj_ent );
1474 : :
1475 : : // For each entity, see if we can find a match
1476 : : // If we find a match, return it
1477 : 373 : bool match_found = false;
1478 [ + - ]: 746 : std::vector< EntityHandle > match_conn;
1479 : : // By default, assume positive sense
1480 : 373 : sense = 1;
1481 [ + + ][ + + ]: 999 : for( unsigned int i = 0; i < adj_ent.size() && match_found == false; i++ )
[ + + ]
1482 : : {
1483 : : // Get the connectivity
1484 [ + - ][ + - ]: 626 : error = mdbImpl->get_connectivity( &( adj_ent[i] ), 1, match_conn );
1485 [ - + ]: 668 : if( error != MB_SUCCESS ) continue;
1486 : :
1487 : : // Make sure they have the same number of vertices (higher order elements ?)
1488 [ + + ]: 626 : if( match_conn.size() != connectivity.size() ) continue;
1489 : :
1490 : : // Find a matching node
1491 : 584 : std::vector< EntityHandle >::iterator iter;
1492 [ + - ]: 584 : std::vector< EntityHandle >::iterator end_corner = match_conn.begin() + nb_corner_nodes;
1493 [ + - ][ + - ]: 584 : iter = std::find( match_conn.begin(), end_corner, connectivity[0] );
1494 [ + - ][ - + ]: 584 : if( iter == end_corner ) continue;
1495 : :
1496 : : // Rotate to match connectivity
1497 : : // rotate only corner nodes, ignore the rest from now on
1498 [ + - ]: 584 : std::rotate( match_conn.begin(), iter, end_corner );
1499 : :
1500 : 584 : bool they_match = true;
1501 [ + + ]: 704 : for( int j = 1; j < nb_corner_nodes; j++ )
1502 : : {
1503 [ + - ][ + - ]: 656 : if( connectivity[j] != match_conn[j] )
[ + + ]
1504 : : {
1505 : 536 : they_match = false;
1506 : 536 : break;
1507 : : }
1508 : : }
1509 : :
1510 : : // If we didn't get a match
1511 [ + + ]: 584 : if( !they_match )
1512 : : {
1513 : : // Try the opposite sense
1514 : 536 : they_match = true;
1515 : :
1516 : 536 : int k = nb_corner_nodes - 1;
1517 [ + - ]: 664 : for( int j = 1; j < nb_corner_nodes; )
1518 : : {
1519 [ + - ][ + - ]: 664 : if( connectivity[j] != match_conn[k] )
[ + + ]
1520 : : {
1521 : 536 : they_match = false;
1522 : 536 : break;
1523 : : }
1524 : 128 : ++j;
1525 : 128 : --k;
1526 : : }
1527 : : // If they matched here, sense is reversed
1528 [ - + ]: 536 : if( they_match ) sense = -1;
1529 : : }
1530 : 584 : match_found = they_match;
1531 [ + + ][ + - ]: 584 : if( match_found ) handle = adj_ent[i];
1532 : : }
1533 : :
1534 : : // If we didn't find a match, create an element
1535 [ + + ][ + - ]: 373 : if( !match_found ) error = mdbImpl->create_element( type, &connectivity[0], connectivity.size(), handle );
[ + - ]
1536 : :
1537 : 373 : return error;
1538 : : }
1539 : :
1540 : 373 : ErrorCode ReadNCDF::find_side_element_type( const int exodus_id, ExoIIElementType& elem_type, ReadBlockData& block_data,
1541 : : int& df_index, int side_id )
1542 : : {
1543 : 373 : std::vector< ReadBlockData >::iterator iter, end_iter;
1544 : 373 : iter = blocksLoading.begin();
1545 : 373 : end_iter = blocksLoading.end();
1546 : 373 : elem_type = EXOII_MAX_ELEM_TYPE;
1547 : :
1548 [ + - ][ + - ]: 781 : for( ; iter != end_iter; ++iter )
[ + - ]
1549 : : {
1550 [ + - ][ + - ]: 781 : if( exodus_id >= ( *iter ).startExoId && exodus_id < ( *iter ).startExoId + ( *iter ).numElements )
[ + - ][ + - ]
[ + + ][ + + ]
1551 : : {
1552 [ + - ]: 373 : elem_type = ( *iter ).elemType;
1553 : :
1554 : : // If we're not reading this block in
1555 [ + - ][ - + ]: 373 : if( !( *iter ).reading_in )
1556 : : {
1557 : : // Offset df_indes according to type
1558 [ # # ][ # # ]: 0 : if( elem_type >= EXOII_HEX && elem_type <= EXOII_HEX27 )
1559 : 0 : df_index += 4;
1560 [ # # ][ # # ]: 0 : else if( elem_type >= EXOII_TETRA && elem_type <= EXOII_TETRA14 )
1561 : 0 : df_index += 3;
1562 [ # # ][ # # ]: 0 : else if( elem_type >= EXOII_QUAD && elem_type <= EXOII_QUAD9 )
1563 : 0 : df_index += 2;
1564 [ # # ][ # # ]: 0 : else if( elem_type >= EXOII_SHELL && elem_type <= EXOII_SHELL9 )
1565 : : {
1566 [ # # ][ # # ]: 0 : if( side_id == 1 || side_id == 2 )
1567 : 0 : df_index += 4;
1568 : : else
1569 : 0 : df_index += 2;
1570 : : }
1571 [ # # ][ # # ]: 0 : else if( elem_type >= EXOII_TRI && elem_type <= EXOII_TRI7 )
1572 : 0 : df_index += 3;
1573 : :
1574 : 0 : return MB_FAILURE;
1575 : : }
1576 : :
1577 [ + - ][ + - ]: 373 : block_data = *iter;
1578 : :
1579 : 373 : return MB_SUCCESS;
1580 : : }
1581 : : }
1582 : 373 : return MB_FAILURE;
1583 : : }
1584 : :
1585 : 2 : ErrorCode ReadNCDF::read_qa_records( EntityHandle file_set )
1586 : : {
1587 [ + - ]: 2 : std::vector< std::string > qa_records;
1588 [ + - ]: 2 : read_qa_information( qa_records );
1589 : :
1590 [ + - ]: 4 : std::vector< char > tag_data;
1591 [ + - ][ + - ]: 10 : for( std::vector< std::string >::iterator i = qa_records.begin(); i != qa_records.end(); ++i )
[ + + ]
1592 : : {
1593 [ + - ][ + - ]: 8 : std::copy( i->begin(), i->end(), std::back_inserter( tag_data ) );
[ + - ][ + - ]
1594 [ + - ]: 8 : tag_data.push_back( '\0' );
1595 : : }
1596 : :
1597 : : // If there were qa_records...tag them to the mCurrentMeshHandle
1598 [ + - ]: 2 : if( !tag_data.empty() )
1599 : : {
1600 [ + - ]: 2 : const void* ptr = &tag_data[0];
1601 : 2 : int size = tag_data.size();
1602 [ + - ][ - + ]: 2 : if( mdbImpl->tag_set_by_ptr( mQaRecordTag, &file_set, 1, &ptr, &size ) != MB_SUCCESS ) { return MB_FAILURE; }
1603 : : }
1604 : :
1605 : 4 : return MB_SUCCESS;
1606 : : }
1607 : :
1608 : 2 : ErrorCode ReadNCDF::read_qa_information( std::vector< std::string >& qa_record_list )
1609 : : {
1610 : : // Inquire on the genesis file to find the number of qa records
1611 : :
1612 : 2 : int number_records = 0;
1613 : :
1614 : : int temp_dim;
1615 [ + - ][ + - ]: 2 : GET_DIM( temp_dim, "num_qa_rec", number_records );
[ + - ][ - + ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1616 [ + - ]: 2 : std::vector< char > data( max_str_length + 1 );
1617 : :
1618 [ + + ]: 4 : for( int i = 0; i < number_records; i++ )
1619 : : {
1620 [ + + ]: 10 : for( int j = 0; j < 4; j++ )
1621 : : {
1622 [ + - ]: 8 : data[max_str_length] = '\0';
1623 [ + - ][ + - ]: 8 : if( read_qa_string( &data[0], i, j ) != MB_SUCCESS ) return MB_FAILURE;
[ - + ]
1624 : :
1625 [ + - ][ + - ]: 8 : qa_record_list.push_back( &data[0] );
[ + - ]
1626 : : }
1627 : : }
1628 : :
1629 : 2 : return MB_SUCCESS;
1630 : : }
1631 : :
1632 : 8 : ErrorCode ReadNCDF::read_qa_string( char* temp_string, int record_number, int record_position )
1633 : : {
1634 [ + - ]: 8 : std::vector< int > dims;
1635 : 8 : int temp_var = -1;
1636 [ + - ][ + - ]: 8 : GET_VAR( "qa_records", temp_var, dims );
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ - + ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1637 [ - + ]: 8 : if( -1 == temp_var )
1638 : : {
1639 [ # # ][ # # ]: 0 : MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting qa record variable" );
[ # # ][ # # ]
[ # # ]
1640 : : return MB_FAILURE;
1641 : : }
1642 : : size_t count[3], start[3];
1643 : 8 : start[0] = record_number;
1644 : 8 : start[1] = record_position;
1645 : 8 : start[2] = 0;
1646 : 8 : count[0] = 1;
1647 : 8 : count[1] = 1;
1648 : 8 : count[2] = max_str_length;
1649 [ + - ]: 8 : int fail = nc_get_vara_text( ncFile, temp_var, start, count, temp_string );
1650 [ - + ]: 8 : if( NC_NOERR != fail )
1651 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem setting current record number variable position" ); }
[ # # ][ # # ]
[ # # ]
1652 : : // Get the variable id in the exodus file
1653 : :
1654 : 8 : return MB_SUCCESS;
1655 : : }
1656 : :
1657 : : // The cub_file_set contains the mesh to be updated. There could exist other
1658 : : // file sets that should be kept separate, such as the geometry file set from
1659 : : // ReadCGM.
1660 : 0 : ErrorCode ReadNCDF::update( const char* exodus_file_name, const FileOptions& opts, const int num_blocks,
1661 : : const int* blocks_to_load, const EntityHandle cub_file_set )
1662 : : {
1663 : : // Function : updating current database from new exodus_file.
1664 : : // Creator: Jane Hu
1665 : : // opts is currently designed as following
1666 : : // tdata = <var_name>[, time][,op][,destination]
1667 : : // where var_name show the tag name to be updated, this version just takes
1668 : : // coord.
1669 : : // time is the optional, and it gives time step of each of the mesh
1670 : : // info in exodus file. It start from 1.
1671 : : // op is the operation that is going to be performed on the var_name info.
1672 : : // currently support 'set'
1673 : : // destination shows where to store the updated info, currently assume it is
1674 : : // stored in the same database by replacing the old info if there's no input
1675 : : // for destination, or the destination data is given in exodus format and just
1676 : : // need to update the coordinates.
1677 : : // Assumptions:
1678 : : // 1. Assume the num_el_blk's in both DB and update exodus file are the same.
1679 : : // 2. Assume num_el_in_blk1...num_el_in_blk(num_el_blk) numbers are matching, may in
1680 : : // different order. example: num_el_in_blk11 = num_el_in_blk22 && num_el_in_blk12 =
1681 : : // num_el_in_blk21.
1682 : : // 3. In exodus file, get node_num_map
1683 : : // 4. loop through the node_num_map, use it to find the node in the cub file.
1684 : : // 5. Replace coord[0][n] with coordx[m] + vals_nod_var1(time_step, m) for all directions for
1685 : : // matching nodes.
1686 : : // Test: try to get hexes
1687 : :
1688 : : // *******************************************************************
1689 : : // Move nodes to their deformed locations.
1690 : : // *******************************************************************
1691 : : ErrorCode rval;
1692 [ # # ]: 0 : std::string s;
1693 [ # # ]: 0 : rval = opts.get_str_option( "tdata", s );
1694 [ # # ][ # # ]: 0 : if( MB_SUCCESS != rval ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem reading file options" ); }
[ # # ][ # # ]
[ # # ][ # # ]
1695 [ # # ]: 0 : std::vector< std::string > tokens;
1696 [ # # ]: 0 : tokenize( s, tokens, "," );
1697 : :
1698 : : // 1. Check for time step to find the match time
1699 : 0 : int time_step = 1;
1700 [ # # ][ # # ]: 0 : if( tokens.size() > 1 && !tokens[1].empty() )
[ # # ][ # # ]
1701 : : {
1702 [ # # ]: 0 : const char* time_s = tokens[1].c_str();
1703 : : char* endptr;
1704 : 0 : long int pval = strtol( time_s, &endptr, 0 );
1705 [ # # ]: 0 : std::string end = endptr;
1706 [ # # ]: 0 : if( !end.empty() ) // Syntax error
1707 : 0 : return MB_TYPE_OUT_OF_RANGE;
1708 : :
1709 : : // Check for overflow (parsing long int, returning int)
1710 : 0 : time_step = pval;
1711 [ # # ]: 0 : if( pval != (long int)time_step ) return MB_TYPE_OUT_OF_RANGE;
1712 [ # # ][ # # ]: 0 : if( time_step <= 0 ) return MB_TYPE_OUT_OF_RANGE;
1713 : : }
1714 : :
1715 : : // 2. Check for the operations, currently support set.
1716 : : const char* op;
1717 [ # # ][ # # ]: 0 : if( tokens.size() < 3 || ( tokens[2] != "set" && tokens[2] != "add" ) )
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
1718 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_TYPE_OUT_OF_RANGE, "ReadNCDF: invalid operation specified for update" ); }
[ # # ][ # # ]
[ # # ]
1719 [ # # ]: 0 : op = tokens[2].c_str();
1720 : :
1721 : : // Open netcdf/exodus file
1722 : 0 : ncFile = 0;
1723 [ # # ]: 0 : int fail = nc_open( exodus_file_name, 0, &ncFile );
1724 [ # # ]: 0 : if( !ncFile )
1725 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FILE_DOES_NOT_EXIST, "ReadNCDF:: problem opening Netcdf/Exodus II file " << exodus_file_name ); }
[ # # ][ # # ]
[ # # ][ # # ]
1726 : :
1727 [ # # ]: 0 : rval = read_exodus_header();
1728 [ # # ]: 0 : if( MB_SUCCESS != rval ) return rval;
1729 : :
1730 : : // Check to make sure that the requested time step exists
1731 : 0 : int ncdim = -1;
1732 : : int max_time_steps;
1733 [ # # ][ # # ]: 0 : GET_DIM( ncdim, "time_step", max_time_steps );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1734 [ # # ]: 0 : if( -1 == ncdim )
1735 : : {
1736 [ # # ][ # # ]: 0 : std::cout << "ReadNCDF: could not get number of time steps" << std::endl;
1737 : 0 : return MB_FAILURE;
1738 : : }
1739 [ # # ][ # # ]: 0 : std::cout << " Maximum time step=" << max_time_steps << std::endl;
[ # # ]
1740 [ # # ]: 0 : if( max_time_steps < time_step )
1741 : : {
1742 [ # # ][ # # ]: 0 : std::cout << "ReadNCDF: time step is greater than max_time_steps" << std::endl;
1743 : 0 : return MB_FAILURE;
1744 : : }
1745 : :
1746 : : // Get the time
1747 [ # # ]: 0 : std::vector< double > times( max_time_steps );
1748 : 0 : int nc_var = -1;
1749 [ # # ][ # # ]: 0 : GET_1D_DBL_VAR( "time_whole", nc_var, times );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
1750 [ # # ][ # # ]: 0 : if( -1 == nc_var ) { std::cout << "ReadNCDF: unable to get time variable" << std::endl; }
[ # # ]
1751 : : else
1752 : : {
1753 [ # # ][ # # ]: 0 : std::cout << " Step " << time_step << " is at " << times[time_step - 1] << " seconds" << std::endl;
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1754 : : }
1755 : :
1756 : : // Read in the node_num_map.
1757 [ # # ]: 0 : std::vector< int > ptr( numberNodes_loading );
1758 : :
1759 : 0 : int varid = -1;
1760 [ # # ][ # # ]: 0 : GET_1D_INT_VAR( "node_num_map", varid, ptr );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
1761 [ # # ][ # # ]: 0 : if( -1 == varid ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting node number map data" ); }
[ # # ][ # # ]
[ # # ][ # # ]
1762 : :
1763 : : // Read in the deformations.
1764 [ # # ]: 0 : std::vector< std::vector< double > > deformed_arrays( 3 );
1765 [ # # ]: 0 : std::vector< std::vector< double > > orig_coords( 3 );
1766 [ # # ][ # # ]: 0 : deformed_arrays[0].reserve( numberNodes_loading );
1767 [ # # ][ # # ]: 0 : deformed_arrays[1].reserve( numberNodes_loading );
1768 [ # # ][ # # ]: 0 : deformed_arrays[2].reserve( numberNodes_loading );
1769 [ # # ][ # # ]: 0 : orig_coords[0].reserve( numberNodes_loading );
1770 [ # # ][ # # ]: 0 : orig_coords[1].reserve( numberNodes_loading );
1771 [ # # ][ # # ]: 0 : orig_coords[2].reserve( numberNodes_loading );
1772 : 0 : size_t start[2] = { static_cast< size_t >( time_step - 1 ), 0 },
1773 : 0 : count[2] = { 1, static_cast< size_t >( numberNodes_loading ) };
1774 [ # # ]: 0 : std::vector< int > dims;
1775 : 0 : int coordx = -1, coordy = -1, coordz = -1;
1776 [ # # ][ # # ]: 0 : GET_VAR( "vals_nod_var1", coordx, dims );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1777 [ # # ][ # # ]: 0 : GET_VAR( "vals_nod_var2", coordy, dims );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1778 [ # # ][ # # ]: 0 : if( numberDimensions_loading == 3 ) GET_VAR( "vals_nod_var3", coordz, dims );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
1779 [ # # ][ # # ]: 0 : if( -1 == coordx || -1 == coordy || ( numberDimensions_loading == 3 && -1 == coordz ) )
[ # # ][ # # ]
1780 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting coords variable" ); }
[ # # ][ # # ]
[ # # ]
1781 : :
1782 [ # # ][ # # ]: 0 : fail = nc_get_vara_double( ncFile, coordx, start, count, &deformed_arrays[0][0] );
[ # # ]
1783 [ # # ][ # # ]: 0 : if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting x deformation array" ); }
[ # # ][ # # ]
[ # # ][ # # ]
1784 : :
1785 [ # # ][ # # ]: 0 : fail = nc_get_vara_double( ncFile, coordy, start, count, &deformed_arrays[1][0] );
[ # # ]
1786 [ # # ][ # # ]: 0 : if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting y deformation array" ); }
[ # # ][ # # ]
[ # # ][ # # ]
1787 [ # # ]: 0 : if( numberDimensions_loading == 3 )
1788 : : {
1789 [ # # ][ # # ]: 0 : fail = nc_get_vara_double( ncFile, coordz, start, count, &deformed_arrays[2][0] );
[ # # ]
1790 [ # # ][ # # ]: 0 : if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting z deformation array" ); }
[ # # ][ # # ]
[ # # ][ # # ]
1791 : : }
1792 : :
1793 : 0 : int coord1 = -1, coord2 = -1, coord3 = -1;
1794 [ # # ][ # # ]: 0 : GET_1D_DBL_VAR( "coordx", coord1, orig_coords[0] );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
1795 [ # # ][ # # ]: 0 : if( -1 == coord1 ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting x coord array" ); }
[ # # ][ # # ]
[ # # ][ # # ]
1796 [ # # ][ # # ]: 0 : GET_1D_DBL_VAR( "coordy", coord2, orig_coords[1] );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
1797 [ # # ][ # # ]: 0 : if( -1 == coord2 ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting y coord array" ); }
[ # # ][ # # ]
[ # # ][ # # ]
1798 [ # # ]: 0 : if( numberDimensions_loading == 3 )
1799 : : {
1800 [ # # ][ # # ]: 0 : GET_1D_DBL_VAR( "coordz", coord3, orig_coords[2] );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
1801 [ # # ][ # # ]: 0 : if( -1 == coord3 ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting z coord array" ); }
[ # # ][ # # ]
[ # # ][ # # ]
1802 : : }
1803 : :
1804 : : // b. Deal with DB file : get node info. according to node_num_map.
1805 [ # # ][ # # ]: 0 : if( tokens[0] != "coord" && tokens[0] != "COORD" ) return MB_NOT_IMPLEMENTED;
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1806 : :
1807 [ # # ][ # # ]: 0 : if( strcmp( op, "set" ) && strcmp( op, " set" ) ) return MB_NOT_IMPLEMENTED;
1808 : :
1809 : : // Two methods of matching nodes (id vs. proximity)
1810 : 0 : const bool match_node_ids = true;
1811 : :
1812 : : // Get nodes in cubit file
1813 [ # # ]: 0 : Range cub_verts;
1814 [ # # ]: 0 : rval = mdbImpl->get_entities_by_type( cub_file_set, MBVERTEX, cub_verts );
1815 [ # # ]: 0 : if( MB_SUCCESS != rval ) return rval;
1816 [ # # ][ # # ]: 0 : std::cout << " cub_file_set contains " << cub_verts.size() << " nodes." << std::endl;
[ # # ][ # # ]
[ # # ]
1817 : :
1818 : : // Some accounting
1819 [ # # ][ # # ]: 0 : std::cout << " exodus file contains " << numberNodes_loading << " nodes." << std::endl;
[ # # ][ # # ]
1820 : 0 : double max_magnitude = 0;
1821 : 0 : double average_magnitude = 0;
1822 : 0 : int found = 0;
1823 : 0 : int lost = 0;
1824 [ # # ]: 0 : std::map< int, EntityHandle > cub_verts_id_map;
1825 [ # # ]: 0 : AdaptiveKDTree kdtree( mdbImpl );
1826 : : EntityHandle root;
1827 : :
1828 : : // Should not use cub verts unless they have been matched. Place in a map
1829 : : // for fast handle_by_id lookup.
1830 [ # # ]: 0 : std::map< int, EntityHandle > matched_cub_vert_id_map;
1831 : :
1832 : : // Place cub verts in a map for searching by id
1833 : : if( match_node_ids )
1834 : : {
1835 [ # # ][ # # ]: 0 : std::vector< int > cub_ids( cub_verts.size() );
1836 [ # # ][ # # ]: 0 : rval = mdbImpl->tag_get_data( mGlobalIdTag, cub_verts, &cub_ids[0] );
1837 [ # # ]: 0 : if( MB_SUCCESS != rval ) return rval;
1838 [ # # ][ # # ]: 0 : for( unsigned i = 0; i != cub_verts.size(); ++i )
[ # # ]
1839 : : {
1840 [ # # ][ # # ]: 0 : cub_verts_id_map.insert( std::pair< int, EntityHandle >( cub_ids[i], cub_verts[i] ) );
[ # # ][ # # ]
1841 : 0 : }
1842 : :
1843 : : // Place cub verts in a kdtree for searching by proximity
1844 : : }
1845 : : else
1846 : : {
1847 : : FileOptions tree_opts( "MAX_PER_LEAF=1;SPLITS_PER_DIR=1;CANDIDATE_PLANE_SET=0" );
1848 : : rval = kdtree.build_tree( cub_verts, &root, &tree_opts );
1849 : : if( MB_SUCCESS != rval ) return rval;
1850 : : AdaptiveKDTreeIter tree_iter;
1851 : : rval = kdtree.get_tree_iterator( root, tree_iter );
1852 : : if( MB_SUCCESS != rval ) return rval;
1853 : : }
1854 : :
1855 : : // For each exo vert, find the matching cub vert
1856 [ # # ]: 0 : for( int i = 0; i < numberNodes_loading; ++i )
1857 : : {
1858 [ # # ]: 0 : int exo_id = ptr[i];
1859 [ # # ][ # # ]: 0 : CartVect exo_coords( orig_coords[0][i], orig_coords[1][i], orig_coords[2][i] );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1860 : 0 : EntityHandle cub_vert = 0;
1861 : 0 : bool found_match = false;
1862 : :
1863 : : // By id
1864 : : if( match_node_ids )
1865 : : {
1866 [ # # ]: 0 : std::map< int, EntityHandle >::iterator i_iter;
1867 [ # # ]: 0 : i_iter = cub_verts_id_map.find( exo_id );
1868 [ # # ][ # # ]: 0 : if( i_iter != cub_verts_id_map.end() )
1869 : : {
1870 : 0 : found_match = true;
1871 [ # # ]: 0 : cub_vert = i_iter->second;
1872 : : }
1873 : :
1874 : : // By proximity
1875 : : }
1876 : : else
1877 : : {
1878 : : // The MAX_NODE_DIST is the farthest distance to search for a node.
1879 : : // For the 1/12th symmetry 85 pin model, the max node dist could not be less
1880 : : // than 1e-1 (March 26, 2010).
1881 : : const double MAX_NODE_DIST = 1e-1;
1882 : :
1883 : : std::vector< EntityHandle > leaves;
1884 : : double min_dist = MAX_NODE_DIST;
1885 : : rval = kdtree.distance_search( exo_coords.array(), MAX_NODE_DIST, leaves );
1886 : : if( MB_SUCCESS != rval ) return rval;
1887 : : for( std::vector< EntityHandle >::const_iterator j = leaves.begin(); j != leaves.end(); ++j )
1888 : : {
1889 : : std::vector< EntityHandle > leaf_verts;
1890 : : rval = mdbImpl->get_entities_by_type( *j, MBVERTEX, leaf_verts );
1891 : : if( MB_SUCCESS != rval ) return rval;
1892 : : for( std::vector< EntityHandle >::const_iterator k = leaf_verts.begin(); k != leaf_verts.end(); ++k )
1893 : : {
1894 : : CartVect orig_cub_coords, difference;
1895 : : rval = mdbImpl->get_coords( &( *k ), 1, orig_cub_coords.array() );
1896 : : if( MB_SUCCESS != rval ) return rval;
1897 : : difference = orig_cub_coords - exo_coords;
1898 : : double dist = difference.length();
1899 : : if( dist < min_dist )
1900 : : {
1901 : : min_dist = dist;
1902 : : cub_vert = *k;
1903 : : }
1904 : : }
1905 : : }
1906 : : if( 0 != cub_vert ) found_match = true;
1907 : : }
1908 : :
1909 : : // If a match is found, update it with the deformed coords from the exo file.
1910 [ # # ]: 0 : if( found_match )
1911 : : {
1912 [ # # ]: 0 : CartVect updated_exo_coords;
1913 [ # # ][ # # ]: 0 : matched_cub_vert_id_map.insert( std::pair< int, EntityHandle >( exo_id, cub_vert ) );
1914 [ # # ][ # # ]: 0 : updated_exo_coords[0] = orig_coords[0][i] + deformed_arrays[0][i];
[ # # ][ # # ]
[ # # ]
1915 [ # # ][ # # ]: 0 : updated_exo_coords[1] = orig_coords[1][i] + deformed_arrays[1][i];
[ # # ][ # # ]
[ # # ]
1916 [ # # ][ # # ]: 0 : if( numberDimensions_loading == 3 ) updated_exo_coords[2] = orig_coords[2][i] + deformed_arrays[2][i];
[ # # ][ # # ]
[ # # ][ # # ]
1917 [ # # ][ # # ]: 0 : rval = mdbImpl->set_coords( &cub_vert, 1, updated_exo_coords.array() );
1918 [ # # ]: 0 : if( MB_SUCCESS != rval ) return rval;
1919 : 0 : ++found;
1920 : : double magnitude =
1921 [ # # ][ # # ]: 0 : sqrt( deformed_arrays[0][i] * deformed_arrays[0][i] + deformed_arrays[1][i] * deformed_arrays[1][i] +
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
1922 [ # # ][ # # ]: 0 : deformed_arrays[2][i] * deformed_arrays[2][i] );
[ # # ][ # # ]
1923 [ # # ]: 0 : if( magnitude > max_magnitude ) max_magnitude = magnitude;
1924 : 0 : average_magnitude += magnitude;
1925 : : }
1926 : : else
1927 : : {
1928 : 0 : ++lost;
1929 [ # # ][ # # ]: 0 : std::cout << "cannot match exo node " << exo_id << " " << exo_coords << std::endl;
[ # # ][ # # ]
[ # # ]
1930 : : }
1931 : : }
1932 : :
1933 : : // Exo verts that could not be matched have already been printed. Now print the
1934 : : // cub verts that could not be matched.
1935 [ # # ][ # # ]: 0 : if( matched_cub_vert_id_map.size() < cub_verts.size() )
1936 : : {
1937 [ # # ]: 0 : Range unmatched_cub_verts = cub_verts;
1938 [ # # ][ # # ]: 0 : for( std::map< int, EntityHandle >::const_iterator i = matched_cub_vert_id_map.begin();
[ # # ][ # # ]
1939 [ # # ]: 0 : i != matched_cub_vert_id_map.end(); ++i )
1940 : : {
1941 [ # # ][ # # ]: 0 : unmatched_cub_verts.erase( i->second );
1942 : : }
1943 [ # # ][ # # ]: 0 : for( Range::const_iterator i = unmatched_cub_verts.begin(); i != unmatched_cub_verts.end(); ++i )
[ # # ][ # # ]
[ # # ]
1944 : : {
1945 : : int cub_id;
1946 [ # # ][ # # ]: 0 : rval = mdbImpl->tag_get_data( mGlobalIdTag, &( *i ), 1, &cub_id );
1947 [ # # ]: 0 : if( MB_SUCCESS != rval ) return rval;
1948 [ # # ]: 0 : CartVect cub_coords;
1949 [ # # ][ # # ]: 0 : rval = mdbImpl->get_coords( &( *i ), 1, cub_coords.array() );
[ # # ]
1950 [ # # ]: 0 : if( MB_SUCCESS != rval ) return rval;
1951 [ # # ][ # # ]: 0 : std::cout << "cannot match cub node " << cub_id << " " << cub_coords << std::endl;
[ # # ][ # # ]
[ # # ]
1952 : : }
1953 [ # # ][ # # ]: 0 : std::cout << " " << unmatched_cub_verts.size() << " nodes from the cub file could not be matched."
[ # # ][ # # ]
1954 [ # # ][ # # ]: 0 : << std::endl;
1955 : : }
1956 : :
1957 : : // Summarize statistics
1958 [ # # ][ # # ]: 0 : std::cout << " " << found << " nodes from the exodus file were matched in the cub_file_set ";
[ # # ]
1959 [ # # ][ # # ]: 0 : if( match_node_ids ) { std::cout << "by id." << std::endl; }
1960 : : else
1961 : : {
1962 : : std::cout << "by proximity." << std::endl;
1963 : : }
1964 : :
1965 : : // Fail if all of the nodes could not be matched.
1966 [ # # ]: 0 : if( 0 != lost )
1967 : : {
1968 [ # # ][ # # ]: 0 : std::cout << "Error: " << lost << " nodes from the exodus file could not be matched." << std::endl;
[ # # ][ # # ]
1969 : : // return MB_FAILURE;
1970 : : }
1971 [ # # ][ # # ]: 0 : std::cout << " maximum node displacement magnitude: " << max_magnitude << " cm" << std::endl;
[ # # ][ # # ]
1972 [ # # ][ # # ]: 0 : std::cout << " average node displacement magnitude: " << average_magnitude / found << " cm" << std::endl;
[ # # ][ # # ]
1973 : :
1974 : : // *******************************************************************
1975 : : // Remove dead elements from the MOAB instance.
1976 : : // *******************************************************************
1977 : :
1978 : : // How many element variables are in the file?
1979 : : int n_elem_var;
1980 [ # # ][ # # ]: 0 : GET_DIM( ncdim, "num_elem_var", n_elem_var );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1981 : :
1982 : : // Get element variable names
1983 : 0 : varid = -1;
1984 [ # # ]: 0 : int cstatus = nc_inq_varid( ncFile, "name_elem_var", &varid );
1985 [ # # ]: 0 : std::vector< char > names_memory( n_elem_var * max_str_length );
1986 [ # # ]: 0 : std::vector< char* > names( n_elem_var );
1987 [ # # ]: 0 : for( int i = 0; i < n_elem_var; ++i )
1988 [ # # ][ # # ]: 0 : names[i] = &names_memory[i * max_str_length];
1989 [ # # ][ # # ]: 0 : if( cstatus != NC_NOERR || varid == -1 )
1990 : : {
1991 [ # # ][ # # ]: 0 : std::cout << "ReadNCDF: name_elem_var does not exist" << std::endl;
1992 : 0 : return MB_FAILURE;
1993 : : }
1994 [ # # ][ # # ]: 0 : int status = nc_get_var_text( ncFile, varid, &names_memory[0] );
1995 [ # # ][ # # ]: 0 : if( NC_NOERR != status ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF: Problem getting element variable names" ); }
[ # # ][ # # ]
[ # # ][ # # ]
1996 : :
1997 : : // Is one of the element variable names "death_status"? If so, get its index
1998 : : // in the element variable array.
1999 : : int death_index;
2000 : 0 : bool found_death_index = false;
2001 [ # # ]: 0 : for( int i = 0; i < n_elem_var; ++i )
2002 : : {
2003 [ # # ][ # # ]: 0 : std::string temp( names[i] );
2004 [ # # ]: 0 : std::string::size_type pos0 = temp.find( "death_status" );
2005 [ # # ]: 0 : std::string::size_type pos1 = temp.find( "Death_Status" );
2006 [ # # ]: 0 : std::string::size_type pos2 = temp.find( "DEATH_STATUS" );
2007 [ # # ][ # # ]: 0 : if( std::string::npos != pos0 || std::string::npos != pos1 || std::string::npos != pos2 )
[ # # ]
2008 : : {
2009 : 0 : found_death_index = true;
2010 : 0 : death_index = i + 1; // NetCDF variables start with 1
2011 [ # # ]: 0 : break;
2012 : : }
2013 : 0 : }
2014 [ # # ][ # # ]: 0 : if( !found_death_index ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF: Problem getting index of death_status variable" ); }
[ # # ][ # # ]
[ # # ][ # # ]
2015 : :
2016 : : // The exodus header has already been read. This contains the number of element
2017 : : // blocks.
2018 : :
2019 : : // Dead elements are listed by block. Read the block headers to determine how
2020 : : // many elements are in each block.
2021 [ # # ]: 0 : rval = read_block_headers( blocks_to_load, num_blocks );
2022 [ # # ][ # # ]: 0 : if( MB_FAILURE == rval ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF: Problem reading block headers" ); }
[ # # ][ # # ]
[ # # ][ # # ]
2023 : :
2024 : : // Dead elements from the Exodus file can be located in the cub_file_set by id
2025 : : // or by connectivity. Currently, finding elements by id requires careful book
2026 : : // keeping when constructing the model in Cubit. To avoid this, one can match
2027 : : // elements by connectivity instead.
2028 : 0 : const bool match_elems_by_connectivity = true;
2029 : :
2030 : : // Get the element id map. The ids in the map are from the elements in the blocks.
2031 : : // elem_num_map(blk1 elem ids, blk2 elem ids, blk3 elem ids, ...)
2032 [ # # ]: 0 : std::vector< int > elem_ids( numberNodes_loading );
2033 : : if( !match_elems_by_connectivity )
2034 : : {
2035 : : GET_1D_INT_VAR( "elem_num_map", varid, elem_ids );
2036 : : if( -1 == varid ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF: Problem getting element number map data" ); }
2037 : : }
2038 : :
2039 : : // For each block
2040 : 0 : int first_elem_id_in_block = 0;
2041 : 0 : int block_count = 1; // NetCDF variables start with 1
2042 : 0 : int total_elems = 0;
2043 : 0 : int total_dead_elems = 0;
2044 [ # # ][ # # ]: 0 : for( std::vector< ReadBlockData >::iterator i = blocksLoading.begin(); i != blocksLoading.end(); ++i )
[ # # ]
2045 : : {
2046 : : // Get the ncdf connect variable
2047 [ # # ]: 0 : std::string temp_string( "connect" );
2048 [ # # ][ # # ]: 0 : std::stringstream temp_ss;
[ # # ]
2049 [ # # ]: 0 : temp_ss << block_count;
2050 [ # # ][ # # ]: 0 : temp_string += temp_ss.str();
2051 [ # # ]: 0 : temp_string += "\0";
2052 [ # # ][ # # ]: 0 : std::vector< int > ldims;
2053 [ # # ][ # # ]: 0 : GET_VAR( temp_string.c_str(), nc_var, ldims );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
2054 [ # # ][ # # ]: 0 : if( !nc_var ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting connect variable" ); }
[ # # ][ # # ]
[ # # ][ # # ]
2055 : : // The element type is an attribute of the connectivity variable
2056 : : nc_type att_type;
2057 : : size_t att_len;
2058 [ # # ]: 0 : fail = nc_inq_att( ncFile, nc_var, "elem_type", &att_type, &att_len );
2059 [ # # ][ # # ]: 0 : if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting elem type attribute" ); }
[ # # ][ # # ]
[ # # ][ # # ]
2060 [ # # ][ # # ]: 0 : std::vector< char > dum_str( att_len + 1 );
2061 [ # # ][ # # ]: 0 : fail = nc_get_att_text( ncFile, nc_var, "elem_type", &dum_str[0] );
2062 [ # # ][ # # ]: 0 : if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting elem type" ); }
[ # # ][ # # ]
[ # # ][ # # ]
2063 [ # # ][ # # ]: 0 : ExoIIElementType elem_type = ExoIIUtil::static_element_name_to_type( &dum_str[0] );
2064 [ # # ]: 0 : ( *i ).elemType = elem_type;
2065 [ # # ]: 0 : const EntityType mb_type = ExoIIUtil::ExoIIElementMBEntity[( *i ).elemType];
2066 : :
2067 : : // Get the number of nodes per element
2068 [ # # ]: 0 : unsigned int nodes_per_element = ExoIIUtil::VerticesPerElement[( *i ).elemType];
2069 : :
2070 : : // Read the connectivity into that memory.
2071 : : // int exo_conn[i->numElements][nodes_per_element];
2072 [ # # ][ # # ]: 0 : int* exo_conn = new int[i->numElements * nodes_per_element];
[ # # ]
2073 : : // NcBool status = temp_var->get(&exo_conn[0][0], i->numElements, nodes_per_element);
2074 : 0 : size_t lstart[2] = { 0, 0 }, lcount[2];
2075 [ # # ]: 0 : lcount[0] = i->numElements;
2076 : 0 : lcount[1] = nodes_per_element;
2077 [ # # ]: 0 : fail = nc_get_vara_int( ncFile, nc_var, lstart, lcount, exo_conn );
2078 [ # # ]: 0 : if( NC_NOERR != fail )
2079 : : {
2080 [ # # ]: 0 : delete[] exo_conn;
2081 [ # # ][ # # ]: 0 : MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting connectivity" );
[ # # ][ # # ]
[ # # ]
2082 : : }
2083 : :
2084 : : // Get the death_status at the correct time step.
2085 [ # # ][ # # ]: 0 : std::vector< double > death_status( i->numElements ); // It seems wrong, but it uses doubles
[ # # ]
2086 [ # # ][ # # ]: 0 : std::string array_name( "vals_elem_var" );
2087 [ # # ][ # # ]: 0 : temp_ss.str( "" ); // stringstream won't clear by temp.clear()
2088 [ # # ]: 0 : temp_ss << death_index;
2089 [ # # ][ # # ]: 0 : array_name += temp_ss.str();
2090 [ # # ]: 0 : array_name += "eb";
2091 [ # # ][ # # ]: 0 : temp_ss.str( "" ); // stringstream won't clear by temp.clear()
2092 [ # # ]: 0 : temp_ss << block_count;
2093 [ # # ][ # # ]: 0 : array_name += temp_ss.str();
2094 [ # # ]: 0 : array_name += "\0";
2095 [ # # ][ # # ]: 0 : GET_VAR( array_name.c_str(), nc_var, ldims );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
2096 [ # # ][ # # ]: 0 : if( !nc_var ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting death status variable" ); }
[ # # ][ # # ]
[ # # ][ # # ]
2097 : 0 : lstart[0] = time_step - 1;
2098 : 0 : lstart[1] = 0;
2099 : 0 : lcount[0] = 1;
2100 [ # # ]: 0 : lcount[1] = i->numElements;
2101 [ # # ][ # # ]: 0 : status = nc_get_vara_double( ncFile, nc_var, lstart, lcount, &death_status[0] );
2102 [ # # ]: 0 : if( NC_NOERR != status )
2103 : : {
2104 [ # # ]: 0 : delete[] exo_conn;
2105 [ # # ][ # # ]: 0 : MB_SET_ERR( MB_FAILURE, "ReadNCDF: Problem setting time step for death_status" );
[ # # ][ # # ]
[ # # ]
2106 : : }
2107 : :
2108 : : // Look for dead elements. If there is too many dead elements and this starts
2109 : : // to take too long, I should put the elems in a kd-tree for more efficient
2110 : : // searching. Alternatively I could get the exo connectivity and match nodes.
2111 : 0 : int dead_elem_counter = 0, missing_elem_counter = 0;
2112 [ # # ][ # # ]: 0 : for( int j = 0; j < i->numElements; ++j )
2113 : : {
2114 [ # # ][ # # ]: 0 : if( 1 != death_status[j] )
2115 : : {
2116 [ # # ][ # # ]: 0 : Range cub_elem, cub_nodes;
[ # # ]
2117 : : if( match_elems_by_connectivity )
2118 : : {
2119 : : // Get exodus nodes for the element
2120 [ # # ]: 0 : std::vector< int > elem_conn( nodes_per_element );
2121 [ # # ]: 0 : for( unsigned int k = 0; k < nodes_per_element; ++k )
2122 : : {
2123 : : // elem_conn[k] = exo_conn[j][k];
2124 [ # # ]: 0 : elem_conn[k] = exo_conn[j * nodes_per_element + k];
2125 : : }
2126 : : // Get the ids of the nodes (assume we are matching by id)
2127 : : // Remember that the exodus array locations start with 1 (not 0).
2128 [ # # ][ # # ]: 0 : std::vector< int > elem_conn_node_ids( nodes_per_element );
2129 [ # # ]: 0 : for( unsigned int k = 0; k < nodes_per_element; ++k )
2130 : : {
2131 [ # # ][ # # ]: 0 : elem_conn_node_ids[k] = ptr[elem_conn[k] - 1];
[ # # ]
2132 : : }
2133 : : // Get the cub_file_set nodes by id
2134 : : // The map is a log search and takes almost no time.
2135 : : // MOAB's linear tag search takes 5-10 minutes.
2136 [ # # ]: 0 : for( unsigned int k = 0; k < nodes_per_element; ++k )
2137 : : {
2138 [ # # ]: 0 : std::map< int, EntityHandle >::iterator k_iter;
2139 [ # # ][ # # ]: 0 : k_iter = matched_cub_vert_id_map.find( elem_conn_node_ids[k] );
2140 : :
2141 [ # # ][ # # ]: 0 : if( k_iter == matched_cub_vert_id_map.end() )
2142 : : {
2143 [ # # ][ # # ]: 0 : std::cout << "ReadNCDF: Found no cub node with id=" << elem_conn_node_ids[k]
[ # # ]
2144 [ # # ][ # # ]: 0 : << ", but expected to find only 1." << std::endl;
2145 : 0 : break;
2146 : : }
2147 [ # # ][ # # ]: 0 : cub_nodes.insert( k_iter->second );
2148 : : }
2149 : :
2150 [ # # ][ # # ]: 0 : if( nodes_per_element != cub_nodes.size() )
2151 : : {
2152 [ # # ][ # # ]: 0 : std::cout << "ReadNCDF: nodes_per_elemenet != cub_nodes.size()" << std::endl;
2153 [ # # ]: 0 : delete[] exo_conn;
2154 : 0 : return MB_INVALID_SIZE;
2155 : : }
2156 : :
2157 : : // Get the cub_file_set element with the same nodes
2158 [ # # ]: 0 : int to_dim = CN::Dimension( mb_type );
2159 [ # # ]: 0 : rval = mdbImpl->get_adjacencies( cub_nodes, to_dim, false, cub_elem );
2160 [ # # ]: 0 : if( MB_SUCCESS != rval )
2161 : : {
2162 [ # # ][ # # ]: 0 : std::cout << "ReadNCDF: could not get dead cub element" << std::endl;
2163 [ # # ]: 0 : delete[] exo_conn;
2164 [ # # ]: 0 : return rval;
2165 : 0 : }
2166 : :
2167 : : // Pronto/Presto renumbers elements, so matching cub and exo elements by
2168 : : // id is not possible at this time.
2169 : : }
2170 : : else
2171 : : {
2172 : : // Get dead element's id
2173 : : int elem_id = elem_ids[first_elem_id_in_block + j];
2174 : : void* id[] = { &elem_id };
2175 : : // Get the element by id
2176 : : rval = mdbImpl->get_entities_by_type_and_tag( cub_file_set, mb_type, &mGlobalIdTag, id, 1, cub_elem,
2177 : : Interface::INTERSECT );
2178 : : if( MB_SUCCESS != rval )
2179 : : {
2180 : : delete[] exo_conn;
2181 : : return rval;
2182 : : }
2183 : : }
2184 : :
2185 [ # # ][ # # ]: 0 : if( 1 == cub_elem.size() )
2186 : : {
2187 : : // Delete the dead element from the cub file. It will be removed from sets
2188 : : // ONLY if they are tracking meshsets.
2189 [ # # ]: 0 : rval = mdbImpl->remove_entities( cub_file_set, cub_elem );
2190 [ # # ]: 0 : if( MB_SUCCESS != rval )
2191 : : {
2192 [ # # ]: 0 : delete[] exo_conn;
2193 : 0 : return rval;
2194 : : }
2195 [ # # ]: 0 : rval = mdbImpl->delete_entities( cub_elem );
2196 [ # # ]: 0 : if( MB_SUCCESS != rval )
2197 : : {
2198 [ # # ]: 0 : delete[] exo_conn;
2199 : 0 : return rval;
2200 : : }
2201 : : }
2202 : : else
2203 : : {
2204 [ # # ][ # # ]: 0 : std::cout << "ReadNCDF: Should have found 1 element with type=" << mb_type
2205 [ # # ][ # # ]: 0 : << " in cub_file_set, but instead found " << cub_elem.size() << std::endl;
[ # # ][ # # ]
2206 [ # # ]: 0 : rval = mdbImpl->list_entities( cub_nodes );
2207 : 0 : ++missing_elem_counter;
2208 [ # # ]: 0 : delete[] exo_conn;
2209 : 0 : return MB_FAILURE;
2210 : : }
2211 [ # # ]: 0 : ++dead_elem_counter;
2212 : : }
2213 : : }
2214 : : // Print some statistics
2215 [ # # ][ # # ]: 0 : temp_ss.str( "" );
2216 [ # # ][ # # ]: 0 : temp_ss << i->blockId;
2217 : 0 : total_dead_elems += dead_elem_counter;
2218 [ # # ]: 0 : total_elems += i->numElements;
2219 [ # # ][ # # ]: 0 : std::cout << " Block " << temp_ss.str() << " has " << dead_elem_counter << "/" << i->numElements
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
2220 [ # # ][ # # ]: 0 : << " dead elements." << std::endl;
2221 [ # # ]: 0 : if( 0 != missing_elem_counter )
2222 : : {
2223 [ # # ][ # # ]: 0 : std::cout << " " << missing_elem_counter
2224 [ # # ][ # # ]: 0 : << " dead elements in this block were not found in the cub_file_set." << std::endl;
2225 : : }
2226 : :
2227 : : // Advance the pointers into element ids and block_count. Memory cleanup.
2228 [ # # ]: 0 : first_elem_id_in_block += i->numElements;
2229 : 0 : ++block_count;
2230 [ # # ][ # # ]: 0 : delete[] exo_conn;
2231 : 0 : }
2232 : :
2233 [ # # ][ # # ]: 0 : std::cout << " Total: " << total_dead_elems << "/" << total_elems << " dead elements." << std::endl;
[ # # ][ # # ]
[ # # ][ # # ]
2234 : :
2235 : : // Close the file
2236 [ # # ]: 0 : fail = nc_close( ncFile );
2237 [ # # ][ # # ]: 0 : if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "Trouble closing file" ); }
[ # # ][ # # ]
[ # # ][ # # ]
2238 : :
2239 : 0 : ncFile = 0;
2240 : 0 : return MB_SUCCESS;
2241 : : }
2242 : :
2243 : 0 : void ReadNCDF::tokenize( const std::string& str, std::vector< std::string >& tokens, const char* delimiters )
2244 : : {
2245 : 0 : std::string::size_type last = str.find_first_not_of( delimiters, 0 );
2246 : 0 : std::string::size_type pos = str.find_first_of( delimiters, last );
2247 [ # # ][ # # ]: 0 : while( std::string::npos != pos && std::string::npos != last )
2248 : : {
2249 [ # # ]: 0 : tokens.push_back( str.substr( last, pos - last ) );
2250 : 0 : last = str.find_first_not_of( delimiters, pos );
2251 : 0 : pos = str.find_first_of( delimiters, last );
2252 [ # # ]: 0 : if( std::string::npos == pos ) pos = str.size();
2253 : : }
2254 : 0 : }
2255 : :
2256 [ + - ][ + - ]: 228 : } // namespace moab
|
