add unref mesher

.gitlab-ci.yml

 makedg :
-  image : immc/dg-valid:v0
+  image : immc/dg-valid:v0.1
   script:
     - mkdir build
     - cd build
@@ -8,7 +8,7 @@ makedg :
     - ctest --output-on-failure
 
 win64 : 
-  image : immc/dg-mingw:v0.2.8
+  image : immc/dg-mingw:v0.2.11
   script :
     - mkdir build
     - cd build

CMakeLists.txt

@@ -197,7 +197,7 @@ add_subdirectory(tools)
 
 macro(dg_add_module_ module FILES)
   add_library(${module} SHARED ${FILES})
-  set_target_properties(${module} PROPERTIES COMPILE_FLAGS "-Wall -pedantic -DFORTIFY_SOURCE=1 -Wno-long-long -Wno-literal-suffix -Wno-unused-result" LINK_FLAGS "-Wall")
+  set_target_properties(${module} PROPERTIES COMPILE_FLAGS "-Wall -pedantic -DFORTIFY_SOURCE=1 -Wno-unknown-pragmas -Wno-long-long -Wno-literal-suffix -Wno-unused-result" LINK_FLAGS "-Wall")
   if(LAPACK_LINKER_FLAGS)
     set_target_properties(${module} PROPERTIES LINK_FLAGS ${LAPACK_LINKER_FLAGS})
   endif(LAPACK_LINKER_FLAGS)

benchmarks/bohai/mesh.py

+import slimPre
+R=6341e3
+
+bath_tiff = slimPre.fetch_ftp('slim_data/bohai/bohai_bathymetry.tif')
+coast_prj = slimPre.fetch_ftp('slim_data/bohai/no-intertidal/lines.prj')
+coast_dbf = slimPre.fetch_ftp('slim_data/bohai/no-intertidal/lines.dbf')
+coast_shx = slimPre.fetch_ftp('slim_data/bohai/no-intertidal/lines.shx')
+coast_shp = slimPre.fetch_ftp('slim_data/bohai/no-intertidal/lines.shp')
+
+bath = slimPre.geo_tiff(bath_tiff)
+dist = slimPre.distance_shp(coast_shp, 1000, [1,2,20])
+
+def fsize(x, y, z) :
+    return 10000/R
+
+physicals = {1:"coast", 2:"coast", 10:"yellow_sea_north", 11:"yellow_sea_center", 12:"yellow_sea_south", 20:"yellow_river"}
+
+slimPre.mesh_shp([coast_shp], 120 , 39, fsize, 500/R, "bohai.msh", "+proj=utm +ellps=WGS84 +zone=51", physicals)

dgpy/scripts/slimPre.py

@@ -971,3 +971,27 @@ def partition_nb():
 def exit(i):
     """Exit properly the program"""
     dgpy.Msg.Exit(i)
+
+def mesh_shp(shp_inputs, lon, lat, mesh_size, mesh_size_min, output, projection = "+latlong +ellps=WGS84",  boundary_names={}):
+    """Generate a mesh file in gmsh .msh file format (http://gmsh.info/). The mesh is generated on a sphere of radius 1. The mesh size
+    should be specified on this sphere (i.e. distances should be divided by the earth radius)
+
+    arguments:
+
+    * shp_inputs : list of string
+        vector of path to files in the ESRI shape file format. The shape files can be composed of point, lines or polygons. The projection
+        used in the shape files should be longitude-latitude (in degrees).
+    * lon, lat : float
+        position of one point inside the domain, the mesh region is the close surface containing this point.
+    * mesh_size : (float x, float y, float z) -> float
+        callback to get the desired element size, (on a sphere of radius 1) at one point on this sphere
+    * mesh_size_min : float
+        minimum value of the mesh size (approximatively the minimum value returned by mesh_size) on a sphere of radius 1
+    * output : string
+        name of the output file
+    * projection : string
+        description of the projection used in the output file in the libproj "+" format (the mesh itself is generated on the sphere)
+    * boundary_name : dictionary of (int, string)
+        to assign physicals names to the boundary elements of the mesh, a field (of type integer64) called "entity" should be present in the shape files. Then this parameter can be used to map the values of this field to physical names
+    """
+    dgpy.meshShp(shp_inputs, lon, lat, mesh_size, mesh_size_min, output, projection, boundary_names)

