intx_rll_cs_sphere_test.cpp

Go to the documentation of this file.

/*
 * IntxRllCssphere_test.cpp
 */
 
#include "moab/IntxMesh/IntxRllCssphere.hpp"
#include "moab/IntxMesh/IntxUtils.hpp"
#include "TestUtil.hpp"
 
using namespace moab;
 
int main( int argc, char* argv[] )
{
 // check command line arg// Euler grid is red, arrival, Lagrangian is blue, departure
 // will will keep the
 const char* filename_mesh1 = STRINGIFY( MESHDIR ) "/mbcslam/outRLLMesh.g";
 const char* filename_mesh2 = STRINGIFY( MESHDIR ) "/mbcslam/outCSMesh.g";
 double R = 1.; // input
 double epsrel = 1.e-8;
 const char* newFile = "intx.vtk";
 if( argc == 6 )
 {
 filename_mesh1 = argv[1];
 filename_mesh2 = argv[2];
 R = atof( argv[3] );
 epsrel = atof( argv[4] );
 newFile = argv[5];
 }
 else
 {
 printf( "Usage: %s <mesh_filename1> <mesh_filename2> <radius> <epsrel> <newFile>\n", argv[0] );
 if( argc != 1 ) return 1;
 printf( "No files specified. Defaulting to: %s %s %f %g %s\n", filename_mesh1, filename_mesh2, R, epsrel,
 newFile );
 }
 
 // read meshes in 2 file sets
 Core moab;
 Interface* mb = &moab; // global
 EntityHandle sf1, sf2;
 ErrorCode rval = mb->create_meshset( MESHSET_SET, sf1 );
 if( MB_SUCCESS != rval ) return 1;
 rval = mb->create_meshset( MESHSET_SET, sf2 );
 if( MB_SUCCESS != rval ) return 1;
 rval = mb->load_file( filename_mesh1, &sf1 );
 if( MB_SUCCESS != rval ) return 1;
 rval = mb->load_file( filename_mesh2, &sf2 );
 if( MB_SUCCESS != rval ) return 1;
 
 EntityHandle outputSet;
 rval = mb->create_meshset( MESHSET_SET, outputSet );
 if( MB_SUCCESS != rval ) return 1;
 
 // IntxUtils
 rval = IntxUtils::fix_degenerate_quads( mb, sf1 );
 if( MB_SUCCESS != rval ) return 1;
 
 IntxAreaUtils areaAdaptor;
 rval = areaAdaptor.positive_orientation( mb, sf1, R );
 if( MB_SUCCESS != rval ) return 1;
 
 rval = areaAdaptor.positive_orientation( mb, sf2, R );
 if( MB_SUCCESS != rval ) return 1;
 
 /*// set the edge tags on all elements sf1
 rval = set_edge_type_flag(mb, sf1); // form all edges, and set on them type 1 if constant
 latitude
 // add them to the set after this, just so we have them
 if (MB_SUCCESS != rval)
 return 1;*/
 
 IntxRllCssphere worker( mb );
 
 worker.set_error_tolerance( R * epsrel );
 // worker.SetEntityType(moab::MBQUAD);
 worker.set_radius( R );
 // worker.enable_debug();
 rval = worker.FindMaxEdges( sf1, sf2 );
 if( MB_SUCCESS != rval ) return 1;
 
 rval = worker.intersect_meshes( sf1, sf2, outputSet );
 // compute total area with 2 methods
 
 if( MB_SUCCESS != rval ) std::cout << " failed to intersect meshes\n";
 rval = mb->write_mesh( newFile, &outputSet, 1 );
 if( MB_SUCCESS != rval ) return 1;
 return 0;
}