[MOAB-dev] r3559 - in MOAB/trunk: . tools/dagmc
bmsmith6 at wisc.edu
bmsmith6 at wisc.edu
Mon Mar 1 11:02:28 CST 2010
Author: bmsmith
Date: 2010-03-01 11:02:28 -0600 (Mon, 01 Mar 2010)
New Revision: 3559
Modified:
MOAB/trunk/ReadNCDF.cpp
MOAB/trunk/tools/dagmc/cub2h5m.cc
Log:
These changes/bug fixes were made to support the 5-pin space reactor model in DagMC.
ReadNCDF:
Create a map for fast handle_by_id for nodes
Locate dead elements by node connectivity.
cub2h5m:
Fix special case of a planar surface through the origin for the signed volume calculation
Remove surfaces (then volumes) if all of their faces are dead elements to avoid OBB problems
Modified: MOAB/trunk/ReadNCDF.cpp
===================================================================
--- MOAB/trunk/ReadNCDF.cpp 2010-02-26 20:57:00 UTC (rev 3558)
+++ MOAB/trunk/ReadNCDF.cpp 2010-03-01 17:02:28 UTC (rev 3559)
@@ -28,6 +28,7 @@
#include <cmath>
#include <memory>
#include <sstream>
+#include <map>
#include "MBCN.hpp"
#include "MBRange.hpp"
@@ -1614,7 +1615,10 @@
//4. loop through the node_num_map, use it to find the node in the cub file.
//5. Replace coord[0][n] with coordx[m]+vals_nod_var1(time_step, m) for all directions for matching nodes.
// Test: try to get hexes
-
+
+ // *******************************************************************
+ // Move nodes to their deformed locations.
+ // *******************************************************************
MBErrorCode rval;
std::string s;
rval = opts.get_str_option("tdata", s );
@@ -1821,10 +1825,18 @@
double max_magnitude = 0;
double average_magnitude = 0;
+ // Find the matching cub node by id or proximity
+ const bool match_node_ids = true;
+
+ // We should not use cub verts unless they have been matched. Place in a map
+ // for fast handle_by_id lookup.
+ std::map<int,MBEntityHandle> matched_cub_vert_id_map;
+
// For each exo vert find the closest cub vert. Use proximity because kd-tree
// search is O(logn) but MOAB tag search is O(n). For 150k nodes this
- // is 5 minutes faster.
- const double GEOMETRY_RESABS = 1e-6;
+ // is 5 minutes faster. The MAX_NODE_DIST is the farthest that we will search
+ // for a node. Note that the exodus file is single precision.
+ const double MAX_NODE_DIST = 1e-3;
std::cout << " exodus file contains " << numberNodes_loading << " nodes."
<< std::endl;
for(int i=0; i<numberNodes_loading; ++i) {
@@ -1832,29 +1844,46 @@
std::vector<MBEntityHandle> leaves;
int exo_id = ptr[i];
MBCartVect exo_coords( orig_coords[0][i], orig_coords[1][i], orig_coords[2][i] );
- rval = kdtree.leaves_within_distance( root, exo_coords.array(), GEOMETRY_RESABS, leaves );
+
+ double min_dist = MAX_NODE_DIST;
+
+ rval = kdtree.leaves_within_distance( root, exo_coords.array(), MAX_NODE_DIST, leaves );
if(MB_SUCCESS != rval) return rval;
for(std::vector<MBEntityHandle>::const_iterator j=leaves.begin();
j!=leaves.end(); ++j) {
- if(0 != cub_vert) break;
+ //if(0 != cub_vert) break;
std::vector<MBEntityHandle> leaf_verts;
rval = mdbImpl->get_entities_by_type( *j, MBVERTEX, leaf_verts );
if(MB_SUCCESS != rval) return rval;
- // Find the matching cub vert by id, then do a sanity check on distance
+
+ // Find the matching cub node by id or proximity
for(std::vector<MBEntityHandle>::const_iterator k=leaf_verts.begin();
k!=leaf_verts.end(); ++k) {
- int cub_id;
- rval = mdbImpl->tag_get_data( mGlobalIdTag, &(*k), 1, &cub_id );
- if(MB_SUCCESS != rval) return rval;
- if(exo_id == cub_id) {
- cub_vert = *k;
- break;
- }
+ if(match_node_ids) {
+ int cub_id;
+ rval = mdbImpl->tag_get_data( mGlobalIdTag, &(*k), 1, &cub_id );
+ if(MB_SUCCESS != rval) return rval;
+ if(exo_id == cub_id) {
+ cub_vert = *k;
+ break;
+ }
+ } else {
+ MBCartVect orig_cub_coords, difference;
+ rval = mdbImpl->get_coords( &(*k), 1, orig_cub_coords.array() );
+ if(MB_SUCCESS != rval) return rval;
+ difference = orig_cub_coords - exo_coords;
+ double dist = difference.length();
+ if(dist < min_dist) {
+ min_dist = dist;
+ cub_vert = *k;
+ }
+ }
}
}
// If a cub vert is found, update it with the deformed coords from the exo file.
if(0 != cub_vert) {
MBCartVect updated_exo_coords;
+ matched_cub_vert_id_map[exo_id] = cub_vert;
updated_exo_coords[0] = orig_coords[0][i] + deformed_arrays[0][i];
updated_exo_coords[1] = orig_coords[1][i] + deformed_arrays[1][i];
if(numberDimensions_loading == 3 )
@@ -1874,14 +1903,28 @@
}
std::cout << " " << found
- << " nodes from the exodus file were matched by id in the cub_file_set."
- << std::endl;
- std::cout << " " << lost << " nodes from the exodus file could not be matched."
- << std::endl;
- std::cout << " maximum node displacement magnitude=" << max_magnitude << std::endl;
- std::cout << " average node displacement magnitude=" << average_magnitude/found
- << std::endl;
+ << " nodes from the exodus file were matched in the cub_file_set ";
+ if(match_node_ids) {
+ std::cout << "by id." << std::endl;
+ } else {
+ std::cout << "by proximity." << std::endl;
+ }
+ // Fail if all of the nodes could not be matched.
+ if(0 != lost) {
+ std::cout << "Error: " << lost << " nodes from the exodus file could not be matched."
+ << std::endl;
+ return MB_FAILURE;
+ }
+ std::cout << " maximum node displacement magnitude: " << max_magnitude
+ << " cm" << std::endl;
+ std::cout << " average node displacement magnitude: " << average_magnitude/found
+ << " cm" << std::endl;
+
+ // *******************************************************************
+ // Remove dead elements from the MOAB instance.
+ // *******************************************************************
+
// How many element variables are in the file?
ncdim = ncFile->get_dim("num_elem_var");
if (!ncdim->is_valid()) {
@@ -1915,15 +1958,17 @@
bool found_death_index = false;
for(int i=0; i<n_elem_var; ++i) {
std::string temp(names[i]);
- std::string::size_type pos = temp.find("death_status");
- if(std::string::npos != pos) {
+ std::string::size_type pos0 = temp.find("death_status");
+ std::string::size_type pos1 = temp.find("Death_Status");
+ std::string::size_type pos2 = temp.find("DEATH_STATUS");
+ if(std::string::npos!=pos0 || std::string::npos!=pos1 || std::string::npos!=pos2) {
found_death_index = true;
death_index = i+1; // NetCDF variables start with 1
break;
}
}
if(!found_death_index) {
- readMeshIface->report_error("ReadNCDF:: Problem getting index of death_status variable.");
+ std::cout << "ReadNCDF: Problem getting index of death_status variable." << std::endl;
return MB_FAILURE;
}
@@ -1933,87 +1978,189 @@
// Dead elements are listed by block. Read the block headers to determine how
// many elements are in each block.
rval = read_block_headers(blocks_to_load, num_blocks);
- if (MB_FAILURE == rval) return rval;
+ if (MB_FAILURE == rval) {
+ std::cout << "ReadNCDF: Problem reading block headers." << std::endl;
+ return rval;
+ }
+ // Dead elements from the Exodus file can be located in the cub_file_set by id
+ // or by connectivity. Currently, finding elements by id requires careful book
+ // keeping when constructing the model in Cubit. To avoid this, one can match
+ // elements by connectivity instead.
+ const bool match_elems_by_connectivity = true;
+
// Get the element id map. The ids in the map are from the elementsin the blocks.
// elem_num_map( blk1 elem ids, blk2 elem ids, blk3 elem ids, ... )
std::vector<int> elem_ids(numberNodes_loading);
- temp_var = ncFile->get_var("elem_num_map");
- if (NULL == temp_var || !temp_var->is_valid()) {
- readMeshIface->report_error("ReadNCDF:: Problem getting element number map variable.");
- return MB_FAILURE;
+ if(!match_elems_by_connectivity) {
+ temp_var = ncFile->get_var("elem_num_map");
+ if (NULL == temp_var || !temp_var->is_valid()) {
+ std::cout << "ReadNCDF: Problem getting element number map variable." << std::endl;
+ return MB_FAILURE;
+ }
+ status = temp_var->get(&elem_ids[0], numberElements_loading);
+ if (0 == status) {
+ std::cout << "ReadNCDF: Problem getting element number map data." << std::endl;
+ return MB_FAILURE;
+ }
}
- status = temp_var->get(&elem_ids[0], numberElements_loading);
- if (0 == status) {
- readMeshIface->report_error("ReadNCDF:: Problem getting element number map data.");
- return MB_FAILURE;
- }
- // For each block, get the death_status at the correct time step.
+ // For each block
int first_elem_id_in_block = 0;
int block_count = 1; // NetCDF variables start with 1
for(std::vector<ReadBlockData>::iterator i=blocksLoading.begin();
i!=blocksLoading.end(); ++i) {
+
+ // get the ncdf connect variable
+ std::string temp_string("connect");
+ std::stringstream temp_ss;
+ temp_ss << block_count;
+ temp_string += temp_ss.str();
+ temp_string += "\0";
+ temp_var = ncFile->get_var(temp_string.c_str());
+ if (NULL == temp_var || !temp_var->is_valid()) {
+ std::cout << "MBCN: Problem getting connect variable." << std::endl;
+ return MB_FAILURE;
+ }
+ // the element type is an attribute of the connectivity variable
+ NcAtt *temp_att = temp_var->get_att("elem_type");
+ if (NULL == temp_att || !temp_att->is_valid()) {
+ std::cout << "MBCN:: Problem getting elem type attribute." << std::endl;
+ return MB_FAILURE;
+ }
+ // Get the MOAB element type from the Exodus attribute type
+ char *dum_str = temp_att->as_string(0);
+ ExoIIElementType elem_type = ExoIIUtil::static_element_name_to_type(dum_str);
+ delete [] dum_str;
+ (*i).elemType = elem_type;
+ const MBEntityType mb_type = ExoIIUtil::ExoIIElementMBEntity[(*i).elemType];
+
+ // Get the number of nodes per element
+ unsigned int nodes_per_element = ExoIIUtil::VerticesPerElement[(*i).elemType];
+
+ // read the connectivity into that memory, this will take only part of the array
+ // 1/2 if sizeof(MBEntityHandle) == 64 bits.
+ int exo_conn[i->numElements][nodes_per_element];
+ NcBool status = temp_var->get( &exo_conn[0][0], i->numElements, nodes_per_element);
+ if (0 == status) {
+ std::cout << "MBCN: Problem getting connectivity." << std::endl;
+ return MB_FAILURE;
+ }
+
+ // get the death_status at the correct time step.
std::vector<double> death_status(i->numElements); // it seems wrong, but it uses doubles
std::string array_name("vals_elem_var");
- std::stringstream temp;
- temp << death_index;
- array_name += temp.str();
+ temp_ss.str(""); // stringstream won't clear by temp.clear()
+ temp_ss << death_index;
+ array_name += temp_ss.str();
array_name += "eb";
- temp.str(""); // stringstream won't clear by temp.clear()
- temp << block_count;
- array_name += temp.str();
+ temp_ss.str(""); // stringstream won't clear by temp.clear()
+ temp_ss << block_count;
+ array_name += temp_ss.str();
array_name += "\0";
temp_var = ncFile->get_var(array_name.c_str());
if (NULL == temp_var || !temp_var->is_valid()) {
- std::cout << "ReadNCDF:: Problem getting death_status variable." << std::endl;;
+ std::cout << "ReadNCDF: Problem getting death_status variable." << std::endl;;
return MB_FAILURE;
}
status = temp_var->set_cur(time_step-1, 0);
if (0 == status) {
- readMeshIface->report_error("ReadNCDF:: Problem setting time step for death_status.");
+ std::cout << "ReadNCDF: Problem setting time step for death_status." << std::endl;
return MB_FAILURE;
}
status = temp_var->get(&death_status[0], 1, i->numElements);
if (0 == status) {
- readMeshIface->report_error("ReadNCDF:: Problem getting death_status array.");
+ std::cout << "ReadNCDF: Problem getting death_status array." << std::endl;
return MB_FAILURE;
}
// Look for dead elements. If there is too many dead elements and this starts
// to take too long, I should put the elems in a kd-tree for more efficient
// searching. Alternatively I could get the exo connectivity and match nodes.
- int dead_elem_counter = 0;
+ int dead_elem_counter = 0, missing_elem_counter = 0;
for (int j=0; j<i->numElements; ++j) {
if (1 != death_status[j]) {
- // get dead element's id
- int elem_id = elem_ids[first_elem_id_in_block+j];
- void *id[] = {&elem_id};
- // assume hex elems!!!!!!!!
- MBRange cub_elem;
- rval = mdbImpl->get_entities_by_type_and_tag( cub_file_set, MBHEX,
- &mGlobalIdTag, id, 1, cub_elem,
- MBInterface::INTERSECT );
- if(MB_SUCCESS != rval) return rval;
- if(1 != cub_elem.size()) {
- std::cout << "ReadNCDF: Should have found 1 element, but instead found "
+
+ MBRange cub_elem, cub_nodes;
+ if(match_elems_by_connectivity) {
+ // get exodus nodes for the element
+ std::vector<int> elem_conn(nodes_per_element);
+ for(unsigned int k=0; k<nodes_per_element; ++k) {
+ elem_conn[k] = exo_conn[j][k];
+ }
+ // get the ids of the nodes (assume we are matching by id)
+ // Remember that the exodus array locations start with 1 (not 0).
+ std::vector<int> elem_conn_node_ids(nodes_per_element);
+ for(unsigned int k=0; k<nodes_per_element; ++k) {
+ elem_conn_node_ids[k] = ptr[ elem_conn[k]-1 ];
+ }
+ // get the cub_file_set nodes by id
+ // The map is a log search and takes almost no time.
+ // MOAB's linear tag search takes 5-10 minutes.
+ for(unsigned int k=0; k<nodes_per_element; ++k) {
+ std::map<int,MBEntityHandle>::iterator k_iter;
+ k_iter = matched_cub_vert_id_map.find( elem_conn_node_ids[k] );
+
+ if(k_iter == matched_cub_vert_id_map.end()) {
+ std::cout << "ReadNCDF: Found no cub node with id=" << elem_conn_node_ids[k]
+ << ", but expected to find only 1." << std::endl;
+ break;
+ }
+ cub_nodes.insert( k_iter->second );
+ }
+
+ if(nodes_per_element != cub_nodes.size()) {
+ std::cout << "ReadNCDF: nodes_per_elemenet != cub_nodes.size()" << std::endl;
+ return MB_INVALID_SIZE;
+ }
+
+ // get the cub_file_set element with the same nodes
+ int to_dim = MBCN::Dimension(mb_type);
+ rval = mdbImpl->get_adjacencies( cub_nodes, to_dim, false, cub_elem);
+ if(MB_SUCCESS != rval) return rval;
+
+ // Pronto/Presto renumbers elements, so matching cub and exo elements by
+ // id is not possible at this time.
+ } else {
+
+ // get dead element's id
+ int elem_id = elem_ids[first_elem_id_in_block+j];
+ void *id[] = {&elem_id};
+ // get the element by id
+ rval = mdbImpl->get_entities_by_type_and_tag( cub_file_set, mb_type,
+ &mGlobalIdTag, id, 1, cub_elem,
+ MBInterface::INTERSECT );
+ if(MB_SUCCESS != rval) return rval;
+ }
+
+ if(1 == cub_elem.size()) {
+ // Delete the dead element from the cub file. It will be removed from sets
+ // ONLY if they are tracking meshsets.
+ rval = mdbImpl->remove_entities( cub_file_set, cub_elem );
+ if(MB_SUCCESS != rval) return rval;
+ rval = mdbImpl->delete_entities( cub_elem );
+ if(MB_SUCCESS != rval) return rval;
+ } else {
+ std::cout << "ReadNCDF: Should have found 1 element with type="
+ << mb_type << " in cub_file_set, but instead found "
<< cub_elem.size() << std::endl;
+ rval = mdbImpl->list_entities( cub_nodes );
+ ++missing_elem_counter;
return MB_FAILURE;
}
- // Delete the dead element from the cub file. It will be removed from sets
- // ONLY if they are tracking meshsets.
- rval = mdbImpl->remove_entities( cub_file_set, cub_elem );
- if(MB_SUCCESS != rval) return rval;
- rval = mdbImpl->delete_entities( cub_elem );
- if(MB_SUCCESS != rval) return rval;
++dead_elem_counter;
}
}
// Print some statistics
- temp.str("");
- temp << i->blockId;
- std::cout << " Block " << temp.str() << " has " << dead_elem_counter << "/"
- << i->numElements << " dead elements." << std::endl;
+ temp_ss.str("");
+ temp_ss << i->blockId;
+ std::cout << " Block " << temp_ss.str() << " has " << dead_elem_counter << "/"
+ << i->numElements << " dead elements." << std::endl;
+ if(0 != missing_elem_counter) {
+ std::cout << " " << missing_elem_counter
+ << " dead elements in this block were not found in the cub_file_set. " << std::endl;
+ }
+
// advance the pointers into element ids and block_count
first_elem_id_in_block += i->numElements;
++block_count;
Modified: MOAB/trunk/tools/dagmc/cub2h5m.cc
===================================================================
--- MOAB/trunk/tools/dagmc/cub2h5m.cc 2010-02-26 20:57:00 UTC (rev 3558)
+++ MOAB/trunk/tools/dagmc/cub2h5m.cc 2010-03-01 17:02:28 UTC (rev 3559)
@@ -10,6 +10,7 @@
#include "MBSkinner.hpp"
#include "quads_to_tris.hpp"
#include <limits>
+#include <cstdlib>
#define GF_CUBIT_FILE_TYPE "CUBIT"
#define GF_STEP_FILE_TYPE "STEP"
@@ -18,6 +19,114 @@
#define GF_ACIS_BIN_FILE_TYPE "ACIS_SAB"
#define GF_OCC_BREP_FILE_TYPE "OCC"
+// DAGMC cannot build an OBB tree if all of a volume's surfaces have no facets.
+// To prevent this, remove the cgm surface set if the cub surface set exists,
+// but had its faced removed (due to dead elements). Remember that the cgm_file_set
+// is not TRACKING.
+MBErrorCode remove_empty_cgm_surfs_and_vols( MBInterface *MBI,
+ const MBEntityHandle cgm_file_set,
+ const MBTag idTag,
+ const MBTag dimTag,
+ const bool debug ) {
+
+ MBErrorCode result;
+ const int two = 2;
+ const void* const two_val[] = {&two};
+ MBRange cgm_surfs;
+ result = MBI->get_entities_by_type_and_tag(cgm_file_set, MBENTITYSET, &dimTag,
+ two_val, 1, cgm_surfs );
+ if(MB_SUCCESS != result) return result;
+
+ for(MBRange::iterator i=cgm_surfs.begin(); i!=cgm_surfs.end(); ++i) {
+ int n_tris;
+ result = MBI->get_number_entities_by_type( *i, MBTRI, n_tris );
+ if(MB_SUCCESS != result) return result;
+
+ if(0==n_tris) {
+ int surf_id;
+ result = MBI->tag_get_data(idTag, &(*i), 1, &surf_id);
+ assert(MB_SUCCESS == result);
+
+ MBRange parent_vols;
+ result = MBI->get_parent_meshsets( *i, parent_vols );
+ assert(MB_SUCCESS == result);
+ for(MBRange::iterator j=parent_vols.begin(); j!=parent_vols.end(); ++j) {
+ result = MBI->remove_parent_child( *j, *i );
+ assert(MB_SUCCESS == result);
+ }
+ MBRange child_curves;
+ result = MBI->get_child_meshsets( *i, child_curves );
+ assert(MB_SUCCESS == result);
+ for(MBRange::iterator j=child_curves.begin(); j!=child_curves.end(); ++j) {
+ result = MBI->remove_parent_child( *i, *j );
+ assert(MB_SUCCESS == result);
+ }
+ result = MBI->remove_entities( cgm_file_set, &(*i), 1 );
+ assert(MB_SUCCESS == result);
+
+ // Is the set contained anywhere else? If the surface is in a CUBIT group,
+ // such as "unmerged_surfs" it will cause write_mesh to fail. This should
+ // be a MOAB bug.
+ MBRange all_sets;
+ result = MBI->get_entities_by_type( 0, MBENTITYSET, all_sets );
+ assert(MB_SUCCESS == result);
+ for(MBRange::iterator j=all_sets.begin(); j!=all_sets.end(); ++j) {
+ if(MBI->contains_entities( *j, &(*i), 1 )) {
+ result = MBI->remove_entities( *j, &(*i), 1 );
+ assert(MB_SUCCESS == result);
+ }
+ }
+
+ result = MBI->delete_entities( &(*i), 1 );
+ assert(MB_SUCCESS == result);
+ std::cout << " Surface " << surf_id
+ << ": removed because all of its mesh faces have been removed"
+ << std::endl;
+ }
+ }
+
+ // get all volumes
+ const int three = 3;
+ const void* const three_val[] = {&three};
+ MBRange cgm_vols;
+ result = MBI->get_entities_by_type_and_tag(cgm_file_set, MBENTITYSET, &dimTag,
+ three_val, 1, cgm_vols );
+ if(MB_SUCCESS != result) return result;
+
+ for(MBRange::iterator i=cgm_vols.begin(); i!=cgm_vols.end(); ++i) {
+ // get the volume's number of surfaces
+ int n_surfs;
+ result = MBI->num_child_meshsets( *i, &n_surfs );
+ assert(MB_SUCCESS == result);
+
+ // if no surfaces remain, remove the volume
+ if(0==n_surfs) {
+ int vol_id;
+ result = MBI->tag_get_data(idTag, &(*i), 1, &vol_id);
+ assert(MB_SUCCESS == result);
+ // Is the set contained anywhere else? If the surface is in a CUBIT group,
+ // such as "unmerged_surfs" it will cause write_mesh to fail. This should
+ // be a MOAB bug.
+ MBRange all_sets;
+ result = MBI->get_entities_by_type( 0, MBENTITYSET, all_sets );
+ assert(MB_SUCCESS == result);
+ for(MBRange::iterator j=all_sets.begin(); j!=all_sets.end(); ++j) {
+ if(MBI->contains_entities( *j, &(*i), 1 )) {
+ result = MBI->remove_entities( *j, &(*i), 1 );
+ assert(MB_SUCCESS == result);
+ }
+ }
+ result = MBI->delete_entities( &(*i), 1 );
+ assert(MB_SUCCESS == result);
+ std::cout << " Volume " << vol_id
+ << ": removed because all of its surfaces have been removed"
+ << std::endl;
+ }
+ }
+ return MB_SUCCESS;
+}
+
+
// Given parent volume senses, an id, and a set handle, this function creates a
// new surface set with dimension, geometry category, id, and sense tags.
MBErrorCode build_new_surface( MBInterface *MBI,
@@ -30,26 +139,27 @@
MBErrorCode result;
result = MBI->create_meshset( 0, new_surf );
- assert(MB_SUCCESS == result);
+ if(MB_SUCCESS != result) return result;
if(0 != forward_parent_vol) {
result = MBI->add_parent_child( forward_parent_vol, new_surf );
- assert(MB_SUCCESS == result);
+ if(MB_SUCCESS != result) return result;
}
if(0 != reverse_parent_vol) {
result = MBI->add_parent_child( reverse_parent_vol, new_surf );
- assert(MB_SUCCESS == result);
+ if(MB_SUCCESS != result) return result;
}
const int two = 2;
result = MBI->tag_set_data( dimTag, &new_surf, 1, &two );
- assert(MB_SUCCESS == result);
+ if(MB_SUCCESS != result) return result;
result = MBI->tag_set_data( idTag, &new_surf, 1, &new_surf_id );
- assert(MB_SUCCESS == result);
+ if(MB_SUCCESS != result) return result;
const char geom_category[CATEGORY_TAG_SIZE] = {"Surface\0"};
result = MBI->tag_set_data( categoryTag, &new_surf, 1, &geom_category );
- assert(MB_SUCCESS == result);
+ if(MB_SUCCESS != result) return result;
MBEntityHandle vols[2] = { forward_parent_vol, reverse_parent_vol };
result = MBI->tag_set_data( senseTag, &new_surf, 1, vols );
- assert(MB_SUCCESS == result);
+ if(MB_SUCCESS != result) return result;
+
return MB_SUCCESS;
}
@@ -63,17 +173,17 @@
for(MBRange::const_iterator i=faces.begin(); i!=faces.end(); ++i) {
MBRange adj_elem;
result = MBI->get_adjacencies( &(*i), 1, 3, false, adj_elem );
- assert(MB_SUCCESS == result);
- assert(1 == adj_elem.size());
+ if(MB_SUCCESS != result) return result;
+ if(1 != adj_elem.size()) return MB_INVALID_SIZE;
// get center of element
const MBEntityHandle *elem_conn;
int n_nodes;
result = MBI->get_connectivity( adj_elem.front(), elem_conn, n_nodes );
- assert(MB_SUCCESS == result);
+ if(MB_SUCCESS != result) return result;
MBCartVect elem_coords[n_nodes];
result = MBI->get_coords( elem_conn, n_nodes, elem_coords[0].array() );
- assert(MB_SUCCESS == result);
+ if(MB_SUCCESS != result) return result;
MBCartVect elem_center(0.0);
for(int j=0; j<n_nodes; ++j) elem_center += elem_coords[j];
elem_center /= n_nodes;
@@ -81,10 +191,11 @@
// get the center of the face
const MBEntityHandle *face_conn;
result = MBI->get_connectivity( *i, face_conn, n_nodes );
- assert(MB_SUCCESS == result);
+ if(MB_SUCCESS != result) return result;
MBCartVect face_coords[n_nodes];
result = MBI->get_coords( face_conn, n_nodes, face_coords[0].array() );
- assert(MB_SUCCESS == result);
+ if(MB_SUCCESS != result) return result;
+ assert(4 == n_nodes);
MBCartVect face_center(0.0);
for(int j=0; j<n_nodes; ++j) face_center += face_coords[j];
face_center /= n_nodes;
@@ -111,7 +222,7 @@
MBEntityHandle new_face_conn[4] = { face_conn[3], face_conn[2],
face_conn[1], face_conn[0] };
result = MBI->set_connectivity( *i, new_face_conn, 4 );
- assert(MB_SUCCESS == result);
+ if(MB_SUCCESS != result) return result;
}
}
return MB_SUCCESS;
@@ -146,8 +257,7 @@
}
}
-// Use this to get quad faces from hex elems. MBSkinner does not produce the
-// expected result.
+// Use this to get quad faces from hex elems.
MBErrorCode skin_hex_elems(MBInterface *MBI, MBRange elems, const int dim,
MBRange &faces ) {
// get faces of each hex
@@ -160,14 +270,14 @@
for(MBRange::iterator i=elems.begin(); i!=elems.end(); ++i) {
MBRange elem_faces;
result = MBI->get_adjacencies( &(*i), 1, 2, true, elem_faces );
- assert(MB_SUCCESS == result);
- assert(faces_per_elem == elem_faces.size());
+ if(MB_SUCCESS != result) return result;
+ if(faces_per_elem != elem_faces.size()) return MB_INVALID_SIZE;
for(MBRange::iterator j=elem_faces.begin(); j!=elem_faces.end(); ++j) {
const MBEntityHandle *conn;
int n_nodes;
MBErrorCode result = MBI->get_connectivity( *j, conn, n_nodes );
- assert(MB_SUCCESS == result);
- assert(nodes_per_face == n_nodes);
+ if(MB_SUCCESS != result) return result;
+ if(nodes_per_face != n_nodes) return MB_INVALID_SIZE;
// Sort the node handles of the face
for(int k=0; k<nodes_per_face; ++k) f[counter][k] = conn[k];
qsort( &f[counter][0], nodes_per_face, sizeof(MBEntityHandle), handle_compare );
@@ -185,8 +295,8 @@
if(n_faces-1 == i) {
MBRange face_handle;
result = MBI->get_adjacencies( &(f[i][0]), nodes_per_face, 2, false, face_handle );
- assert(MB_SUCCESS == result);
- assert(1 == face_handle.size());
+ if(MB_SUCCESS != result) return result;
+ if(1 != face_handle.size()) return MB_INVALID_SIZE;
faces.insert( face_handle.front() );
// Due to sort, if a duplicate exists it must be next
} else if( f[i][0]==f[i+1][0] && f[i][1]==f[i+1][1] &&
@@ -201,8 +311,8 @@
} else {
MBRange face_handle;
result = MBI->get_adjacencies( &(f[i][0]), nodes_per_face, 2, false, face_handle );
- assert(MB_SUCCESS == result);
- assert(1 == face_handle.size());
+ if(MB_SUCCESS != result) return result;
+ if(1 != face_handle.size()) return MB_INVALID_SIZE;
faces.insert( face_handle.front() );
}
}
@@ -279,7 +389,7 @@
double ratio[n_elems];
for(int i=0; i<n_elems; ++i) ratio[i] = (defo[i]-orig[i])/orig[i];
- plot_histogram( "Predeformed Element Volume", "Num_Elems", "Volume", 10, orig, n_elems );
+ plot_histogram( "Predeformed Element Volume", "Num_Elems", "Volume [cc]", 10, orig, n_elems );
std::string title = "Element Volume Change Ratio at Time Step " + time_step;
plot_histogram( title, "Num_Elems", "Volume Ratio", 10, ratio, n_elems );
@@ -301,11 +411,11 @@
const MBEntityHandle *conn;
int num_vertices;
MBErrorCode result = MBI->get_connectivity( element, conn, num_vertices );
- assert(MB_SUCCESS == result);
+ if(MB_SUCCESS != result) return result;
MBCartVect coords[num_vertices];
result = MBI->get_coords( conn, num_vertices, coords[0].array() );
- assert(MB_SUCCESS == result);
+ if(MB_SUCCESS != result) return result;
switch( type )
{
@@ -352,9 +462,14 @@
/* Calculate the signed volumes beneath the surface (x 6.0). Use the triangle's
cannonical sense. Do not take sense tags into account. Code taken from
- DagMC::measure_volume. */
+ DagMC::measure_volume.
+
+ Special Case: If the surface is planar, and the plane includes the origin,
+ the signed volume will be ~0. If the signed volume is ~0 then offset everything
+ by a random amount and try again. */
MBErrorCode get_signed_volume( MBInterface *MBI,
- const MBEntityHandle surf_set,
+ const MBEntityHandle surf_set,
+ const MBCartVect offset,
double &signed_volume) {
MBErrorCode rval;
MBRange tris;
@@ -366,11 +481,19 @@
MBCartVect coords[3];
for (MBRange::iterator j = tris.begin(); j != tris.end(); ++j) {
rval = MBI->get_connectivity( *j, conn, len, true );
- if (MB_SUCCESS != rval) return rval;
- assert(3 == len);
+ if(MB_SUCCESS != rval) return rval;
+ if(3 != len) return MB_INVALID_SIZE;
rval = MBI->get_coords( conn, 3, coords[0].array() );
- if (MB_SUCCESS != rval) return rval;
-
+ if(MB_SUCCESS != rval) return rval;
+
+ // Apply offset to avoid calculating 0 for cases when the origin is in the
+ // plane of the surface.
+ for(unsigned int k=0; k<3; ++k) {
+ coords[k][0] += offset[0];
+ coords[k][1] += offset[1];
+ coords[k][2] += offset[2];
+ }
+
coords[1] -= coords[0];
coords[2] -= coords[0];
signed_volume += (coords[0] % (coords[1] * coords[2]));
@@ -405,9 +528,10 @@
const void* const tag_vals[] = { &surf_id, &two };
result = MBI->get_entities_by_type_and_tag(cub_file_set, MBENTITYSET, tags,
tag_vals, 2, cub_surf );
- assert(MB_SUCCESS == result);
+ if(MB_SUCCESS != result) return result;
if(1 != cub_surf.size()) {
- std::cout << " surf_id=" << surf_id << " no meshed surface found"
+ std::cout << " Surface " << surf_id
+ << ": no meshed representation found, using CAD representation instead"
<< std::endl;
continue;
}
@@ -415,22 +539,38 @@
// Get tris that represent the quads of the cub surf
MBRange quads;
result = MBI->get_entities_by_type( cub_surf.front(), MBQUAD, quads );
- assert(MB_SUCCESS == result);
+ if(MB_SUCCESS != result) return result;
MBRange cub_tris;
result = make_tris_from_quads( MBI, quads, cub_tris );
// Add the tris to the same surface meshset as the quads are inside.
result = MBI->add_entities( cub_surf.front(), cub_tris );
- assert( MB_SUCCESS == result);
-
- // get the signed volume for each surface representation
+ if(MB_SUCCESS != result) return result;
+
+ // get the signed volume for each surface representation. Keep trying until
+ // The signed volumes are much greater than numerical precision. Planar
+ // surfaces will have a signed volume of zero if the plane goes through the
+ // origin, unless we apply an offset.
+ const int n_attempts = 10;
+ const int max_random = 10;
+ const double min_signed_vol = 0.1;
double cgm_signed_vol, cub_signed_vol;
- result = get_signed_volume( MBI, *i, cgm_signed_vol );
- assert( MB_SUCCESS == result);
- result = get_signed_volume( MBI, cub_surf.front(), cub_signed_vol );
- assert( MB_SUCCESS == result);
- if(debug) std::cout << " surf_id=" << surf_id << " cgm_signed_vol="
- << cgm_signed_vol << " cub_signed_vol=" << cub_signed_vol << std::endl;
+ for(int j=0; j<n_attempts; ++j) {
+ cgm_signed_vol = 0;
+ cub_signed_vol = 0;
+ MBCartVect offset(std::rand()%max_random, std::rand()%max_random, std::rand()%max_random);
+ result = get_signed_volume( MBI, *i, offset, cgm_signed_vol );
+ if(MB_SUCCESS != result) return result;
+ result = get_signed_volume( MBI, cub_surf.front(), offset, cub_signed_vol );
+ if(MB_SUCCESS != result) return result;
+ if(debug) std::cout << " surf_id=" << surf_id << " cgm_signed_vol="
+ << cgm_signed_vol << " cub_signed_vol=" << cub_signed_vol << std::endl;
+ if(fabs(cgm_signed_vol)>min_signed_vol && fabs(cub_signed_vol)>min_signed_vol) break;
+ if(n_attempts == j+1) {
+ std::cout << "error: signed volume could not be calculated unambiguously" << std::cout;
+ return MB_FAILURE;
+ }
+ }
// If the sign is different, reverse the cgm senses so that both
// representations have the same signed volume.
@@ -438,14 +578,14 @@
(cgm_signed_vol>0 && cub_signed_vol<0) ) {
MBEntityHandle cgm_surf_volumes[2], reversed_cgm_surf_volumes[2];
result = MBI->tag_get_data( senseTag, &(*i), 1, cgm_surf_volumes );
- assert(MB_SUCCESS == result);
+ if(MB_SUCCESS != result) return result;
if(MB_SUCCESS != result) return result;
reversed_cgm_surf_volumes[0] = cgm_surf_volumes[1];
reversed_cgm_surf_volumes[1] = cgm_surf_volumes[0];
result = MBI->tag_set_data( senseTag, &(*i), 1, reversed_cgm_surf_volumes );
- assert(MB_SUCCESS == result);
+ if(MB_SUCCESS != result) return result;
}
}
@@ -481,34 +621,38 @@
const void* const tag_vals[] = { &surf_id, &two };
result = MBI->get_entities_by_type_and_tag(cub_file_set, MBENTITYSET, tags,
tag_vals, 2, cub_surf );
- assert(MB_SUCCESS == result);
+ if(MB_SUCCESS != result) return result;
if(1 != cub_surf.size()) {
- std::cout << " Surface " << surf_id << ": no meshed representation found" << std::endl;
+ std::cout << " Surface " << surf_id
+ << ": no meshed representation found, using CAD representation instead"
+ << std::endl;
continue;
}
// Get tris that represent the quads of the cub surf
MBRange quads;
result = MBI->get_entities_by_type( cub_surf.front(), MBQUAD, quads );
- assert(MB_SUCCESS == result);
+ if(MB_SUCCESS != result) return result;
+
MBRange cub_tris;
result = make_tris_from_quads( MBI, quads, cub_tris );
-
+ if(MB_SUCCESS != result) return result;
+
// Remove the tris from the cgm surf. Don't forget to remove them from the
// cgm_file_set because it is not TRACKING.
MBRange cgm_tris;
result = MBI->get_entities_by_type( *i, MBTRI, cgm_tris );
- assert( MB_SUCCESS == result);
+ if(MB_SUCCESS != result) return result;
result = MBI->remove_entities( *i, cgm_tris );
- assert( MB_SUCCESS == result);
+ if(MB_SUCCESS != result) return result;
result = MBI->remove_entities( cgm_file_set, cgm_tris );
- assert( MB_SUCCESS == result);
+ if(MB_SUCCESS != result) return result;
result = MBI->delete_entities( cgm_tris );
- assert( MB_SUCCESS == result);
+ if(MB_SUCCESS != result) return result;
// Add the cub_tris to the cgm_surf
result = MBI->add_entities( *i, cub_tris );
- assert(MB_SUCCESS == result);
+ if(MB_SUCCESS != result) return result;
}
return MB_SUCCESS;
@@ -536,13 +680,13 @@
two_val, 1, cgm_surfs );
if(MB_SUCCESS != result) return result;
- // Get the maximum surface id. This is so that new surfaces to do have
+ // Get the maximum surface id. This is so that new surfaces do not have
// duplicate ids.
int max_surf_id = -1;
for(MBRange::const_iterator i=cgm_surfs.begin(); i!=cgm_surfs.end(); ++i) {
int surf_id;
result = MBI->tag_get_data( idTag, &(*i), 1, &surf_id );
- assert(MB_SUCCESS == result);
+ if(MB_SUCCESS != result) return result;
if(max_surf_id < surf_id) max_surf_id = surf_id;
}
std::cout << " Maximum surface id=" << max_surf_id << std::endl;
@@ -553,13 +697,15 @@
MBRange cgm_vols;
result = MBI->get_entities_by_type_and_tag(cgm_file_set, MBENTITYSET, &dimTag,
three_val, 1, cgm_vols );
- assert(MB_SUCCESS == result);
+ if(MB_SUCCESS != result) return result;
+
// get the corresponding cub volume
for(MBRange::iterator i=cgm_vols.begin(); i!=cgm_vols.end(); i++) {
int vol_id;
result = MBI->tag_get_data(idTag, &(*i), 1, &vol_id);
+ assert(MB_SUCCESS == result);
if(MB_SUCCESS != result) return result;
- std::cout << " Volume " << vol_id << std::endl;
+ std::cout << " Volume " << vol_id;
// Find the meshed vol set with the same id
MBRange cub_vol;
@@ -568,68 +714,94 @@
result = MBI->get_entities_by_type_and_tag(cub_file_set, MBENTITYSET, tags,
tag_vals, 2, cub_vol );
assert(MB_SUCCESS == result);
+ if(MB_SUCCESS != result) return result;
if(1 != cub_vol.size()) {
- std::cout << " no meshed volume found" << std::endl;
+ std::cout << ": no meshed representation found" << std::endl;
continue;
+ } else {
+ std::cout << std::endl;
}
// get the mesh elements of the volume.
MBRange elems;
result = MBI->get_entities_by_type( cub_vol.front(), MBHEX, elems );
- assert(MB_SUCCESS == result);
+ assert(MB_SUCCESS == result);
+ if(MB_SUCCESS != result) return result;
if(debug) std::cout << " found " << elems.size() << " hex elems" << std::endl;
// skin the volumes
- // BEWARE: THE MBSKINNER DOES NOT PRODUCE THE EXPECTED RESULT. CONFIRMED IN VISIT.
MBSkinner tool(MBI);
- MBRange moab_faces, my_faces;
- result = tool.find_skin( elems, 2, moab_faces, true );
+ MBRange skin_faces;
+ result = tool.find_skin( elems, 2, skin_faces, true );
assert(MB_SUCCESS == result);
- skin_hex_elems( MBI, elems, 2, my_faces );
- assert(MB_SUCCESS == result);
- if(debug) std::cout << " moab_faces=" << moab_faces.size() << " my_faces="
- << my_faces.size() << std::endl;
+ if(MB_SUCCESS != result) return result;
+ // Reconcile the difference between faces of the cub file surfaces and skin
+ // faces of the 3D mesh. The difference exists because dead elements have been
+ // removed. Faces are divided into:
+ // cub_faces (in) - the faces in the cub file surface
+ // skin_faces (in) - the faces of the 3D mesh elements
+ // common_faces (out)- the faces common to both the cub file surface and 3D mesh
+ // they are still adjacent to this volume
+ // old_faces (out)- the faces of the cub surface not on the 3D mesh skin
+ // they are no longer adjacent to this vol
+ // new_faces (out)- the faces of the 3D mesh skin not on the cub surface
+ // they are now adjacent to this volume
+
// get cub child surfaces.
MBRange cub_surfs;
result = MBI->get_child_meshsets( cub_vol.front(), cub_surfs );
assert(MB_SUCCESS == result);
+ if(MB_SUCCESS != result) return result;
for(MBRange::iterator j=cub_surfs.begin(); j!=cub_surfs.end(); ++j) {
+ int surf_id;
+ result = MBI->tag_get_data( idTag, &(*j), 1, &surf_id );
+ assert(MB_SUCCESS == result);
+ if(MB_SUCCESS != result) return result;
// get the quads on each surface
MBRange cub_faces;
result = MBI->get_entities_by_type( *j, MBQUAD, cub_faces );
assert(MB_SUCCESS == result);
+ if(MB_SUCCESS != result) return result;
// Meshed volumes must have meshed surfaces
- assert(!cub_faces.empty());
+ if(cub_faces.empty()) {
+ std::cout << " Surface " << surf_id << ": contains no meshed faces" << std::endl;
+ // return MB_ENTITY_NOT_FOUND;
+ }
// get the faces common to both the skin and this surface
- MBRange common_faces = intersect( cub_faces, my_faces );
- // find the surfaces faces not on the skin - these are orphaned and need removed
- MBRange orphaned_faces = subtract( cub_faces, common_faces );
- result = MBI->remove_entities( *j, orphaned_faces );
+ MBRange common_faces = intersect( cub_faces, skin_faces );
+ // find the surface faces not on the skin - these are old and need removed
+ MBRange old_faces = subtract( cub_faces, common_faces );
+ result = MBI->remove_entities( *j, old_faces );
assert(MB_SUCCESS == result);
- int surf_id;
- result = MBI->tag_get_data( idTag, &(*j), 1, &surf_id );
- assert(MB_SUCCESS == result);
+ if(MB_SUCCESS != result) return result;
+
// remove the common faces from the skin faces
- my_faces = subtract( my_faces, common_faces );
- // If no orphaned faces exist we are done
- if(orphaned_faces.empty()) continue;
- std::cout << " Surface " << surf_id << " had " << orphaned_faces.size()
- << " orphaned faces removed" << std::endl;
- // place the orphaned faces in a new surface. Get the parent vols of the
- // surf.
+ skin_faces = subtract( skin_faces, common_faces );
+ // If no old faces exist we are done
+ if(old_faces.empty()) continue;
+ std::cout << " Surface " << surf_id << ": " << old_faces.size()
+ << " old faces removed" << std::endl;
+ // Place the old faces in a new surface, because they may still be adjacent
+ // to 3D mesh in another volume. Get the parent vols of the surface.
MBRange cgm_surf;
const MBTag tags[] = {idTag, dimTag};
const void* const tag_vals[] = { &surf_id, &two };
result = MBI->get_entities_by_type_and_tag(cgm_file_set, MBENTITYSET, tags,
tag_vals, 2, cgm_surf );
assert(MB_SUCCESS == result);
- assert(1 == cgm_surf.size());
+ if(MB_SUCCESS != result) return result;
+ if(1 != cgm_surf.size()) {
+ std::cout << "invalid size" << std::endl;
+ return MB_INVALID_SIZE;
+ }
MBEntityHandle cgm_vols[2], cub_vols[2];
result = MBI->tag_get_data( senseTag, cgm_surf, &cgm_vols );
assert(MB_SUCCESS == result);
+ if(MB_SUCCESS != result) return result;
result = MBI->tag_get_data( senseTag, &(*j), 1, &cub_vols );
assert(MB_SUCCESS == result);
+ if(MB_SUCCESS != result) return result;
// for the new surf, replace the current volume with the impl compl vol.
// This is because the faces that no longer exist will become adjacent to
// the impl compl
@@ -643,10 +815,10 @@
}
// If both sides of the surface are the impl comp, do not create the surface.
if(0==cgm_vols[0] && 0==cgm_vols[1]) {
- std::cout << " New surface was not created for orphaned faces because both parents are impl_compl volume " << std::endl;
+ std::cout << " New surface was not created for old faces because both parents are impl_compl volume " << std::endl;
continue;
}
- // build the new surface, convert quads to tris, and add the faces.
+ // build the new surface.
MBEntityHandle new_cgm_surf, new_cub_surf;
++max_surf_id;
result = build_new_surface( MBI, new_cgm_surf,
@@ -654,30 +826,41 @@
dimTag, idTag,
categoryTag, senseTag );
assert(MB_SUCCESS == result);
+ if(MB_SUCCESS != result) return result;
result = build_new_surface( MBI, new_cub_surf,
cub_vols[0], cub_vols[1], max_surf_id,
dimTag, idTag,
categoryTag, senseTag );
assert(MB_SUCCESS == result);
- // add the new surface to the file set and populate it with faces
+ if(MB_SUCCESS != result) return result;
+ // add the new surface to the file set and populate it with the old faces
result = MBI->add_entities( cgm_file_set, &new_cgm_surf, 1 );
assert(MB_SUCCESS == result);
+ if(MB_SUCCESS != result) return result;
+ assert(MB_SUCCESS == result);
result = MBI->add_entities( cub_file_set, &new_cub_surf, 1 );
- assert(MB_SUCCESS == result);
- result = MBI->add_entities( new_cub_surf, orphaned_faces );
+ if(MB_SUCCESS != result) return result;
assert(MB_SUCCESS == result);
- std::cout << " Surface " << max_surf_id << " was created for "
- << orphaned_faces.size() << " orphaned faces" << std::endl;
+ result = MBI->add_entities( new_cub_surf, old_faces );
+ assert(MB_SUCCESS == result);
+ if(MB_SUCCESS != result) return result;
+
+ std::cout << " Surface " << max_surf_id << ": created for "
+ << old_faces.size() << " old faces" << std::endl;
}
- // remaining skin faces must be assigned to a surface
- if(my_faces.empty()) continue;
- std::cout << " Surface " << max_surf_id+1 << " will be created for "
- << my_faces.size() << " dead skin faces" << std::endl;
+ // the remaining skin faces are newly exposed faces
+ MBRange new_faces = skin_faces;
+ // new skin faces must be assigned to a surface
+ if(new_faces.empty()) continue;
+ std::cout << " Surface " << max_surf_id+1 << ": created for "
+ << new_faces.size() << " new faces" << std::endl;
+
// Ensure that faces are oriented outwards
- result = orient_faces_outward( MBI, my_faces, debug );
+ result = orient_faces_outward( MBI, new_faces, debug );
assert(MB_SUCCESS == result);
+ if(MB_SUCCESS != result) return result;
// Create the new surface.
MBEntityHandle new_cgm_surf, new_cub_surf;
@@ -687,18 +870,23 @@
dimTag, idTag,
categoryTag, senseTag );
assert(MB_SUCCESS == result);
+ if(MB_SUCCESS != result) return result;
result = build_new_surface( MBI, new_cub_surf,
cub_vol.front(), 0, max_surf_id,
dimTag, idTag,
categoryTag, senseTag );
assert(MB_SUCCESS == result);
+ if(MB_SUCCESS != result) return result;
// Insert the new surf into file sets and populate it with faces.
result = MBI->add_entities( cgm_file_set, &new_cgm_surf, 1 );
assert(MB_SUCCESS == result);
+ if(MB_SUCCESS != result) return result;
result = MBI->add_entities( cub_file_set, &new_cub_surf, 1 );
- assert(MB_SUCCESS == result);
- result = MBI->add_entities( new_cub_surf, my_faces );
assert(MB_SUCCESS == result);
+ if(MB_SUCCESS != result) return result;
+ result = MBI->add_entities( new_cub_surf, new_faces );
+ assert(MB_SUCCESS == result);
+ if(MB_SUCCESS != result) return result;
}
return MB_SUCCESS;
@@ -729,7 +917,6 @@
// get signed volume from cgm and cub surf MUST COME BEFORE COORD UPDATE, NEEDS TRIS
// reverse cgm surface sense if needed
// replace cgm surf tris with cub surf tris
-// X remove quads from cub surf
// measure volume of predeformed cub elements
// convert cub volumes sets to tracking so that dead elems are removed from vol sets
// update coordinates and delete dead elems
@@ -738,12 +925,13 @@
// for each cub volume
// skin volume elems to get faces
// for each child cub surface
-// remove surface faces that are not in skin
+// assign old skin faces to a new surface in case they are adjacent to another volume
// orient each skin face outward
-// assign new (leftover) skin faces to a surface
+// assign new skin faces to a new surface
// for each surface
// remove existing tris (from before the update)
// convert quads to tris
+// remove empty surfaces and volumes due to dead elements
int main( int argc, char* argv[] )
{
const bool debug = false;
@@ -822,7 +1010,7 @@
MBErrorCode result;
MBEntityHandle cub_file_set;
result = MBI->create_meshset( 0, cub_file_set );
- assert(MB_SUCCESS == result);
+ if(MB_SUCCESS != result) return result;
//char cub_options[256] = "120";
//result = MBI->load_file(cub_name, cub_file_set, cub_options, NULL, 0, 0);
result = MBI->load_file(cub_name, &cub_file_set, 0, NULL, 0, 0);
@@ -836,10 +1024,10 @@
dist_tol,norm_tol,len_tol);
MBEntityHandle cgm_file_set;
result = MBI->create_meshset( 0, cgm_file_set );
- assert(MB_SUCCESS == result);
+ if(MB_SUCCESS != result) return result;
result = MBI->load_file(sat_name, &cgm_file_set, cgm_options,NULL,0,0);
if(MB_SUCCESS != result) return result;
- std::cout << "Geometry file read." << std::endl;
+ std::cout << "CAD file read." << std::endl;
// Create tags
MBTag dimTag, idTag, categoryTag, senseTag;
@@ -861,15 +1049,15 @@
// cub representations. Change the sense of the cgm representation to match
// the cub representation.
result = fix_surface_senses( MBI, cgm_file_set, cub_file_set, idTag, dimTag, senseTag, debug );
- assert(MB_SUCCESS == result);
- std::cout << "Fixed geometry surface senses to match meshed surface senses." << std::endl;
+ if(MB_SUCCESS != result) return result;
+ std::cout << "Fixed CAD surface senses to match meshed surface senses." << std::endl;
// Get the 3D elements in the cub file and measure their volume.
MBRange orig_elems;
std::vector<double> orig_size;
if(determine_volume_change) {
result = MBI->get_entities_by_dimension( 0, 3, orig_elems );
- assert(MB_SUCCESS == result);
+ if(MB_SUCCESS != result) return result;
orig_size.resize(orig_elems.size());
for(unsigned int i=0; i<orig_elems.size(); ++i) {
orig_size[i] = measure( MBI, orig_elems[i] );
@@ -883,10 +1071,10 @@
MBRange cub_vols;
result = MBI->get_entities_by_type_and_tag(cub_file_set, MBENTITYSET, &dimTag,
three_val, 1, cub_vols );
- assert(MB_SUCCESS == result);
+ if(MB_SUCCESS != result) return result;
for(MBRange::const_iterator i=cub_vols.begin(); i!=cub_vols.end(); ++i) {
result = MBI->set_meshset_options( *i, MESHSET_TRACK_OWNER );
- assert(MB_SUCCESS == result);
+ if(MB_SUCCESS != result) return result;
}
// Update the coordinates if needed. Do not do this before checking surface
@@ -895,7 +1083,7 @@
// The cub node ids are unique because cgm vertex ids are tagged on the vertex
// meshset and not the vertex itself.
//result = MBI->delete_entities( &cub_file_set, 1 );
- //assert(MB_SUCCESS == result);
+ //if(MB_SUCCESS != result) return result;
// Assume dead elements exist until I think of something better.
bool dead_elements_exist = true;
if(update_coords) {
@@ -909,7 +1097,9 @@
//result = my_ex_reader.load_file(exo_name, cub_file_set, exo_opts, NULL, 0 , 0);
result = my_ex_reader.load_file(exo_name, &cub_file_set, exo_opts, NULL, 0 , 0);
if(MB_SUCCESS != result) {
- std::cout << "coordinate update failed" << std::endl;
+ std::string last_error;
+ MBI->get_last_error(last_error);
+ std::cout << "coordinate update failed, " << last_error << std::endl;
return result;
}
std::cout << "Updated mesh nodes with deformed coordinates from exodus file." << std::endl;
@@ -920,7 +1110,7 @@
// still exist.
MBRange defo_elems;
result = MBI->get_entities_by_dimension( 0, 3, defo_elems );
- assert(MB_SUCCESS == result);
+ if(MB_SUCCESS != result) return result;
// Determine the volume of the elements now that a deformation has been
// applied. Condense the original array by removing dead elements.
@@ -943,8 +1133,9 @@
if(update_coords && dead_elements_exist) {
result = add_dead_elems_to_impl_compl( MBI, cgm_file_set, cub_file_set,
idTag, dimTag, categoryTag, senseTag, debug );
- assert(MB_SUCCESS == result);
- std::cout << "Placed dead elements to implicit complement volume and added required surfaces." << std::endl;
+ if(MB_SUCCESS != result) return result;
+ std::cout << "Placed dead elements in implicit complement volume and added required surfaces."
+ << std::endl;
}
// The quads in the cub_file_set have been updated for dead elements. For each
@@ -952,9 +1143,15 @@
// with cub_tris (created from the quads). Note the a surface that is not
// meshed (in cub file) will not be effected.
result = replace_faceted_cgm_surfs( MBI, cgm_file_set, cub_file_set, idTag, dimTag, debug );
- assert(MB_SUCCESS == result);
- std::cout << "Replaced faceted geometry surfaces with meshed surfaces of triangles." << std::endl;
+ if(MB_SUCCESS != result) return result;
+ std::cout << "Replaced faceted CAD surfaces with meshed surfaces of triangles."
+ << std::endl;
+ result = remove_empty_cgm_surfs_and_vols( MBI, cgm_file_set, idTag, dimTag, debug );
+ if(MB_SUCCESS != result) return result;
+ std::cout << "Removed surfaces and volumes that no longer have any triangles."
+ << std::endl;
+
result = MBI->write_mesh( out_name, &cgm_file_set, 1 );
if(MB_SUCCESS != result) {
std::cout << "write mesh failed" << std::endl;
More information about the moab-dev
mailing list