moab
|
00001 00017 #ifdef WIN32 00018 #ifdef _DEBUG 00019 // turn off warnings that say they debugging identifier has been truncated 00020 // this warning comes up when using some STL containers 00021 #pragma warning(disable : 4786) 00022 #endif 00023 #endif 00024 00025 #ifndef NETCDF_FILE 00026 # error Attempt to compile WriteNCDF with NetCDF support disabled 00027 #endif 00028 00029 #include "WriteNCDF.hpp" 00030 00031 #include "netcdf.h" 00032 #include <utility> 00033 #include <algorithm> 00034 #include <time.h> 00035 #include <string> 00036 #include <vector> 00037 #include <stdio.h> 00038 #include <string.h> 00039 #include <assert.h> 00040 00041 #include "moab/Interface.hpp" 00042 #include "moab/Range.hpp" 00043 #include "moab/CN.hpp" 00044 #include "MBTagConventions.hpp" 00045 #include "Internals.hpp" 00046 #include "ExoIIUtil.hpp" 00047 #include "moab/WriteUtilIface.hpp" 00048 #include "exodus_order.h" 00049 00050 namespace moab { 00051 00052 const int TIME_STR_LEN = 11; 00053 00054 #define INS_ID(stringvar, prefix, id) \ 00055 sprintf(stringvar, prefix, id) 00056 00057 #define GET_DIM(ncdim, name, val)\ 00058 { \ 00059 int gdfail = nc_inq_dimid(ncFile, name, &ncdim); \ 00060 if (NC_NOERR == gdfail) { \ 00061 size_t tmp_val; \ 00062 gdfail = nc_inq_dimlen(ncFile, ncdim, &tmp_val); \ 00063 if (NC_NOERR != gdfail) { \ 00064 readMeshIface->report_error("ReadNCDF:: couldn't get dimension length."); \ 00065 return MB_FAILURE; \ 00066 } \ 00067 else val = tmp_val; \ 00068 } else val = 0;} 00069 00070 #define GET_DIMB(ncdim, name, varname, id, val) \ 00071 INS_ID(name, varname, id); \ 00072 GET_DIM(ncdim, name, val); 00073 00074 #define GET_VAR(name, id, dims) \ 00075 { \ 00076 id = -1;\ 00077 int gvfail = nc_inq_varid(ncFile, name, &id); \ 00078 if (NC_NOERR == gvfail) { \ 00079 int ndims;\ 00080 gvfail = nc_inq_varndims(ncFile, id, &ndims);\ 00081 if (NC_NOERR == gvfail) {\ 00082 dims.resize(ndims); \ 00083 gvfail = nc_inq_vardimid(ncFile, id, &dims[0]);}}} 00084 00085 WriterIface* WriteNCDF::factory( Interface* iface ) 00086 { return new WriteNCDF( iface ); } 00087 00088 WriteNCDF::WriteNCDF(Interface *impl) 00089 : mdbImpl(impl), ncFile(0), mCurrentMeshHandle(0) 00090 { 00091 assert(impl != NULL); 00092 00093 impl->query_interface( mWriteIface ); 00094 00095 // initialize in case tag_get_handle fails below 00097 int zero = 0, negone = -1; 00098 impl->tag_get_handle(MATERIAL_SET_TAG_NAME, 1, MB_TYPE_INTEGER, 00099 mMaterialSetTag, MB_TAG_SPARSE|MB_TAG_CREAT, &negone); 00100 00101 impl->tag_get_handle(DIRICHLET_SET_TAG_NAME, 1, MB_TYPE_INTEGER, 00102 mDirichletSetTag, MB_TAG_SPARSE|MB_TAG_CREAT, &negone); 00103 00104 impl->tag_get_handle(NEUMANN_SET_TAG_NAME, 1, MB_TYPE_INTEGER, 00105 mNeumannSetTag, MB_TAG_SPARSE|MB_TAG_CREAT, &negone); 00106 00107 impl->tag_get_handle(GLOBAL_ID_TAG_NAME, 1, MB_TYPE_INTEGER, 00108 mGlobalIdTag, MB_TAG_SPARSE|MB_TAG_CREAT, &zero); 00109 00110 int dum_val_array[] = {-1, -1, -1, -1}; 00111 impl->tag_get_handle(HAS_MID_NODES_TAG_NAME, 4, MB_TYPE_INTEGER, 00112 mHasMidNodesTag, MB_TAG_SPARSE|MB_TAG_CREAT, dum_val_array); 00113 00114 impl->tag_get_handle( "distFactor", 0, MB_TYPE_DOUBLE, mDistFactorTag, 00115 MB_TAG_SPARSE|MB_TAG_VARLEN|MB_TAG_CREAT ); 00116 00117 impl->tag_get_handle( "qaRecord", 0, MB_TYPE_OPAQUE, mQaRecordTag, 00118 MB_TAG_SPARSE|MB_TAG_VARLEN|MB_TAG_CREAT ); 00119 00120 impl->tag_get_handle("WriteNCDF element mark", 1, MB_TYPE_BIT, mEntityMark, MB_TAG_CREAT); 00121 00122 } 00123 00124 WriteNCDF::~WriteNCDF() 00125 { 00126 mdbImpl->release_interface(mWriteIface); 00127 00128 mdbImpl->tag_delete(mEntityMark); 00129 00130 if (ncFile) ncFile = 0; 00131 } 00132 00133 void WriteNCDF::reset_block(std::vector<MaterialSetData> &block_info) 00134 { 00135 std::vector<MaterialSetData>::iterator iter; 00136 00137 for (iter = block_info.begin(); iter != block_info.end(); iter++) 00138 { 00139 iter->elements.clear(); 00140 } 00141 } 00142 00143 void WriteNCDF::time_and_date(char* time_string, char* date_string) 00144 { 00145 struct tm* local_time; 00146 time_t calendar_time; 00147 00148 calendar_time = time(NULL); 00149 local_time = localtime(&calendar_time); 00150 00151 assert(NULL != time_string && NULL != date_string); 00152 00153 strftime(time_string, TIME_STR_LEN, "%H:%M:%S", local_time); 00154 strftime(date_string, TIME_STR_LEN, "%m/%d/%Y", local_time); 00155 00156 // terminate with NULL character 00157 time_string[10] = (char)NULL; 00158 date_string[10] = (char)NULL; 00159 } 00160 00161 ErrorCode WriteNCDF::write_file(const char *exodus_file_name, 00162 const bool overwrite, 00163 const FileOptions&, 00164 const EntityHandle *ent_handles, 00165 const int num_sets, 00166 const std::vector<std::string> &qa_records, 00167 const Tag*, 00168 int, 00169 int user_dimension) 00170 { 00171 assert(0 != mMaterialSetTag && 00172 0 != mNeumannSetTag && 00173 0 != mDirichletSetTag); 00174 00175 if (user_dimension == 0) 00176 mdbImpl->get_dimension( user_dimension ); 00177 00178 00179 std::vector<EntityHandle> blocks, nodesets, sidesets, entities; 00180 00181 // separate into blocks, nodesets, sidesets 00182 00183 if (num_sets == 0) { 00184 // default to all defined block, nodeset and sideset-type sets 00185 Range this_range; 00186 mdbImpl->get_entities_by_type_and_tag(0, MBENTITYSET, &mMaterialSetTag, NULL, 1, this_range); 00187 std::copy(this_range.begin(), this_range.end(), std::back_inserter(blocks)); 00188 this_range.clear(); 00189 mdbImpl->get_entities_by_type_and_tag(0, MBENTITYSET, &mDirichletSetTag, NULL, 1, this_range); 00190 std::copy(this_range.begin(), this_range.end(), std::back_inserter(nodesets)); 00191 this_range.clear(); 00192 mdbImpl->get_entities_by_type_and_tag(0, MBENTITYSET, &mNeumannSetTag, NULL, 1, this_range); 00193 std::copy(this_range.begin(), this_range.end(), std::back_inserter(sidesets)); 00194 00195 // If there is nothing to write, write everything as one block. 00196 if (blocks.empty() && nodesets.empty() && sidesets.empty()) 00197 { 00198 this_range.clear(); 00199 for (int d = user_dimension; d > 0 && this_range.empty(); --d) 00200 mdbImpl->get_entities_by_dimension( 0, d, this_range, false ); 00201 if (this_range.empty()) 00202 return MB_FILE_WRITE_ERROR; 00203 00204 EntityHandle block_handle; 00205 int block_id = 1; 00206 mdbImpl->create_meshset( MESHSET_SET, block_handle ); 00207 mdbImpl->tag_set_data( mMaterialSetTag, &block_handle, 1, &block_id ); 00208 mdbImpl->add_entities( block_handle, this_range ); 00209 blocks.push_back( block_handle ); 00210 } 00211 } 00212 else { 00213 int dummy; 00214 for (const EntityHandle *iter = ent_handles; iter < ent_handles+num_sets; iter++) 00215 { 00216 if (MB_SUCCESS == mdbImpl->tag_get_data(mMaterialSetTag, &(*iter), 1, &dummy) && 00217 -1 != dummy) 00218 blocks.push_back(*iter); 00219 else if (MB_SUCCESS == mdbImpl->tag_get_data(mDirichletSetTag, &(*iter), 1, &dummy) && 00220 -1 != dummy) 00221 nodesets.push_back(*iter); 00222 else if (MB_SUCCESS == mdbImpl->tag_get_data(mNeumannSetTag, &(*iter), 1, &dummy) && 00223 -1 != dummy) 00224 sidesets.push_back(*iter); 00225 } 00226 } 00227 00228 00229 // if there is nothing to write just return. 00230 if (blocks.empty() && nodesets.empty() && sidesets.empty()) 00231 return MB_FILE_WRITE_ERROR; 00232 00233 // try to get mesh information 00234 ExodusMeshInfo mesh_info; 00235 00236 std::vector<MaterialSetData> block_info; 00237 std::vector<NeumannSetData> sideset_info; 00238 std::vector<DirichletSetData> nodeset_info; 00239 00240 mesh_info.num_dim = user_dimension; 00241 00242 if (qa_records.empty()) { 00243 // qa records are empty - initialize some MB-standard ones 00244 mesh_info.qaRecords.push_back("MB"); 00245 mesh_info.qaRecords.push_back("0.99"); 00246 char string1[80], string2[80]; 00247 time_and_date(string2, string1); 00248 mesh_info.qaRecords.push_back(string2); 00249 mesh_info.qaRecords.push_back(string1); 00250 } 00251 else { 00252 // constrained to multiples of 4 qa records 00253 assert(qa_records.size()%4 == 0); 00254 00255 std::copy(qa_records.begin(), qa_records.end(), 00256 std::back_inserter(mesh_info.qaRecords)); 00257 } 00258 00259 block_info.clear(); 00260 if(gather_mesh_information(mesh_info, block_info, sideset_info, nodeset_info, 00261 blocks, sidesets, nodesets) != MB_SUCCESS) 00262 { 00263 reset_block(block_info); 00264 return MB_FAILURE; 00265 } 00266 00267 00268 // try to open the file after gather mesh info succeeds 00269 int fail = nc_create(exodus_file_name, overwrite ? NC_CLOBBER : NC_NOCLOBBER, &ncFile); 00270 if (NC_NOERR != fail) { 00271 reset_block(block_info); 00272 return MB_FAILURE; 00273 } 00274 00275 if( write_header(mesh_info, block_info, sideset_info, 00276 nodeset_info, exodus_file_name) != MB_SUCCESS) 00277 { 00278 reset_block(block_info); 00279 return MB_FAILURE; 00280 } 00281 00282 if( write_nodes(mesh_info.num_nodes, mesh_info.nodes, mesh_info.num_dim) != MB_SUCCESS ) 00283 { 00284 reset_block(block_info); 00285 return MB_FAILURE; 00286 } 00287 00288 if( write_elementblocks(block_info) ) 00289 { 00290 reset_block(block_info); 00291 return MB_FAILURE; 00292 } 00293 00294 // write the three maps 00295 if( write_global_node_order_map(mesh_info.num_nodes, mesh_info.nodes) != MB_SUCCESS ) 00296 { 00297 reset_block(block_info); 00298 return MB_FAILURE; 00299 } 00300 00301 if(write_global_element_order_map(mesh_info.num_elements) != MB_SUCCESS ) 00302 { 00303 reset_block(block_info); 00304 return MB_FAILURE; 00305 } 00306 00307 if(write_element_order_map(mesh_info.num_elements) != MB_SUCCESS ) 00308 { 00309 reset_block(block_info); 00310 return MB_FAILURE; 00311 } 00312 00313 /* 00314 if(write_elementmap(mesh_info) != MB_SUCCESS) 00315 return MB_FAILURE; 00316 */ 00317 00318 if(write_BCs(sideset_info, nodeset_info) != MB_SUCCESS) 00319 { 00320 reset_block(block_info); 00321 return MB_FAILURE; 00322 } 00323 00324 if( write_qa_records( mesh_info.qaRecords) != MB_SUCCESS ) 00325 return MB_FAILURE; 00326 00327 // copy the qa records into the argument 00328 //mesh_info.qaRecords.swap(qa_records); 00329 // close the file 00330 fail = nc_close(ncFile); 00331 if (NC_NOERR != fail) { 00332 mWriteIface->report_error("Trouble closing file."); 00333 return MB_FAILURE; 00334 } 00335 00336 return MB_SUCCESS; 00337 } 00338 00339 ErrorCode WriteNCDF::gather_mesh_information( 00340 ExodusMeshInfo &mesh_info, 00341 std::vector<MaterialSetData> &block_info, 00342 std::vector<NeumannSetData> &sideset_info, 00343 std::vector<DirichletSetData> &nodeset_info, 00344 std::vector<EntityHandle> &blocks, 00345 std::vector<EntityHandle> &sidesets, 00346 std::vector<EntityHandle> &nodesets) 00347 { 00348 ErrorCode rval; 00349 std::vector<EntityHandle>::iterator vector_iter, end_vector_iter; 00350 00351 mesh_info.num_nodes = 0; 00352 mesh_info.num_elements = 0; 00353 mesh_info.num_elementblocks = 0; 00354 00355 int id = 0; 00356 00357 vector_iter= blocks.begin(); 00358 end_vector_iter = blocks.end(); 00359 00360 std::vector<EntityHandle> parent_meshsets; 00361 00362 // clean out the bits for the element mark 00363 rval = mdbImpl->tag_delete(mEntityMark); 00364 if (MB_SUCCESS != rval) 00365 return rval; 00366 rval = mdbImpl->tag_get_handle("WriteNCDF element mark", 1, MB_TYPE_BIT, mEntityMark, MB_TAG_CREAT); 00367 if (MB_SUCCESS != rval) 00368 return rval; 00369 00370 int highest_dimension_of_element_blocks = 0; 00371 00372 for(vector_iter = blocks.begin(); vector_iter != blocks.end(); vector_iter++) 00373 { 00374 00375 MaterialSetData block_data; 00376 00377 //for the purpose of qa records, get the parents of these blocks 00378 if( mdbImpl->get_parent_meshsets( *vector_iter, parent_meshsets ) != MB_SUCCESS ) 00379 return MB_FAILURE; 00380 00381 // get all Entity Handles in the mesh set 00382 Range dummy_range; 00383 rval = mdbImpl->get_entities_by_handle(*vector_iter, dummy_range, true ); 00384 if (MB_SUCCESS != rval) 00385 return rval; 00386 00387 // skip empty blocks 00388 if (dummy_range.empty()) 00389 continue; 00390 00391 // get the block's id 00392 if(mdbImpl->tag_get_data(mMaterialSetTag, &(*vector_iter), 1, &id) != MB_SUCCESS ) { 00393 mWriteIface->report_error("Couldn't get block id from a tag for an element block."); 00394 return MB_FAILURE; 00395 } 00396 00397 block_data.id = id; 00398 block_data.number_attributes = 0; 00399 00400 // wait a minute, we are doing some filtering here that doesn't make sense at this level CJS 00401 00402 // find the dimension of the last entity in this range 00403 int this_dim = CN::Dimension(TYPE_FROM_HANDLE(dummy_range.back())); 00404 if (this_dim > 3) { 00405 mWriteIface->report_error("Block %d contains entity sets.", id); 00406 return MB_TYPE_OUT_OF_RANGE; 00407 } 00408 block_data.elements = dummy_range.subset_by_dimension( this_dim ); 00409 00410 // end of -- wait a minute, we are doing some filtering here that doesn't make sense at this level CJS 00411 00412 00413 // get the entity type for this block, verifying that it's the same for all elements 00414 EntityType entity_type = TYPE_FROM_HANDLE(block_data.elements.front()); 00415 if (!block_data.elements.all_of_type(entity_type)) { 00416 mWriteIface->report_error("Entities in block %i not of common type", id); 00417 return MB_FAILURE; 00418 } 00419 00420 int dimension = -1; 00421 if(entity_type == MBQUAD || entity_type == MBTRI) 00422 dimension = 3; // output shells by default 00423 else if(entity_type == MBEDGE) 00424 dimension = 2; // SHOULD THIS BE 1?? -- J.Kraftcheck, August, 2011 00425 else 00426 dimension = CN::Dimension(entity_type); 00427 00428 if( dimension > highest_dimension_of_element_blocks ) 00429 highest_dimension_of_element_blocks = dimension; 00430 00431 std::vector<EntityHandle> tmp_conn; 00432 rval = mdbImpl->get_connectivity(&(block_data.elements.front()), 1, tmp_conn); 00433 if (MB_SUCCESS != rval) 00434 return rval; 00435 block_data.element_type = ExoIIUtil::get_element_type_from_num_verts(tmp_conn.size(), entity_type, dimension); 00436 00437 if (block_data.element_type == EXOII_MAX_ELEM_TYPE) { 00438 mWriteIface->report_error("Element type in block %i didn't get set correctly", id); 00439 return MB_FAILURE; 00440 } 00441 00442 block_data.number_nodes_per_element = ExoIIUtil::VerticesPerElement[block_data.element_type]; 00443 00444 // number of nodes for this block 00445 block_data.number_elements = block_data.elements.size(); 00446 00447 // total number of elements 00448 mesh_info.num_elements += block_data.number_elements; 00449 00450 // get the nodes for the elements 00451 rval = mWriteIface->gather_nodes_from_elements(block_data.elements, mEntityMark, mesh_info.nodes); 00452 if (MB_SUCCESS != rval) 00453 return rval; 00454 00455 if(!sidesets.empty()) 00456 { 00457 // if there are sidesets, keep track of which elements are being written out 00458 for(Range::iterator iter = block_data.elements.begin(); 00459 iter != block_data.elements.end(); ++iter) 00460 { 00461 unsigned char bit = 0x1; 00462 rval = mdbImpl->tag_set_data(mEntityMark, &(*iter), 1, &bit); 00463 if (MB_SUCCESS != rval) 00464 return rval; 00465 } 00466 } 00467 00468 block_info.push_back( block_data ); 00469 00470 std::vector<char*> *qa_rec = NULL; 00471 00472 if (mdbImpl->tag_get_data(mQaRecordTag, &(*vector_iter), 1, &qa_rec) == MB_SUCCESS && 00473 NULL != qa_rec) 00474 { 00475 // there are qa records on this block - copy over to mesh qa records 00476 00477 // constrained to multiples of 4 qa records 00478 assert(qa_rec->size()%4 == 0); 00479 00480 for(unsigned int k=0; k<qa_rec->size(); k++) 00481 mesh_info.qaRecords.push_back( (*qa_rec)[k] ); 00482 } 00483 } 00484 00485 00486 mesh_info.num_elementblocks = block_info.size(); 00487 00488 //if user hasn't entered dimension, we figure it out 00489 if( mesh_info.num_dim == 0 ) 00490 { 00491 //never want 1 or zero dimensions 00492 if( highest_dimension_of_element_blocks < 2 ) 00493 mesh_info.num_dim = 3; 00494 else 00495 mesh_info.num_dim = highest_dimension_of_element_blocks; 00496 } 00497 00498 00499 Range::iterator range_iter, end_range_iter; 00500 range_iter = mesh_info.nodes.begin(); 00501 end_range_iter = mesh_info.nodes.end(); 00502 00503 mesh_info.num_nodes = mesh_info.nodes.size(); 00504 00505 //------nodesets-------- 00506 00507 vector_iter= nodesets.begin(); 00508 end_vector_iter = nodesets.end(); 00509 00510 for(; vector_iter != end_vector_iter; vector_iter++) 00511 { 00512 00513 DirichletSetData nodeset_data; 00514 nodeset_data.id = 0; 00515 nodeset_data.number_nodes = 0; 00516 00517 // get the nodeset's id 00518 if(mdbImpl->tag_get_data(mDirichletSetTag,&(*vector_iter), 1,&id) != MB_SUCCESS) { 00519 mWriteIface->report_error("Couldn't get id tag for nodeset %i", id); 00520 return MB_FAILURE; 00521 } 00522 00523 nodeset_data.id = id; 00524 00525 std::vector<EntityHandle> node_vector; 00526 //get the nodes of the nodeset that are in mesh_info.nodes 00527 if( mdbImpl->get_entities_by_handle(*vector_iter, node_vector, true) != MB_SUCCESS ) { 00528 mWriteIface->report_error("Couldn't get nodes in nodeset %i", id); 00529 return MB_FAILURE; 00530 } 00531 00532 00533 00534 //get the tag for distribution factors 00535 const double *dist_factor_vector; 00536 int dist_factor_size; 00537 const void* ptr = 0; 00538 00539 int has_dist_factors = 0; 00540 if(mdbImpl->tag_get_by_ptr(mDistFactorTag,&(*vector_iter), 1, &ptr, &dist_factor_size) == MB_SUCCESS && 00541 dist_factor_size) 00542 has_dist_factors = 1; 00543 dist_factor_size /= sizeof(double); 00544 dist_factor_vector = reinterpret_cast<const double*>(ptr); 00545 std::vector<EntityHandle>::iterator iter, end_iter; 00546 iter = node_vector.begin(); 00547 end_iter= node_vector.end(); 00548 00549 int j=0; 00550 unsigned char node_marked = 0; 00551 ErrorCode result; 00552 for(; iter != end_iter; iter++) 00553 { 00554 if (TYPE_FROM_HANDLE(*iter) != MBVERTEX) continue; 00555 result = mdbImpl->tag_get_data(mEntityMark, &(*iter), 1, &node_marked); 00556 if (MB_SUCCESS != result) { 00557 mWriteIface->report_error("Couldn't get mark data."); 00558 return result; 00559 } 00560 00561 if(node_marked == 0x1) 00562 { 00563 nodeset_data.nodes.push_back( *iter ); 00564 if( has_dist_factors != 0) 00565 nodeset_data.node_dist_factors.push_back( dist_factor_vector[j] ); 00566 else 00567 nodeset_data.node_dist_factors.push_back( 1.0 ); 00568 } 00569 j++; 00570 } 00571 00572 nodeset_data.number_nodes = nodeset_data.nodes.size(); 00573 nodeset_info.push_back( nodeset_data ); 00574 00575 } 00576 00577 //------sidesets-------- 00578 vector_iter= sidesets.begin(); 00579 end_vector_iter = sidesets.end(); 00580 00581 for(; vector_iter != end_vector_iter; vector_iter++) 00582 { 00583 00584 NeumannSetData sideset_data; 00585 00586 // get the sideset's id 00587 if(mdbImpl->tag_get_data(mNeumannSetTag,&(*vector_iter), 1,&id) != MB_SUCCESS) 00588 return MB_FAILURE; 00589 00590 sideset_data.id = id; 00591 sideset_data.mesh_set_handle = *vector_iter; 00592 00593 //get the sides in two lists, one forward the other reverse; starts with forward sense 00594 // by convention 00595 Range forward_elems, reverse_elems; 00596 if(get_sideset_elems(*vector_iter, 0, forward_elems, reverse_elems) == MB_FAILURE) 00597 return MB_FAILURE; 00598 00599 ErrorCode result = get_valid_sides(forward_elems, mesh_info, 1, sideset_data); 00600 if (MB_SUCCESS != result) { 00601 mWriteIface->report_error("Couldn't get valid sides data."); 00602 return result; 00603 } 00604 result = get_valid_sides(reverse_elems, mesh_info, -1, sideset_data); 00605 if (MB_SUCCESS != result) { 00606 mWriteIface->report_error("Couldn't get valid sides data."); 00607 return result; 00608 } 00609 00610 sideset_data.number_elements = sideset_data.elements.size(); 00611 sideset_info.push_back( sideset_data ); 00612 00613 } 00614 00615 return MB_SUCCESS; 00616 00617 } 00618 00619 ErrorCode WriteNCDF::get_valid_sides(Range &elems, ExodusMeshInfo& /*mesh_info*/, 00620 const int sense, 00621 NeumannSetData &sideset_data) 00622 { 00623 // this is where we see if underlying element of side set element is included in output 00624 00625 //get the sideset-based info for distribution factors 00626 const double *dist_factor_vector = 0; 00627 int dist_factor_size = 0; 00628 00629 // initialize dist_fac_iter to get rid of compiler warning 00630 const double* dist_fac_iter = 0; 00631 const void* ptr = 0; 00632 bool has_dist_factors = false; 00633 if(mdbImpl->tag_get_by_ptr(mDistFactorTag, 00634 &(sideset_data.mesh_set_handle), 1, &ptr, &dist_factor_size) == MB_SUCCESS && 00635 dist_factor_size) 00636 { 00637 has_dist_factors = true; 00638 dist_factor_vector = reinterpret_cast<const double*>(ptr); 00639 dist_fac_iter = dist_factor_vector; 00640 dist_factor_size /= sizeof(double); 00641 } 00642 00643 unsigned char element_marked = 0; 00644 ErrorCode result; 00645 for(Range::iterator iter = elems.begin(); iter != elems.end(); iter++) 00646 { 00647 // should insert here if "side" is a quad/tri on a quad/tri mesh 00648 result = mdbImpl->tag_get_data(mEntityMark, &(*iter), 1, &element_marked); 00649 if (MB_SUCCESS != result) { 00650 mWriteIface->report_error("Couldn't get mark data."); 00651 return result; 00652 } 00653 00654 if(element_marked == 0x1) 00655 { 00656 sideset_data.elements.push_back( *iter ); 00657 00658 // TJT TODO: the sense should really be # edges + 1or2 00659 sideset_data.side_numbers.push_back((sense == 1 ? 1 : 2)); 00660 } 00661 else //then "side" is probably a quad/tri on a hex/tet mesh 00662 { 00663 std::vector<EntityHandle> parents; 00664 int dimension = CN::Dimension( TYPE_FROM_HANDLE(*iter)); 00665 00666 //get the adjacent parent element of "side" 00667 if( mdbImpl->get_adjacencies( &(*iter), 1, dimension+1, false, parents) != MB_SUCCESS ) { 00668 mWriteIface->report_error("Couldn't get adjacencies for sideset."); 00669 return MB_FAILURE; 00670 } 00671 00672 if(!parents.empty()) 00673 { 00674 //make sure the adjacent parent element will be output 00675 for(unsigned int k=0; k<parents.size(); k++) 00676 { 00677 result = mdbImpl->tag_get_data(mEntityMark, &(parents[k]), 1, &element_marked); 00678 if (MB_SUCCESS != result) { 00679 mWriteIface->report_error("Couldn't get mark data."); 00680 return result; 00681 } 00682 00683 int side_no, this_sense, this_offset; 00684 if(element_marked == 0x1 && 00685 mdbImpl->side_number(parents[k], *iter, side_no, 00686 this_sense, this_offset) == MB_SUCCESS && 00687 this_sense == sense) { 00688 sideset_data.elements.push_back(parents[k]); 00689 sideset_data.side_numbers.push_back(side_no+1); 00690 break; 00691 } 00692 } 00693 } 00694 else 00695 { 00696 mWriteIface->report_error("No parent element exists for element in sideset %i", sideset_data.id); 00697 return MB_FAILURE; 00698 } 00699 } 00700 00701 if( sideset_data.elements.size() != 0 ) 00702 { 00703 // distribution factors 00704 int num_nodes = CN::VerticesPerEntity(TYPE_FROM_HANDLE(*iter)); 00705 if( has_dist_factors ) 00706 { 00707 std::copy(dist_fac_iter, dist_fac_iter + num_nodes, 00708 std::back_inserter(sideset_data.ss_dist_factors) ); 00709 dist_fac_iter += num_nodes; 00710 } 00711 else 00712 { 00713 for(int j=0; j<num_nodes; j++) 00714 sideset_data.ss_dist_factors.push_back( 1.0 ); 00715 } 00716 } 00717 } 00718 00719 return MB_SUCCESS; 00720 } 00721 00722 ErrorCode WriteNCDF::write_qa_records(std::vector<std::string> &qa_record_list) 00723 { 00724 int i = 0; 00725 00726 for(std::vector<std::string>::iterator string_it = qa_record_list.begin(); 00727 string_it != qa_record_list.end(); ) 00728 { 00729 for (int j = 0; j < 4; j++) 00730 { 00731 write_qa_string((*string_it++).c_str(), i, j); 00732 } 00733 i++; 00734 } 00735 00736 return MB_SUCCESS; 00737 } 00738 00739 ErrorCode WriteNCDF::write_qa_string(const char *string, 00740 int record_number, 00741 int record_position) 00742 { 00743 // get the variable id in the exodus file 00744 00745 std::vector<int> dims; 00746 int temp_var = -1; 00747 GET_VAR("qa_records", temp_var, dims); 00748 if (-1 == temp_var) { 00749 mWriteIface->report_error("WriteNCDF:: Problem getting qa record variable."); 00750 return MB_FAILURE; 00751 } 00752 size_t count[3], start[3]; 00753 00754 // write out the record 00755 start[0] = record_number; 00756 start[1] = record_position; 00757 start[2] = 0; 00758 00759 count[0] = 1; 00760 count[1] = 1; 00761 count[2] = (long)strlen(string) +1; 00762 int fail = nc_put_vara_text(ncFile, temp_var, start, count, string); 00763 if (NC_NOERR != fail) { 00764 mWriteIface->report_error("Failed to position qa string variable."); 00765 return MB_FAILURE; 00766 } 00767 00768 return MB_SUCCESS; 00769 } 00770 00771 00772 ErrorCode WriteNCDF::write_nodes(int num_nodes, Range& nodes, int dimension) 00773 { 00774 // write coordinates names 00775 int nc_var = -1; 00776 std::vector<int> dims; 00777 GET_VAR("coor_names", nc_var, dims); 00778 if (-1 == nc_var) { 00779 mWriteIface->report_error("Trouble getting coordinate name variable."); 00780 return MB_FAILURE; 00781 } 00782 00783 size_t start[2] = {0, 0}, count[2] = {1, ExoIIInterface::MAX_STR_LENGTH}; 00784 char dum_str[ExoIIInterface::MAX_STR_LENGTH]; 00785 strcpy(dum_str, "x"); 00786 int fail = nc_put_vara_text(ncFile, nc_var, start, count, dum_str); 00787 if (NC_NOERR != fail) { 00788 mWriteIface->report_error("Trouble adding x coordinate name; netcdf message:"); 00789 const char *err = nc_strerror(fail); 00790 mWriteIface->report_error("%s", err); 00791 return MB_FAILURE; 00792 } 00793 00794 start[0] = 1; 00795 strcpy(dum_str, "y"); 00796 fail = nc_put_vara_text(ncFile, nc_var, start, count, dum_str); 00797 if (NC_NOERR != fail) { 00798 mWriteIface->report_error("Trouble adding y coordinate name."); 00799 const char *err = nc_strerror(fail); 00800 mWriteIface->report_error("%s", err); 00801 return MB_FAILURE; 00802 } 00803 00804 start[0] = 2; 00805 strcpy(dum_str, "z"); 00806 fail = nc_put_vara_text(ncFile, nc_var, start, count, dum_str); 00807 if (NC_NOERR != fail) { 00808 mWriteIface->report_error("Trouble adding z coordinate name."); 00809 const char *err = nc_strerror(fail); 00810 mWriteIface->report_error("%s", err); 00811 return MB_FAILURE; 00812 } 00813 00814 //see if should transform coordinates 00815 ErrorCode result; 00816 Tag trans_tag; 00817 result = mdbImpl->tag_get_handle( MESH_TRANSFORM_TAG_NAME, 16, MB_TYPE_DOUBLE, trans_tag); 00818 bool transform_needed = true; 00819 if( result == MB_TAG_NOT_FOUND ) 00820 transform_needed = false; 00821 00822 int num_coords_to_fill = transform_needed ? 3 : dimension; 00823 00824 std::vector<double*> coord_arrays(3); 00825 coord_arrays[0] = new double[num_nodes]; 00826 coord_arrays[1] = new double[num_nodes]; 00827 coord_arrays[2] = NULL; 00828 00829 if( num_coords_to_fill == 3 ) 00830 coord_arrays[2] = new double[num_nodes]; 00831 00832 result = mWriteIface->get_node_coords(dimension, num_nodes, nodes, mGlobalIdTag, 1, coord_arrays); 00833 if(result != MB_SUCCESS) 00834 { 00835 delete [] coord_arrays[0]; 00836 delete [] coord_arrays[1]; 00837 if(coord_arrays[2]) delete [] coord_arrays[2]; 00838 return result; 00839 } 00840 00841 if( transform_needed ) 00842 { 00843 double trans_matrix[16]; 00844 const EntityHandle mesh = 0; 00845 result = mdbImpl->tag_get_data( trans_tag, &mesh, 0, trans_matrix ); 00846 if (MB_SUCCESS != result) { 00847 mWriteIface->report_error("Couldn't get transform data."); 00848 return result; 00849 } 00850 00851 for( int i=0; i<num_nodes; i++) 00852 { 00853 00854 double vec1[3]; 00855 double vec2[3]; 00856 00857 vec2[0] = coord_arrays[0][i]; 00858 vec2[1] = coord_arrays[1][i]; 00859 vec2[2] = coord_arrays[2][i]; 00860 00861 for( int row=0; row<3; row++ ) 00862 { 00863 vec1[row] = 0.0; 00864 for( int col = 0; col<3; col++ ) 00865 { 00866 vec1[row] += ( trans_matrix[ (row*4)+col ] * vec2[col] ); 00867 } 00868 } 00869 00870 coord_arrays[0][i] = vec1[0]; 00871 coord_arrays[1][i] = vec1[1]; 00872 coord_arrays[2][i] = vec1[2]; 00873 00874 } 00875 } 00876 00877 00878 // write the nodes 00879 nc_var = -1; 00880 GET_VAR("coord", nc_var, dims); 00881 if (-1 == nc_var) { 00882 mWriteIface->report_error("Trouble getting coordinate variable."); 00883 return MB_FAILURE; 00884 } 00885 start[0] = 0; 00886 count[1] = num_nodes; 00887 fail = nc_put_vara_double(ncFile, nc_var, start, count, &(coord_arrays[0][0])); 00888 if (NC_NOERR != fail) { 00889 mWriteIface->report_error("Trouble writing x coordinate."); 00890 return MB_FAILURE; 00891 } 00892 00893 start[0] = 1; 00894 fail = nc_put_vara_double(ncFile, nc_var, start, count, &(coord_arrays[1][0])); 00895 if (NC_NOERR != fail) { 00896 mWriteIface->report_error("Trouble writing y coordinate."); 00897 return MB_FAILURE; 00898 } 00899 00900 start[0] = 2; 00901 fail = nc_put_vara_double(ncFile, nc_var, start, count, &(coord_arrays[2][0])); 00902 if (NC_NOERR != fail) { 00903 mWriteIface->report_error("Trouble writing z coordinate."); 00904 return MB_FAILURE; 00905 } 00906 00907 delete [] coord_arrays[0]; 00908 delete [] coord_arrays[1]; 00909 if(coord_arrays[2]) delete [] coord_arrays[2]; 00910 00911 return MB_SUCCESS; 00912 00913 } 00914 00915 00916 ErrorCode WriteNCDF::write_header(ExodusMeshInfo& mesh_info, 00917 std::vector<MaterialSetData> &block_info, 00918 std::vector<NeumannSetData> &sideset_info, 00919 std::vector<DirichletSetData> &nodeset_info, 00920 const char *filename) 00921 { 00922 00923 // get the date and time 00924 char time[TIME_STR_LEN]; 00925 char date[TIME_STR_LEN]; 00926 time_and_date(time,date); 00927 00928 std::string title_string = "MOAB"; 00929 title_string.append( "(" ); 00930 title_string.append( filename ); 00931 title_string.append( "): "); 00932 title_string.append( date ); 00933 title_string.append( ": " ); 00934 title_string.append( "time " ); 00935 00936 if(title_string.length() > ExoIIInterface::MAX_LINE_LENGTH) 00937 title_string.resize( ExoIIInterface::MAX_LINE_LENGTH ); 00938 00939 // initialize the exodus file 00940 00941 int result = initialize_exodus_file(mesh_info, 00942 block_info, sideset_info, 00943 nodeset_info, title_string.c_str()); 00944 00945 if(result == MB_FAILURE) 00946 return MB_FAILURE; 00947 00948 return MB_SUCCESS; 00949 } 00950 00951 ErrorCode WriteNCDF::write_elementblocks(std::vector<MaterialSetData> &block_data ) 00952 { 00953 00954 unsigned int i; 00955 int block_index = 0; // index into block list, may != 1 if there are inactive blocks 00956 int exodus_id = 1; 00957 00958 for(i=0; i< block_data.size(); i++) 00959 { 00960 MaterialSetData& block = block_data[i]; 00961 00962 unsigned int num_nodes_per_elem = block.number_nodes_per_element; 00963 00964 // write out the id for the block 00965 00966 int id = block.id; 00967 int num_values = 1; 00968 00969 if( write_exodus_integer_variable("eb_prop1", &id, block_index, num_values) != MB_SUCCESS ) 00970 { 00971 mWriteIface->report_error("Problem writing element block id %i", id); 00972 } 00973 00974 // write out the block status 00975 00976 int status = 1; 00977 if (!block.number_elements) 00978 { 00979 mWriteIface->report_error("Warning: No elements in block %i", id); 00980 } 00981 00982 if( write_exodus_integer_variable("eb_status", &status, block_index, num_values) != MB_SUCCESS ) 00983 { 00984 mWriteIface->report_error("Problem writing element block status"); 00985 } 00986 00987 // map the connectivity to the new nodes 00988 const unsigned int num_elem = block.number_elements; 00989 const unsigned int num_nodes = num_nodes_per_elem * num_elem; 00990 int* connectivity = new int[num_nodes]; 00991 00992 ErrorCode result = mWriteIface->get_element_connect( 00993 num_elem, num_nodes_per_elem, mGlobalIdTag, block.elements, mGlobalIdTag, exodus_id ,connectivity); 00994 00995 if(result != MB_SUCCESS) { 00996 mWriteIface->report_error("Couldn't get element array to write from."); 00997 delete [] connectivity; 00998 return result; 00999 } 01000 01001 // if necessary, convert from EXODUS to CN node order 01002 const EntityType elem_type = ExoIIUtil::ExoIIElementMBEntity[block.element_type]; 01003 assert( block.elements.all_of_type( elem_type ) ); 01004 const int* reorder = exodus_elem_order_map[elem_type][block.number_nodes_per_element]; 01005 if (reorder) 01006 WriteUtilIface::reorder( reorder, connectivity, 01007 block.number_elements, 01008 block.number_nodes_per_element ); 01009 01010 exodus_id += num_elem; 01011 01012 char wname[80]; 01013 INS_ID(wname, "connect%u", i+1); 01014 std::vector<int> dims; 01015 int nc_var = -1; 01016 GET_VAR(wname, nc_var, dims); 01017 if (-1 == nc_var) { 01018 mWriteIface->report_error("Couldn't get connectivity variable."); 01019 delete [] connectivity; 01020 return MB_FAILURE; 01021 } 01022 01023 size_t start[2] = {0, 0}, count[2] = {num_elem, num_nodes_per_elem}; 01024 int fail = nc_put_vara_int(ncFile, nc_var, start, count, connectivity); 01025 if (NC_NOERR != fail) { 01026 mWriteIface->report_error("Couldn't write connectivity variable."); 01027 delete [] connectivity; 01028 return MB_FAILURE; 01029 } 01030 01031 block_index++; 01032 01033 } 01034 01035 return MB_SUCCESS; 01036 } 01037 01038 01039 ErrorCode WriteNCDF::write_global_node_order_map(int num_nodes, Range& nodes) 01040 { 01041 // note: this routine bypasses the standard exodusII interface for efficiency! 01042 01043 // node order map 01044 int* map = new int[num_nodes]; 01045 01046 // for now, output a dummy map! 01047 01048 Range::iterator range_iter, end_iter; 01049 range_iter = nodes.begin(); 01050 end_iter = nodes.end(); 01051 01052 int i=0; 01053 01054 for(; range_iter != end_iter; range_iter++) 01055 { 01056 // TODO -- do we really want to cast this to an int? 01057 map[i++] = (int)ID_FROM_HANDLE(*range_iter); 01058 } 01059 01060 // output array and cleanup 01061 01062 int error = write_exodus_integer_variable("node_num_map", 01063 map, 01064 0, 01065 num_nodes); 01066 01067 01068 if(map) 01069 delete [] map; 01070 01071 if( error < 0 ) 01072 { 01073 mWriteIface->report_error("Failed writing global node order map"); 01074 return MB_FAILURE; 01075 } 01076 01077 return MB_SUCCESS; 01078 } 01079 01080 ErrorCode WriteNCDF::write_global_element_order_map(int num_elements) 01081 { 01082 01083 // allocate map array 01084 int* map = new int[num_elements]; 01085 01086 01087 // Many Sandia codes assume this map is unique, and CUBIT does not currently 01088 // have unique ids for all elements. Therefore, to make sure nothing crashes, 01089 // insert a dummy map... 01090 01091 for(int i=0; i<num_elements; i++) 01092 { 01093 map[i] = i+1; 01094 } 01095 01096 01097 // output array and cleanup 01098 01099 int error = write_exodus_integer_variable("elem_num_map", 01100 map, 01101 0, 01102 num_elements); 01103 01104 01105 if(map) 01106 delete [] map; 01107 01108 if( error < 0 ) 01109 { 01110 mWriteIface->report_error("Failed writing global element order map"); 01111 return MB_FAILURE; 01112 } 01113 01114 return MB_SUCCESS; 01115 } 01116 01117 01118 ErrorCode WriteNCDF::write_element_order_map(int num_elements) 01119 { 01120 // note: this routine bypasses the standard exodusII interface for efficiency! 01121 01122 // element order map 01123 int* map = new int[num_elements]; 01124 01125 // for now, output a dummy map! 01126 01127 for(int i=0; i<num_elements; i++) 01128 { 01129 map[i] = i+1; 01130 } 01131 01132 // output array and cleanup 01133 01134 int error = write_exodus_integer_variable("elem_map", 01135 map, 01136 0, 01137 num_elements); 01138 01139 01140 if(map) 01141 delete [] map; 01142 01143 if( error < 0 ) 01144 { 01145 mWriteIface->report_error("Failed writing element map"); 01146 return MB_FAILURE; 01147 } 01148 01149 return MB_SUCCESS; 01150 } 01151 01152 01153 01154 01155 ErrorCode WriteNCDF::write_exodus_integer_variable(const char* variable_name, 01156 int *variable_array, 01157 int start_position, 01158 int number_values) 01159 { 01160 01161 // note: this routine bypasses the standard exodusII interface for efficiency! 01162 01163 // write directly to netcdf interface for efficiency 01164 01165 // get the variable id of the element map 01166 int nc_var = -1; 01167 std::vector<int> dims; 01168 GET_VAR(variable_name, nc_var, dims); 01169 if (-1 == nc_var) 01170 { 01171 mWriteIface->report_error("WriteNCDF: failed to locate variable %s in file.", variable_name); 01172 return MB_FAILURE; 01173 } 01174 // this contortion is necessary because netCDF is expecting nclongs; 01175 // fortunately it's necessary only when ints and nclongs aren't the same size 01176 01177 size_t start[1], count[1]; 01178 start[0] = start_position; 01179 count[0] = number_values; 01180 01181 int fail = NC_NOERR; 01182 if (sizeof(int) == sizeof(long)) { 01183 fail = nc_put_vara_int(ncFile, nc_var, start, count, variable_array); 01184 } else { 01185 long *lptr = new long[number_values]; 01186 for (int jj = 0; jj < number_values; jj++) 01187 lptr[jj] = variable_array[jj]; 01188 fail = nc_put_vara_long(ncFile, nc_var, start, count, lptr); 01189 delete [] lptr; 01190 } 01191 if (NC_NOERR != fail) 01192 { 01193 mWriteIface->report_error("Failed to store variable %s", variable_name); 01194 return MB_FAILURE; 01195 } 01196 01197 return MB_SUCCESS; 01198 } 01199 01200 01201 01202 ErrorCode WriteNCDF::write_BCs(std::vector<NeumannSetData> &sidesets, 01203 std::vector<DirichletSetData> &nodesets) 01204 { 01205 unsigned int i,j; 01206 int id; 01207 int ns_index = -1; 01208 01209 for(std::vector<DirichletSetData>::iterator ns_it = nodesets.begin(); 01210 ns_it != nodesets.end(); ns_it++) 01211 { 01212 01213 //get number of nodes in set 01214 int number_nodes = (*ns_it).number_nodes; 01215 if (0 == number_nodes) continue; 01216 01217 // if we're here, we have a non-empty nodeset; increment the index 01218 ns_index++; 01219 01220 //get the node set id 01221 id = (*ns_it).id; 01222 01223 //build new array to old exodus ids 01224 int * exodus_id_array = new int[number_nodes]; 01225 double * dist_factor_array = new double[number_nodes]; 01226 01227 std::vector<EntityHandle>::iterator begin_iter, end_iter; 01228 std::vector<double>::iterator other_iter; 01229 begin_iter = (*ns_it).nodes.begin(); 01230 end_iter = (*ns_it).nodes.end(); 01231 other_iter = (*ns_it).node_dist_factors.begin(); 01232 01233 j=0; 01234 int exodus_id; 01235 ErrorCode result; 01236 //fill up node array and dist. factor array at the same time 01237 for(; begin_iter != end_iter; begin_iter++) 01238 { 01239 result = mdbImpl->tag_get_data(mGlobalIdTag, &(*begin_iter), 1, &exodus_id); 01240 if (MB_SUCCESS != result) { 01241 mWriteIface->report_error("Problem getting id tag data."); 01242 return result; 01243 } 01244 01245 exodus_id_array[j] = exodus_id; 01246 dist_factor_array[j] = *(other_iter); 01247 other_iter++; 01248 j++; 01249 } 01250 01251 // write out the id for the nodeset 01252 01253 int num_values = 1; 01254 01255 result = write_exodus_integer_variable("ns_prop1", 01256 &id, ns_index, num_values); 01257 01258 if (result != MB_SUCCESS) 01259 { 01260 mWriteIface->report_error("Problem writing node set id %d", id); 01261 return MB_FAILURE; 01262 } 01263 01264 // write out the nodeset status 01265 01266 int status = 1; 01267 if (!number_nodes) 01268 status = 0; 01269 01270 result = write_exodus_integer_variable("ns_status", 01271 &status, ns_index, num_values); 01272 01273 if (result != MB_SUCCESS) 01274 { 01275 mWriteIface->report_error("Problem writing node set status"); 01276 return MB_FAILURE; 01277 } 01278 01279 //write it out 01280 char wname[80]; 01281 int nc_var = -1; 01282 std::vector<int> dims; 01283 INS_ID(wname, "node_ns%d", ns_index+1); 01284 GET_VAR(wname, nc_var, dims); 01285 if (-1 == nc_var) { 01286 mWriteIface->report_error("Failed to get node_ns variable."); 01287 return MB_FAILURE; 01288 } 01289 01290 size_t start = 0, count = number_nodes; 01291 int fail = nc_put_vara_int(ncFile, nc_var, &start, &count, exodus_id_array); 01292 if(NC_NOERR != fail) 01293 { 01294 mWriteIface->report_error("Failed writing exodus id array"); 01295 return MB_FAILURE; 01296 } 01297 01298 // write out nodeset distribution factors 01299 INS_ID(wname, "dist_fact_ns%d", ns_index+1); 01300 nc_var = -1; 01301 GET_VAR(wname, nc_var, dims); 01302 if (-1 == nc_var) { 01303 mWriteIface->report_error("Failed to get dist_fact variable."); 01304 return MB_FAILURE; 01305 } 01306 fail = nc_put_vara_double(ncFile, nc_var, &start, &count, dist_factor_array); 01307 if(NC_NOERR != fail) 01308 { 01309 mWriteIface->report_error("Failed writing dist factor array"); 01310 return MB_FAILURE; 01311 } 01312 01313 delete [] dist_factor_array; 01314 delete [] exodus_id_array; 01315 } 01316 01317 01318 //now do sidesets 01319 int ss_index = 0; // index of sideset - not the same as 'i' because 01320 // only writing non-empty side sets 01321 for( i=0; i<sidesets.size(); i++) 01322 { 01323 01324 NeumannSetData sideset_data = sidesets[i]; 01325 01326 //get the side set id 01327 int side_set_id = sideset_data.id; 01328 01329 //get number of elements in set 01330 int number_elements = sideset_data.number_elements; 01331 if( number_elements == 0 ) 01332 continue; 01333 01334 //build new array to old exodus ids 01335 int * output_element_ids = new int[number_elements]; 01336 int * output_element_side_numbers = new int[number_elements]; 01337 01338 std::vector<EntityHandle>::iterator begin_iter, end_iter; 01339 begin_iter = sideset_data.elements.begin(); 01340 end_iter = sideset_data.elements.end(); 01341 std::vector<int>::iterator side_iter = sideset_data.side_numbers.begin(); 01342 01343 //get the tag handle 01344 j=0; 01345 int exodus_id; 01346 01347 //for each "side" 01348 for(; begin_iter != end_iter; begin_iter++, side_iter++) 01349 { 01350 ErrorCode result = mdbImpl->tag_get_data(mGlobalIdTag, 01351 &(*begin_iter), 01352 1, &exodus_id); 01353 if (MB_FAILURE == result) { 01354 mWriteIface->report_error("Problem getting exodus id for sideset element %lu", 01355 (long unsigned int) ID_FROM_HANDLE(*begin_iter)); 01356 return result; 01357 } 01358 01359 output_element_ids[j] = exodus_id; 01360 output_element_side_numbers[j++] = *side_iter; 01361 } 01362 01363 if(number_elements) 01364 { 01365 // write out the id for the nodeset 01366 01367 int num_values = 1; 01368 01369 // ss_prop1[ss_index] = side_set_id 01370 ErrorCode result = write_exodus_integer_variable("ss_prop1", 01371 &side_set_id, 01372 ss_index,num_values); 01373 01374 if (result != MB_SUCCESS) 01375 { 01376 mWriteIface->report_error("Problem writing node set id %d", id); 01377 return MB_FAILURE; 01378 } 01379 01380 // FIXME : Something seems wrong here. The we are within a block 01381 // started with if(number_elements), so this condition is always 01382 // false. This code seems to indicate that we want to write all 01383 // sidesets, not just empty ones. But the code both here and in 01384 // initialize_exodus_file() skip empty side sets. 01385 // - j.k. 2007-03-09 01386 int status = 1; 01387 if (!number_elements) 01388 status = 0; 01389 01390 // ss_status[ss_index] = status 01391 result = write_exodus_integer_variable("ss_status", 01392 &status, 01393 ss_index, num_values); 01394 if (result != MB_SUCCESS) 01395 { 01396 mWriteIface->report_error("Problem writing side set status"); 01397 return MB_FAILURE; 01398 } 01399 01400 // Increment ss_index now because we want a) we need to 01401 // increment it somewhere within the if(number_elements) 01402 // block and b) the above calls need a zero-based index 01403 // while the following use a one-based index. 01404 ++ss_index; 01405 01406 char wname[80]; 01407 int nc_var; 01408 std::vector<int> dims; 01409 INS_ID(wname, "elem_ss%d", ss_index); 01410 GET_VAR(wname, nc_var, dims); 01411 if (-1 == nc_var) { 01412 mWriteIface->report_error("Failed to get elem_ss variable."); 01413 return MB_FAILURE; 01414 } 01415 size_t start = 0, count = number_elements; 01416 int fail = nc_put_vara_int(ncFile, nc_var, &start, &count, output_element_ids); 01417 if(NC_NOERR != fail) 01418 { 01419 mWriteIface->report_error("Failed writing sideset element array"); 01420 return MB_FAILURE; 01421 } 01422 01423 INS_ID(wname, "side_ss%d", ss_index); 01424 nc_var = -1; 01425 GET_VAR(wname, nc_var, dims); 01426 if (-1 == nc_var) { 01427 mWriteIface->report_error("Failed to get side_ss variable."); 01428 return MB_FAILURE; 01429 } 01430 fail = nc_put_vara_int(ncFile, nc_var, &start, &count, output_element_side_numbers); 01431 if(NC_NOERR != fail) 01432 { 01433 mWriteIface->report_error("Failed writing sideset side array"); 01434 return MB_FAILURE; 01435 } 01436 01437 INS_ID(wname, "dist_fact_ss%d", ss_index); 01438 nc_var = -1; 01439 GET_VAR(wname, nc_var, dims); 01440 if (-1 == nc_var) { 01441 mWriteIface->report_error("Failed to get sideset dist factors variable."); 01442 return MB_FAILURE; 01443 } 01444 count = sideset_data.ss_dist_factors.size(); 01445 fail = nc_put_vara_double(ncFile, nc_var, &start, &count, &(sideset_data.ss_dist_factors[0])); 01446 if(NC_NOERR != fail) 01447 { 01448 mWriteIface->report_error("Failed writing sideset dist factors array"); 01449 return MB_FAILURE; 01450 } 01451 } 01452 delete [] output_element_ids; 01453 delete [] output_element_side_numbers; 01454 01455 } 01456 01457 return MB_SUCCESS; 01458 } 01459 01460 ErrorCode WriteNCDF::initialize_exodus_file(ExodusMeshInfo &mesh_info, 01461 std::vector<MaterialSetData> &block_data, 01462 std::vector<NeumannSetData> & sideset_data, 01463 std::vector<DirichletSetData> & nodeset_data, 01464 const char* title_string, 01465 bool write_maps, 01466 bool /* write_sideset_distribution_factors */) 01467 { 01468 // This routine takes the place of the exodusii routine ex_put_init, 01469 // and additionally pre-defines variables such as qa, element blocks, 01470 // sidesets and nodesets in a single pass. 01471 // 01472 // This is necessary because of the way exodusII works. Each time the 01473 // netcdf routine endef is called to take the file out of define mode, 01474 // the entire file is copied before the new information is added. 01475 // 01476 // With very large files, this is simply not workable. This routine takes 01477 // the definition portions of all applicable exodus routines and puts them 01478 // in a single definition, preventing repeated copying of the file. 01479 // 01480 // Most of the code is copied directly from the applicable exodus routine, 01481 // and thus this routine may not seem very consistent in usage or variable 01482 // naming, etc. 01483 01484 // perform the initializations 01485 01486 int element_block_index; 01487 01488 // inquire on defined string dimension and general dimension for qa 01489 01490 int dim_str, dim_four, dim_line, dim_time; 01491 if (nc_def_dim(ncFile, "len_string", ExoIIInterface::MAX_STR_LENGTH, &dim_str) != NC_NOERR) 01492 { 01493 mWriteIface->report_error("WriteNCDF: failed to get string length in file"); 01494 return (MB_FAILURE); 01495 } 01496 if (nc_def_dim(ncFile, "len_line", ExoIIInterface::MAX_STR_LENGTH, &dim_line) != NC_NOERR) 01497 { 01498 mWriteIface->report_error("WriteNCDF: failed to get line length in file"); 01499 return (MB_FAILURE); 01500 } 01501 if (nc_def_dim(ncFile, "four", 4, &dim_four) != NC_NOERR) 01502 { 01503 mWriteIface->report_error("WriteNCDF: failed to locate four in file"); 01504 return (MB_FAILURE); 01505 } 01506 if (nc_def_dim(ncFile, "time_step", 1, &dim_time) != NC_NOERR) 01507 { 01508 mWriteIface->report_error("WriteNCDF: failed to locate time step in file"); 01509 return (MB_FAILURE); 01510 } 01511 01512 01513 01514 /* put file into define mode */ 01515 01516 // it is possible that an NT filename using backslashes is in the title string 01517 // this can give fits to unix codes where the backslash is assumed to be an escape 01518 // sequence. For the exodus file, just flip the backslash to a slash to prevent 01519 // this problem 01520 01521 // get a working copy of the title_string; 01522 01523 char working_title[80]; 01524 strncpy(working_title,title_string, 79); 01525 01526 int length = strlen(working_title); 01527 for(int pos = 0; pos < length; pos++) 01528 { 01529 if (working_title[pos] == '\\') 01530 strncpy(&working_title[pos],"/",1); 01531 } 01532 01533 if (NC_NOERR != nc_put_att_text(ncFile, NC_GLOBAL, "title", length, working_title)) 01534 { 01535 mWriteIface->report_error("WriteNCDF: failed to define title attribute"); 01536 return (MB_FAILURE); 01537 } 01538 01539 // add other attributes while we're at it 01540 float dum_vers = 3.04F; 01541 if (NC_NOERR != nc_put_att_float(ncFile, NC_GLOBAL, "api_version", NC_FLOAT, 1, &dum_vers)) { 01542 mWriteIface->report_error("WriteNCDF: failed to define api_version attribute"); 01543 return (MB_FAILURE); 01544 } 01545 dum_vers = 2.05F; 01546 if (NC_NOERR != nc_put_att_float(ncFile, NC_GLOBAL, "version", NC_FLOAT, 1, &dum_vers)) { 01547 mWriteIface->report_error("WriteNCDF: failed to define version attribute"); 01548 return (MB_FAILURE); 01549 } 01550 int dum_siz = sizeof(double); 01551 if (NC_NOERR != nc_put_att_int(ncFile, NC_GLOBAL, "floating_point_word_size", NC_INT, 1, &dum_siz)) { 01552 mWriteIface->report_error("WriteNCDF: failed to define floating pt word size attribute"); 01553 return (MB_FAILURE); 01554 } 01555 01556 // set up number of dimensions 01557 01558 int num_el_blk, num_elem, num_nodes, num_dim; 01559 if (nc_def_dim(ncFile, "num_dim", (size_t)mesh_info.num_dim, &num_dim) != NC_NOERR) 01560 { 01561 mWriteIface->report_error("WriteNCDF: failed to define number of dimensions"); 01562 return (MB_FAILURE); 01563 } 01564 01565 if (nc_def_dim(ncFile, "num_nodes", mesh_info.num_nodes, &num_nodes) != NC_NOERR) 01566 { 01567 mWriteIface->report_error("WriteNCDF: failed to define number of nodes"); 01568 return (MB_FAILURE); 01569 } 01570 01571 if (nc_def_dim(ncFile, "num_elem", mesh_info.num_elements, &num_elem) != NC_NOERR) 01572 { 01573 mWriteIface->report_error("WriteNCDF: failed to define number of elements"); 01574 return (MB_FAILURE); 01575 } 01576 01577 if (nc_def_dim(ncFile, "num_el_blk", mesh_info.num_elementblocks, &num_el_blk) != NC_NOERR) 01578 { 01579 mWriteIface->report_error("WriteNCDF: failed to define number of element blocks"); 01580 return (MB_FAILURE); 01581 } 01582 01583 /* ...and some variables */ 01584 01585 /* element block id status array */ 01586 int idstat = -1; 01587 if (NC_NOERR != nc_def_var(ncFile, "eb_status", NC_LONG, 1, &num_el_blk, &idstat)) 01588 { 01589 mWriteIface->report_error("WriteNCDF: failed to define element block status array"); 01590 return (MB_FAILURE); 01591 } 01592 01593 /* element block id array */ 01594 01595 int idarr = -1; 01596 if (NC_NOERR != nc_def_var(ncFile, "eb_prop1", NC_LONG, 1, &num_el_blk, &idarr)) 01597 { 01598 mWriteIface->report_error("WriteNCDF: failed to define element block id array"); 01599 return (MB_FAILURE); 01600 } 01601 01602 /* store property name as attribute of property array variable */ 01603 if (NC_NOERR != nc_put_att_text(ncFile, idarr, "name", strlen("ID"), "ID")) 01604 { 01605 mWriteIface->report_error("WriteNCDF: failed to store element block property name ID"); 01606 return (MB_FAILURE); 01607 } 01608 01609 // define element blocks 01610 01611 01612 char wname[80]; 01613 for (unsigned int i = 0; i < block_data.size(); i++) 01614 { 01615 MaterialSetData block = block_data[i]; 01616 01617 element_block_index = i+1; 01618 int num_el_in_blk = -1, num_att_in_blk = -1; 01619 int blk_attrib, connect; 01620 01621 /* define number of elements in this block */ 01622 01623 INS_ID(wname, "num_el_in_blk%d", element_block_index); 01624 if (nc_def_dim(ncFile, wname, (size_t)block.number_elements, &num_el_in_blk) != NC_NOERR) 01625 { 01626 mWriteIface->report_error("WriteNCDF: failed to define number of elements/block for block %d", 01627 i+1); 01628 return (MB_FAILURE); 01629 } 01630 01631 01632 /* define number of nodes per element for this block */ 01633 INS_ID(wname, "num_nod_per_el%d", element_block_index); 01634 int num_nod_per_el = -1; 01635 if (nc_def_dim(ncFile, wname, (size_t)block.number_nodes_per_element, &num_nod_per_el) != NC_NOERR) 01636 { 01637 mWriteIface->report_error("WriteNCDF: failed to define number of nodes/element for block %d", 01638 block.id); 01639 return (MB_FAILURE); 01640 } 01641 01642 /* define element attribute array for this block */ 01643 int dims[3]; 01644 if (block.number_attributes > 0) 01645 { 01646 INS_ID(wname, "num_att_in_blk%d", element_block_index); 01647 if (nc_def_dim(ncFile, wname, (size_t)block.number_attributes, &num_att_in_blk) != NC_NOERR) 01648 { 01649 mWriteIface->report_error("WriteNCDF: failed to define number of attributes in block %d", block.id); 01650 return (MB_FAILURE); 01651 } 01652 01653 INS_ID(wname, "attrib%d", element_block_index); 01654 dims[0] = num_el_in_blk; 01655 dims[1] = num_att_in_blk; 01656 if (NC_NOERR != nc_def_var(ncFile, wname, NC_DOUBLE, 2, dims, &blk_attrib)) 01657 { 01658 mWriteIface->report_error("WriteNCDF: failed to define attributes for element block %d", block.id); 01659 return (MB_FAILURE); 01660 } 01661 } 01662 01663 01664 /* define element connectivity array for this block */ 01665 01666 INS_ID(wname, "connect%d", element_block_index); 01667 dims[0] = num_el_in_blk; 01668 dims[1] = num_nod_per_el; 01669 if (NC_NOERR != nc_def_var(ncFile, wname, NC_LONG, 2, dims, &connect)) 01670 { 01671 mWriteIface->report_error("WriteNCDF: failed to create connectivity array for block %d", i+1); 01672 return (MB_FAILURE); 01673 } 01674 01675 /* store element type as attribute of connectivity variable */ 01676 01677 std::string element_type_string(ExoIIUtil::ElementTypeNames[ block.element_type ]); 01678 if (NC_NOERR != nc_put_att_text(ncFile, connect, "elem_type", element_type_string.length(), element_type_string.c_str())) 01679 { 01680 mWriteIface->report_error("WriteNCDF: failed to store element type name %d", (int)block.element_type); 01681 return (MB_FAILURE); 01682 } 01683 } 01684 01685 01686 /* node set id array: */ 01687 01688 int non_empty_nss = 0; 01689 // need to go through nodesets to compute # nodesets, some might be empty 01690 std::vector<DirichletSetData>::iterator ns_it; 01691 for( ns_it = nodeset_data.begin(); 01692 ns_it != nodeset_data.end(); ns_it++) { 01693 if (0 != (*ns_it).number_nodes) non_empty_nss++; 01694 } 01695 01696 int num_ns = -1; 01697 int ns_idstat = -1, ns_idarr = -1; 01698 if (non_empty_nss > 0) 01699 { 01700 if (nc_def_dim(ncFile, "num_node_sets", (size_t)(non_empty_nss), &num_ns) != NC_NOERR) 01701 { 01702 mWriteIface->report_error("WriteNCDF: failed to define number of node sets"); 01703 return (MB_FAILURE); 01704 } 01705 01706 /* node set id status array: */ 01707 01708 if (NC_NOERR != nc_def_var(ncFile, "ns_status", NC_LONG, 1, &num_ns, &ns_idstat)) 01709 { 01710 mWriteIface->report_error("WriteNCDF: failed to create node sets status array"); 01711 return (MB_FAILURE); 01712 } 01713 01714 /* node set id array: */ 01715 if (NC_NOERR != nc_def_var(ncFile, "ns_prop1", NC_LONG, 1, &num_ns, &ns_idarr)) 01716 { 01717 mWriteIface->report_error("WriteNCDF: failed to create node sets property array"); 01718 return (MB_FAILURE); 01719 } 01720 01721 01722 /* store property name as attribute of property array variable */ 01723 if (NC_NOERR != nc_put_att_text(ncFile, NC_GLOBAL, "name", strlen("ID"), "ID")) 01724 { 01725 mWriteIface->report_error("WriteNCDF: failed to store node set property name ID"); 01726 return (MB_FAILURE); 01727 } 01728 01729 // now, define the arrays needed for each node set 01730 01731 01732 int index = 0; 01733 01734 for(unsigned i = 0; i <nodeset_data.size(); i++) 01735 { 01736 DirichletSetData node_set = nodeset_data[i]; 01737 01738 if (node_set.number_nodes == 0) { 01739 mWriteIface->report_error("WriteNCDF: empty nodeset %d", node_set.id); 01740 continue; 01741 } 01742 index++; 01743 01744 int num_nod_ns = -1; 01745 INS_ID(wname, "num_nod_ns%d", index); 01746 if (nc_def_dim(ncFile, wname, (size_t)node_set.number_nodes, &num_nod_ns) != NC_NOERR) 01747 { 01748 mWriteIface->report_error("WriteNCDF: failed to define number of nodes for set %d", node_set.id); 01749 return (MB_FAILURE); 01750 } 01751 01752 /* create variable array in which to store the node set node list */ 01753 int node_ns = -1; 01754 INS_ID(wname, "node_ns%d", index); 01755 if (NC_NOERR != nc_def_var(ncFile, wname, NC_LONG, 1, &num_nod_ns, &node_ns)) 01756 { 01757 mWriteIface->report_error("WriteNCDF: failed to create node set %d node list", node_set.id); 01758 return (MB_FAILURE); 01759 } 01760 01761 // create distribution factor array 01762 int fact_ns = -1; 01763 INS_ID(wname, "dist_fact_ns%d", index); 01764 if (NC_NOERR != nc_def_var(ncFile, wname, NC_DOUBLE, 1, &num_nod_ns, &fact_ns)) 01765 { 01766 mWriteIface->report_error("WriteNCDF: failed to create node set %d distribution factor list", 01767 node_set.id); 01768 return (MB_FAILURE); 01769 } 01770 01771 } 01772 01773 } 01774 01775 /* side set id array: */ 01776 01777 long non_empty_ss = 0; 01778 // need to go through nodesets to compute # nodesets, some might be empty 01779 std::vector<NeumannSetData>::iterator ss_it; 01780 for( ss_it = sideset_data.begin(); 01781 ss_it != sideset_data.end(); ss_it++) { 01782 if (0 != (*ss_it).number_elements) non_empty_ss++; 01783 } 01784 01785 if (non_empty_ss > 0) { 01786 int num_ss = -1; 01787 if (nc_def_dim(ncFile, "num_side_sets", non_empty_ss, &num_ss) != NC_NOERR) 01788 { 01789 mWriteIface->report_error("WriteNCDF: failed to define number of side sets"); 01790 return (MB_FAILURE); 01791 } 01792 01793 /* side set id status array: */ 01794 int ss_idstat = -1, ss_idarr = -1; 01795 if (NC_NOERR != nc_def_var(ncFile, "ss_status", NC_LONG, 1, &num_ss, &ss_idstat)) 01796 { 01797 mWriteIface->report_error("WriteNCDF: failed to define side set status"); 01798 return (MB_FAILURE); 01799 } 01800 01801 /* side set id array: */ 01802 if (NC_NOERR != nc_def_var(ncFile, "ss_prop1", NC_LONG, 1, &num_ss, &ss_idarr)) 01803 { 01804 mWriteIface->report_error( "WriteNCDF: failed to define side set property"); 01805 return (MB_FAILURE); 01806 } 01807 01808 /* store property name as attribute of property array variable */ 01809 if (NC_NOERR != nc_put_att_text(ncFile, ss_idarr, "name", strlen("ID"), "ID")) 01810 { 01811 mWriteIface->report_error("WriteNCDF: failed to store side set property name ID"); 01812 return (MB_FAILURE); 01813 } 01814 01815 // now, define the arrays needed for each side set 01816 01817 int index = 0; 01818 for(unsigned int i = 0; i < sideset_data.size(); i++) 01819 { 01820 NeumannSetData side_set = sideset_data[i]; 01821 01822 // dont define an empty set 01823 if (side_set.number_elements == 0) 01824 continue; 01825 01826 index++; 01827 01828 int num_side_ss = -1; 01829 int elem_ss = -1, side_ss = -1; 01830 INS_ID(wname, "num_side_ss%d", index); 01831 if (nc_def_dim(ncFile, wname, (size_t)side_set.number_elements, &num_side_ss) != NC_NOERR) 01832 { 01833 mWriteIface->report_error("WriteNCDF: failed to define number of sides in side set %d", 01834 side_set.id); 01835 return(MB_FAILURE); 01836 } 01837 01838 INS_ID(wname, "elem_ss%d", index); 01839 if (NC_NOERR != nc_def_var(ncFile, wname, NC_LONG, 1, &num_side_ss, &elem_ss)) 01840 { 01841 mWriteIface->report_error("WriteNCDF: failed to create element list for side set %d", 01842 side_set.id); 01843 return(MB_FAILURE); /* exit define mode and return */ 01844 } 01845 INS_ID(wname, "side_ss%d", index); 01846 if (NC_NOERR != nc_def_var(ncFile, wname, NC_LONG, 1, &num_side_ss, &side_ss)) 01847 { 01848 mWriteIface->report_error("WriteNCDF: failed to create side list for side set %d", 01849 side_set.id); 01850 return(MB_FAILURE); /* exit define mode and return */ 01851 01852 } 01853 01854 // sideset distribution factors 01855 int num_df_ss = -1; 01856 INS_ID(wname, "num_df_ss%d", index); 01857 if (nc_def_dim(ncFile, wname, (size_t)side_set.ss_dist_factors.size(), &num_df_ss) != NC_NOERR) 01858 { 01859 mWriteIface->report_error("WriteNCDF: failed to define number of dist factors in side set %d", 01860 side_set.id); 01861 return(MB_FAILURE); /* exit define mode and return */ 01862 } 01863 01864 /* create variable array in which to store the side set distribution factors */ 01865 01866 int fact_ss = -1; 01867 INS_ID(wname, "dist_fact_ss%d", index); 01868 if (NC_NOERR != nc_def_var(ncFile, wname, NC_LONG, 1, &num_df_ss, &fact_ss)) 01869 { 01870 mWriteIface->report_error("WriteNCDF: failed to create dist factors list for side set %d", 01871 side_set.id); 01872 return(MB_FAILURE); /* exit define mode and return */ 01873 } 01874 } 01875 } 01876 01877 /* node coordinate arrays: */ 01878 01879 int coord, name_coord, dims[3]; 01880 dims[0] = num_dim; 01881 dims[1] = num_nodes; 01882 if (NC_NOERR != nc_def_var(ncFile, "coord", NC_DOUBLE, 2, dims, &coord)) 01883 { 01884 mWriteIface->report_error("WriteNCDF: failed to define node coordinate array"); 01885 return (MB_FAILURE); 01886 } 01887 /* coordinate names array */ 01888 01889 dims[0] = num_dim; 01890 dims[1] = dim_str; 01891 if (NC_NOERR != nc_def_var(ncFile, "coor_names", NC_CHAR, 2, dims, &name_coord)) 01892 { 01893 mWriteIface->report_error("WriteNCDF: failed to define coordinate name array"); 01894 return (MB_FAILURE); 01895 } 01896 01897 // define genesis maps if required 01898 01899 if (write_maps) 01900 { 01901 // element map 01902 int elem_map = -1, elem_map2 = -1, node_map = -1; 01903 if (NC_NOERR != nc_def_var(ncFile, "elem_map", NC_LONG, 1, &num_elem, &elem_map)) 01904 { 01905 mWriteIface->report_error("WriteNCDF: failed to create element map array"); 01906 return (MB_FAILURE); /* exit define mode and return */ 01907 } 01908 01909 // create the element number map 01910 if (NC_NOERR != nc_def_var(ncFile, "elem_num_map", NC_LONG, 1, &num_elem, &elem_map2)) 01911 { 01912 mWriteIface->report_error("WriteNCDF: failed to create element numbering map"); 01913 } 01914 // create node number map 01915 if (NC_NOERR != nc_def_var(ncFile, "node_num_map", NC_LONG, 1, &num_nodes, &node_map)) 01916 { 01917 mWriteIface->report_error("WriteNCDF: failed to create node numbering map array"); 01918 return(MB_FAILURE); /* exit define mode and return */ 01919 } 01920 } 01921 01922 // define qa records to be used 01923 01924 int num_qa_rec = mesh_info.qaRecords.size()/4; 01925 int num_qa = -1; 01926 01927 if (nc_def_dim(ncFile, "num_qa_rec",(long)num_qa_rec, &num_qa) != NC_NOERR) 01928 { 01929 mWriteIface->report_error("WriteNCDF: failed to define qa record array size"); 01930 return (MB_FAILURE); 01931 } 01932 01933 // define qa array 01934 int qa_title; 01935 dims[0] = num_qa; 01936 dims[1] = dim_four; 01937 dims[2] = dim_str; 01938 if (NC_NOERR != nc_def_var(ncFile, "qa_records", NC_CHAR, 3, dims, &qa_title)) 01939 { 01940 mWriteIface->report_error("WriteNCDF: failed to define qa record array"); 01941 return (MB_FAILURE); 01942 } 01943 01944 // take it out of define mode 01945 if (NC_NOERR != nc_enddef(ncFile)) { 01946 mWriteIface->report_error("WriteNCDF: Trouble leaving define mode."); 01947 return (MB_FAILURE); 01948 } 01949 01950 return MB_SUCCESS; 01951 } 01952 01953 01954 ErrorCode WriteNCDF::open_file(const char* filename) 01955 { 01956 // not a valid filname 01957 if(strlen((const char*)filename) == 0) 01958 { 01959 mWriteIface->report_error("Output Exodus filename not specified"); 01960 return MB_FAILURE; 01961 } 01962 01963 int fail = nc_create(filename, NC_CLOBBER, &ncFile); 01964 01965 // file couldn't be opened 01966 if (NC_NOERR != fail) 01967 { 01968 mWriteIface->report_error("Cannot open %s", filename); 01969 return MB_FAILURE; 01970 } 01971 return MB_SUCCESS; 01972 } 01973 01974 ErrorCode WriteNCDF::get_sideset_elems(EntityHandle sideset, int current_sense, 01975 Range &forward_elems, Range &reverse_elems) 01976 { 01977 Range ss_elems, ss_meshsets; 01978 01979 // get the sense tag; don't need to check return, might be an error if the tag 01980 // hasn't been created yet 01981 Tag sense_tag = 0; 01982 mdbImpl->tag_get_handle("SENSE", 1, MB_TYPE_INTEGER, sense_tag); 01983 01984 // get the entities in this set 01985 ErrorCode result = mdbImpl->get_entities_by_handle(sideset, ss_elems, true); 01986 if (MB_FAILURE == result) return result; 01987 01988 // now remove the meshsets into the ss_meshsets; first find the first meshset, 01989 Range::iterator range_iter = ss_elems.begin(); 01990 while (TYPE_FROM_HANDLE(*range_iter) != MBENTITYSET && range_iter != ss_elems.end()) 01991 range_iter++; 01992 01993 // then, if there are some, copy them into ss_meshsets and erase from ss_elems 01994 if (range_iter != ss_elems.end()) { 01995 std::copy(range_iter, ss_elems.end(), range_inserter(ss_meshsets)); 01996 ss_elems.erase(range_iter, ss_elems.end()); 01997 } 01998 01999 02000 // ok, for the elements, check the sense of this set and copy into the right range 02001 // (if the sense is 0, copy into both ranges) 02002 02003 // need to step forward on list until we reach the right dimension 02004 Range::iterator dum_it = ss_elems.end(); 02005 dum_it--; 02006 int target_dim = CN::Dimension(TYPE_FROM_HANDLE(*dum_it)); 02007 dum_it = ss_elems.begin(); 02008 while (target_dim != CN::Dimension(TYPE_FROM_HANDLE(*dum_it)) && 02009 dum_it != ss_elems.end()) 02010 dum_it++; 02011 02012 if (current_sense == 1 || current_sense == 0) 02013 std::copy(dum_it, ss_elems.end(), range_inserter(forward_elems)); 02014 if (current_sense == -1 || current_sense == 0) 02015 std::copy(dum_it, ss_elems.end(), range_inserter(reverse_elems)); 02016 02017 // now loop over the contained meshsets, getting the sense of those and calling this 02018 // function recursively 02019 for (range_iter = ss_meshsets.begin(); range_iter != ss_meshsets.end(); range_iter++) { 02020 02021 // first get the sense; if it's not there, by convention it's forward 02022 int this_sense; 02023 if (0 == sense_tag || 02024 MB_FAILURE == mdbImpl->tag_get_data(sense_tag, &(*range_iter), 1, &this_sense)) 02025 this_sense = 1; 02026 02027 // now get all the entities on this meshset, with the proper (possibly reversed) sense 02028 get_sideset_elems(*range_iter, this_sense*current_sense, 02029 forward_elems, reverse_elems); 02030 } 02031 02032 return result; 02033 } 02034 02035 02036 } // namespace moab 02037 02038