Example 7 - Refining a triangulation

Contents

Example 7 - Refining a triangulation¶

We have seen how the standard meshes can be uniformly refined to finer resolution. The routines used for this task are available to the stripy user for non-uniform refinement as well.

Notebook contents¶

Uniform meshes
Refinement strategies
Visualisation
Targetted refinement
Visualisation

import stripy as stripy
import numpy as np

Visualisation of targetted refinement¶

import k3d


## The four different triangulation strategies

t = [ edge_triangulations[-1],
      edge_t_triangulations[-1],
      centroid_triangulations[-1],
      centroid_t_triangulations[-1] ]

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


mesh_viewer = []
wire_viewer = []


for i in range(0,4):
    
    
    indices = t[i].simplices.astype(np.uint32)
    points  = np.column_stack(t[i].points.T).astype(np.float32)


    mesh_viewer.append(k3d.mesh(points, indices, wireframe=False, 
                   color=0x99AABB, 
                   name="mesh viewer {}".format(i+1),
                   flat_shading=True, opacity=1.0  ))

    wire_viewer.append(k3d.mesh(points, indices, wireframe=True, 
                   color=0x002244, 
                   name="wire frame viewer {}".format(i+1),
                   flat_shading=True, opacity=1.0  ))

    plot += mesh_viewer[i]
    plot += wire_viewer[i]
    

## This helps to manage the wireframe / transparency


indices = ico3.simplices.astype(np.uint32)
points  = np.column_stack(ico3.points.T).astype(np.float32)


background = k3d.mesh(points*0.9, indices, wireframe=False, 
                   color=0xBBBBBB, opacity=1.0, flat_shading=False  )

plot += background
plot.display()

## ## ## 

from ipywidgets import interact, interactive
import ipywidgets as widgets

choices = {  "edge triangulation": 0,
             "edge t triangulation": 1, 
             "centroid triangulation": 2, 
             "centroid t triangulation": 3  }

@interact(choice=choices.keys())
def chooser(choice):

    for i in range(0,4):
        mesh_viewer[i].visible = False
        wire_viewer[i].visible = False

    mesh_viewer[choices[choice]].visible = True
    wire_viewer[choices[choice]].visible = True
    
    return 

import matplotlib.pyplot as plt

%matplotlib inline

import cartopy
import cartopy.crs as ccrs

def mesh_fig(mesh, meshR, name):

    fig = plt.figure(figsize=(10, 10), facecolor="none")
    ax  = plt.subplot(111, projection=ccrs.Orthographic(central_longitude=0.0, central_latitude=0.0, globe=None))
    ax.coastlines(color="lightgrey")
    ax.set_global()

    generator = mesh
    refined   = meshR

    lons0 = np.degrees(generator.lons)
    lats0 = np.degrees(generator.lats)

    lonsR = np.degrees(refined.lons)
    latsR = np.degrees(refined.lats)


    ax.scatter(lons0, lats0, color="Red",
                marker="o", s=150.0, transform=ccrs.PlateCarree())

    ax.scatter(lonsR, latsR, color="DarkBlue",
                marker="o", s=50.0, transform=ccrs.PlateCarree())

    ax.scatter(np.degrees(points[:,0]), np.degrees(points[:,1]), marker="s", s=50, 
               color="#885500", transform=ccrs.PlateCarree())

    segs = refined.identify_segments()

    for s1, s2 in segs:
        ax.plot( [lonsR[s1], lonsR[s2]],
                 [latsR[s1], latsR[s2]], 
                 linewidth=0.5, color="black", transform=ccrs.Geodetic())

    fig.savefig(name, dpi=250, transparent=True)
    
    return



mesh_fig(edge_triangulations[0],        edge_triangulations[-1],     "EdgeByVertex" )

T = edge_triangulations[-1]
E = np.array(T.edge_lengths()).T
A = np.array(T.areas()).T
equant = np.max(E, axis=1) / np.min(E, axis=1)
size_ratio = np.sqrt(np.max(A) / np.min(A))
print ("EBV", T.simplices.shape[0], equant.max(), equant.min(), size_ratio)

mesh_fig(edge_t_triangulations[0],      edge_t_triangulations[-1],     "EdgeByTriangle" )


T = edge_t_triangulations[-1]
E = np.array(T.edge_lengths()).T
A = np.array(T.areas()).T
equant = np.max(E, axis=1) / np.min(E, axis=1)
size_ratio = np.sqrt(np.max(A) / np.min(A))
print ("EBT", T.simplices.shape[0], equant.max(), equant.min(), size_ratio)


mesh_fig(centroid_triangulations[0],    centroid_triangulations[-1],   "CentroidByVertex" )

T = centroid_triangulations[-1]
E = np.array(T.edge_lengths()).T
A = np.array(T.areas()).T
equant = np.max(E, axis=1) / np.min(E, axis=1)
size_ratio = np.sqrt(np.max(A) / np.min(A))
print ("CBV", T.simplices.shape[0], equant.max(), equant.min(), size_ratio)



mesh_fig(centroid_t_triangulations[0],  centroid_t_triangulations[-1], "CentroidByTriangle" )

T = centroid_t_triangulations[-1]
E = np.array(T.edge_lengths()).T
A = np.array(T.areas()).T
equant = np.max(E, axis=1) / np.min(E, axis=1)
size_ratio = np.sqrt(np.max(A) / np.min(A))
print ("CBT", T.simplices.shape[0], equant.max(), equant.min(), size_ratio)

The next example is Ex8-Spline-Tension

previous

Example 6 - Scattered Data and ‘Heat Maps’

next

Example 8 - Spline Tension