unrefMesher.cpp 2.02 KB
Newer Older
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
#ifdef HAVE_UNREF
#include <cmath>
#include "unref.h"
#include "proj_api.h"
#include "unrefMesher.h"
#include "dgMessage.h"

static double evalPythonCallback(PyObject *callback, double x, double y, double z) {
  PyObject *args = Py_BuildValue("(ddd)", x, y, z);
  PyObject *result = PyEval_CallObject(callback, args);
  double r = 1;
  if (result) {
    if (!PyNumber_Check(result))
      Msg::Fatal("python callback does not return a floating point value\n");
    r = PyFloat_AsDouble(result);
    Py_DECREF(result);
  }
  else{
    PyErr_Print();
    Msg::Fatal("error in python callback\n");
  }
  Py_DECREF(args);
  return r;
}

static Vertex project(projPJ pj_in, projPJ pj_out, const Vertex &v) {
  Vertex vout(0, v.x(), v.y(), v.z(), 0.);
  pj_transform(pj_in, pj_out, 1, 1, &vout.x(), &vout.y(), &vout.z());
  return vout;
}

void meshShp(std::vector<std::string> shp_inputs, double lon, double lat, PyObject *callback, double lcmin, const char *filename, const char *projection,  std::map<int, std::string> boundary_names) {
  int NUM_THREADS=1;
  Py_INCREF(callback);
  std::map<int, std::string> physicals[4];
  physicals[1] = boundary_names;
  std::vector<unref::TaggedVertex*> vertices;
  for (auto & filename : shp_inputs){
    auto vs = unref::readSHP(filename.c_str(), lcmin);
    vertices.insert(vertices.end(), vs.begin(), vs.end());
  }
  auto f = std::bind(&evalPythonCallback, callback, std::placeholders::_1, std::placeholders::_2, std::placeholders::_3);
  auto lineloops = unref::genBoundaries(vertices, lon*M_PI/180, lat*M_PI/180, f, NUM_THREADS);
  unref::mesh mesh(lineloops, f, lcmin, NUM_THREADS);
  if(projection) {
    projPJ pj_ll = pj_init_plus("+proj=geocent +ellps=sphere +a=1 +b=1");
    projPJ pj_out = pj_init_plus(projection);
    auto fp = std::bind(&project, pj_ll, pj_out, std::placeholders::_1);
    mesh.writeMsh(filename, physicals, fp, NUM_THREADS);
    pj_free(pj_ll);
    pj_free(pj_out);
  }
  else {
    mesh.writeMsh(filename, physicals, [](const Vertex &x){return x;}, NUM_THREADS);
  }
  Py_DECREF(callback);
}
#endif