Example 4 - stripy gradients

Example 4 - `stripy` gradients#

SRFPACK is a Fortran 77 software package that constructs a smooth interpolatory or approximating surface to data values associated with arbitrarily distributed points. It employs automatically selected tension factors to preserve shape properties of the data and avoid overshoot and undershoot associated with steep gradients.

The next example is Ex5-Smoothing

Define a computational mesh#

Use the (usual) icosahedron with face points included.

import stripy as stripy

xmin = 0.0
xmax = 10.0
ymin = 0.0
ymax = 10.0
extent = [xmin, xmax, ymin, ymax]

spacingX = 0.2
spacingY = 0.2

mesh = stripy.cartesian_meshes.elliptical_mesh(extent, spacingX, spacingY, refinement_levels=3)
mesh = stripy.Triangulation(mesh.x, mesh.y, permute=True)

print("number of points = {}".format(mesh.npoints))

number of points = 15523

Analytic function#

Define a relatively smooth function that we can interpolate from the coarse mesh to the fine mesh and analyse

import numpy as np

def analytic(xs, ys, k1, k2):
     return np.cos(k1*xs) * np.sin(k2*ys)

def analytic_ddx(xs, ys, k1, k2):
     return -k1 * np.sin(k1*xs) * np.sin(k2*ys) / np.cos(ys)

def analytic_ddy(xs, ys, k1, k2):
     return k2 * np.cos(k1*xs) * np.cos(k2*ys) 

analytic_sol = analytic(mesh.x, mesh.y, 0.1, 1.0)
analytic_sol_ddx = analytic_ddx(mesh.x, mesh.y, 0.1, 1.0)
analytic_sol_ddy = analytic_ddy(mesh.x, mesh.y, 0.1, 1.0)

%matplotlib inline
import matplotlib.pyplot as plt

def axis_mesh_field(fig, ax, mesh, field, label):

    ax.axis('off')

    x0 = mesh.x
    y0 = mesh.y
    
    trip = ax.tripcolor(x0, y0, mesh.simplices, field, cmap=plt.cm.RdBu)
    fig.colorbar(trip, ax=ax)
    
    ax.set_title(str(label))
    return

    
fig = plt.figure(figsize=(10, 8), facecolor="none")
ax = fig.add_subplot(111)
axis_mesh_field(fig, ax, mesh, analytic_sol, "analytic solution")

---------------------------------------------------------------------------
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
      4 def axis_mesh_field(fig, ax, mesh, field, label):

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'

Derivatives of solution compared to analytic values#

The gradient method of Triangulation takes a data array f representing values on the mesh vertices and returns the x,y derivatives.

Triangulation.gradient(f, nit=3, tol=0.001)

Derivatives of higher accuracy can be obtained by tweaking tol, which controls the convergence tolerance, or nit which controls the number of iterations to a solution. The default values are set to an optimal trade-off between speed and accuracy.

stripy_ddx, stripy_ddy = mesh.gradient(analytic_sol)
stripy_d2x, _ = mesh.gradient(stripy_ddx)
_, stripy_d2y = mesh.gradient(stripy_ddy)

fig, ax = plt.subplots(3,2, figsize=(12, 15), facecolor="none")

axis_mesh_field(fig, ax[0,0], mesh, analytic_sol, label="original")
axis_mesh_field(fig, ax[1,0], mesh, stripy_ddx, label="ddy")
axis_mesh_field(fig, ax[1,1], mesh, stripy_ddy, label="ddx")
axis_mesh_field(fig, ax[2,0], mesh, stripy_d2x, label="d2x")
axis_mesh_field(fig, ax[2,1], mesh, stripy_d2y, label="d2y")

ax[0,1].axis('off')

plt.show()

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[5], line 1
----> 1 fig, ax = plt.subplots(3,2, figsize=(12, 15), facecolor="none")
      3 axis_mesh_field(fig, ax[0,0], mesh, analytic_sol, label="original")
      4 axis_mesh_field(fig, ax[1,0], mesh, stripy_ddx, label="ddy")

NameError: name 'plt' is not defined

from stripy import _srfpack
def second_gradient_local(self, f, index):
    """
    Return the gradient / 2nd partials of an n-dimensional array.

    The method consists of minimizing a quadratic functional Q(G) over
    gradient vectors (in x and y directions), where Q is an approximation
    to the linearized curvature over the triangulation of a C-1 bivariate
    function \\(F(x,y)\\) which interpolates the nodal values and gradients.

    Args:
        f : array of floats, shape (n,)
            field over which to evaluate the gradient
        nit : int (default: 3)
            number of iterations to reach a convergence tolerance,
            tol nit >= 1
        tol: float (default: 1e-3)
            maximum change in gradient between iterations.
            convergence is reached when this condition is met.

    Returns:
        dfdx : array of floats, shape (n,)
            derivative of f in the x direction
        dfdy : array of floats, shape (n,)
            derivative of f in the y direction

    Notes:
        For SIGMA = 0, optimal efficiency was achieved in testing with
        tol = 0, and nit = 3 or 4.

        The restriction of F to an arc of the triangulation is taken to be
        the Hermite interpolatory tension spline defined by the data values
        and tangential gradient components at the endpoints of the arc, and
        Q is the sum over the triangulation arcs, excluding interior
        constraint arcs, of the linearized curvatures of F along the arcs --
        the integrals over the arcs of \\( (d^2 F / dT^2)^2\\), where \\( d^2 F / dT^2\\)is the second
        derivative of \\(F\\) with respect to distance \\(T\\) along the arc.
    """

    if f.size != self.npoints:
        raise ValueError('f should be the same size as mesh')

    sigma = self.sigma
    iflgs = self.iflgs

    f = self._shuffle_field(f)
    index = self._shuffle_simplices(index)

    ## wrapping: 

    # subroutine gradc(k,ncc,lcc,n,x,y,z,list,lptr,lend,dx,dy,dxx,dxy,dyy,ier) ! in :_srfpack:srfpack.f
    # subroutine gradg(  ncc,lcc,n,x,y,z,list,lptr,lend,iflgs,sigma,nit,dgmax,grad,ier) ! in :_srfpack:srfpack.f

    dx, dy, dxx, dxy, dyy, ierr = _srfpack.gradcs(index+1, self._x, self._y, f, self.lst, self.lptr, self.lend)

    if ierr < 0:
        raise ValueError('ierr={} in gradc\n{}'.format(ierr, _ier_codes[ierr]))

    return dx, dy, dxx, dxy, dyy

def gradient_local(self, f, index):
    """
    Return the gradient at a specified node.

    This routine employs a local method, in which values depend only on nearby
    data points, to compute an estimated gradient at a node.

    `gradient_local()` is more efficient than `gradient()` only if it is unnecessary
    to compute gradients at all of the nodes. Both routines have similar accuracy.
    """
    if f.size != self.npoints:
        raise ValueError('f should be the same size as mesh')

    f = self._shuffle_field(f)
    index = self._shuffle_simplices(index)

    gradX, gradY, l = _srfpack.gradls(index + 1, self._x, self._y, f,\
                                     self.lst, self.lptr, self.lend)

    return gradX, gradY

indices = np.arange(0, mesh.npoints)

dx, dy = gradient_local(mesh, analytic_sol, indices)

stripy_ddy.shape, dx.shape

((15523,), (15523,))

fig, ax = plt.subplots(3,2, figsize=(12, 15), facecolor="none")

axis_mesh_field(fig, ax[0,0], mesh, analytic_sol, label="original")
axis_mesh_field(fig, ax[1,0], mesh, dx, label="ddx")
axis_mesh_field(fig, ax[1,1], mesh, dy, label="ddy")
axis_mesh_field(fig, ax[2,0], mesh, stripy_ddx, label="d2x")
axis_mesh_field(fig, ax[2,1], mesh, stripy_ddy, label="d2y")

ax[0,1].axis('off')

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[10], line 1
----> 1 fig, ax = plt.subplots(3,2, figsize=(12, 15), facecolor="none")
      3 axis_mesh_field(fig, ax[0,0], mesh, analytic_sol, label="original")
      4 axis_mesh_field(fig, ax[1,0], mesh, dx, label="ddx")

NameError: name 'plt' is not defined

indices = np.arange(0, mesh.npoints)
dx, dy, dxx, dxy, dyy = second_gradient_local(mesh, analytic_sol, indices)

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[11], line 2
      1 indices = np.arange(0, mesh.npoints)
----> 2 dx, dy, dxx, dxy, dyy = second_gradient_local(mesh, analytic_sol, indices)

Cell In[6], line 43, in second_gradient_local(self, f, index)
     40 if f.size != self.npoints:
     41     raise ValueError('f should be the same size as mesh')
---> 43 sigma = self.sigma
     44 iflgs = self.iflgs
     46 f = self._shuffle_field(f)

AttributeError: 'Triangulation' object has no attribute 'sigma'

d2x = np.zeros_like(stripy_d2x)
d2y = np.zeros_like(stripy_d2y)
d1x = np.zeros_like(stripy_d2y)
d1y = np.zeros_like(stripy_d2y)
for index in range(0, mesh.npoints):
    dx, dy, dxx, dxy, dyy = second_gradient_local(mesh, analytic_sol, index)
    d2x[index] = dxx
    d2y[index] = dyy
    d1x[index] = dx
    d1y[index] = dy

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[12], line 6
      4 d1y = np.zeros_like(stripy_d2y)
      5 for index in range(0, mesh.npoints):
----> 6     dx, dy, dxx, dxy, dyy = second_gradient_local(mesh, analytic_sol, index)
      7     d2x[index] = dxx
      8     d2y[index] = dyy

Cell In[6], line 43, in second_gradient_local(self, f, index)
     40 if f.size != self.npoints:
     41     raise ValueError('f should be the same size as mesh')
---> 43 sigma = self.sigma
     44 iflgs = self.iflgs
     46 f = self._shuffle_field(f)

AttributeError: 'Triangulation' object has no attribute 'sigma'

fig, ax = plt.subplots(3,2, figsize=(12, 15), facecolor="none")

axis_mesh_field(fig, ax[0,0], mesh, analytic_sol, label="original")
axis_mesh_field(fig, ax[1,0], mesh, d1x, label="ddy")
axis_mesh_field(fig, ax[1,1], mesh, d1y, label="ddx")
axis_mesh_field(fig, ax[2,0], mesh, d2x, label="d2x")
axis_mesh_field(fig, ax[2,1], mesh, d2y, label="d2y")

ax[0,1].axis('off')

plt.show()

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[13], line 1
----> 1 fig, ax = plt.subplots(3,2, figsize=(12, 15), facecolor="none")
      3 axis_mesh_field(fig, ax[0,0], mesh, analytic_sol, label="original")
      4 axis_mesh_field(fig, ax[1,0], mesh, d1x, label="ddy")

NameError: name 'plt' is not defined