removing get_distance from scipy. features to read the distance to a groupOfFaces

dgpy/scripts/slimPre.py

@@ -6,7 +6,7 @@
     SLIM solves the governing equations on adaptative unstructured meshes using the discontinuous Galerkin finite element method.
 
     slimPre contains the tools writing the input files necessary for the simulations (for simulations, see slim.py file), such as partitioned mesh, bathymetry, coriolis, wind, tides, ...
-    
+
     To generate the input files, run the script my_prepro_script.py:
         - On a single core: rundgpy my_prepro_script.py
         - On multiple cores: mpirun -np N rundgpy my_prepro_script.py (replace N by the number of cores)
@@ -28,7 +28,7 @@
         The write_file function writes a netcdf file, using a Region object (giving the nodes for which there is data), a time object (giving the times of the data), and a list of data tuples.
         The first variable of a data tuple is the data name.
         The secod variable of a data tuple is the data array, with as many lines as the number of nodes, and as many columns as the number of time steps
-    
+
     netcdf_to_msh:
         The function transforms a netcdf file to a .msh file.
 
@@ -40,13 +40,13 @@
 
     smooth_bathymetry:
         The function smooths the given bathymetry, following a smoothing criterion, and writes it to a file. It can also write a bathymetry constant per element (P0), for a better stability. If the bathymetry is P0, it will be written in a idx file. Otherwise, it will be written in a netcdf file.
-    
+
     interpolate_from_structured_grid
         The function interpolates from structured data 
-    
+
     interpolate_from_structured_file
         The function interpolates structured data from a GMSH-file
-    
+
     get_data_from_netcdf:
         The function returns the data array contained in the netcdf file given as an argument
 
@@ -61,16 +61,16 @@
 
     make_directory:
         This function makes a directory if it does not exist
-        
+
     partition_mesh:
         This function partitions the mesh (if run in multiprocessing) 
-    
+
     partition_id:
         This function returns the number of the current partition as a string
-    
+
     partition_nb:
         This function returns the number of partitions used as a string
-    
+
     exit:
         This function closes the program
 """
@@ -840,9 +840,9 @@ def get_distance(region, physical_tags):
     * physical_tags
         vector containing the physical tags from which the minimum distance will be computed
     """
-    xyz = region.coordinates
-    dist = slim_private.distance.distance(region._mesh, physical_tags, evaluationPoints=xyz)
-    return np.array(dist)
+    dist = distance_msh(region, 1000, physical_tags)
+    proj = Coordinate_system(region, data_proj="+proj=geocent +ellps=sphere +a=1 +b=1")
+    return dist.evaluate_matrix(proj.coordinates)
 
 def extrude(mesh_file_name, bath_file_name, nb_layers=None, z_layers=None, layers_function=None, mesh_file_name_out='', factor_show=0, periodicity=None):
 
@@ -1012,6 +1012,20 @@ def mesh_shp(shp_inputs, lon, lat, mesh_size, mesh_size_min, output, projection
     """
     dgpy.meshShp(shp_inputs, lon, lat, mesh_size, mesh_size_min, output, projection, boundary_names)
 
+class distance_msh(slim_private.pre_evaluator) :
+    """ compute the distance from an interface region."""
+    def __init__(self, region, dist) :
+        """
+        arguments:
+        * region : slim region 
+            (only implemented for interfaces region)
+        * dmax : float
+             points on line will be inserted every dmax interval
+        """
+        if region._type is not "interface":
+            dgpy.Msg.Fatal('distance_msh is only implemented for region of interfaces')
+        super(distance_msh, self).__init__(dgpy.slimPreEvaluatorDistance(region._groups, region._mesh._mesh_proj, dist))
+
 class distance_shp(slim_private.pre_evaluator) :
     """ compute the distance from a set of points given in an esri shape file (shp). Distances are evaluated on an ellipsoid."""
     def __init__(self, shp_filename,  dmax, tags=[], ellps="+ellps=WGSG4") :
@@ -1026,7 +1040,7 @@ class distance_shp(slim_private.pre_evaluator) :
         * ellps : string (default: "+ellps=WGS84")
             libproj string describing the ellipsoid
         """
-        super(distance_shp, self).__init__(dgpy.slimPreEvaluatorDistanceSHP(shp_filename, dmax, tags))
+        super(distance_shp, self).__init__(dgpy.slimPreEvaluatorDistance(shp_filename, dmax, tags))
 
 class geo_tiff(slim_private.pre_evaluator) :
     """ extract data from a geo TIFF file """

mesh/mesh.i

@@ -59,6 +59,7 @@
 %include "dgMeshPartition.h"
 
 %template(VectorGroupCollection) std::vector<const dgGroupCollection*>;
+%template(FaceGroupVector) std::vector<dgGroupOfFaces*>;
 
 %typemap(typecheck) string = char *;
 %typemap(typecheck) const string & = char *;

modules/slimFunction/CMakeLists.txt

@@ -12,7 +12,7 @@ set(SRC
   slimExportNetCDF.cpp
   unrefMesher.cpp
   slimPreEvaluator.cpp
-  slimPreEvaluatorDistanceSHP.cpp
+  slimPreEvaluatorDistance.cpp
   slimPreEvaluatorGeoTIFF.cpp
 )
 
@@ -83,7 +83,7 @@ if(ENABLE_PROJ)
     ExternalProject_Add(
       unref 
       GIT_REPOSITORY https://git.immc.ucl.ac.be/unref/unref.git
-      GIT_TAG 2d5baab40e344cd8c2677ec05b2f3763dece6f3b
+      GIT_TAG a6b4ae1751c44082692ac350040d31907457a669
       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

modules/slimFunction/slimFunction.i

@@ -14,13 +14,15 @@
   #include "slimExportNetCDF.h"
   #include "unrefMesher.h"
   #include "slimPreEvaluator.h"
-  #include "slimPreEvaluatorDistanceSHP.h"
+  #include "slimPreEvaluatorDistance.h"
   #include "slimPreEvaluatorGeoTIFF.h"
 %}
 
 %import(module="dgpy.dgCommon") "dgCommon.i"
 %import(module="dgpy.dgFunction") "functor.h"
 %import(module="dgpy.dgFunction") "function.h"
+//%import(module="dgpy.dgMesh") "mesh.i"
+
 
 %include "slimFunctionConfig.h"
 %include "unrefMesher.h"
@@ -36,5 +38,5 @@
 %include "slimFES.h"
 %include "slimExportNetCDF.h"
 %include "slimPreEvaluator.h"
-%include "slimPreEvaluatorDistanceSHP.h"
+%include "slimPreEvaluatorDistance.h"
 %include "slimPreEvaluatorGeoTIFF.h"

modules/slimFunction/slimPreEvaluatorDistanceSHP.cpp → modules/slimFunction/slimPreEvaluatorDistance.cpp

 #include "slimFunctionConfig.h"
 #if defined(HAVE_UNREF) and defined(HAVE_FLANN)
-#include "slimPreEvaluatorDistanceSHP.h"
+#include "slimPreEvaluatorDistance.h"
+#include "proj_api.h"
+#include "dgGroupOfElements.h"
+#include "dgDofContainer.h"
 #include <cmath> // required for flann
 #include "flann/flann.hpp"
 namespace flann {namespace serialization {BASIC_TYPE_SERIALIZER(unsigned long long);}} // hack to build build on mingw_w64 x86_64
 #include "unref.h"
-slimPreEvaluatorDistanceSHP::slimPreEvaluatorDistanceSHP(std::string shp, double dmax, std::vector<int> tags, std::string ellipsoid) {
+slimPreEvaluatorDistance::slimPreEvaluatorDistance(std::string shp, double dmax, std::vector<int> tags, std::string ellipsoid) {
   _ellipsoid = ellipsoid;
   auto vertices = unref::readSHP(shp.c_str(), dmax, ellipsoid);
   int n = 0;
@@ -15,7 +18,6 @@ slimPreEvaluatorDistanceSHP::slimPreEvaluatorDistanceSHP(std::string shp, double
   }
   _points = new flann::Matrix<float>(new float[3*n],n,3);
   n = 0;
-
   for (const auto &v : vertices) {
     if ((!tags.empty()) && std::find(tags.begin(), tags.end(), v->tag()) == tags.end()) continue;
     (*_points)[n][0] = v->x();
@@ -27,7 +29,56 @@ slimPreEvaluatorDistanceSHP::slimPreEvaluatorDistanceSHP(std::string shp, double
   _index->buildIndex();
 }
 
-slimPreEvaluatorDistanceSHP::~slimPreEvaluatorDistanceSHP() {
+slimPreEvaluatorDistance::slimPreEvaluatorDistance(std::vector<dgGroupOfFaces*> groups, std::string meshProj, double dmax) {
+  fullMatrix<double> coord;
+  std::vector<std::vector<double> > tmpcoord;
+
+  projPJ pjLocal = pj_init_plus(meshProj.c_str());
+  projPJ pjGlobal = pj_init_plus("+proj=geocent +ellps=sphere +a=1 +b=1");
+  int n = 0;
+  for (size_t iFaceGroup = 0; iFaceGroup < groups.size(); ++iFaceGroup){
+    dgGroupOfFaces *faceGroup= groups[iFaceGroup];
+    if (faceGroup->getNbNodes() != 2){
+      Msg::Warning("Only first order 1D interfaces are currently taken into account for evaluate distance.");
+      continue;
+    }
+    const dgFullMatrix<double> &groupCoord = faceGroup->elementGroup(0).getGroupCollection()->coordinatesDof().getGroupProxy(&faceGroup->elementGroup(0));
+    for (size_t iFace = 0; iFace < faceGroup->size(); ++iFace){
+      const std::vector<int> &cl = faceGroup->closure(iFace, 0);
+      size_t iElem = faceGroup->elementId(iFace, 0);
+      groupCoord.getBlockProxy(iElem, coord);
+      double x0 = coord(cl[0], 0);
+      double y0 = coord(cl[0], 1);
+      double z0 = coord(cl[0], 2);
+      double x1 = coord(cl[1], 0);
+      double y1 = coord(cl[1], 1);
+      double z1 = coord(cl[1], 2);
+      pj_transform(pjLocal,pjGlobal, 1, 1, &x0, &y0, &z0);
+      pj_transform(pjLocal,pjGlobal, 1, 1, &x1, &y1, &z1);
+      double dist = sqrt((x0-x1)*(x0-x1)+(y0-y1)*(y0-y1)+(z0-z1)*(z0-z1));
+      int nPoints = ceil(dist/dmax) + 2; 
+      for (int i = 0; i < nPoints; ++i){
+        double px = x0 * (1-i/(nPoints-1)) + x1 * (i/(nPoints-1));
+        double py = y0 * (1-i/(nPoints-1)) + y1 * (i/(nPoints-1));
+        double pz = z0 * (1-i/(nPoints-1)) + z1 * (i/(nPoints-1));
+        std::vector<double> tmp = {px, py, pz};
+        tmpcoord.push_back(tmp);
+        n +=1;
+      }
+    }
+  }
+  pj_free(pjLocal);
+  pj_free(pjGlobal);
+
+  _points = new flann::Matrix<float>(new float[3*n],n,3);
+  for (int i = 0; i < n; ++i)
+    for (int j = 0; j < 3; ++j)
+      (*_points)[i][j] = tmpcoord[i][j];
+  _index = new flann::Index<flann::L2<float> >(*_points, flann::KDTreeSingleIndexParams());
+  _index->buildIndex();
+}
+
+slimPreEvaluatorDistance::~slimPreEvaluatorDistance() {
   if (_index) {
     delete _index;
   }
@@ -37,7 +88,7 @@ slimPreEvaluatorDistanceSHP::~slimPreEvaluatorDistanceSHP() {
   }
 }
 
-double slimPreEvaluatorDistanceSHP::evaluate(double x, double y, double z) {
+double slimPreEvaluatorDistance::evaluate(double x, double y, double z) {
   float query[3] = {(float)x,(float)y,(float)z};
   int id;
   float dist;

modules/slimFunction/slimPreEvaluatorDistanceSHP.h → modules/slimFunction/slimPreEvaluatorDistance.h

-#ifndef SLIM_PRE_EVLUATOR_DISTANCE_SHP_H_
-#define SLIM_PRE_EVLUATOR_DISTANCE_SHP_H_
+#ifndef SLIM_PRE_EVLUATOR_DISTANCE_H_
+#define SLIM_PRE_EVLUATOR_DISTANCE_H_
 #include "slimFunctionConfig.h"
+class dgGroupOfFaces;
 #ifdef HAVE_PROJ
 #ifdef HAVE_UNREF
 #ifdef HAVE_FLANN
@@ -9,16 +10,17 @@
 #include <cmath> // required for flann
 #include "flann/flann.h"
 #endif
-class slimPreEvaluatorDistanceSHP:public slimPreEvaluator{
+class slimPreEvaluatorDistance:public slimPreEvaluator{
 #ifndef SWIG
   flann::Index<flann::L2<float> > *_index;
   flann::Matrix<float> *_points;
   std::string _ellipsoid;
 #endif
  public:
-  slimPreEvaluatorDistanceSHP(std::string shp, double dmax, std::vector<int> tags = std::vector<int>(), std::string ellipsoid="+ellps=WGS84");
   std::string coordinateSystem() {return "+proj=geocent " + _ellipsoid;}
-  ~slimPreEvaluatorDistanceSHP();
+  slimPreEvaluatorDistance(std::string shp, double dmax, std::vector<int> tags = std::vector<int>(), std::string ellipsoid="+ellps=WGS84");
+  slimPreEvaluatorDistance(std::vector<dgGroupOfFaces*> groups, std::string meshProj, double dmax);
+  ~slimPreEvaluatorDistance();
   double evaluate(double x, double y, double z) ;
 };
 #endif