#include "moab/Core.hpp"
#include "moab/Interface.hpp"
#include "moab/Range.hpp"
#include "moab/ParallelComm.hpp"
#include "MBParallelConventions.h"

/* ========================================================== */
#define PRINT_LAST_ERROR \
    if (MB_SUCCESS != result) {\
      std::string tmp_str;\
      mb->get_last_error(tmp_str);\
      std::cout << "Failure; message:" << tmp_str << std::endl;\
      delete mb;\
      return result;\
    }
    
/* ========================================================== */

using namespace moab;
using namespace std;


int main(int argc, char **argv) {

  EntityHandle mesh;
  EntityHandle root;

  MPI_Init(&argc, &argv);

  Interface* mb = new Core;
  if (NULL == mb)
    return 1;
  
  ErrorCode result = MB_SUCCESS;

  string options;
  string mesh_file_name;
  
  mesh_file_name = argv[1];

  options = "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARALLEL_RESOLVE_SHARED_ENTS";

  ParallelComm* pcomm = new ParallelComm(mb, MPI_COMM_WORLD);
  int nprocs = pcomm->proc_config().proc_size();
  int rank = pcomm->proc_config().proc_rank();
  MPI_Comm comm = pcomm->proc_config().proc_comm();


  result = mb->load_file(mesh_file_name.c_str(), 0, options.c_str());
  PRINT_LAST_ERROR; 
  
  result = pcomm->exchange_ghost_cells(3, // int ghost_dim,
                                     2, // int bridge_dim,
                                     1, //int num_layers,
                                     2, //int addl_ents,
                                     true); // bool store_remote_handles);
  PRINT_LAST_ERROR;
     
  
  Range ents;
  result = mb->get_entities_by_dimension(0, 3, ents);
  PRINT_LAST_ERROR;
  
  Range adjs;
  result = mb->get_adjacencies(ents, 2, true, adjs, Interface::UNION);
  PRINT_LAST_ERROR;

  result = mb->create_meshset(0,mesh,0);
  result = mb->add_entities(mesh,ents);
  result = mb->add_entities(mesh,adjs);

  root = mb->get_root_set();

  result = mb->add_entities(mesh,&root,1);
  
  Range edges;
  result = mb->get_adjacencies(ents, 1, true, edges, Interface::UNION);
  PRINT_LAST_ERROR;  
  result = mb->delete_entities(edges);
  
  std::string writeOpts = "PARALLEL=WRITE_PART";
  std::string outFile = "out.h5m";
  result = mb->write_file(outFile.c_str(), NULL, writeOpts.c_str(),&mesh);
  PRINT_LAST_ERROR;
    
	if (rank == 0){
	 std::cout << "Wrote " << outFile << std::endl;
	}
    
 delete mb;

 return 0;
 }