xmsmesh "hello world"

This is an introduction to the xmsmesh library. This library generates unstructured triangular or mixed triangle, quad meshes. The meshes are generated from input polygons. The resulting mesh will honor the input polygon boundaries. The density of mesh elements is controlled by the spacing of points along the polygon boundaries.

In addition, the density of elements can be controlled by specifying a size function (this is demonstrated in another notebook {see the examples folder}).

This example shows how to generate a mesh from a simple, square polygon with coordinates (0,0), (100,100) and a boundary spacing of 10.0.

In [1]:
import numpy as np
import geoviews as gv
import holoviews as hv
import pandas as pd
import xarray as xr
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import panel as pn

from shapely.geometry import Polygon
from holoviews.streams import Params

import hvplot.pandas
import hvplot.xarray

import meshing_tools
import xmsmesh

hv.extension('bokeh')

Define the Polygon Points and Display

In [2]:
# Define the polygon points
polygon_points = np.array([
    (0, 10,0),   (0, 20,0),   (0, 30,0),   (0, 40,0),    (0, 50,0),   (0, 60,0),   (0, 70,0),   (0, 80,0),
    (0, 90,0),   (0, 100,0),  (10, 100,0), (20, 100,0),  (30, 100,0), (40, 100,0), (50, 100,0), (60, 100,0),
    (70, 100,0), (80, 100,0), (90, 100,0), (100, 100,0), (100, 90,0), (100, 80,0), (100, 70,0), (100, 60,0),
    (100, 50,0), (100, 40,0), (100, 30,0), (100, 20,0),  (100, 10,0), (100, 0,0),  (90, 0,0),   (80, 0,0),
    (70, 0,0),   (60, 0,0),   (50, 0,0),   (40, 0,0),    (30, 0,0),   (20, 0,0),   (10, 0,0),   (0, 0,0)
])

# Display the Polygon
# We are using holoviews instead of geoviews here because we don't need projections for this.
%opts Polygons Points [width=500, height=500]
polygon = Polygon(polygon_points)
polygons = hv.Polygons([polygon,]).redim.range(x=(-15, 115), y=(-15, 115))
points = hv.Points(polygon_points).options(color='black', size=6)
polygons * points
Out[2]:

Create the Mesh Inputs and Generate the Mesh

In [3]:
# Create Meshing Inputs
input_polygon = xmsmesh.meshing.PolyInput(out_poly=polygon_points)
mesh_io = xmsmesh.meshing.MultiPolyMesherIo(poly_inputs=[input_polygon])

# Generate Mesh
succeded, errors = xmsmesh.meshing.mesh_utils.generate_mesh(mesh_io=mesh_io)
if succeded:
    print("Meshing was successful")
else:
    print("Meshing errors found:")
    for err in errors:
        print("\t{}".format(err))

Meshing was successful

Display the Mesh

In [4]:
pts, cells = meshing_tools.xmsmesh_to_dataframe(mesh_io.points, mesh_io.cells)
vert_points = hv.Points(pts, vdims=['z'])
trimesh = hv.TriMesh((cells, vert_points))
polygons * points * trimesh.edgepaths.options(line_width=0.5)
Out[4]:

Read polygon from shapefile

In [5]:
from shapely.geometry import shape
import fiona
shape_file = fiona.open("square_holes_polys.shp")
first = next(iter(shape_file))
shp_geom = shape(first['geometry'])

x, y = shp_geom.exterior.coords.xy
coords = np.array(list(zip(x, y)))

holes = []
for hole in shp_geom.interiors:
    xh, yh = hole.coords.xy
    hole_coords = np.array(list(zip(xh, yh)))
    holes.append(hole_coords)

# Create Meshing Inputs
input_polygon = xmsmesh.meshing.PolyInput(out_poly=coords[:-1], inside_polys=[x[:-1] for x in holes])
mesh_io = xmsmesh.meshing.MultiPolyMesherIo(poly_inputs=[input_polygon])

# Generate Mesh
succeded, errors = xmsmesh.meshing.mesh_utils.generate_mesh(mesh_io=mesh_io)
if succeded:
    print("Meshing was successful")
else:
    print("Meshing errors found:")
    print("{}".format(errors))
    
polygon_w_hole = Polygon(coords, holes)
polygons_holes = hv.Polygons([polygon_w_hole,]).redim.range(x=(-15, 115), y=(-15, 115))

pts, cells = meshing_tools.xmsmesh_to_dataframe(mesh_io.points, mesh_io.cells)
vert_points = hv.Points(pts, vdims=['z'])
trimesh = hv.TriMesh((cells, vert_points))
polygons_holes * trimesh.edgepaths.options(line_width=0.5)

Meshing was successful
Out[5]: