1.1 - Visualisation¶

from quagmire.tools import meshtools
from quagmire import QuagMesh
import numpy as np

Structured grids¶

minX, maxX = -5.0, 5.0
minY, maxY = -5.0, 5.0

resX = 75
resY = 75

DM = meshtools.create_DMDA(minX, maxX, minY, maxY, resX, resY)

print(type(DM))

<class 'petsc4py.PETSc.DMDA'>

This is a native PETSc data management object for structured grids (DMDA). This object has a number of useful methods and attached data which can be listed with

help(DM)

We hand this to QuagMesh to generate the necessary data structures for gradient operations, smoothing, neighbour allocation, etc.

mesh = QuagMesh(DM)

Underlying Mesh type: PixMesh
- cKDTree 0.0020817920012632385s
- Find boundaries 0.0005966669996269047s
- Construct neighbour cloud arrays 0.04875845900096465s, (0.028369209001539275s + 0.02036183300515404s)
- Construct rbf weights 0.004284041999198962s

We attach data to a mesh solely through mesh variables (see Example notebook for details)

    mesh_variable = mesh.add_variable(name="data1")
    mesh_variable.data = np.sin(mesh.coords[:,0] * np.pi)
    mesh_variable.sync()

The sync operation ensures data is coherent across processors - it is harmless and relatively inexpensive so is safe to use even in cases like this where there is no way for information to be out of sync between domains.

mesh_variable = mesh.add_variable(name="data1")
mesh_variable.data = np.sin(mesh.coords[:,0] * np.pi)
mesh_variable.sync()

mesh_variable2 = mesh.add_variable(name="data2")
mesh_variable2.data = np.sin(mesh.coords[:,0] * 0.2 * np.pi) * np.cos(mesh.coords[:,1] * 0.2 * np.pi) 
mesh_variable2.sync()

XY = mesh.coords
XYZ = np.vstack((XY[:,0],XY[:,1], mesh_variable2.data)).T

# import ipyvolume as ipv
# from matplotlib import cm
# norm = cm.colors.Normalize(Z.min(), Z.max())
# color = cm.coolwarm(norm(Z))

# ipv.figure(width=800, height=500)
# #surf = ipv.plot_surface(X, Y, Z, color=color)
# #wire = ipv.plot_wireframe(X, Y, Z, color=color)
# balls = ipv.scatter(X.reshape(-1),Y.reshape(-1),mesh_variable2.data, marker="sphere")
# ipv.zlim(-2,2)
# ipv.show()

import numpy as np
import k3d

plot = k3d.plot(camera_auto_fit=True, grid_visible=False)

plot += k3d.points(XYZ, point_size=0.1)

plot += k3d.surface(mesh_variable2.data.reshape((mesh.nx,mesh.ny)), 
                   xmin=mesh.minX, xmax=mesh.maxX, ymin=mesh.minY, ymax=mesh.maxY,
                   color=0x995500
                   )
plot.display()

import lavavu

lv = lavavu.Viewer(border=False, background="#FFFFFF", resolution=[500,500], near=-10.0)

lvmesh = lv.quads(dims=(mesh.nx, mesh.ny), wireframe=True)
lvmesh.vertices(mesh.coords)
lvmesh.values( mesh_variable.data, "sinx")
lvmesh.colourmap("#FF0000, #555555 #0000FF", range=[-1.0,1.0])

# The mesh can be given a height mapping like this

vertices = np.zeros((mesh.coords.shape[0],3))
vertices[:,0:2] = mesh.coords
vertices[:,2]   = mesh_variable2.data * 0.5

lvmesh2 = lv.quads(dims=(mesh.nx, mesh.ny), wireframe=False)
lvmesh2.vertices(vertices)
lvmesh2.values(mesh_variable2.data,"sinxcosy")
lvmesh2.colourmap("#FF0000, #FFFFFF:0.5 #0000FF", range=[-1.0,1.0])


lv.control.Panel()
lv.control.ObjectList()
lv.control.show()

Unstructured meshes¶

This is handled by PETSc’s DMPlex object, which requires the connectivity of a set of points. The connectivity between points can be triangulated using the built-in mesh creation tools:

x, y, simplices = square_mesh(minX, maxX, minY, maxY, spacingX, spacingY)
x, y, simplices = elliptical_mesh(minX, maxX, minY, maxY, spacingX, spacingY)

and handed to DMPlex using:

DM = meshtools.create_DMPlex(x, y, simplices, boundary_vertices=None, refinement_levels=0)

Alternatively, an arbitrary set of points (without duplicates) can be triangulated and processed as a DMPlex object using:

meshtools.create_DMPlex_from_points(x, y, bmask=None, refinement_levels=0)

If no boundary information is provided, the boundary is assumed to be the convex hull of points.

Parallel notes¶

The triangulation from the root processor is distributed to other processors using the DM object, including boundary points and boundary edges. The mesh can be refined efficiently in parallel using the refine_dm method. The order of this operation is important:

Triangulate points
Mark boundary edges
Distribute DMPlex to other processors
Refine the mesh

Elliptical mesh¶

spacingX = 0.1
spacingY = 0.1

x, y, simplices = meshtools.elliptical_mesh(minX, maxX, minY, maxY, spacingX, spacingY)
DM = meshtools.create_DMPlex(x, y, simplices)

mesh = QuagMesh(DM)

mesh_equant = mesh.neighbour_cloud_distances.mean(axis=1) / ( np.sqrt(mesh.area))

# import ipyvolume as ipv

# a = np.arange(-5, 5)
# X = mesh.tri.x
# Y = mesh.tri.y
# Z = mesh.tri.x * 0.0

# from matplotlib import cm
# norm = cm.colors.Normalize(mesh_equant.min(), mesh_equant.max())
# color = cm.coolwarm(norm(mesh_equant))

# ipv.figure(width=800, height=500)
# ipv.plot_trisurf(X, Y, Z, triangles=mesh.tri.simplices, color=color )
# ipv.show()

import lavavu

lv = lavavu.Viewer(border=False, background="#FFFFFF", resolution=[1000,600], near=-10.0)

# lavavu also works in 3D - so need to stitch in a Z component  (zero or a height)

vertices = np.zeros((mesh.tri.points.shape[0],3))
vertices[:,0:2] = mesh.tri.points
# vertices[:,2] = heights

bnodes = lv.points("Boundary Points", pointsize=10.0, pointtype="shiny", colour="red", opacity=0.75)
bnodes.vertices(vertices[~mesh.bmask])

nodes = lv.points("All Points", pointsize=10.0, pointtype="shiny", colour="blue", opacity=0.75)
nodes.vertices(vertices)

simp = lv.triangles("Triangle Edges", wireframe=True, colour="#442222", opacity=0.75)
simp.vertices(vertices)
simp.indices(mesh.tri.simplices)

tris = lv.triangles("Triangle Areas",  wireframe=False, colour="#77ff88", opacity=1.0)
tris.vertices(vertices-(0.0,0.0,0.01))
tris.indices(mesh.tri.simplices)
tris.values(mesh_equant, label="pointwise_area")

tris.colourmap("#000000, #FFFFFF")
cb = tris.colourbar()

lv.control.Panel()
lv.control.ObjectList()
lv.control.show()

lv.show()

Mesh improvement¶

Applies Lloyd’s algorithm of iterated voronoi construction to improve the mesh point locations. This distributes the points to a more uniform spacing with more equant triangles. It can be very slow for anything but a small mesh. Refining the mesh a few times will produce a large, well-spaced mesh.

bmask = mesh.bmask.copy()

x1, y1 = meshtools.lloyd_mesh_improvement(x, y, bmask, iterations=3)
DM = meshtools.create_DMPlex_from_points(x1, y1, bmask)

mesh1 = QuagMesh(DM)

mesh1_equant = mesh1.neighbour_cloud_distances.mean(axis=1) / ( np.sqrt(mesh1.area))

import lavavu

lv = lavavu.Viewer(border=False, background="#FFFFFF", resolution=[1000,600], near=-10.0)

# lavavu assumes 3D - so need to stitch in a Z component  (zero or a height)

vertices = np.zeros((mesh1.tri.points.shape[0],3))
vertices[:,0:2] = mesh1.tri.points

bnodes = lv.points("Boundary Points", pointsize=10.0, pointtype="shiny", colour="red", opacity=0.75)
bnodes.vertices(vertices[~mesh1.bmask])

simp = lv.triangles("Triangle Edges", wireframe=True, colour="#442222", opacity=0.75)
simp.vertices(vertices)
simp.indices(mesh1.tri.simplices)

tris = lv.triangles("Triangle Areas",  wireframe=False, colour="#77ff88", opacity=1.0)
tris.vertices(vertices-(0.0,0.0,0.01))
tris.indices(mesh1.tri.simplices)
tris.values(mesh1_equant, label="pointwise_area")

tris.colourmap("#000000, #FFFFFF", range=[1.0,0.9*mesh1_equant.mean()])

cb = tris.colourbar()

lv.control.Panel()
lv.control.ObjectList()
lv.control.show()

lv.show()

# Comparison of point-wise area for original and improved mesh

import matplotlib.pyplot as plt
%matplotlib inline

fig, (ax1, ax2) = plt.subplots(1,2, figsize=(10,4))

ax1.hist(mesh_equant, density=True)
ax2.hist(mesh1_equant, density=True)

ax1.set_title('original mesh')
ax2.set_title('improved mesh')

plt.show()

Spherical meshes¶

This unstructed mesh uses PETSc’s DMPlex object, and uses stripy to triangulate on the unit sphere. Multiple meshes may be created, including:

DM = meshtools.create_spherical_DMPlex(lons, lats, simplices, boundary_vertices=None)
DM = meshtools.create_DMPlex_from_spherical_points(lons, lats, simplices, bmask=None, refinement_levels=0)

If no boundary information is provided, the boundary is calculated from any line segments that do not share a triangle with another.a

import stripy

lons, lats, bmask = meshtools.generate_elliptical_points(-40, 40, -80, 80, 0.1, 0.1, 1500, 200)

icomesh = stripy.spherical_meshes.icosahedral_mesh(refinement_levels=5, include_face_points=True)

lons = np.degrees(icomesh.lons)
lats = np.degrees(icomesh.lats)
bmask = None

DM = meshtools.create_DMPlex_from_spherical_points(lons, lats, bmask, refinement_levels=0)
mesh0 = QuagMesh(DM)

Origin =  0.0 0.0 Radius =  40.0 Aspect =  2.0
Underlying Mesh type: sTriMesh
- Delaunay triangulation 0.04886054199994305s
- Calculate node weights and area 0.007725415999971119s
- Find boundaries 0.0037982080000347196s
- cKDTree 0.00373012500006098s
- Construct neighbour cloud arrays 0.1836289579999857s, (0.14000704199997926s + 0.0435912500000768s)
- Construct rbf weights 0.011983624999970743s

# lv = lavavu.Viewer(border=False, background="#FFFFFF", resolution=[1000,600], near=-10.0)

# bnodes0 = plot_points(lv, mesh0.data[~mesh0.bmask], "boundary_points_r0", colour="red", pointsize=10)
# bnodes1 = plot_points(lv, mesh1.data[~mesh1.bmask], "boundary_points_r1", colour="blue", pointsize=10)
# bnodes2 = plot_points(lv, mesh2.data[~mesh2.bmask], "boundary_points_r2", colour="#336611", pointsize=10)

# tri0 = plot_triangles(lv, mesh0.data, mesh0.tri.simplices, "mesh_r0", wireframe=True, linewidth=1.5, colour="red")
# tri1 = plot_triangles(lv, mesh1.data, mesh1.tri.simplices, "mesh_r1", wireframe=True, linewidth=1.0, colour="blue")
# tri2 = plot_triangles(lv, mesh2.data, mesh2.tri.simplices, "mesh_r2", wireframe=True, linewidth=1.0, colour="#336611")

# lv.control.Panel()
# lv.control.ObjectList()
# lv.control.show()

# lv.show()

# import ipyvolume as ipv

# a = np.arange(-5, 5)
# X = mesh0.tri.x
# Y = mesh0.tri.y
# Z = mesh0.tri.z
# C = mesh0.tri.lats

# from matplotlib import cm
# colormap = cm.coolwarm
# znorm = Z - Z.min()
# znorm /= znorm.ptp()
# znorm.min(), znorm.max()
# color = colormap(znorm)

# ipv.figure()
# globe = ipv.plot_trisurf(X, Y, Z, triangles=mesh0.tri.simplices)
# ipv.show()

from PIL import Image
# image = PIL.Image.open("/Users/lmoresi/Desktop/color_etopo1_ice_low_400.gif")
image = PIL.Image.open("/Users/lmoresi/cloudstor/Datasets4Teaching/Venus/RadarBrightness.jpg")
# image = PIL.Image.open("/Users/lmoresi/Downloads/GRAY_LR_SR_OB_DR/GRAY_LR_SR_OB_DR.tif")

# import ipyvolume as ipv

a = np.arange(-5, 5)
X = mesh0.tri.x
Y = mesh0.tri.y
Z = mesh0.tri.z
C = mesh0.tri.lats

u = mesh0.tri.lons / (2.0*np.pi)
v = mesh0.tri.lats / (1.0*np.pi) + 0.5

from matplotlib import cm
colormap = cm.coolwarm
znorm = Z - Z.min()
znorm /= znorm.ptp()
znorm.min(), znorm.max()
color = colormap(znorm)
color[:,3] = v

# fig = ipv.figure(width=1000)
# ipv.pylab.style.axes_off()
# ipv.pylab.style.box_off()


# globe = ipv.plot_trisurf(X, Y, Z, u=u, v=v, texture=image,
#                  triangles=mesh0.tri.simplices)
# ipv.show()

import io
image = PIL.Image.open("/Users/lmoresi/Downloads/GRAY_LR_SR_OB_DR/GRAY_LR_SR_OB_DR.tif")
image2 = image.resize((2050,1025))
buf = io.BytesIO()
image2.save(buf, format="tiff")
bpil = buf.getvalue()

# import ipyvolume as ipv

# u = np.linspace(0.0,    2.0*np.pi, num=180, endpoint=True) 
# v = np.linspace(-np.pi/2, np.pi/2, num=90,  endpoint=True)
# u, v = np.meshgrid(u, v)

# X,Y,Z = stripy.spherical.lonlat2xyz(u.reshape(-1),v.reshape(-1))
# X = X.reshape(u.shape)
# Y = Y.reshape(u.shape)
# Z = Z.reshape(u.shape)

# u /= 2.0 * np.pi
# v = v / np.pi + 0.5

# # v = st_mesh.y.reshape(180,90)

# from matplotlib import cm
# colormap = cm.coolwarm
# znorm = Z - Z.min()
# znorm /= znorm.ptp()
# znorm.min(), znorm.max()

# fig = ipv.figure()
# fig.style['axes']['visible'] = False
# fig.style['box']['visible'] = False

# globe = ipv.plot_mesh(X, Y, Z,  u=u, v=v, texture=image, wireframe=False )
# ipv.show()

# b = open("/Users/lmoresi/cloudstor/Datasets4Teaching/Venus/RadarBrightness.jpg", "rb").read()
# b1 = open("/Users/lmoresi/Desktop/color_etopo1_ice_low_400.gif", "rb").read()
b = open("GreyEarth.tif", "rb").read()
b1 = open("/Users/lmoresi/Downloads/GRAY_LR_SR_OB_DR/GRAY_LR_SR_OB_DR.tif", "rb").read()
b2 = open("/Users/lmoresi/cloudstor/Datasets4Teaching/Venus/RadarBrightness.jpg", "rb").read()
len(b1)

131327272

U = 0.5 + mesh0.tri.lons_map_to_wrapped(mesh0.tri.lons) / (2.0*np.pi) 
V = 0.5 - mesh0.tri.lats / np.pi

import k3d
plot = k3d.plot(camera_auto_fit=False)

indices = mesh0.tri.simplices.astype(np.uint32)
points  = np.column_stack(mesh0.tri.points.T).astype(np.float32)
uv      = np.column_stack((U,V)).astype(np.float32)

plt_mesh = k3d.mesh(points, indices, wireframe=False, color=0xFFFFFF,
                    texture=bpil, texture_file_format="tif" , uvs=uv, 
                    opacity=0.33
                   )

plt_mesh2 = k3d.mesh(points*0.9, indices, wireframe=False, color=0xFFFFFF,
                    texture=b2, texture_file_format="." , uvs=uv, 
                    opacity=1.0
                   )


plot += plt_mesh
plot += plt_mesh2
plot.display()

Save mesh to file¶

A mesh can be saved and imported later. The QuagMesh object has the save_mesh_to_hdf5 method for this, as does meshtools.

Note: Requires PETSc 3.8 or higher

filename = "Ex1-refined_mesh.h5"

# save from QuagMesh object:
# mesh2.save_mesh_to_hdf5(filename)

# save from meshtools:
meshtools.save_DM_to_hdf5(DM_r2, filename)

# load DM from file
DM_r2 = meshtools.create_DMPlex_from_hdf5(filename)

mesh2 = QuagMesh(DM_r2)

print(mesh2.npoints)
print(mesh2.area)

The next example is Ex2-Topography-Meshes

Quagmire.tex