#include "moab/Core.hpp"
#include "moab/Range.hpp"
#include "moab/ParallelComm.hpp"
#include "moab/Skinner.hpp"
#include "moab/CN.hpp"
#include "MBTagConventions.hpp"
//#include "parallel.h"

#include <iostream>
#include <iterator>
#include <list>
#include <algorithm>
#include <sstream>
#include <limits>
#include <string>
#include <time.h>

#include "moab_mpi.h"

const char PARTTAG[] = "PARALLEL_PARTITION";
const char READ_DEL_OPT[] = "READ_DELETE";
const char BCAST_DEL_OPT[] = "BCAST_DELETE";
const char READ_PART_OPT[] = "READ_PART";
const char* const DEFAULT_MODE = READ_PART_OPT;

const char USAGE[] = " [-v <n>] [-R] [-p <parttag>[=val]] [-D|-B|-P] <input_file>";

static void usage( const char* argv0 ) {
  std::cerr << "Usage: " << argv0 << USAGE << std::endl
            << "       " << argv0 << " -h" << std::endl;
  exit(1);
}

static const char* defstr( const char* mode )
{
  static const char s[] = " (default)";
  int len = (mode == DEFAULT_MODE) ? 0 : sizeof(s);
  return s + len;
}    

static void help( const char* argv0 ) {
  std::cout << "Usage: " << argv0 << USAGE << std::endl
            << "  -v  : debug output level (DEBUG_IO option)" << std::endl
            << "  -R  : do not resolve shared entities" << std::endl
            << "  -p  : partition tag, with optional '=' value" << std::endl
            << "          (default = \"" << PARTTAG << "\")" << std::endl
            << "  -D  : read mode as \"" << READ_DEL_OPT << "\"" << defstr(READ_DEL_OPT) << std::endl
            << "  -B  : read mode as \"" << BCAST_DEL_OPT << "\"" << defstr(READ_DEL_OPT) << std::endl
            << "  -P  : read mode as \"" << READ_PART_OPT << "\"" << defstr(READ_DEL_OPT) << std::endl
            ;
  exit(0);
}

using namespace moab;

int check_parallel_read( Interface& mb, ParallelComm& pcomm, bool shared_ents );

int main( int argc, char* argv[] )
{
  const char* read_mode_opt = DEFAULT_MODE;
  const char* part_tag_opt = PARTTAG;
  bool resolve_shared = true;
  const char* input_file = 0;
  const char* debug_flag_str = 0;
  bool no_more_flags = false;

  MPI_Init( &argc, &argv );
  
  for (int i = 1; i < argc; ++i) {
    if (!no_more_flags && argv[i][0] == '-') {
      const char* opts = argv[i]+1;
      for (int j = 0; opts[j]; ++j) {
        switch (opts[j]) {
          case '-':
            no_more_flags = true;
            break;
          case 'v':
            if (++i == argc)
              usage(argv[0]);
            debug_flag_str = argv[i];
            break;
          case 'R':
            resolve_shared = false;
            break;
          case 'p':
            if (++i == argc)
              usage(argv[0]);
            part_tag_opt = argv[i];
            break;
          case 'D':
            read_mode_opt = READ_DEL_OPT;
            break;
          case 'B':
            read_mode_opt = BCAST_DEL_OPT;
            break;
          case 'P':
            read_mode_opt = READ_PART_OPT;
            break;
          case 'h':
            help(argv[0]);
            break;
          default:
            usage(argv[0]);
        }
      }
    }
    else if (!input_file) {
      input_file = argv[i];
    }
    else {
      usage(argv[0]);
    }
  }
  
  if (!input_file) {
    std::cerr << argv[0] << ": no input file specified" << std::endl;
    usage(argv[0]);
  }
  
  
    // build options string
  
  std::ostringstream opts;
  opts << "PARALLEL=" << read_mode_opt;
  
  std::string part_opt( part_tag_opt );
  size_t p = part_opt.find_last_of('=');
  if (p == std::string::npos) {
    opts << ";PARTITION=" << part_opt;
  }
  else {
    char* endptr = 0;
    long n = strtol( part_opt.c_str() + p + 1, &endptr, 0 );
    if (*endptr || p == part_opt.size()-1) {
      std::cerr << "Warning: partition tag option contains an '=' followed "
                   "         by a non-integer value.  Assuming tag name is "
                   "         \"" << part_opt << "\"" << std::endl;
      opts << ";PARTITION=" << part_opt;
    }
    else {
      opts << ";PARTITION=" << part_opt.substr( 0, p );
      opts << ";PARTITION_VAL=" << n;
    }
  }
  
  if (resolve_shared) {
    opts << ";PARALLEL_RESOLVE_SHARED_ENTS";
  }
  
  if (debug_flag_str) {
    char* endptr = 0;
    long n = strtol( debug_flag_str, &endptr, 0 );
    if (*endptr || n < 0 || !*debug_flag_str)
      usage(argv[0]);
    opts << ";DEBUG_IO=" << n;
  }
  
  Core moab;
  Interface& mb = moab;
  ParallelComm pcomm( &mb, MPI_COMM_WORLD );
  if (pcomm.rank() == 0) 
    std::cout << "Loading file: \"" << input_file << "\" with options \"" << opts.str() << '"' << std::endl;
  
  ErrorCode rval = mb.load_file( input_file, 0, opts.str().c_str() );
  if (MB_SUCCESS != rval) {
    std::cerr << input_file << " : file read failed (" << mb.get_error_string(rval) << ")" << std::endl;
    return 1;
  }
  
  int result1 = 0;
  result1 = check_parallel_read( mb, pcomm, resolve_shared );
  
  char fname[1024];
  sprintf(fname,"pparts1_%01d.vtk",pcomm.rank());
  rval = mb.write_file(fname);

  MPI_Finalize();
  return result1;
}