doc/docker-mingw/Dockerfile

-# docker build -t immc/dg-mingw:v0.2.8 .
+# docker build -t immc/dg-mingw:v0.2.11 .
 FROM ubuntu:16.04
 
 RUN apt-get update
@@ -20,6 +20,14 @@ RUN apt-get install -y m4 && \
  make -j8 && \
  make install
 
+RUN apt-get install -y m4 && \
+ wget http://dl.maptools.org/dl/libtiff/tiff-3.8.2.tar.gz && \
+ tar xf tiff-3.8.2.tar.gz && \
+ cd tiff-3.8.2 && \
+ ./configure --prefix=/mingw64 --host=x86_64-w64-mingw32 && \
+ make -j8 && \
+ make install
+
 # python-numpy-scipy-openblas
 WORKDIR /mingw64/src
 RUN wget http://launchpadlibrarian.net/279142124/p7zip-full_16.02+dfsg-1_amd64.deb && \
@@ -56,6 +64,9 @@ RUN tar xf petsc-lite-3.7.3.tar.gz && cd petsc-3.7.3 && \
   make install && \
   cp /mingw64/lib/libpetsc.so.3.7.3 /mingw64/Python34
 
+# missing file
+RUN cp /mingw64/bin/libtiff-3.dll /mingw64/Python34/
+
 # cmake toolchain file
 RUN echo "\ 
 SET(CMAKE_SYSTEM_NAME Windows)\n\
@@ -74,6 +85,7 @@ SET(PYTHON_EXECUTABLE /mingw64/Python34/python.exe)\n\
 SET(PYTHON_INCLUDE_DIR /mingw64/Python34/include/)\n\
 SET(PYTHON_LIBRARY /mingw64/lib/libpython3.4.a)\n" > /mingw64/cmake-mingw
 
+RUN apt-get install -y git
 RUN useradd -ms /bin/bash builder
 USER builder
 WORKDIR /home/builder

doc/docker-validator/Dockerfile

@@ -15,6 +15,7 @@ RUN /bin/bash -c 'cd petsc4py-3.6.0/; PETSC_DIR=/usr/lib/petscdir/3.6.2/x86_64-l
 
 RUN apt-get install -y gmsh
 RUN apt-get install -y python3-netcdf4
+RUN apt-get install -y libtiff-dev
 
 VOLUME ["/etc/gitlab-runner"]
 

modules/slimFunction/CMakeLists.txt

@@ -10,9 +10,11 @@ set(SRC
   slimTpxo.cpp
   slimFES.cpp
   slimExportNetCDF.cpp
+  unrefMesher.cpp
 )
 
 set(LIBS "")
+set(EXTERNAL_MODULE_DEPS "")
 
 option(ENABLE_NETCDF "Enable NETCDF" ON)
 if(ENABLE_NETCDF)
@@ -46,6 +48,7 @@ if(ENABLE_FES)
 endif(ENABLE_FES)
 
 option(ENABLE_PROJ "Enable PROJ" ON)
+option(ENABLE_UNREF "Enable UNREF mesh generator (requires ENABLE_PROJ)" ON)
 if(ENABLE_PROJ)
   if(WIN32)
   set (PROJ_LIB ${CMAKE_CURRENT_BINARY_DIR}/libproj/build/lib/libproj_4_9.a)
@@ -68,14 +71,31 @@ if(ENABLE_PROJ)
   message("slimFunction configured with PROJ")
   add_definitions("-DHAVE_PROJ")
   include_directories(${PROJ_INC})
+  list(APPEND EXTERNAL_MODULE_DEPS libproj)
   list(APPEND LIBS ${PROJ_LIB})
+  if (ENABLE_UNREF)
+    include(ExternalProject)
+    ExternalProject_Add(
+      unref 
+      GIT_REPOSITORY https://git.immc.ucl.ac.be/unref/unref.git
+      GIT_TAG 587ee41
+      INSTALL_COMMAND ""
+      CMAKE_ARGS -DCMAKE_BUILD_TYPE=${CMAKE_BUILD_TYPE} -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TOOLCHAIN_FILE} -DEXTERNAL_LIBPROJ=ON -DPROJ_INCLUDE_PATH=${PROJ_INC} -DPROJ_LIBRARY=${PROJ_LIB} -DBUILD_EXECUTABLE=OFF -DCMAKE_CXX_FLAGS=-fPIC
+      PREFIX unref
+      BINARY_DIR unref/build
+    )
+    list(APPEND EXTERNAL_MODULE_DEPS unref)
+    add_definition(-DHAVE_UNREF)
+    list(APPEND LIBS ${CMAKE_CURRENT_BINARY_DIR}/unref/build/libunref.a)
+    include_directories(${CMAKE_CURRENT_BINARY_DIR}/unref/src/unref)
+  endif(ENABLE_UNREF)
 endif(ENABLE_PROJ)
 
-dg_add_module(dgSlimFunction "${SRC}" "${LIBS}")
+dg_add_module(dgSlimFunction "${SRC}" "${LIBS};${PYTHON_LIBRARY}")
 dg_add_swig_module(slimFunction slimFunction.i dgSlimFunction)
 if(ENABLE_PROJ)
   set(SWIG_MODULE_proj_EXTRA_DEPS "libproj;dg")
   dg_add_swig_module(proj proj.i "${PROJ_LIB};dg")
-  add_dependencies(dgSlimFunction libproj)
+  add_dependencies(dgSlimFunction ${EXTERNAL_MODULE_DEPS})
 endif(ENABLE_PROJ)
 dg_add_test_directory(tests valentin.vallaeys@uclouvain.be)

modules/slimFunction/slimFunction.i

@@ -12,12 +12,14 @@
   #include "slimIterateHelper.h"
   #include "slimFES.h"
   #include "slimExportNetCDF.h"
+  #include "unrefMesher.h"
 %}
 
 %import(module="dgpy.dgCommon") "dgCommon.i"
 %import(module="dgpy.dgFunction") "functor.h"
 %import(module="dgpy.dgFunction") "function.h"
 
+%include "unrefMesher.h"
 %rename (_print) print;
 %include "ncDataDz.h"
 %include "slimFunction.h"

modules/slimFunction/unrefMesher.cpp

+#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

modules/slimFunction/unrefMesher.h

+#ifndef UNREF_H_
+#define UNREF_H_
+#ifdef HAVE_UNREF
+#include <string>
+#include <map>
+#include <vector>
+
+#include "Python.h"
+
+void meshShp(std::vector<std::string> shp_inputs, double lon, double lat, PyObject *callback, double lcmin, const char *filename, const char *projection=NULL, std::map<int, std::string> boundary_names = std::map<int,std::string>()) ;
+#endif
+#endif

tools/dgCommon.i

@@ -35,6 +35,7 @@ namespace std {
   %template(StringVector) vector<std::string, std::allocator<std::string> >;
   %template(MapStringVectorDouble) map<std::string, std::vector<double> >;
   %template(MapStringInt) map<std::string, int >;
+  %template(MapIntString) map<int, std::string >;
   %template(MapStringDouble) map<std::string, double>;
   %template(MapIntInt) std::map< int,int,std::less< int >,std::allocator< std::pair< int const,int > > >;
   %template(PairIntInt) pair<int, int>;