MOAB: Mesh Oriented datABase  (version 5.2.1)
intx_rll_cs_sphere_test.cpp
Go to the documentation of this file.
00001 /*
00002  * IntxRllCssphere_test.cpp
00003  */
00004 
00005 #include "moab/IntxMesh/IntxRllCssphere.hpp"
00006 #include "moab/IntxMesh/IntxUtils.hpp"
00007 #include "TestUtil.hpp"
00008 
00009 using namespace moab;
00010 
00011 int main( int argc, char* argv[] )
00012 {
00013     // check command line arg// Euler grid is red, arrival, Lagrangian is blue, departure
00014     // will will keep the
00015     const char* filename_mesh1 = STRINGIFY( MESHDIR ) "/mbcslam/outRLLMesh.g";
00016     const char* filename_mesh2 = STRINGIFY( MESHDIR ) "/mbcslam/outCSMesh.g";
00017     double R                   = 1.;  // input
00018     double epsrel              = 1.e-8;
00019     const char* newFile        = "intx.vtk";
00020     if( argc == 6 )
00021     {
00022         filename_mesh1 = argv[1];
00023         filename_mesh2 = argv[2];
00024         R              = atof( argv[3] );
00025         epsrel         = atof( argv[4] );
00026         newFile        = argv[5];
00027     }
00028     else
00029     {
00030         printf( "Usage: %s <mesh_filename1> <mesh_filename2> <radius> <epsrel> <newFile>\n", argv[0] );
00031         if( argc != 1 ) return 1;
00032         printf( "No files specified.  Defaulting to: %s  %s  %f %g %s\n", filename_mesh1, filename_mesh2, R, epsrel,
00033                 newFile );
00034     }
00035 
00036     // read meshes in 2 file sets
00037     Core moab;
00038     Interface* mb = &moab;  // global
00039     EntityHandle sf1, sf2;
00040     ErrorCode rval = mb->create_meshset( MESHSET_SET, sf1 );
00041     if( MB_SUCCESS != rval ) return 1;
00042     rval = mb->create_meshset( MESHSET_SET, sf2 );
00043     if( MB_SUCCESS != rval ) return 1;
00044     rval = mb->load_file( filename_mesh1, &sf1 );
00045     if( MB_SUCCESS != rval ) return 1;
00046     rval = mb->load_file( filename_mesh2, &sf2 );
00047     if( MB_SUCCESS != rval ) return 1;
00048 
00049     EntityHandle outputSet;
00050     rval = mb->create_meshset( MESHSET_SET, outputSet );
00051     if( MB_SUCCESS != rval ) return 1;
00052 
00053     // IntxUtils
00054     rval = IntxUtils::fix_degenerate_quads( mb, sf1 );
00055     if( MB_SUCCESS != rval ) return 1;
00056 
00057     IntxAreaUtils areaAdaptor;
00058     rval = areaAdaptor.positive_orientation( mb, sf1, R );
00059     if( MB_SUCCESS != rval ) return 1;
00060 
00061     rval = areaAdaptor.positive_orientation( mb, sf2, R );
00062     if( MB_SUCCESS != rval ) return 1;
00063 
00064     /*// set the edge tags on all elements sf1
00065     rval = set_edge_type_flag(mb, sf1); // form all edges, and set on them type 1 if constant
00066     latitude
00067     // add them to the set after this, just so we have them
00068     if (MB_SUCCESS != rval)
00069       return 1;*/
00070 
00071     IntxRllCssphere worker( mb );
00072 
00073     worker.set_error_tolerance( R * epsrel );
00074     // worker.SetEntityType(moab::MBQUAD);
00075     worker.set_radius( R );
00076     // worker.enable_debug();
00077     rval = worker.FindMaxEdges( sf1, sf2 );
00078     if( MB_SUCCESS != rval ) return 1;
00079 
00080     rval = worker.intersect_meshes( sf1, sf2, outputSet );
00081     // compute total area with 2 methods
00082 
00083     if( MB_SUCCESS != rval ) std::cout << " failed to intersect meshes\n";
00084     rval = mb->write_mesh( newFile, &outputSet, 1 );
00085     if( MB_SUCCESS != rval ) return 1;
00086     return 0;
00087 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines