changed names for better clarity

used broadcasting for cleaner code

docs/source/reference/tools.rst

@@ -68,7 +68,4 @@ Mesh Converter
 The following functions can be used to convert a mesh file supported by
 'meshio<https://github.com/nschloe/meshio>'_ to a set of surface points.
 
-.. autofunction:: pysph.tools.gen_geom_surf.gen_surf_points
-
-.. autofunction:: pysph.tools.gen_geom_surf.gen_surf_points_uniform
-
+.. autofunction:: pysph.tools.read_mesh.mesh2points

pysph/tools/geom_surf_points.pyx → pysph/tools/mesh_tools.pyx

@@ -276,9 +276,10 @@ cdef double dot(double[:] normal, double[:] point, double[:] face_centre):
            normal[1]*(point[1]-face_centre[1]) + \
            normal[2]*(point[2]-face_centre[2])
 
-def get_surface_points_uniform(x, y, z, cells, normals, dx_sph=1,
-                               h_sph=1, radius_scale=1.0,
-                               dx_triangle=None):
+
+def surf_points_uniform(x, y, z, cells, normals, dx_sph=1,
+                        h_sph=1, radius_scale=1.0,
+                        dx_triangle=None):
     """Generates points to cover surface described by a set of points
     The function generates a grid with a spacing of dx_sph and keeps points
     on the grid which lie within the object.
@@ -344,7 +345,7 @@ def get_surface_points_uniform(x, y, z, cells, normals, dx_sph=1,
     return xf, yf, zf
 
 
-def get_surface_points(x, y, z, cells, dx_triangle):
+def surface_points(x, y, z, cells, dx_triangle):
     """ Generates points to cover surface described by given set of connected points
     Parameters
     ----------

pysph/tools/read_mesh.py

@@ -0,0 +1,92 @@
+'''
+Functions can be used to generate points to describe a given input
+mesh.
+
+Supported mesh formats: All file formats supported by meshio.
+(https://github.com/nschloe/meshio)
+'''
+
+import numpy as np
+import meshio
+from pysph.tools.mesh_tools import surface_points, surf_points_uniform
+
+
+class Mesh:
+    def __init__(self, file_name, file_type=None):
+        if file_type is None:
+            self.mesh = meshio.read(file_name)
+        else:
+            self.mesh = meshio.read(file_name, file_type)
+
+        self.cells = np.array([], dtype=int).reshape(0, 3)
+
+    def extract_connectivity_info(self):
+        cell_blocks = self.mesh.cells
+        for block in cell_blocks:
+            self.cells = np.concatenate((self.cells, block.data))
+
+        return self.cells
+
+    def extract_coordinates(self):
+        x, y, z = self.mesh.points.T
+        self.x, self.y, self.z = x, y, z
+
+        return x, y, z
+
+    def compute_normals(self):
+        n = self.cells.shape[0]
+        self.normals = np.zeros((n, 3))
+        points = self.mesh.points
+
+        idx = self.cells
+        pts = points[idx]
+        a = pts[:, 1] - pts[:, 0]
+        b = pts[:, 2] - pts[:, 0]
+
+        normals = np.cross(a, b)
+        mag = np.linalg.norm(normals, axis=1)
+        mag.shape = (n, 1)
+        self.normals = normals/mag
+
+        return self.normals
+
+
+def mesh2points(file_name, dx, file_format=None, uniform=False):
+    '''
+    Generates points with a spacing dx to describe the surface of the
+    input mesh file.
+
+    Supported file formats: Refer to https://github.com/nschloe/meshio
+
+    Only works with triangle meshes.
+    Parameters
+    ----------
+    file_name : string
+        Mesh file name
+    dx : float
+        Required spacing between generated particles
+    file_format : str
+        Mesh file format
+    uniform : bool
+        If True generates points on a grid of spacing dx
+    Returns
+    -------
+    xf, yf, zf : ndarray
+        1d numpy arrays with x, y, z coordinates of covered surface
+    '''
+    mesh = Mesh(file_name)
+    cells = mesh.extract_connectivity_info()
+    x, y, z = mesh.extract_coordinates()
+
+    if uniform is False:
+        xf, yf, zf = surface_points(x, y, z, cells, dx)
+
+    else:
+        if file_format is 'stl':
+            normals = mesh.mesh.cell_data['facet_normals'][0]
+        else:
+            normals = mesh.compute_normals()
+
+        xf, yf, zf = surf_points_uniform(x, y, z, cells, normals, dx, dx)
+
+    return xf, yf, zf

pysph/tools/tests/test_geom_surf_points.py → pysph/tools/tests/test_mesh_tools.py

@@ -2,55 +2,59 @@
 import unittest
 import pytest
 from pysph.base.particle_array import ParticleArray
-import pysph.tools.geom_surf_points as G
+import pysph.tools.mesh_tools as G
 from pysph.base.utils import get_particle_array
 
+
 # Data of a unit length cube
-points = np.array([[0., 0., 0.],
-                   [0., 1., 0.],
-                   [1., 1., 0.],
-                   [1., 0., 0.],
-                   [0., 0., 1.],
-                   [0., 1., 1.],
-                   [1., 0., 1.],
-                   [1., 1., 1.]])
-
-x_cube, y_cube, z_cube = points.T
-
-cells = np.array([[0, 1, 2],
-                  [0, 2, 3],
-                  [0, 4, 5],
-                  [0, 5, 1],
-                  [0, 3, 6],
-                  [0, 6, 4],
-                  [4, 6, 7],
-                  [4, 7, 5],
-                  [3, 2, 7],
-                  [3, 7, 6],
-                  [1, 5, 7],
-                  [1, 7, 2]])
-
-normals = np.array([[0.,  0., -1.],
-                    [0.,  0., -1.],
-                    [-1.,  0.,  0.],
-                    [-1.,  0.,  0.],
-                    [0., -1.,  0.],
-                    [0., -1.,  0.],
-                    [0.,  0.,  1.],
-                    [0.,  0.,  1.],
-                    [1.,  0.,  0.],
-                    [1.,  0.,  0.],
-                    [0.,  1.,  0.],
-                    [0.,  1.,  0.]])
-
-vectors = np.zeros((len(cells), 3, 3))
-for i, cell in enumerate(cells):
-    idx1, idx2, idx3 = cell
-
-    vector = np.array([[x_cube[idx1], y_cube[idx1], z_cube[idx1]],
-                       [x_cube[idx2], y_cube[idx2], z_cube[idx2]],
-                       [x_cube[idx3], y_cube[idx3], z_cube[idx3]]])
-    vectors[i] = vector
+def cube_data():
+    points = np.array([[0., 0., 0.],
+                       [0., 1., 0.],
+                       [1., 1., 0.],
+                       [1., 0., 0.],
+                       [0., 0., 1.],
+                       [0., 1., 1.],
+                       [1., 0., 1.],
+                       [1., 1., 1.]])
+
+    x_cube, y_cube, z_cube = points.T
+
+    cells = np.array([[0, 1, 2],
+                      [0, 2, 3],
+                      [0, 4, 5],
+                      [0, 5, 1],
+                      [0, 3, 6],
+                      [0, 6, 4],
+                      [4, 6, 7],
+                      [4, 7, 5],
+                      [3, 2, 7],
+                      [3, 7, 6],
+                      [1, 5, 7],
+                      [1, 7, 2]])
+
+    normals = np.array([[0.,  0., -1.],
+                        [0.,  0., -1.],
+                        [-1.,  0.,  0.],
+                        [-1.,  0.,  0.],
+                        [0., -1.,  0.],
+                        [0., -1.,  0.],
+                        [0.,  0.,  1.],
+                        [0.,  0.,  1.],
+                        [1.,  0.,  0.],
+                        [1.,  0.,  0.],
+                        [0.,  1.,  0.],
+                        [0.,  1.,  0.]])
+
+    vectors = np.zeros((len(cells), 3, 3))
+    for i, cell in enumerate(cells):
+        idx1, idx2, idx3 = cell
+
+        vector = np.array([[x_cube[idx1], y_cube[idx1], z_cube[idx1]],
+                           [x_cube[idx2], y_cube[idx2], z_cube[idx2]],
+                           [x_cube[idx3], y_cube[idx3], z_cube[idx3]]])
+        vectors[i] = vector
+
+    return x_cube, y_cube, z_cube, cells, normals, vectors
 
 
 class TestGeometry(unittest.TestCase):
@@ -92,6 +96,7 @@ def test_fill_triangle_throws_polygon_mesh_error(self):
     def test_get_points_from_mgrid(self):
         """Find neighbouring particles around a unit cube"""
         h = 0.1
+        x_cube, y_cube, z_cube, cells, normals, vectors = cube_data()
 
         x, y, z, x_list, y_list, z_list, vectors = \
         G._get_surface_mesh(x_cube, y_cube, z_cube, cells, h, uniform=True)
@@ -123,29 +128,32 @@ def on_surface(x, y, z): return surface1(x, y, z) or \
     def test_get_surface_mesh(self):
         """Check if mesh is generated correctly for unit cube"""
 
+        x_cube, y_cube, z_cube, cells, normals, vectors = cube_data()
         x, y, z = G._get_surface_mesh(x_cube, y_cube, z_cube, cells, 0.1)
         h = np.finfo(float).eps
         self._cube_assert(x, y, z, h)
 
     def test_get_surface_points(self):
         """Check if surface is generated correctly for unit cube"""
         h = 0.1
-        x, y, z = G.get_surface_points(x_cube, y_cube, z_cube, cells, h)
+        x_cube, y_cube, z_cube, cells, normals, vectors = cube_data()
+        x, y, z = G.surface_points(x_cube, y_cube, z_cube, cells, h)
         self._cube_assert(x, y, z, h)
 
     def test_get_surface_points_uniform(self):
         """Check if uniform surface is generated correctly for unit cube"""
         h = 0.1
-        x, y, z = G.get_surface_points_uniform(x_cube, y_cube, z_cube,
-                                               cells, normals, 1.0, 1.0)
+        x_cube, y_cube, z_cube, cells, normals, vectors = cube_data()
+        x, y, z = G.surf_points_uniform(x_cube, y_cube, z_cube,
+                                        cells, normals, 1.0, 1.0)
         self._cube_assert(x, y, z, h)
 
     def test_prism(self):
         tri_normal = np.array([0, -1, 0])
         tri_points = np.array([[0, 0, 0], [1, 0, 0], [0, 0, 1]])
         h = 1/1.5
         prism_normals, prism_points, prism_face_centres = \
-          G.prism(tri_normal, tri_points, h)
+            G.prism(tri_normal, tri_points, h)
         assert np.array([-1, 0, 0]) in prism_normals
         assert np.array([0, 1, 0]) in prism_points
         assert np.array([0.5, 0.5, 0]) in prism_face_centres