# How to create random, unstructured mesh of surface

I am working on a supersonic boundary element method (https://github.com/usuaero/MachLine). In order to test its sensitivity to the boundary discretization, I'd like to be able to create "random" meshes of simple shapes (e.g. a cone) like what's shown here for a sphere, except using triangles.

Generating a set of random vertices shouldn't be difficult. What algorithm/tool could I use for then piecing those vertices together into a set of panels?

• Random points on the surface + Delaunay triangulation; but after that, if you don't want triangular elements, you'd have to combine adjacent triangles into other polygons. Commented Jul 26, 2022 at 21:13
• What about using gmsh or gid? gmsh is free, and gid is charged. Commented Jul 26, 2022 at 22:19
• Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.
– Community Bot
Commented Jul 26, 2022 at 22:37
• The fluid-dynamics part of geophysics has a lot of simulations which are on planet-scale and consequently they have existing literature about which meshes of sphere surfaces are well suited. Commented Jul 27, 2022 at 12:16

One option could be to compute the 'restricted' Delaunay triangulation of the set of scattered points you place on the surface of the sphere/cone/object you're working with. This will determine the subset of triangular faces of the 3-dimensional Delaunay tessellation (tetrahedrons) induced by the points. The CGAL library, INRIA's geogram, or my jigsaw package are all options for computing rDTs.

Another option could be to use a surface meshing library to generate a high quality mesh for your surfaces, and then perturb the vertices to degrade its quality. I've had some success using this approach to test the robustness of order-of-accuracy estimates to mesh distortion.

Turns out scipy has an algorithm for generating convex hulls for a set of points, described here https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.ConvexHull.html. I'm not sure what the underlying algorithm is, but it did what I needed it to. Here's my solution coded up:

def _get_random_points_on_surface_of_sphere(N, r):
"""Generates a set of points randomly distributed across the surface of a sphere.

Parameters
----------
N : integer
Number of points.

r : float

Returns
-------
ndarray
Array of points.
"""

# Initialize
points = np.zeros((N,3))

# Get azimuth angles
psi = np.random.random(N)*2.0*np.pi

# Get elevation angles (weighted towards equator, so the distribution remains even)
theta = np.arccos(1.0-2.0*np.random.random(N))

# Get coordinates
points[:,0] = r*np.sin(theta)*np.cos(psi)
points[:,1] = r*np.sin(theta)*np.sin(psi)
points[:,2] = r*np.cos(theta)

return points

def generate_random_sphere(filename, N, r):
"""Generates a random unstructured mesh of a sphere.

Parameters
----------
filename : str
Name of the file to write the mesh to. Must have '.vtk' extension.
N : integer
Number of vertices.

r : float
Array of points.
"""

# Get vertices
verts = _get_random_points_on_surface_of_sphere(N, r)

# Create convex hull
hull = ConvexHull(verts)

# Make sure normal vector points inward
fixed_simplices = []
for simplex in hull.simplices:

# Calculate normal and centroid
centroid = np.sum(verts[simplex,:], axis=0) / 3.0
normal = np.cross(verts[simplex[1],:] - verts[simplex[0],:], verts[simplex[2],:] - verts[simplex[1],:])

# Check the normal points outward
if np.dot(centroid, normal) > 0.0:
fixed_simplices.append(simplex)
else:
fixed_simplices.append(simplex[::-1])

# Export mesh
_export_vtk(filename, verts, np.array(fixed_simplices))


Here's what one looks like with 500 vertices.