Example 2 - stripy predefined meshes#
One common use of stripy is in meshing x,y coordinates and, to this end, we provide pre-defined meshes for square and elliptical triangulations. A random mesh is included as a counterpoint to the regular meshes. Each of these meshes is also a Triangulation object.
The Cartesian mesh classes in stripy are:
stripy.cartesian_meshes.square_mesh(extent, spacingX, spacingY)
stripy.cartesian_meshes.elliptical_mesh(extent, spacingX, spacingY)
stripy.cartesian_meshes.random_mesh(extent, number_of_points=5000)
Any of the above meshes can be uniformly refined by specifying the refinement_levels parameter. The square and elliptical meshes come with a random_scale parameter that specifies the amount of random noise to be added to the mesh (random_scale=0 by default).
Notebook contents#
The next example is Ex3-Interpolation
Sample meshes#
We create a number of meshes from the basic types available in stripy with approximately similar numbers of vertices.
import stripy
## A bunch of meshes with roughly similar overall numbers of points / triangles
xmin = 0.0
xmax = 10.0
ymin = 0.0
ymax = 10.0
extent = [xmin, xmax, ymin, ymax]
spacingX = 1.0
spacingY = 1.0
nsamples = 5000
str_fmt = "{:25} {:3}\t{:6}"
square0 = stripy.cartesian_meshes.square_mesh(extent, spacingX, spacingY, refinement_levels=0)
square2 = stripy.cartesian_meshes.square_mesh(extent, spacingX, spacingY, refinement_levels=2)
squareR = stripy.cartesian_meshes.square_mesh(extent, spacingX, spacingY, refinement_levels=4)
print(str_fmt.format('Square mesh', square0.npoints, squareR.npoints))
ellip0 = stripy.cartesian_meshes.elliptical_mesh(extent, spacingX, spacingY, refinement_levels=0)
ellip2 = stripy.cartesian_meshes.elliptical_mesh(extent, spacingX, spacingY, refinement_levels=2)
ellipR = stripy.cartesian_meshes.elliptical_mesh(extent, spacingX, spacingY, refinement_levels=4)
print(str_fmt.format('Elliptical mesh', ellip0.npoints, ellipR.npoints))
randR = stripy.cartesian_meshes.random_mesh(extent, nsamples)
rand0 = stripy.Triangulation(randR.x[::50], randR.y[::50])
rand2 = stripy.Triangulation(randR.x[::25], randR.y[::25])
print(str_fmt.format('Random mesh', rand0.npoints, randR.npoints))
Square mesh 121 25921
Elliptical mesh 21 1301
(276,) (276,)
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[1], line 31
28 print(str_fmt.format('Elliptical mesh', ellip0.npoints, ellipR.npoints))
30 randR = stripy.cartesian_meshes.random_mesh(extent, nsamples)
---> 31 rand0 = stripy.Triangulation(randR.x[::50], randR.y[::50])
32 rand2 = stripy.Triangulation(randR.x[::25], randR.y[::25])
35 print(str_fmt.format('Random mesh', rand0.npoints, randR.npoints))
File ~/.local/lib/python3.12/site-packages/stripy/cartesian.py:141, in Triangulation.__init__(self, x, y, refinement_levels, permute, tree)
138 self.permute = permute
139 self.tree = tree
--> 141 self._update_triangulation(x, y)
143 for r in range(0,refinement_levels):
144 x, y = self.uniformly_refine_triangulation(faces=False, trisect=False)
File ~/.local/lib/python3.12/site-packages/stripy/cartesian.py:196, in Triangulation._update_triangulation(self, x, y)
194 raise ValueError('ierr={} in trmesh\n{}'.format(ierr, _ier_codes[9999]))
195 elif ierr != 0:
--> 196 raise ValueError('ierr={} in trmesh\n{}'.format(ierr, _ier_codes[ierr]))
198 self.npoints = npoints
199 self._x = x
ValueError: ierr=-2 in trmesh
the first three nodes are collinear.
Set permute to True or reorder nodes manually.
print("Square: {}".format(square0.__doc__))
print("Elliptical: {}".format(ellip0.__doc__))
print("Random: {}".format(randR.__doc__))
Square:
A square mesh where points are uniformly populated
along x and y directions defined by extent.
Elliptical:
An elliptical mesh where points are successively populated at an
increasing radius from the midpoint of the extent.
Caution in parallel and for reproducibility - random noise in point locations !
Random:
A mesh of random points. Take care if you use this is parallel
as the location of points will not be the same on all processors.
Analysis of the characteristics of the triangulations#
We plot a histogram of the (spherical) areas of the triangles in each of the triangulations normalised by the average area. This is one measure of the uniformity of each mesh.
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
def area_histo(mesh):
freq, area_bin = np.histogram(mesh.areas(), bins=20, normed=True)
area = 0.5 * (area_bin[1:] + area_bin[:-1])
norm_area = area / mesh.areas().mean()
return norm_area, 0.25 * freq*area / np.pi**2
def add_plot(axis, mesh, xlim, ylim):
u, v = area_histo(mesh)
width = (u.max() - u.min()) / 30.
axis.bar(u, v, width=width)
axis.set_ylim(ylim)
axis.plot([1.0,1.0], [0.0,1.5], linewidth=1.0, linestyle="-.", color="Black")
return
fig, ax = plt.subplots(1,3, figsize=(14,4))
xlim=(0.75,1.5)
ylim=(0.0,0.05)
# square
add_plot(ax[0], squareR, xlim, ylim)
# elliptical
add_plot(ax[1], ellipR, xlim, ylim)
# random (this one is very different from the others ... )
add_plot(ax[2], randR, xlim, ylim)
fig.savefig("AreaDistributionsByMesh.png", dpi=250, transparent=True)
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
Cell In[3], line 1
----> 1 get_ipython().run_line_magic('matplotlib', 'inline')
2 import matplotlib.pyplot as plt
3 import numpy as np
File ~/.local/lib/python3.12/site-packages/IPython/core/interactiveshell.py:2482, in InteractiveShell.run_line_magic(self, magic_name, line, _stack_depth)
2480 kwargs['local_ns'] = self.get_local_scope(stack_depth)
2481 with self.builtin_trap:
-> 2482 result = fn(*args, **kwargs)
2484 # The code below prevents the output from being displayed
2485 # when using magics with decorator @output_can_be_silenced
2486 # when the last Python token in the expression is a ';'.
2487 if getattr(fn, magic.MAGIC_OUTPUT_CAN_BE_SILENCED, False):
File ~/.local/lib/python3.12/site-packages/IPython/core/magics/pylab.py:103, in PylabMagics.matplotlib(self, line)
98 print(
99 "Available matplotlib backends: %s"
100 % _list_matplotlib_backends_and_gui_loops()
101 )
102 else:
--> 103 gui, backend = self.shell.enable_matplotlib(args.gui)
104 self._show_matplotlib_backend(args.gui, backend)
File ~/.local/lib/python3.12/site-packages/IPython/core/interactiveshell.py:3667, in InteractiveShell.enable_matplotlib(self, gui)
3664 import matplotlib_inline.backend_inline
3666 from IPython.core import pylabtools as pt
-> 3667 gui, backend = pt.find_gui_and_backend(gui, self.pylab_gui_select)
3669 if gui != None:
3670 # If we have our first gui selection, store it
3671 if self.pylab_gui_select is None:
File ~/.local/lib/python3.12/site-packages/IPython/core/pylabtools.py:338, in find_gui_and_backend(gui, gui_select)
321 def find_gui_and_backend(gui=None, gui_select=None):
322 """Given a gui string return the gui and mpl backend.
323
324 Parameters
(...)
335 'WXAgg','Qt4Agg','module://matplotlib_inline.backend_inline','agg').
336 """
--> 338 import matplotlib
340 if _matplotlib_manages_backends():
341 backend_registry = matplotlib.backends.registry.backend_registry
ModuleNotFoundError: No module named 'matplotlib'
Plot and compare the predefined meshes#
def mesh_fig(mesh, meshR, name):
fig = plt.figure(figsize=(10, 10), facecolor="none")
ax = plt.subplot(111)
ax.axis('off')
generator = mesh
refined = meshR
x0 = generator.x
y0 = generator.y
xR = refined.x
yR = refined.y
ax.scatter(x0, y0, color="Red", marker="o", s=150.0)
ax.scatter(xR, yR, color="DarkBlue", marker="o", s=50.0)
ax.triplot(xR, yR, refined.simplices, color="black", linewidth=0.5)
fig.savefig(name, dpi=250, transparent=True)
return
mesh_fig(square0, square2, "Square" )
mesh_fig(ellip0, ellip2, "Elliptical" )
mesh_fig(rand0, rand2, "Random" )
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[4], line 26
22 fig.savefig(name, dpi=250, transparent=True)
24 return
---> 26 mesh_fig(square0, square2, "Square" )
28 mesh_fig(ellip0, ellip2, "Elliptical" )
30 mesh_fig(rand0, rand2, "Random" )
Cell In[4], line 3, in mesh_fig(mesh, meshR, name)
1 def mesh_fig(mesh, meshR, name):
----> 3 fig = plt.figure(figsize=(10, 10), facecolor="none")
4 ax = plt.subplot(111)
5 ax.axis('off')
NameError: name 'plt' is not defined
The next example is Ex3-Interpolation