Native Mesh Contouring for Vertex Fields

[andrewdnolan]

Vertex contours on the native mesh can produce saddle points, which violate the cycle/path assumptions of cell contours.

More work will be required to handle these saddle points.

[andrewdnolan]

import io
from contextlib import redirect_stdout

import matplotlib.pyplot as plt
import numpy as np

from mpas_tools.mesh.conversion import convert, cull
from mpas_tools.planar_hex import make_planar_hex_mesh
from mpas_tools.translate import center

import mosaic
from mosaic.contour import MPASContourGenerator


def sinusoid(x, y, A=2.0, N=1.5):
    """ """
    x_min = np.min(x)
    y_min = np.min(y)
    x_max = np.max(x)
    y_max = np.max(y)

    L_x = (x_max - x_min) / 2
    L_y = (y_max - y_min) / 2
    return A * np.sin(N * np.pi * (x / L_x)) * np.sin(N * np.pi * (y / L_y))


def planar_hex_mesh(N=100):
    """Wrapper for MPAS-Tools functions"""
    output_capture = io.StringIO()

    # Redirect stdout to the buffer within the 'with' block
    with redirect_stdout(output_capture):
        ds = make_planar_hex_mesh(
            nx=N, ny=N, dc=1 / N, nonperiodic_x=True, nonperiodic_y=True
        )
        ds = cull(ds)
        ds = convert(ds)
        # returns None, does centering inplace
        center(ds)

    return ds

if __name__ == "__main__":

    ds = planar_hex_mesh(50)
    
    field = sinusoid(ds.xVertex, ds.yVertex)
    descriptor = mosaic.Descriptor(ds)
    
    mask = field < 0.0
    
    generator = MPASContourGenerator(descriptor, field)
    
    graph = generator._create_vertex_contour_graph(mask, filled=False)
    
    fig, ax = plt.subplots(layout='constrained')
    
    mosaic.polypcolor(ax, descriptor, field, alpha=1 / 3, ec="k", lw=1 / 5)
    
    
    for i, j in graph.edges:
        ax.plot(
            descriptor.ds.xCell[[i,j]],
            descriptor.ds.yCell[[i,j]],
            c='k'
        )
    
    ax.scatter(
        generator.ds.xEdge[generator.boundary_edges],
        generator.ds.yEdge[generator.boundary_edges],
        c='r',
        s=2.5,
        alpha=0.5
    )

    plt.show()