MOAB: Mesh Oriented datABase  (version 5.2.1)
cub2h5m.cpp
Go to the documentation of this file.
00001 #include "moab/GeomQueryTool.hpp"
00002 #include "moab/GeomTopoTool.hpp"
00003 #include "InitCGMA.hpp"
00004 #include "CGMApp.hpp"
00005 #include "moab/CN.hpp"
00006 #include "moab/Core.hpp"
00007 #include "moab/CartVect.hpp"
00008 #include "moab/FileOptions.hpp"
00009 #include "moab/Skinner.hpp"
00010 #include "quads_to_tris.hpp"
00011 #include <limits>
00012 #include <cstdlib>
00013 #include <sstream>
00014 #include <time.h>
00015 
00016 #define GF_CUBIT_FILE_TYPE    "CUBIT"
00017 #define GF_STEP_FILE_TYPE     "STEP"
00018 #define GF_IGES_FILE_TYPE     "IGES"
00019 #define GF_ACIS_TXT_FILE_TYPE "ACIS_SAT"
00020 #define GF_ACIS_BIN_FILE_TYPE "ACIS_SAB"
00021 #define GF_OCC_BREP_FILE_TYPE "OCC"
00022 
00023 using namespace moab;
00024 
00025 void tokenize( const std::string& str, std::vector< std::string >& tokens, const char* delimiters )
00026 {
00027     std::string::size_type last = str.find_first_not_of( delimiters, 0 );
00028     std::string::size_type pos  = str.find_first_of( delimiters, last );
00029     if( std::string::npos == pos )
00030         tokens.push_back( str );
00031     else
00032         while( std::string::npos != pos && std::string::npos != last )
00033         {
00034             tokens.push_back( str.substr( last, pos - last ) );
00035             last = str.find_first_not_of( delimiters, pos );
00036             pos  = str.find_first_of( delimiters, last );
00037             if( std::string::npos == pos ) pos = str.size();
00038         }
00039 }
00040 
00041 ErrorCode get_group_names( Interface* MBI, const EntityHandle group_set, const Tag nameTag,
00042                            std::vector< std::string >& grp_names )
00043 {
00044     // get names
00045     char name0[NAME_TAG_SIZE];
00046     std::fill( name0, name0 + NAME_TAG_SIZE, '\0' );
00047     ErrorCode result = MBI->tag_get_data( nameTag, &group_set, 1, &name0 );
00048     if( MB_SUCCESS != result && MB_TAG_NOT_FOUND != result ) return MB_FAILURE;
00049 
00050     if( MB_TAG_NOT_FOUND != result ) grp_names.push_back( std::string( name0 ) );
00051 
00052     return MB_SUCCESS;
00053 }
00054 
00055 // For each material, sum the volume. If the coordinates were updated for
00056 // deformation, summarize the volume change.
00057 ErrorCode summarize_cell_volume_change( Interface* MBI, const EntityHandle cgm_file_set, const Tag categoryTag,
00058                                         const Tag dimTag, const Tag sizeTag, const Tag nameTag, const Tag idTag,
00059                                         const bool conserve_mass, const bool debug )
00060 {
00061     // get groups
00062     ErrorCode rval;
00063     const char group_category[CATEGORY_TAG_SIZE] = { "Group\0" };
00064     const void* const group_val[]                = { &group_category };
00065     Range groups;
00066     rval = MBI->get_entities_by_type_and_tag( 0, MBENTITYSET, &categoryTag, group_val, 1, groups );
00067     if( MB_SUCCESS != rval ) return rval;
00068 
00069     // Get the maximum group id. This is so that new groups do not have
00070     // duplicate ids.
00071     int max_grp_id = -1;
00072     for( Range::const_iterator i = groups.begin(); i != groups.end(); ++i )
00073     {
00074         int grp_id;
00075         rval = MBI->tag_get_data( idTag, &( *i ), 1, &grp_id );
00076         if( MB_SUCCESS != rval ) return rval;
00077         if( max_grp_id < grp_id ) max_grp_id = grp_id;
00078     }
00079     if( conserve_mass )
00080     {
00081         std::cout << "  Altering group densities to conserve mass for each volume." << std::endl;
00082         std::cout << "    maximum group id=" << max_grp_id << std::endl;
00083     }
00084 
00085     for( Range::const_iterator i = groups.begin(); i != groups.end(); ++i )
00086     {
00087         // get group names
00088         std::vector< std::string > grp_names;
00089         rval = get_group_names( MBI, *i, nameTag, grp_names );
00090         if( MB_SUCCESS != rval ) return MB_FAILURE;
00091 
00092         // determine if it is a material group
00093         bool material_grp = false;
00094         int mat_id        = -1;
00095         double rho        = 0;
00096         for( std::vector< std::string >::const_iterator j = grp_names.begin(); j != grp_names.end(); ++j )
00097         {
00098             if( std::string::npos != ( *j ).find( "mat" ) && std::string::npos != ( *j ).find( "rho" ) )
00099             {
00100                 material_grp = true;
00101                 std::cout << "  material group: " << *j << std::endl;
00102 
00103                 // get the density and material id
00104                 std::vector< std::string > tokens;
00105                 tokenize( *j, tokens, "_" );
00106                 mat_id = atoi( tokens[1].c_str() );
00107                 rho    = atof( tokens[3].c_str() );
00108             }
00109         }
00110         if( !material_grp ) continue;
00111 
00112         // get the volume sets of the material group
00113         const int three               = 3;
00114         const void* const three_val[] = { &three };
00115         Range vols;
00116         rval = MBI->get_entities_by_type_and_tag( *i, MBENTITYSET, &dimTag, three_val, 1, vols );
00117         if( MB_SUCCESS != rval ) return rval;
00118 
00119         // for each volume, sum predeformed and deformed volume
00120         double orig_grp_volume = 0, defo_grp_volume = 0;
00121 
00122         moab::GeomTopoTool gtt  = moab::GeomTopoTool( MBI, false );
00123         moab::GeomQueryTool gqt = moab::GeomQueryTool( &gtt );
00124         for( Range::const_iterator j = vols.begin(); j != vols.end(); ++j )
00125         {
00126             double defo_size = 0, orig_size = 0;
00127             rval = gqt.measure_volume( *j, defo_size );
00128             if( MB_SUCCESS != rval ) return rval;
00129             defo_grp_volume += defo_size;
00130             rval = MBI->tag_get_data( sizeTag, &( *j ), 1, &orig_size );
00131             if( MB_SUCCESS != rval ) return rval;
00132             orig_grp_volume += orig_size;
00133 
00134             // calculate a new density to conserve mass through the deformation
00135             if( !conserve_mass ) continue;
00136             double new_rho = rho * orig_size / defo_size;
00137 
00138             // create a group for the volume with modified density
00139             EntityHandle new_grp;
00140             rval = MBI->create_meshset( MESHSET_SET, new_grp );
00141             if( MB_SUCCESS != rval ) return rval;
00142             std::stringstream new_name_ss;
00143             new_name_ss << "mat_" << mat_id << "_rho_" << new_rho << "\0";
00144             std::string new_name;
00145             new_name_ss >> new_name;
00146             rval = MBI->tag_set_data( nameTag, &new_grp, 1, new_name.c_str() );
00147             if( MB_SUCCESS != rval ) return rval;
00148             max_grp_id++;
00149             rval = MBI->tag_set_data( idTag, &new_grp, 1, &max_grp_id );
00150             if( MB_SUCCESS != rval ) return rval;
00151             const char group_category2[CATEGORY_TAG_SIZE] = "Group\0";
00152             rval = MBI->tag_set_data( categoryTag, &new_grp, 1, group_category2 );
00153             if( MB_SUCCESS != rval ) return rval;
00154 
00155             // add the volume to the new group
00156             rval = MBI->add_entities( new_grp, &( *j ), 1 );
00157             if( MB_SUCCESS != rval ) return rval;
00158 
00159             // add the new grp to the cgm_file_set
00160             rval = MBI->add_entities( cgm_file_set, &new_grp, 1 );
00161             if( MB_SUCCESS != rval ) return rval;
00162 
00163             // remove the volume from the old group
00164             rval = MBI->remove_entities( *i, &( *j ), 1 );
00165             if( MB_SUCCESS != rval ) return rval;
00166             if( debug ) std::cout << "    new group: " << new_name << " id=" << max_grp_id << std::endl;
00167         }
00168 
00169         std::cout << "    orig_volume=" << orig_grp_volume << " defo_volume=" << defo_grp_volume
00170                   << " defo/orig=" << defo_grp_volume / orig_grp_volume << std::endl;
00171     }
00172 
00173     return MB_SUCCESS;
00174 }
00175 
00176 // We cannot build an OBB tree if all of a volume's surfaces have no facets.
00177 // To prevent this, remove the cgm surface set if the cub surface set exists,
00178 // but had its faced removed (due to dead elements). Remember that the cgm_file_set
00179 // is not TRACKING.
00180 ErrorCode remove_empty_cgm_surfs_and_vols( Interface* MBI, const EntityHandle cgm_file_set, const Tag idTag,
00181                                            const Tag dimTag, const bool /*debug */ )
00182 {
00183 
00184     ErrorCode result;
00185     const int two               = 2;
00186     const void* const two_val[] = { &two };
00187     Range cgm_surfs;
00188     result = MBI->get_entities_by_type_and_tag( cgm_file_set, MBENTITYSET, &dimTag, two_val, 1, cgm_surfs );
00189     if( MB_SUCCESS != result ) return result;
00190 
00191     for( Range::iterator i = cgm_surfs.begin(); i != cgm_surfs.end(); ++i )
00192     {
00193         int n_tris;
00194         result = MBI->get_number_entities_by_type( *i, MBTRI, n_tris );
00195         if( MB_SUCCESS != result ) return result;
00196 
00197         if( 0 == n_tris )
00198         {
00199             int surf_id;
00200             result = MBI->tag_get_data( idTag, &( *i ), 1, &surf_id );
00201             assert( MB_SUCCESS == result );
00202 
00203             Range parent_vols;
00204             result = MBI->get_parent_meshsets( *i, parent_vols );
00205             assert( MB_SUCCESS == result );
00206             for( Range::iterator j = parent_vols.begin(); j != parent_vols.end(); ++j )
00207             {
00208                 result = MBI->remove_parent_child( *j, *i );
00209                 assert( MB_SUCCESS == result );
00210             }
00211             Range child_curves;
00212             result = MBI->get_child_meshsets( *i, child_curves );
00213             assert( MB_SUCCESS == result );
00214             for( Range::iterator j = child_curves.begin(); j != child_curves.end(); ++j )
00215             {
00216                 result = MBI->remove_parent_child( *i, *j );
00217                 assert( MB_SUCCESS == result );
00218             }
00219             result = MBI->remove_entities( cgm_file_set, &( *i ), 1 );
00220             assert( MB_SUCCESS == result );
00221 
00222             // Is the set contained anywhere else? If the surface is in a CUBIT group,
00223             // such as "unmerged_surfs" it will cause write_mesh to fail. This should
00224             // be a MOAB bug.
00225             Range all_sets;
00226             result = MBI->get_entities_by_type( 0, MBENTITYSET, all_sets );
00227             assert( MB_SUCCESS == result );
00228             for( Range::iterator j = all_sets.begin(); j != all_sets.end(); ++j )
00229             {
00230                 if( MBI->contains_entities( *j, &( *i ), 1 ) )
00231                 {
00232                     result = MBI->remove_entities( *j, &( *i ), 1 );
00233                     assert( MB_SUCCESS == result );
00234                 }
00235             }
00236 
00237             result = MBI->delete_entities( &( *i ), 1 );
00238             assert( MB_SUCCESS == result );
00239             std::cout << "  Surface " << surf_id << ": removed because all of its mesh faces have been removed"
00240                       << std::endl;
00241         }
00242     }
00243 
00244     // get all volumes
00245     const int three               = 3;
00246     const void* const three_val[] = { &three };
00247     Range cgm_vols;
00248     result = MBI->get_entities_by_type_and_tag( cgm_file_set, MBENTITYSET, &dimTag, three_val, 1, cgm_vols );
00249     if( MB_SUCCESS != result ) return result;
00250 
00251     for( Range::iterator i = cgm_vols.begin(); i != cgm_vols.end(); ++i )
00252     {
00253         // get the volume's number of surfaces
00254         int n_surfs;
00255         result = MBI->num_child_meshsets( *i, &n_surfs );
00256         assert( MB_SUCCESS == result );
00257 
00258         // if no surfaces remain, remove the volume
00259         if( 0 == n_surfs )
00260         {
00261             int vol_id;
00262             result = MBI->tag_get_data( idTag, &( *i ), 1, &vol_id );
00263             assert( MB_SUCCESS == result );
00264             // Is the set contained anywhere else? If the surface is in a CUBIT group,
00265             // such as "unmerged_surfs" it will cause write_mesh to fail. This should
00266             // be a MOAB bug.
00267             Range all_sets;
00268             result = MBI->get_entities_by_type( 0, MBENTITYSET, all_sets );
00269             assert( MB_SUCCESS == result );
00270             for( Range::iterator j = all_sets.begin(); j != all_sets.end(); ++j )
00271             {
00272                 if( MBI->contains_entities( *j, &( *i ), 1 ) )
00273                 {
00274                     result = MBI->remove_entities( *j, &( *i ), 1 );
00275                     assert( MB_SUCCESS == result );
00276                 }
00277             }
00278             result = MBI->delete_entities( &( *i ), 1 );
00279             assert( MB_SUCCESS == result );
00280             std::cout << "  Volume " << vol_id << ": removed because all of its surfaces have been removed"
00281                       << std::endl;
00282         }
00283     }
00284     return MB_SUCCESS;
00285 }
00286 
00287 // Given parent volume senses, an id, and a set handle, this function creates a
00288 // new surface set with dimension, geometry category, id, and sense tags.
00289 ErrorCode build_new_surface( Interface* MBI, EntityHandle& new_surf, const EntityHandle forward_parent_vol,
00290                              const EntityHandle reverse_parent_vol, const int new_surf_id, const Tag dimTag,
00291                              const Tag idTag, const Tag categoryTag, const Tag senseTag )
00292 {
00293 
00294     ErrorCode result;
00295     result = MBI->create_meshset( 0, new_surf );
00296     if( MB_SUCCESS != result ) return result;
00297     if( 0 != forward_parent_vol )
00298     {
00299         result = MBI->add_parent_child( forward_parent_vol, new_surf );
00300         if( MB_SUCCESS != result ) return result;
00301     }
00302     if( 0 != reverse_parent_vol )
00303     {
00304         result = MBI->add_parent_child( reverse_parent_vol, new_surf );
00305         if( MB_SUCCESS != result ) return result;
00306     }
00307     const int two = 2;
00308     result        = MBI->tag_set_data( dimTag, &new_surf, 1, &two );
00309     if( MB_SUCCESS != result ) return result;
00310     result = MBI->tag_set_data( idTag, &new_surf, 1, &new_surf_id );
00311     if( MB_SUCCESS != result ) return result;
00312     const char geom_category[CATEGORY_TAG_SIZE] = { "Surface\0" };
00313     result                                      = MBI->tag_set_data( categoryTag, &new_surf, 1, &geom_category );
00314     if( MB_SUCCESS != result ) return result;
00315     EntityHandle vols[2] = { forward_parent_vol, reverse_parent_vol };
00316     result               = MBI->tag_set_data( senseTag, &new_surf, 1, vols );
00317     if( MB_SUCCESS != result ) return result;
00318 
00319     return MB_SUCCESS;
00320 }
00321 
00322 // Given a face, orient it outward wrt its adjacent mesh element.
00323 // Each face must be adjacent to exactly one mesh element.
00324 ErrorCode orient_faces_outward( Interface* MBI, const Range faces, const bool /*debug*/ )
00325 {
00326 
00327     ErrorCode result;
00328     for( Range::const_iterator i = faces.begin(); i != faces.end(); ++i )
00329     {
00330         Range adj_elem;
00331         result = MBI->get_adjacencies( &( *i ), 1, 3, false, adj_elem );
00332         if( MB_SUCCESS != result ) return result;
00333         if( 1 != adj_elem.size() ) return MB_INVALID_SIZE;
00334 
00335         // get connectivity for element and face
00336         const EntityHandle* elem_conn;
00337         int elem_n_nodes;
00338         result = MBI->get_connectivity( adj_elem.front(), elem_conn, elem_n_nodes );
00339         if( MB_SUCCESS != result ) return result;
00340         const EntityHandle* face_conn;
00341         int face_n_nodes;
00342         result = MBI->get_connectivity( *i, face_conn, face_n_nodes );
00343         if( MB_SUCCESS != result ) return result;
00344 
00345         // Get the sense of the face wrt the element
00346         EntityType elem_type = MBI->type_from_handle( adj_elem.front() );
00347         EntityType face_type = MBI->type_from_handle( *i );
00348         int face_num, offset;
00349         int sense          = 0;
00350         const int face_dim = CN::Dimension( face_type );
00351         int rval = CN::SideNumber( elem_type, elem_conn, face_conn, face_n_nodes, face_dim, face_num, sense, offset );
00352         if( 0 != rval ) return MB_FAILURE;
00353 
00354         // If the face is not oriented outward wrt the element, reverse it
00355         if( -1 == sense )
00356         {
00357             EntityHandle new_face_conn[4] = { face_conn[3], face_conn[2], face_conn[1], face_conn[0] };
00358             result                        = MBI->set_connectivity( *i, new_face_conn, 4 );
00359             if( MB_SUCCESS != result ) return result;
00360         }
00361     }
00362     return MB_SUCCESS;
00363 }
00364 
00365 /* Isolate dead code in a preproc-killed block - sjackson 11/22/10 */
00366 #if 0
00367 
00368 /* qsort int comparison function */
00369 int handle_compare(const void *a, const void *b)
00370 {
00371   const EntityHandle *ia = (const EntityHandle *)a; // casting pointer types
00372   const EntityHandle *ib = (const EntityHandle *)b;
00373   return *ia - *ib;
00374   /* integer comparison: returns negative if b > a
00375    and positive if a > b */
00376 }
00377 
00378 // qsort face comparison function. assume each face has 4 nodes
00379 int compare_face(const void *a, const void *b)
00380 {
00381   EntityHandle *ia = (EntityHandle *)a;
00382   EntityHandle *ib = (EntityHandle *)b;
00383   if(*ia == *ib)
00384   {
00385     if(*(ia+1) == *(ib+1))
00386     {
00387       if(*(ia+2) == *(ib+2))
00388       {
00389         return (int)(*(ia+3) - *(ib+3));
00390       }
00391       else
00392       {
00393         return (int)(*(ia+2) - *(ib+2));
00394       }
00395     }
00396     else
00397     {
00398       return (int)(*(ia+1) - *(ib+1));
00399     }
00400   }
00401   else
00402   {
00403     return (int)(*ia - *ib);
00404   }
00405 }
00406 
00407 // Use this to get quad faces from hex elems.
00408 ErrorCode skin_hex_elems(Interface *MBI, Range elems, const int dim,
00409     Range &faces )
00410 {
00411   // get faces of each hex
00412   const int nodes_per_face = 4;
00413   const unsigned int faces_per_elem = 6;
00414   unsigned int n_faces = faces_per_elem*elems.size();
00415   EntityHandle f[n_faces][nodes_per_face];
00416   ErrorCode result;
00417   int counter = 0;
00418   for(Range::iterator i=elems.begin(); i!=elems.end(); ++i)
00419   {
00420     Range elem_faces;
00421     result = MBI->get_adjacencies( &(*i), 1, 2, true, elem_faces );
00422     if(MB_SUCCESS != result) return result;
00423     if(faces_per_elem != elem_faces.size()) return MB_INVALID_SIZE;
00424     for(Range::iterator j=elem_faces.begin(); j!=elem_faces.end(); ++j)
00425     {
00426       const EntityHandle *conn;
00427       int n_nodes;
00428       ErrorCode result = MBI->get_connectivity( *j, conn, n_nodes );
00429       if(MB_SUCCESS != result) return result;
00430       if(nodes_per_face != n_nodes) return MB_INVALID_SIZE;
00431       // Sort the node handles of the face
00432       for(int k=0; k<nodes_per_face; ++k) f[counter][k] = conn[k];
00433       qsort( &f[counter][0], nodes_per_face, sizeof(EntityHandle), handle_compare );
00434       ++counter;
00435     }
00436   }
00437 
00438   // Sort the faces by the first node handle, then second node, then third node...
00439   qsort( &f[0][0], n_faces, nodes_per_face*sizeof(EntityHandle), compare_face );
00440 
00441   // if a face has 1 or more duplicates, it is not on the skin
00442   faces.clear();
00443   for(unsigned int i=0; i<n_faces; ++i)
00444   {
00445     // if the last face is tested, it must be on the skin
00446     if(n_faces-1 == i)
00447     {
00448       Range face_handle;
00449       result = MBI->get_adjacencies( &(f[i][0]), nodes_per_face, 2, false, face_handle );
00450       if(MB_SUCCESS != result) return result;
00451       if(1 != face_handle.size()) return MB_INVALID_SIZE;
00452       faces.insert( face_handle.front() );
00453       // Due to sort, if a duplicate exists it must be next
00454     }
00455     else if( f[i][0]==f[i+1][0] && f[i][1]==f[i+1][1] &&
00456         f[i][2]==f[i+1][2] && f[i][3]==f[i+1][3] )
00457     {
00458       ++i;
00459       while( f[i][0]==f[i+1][0] && f[i][1]==f[i+1][1] &&
00460           f[i][2]==f[i+1][2] && f[i][3]==f[i+1][3] )
00461       {
00462         std::cout << "    skin WARNING: non-manifold face" << std::endl;
00463         ++i;
00464       }
00465       // otherwise it is on the skin
00466     }
00467     else
00468     {
00469       Range face_handle;
00470       result = MBI->get_adjacencies( &(f[i][0]), nodes_per_face, 2, false, face_handle );
00471       if(MB_SUCCESS != result) return result;
00472       if(1 != face_handle.size()) return MB_INVALID_SIZE;
00473       faces.insert( face_handle.front() );
00474     }
00475   }
00476 
00477   return MB_SUCCESS;
00478 }
00479 
00480 #endif  // dead code isolation
00481 // Given a 1D array of data, axis labels, title, and number of bins, create a
00482 // histogram.
00483 void plot_histogram( const std::string& title, const std::string& x_axis_label, const std::string& y_axis_label,
00484                      const int n_bins, const double data[], const int n_data )
00485 {
00486     // find max and min
00487     double min = std::numeric_limits< double >::max();
00488     double max = -std::numeric_limits< double >::max();
00489     for( int i = 0; i < n_data; ++i )
00490     {
00491         if( min > data[i] ) min = data[i];
00492         if( max < data[i] ) max = data[i];
00493     }
00494 
00495     // make bins for histogram
00496     double bin_width = ( max - min ) / n_bins;
00497     std::vector< int > bins( n_bins );
00498     for( int i = 0; i < n_bins; ++i )
00499         bins[i] = 0;
00500 
00501     // fill the bins
00502     for( int i = 0; i < n_data; ++i )
00503     {
00504         double diff = data[i] - min;
00505         int bin     = diff / bin_width;
00506         if( 9 < bin ) bin = 9;  // cheap fix for numerical precision error
00507         if( 0 > bin ) bin = 0;  // cheap fix for numerical precision error
00508         ++bins[bin];
00509     }
00510 
00511     // create bars
00512     int max_bin = 0;
00513     for( int i = 0; i < n_bins; ++i )
00514         if( max_bin < bins[i] ) max_bin = bins[i];
00515     int bar_height;
00516     int max_bar_chars = 72;
00517     std::vector< std::string > bars( n_bins );
00518     for( int i = 0; i < n_bins; ++i )
00519     {
00520         bar_height = ( max_bar_chars * bins[i] ) / max_bin;
00521         for( int j = 0; j < bar_height; ++j )
00522             bars[i] += "*";
00523     }
00524 
00525     // print histogram header
00526     std::cout << std::endl;
00527     std::cout << "                                 " << title << std::endl;
00528 
00529     // print results
00530     std::cout.width( 15 );
00531     std::cout << min << std::endl;
00532     for( int i = 0; i < n_bins; ++i )
00533     {
00534         std::cout.width( 15 );
00535         std::cout << min + ( ( i + 1 ) * bin_width );
00536         std::cout.width( max_bar_chars );
00537         std::cout << bars[i] << bins[i] << std::endl;
00538     }
00539 
00540     // print histogram footer
00541     std::cout.width( 15 );
00542     std::cout << y_axis_label;
00543     std::cout.width( max_bar_chars );
00544     std::cout << " " << x_axis_label << std::endl;
00545     std::cout << std::endl;
00546 }
00547 
00548 // This is a helper function that creates data and labels for histograms.
00549 void generate_plots( const double orig[], const double defo[], const int n_elems, const std::string& time_step )
00550 {
00551 
00552     // find volume ratio then max and min
00553     double* ratio = new double[n_elems];
00554     for( int i = 0; i < n_elems; ++i )
00555         ratio[i] = ( defo[i] - orig[i] ) / orig[i];
00556 
00557     plot_histogram( "Predeformed Element Volume", "Num_Elems", "Volume [cc]", 10, orig, n_elems );
00558     std::string title = "Element Volume Change Ratio at Time Step " + time_step;
00559     plot_histogram( title, "Num_Elems", "Volume Ratio", 10, ratio, n_elems );
00560     delete[] ratio;
00561 }
00562 
00563 // Given four nodes, calculate the tet volume.
00564 inline static double tet_volume( const CartVect& v0, const CartVect& v1, const CartVect& v2, const CartVect& v3 )
00565 {
00566     return 1. / 6. * ( ( ( v1 - v0 ) * ( v2 - v0 ) ) % ( v3 - v0 ) );
00567 }
00568 
00569 // Measure and tet volume are taken from measure.cpp
00570 double measure( Interface* MBI, const EntityHandle element )
00571 {
00572     EntityType type = MBI->type_from_handle( element );
00573 
00574     const EntityHandle* conn;
00575     int num_vertices;
00576     ErrorCode result = MBI->get_connectivity( element, conn, num_vertices );
00577     if( MB_SUCCESS != result ) return result;
00578 
00579     std::vector< CartVect > coords( num_vertices );
00580     result = MBI->get_coords( conn, num_vertices, coords[0].array() );
00581     if( MB_SUCCESS != result ) return result;
00582 
00583     switch( type )
00584     {
00585         case MBEDGE:
00586             return ( coords[0] - coords[1] ).length();
00587         case MBTRI:
00588             return 0.5 * ( ( coords[1] - coords[0] ) * ( coords[2] - coords[0] ) ).length();
00589         case MBQUAD:
00590             num_vertices = 4;
00591         case MBPOLYGON: {
00592             CartVect mid( 0, 0, 0 );
00593             for( int i = 0; i < num_vertices; ++i )
00594                 mid += coords[i];
00595             mid /= num_vertices;
00596 
00597             double sum = 0.0;
00598             for( int i = 0; i < num_vertices; ++i )
00599             {
00600                 int j = ( i + 1 ) % num_vertices;
00601                 sum += ( ( mid - coords[i] ) * ( mid - coords[j] ) ).length();
00602             }
00603             return 0.5 * sum;
00604         }
00605         case MBTET:
00606             return tet_volume( coords[0], coords[1], coords[2], coords[3] );
00607         case MBPYRAMID:
00608             return tet_volume( coords[0], coords[1], coords[2], coords[4] ) +
00609                    tet_volume( coords[0], coords[2], coords[3], coords[4] );
00610         case MBPRISM:
00611             return tet_volume( coords[0], coords[1], coords[2], coords[5] ) +
00612                    tet_volume( coords[3], coords[5], coords[4], coords[0] ) +
00613                    tet_volume( coords[1], coords[4], coords[5], coords[0] );
00614         case MBHEX:
00615             return tet_volume( coords[0], coords[1], coords[3], coords[4] ) +
00616                    tet_volume( coords[7], coords[3], coords[6], coords[4] ) +
00617                    tet_volume( coords[4], coords[5], coords[1], coords[6] ) +
00618                    tet_volume( coords[1], coords[6], coords[3], coords[4] ) +
00619                    tet_volume( coords[2], coords[6], coords[3], coords[1] );
00620         default:
00621             return 0.0;
00622     }
00623 }
00624 
00625 /* Calculate the signed volumes beneath the surface (x 6.0). Use the triangle's
00626  canonical sense. Do not take sense tags into account. Code taken from
00627  DagMC::measure_volume.
00628 
00629  Special Case: If the surface is planar, and the plane includes the origin,
00630  the signed volume will be ~0. If the signed volume is ~0 then offset everything
00631  by a random amount and try again. */
00632 ErrorCode get_signed_volume( Interface* MBI, const EntityHandle surf_set, const CartVect& offset,
00633                              double& signed_volume )
00634 {
00635     ErrorCode rval;
00636     Range tris;
00637     rval = MBI->get_entities_by_type( surf_set, MBTRI, tris );
00638     if( MB_SUCCESS != rval ) return rval;
00639     signed_volume = 0.0;
00640     const EntityHandle* conn;
00641     int len;
00642     CartVect coords[3];
00643     for( Range::iterator j = tris.begin(); j != tris.end(); ++j )
00644     {
00645         rval = MBI->get_connectivity( *j, conn, len, true );
00646         if( MB_SUCCESS != rval ) return rval;
00647         if( 3 != len ) return MB_INVALID_SIZE;
00648         rval = MBI->get_coords( conn, 3, coords[0].array() );
00649         if( MB_SUCCESS != rval ) return rval;
00650 
00651         // Apply offset to avoid calculating 0 for cases when the origin is in the
00652         // plane of the surface.
00653         for( unsigned int k = 0; k < 3; ++k )
00654         {
00655             coords[k][0] += offset[0];
00656             coords[k][1] += offset[1];
00657             coords[k][2] += offset[2];
00658         }
00659 
00660         coords[1] -= coords[0];
00661         coords[2] -= coords[0];
00662         signed_volume += ( coords[0] % ( coords[1] * coords[2] ) );
00663     }
00664     return MB_SUCCESS;
00665 }
00666 
00667 // The cgm and cub surfaces may not have the same sense. Create triangles that
00668 // represent the quads in the cub surface. Calculate the signed volume of both
00669 // the cgm and cub surface. If they are different, change the cgm sense so that
00670 // it matches the sense of the cub surface.
00671 ErrorCode fix_surface_senses( Interface* MBI, const EntityHandle cgm_file_set, const EntityHandle cub_file_set,
00672                               const Tag idTag, const Tag dimTag, const Tag senseTag, const bool debug )
00673 {
00674     ErrorCode result;
00675     const int two               = 2;
00676     const void* const two_val[] = { &two };
00677     Range cgm_surfs;
00678     result = MBI->get_entities_by_type_and_tag( cgm_file_set, MBENTITYSET, &dimTag, two_val, 1, cgm_surfs );
00679     if( MB_SUCCESS != result ) return result;
00680     for( Range::iterator i = cgm_surfs.begin(); i != cgm_surfs.end(); ++i )
00681     {
00682         int surf_id;
00683         result = MBI->tag_get_data( idTag, &( *i ), 1, &surf_id );
00684         if( MB_SUCCESS != result ) return result;
00685 
00686         // Find the meshed surface set with the same id
00687         Range cub_surf;
00688         const Tag tags[]             = { idTag, dimTag };
00689         const void* const tag_vals[] = { &surf_id, &two };
00690         result = MBI->get_entities_by_type_and_tag( cub_file_set, MBENTITYSET, tags, tag_vals, 2, cub_surf );
00691         if( MB_SUCCESS != result ) return result;
00692         if( 1 != cub_surf.size() )
00693         {
00694             std::cout << "  Surface " << surf_id << ": no meshed representation found, using CAD representation instead"
00695                       << std::endl;
00696             continue;
00697         }
00698 
00699         // Get tris that represent the quads of the cub surf
00700         Range quads;
00701         result = MBI->get_entities_by_type( cub_surf.front(), MBQUAD, quads );
00702         if( MB_SUCCESS != result ) return result;
00703         Range cub_tris;
00704         result = make_tris_from_quads( MBI, quads, cub_tris );
00705         if( MB_SUCCESS != result ) return result;
00706 
00707         // Add the tris to the same surface meshset as the quads are inside.
00708         result = MBI->add_entities( cub_surf.front(), cub_tris );
00709         if( MB_SUCCESS != result ) return result;
00710 
00711         // get the signed volume for each surface representation. Keep trying until
00712         // The signed volumes are much greater than numerical precision. Planar
00713         // surfaces will have a signed volume of zero if the plane goes through the
00714         // origin, unless we apply an offset.
00715         const int n_attempts        = 100;
00716         const int max_random        = 10;
00717         const double min_signed_vol = 0.1;
00718         double cgm_signed_vol, cub_signed_vol;
00719         for( int j = 0; j < n_attempts; ++j )
00720         {
00721             cgm_signed_vol = 0;
00722             cub_signed_vol = 0;
00723             CartVect offset( std::rand() % max_random, std::rand() % max_random, std::rand() % max_random );
00724             result = get_signed_volume( MBI, *i, offset, cgm_signed_vol );
00725             if( MB_SUCCESS != result ) return result;
00726             result = get_signed_volume( MBI, cub_surf.front(), offset, cub_signed_vol );
00727             if( MB_SUCCESS != result ) return result;
00728             if( debug )
00729                 std::cout << "  surf_id=" << surf_id << " cgm_signed_vol=" << cgm_signed_vol
00730                           << " cub_signed_vol=" << cub_signed_vol << std::endl;
00731             if( fabs( cgm_signed_vol ) > min_signed_vol && fabs( cub_signed_vol ) > min_signed_vol ) break;
00732             if( n_attempts == j + 1 )
00733             {
00734                 std::cout << "error: signed volume could not be calculated unambiguously" << std::endl;
00735                 return MB_FAILURE;
00736             }
00737         }
00738 
00739         // If the sign is different, reverse the cgm senses so that both
00740         // representations have the same signed volume.
00741         if( ( cgm_signed_vol < 0 && cub_signed_vol > 0 ) || ( cgm_signed_vol > 0 && cub_signed_vol < 0 ) )
00742         {
00743             EntityHandle cgm_surf_volumes[2], reversed_cgm_surf_volumes[2];
00744             result = MBI->tag_get_data( senseTag, &( *i ), 1, cgm_surf_volumes );
00745             if( MB_SUCCESS != result ) return result;
00746             if( MB_SUCCESS != result ) return result;
00747 
00748             reversed_cgm_surf_volumes[0] = cgm_surf_volumes[1];
00749             reversed_cgm_surf_volumes[1] = cgm_surf_volumes[0];
00750 
00751             result = MBI->tag_set_data( senseTag, &( *i ), 1, reversed_cgm_surf_volumes );
00752             if( MB_SUCCESS != result ) return result;
00753         }
00754     }
00755 
00756     return MB_SUCCESS;
00757 }
00758 
00759 // The quads in the cub_file_set have been updated for dead elements. For each
00760 // cgm_surf, if there exists a cub_surf with the same id, replace the cgm tris
00761 // with cub_tris (created from the quads). Note the a surface that is not
00762 // meshed (in cub file) will not be effected.
00763 ErrorCode replace_faceted_cgm_surfs( Interface* MBI, const EntityHandle cgm_file_set, const EntityHandle cub_file_set,
00764                                      const Tag idTag, const Tag dimTag, const bool debug )
00765 {
00766     ErrorCode result;
00767     const int two               = 2;
00768     const void* const two_val[] = { &two };
00769     Range cgm_surfs;
00770     result = MBI->get_entities_by_type_and_tag( cgm_file_set, MBENTITYSET, &dimTag, two_val, 1, cgm_surfs );
00771     if( MB_SUCCESS != result ) return result;
00772 
00773     for( Range::iterator i = cgm_surfs.begin(); i != cgm_surfs.end(); ++i )
00774     {
00775         int surf_id;
00776         result = MBI->tag_get_data( idTag, &( *i ), 1, &surf_id );
00777         if( MB_SUCCESS != result ) return result;
00778         if( debug ) std::cout << "surf_id=" << surf_id << std::endl;
00779 
00780         // Find the meshed surface set with the same id
00781         Range cub_surf;
00782         const Tag tags[]             = { idTag, dimTag };
00783         const void* const tag_vals[] = { &surf_id, &two };
00784         result = MBI->get_entities_by_type_and_tag( cub_file_set, MBENTITYSET, tags, tag_vals, 2, cub_surf );
00785         if( MB_SUCCESS != result ) return result;
00786         if( 1 != cub_surf.size() )
00787         {
00788             std::cout << "  Surface " << surf_id << ": no meshed representation found, using CAD representation instead"
00789                       << std::endl;
00790             continue;
00791         }
00792 
00793         // Get tris that represent the quads of the cub surf
00794         Range quads;
00795         result = MBI->get_entities_by_type( cub_surf.front(), MBQUAD, quads );
00796         if( MB_SUCCESS != result ) return result;
00797 
00798         Range cub_tris;
00799         result = make_tris_from_quads( MBI, quads, cub_tris );
00800         if( MB_SUCCESS != result ) return result;
00801 
00802         // Remove the tris from the cgm surf. Don't forget to remove them from the
00803         // cgm_file_set because it is not TRACKING.
00804         Range cgm_tris;
00805         result = MBI->get_entities_by_type( *i, MBTRI, cgm_tris );
00806         if( MB_SUCCESS != result ) return result;
00807         result = MBI->remove_entities( *i, cgm_tris );
00808         if( MB_SUCCESS != result ) return result;
00809         result = MBI->remove_entities( cgm_file_set, cgm_tris );
00810         if( MB_SUCCESS != result ) return result;
00811         result = MBI->delete_entities( cgm_tris );
00812         if( MB_SUCCESS != result ) return result;
00813 
00814         // Add the cub_tris to the cgm_surf
00815         result = MBI->add_entities( *i, cub_tris );
00816         if( MB_SUCCESS != result ) return result;
00817     }
00818 
00819     return MB_SUCCESS;
00820 }
00821 
00822 // Dead elements need removed from the simulation. This is done by removing them
00823 // from their volume set and adding them to the implicit complement. New surfaces
00824 // must be created for this.
00825 // IF MODIFYING THIS CODE, BE AWARE THAT DEAD ELEMENTS CAN BE ADJACENT TO MORE
00826 // THAN ONE SURFACE, MAKING THE ASSOCIATION BETWEEN NEWLY EXPOSED AND EXISTING
00827 // SURFACES AMBIGUOUS.
00828 ErrorCode add_dead_elems_to_impl_compl( Interface* MBI, const EntityHandle cgm_file_set,
00829                                         const EntityHandle cub_file_set, const Tag idTag, const Tag dimTag,
00830                                         const Tag categoryTag, const Tag senseTag, const bool debug )
00831 {
00832 
00833     // Get the cgm surfaces
00834     ErrorCode result;
00835     const int two               = 2;
00836     const void* const two_val[] = { &two };
00837     Range cgm_surfs;
00838     result = MBI->get_entities_by_type_and_tag( cgm_file_set, MBENTITYSET, &dimTag, two_val, 1, cgm_surfs );
00839     if( MB_SUCCESS != result ) return result;
00840 
00841     // Get the maximum surface id. This is so that new surfaces do not have
00842     // duplicate ids.
00843     int max_surf_id = -1;
00844     for( Range::const_iterator i = cgm_surfs.begin(); i != cgm_surfs.end(); ++i )
00845     {
00846         int surf_id;
00847         result = MBI->tag_get_data( idTag, &( *i ), 1, &surf_id );
00848         if( MB_SUCCESS != result ) return result;
00849         if( max_surf_id < surf_id ) max_surf_id = surf_id;
00850     }
00851     std::cout << "  Maximum surface id=" << max_surf_id << std::endl;
00852 
00853     // For each cgm volume, does a cub volume with the same id exist?
00854     const int three               = 3;
00855     const void* const three_val[] = { &three };
00856     Range cgm_vols;
00857     result = MBI->get_entities_by_type_and_tag( cgm_file_set, MBENTITYSET, &dimTag, three_val, 1, cgm_vols );
00858     if( MB_SUCCESS != result ) return result;
00859 
00860     // get the corresponding cub volume
00861     for( Range::iterator i = cgm_vols.begin(); i != cgm_vols.end(); ++i )
00862     {
00863         int vol_id;
00864         result = MBI->tag_get_data( idTag, &( *i ), 1, &vol_id );
00865         assert( MB_SUCCESS == result );
00866         if( MB_SUCCESS != result ) return result;
00867         std::cout << "  Volume " << vol_id;
00868 
00869         // Find the meshed vol set with the same id
00870         Range cub_vol;
00871         const Tag tags[]             = { idTag, dimTag };
00872         const void* const tag_vals[] = { &vol_id, &three };
00873         result = MBI->get_entities_by_type_and_tag( cub_file_set, MBENTITYSET, tags, tag_vals, 2, cub_vol );
00874         assert( MB_SUCCESS == result );
00875         if( MB_SUCCESS != result ) return result;
00876         if( 1 != cub_vol.size() )
00877         {
00878             std::cout << ": no meshed representation found" << std::endl;
00879             continue;
00880         }
00881         else
00882         {
00883             std::cout << std::endl;
00884         }
00885 
00886         // get the mesh elements of the volume.
00887         Range elems;
00888         result = MBI->get_entities_by_type( cub_vol.front(), MBHEX, elems );
00889         assert( MB_SUCCESS == result );
00890         if( MB_SUCCESS != result ) return result;
00891         if( debug ) std::cout << "    found " << elems.size() << " hex elems" << std::endl;
00892 
00893         // skin the volumes
00894         Skinner tool( MBI );
00895         Range skin_faces;
00896         result = tool.find_skin( 0, elems, 2, skin_faces, true );
00897         assert( MB_SUCCESS == result );
00898         if( MB_SUCCESS != result ) return result;
00899 
00900         // Reconcile the difference between faces of the cub file surfaces and skin
00901         // faces of the 3D mesh. The difference exists because dead elements have been
00902         // removed. Faces are divided into:
00903         //   cub_faces    (in) - the faces in the cub file surface
00904         //   skin_faces   (in) - the faces of the 3D mesh elements
00905         //   common_faces (out)- the faces common to both the cub file surface and 3D mesh
00906         //                       they are still adjacent to this volume
00907         //   old_faces    (out)- the faces of the cub surface not on the 3D mesh skin
00908         //                       they are no longer adjacent to this vol
00909         //   new_faces    (out)- the faces of the 3D mesh skin not on the cub surface
00910         //                       they are now adjacent to this volume
00911 
00912         // get cub child surfaces.
00913         Range cub_surfs;
00914         result = MBI->get_child_meshsets( cub_vol.front(), cub_surfs );
00915         assert( MB_SUCCESS == result );
00916         if( MB_SUCCESS != result ) return result;
00917         for( Range::iterator j = cub_surfs.begin(); j != cub_surfs.end(); ++j )
00918         {
00919             int surf_id;
00920             result = MBI->tag_get_data( idTag, &( *j ), 1, &surf_id );
00921             assert( MB_SUCCESS == result );
00922             if( MB_SUCCESS != result ) return result;
00923             // get the quads on each surface
00924             Range cub_faces;
00925             result = MBI->get_entities_by_type( *j, MBQUAD, cub_faces );
00926             assert( MB_SUCCESS == result );
00927             if( MB_SUCCESS != result ) return result;
00928             // Meshed volumes must have meshed surfaces
00929             if( cub_faces.empty() )
00930             {
00931                 std::cout << "    Surface " << surf_id << ": contains no meshed faces" << std::endl;
00932                 // return MB_ENTITY_NOT_FOUND;
00933             }
00934             // get the faces common to both the skin and this surface
00935             Range common_faces = intersect( cub_faces, skin_faces );
00936             // find the surface faces not on the skin - these are old and need removed
00937             Range old_faces = subtract( cub_faces, common_faces );
00938             result          = MBI->remove_entities( *j, old_faces );
00939             assert( MB_SUCCESS == result );
00940             if( MB_SUCCESS != result ) return result;
00941 
00942             // remove the common faces from the skin faces
00943             skin_faces = subtract( skin_faces, common_faces );
00944             // If no old faces exist we are done
00945             if( old_faces.empty() ) continue;
00946             std::cout << "    Surface " << surf_id << ": " << old_faces.size() << " old faces removed" << std::endl;
00947             // Place the old faces in a new surface, because they may still be adjacent
00948             // to 3D mesh in another volume. Get the parent vols of the surface.
00949             Range cgm_surf;
00950             const Tag tags2[]             = { idTag, dimTag };
00951             const void* const tag_vals2[] = { &surf_id, &two };
00952             result = MBI->get_entities_by_type_and_tag( cgm_file_set, MBENTITYSET, tags2, tag_vals2, 2, cgm_surf );
00953             assert( MB_SUCCESS == result );
00954             if( MB_SUCCESS != result ) return result;
00955             if( 1 != cgm_surf.size() )
00956             {
00957                 std::cout << "invalid size" << std::endl;
00958                 return MB_INVALID_SIZE;
00959             }
00960             EntityHandle cgm_vols2[2], cub_vols[2];
00961             result = MBI->tag_get_data( senseTag, cgm_surf, &cgm_vols2 );
00962             assert( MB_SUCCESS == result );
00963             if( MB_SUCCESS != result ) return result;
00964             result = MBI->tag_get_data( senseTag, &( *j ), 1, &cub_vols );
00965             assert( MB_SUCCESS == result );
00966             if( MB_SUCCESS != result ) return result;
00967             // for the new surf, replace the current volume with the impl compl vol.
00968             // This is because the faces that no longer exist will become adjacent to
00969             // the impl compl
00970             if( *i == cgm_vols2[0] )
00971             {
00972                 cgm_vols2[0] = 0;
00973                 cub_vols[0]  = 0;
00974             }
00975             if( *i == cgm_vols2[1] )
00976             {
00977                 cgm_vols2[1] = 0;
00978                 cub_vols[1]  = 0;
00979             }
00980             // If both sides of the surface are the impl comp, do not create the surface.
00981             if( 0 == cgm_vols2[0] && 0 == cgm_vols2[1] )
00982             {
00983                 std::cout << "    New surface was not created for old faces because both parents "
00984                              "are impl_compl volume "
00985                           << std::endl;
00986                 continue;
00987             }
00988             // build the new surface.
00989             EntityHandle new_cgm_surf, new_cub_surf;
00990             ++max_surf_id;
00991             result = build_new_surface( MBI, new_cgm_surf, cgm_vols2[0], cgm_vols2[1], max_surf_id, dimTag, idTag,
00992                                         categoryTag, senseTag );
00993             assert( MB_SUCCESS == result );
00994             if( MB_SUCCESS != result ) return result;
00995             result = build_new_surface( MBI, new_cub_surf, cub_vols[0], cub_vols[1], max_surf_id, dimTag, idTag,
00996                                         categoryTag, senseTag );
00997             assert( MB_SUCCESS == result );
00998             if( MB_SUCCESS != result ) return result;
00999             // add the new surface to the file set and populate it with the old faces
01000             result = MBI->add_entities( cgm_file_set, &new_cgm_surf, 1 );
01001             assert( MB_SUCCESS == result );
01002             if( MB_SUCCESS != result ) return result;
01003             assert( MB_SUCCESS == result );
01004             result = MBI->add_entities( cub_file_set, &new_cub_surf, 1 );
01005             if( MB_SUCCESS != result ) return result;
01006             assert( MB_SUCCESS == result );
01007             result = MBI->add_entities( new_cub_surf, old_faces );
01008             assert( MB_SUCCESS == result );
01009             if( MB_SUCCESS != result ) return result;
01010 
01011             std::cout << "    Surface " << max_surf_id << ": created for " << old_faces.size() << " old faces"
01012                       << std::endl;
01013         }
01014 
01015         // the remaining skin faces are newly exposed faces
01016         Range new_faces = skin_faces;
01017 
01018         // new skin faces must be assigned to a surface
01019         if( new_faces.empty() ) continue;
01020         std::cout << "    Surface " << max_surf_id + 1 << ": created for " << new_faces.size() << " new faces"
01021                   << std::endl;
01022 
01023         // Ensure that faces are oriented outwards
01024         result = orient_faces_outward( MBI, new_faces, debug );
01025         assert( MB_SUCCESS == result );
01026         if( MB_SUCCESS != result ) return result;
01027 
01028         // Create the new surface.
01029         EntityHandle new_cgm_surf, new_cub_surf;
01030         ++max_surf_id;
01031         result = build_new_surface( MBI, new_cgm_surf, *i, 0, max_surf_id, dimTag, idTag, categoryTag, senseTag );
01032         assert( MB_SUCCESS == result );
01033         if( MB_SUCCESS != result ) return result;
01034         result = build_new_surface( MBI, new_cub_surf, cub_vol.front(), 0, max_surf_id, dimTag, idTag, categoryTag,
01035                                     senseTag );
01036         assert( MB_SUCCESS == result );
01037         if( MB_SUCCESS != result ) return result;
01038         // Insert the new surf into file sets and populate it with faces.
01039         result = MBI->add_entities( cgm_file_set, &new_cgm_surf, 1 );
01040         assert( MB_SUCCESS == result );
01041         if( MB_SUCCESS != result ) return result;
01042         result = MBI->add_entities( cub_file_set, &new_cub_surf, 1 );
01043         assert( MB_SUCCESS == result );
01044         if( MB_SUCCESS != result ) return result;
01045         result = MBI->add_entities( new_cub_surf, new_faces );
01046         assert( MB_SUCCESS == result );
01047         if( MB_SUCCESS != result ) return result;
01048     }
01049 
01050     return MB_SUCCESS;
01051 }
01052 
01053 /* Get the type of a file.
01054  Return value is one of the above constants
01055  */
01056 const char* get_geom_file_type( const char* filename );
01057 const char* get_geom_fptr_type( FILE* file );
01058 
01059 int is_cubit_file( FILE* file );
01060 int is_step_file( FILE* file );
01061 int is_iges_file( FILE* file );
01062 int is_acis_txt_file( FILE* file );
01063 int is_acis_bin_file( FILE* file );
01064 int is_occ_brep_file( FILE* file );
01065 
01066 double DEFAULT_DISTANCE = 0.001;
01067 double DEFAULT_LEN      = 0.0;
01068 int DEFAULT_NORM        = 5;
01069 
01070 // load cub file
01071 // load cgm file
01072 // for each surface
01073 //   convert cub surf quads to tris
01074 //   get signed volume from cgm and cub surf         MUST COME BEFORE COORD UPDATE, NEEDS MBTRIS
01075 //   reverse cgm surface sense if needed
01076 //   replace cgm surf tris with cub surf tris
01077 // measure volume of predeformed cub elements
01078 // convert cub volumes sets to tracking so that dead elems are removed from vol sets
01079 // update coordinates and delete dead elems
01080 // measure volume of deformed cub elems
01081 // print histogram of volume change
01082 // for each cub volume
01083 //   skin volume elems to get faces
01084 //   for each child cub surface
01085 //     assign old skin faces to a new surface in case they are adjacent to another volume
01086 //   orient each skin face outward
01087 //   assign new skin faces to a new surface
01088 // for each surface
01089 //   remove existing tris (from before the update)
01090 //   convert quads to tris
01091 // remove empty surfaces and volumes due to dead elements
01092 int main( int argc, char* argv[] )
01093 {
01094     clock_t start_time = clock();
01095 
01096     const bool debug      = false;
01097     const char* file_type = NULL;
01098 
01099     const char* cub_name  = 0;
01100     const char* exo_name  = 0;
01101     const char* out_name  = 0;
01102     const char* time_step = 0;
01103     const char* sat_name  = 0;
01104     double dist_tol = 0.001, len_tol = 0.0;
01105     int norm_tol = 5;
01106 
01107     if( 6 != argc && 9 != argc )
01108     {
01109         std::cerr << "To read meshed geometry for MOAB:" << std::endl;
01110         std::cerr << "$> <cub_file.cub> <acis_file.sat> <facet_tol> <output_file.h5m> conserve_mass<bool>" << std::endl;
01111         std::cerr << "To read meshed geometry for MOAB and update node coordinates:" << std::endl;
01112         std::cerr << "$> <cub_file.cub> <acis_file.sat> <facet_tol> <output_file.h5m> "
01113                      "<deformed_exo_file.e> time_step<int> check_vol_change<bool> conserve_mass<bool>"
01114                   << std::endl;
01115         exit( 4 );
01116     }
01117 
01118     // check filenames for proper suffix
01119     std::string temp;
01120     cub_name = argv[1];
01121     temp.assign( cub_name );
01122     if( std::string::npos == temp.find( ".cub" ) )
01123     {
01124         std::cerr << "cub_file does not have correct suffix" << std::endl;
01125         return 1;
01126     }
01127     sat_name = argv[2];  // needed because the cub file's embedded sat file does not have groups
01128     temp.assign( sat_name );
01129     if( std::string::npos == temp.find( ".sat" ) )
01130     {
01131         std::cerr << "sat_file does not have correct suffix" << std::endl;
01132         return 1;
01133     }
01134     out_name = argv[4];
01135     temp.assign( out_name );
01136     if( std::string::npos == temp.find( ".h5m" ) )
01137     {
01138         std::cerr << "out_file does not have correct suffix" << std::endl;
01139         return 1;
01140     }
01141 
01142     // get the facet tolerance
01143     dist_tol = atof( argv[3] );
01144     if( 0 > dist_tol || 1 < dist_tol )
01145     {
01146         std::cout << "error: facet_tolerance=" << dist_tol << std::endl;
01147         return 1;
01148     }
01149 
01150     // Should the nodes be updated?
01151     bool update_coords = false;
01152     if( 9 == argc )
01153     {
01154         exo_name = argv[5];
01155         temp.assign( exo_name );
01156         if( std::string::npos == temp.find( ".e" ) )
01157         {
01158             std::cerr << "e_file does not have correct suffix" << std::endl;
01159             return 1;
01160         }
01161         time_step     = argv[6];
01162         update_coords = true;
01163     }
01164 
01165     // Should the volume change be determined?
01166     bool determine_volume_change = false;
01167     if( 9 == argc )
01168     {
01169         temp.assign( argv[7] );
01170         if( std::string::npos != temp.find( "true" ) ) determine_volume_change = true;
01171     }
01172 
01173     // Should densities be changed to conserve mass?
01174     bool conserve_mass = false;
01175     if( 9 == argc )
01176     {
01177         temp.assign( argv[8] );
01178         if( std::string::npos != temp.find( "true" ) ) conserve_mass = true;
01179     }
01180     else if( 6 == argc )
01181     {
01182         temp.assign( argv[5] );
01183         if( std::string::npos != temp.find( "true" ) ) conserve_mass = true;
01184     }
01185 
01186     // Get CGM file type
01187     if( !file_type )
01188     {
01189         file_type = get_geom_file_type( cub_name );
01190         if( !file_type )
01191         {
01192             std::cerr << cub_name << " : unknown file type, try '-t'" << std::endl;
01193             exit( 1 );
01194         }
01195     }
01196 
01197     // Read the mesh from the cub file with Tqcdfr
01198     Core* MBI = new Core();
01199     ErrorCode result;
01200     EntityHandle cub_file_set;
01201     result = MBI->create_meshset( 0, cub_file_set );
01202     if( MB_SUCCESS != result ) return result;
01203     // Do not ignore the Cubit file version. In testing, a cub file from Cubit12
01204     // did not work.
01205     // char cub_options[256] = "120";
01206     // char cub_options[256] = "IGNORE_VERSION";
01207     // result = MBI->load_file(cub_name, &cub_file_set, cub_options, NULL, 0, 0);
01208     result = MBI->load_file( cub_name, &cub_file_set, 0, NULL, 0, 0 );
01209     if( MB_SUCCESS != result )
01210     {
01211         std::cout << "error: problem reading cub file" << std::endl;
01212         return result;
01213     }
01214     std::cout << "Mesh file read." << std::endl;
01215 
01216     // Read the ACIS file with ReadCGM
01217     char cgm_options[256];
01218     std::cout << "  facet tolerance=" << dist_tol << std::endl;
01219     sprintf( cgm_options,
01220              "CGM_ATTRIBS=yes;FACET_DISTANCE_TOLERANCE=%g;FACET_NORMAL_TOLERANCE=%d;MAX_FACET_EDGE_"
01221              "LENGTH=%g;",
01222              dist_tol, norm_tol, len_tol );
01223     EntityHandle cgm_file_set;
01224     result = MBI->create_meshset( 0, cgm_file_set );
01225     if( MB_SUCCESS != result ) return result;
01226     result = MBI->load_file( sat_name, &cgm_file_set, cgm_options, NULL, 0, 0 );
01227     if( MB_SUCCESS != result )
01228     {
01229         std::cout << "error: problem reading sat file" << std::endl;
01230         return result;
01231     }
01232     std::cout << "CAD file read." << std::endl;
01233 
01234     // Create tags
01235     Tag dimTag, idTag, categoryTag, senseTag, sizeTag, nameTag;
01236     result = MBI->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, dimTag, MB_TAG_SPARSE | MB_TAG_CREAT );
01237     if( MB_SUCCESS != result ) return result;
01238     idTag  = MBI->globalId_tag();
01239     result = MBI->tag_get_handle( CATEGORY_TAG_NAME, CATEGORY_TAG_SIZE, MB_TYPE_OPAQUE, categoryTag,
01240                                   MB_TAG_SPARSE | MB_TAG_CREAT );
01241     if( MB_SUCCESS != result ) return result;
01242     result = MBI->tag_get_handle( "GEOM_SENSE_2", 2, MB_TYPE_HANDLE, senseTag, MB_TAG_SPARSE | MB_TAG_CREAT );
01243     if( MB_SUCCESS != result ) return result;
01244     result = MBI->tag_get_handle( "GEOM_SIZE", 1, MB_TYPE_DOUBLE, sizeTag, MB_TAG_DENSE | MB_TAG_CREAT );
01245     if( MB_SUCCESS != result ) return result;
01246     result = MBI->tag_get_handle( NAME_TAG_NAME, NAME_TAG_SIZE, MB_TYPE_OPAQUE, nameTag, MB_TAG_SPARSE | MB_TAG_CREAT );
01247     if( MB_SUCCESS != result ) return result;
01248 
01249     // Create triangles from the quads of the cub surface sets and add them to the
01250     // cub surface sets. Get the signed volume of each surface for both cgm and
01251     // cub representations. Change the sense of the cgm representation to match
01252     // the cub representation.
01253     result = fix_surface_senses( MBI, cgm_file_set, cub_file_set, idTag, dimTag, senseTag, debug );
01254     if( MB_SUCCESS != result ) return result;
01255     std::cout << "Fixed CAD surface senses to match meshed surface senses." << std::endl;
01256 
01257     // Get the 3D elements in the cub file and measure their volume.
01258     Range orig_elems;
01259     std::vector< double > orig_size;
01260     if( determine_volume_change )
01261     {
01262         result = MBI->get_entities_by_dimension( 0, 3, orig_elems );
01263         if( MB_SUCCESS != result ) return result;
01264         orig_size.resize( orig_elems.size() );
01265         for( unsigned int i = 0; i < orig_elems.size(); ++i )
01266         {
01267             orig_size[i] = measure( MBI, orig_elems[i] );
01268         }
01269     }
01270 
01271     // Before updating the nodes and removing dead elements, force the cub volume
01272     // sets to track ownership so that dead elements will be deleted from the sets.
01273     const int three               = 3;
01274     const void* const three_val[] = { &three };
01275     Range cub_vols;
01276     result = MBI->get_entities_by_type_and_tag( cub_file_set, MBENTITYSET, &dimTag, three_val, 1, cub_vols );
01277     if( MB_SUCCESS != result ) return result;
01278     for( Range::const_iterator i = cub_vols.begin(); i != cub_vols.end(); ++i )
01279     {
01280         result = MBI->set_meshset_options( *i, MESHSET_TRACK_OWNER );
01281         if( MB_SUCCESS != result ) return result;
01282     }
01283 
01284     // Tag volume sets with the undeformed size of each volume.
01285     Range vol_sets;
01286     result = MBI->get_entities_by_type_and_tag( cgm_file_set, MBENTITYSET, &dimTag, three_val, 1, vol_sets );
01287     if( MB_SUCCESS != result ) return result;
01288 
01289     moab::GeomTopoTool gtt  = moab::GeomTopoTool( MBI, false );
01290     moab::GeomQueryTool gqt = moab::GeomQueryTool( &gtt );
01291 
01292     for( Range::const_iterator i = vol_sets.begin(); i != vol_sets.end(); ++i )
01293     {
01294         double size;
01295         result = gqt.measure_volume( *i, size );
01296         if( MB_SUCCESS != result ) return result;
01297         result = MBI->tag_set_data( sizeTag, &( *i ), 1, &size );
01298         if( MB_SUCCESS != result ) return result;
01299     }
01300 
01301     // Update the coordinates if needed. Do not do this before checking surface
01302     // sense, because the coordinate update could deform the surfaces too much
01303     // to make an accurate comparison.
01304     // The cub node ids are unique because cgm vertex ids are tagged on the vertex
01305     // meshset and not the vertex itself.
01306     // result = MBI->delete_entities( &cub_file_set, 1 );
01307     // if(MB_SUCCESS != result) return result;
01308     // Assume dead elements exist until I think of something better.
01309     bool dead_elements_exist = true;
01310     if( update_coords )
01311     {
01312         // ReadNCDF my_ex_reader(MBI);
01313         char exo_options[120] = "tdata=coord,";
01314         strcat( exo_options, time_step );
01315         strcat( exo_options, ",set" );
01316         // FileOptions exo_opts(exo_options)  ;
01317         // opts = "tdata=coord, 100, sum, temp.exo";
01318         // result =  my_ex_reader.load_file(exo_name, cgm_file_set, exo_opts, NULL, 0 , 0);
01319         // result =  my_ex_reader.load_file(exo_name, cub_file_set, exo_opts, NULL, 0 , 0);
01320         // result = my_ex_reader.load_file(exo_name, &cub_file_set, exo_opts, NULL, 0 , 0);
01321         MBI->load_file( exo_name, &cub_file_set, exo_options );
01322         if( MB_SUCCESS != result )
01323         {
01324             std::string last_error;
01325             MBI->get_last_error( last_error );
01326             std::cout << "coordinate update failed, " << last_error << std::endl;
01327             return result;
01328         }
01329         std::cout << "Updated mesh nodes with deformed coordinates from exodus file." << std::endl;
01330     }
01331 
01332     if( determine_volume_change )
01333     {
01334         // Dead elements have been removed by the deformation. Get the elements that
01335         // still exist.
01336         Range defo_elems;
01337         result = MBI->get_entities_by_dimension( 0, 3, defo_elems );
01338         if( MB_SUCCESS != result ) return result;
01339 
01340         // Determine the volume of the elements now that a deformation has been
01341         // applied. Condense the original array by removing dead elements.
01342         double* orig_size_condensed = new double[defo_elems.size()];
01343         double* defo_size_condensed = new double[defo_elems.size()];
01344         int j                       = 0;
01345         for( unsigned int i = 0; i < orig_elems.size(); ++i )
01346         {
01347             if( orig_elems[i] == defo_elems[j] )
01348             {
01349                 orig_size_condensed[j] = orig_size[i];
01350                 defo_size_condensed[j] = measure( MBI, defo_elems[j] );
01351                 ++j;
01352             }
01353         }
01354         generate_plots( orig_size_condensed, defo_size_condensed, defo_elems.size(), std::string( time_step ) );
01355         delete[] orig_size_condensed;  // can't use the same delete[] for both
01356         delete[] defo_size_condensed;
01357     }
01358 
01359     // Deal with dead elements. For now, add them to the impl_compl volume.
01360     // Extra surfaces are created to do this.
01361     if( update_coords && dead_elements_exist )
01362     {
01363         result = add_dead_elems_to_impl_compl( MBI, cgm_file_set, cub_file_set, idTag, dimTag, categoryTag, senseTag,
01364                                                debug );
01365         if( MB_SUCCESS != result ) return result;
01366         std::cout << "Placed dead elements in implicit complement volume and added required surfaces." << std::endl;
01367     }
01368 
01369     // The quads in the cub_file_set have been updated for dead elements. For each
01370     // cgm_surf, if there exists a cub_surf with the same id, replace the cgm tris
01371     // with cub_tris (created from the quads). Note the a surface that is not
01372     // meshed (in cub file) will not be effected.
01373     result = replace_faceted_cgm_surfs( MBI, cgm_file_set, cub_file_set, idTag, dimTag, debug );
01374     if( MB_SUCCESS != result ) return result;
01375     std::cout << "Replaced faceted CAD surfaces with meshed surfaces of triangles." << std::endl;
01376 
01377     result = remove_empty_cgm_surfs_and_vols( MBI, cgm_file_set, idTag, dimTag, debug );
01378     if( MB_SUCCESS != result ) return result;
01379     std::cout << "Removed surfaces and volumes that no longer have any triangles." << std::endl;
01380 
01381     // For each material, sum the volume. If the coordinates were updated for
01382     // deformation, summarize the volume change.
01383     result = summarize_cell_volume_change( MBI, cgm_file_set, categoryTag, dimTag, sizeTag, nameTag, idTag,
01384                                            conserve_mass, debug );
01385     if( MB_SUCCESS != result ) return result;
01386     std::cout << "Summarized the volume change of each material, with respect to the solid model." << std::endl;
01387 
01388     result = MBI->write_mesh( out_name, &cgm_file_set, 1 );
01389     if( MB_SUCCESS != result )
01390     {
01391         std::cout << "write mesh failed" << std::endl;
01392         return result;
01393     }
01394     std::cout << "Saved output file for mesh-based analysis." << std::endl;
01395 
01396     clock_t end_time = clock();
01397     std::cout << "  " << (double)( end_time - start_time ) / CLOCKS_PER_SEC << " seconds" << std::endl;
01398     std::cout << std::endl;
01399 
01400     return 0;
01401 }
01402 
01403 const char* get_geom_file_type( const char* name )
01404 {
01405     FILE* file;
01406     const char* result = 0;
01407 
01408     file = fopen( name, "r" );
01409     if( file )
01410     {
01411         result = get_geom_fptr_type( file );
01412         fclose( file );
01413     }
01414 
01415     return result;
01416 }
01417 
01418 const char* get_geom_fptr_type( FILE* file )
01419 {
01420     static const char* CUBIT_NAME = GF_CUBIT_FILE_TYPE;
01421     static const char* STEP_NAME  = GF_STEP_FILE_TYPE;
01422     static const char* IGES_NAME  = GF_IGES_FILE_TYPE;
01423     static const char* SAT_NAME   = GF_ACIS_TXT_FILE_TYPE;
01424     static const char* SAB_NAME   = GF_ACIS_BIN_FILE_TYPE;
01425     static const char* BREP_NAME  = GF_OCC_BREP_FILE_TYPE;
01426 
01427     if( is_cubit_file( file ) )
01428         return CUBIT_NAME;
01429     else if( is_step_file( file ) )
01430         return STEP_NAME;
01431     else if( is_iges_file( file ) )
01432         return IGES_NAME;
01433     else if( is_acis_bin_file( file ) )
01434         return SAB_NAME;
01435     else if( is_acis_txt_file( file ) )
01436         return SAT_NAME;
01437     else if( is_occ_brep_file( file ) )
01438         return BREP_NAME;
01439     else
01440         return 0;
01441 }
01442 
01443 int is_cubit_file( FILE* file )
01444 {
01445     unsigned char buffer[4];
01446     return !fseek( file, 0, SEEK_SET ) && fread( buffer, 4, 1, file ) && !memcmp( buffer, "CUBE", 4 );
01447 }
01448 
01449 int is_step_file( FILE* file )
01450 {
01451     unsigned char buffer[9];
01452     return !fseek( file, 0, SEEK_SET ) && fread( buffer, 9, 1, file ) && !memcmp( buffer, "ISO-10303", 9 );
01453 }
01454 
01455 int is_iges_file( FILE* file )
01456 {
01457     unsigned char buffer[10];
01458     return !fseek( file, 72, SEEK_SET ) && fread( buffer, 10, 1, file ) && !memcmp( buffer, "S      1\r\n", 10 );
01459 }
01460 
01461 int is_acis_bin_file( FILE* file )
01462 {
01463     char buffer[15];
01464     return !fseek( file, 0, SEEK_SET ) && fread( buffer, 15, 1, file ) && !memcmp( buffer, "ACIS BinaryFile", 9 );
01465 }
01466 
01467 int is_acis_txt_file( FILE* file )
01468 {
01469     char buffer[5];
01470     int version, length;
01471 
01472     if( fseek( file, 0, SEEK_SET ) || 2 != fscanf( file, "%d %*d %*d %*d %d ", &version, &length ) ) return 0;
01473 
01474     if( version < 1 || version > 0xFFFF ) return 0;
01475 
01476     // Skip application name
01477     if( fseek( file, length, SEEK_CUR ) ) return 0;
01478 
01479     // Read length of version string followed by first 5 characters
01480     if( 2 != fscanf( file, "%d %4s", &length, buffer ) ) return 0;
01481 
01482     return !strcmp( buffer, "ACIS" );
01483 }
01484 
01485 int is_occ_brep_file( FILE* file )
01486 {
01487     unsigned char buffer[6];
01488     return !fseek( file, 0, SEEK_SET ) && fread( buffer, 6, 1, file ) && !memcmp( buffer, "DBRep_", 6 );
01489 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines