dual_test.cpp

Go to the documentation of this file.

// tests dual construction code
 
#include "moab/Core.hpp"
#include "moab/Range.hpp"
#include "moab/MeshTopoUtil.hpp"
#include "moab/DualTool.hpp"
#include "../TestUtil.hpp"
#include <iostream>
#include <string>
 
using namespace moab;
 
std::string default_input = TestDir + "unittest/hex01.vtk";
const char default_output[] = "dual.vtk";
 
Interface* gMB;
 
int main( int argc, char* argv[] )
{
 const char* f1 = default_input.c_str();
 const char* f2 = default_output;
 if( argc <= 2 )
 {
 std::cout << "Usage: dual_test <mesh_file_name> <out file>" << std::endl;
 std::cout << "using default input : " << default_input << " and output " << default_output << "\n";
 }
 else
 {
 f1 = argv[1];
 f2 = argv[2];
 }
 
 gMB = new Core();
 
 // read the mesh file
 ErrorCode result = gMB->load_mesh( f1 );
 if( MB_SUCCESS != result )
 {
 std::cout << "Problems reading file " << argv[1] << "." << std::endl;
 delete gMB;
 return 1;
 }
 
 // make sure aentities are created
 Range all_verts;
 result = gMB->get_entities_by_dimension( 0, 0, all_verts );
 if( MB_SUCCESS != result ) std::cout << "Problem getting vertices." << std::endl;
 
 MeshTopoUtil( gMB ).construct_aentities( all_verts );
 
 // get counts before constructing dual
 int num_edges;
 result = gMB->get_number_entities_by_dimension( 0, 1, num_edges );
 if( MB_SUCCESS != result ) std::cout << "Problem getting number of edges." << std::endl;
 
 // see if it's all-hex or all-quad, and construct hex/quad dual if so
 int num_hex, num_quad, num_3d, num_2d;
 result = gMB->get_number_entities_by_dimension( 0, 2, num_2d );
 if( MB_SUCCESS != result ) return 0;
 result = gMB->get_number_entities_by_dimension( 0, 3, num_3d );
 if( MB_SUCCESS != result ) return 0;
 result = gMB->get_number_entities_by_type( 0, MBQUAD, num_quad );
 if( MB_SUCCESS != result ) return 0;
 result = gMB->get_number_entities_by_type( 0, MBHEX, num_hex );
 if( MB_SUCCESS != result ) return 0;
 
 // construct the dual for this mesh
 DualTool dt( gMB );
 if( num_quad == num_2d && num_hex == num_3d )
 result = dt.construct_hex_dual( 0, 0 );
 else
 result = dt.construct_dual( 0, 0 );
 if( MB_SUCCESS != result )
 {
 std::cout << "Problems constructing dual." << std::endl;
 return 0;
 }
 
 // print information about the dual
 Range dual_cells, dual_faces;
 result = dt.get_dual_entities( 0, 0, 2, dual_faces );
 if( MB_SUCCESS != result )
 std::cout << "Problem getting dual faces." << std::endl;
 else
 std::cout << "Found " << dual_faces.size() << "/" << num_edges << " dual faces." << std::endl;
 
 result = dt.get_dual_entities( 0, 0, 3, dual_cells );
 if( MB_SUCCESS != result )
 std::cout << "Problem getting dual cells." << std::endl;
 else
 std::cout << "Found " << dual_cells.size() << "/" << all_verts.size() << " dual cells." << std::endl;
 
 // print information about dual hyperplanes, if any
 Tag hp_tag;
 Range hp_sets;
 if( num_2d == num_quad )
 {
 // get sets with the right tag
 result = gMB->tag_get_handle( DualTool::DUAL_CURVE_TAG_NAME, 1, MB_TYPE_HANDLE, hp_tag );
 if( MB_SUCCESS == result )
 {
 result = gMB->get_entities_by_type_and_tag( 0, MBENTITYSET, &hp_tag, NULL, 1, hp_sets );
 if( MB_SUCCESS == result && !hp_sets.empty() )
 std::cout << "Found " << hp_sets.size() << " 1d dual hyperplanes (chords)." << std::endl;
 }
 }
 
 if( num_3d == num_hex )
 {
 // get sets with the right tag
 result = gMB->tag_get_handle( DualTool::DUAL_SURFACE_TAG_NAME, 1, MB_TYPE_HANDLE, hp_tag );
 if( MB_SUCCESS == result )
 {
 hp_sets.clear();
 result = gMB->get_entities_by_type_and_tag( 0, MBENTITYSET, &hp_tag, NULL, 1, hp_sets );
 if( MB_SUCCESS == result && !hp_sets.empty() )
 std::cout << "Found " << hp_sets.size() << " 2d dual hyperplanes (sheets)." << std::endl;
 }
 }
 
 // write GMV file
 gMB->write_file( f2 );
 delete gMB;
 
 return 0;
}