int check_parallel_read( Interface& mb, ParallelComm& pcomm, bool /*shared_ents*/ )
{
  int error_count = 0;

  ErrorCode rval = pcomm.resolve_shared_ents( 0, 3); 
  
  rval = pcomm.exchange_ghost_cells(3,2,1,0,true,true);

  pcomm.assign_global_ids(0, 3, 0, true, true, false);

  Range faces;
  Range ghosts;
  if (pcomm.rank() == 0)
  {
    pcomm.get_iface_entities(1,2,faces);
    int ifaces_size = faces.size(); 
    std::cout<<"iface faces size "<<ifaces_size<<std::endl;
    int numfaces0;
    mb.get_number_entities_by_dimension(0,2,numfaces0);
    rval = mb.get_adjacencies(faces, 3, false, ghosts, Interface::UNION);
    faces.clear();
    rval = mb.get_adjacencies(ghosts, 0, true, faces, Interface::UNION);
    int numfaces1;
    mb.get_number_entities_by_dimension(0,2,numfaces1);
    std::cout<<"ghost cells size "<<ghosts.size();
    unsigned char status;
    int count = 0;
    for (unsigned int i = 0; i < ghosts.size(); i++)
    {
      rval = pcomm.get_pstatus(ghosts[i], status);
      if (status == PSTATUS_SHARED)
	count++;
    }
    std::cout<<" shared = "<<count<<std::endl;
  
    std::vector<unsigned char> pstatus;
    Range::iterator rit;
    std::vector<EntityHandle> dum_range_ghost;
    ErrorCode result;
    count=0;
    for (rit = ghosts.begin(); rit != ghosts.end(); rit++) 
    {
      EntityHandle elem = *rit;
      pstatus.resize(1);
      rval = mb.tag_get_data(pcomm.pstatus_tag(), &elem, 1, &pstatus[0]);
      if (pstatus[0] != PSTATUS_SHARED && pstatus[0] != PSTATUS_NOT_OWNED && pstatus[0] != PSTATUS_GHOST) 
      {
        count++;	
        short int num_sub_ent = CN::NumSubEntities(mb.type_from_handle(elem), 2);
        for (unsigned int side_number = 0; side_number < num_sub_ent; side_number++)
        {
          EntityHandle neighbor = 0;
          std::vector<EntityHandle> conn;
          result = mb.get_connectivity(&elem, 1, conn);
          EntityHandle bridge_verts[MAX_SUB_ENTITIES];
          int bridge_indices[MAX_SUB_ENTITIES];

          Range vertex_neighbors;
          int num_sub_verts = CN::VerticesPerEntity( CN::SubEntityType(mb.type_from_handle(elem), 2, side_number ) );
          CN::SubEntityVertexIndices(mb.type_from_handle(elem), 2, side_number, bridge_indices );
          for (int j = 0; j < num_sub_verts; ++j)
            bridge_verts[j]= conn[bridge_indices[j]];

          result = mb.get_adjacencies(&bridge_verts[0], num_sub_verts, 3, true, vertex_neighbors);

          if (vertex_neighbors.size()==1)
          {
             neighbor=elem;
             if (vertex_neighbors[0] != (elem))
             {
                std::cout<<"error"<<std::endl;
                return MB_MULTIPLE_ENTITIES_FOUND;
             }
          }
          else if (vertex_neighbors.size()==2)
          {
            moab::Range::iterator iter;
            for(iter = vertex_neighbors.begin(); iter != vertex_neighbors.end(); ++iter)
            {
              if ((*iter) != (elem))
              {
                 neighbor=(*iter);
              }
            }
          }
          else if (vertex_neighbors.size() == 0)
          {
	    std::cout<<"error"<<std::endl;
          }
          else if (vertex_neighbors.size() > 2)
          {
	    std::cout<<"error"<<std::endl;
          }
          if (neighbor==elem) 
          {
             dum_range_ghost.push_back(neighbor);
          }
        }
      }
    }
    std::cout<<"owned entities = "<<count<<" adjacent entities not found = "<<dum_range_ghost.size()<<" the number of adjacencies not found should be "<<ifaces_size<<" smaller"<<std::endl;
  }

  return error_count;
}