/** @example HelloMOAB.cpp
 * Description: read a mesh, get the entities.\n
 * HelloMOAB is a simple test file which is used to read meshes from VTK file and test how many entities there are.\n
 *
 * To run: ./HelloMOAB [meshfile]\n
 * (default values can run if users don't specify a mesh file)
 */


#include "moab/Core.hpp"
#include <iostream>

using namespace moab;
using namespace std;

#ifndef MESH_DIR
#define MESH_DIR "."
#endif

// Note: change the file name below to test a trivial "No such file or directory" error
string test_file_name = string(MESH_DIR) + string("/64bricks_12ktet.h5m");

int main(int argc, char **argv)
{
  // Get MOAB instance
  Interface* mb = new (std::nothrow) Core;
  if (NULL == mb)
    return 1;

  // Need option handling here for input filename
  if (argc > 1) {
    // User has input a mesh file
    test_file_name = argv[1];
  }

  // Load the mesh from file
  ErrorCode rval = mb->load_mesh(test_file_name.c_str());MB_CHK_ERR(rval);

  // Get verts entities, by type
  Range verts;
  rval = mb->get_entities_by_type(0, MBVERTEX, verts);MB_CHK_ERR(rval);

  // Get edge entities, by type
  Range edges;
  rval = mb->get_entities_by_type(0, MBEDGE, edges);MB_CHK_ERR(rval);

  // Get faces, by dimension, so we stay generic to entity type
  Range faces;
  rval = mb->get_entities_by_dimension(0, 2, faces);MB_CHK_ERR(rval);

  // Get regions, by dimension, so we stay generic to entity type
  Range elems;
  rval = mb->get_entities_by_dimension(0, 3, elems);MB_CHK_ERR(rval);

  // Output the number of entities
  cout << "Number of vertices is " << verts.size() << endl;
  cout << "Number of edges is " << edges.size() << endl;
  cout << "Number of faces is " << faces.size() << endl;
  cout << "Number of elements is " << elems.size() << endl;

  // get and create all faces adjacent to cells
  faces.clear();
  rval = mb->get_adjacencies(elems, 2, true, faces, Interface::UNION); MB_CHK_ERR(rval);
  cout << "Number of faces adjacent to cells: " << faces.size() << endl;

  // get and create all edges adjacent to cells
  edges.clear();
  rval = mb->get_adjacencies(elems, 1, true, edges, Interface::UNION); MB_CHK_ERR(rval);
  cout << "Number of edges adjacent to cells: " << edges.size() << endl;

  // put in a new set all faces
  EntityHandle faces_set;
  rval = mb->create_meshset(MESHSET_SET, faces_set);MB_CHK_ERR(rval);
  rval = mb->add_entities(faces_set, faces);MB_CHK_ERR(rval);

  // put in a new set all edges
  EntityHandle edges_set;
  rval = mb->create_meshset(MESHSET_SET, edges_set);MB_CHK_ERR(rval);
  rval = mb->add_entities(edges_set, edges);MB_CHK_ERR(rval);

  // save a new file with all faces and a new file with all edges
  rval = mb->write_file("faces.vtk", 0, 0, &faces_set, 1);MB_CHK_ERR(rval);

  rval = mb->write_file("edges.vtk", 0, 0, &edges_set, 1);MB_CHK_ERR(rval);

  // get the first cell
  if (!elems.empty())
  {
    EntityHandle firstCell=elems[0];
    // get faces adjacent to it
    std::vector<EntityHandle> fs;
    rval = mb->get_adjacencies(&firstCell, 1, 2, false, fs); MB_CHK_ERR(rval);
    cout<<" faces adjacent to first cell " << fs.size() << endl;

    mb->list_entity(firstCell);

    const EntityHandle * conn = NULL;
    int numV=0;
    rval = mb->get_connectivity(firstCell, conn , numV);  MB_CHK_ERR(rval);
    std::vector<double> coords;
    coords.resize(3*numV);
    rval = mb->get_coords(conn, numV, &coords[0]);
    for (int i=0; i<numV; i++)
    {
      cout <<" vertex id " << mb->id_from_handle(conn[i]) << " xyz: " << coords[3*i] << " "
           << coords[3*i+1] << " " << coords[3*i+2] << endl;
    }

  }
  delete mb;

  return 0;
}