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