#include <assert.h>
#include <math.h>
#include <sstream>
#include <string>
#include <moab/Skinner.hpp>
#include <moab/Range.hpp>
#include <moab/Core.hpp>
#include <moab/MeshTopoUtil.hpp>
#include <X11/X.h>
#include <iostream>
#include <vector>
#include <list>
#include <stdio.h>
#include <stdlib.h>
#include <ostream>

/** Macro used to return if MOAB operation failed 
 */
#define RINS(imesh, ierr) \
do { \
  if (moab::MB_SUCCESS != ierr)  {\
    std::cerr << "MOAB failure: " \
              << (imesh)->get_error_string(ierr) << "\n"; \
    abort(); \
  }\
} while(0)


double entity_length(moab::Interface *imesh, const moab::Range &entities);

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

  moab::Core mesh;
  moab::ErrorCode ierr = moab::MB_FAILURE;
  char *filename = NULL;
  if (argc > 1) {
    filename = argv[1];
    ierr = mesh.load_file(filename); 
    if (moab::MB_SUCCESS != ierr) {
      std::cout << "Cannot read file " << filename << "\n";
    }
  } else {
    std::cout << "Usage: " << argv[0] << " FILE\n"; 
    return 1;
  }
 
  moab::Range elements;
  moab::Range edges;
  moab::Skinner skinner(&mesh);

  ierr = mesh.get_entities_by_dimension(0, 2, elements);
  RINS(&mesh, ierr);
  ierr =  skinner.find_skin(elements, 1, edges);
  RINS(&mesh, ierr);

  double length = entity_length(&mesh, edges);

  std::cout << "Perim length: " << length << "\n";
  return EXIT_SUCCESS;
}

static 
double edge_length(moab::Interface *imesh, int geodim, 
                   moab::EntityHandle handle) {
  int dim = imesh->dimension_from_handle(handle);  
  if (1 != dim) {
    std::cout << "Entity is not edge\n";
    abort();
  }
  moab::ErrorCode ierr;
  moab::Range vertices;
  ierr = imesh->get_connectivity(&handle, 1, vertices);
  RINS(imesh, ierr);
  double coords[6];
  ierr = imesh->get_coords(vertices, coords);

  std::cout << "V: " << vertices[0] << " " << vertices[1] << "\n";
  for (int i =0; i<6; i++) {
    std::cout << coords[i] << " " ;
  }
  std::cout << "\n";

  RINS(imesh, ierr);
  double length=0.0;
  for (int j=0; j<geodim; j++) {
     double d;
     d = coords[j] - coords[j+geodim];
     length += d*d;
  }
  return sqrt(length);
}

double entity_length(moab::Interface *imesh, const moab::Range &entities) {
  double length=0;
  moab::ErrorCode ierr;
  int geodim=-1;
  ierr = imesh->get_dimension(geodim);
  RINS(imesh, ierr);
  for (moab::Range::const_iterator ei=entities.begin(); 
                                    ei != entities.end();
                                   ++ei) {
    moab::EntityHandle handle = *ei;
    length += edge_length(imesh, geodim, handle);
  }
  return length;
}