crop_vol_test.cpp

Go to the documentation of this file.

/**
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 *
 */
/*
 * crop_vol_test.cpp
 * \brief Driver to create a volume from 2 surfaces (top and bottom), a user defined polyline
 * (contour) and a direction
 *
 */
 
#include "moab/Core.hpp"
#include <iostream>
#include <fstream>
#include <cstring>
#include "moab/FBEngine.hpp"
#include "moab/GeomTopoTool.hpp"
#include "TestUtil.hpp"
 
std::string filename_top;
std::string filename_bot;
std::string polygon_file_name;
double min_dot = 0.8;
std::string vol_file;
 
int number_tests_successful = 0;
int number_tests_failed = 0;
 
using namespace moab;
 
ErrorCode volume_test( FBEngine* pFacet );
 
void handle_error_code( ErrorCode rv, int& number_failed, int& number_successful )
{
 if( rv == MB_SUCCESS )
 {
 std::cout << "Success";
 number_successful++;
 }
 else
 {
 std::cout << "Failure";
 number_failed++;
 }
}
 
int main( int argc, char* argv[] )
{
 filename_bot = TestDir + "unittest/BedCrop2.h5m";
 filename_top = TestDir + "unittest/SECrop2.h5m"; //"/polyPB.txt";
 polygon_file_name = TestDir + "unittest/poly14.txt";
 vol_file = "volIce.h5m"; // output
 
 number_tests_successful = 0;
 bool remove_output = true;
 
 if( argc == 1 )
 {
 std::cout << "Using default input files: " << filename_bot << " " << filename_top << " " << polygon_file_name
 << " " << min_dot << " " << vol_file << std::endl;
 std::cout << " default output file: " << vol_file << " will be deleted \n";
 }
 else if( argc == 6 )
 {
 remove_output = false;
 filename_bot = argv[1];
 filename_top = argv[2];
 polygon_file_name = argv[3];
 min_dot = atof( argv[4] );
 vol_file = argv[5];
 }
 else
 {
 std::cout << " wrong number of arguments\n";
 return 1; // fail
 }
 Core mbcore;
 Interface* mb = &mbcore;
 
 ErrorCode rval = mb->load_file( filename_bot.c_str() );MB_CHK_SET_ERR( rval, "failed to load bed file" );
 
 rval = mb->load_file( filename_top.c_str() ); // load the second face (we know we have one face in each file)
 MB_CHK_SET_ERR( rval, "failed to load top file" );
 
 FBEngine* pFacet = new FBEngine( mb, NULL, true ); // smooth facetting, no OBB tree passed
 
 if( !pFacet ) return 1; // error
 
 // should the init be part of constructor or not?
 // this is where the obb tree is constructed, and smooth faceting initialized, too.
 rval = pFacet->Init();MB_CHK_SET_ERR( rval, "failed to initialize smoothing" );
 
 std::cout << "volume creation test: ";
 rval = volume_test( pFacet );
 handle_error_code( rval, number_tests_failed, number_tests_successful );
 std::cout << "\n";
 
 if( remove_output )
 {
 remove( vol_file.c_str() );
 }
 return number_tests_failed;
}
 
ErrorCode volume_test( FBEngine* pFacet )
{
 // there should be 2 faces in the model
 std::ifstream datafile( polygon_file_name.c_str(), std::ifstream::in );
 if( !datafile )
 {
 std::cout << "can't read file\n";
 return MB_FAILURE;
 }
 //
 char temp[100];
 double direction[3]; // normalized
 double gridSize;
 datafile.getline( temp, 100 ); // first line
 
 // get direction and mesh size along polygon segments, from file
 sscanf( temp, " %lf %lf %lf %lf ", direction, direction + 1, direction + 2, &gridSize );
 // NORMALIZE(direction);// just to be sure
 
 std::vector< double > xyz;
 while( !datafile.eof() )
 {
 datafile.getline( temp, 100 );
 // int id = 0;
 double x, y, z;
 int nr = sscanf( temp, "%lf %lf %lf", &x, &y, &z );
 if( nr == 3 )
 {
 xyz.push_back( x );
 xyz.push_back( y );
 xyz.push_back( z );
 }
 }
 int sizePolygon = (int)xyz.size() / 3;
 if( sizePolygon < 3 )
 {
 std::cerr << " Not enough points in the polygon" << std::endl;
 return MB_FAILURE;
 }
 EntityHandle root_set;
 ErrorCode rval = pFacet->getRootSet( &root_set );MB_CHK_SET_ERR( rval, "ERROR : getRootSet failed!" );
 
 int top = 2; // iBase_FACE;
 
 Range faces;
 rval = pFacet->getEntities( root_set, top, faces );MB_CHK_SET_ERR( rval, "Failed to get faces in volume_test." );
 
 if( 2 != faces.size() )
 {
 std::cerr << "expected 2 faces, top and bottom\n";
 return MB_FAILURE;
 }
 EntityHandle newFace1; // first test is with closed surface
 rval = pFacet->split_surface_with_direction( faces[0], xyz, direction, /*closed*/ 1, min_dot, newFace1 );MB_CHK_SET_ERR( rval, "Failed to crop first face." );
 
 EntityHandle newFace2; // first test is with closed surface
 rval = pFacet->split_surface_with_direction( faces[1], xyz, direction, /*closed*/ 1, min_dot, newFace2 );MB_CHK_SET_ERR( rval, "Failed to crop second face." );
 
 // here we make the assumption that the edges, vertices are matching fine...
 EntityHandle volume;
 rval = pFacet->create_volume_with_direction( newFace1, newFace2, direction, volume );MB_CHK_SET_ERR( rval, "Failed to create volume." );
 
 Interface* mb = pFacet->moab_instance();
 pFacet->delete_smooth_tags();
 // duplicate -- extract a new geom topo tool
 GeomTopoTool* gtt = pFacet->get_gtt();
 GeomTopoTool* duplicate = NULL;
 std::vector< EntityHandle > gents;
 gents.push_back( volume );
 rval = gtt->duplicate_model( duplicate, &gents );MB_CHK_SET_ERR( rval, "Failed to extract volume." );
 EntityHandle newRootSet = duplicate->get_root_model_set();
 delete pFacet;
 pFacet = NULL; // try not to write the obb tree
 delete duplicate;
 rval = mb->write_file( vol_file.c_str(), NULL, NULL, &newRootSet, 1 );
 
 return rval;
}