MOAB: Mesh Oriented datABase  (version 5.3.0)
proj1.cpp
Go to the documentation of this file.
00001 /*
00002  * proj1.cpp
00003  *
00004  *  project on a sphere of radius R, delete sets if needed, and delete edges between parts
00005  *  (created by resolve shared ents)
00006  */
00007 
00008 #include "moab/Core.hpp"
00009 #include "moab/Interface.hpp"
00010 #include <iostream>
00011 #include <cmath>
00012 
00013 #include "moab/IntxMesh/IntxUtils.hpp"
00014 #include <cassert>
00015 using namespace moab;
00016 
00017 double radius = 1.;  // in m:  6371220.
00018 
00019 int main( int argc, char** argv )
00020 {
00021 
00022     bool delete_partition_sets = false;
00023 
00024     if( argc < 3 ) return 1;
00025 
00026     int index         = 1;
00027     char* input_mesh1 = argv[1];
00028     char* output      = argv[2];
00029     while( index < argc )
00030     {
00031         if( !strcmp( argv[index], "-R" ) )  // this is for radius to project
00032         {
00033             radius = atof( argv[++index] );
00034         }
00035         if( !strcmp( argv[index], "-DS" ) )  // delete partition sets
00036         {
00037             delete_partition_sets = true;
00038         }
00039 
00040         if( !strcmp( argv[index], "-h" ) )
00041         {
00042             std::cout << " usage: proj1 <input> <output> -R <value>  -DS (delete partition sets)\n";
00043             return 1;
00044         }
00045         index++;
00046     }
00047 
00048     Core moab;
00049     Interface& mb = moab;
00050 
00051     ErrorCode rval = mb.load_mesh( input_mesh1 );
00052 
00053     std::cout << " -R " << radius << " input: " << input_mesh1 << "  output: " << output << "\n";
00054 
00055     Range verts;
00056     rval = mb.get_entities_by_dimension( 0, 0, verts );
00057     if( MB_SUCCESS != rval ) return 1;
00058 
00059     double *x_ptr, *y_ptr, *z_ptr;
00060     int count;
00061     rval = mb.coords_iterate( verts.begin(), verts.end(), x_ptr, y_ptr, z_ptr, count );
00062     if( MB_SUCCESS != rval ) return 1;
00063     assert( count == (int)verts.size() );  // should end up with just one contiguous chunk of vertices
00064 
00065     for( int v = 0; v < count; v++ )
00066     {
00067         // EntityHandle v = verts[v];
00068         CartVect pos( x_ptr[v], y_ptr[v], z_ptr[v] );
00069         pos      = pos / pos.length();
00070         pos      = radius * pos;
00071         x_ptr[v] = pos[0];
00072         y_ptr[v] = pos[1];
00073         z_ptr[v] = pos[2];
00074     }
00075 
00076     Range edges;
00077     rval = mb.get_entities_by_dimension( 0, 1, edges );
00078     if( MB_SUCCESS != rval ) return 1;
00079     // write edges to a new set, and after that, write the set, delete the edges and the set
00080     EntityHandle sf1;
00081     rval = mb.create_meshset( MESHSET_SET, sf1 );
00082     if( MB_SUCCESS != rval ) return 1;
00083     rval = mb.add_entities( sf1, edges );
00084     if( MB_SUCCESS != rval ) return 1;
00085     rval = mb.write_mesh( "edgesOnly.h5m", &sf1, 1 );
00086     if( MB_SUCCESS != rval ) return 1;
00087     rval = mb.delete_entities( &sf1, 1 );
00088     if( MB_SUCCESS != rval ) return 1;
00089     mb.delete_entities( edges );
00090     mb.write_file( output );
00091 
00092     if( delete_partition_sets )
00093     {
00094         Tag par_tag;
00095         rval = mb.tag_get_handle( "PARALLEL_PARTITION", par_tag );
00096         if( MB_SUCCESS == rval )
00097 
00098         {
00099             Range par_sets;
00100             rval =
00101                 mb.get_entities_by_type_and_tag( 0, MBENTITYSET, &par_tag, NULL, 1, par_sets, moab::Interface::UNION );
00102             if( !par_sets.empty() ) mb.delete_entities( par_sets );
00103             mb.tag_delete( par_tag );
00104         }
00105     }
00106 
00107     return 0;
00108 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines