Scalar Paving

This example illustrates how to influence the size of elements in the generated 2d mesh by specifying a size function. This process is referred to as scalar paving. The size function is specified using the InterpLinear or InterpIdw classes. The InterpLinear class performs spatial interpolation from points and triangles. This example uses a simple polygon with a set of 5 points and 4 triangles to define a linear size function.

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 param
import hvplot.pandas
import hvplot.xarray

import meshing_tools
import xmsmesh
import xmsinterp

hv.extension('bokeh')

Some Helper Functions

In [2]:
def pts_to_dataframe(pts):
    np_pts = np.array(pts)
    return pd.DataFrame(np_pts, columns=['x', 'y', 'z'])

def tri_list_to_dataframe(tris):
    np_tris = np.array(tris)
    return pd.DataFrame([(np_tris[x], np_tris[x+1], np_tris[x+2]) for x in range(0, len(tris), 3)], columns=['v0', 'v1', 'v2'])

Define the Polygon Points and Display

In [3]:
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[3]:

Setup Linear Interpolation Object for Size Function

Create linear interpolator for the size function.

In [4]:
df_pts = pts_to_dataframe(((-10, -10, 10), (-10, 110, 10), (110, 110, 10), (110, -10, 10), (60, 70, 1)))
df_tris = tri_list_to_dataframe((0, 4, 1, 1, 4, 2, 2, 4, 3, 3, 4, 0))
vert_points = hv.Points(df_pts)
trimesh = hv.TriMesh((df_tris, vert_points))
polygons * points * trimesh
Out[4]:
In [5]:
linear = xmsinterp.interpolate.InterpLinear(pts=df_pts.values)
_ = linear.interp_to_pt((0, 0, 0)) # TODO: Fix constructor... This will build the triangles

Create the Mesh Inputs and Generate the Mesh

In [6]:
input_polygon = xmsmesh.meshing.PolyInput(out_poly=polygon_points, size_function=linear)
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 [7]:
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[7]:

Setup IDW Interpolation Object for Size Function and Generate Mesh

Create idw interpolator for the size function.

In [8]:
# use idw as the size function
idw = xmsinterp.interpolate.InterpIdw(pts=df_pts.values)
#mesh_io.poly_inputs[0].size_function = idw
input_polygon = xmsmesh.meshing.PolyInput(out_poly=polygon_points, size_function=idw)
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 [9]:
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[9]: