Interpolate numpy ndarray to Function

1

I have an analytical solution for a problem given as a numpy ndarray and want to do some error analysis. The analytical solution is not explicit in (x,y) and is expensive to compute so I can’t really use something like UserExpression for example. Instead I have three arrays, all N x N, representing the solution values and the corresponding x and y coordinates (non-uniform), mocked up below. What’s a good way to interpolate the solution values onto a FunctionSpace? The complication is that I want to minimize interpolation errors as I want to use this for a convergence analysis. In the example below I use a KDTree to lookup the closest known solution for each mesh coordinate. Would this be the right approach?

from dolfin import *
import numpy as np
import matplotlib.pyplot as plt
from scipy import spatial

x = np.linspace(-5, 5, 100)
y = np.linspace(-10, 10, 100)
X, Y = np.meshgrid(x, y)
# dummy func, note my function is not explicit in x, y
Z = np.sin(X*2 +Y**2) / (X**2 + Y**2)
plt.contourf(X, Y, Z)

mesh = RectangleMesh(Point(-5,-10), Point(5,10), 20, 20)
V = FunctionSpace(mesh, 'CG', 1)
N = mesh.geometry().dim()
coord = V.tabulate_dof_coordinates().reshape(V.dim(), N)

# Shoe-horn function coordinate data for entry into KDTree routine
combined_X_Y_arrays = np.dstack([X.ravel(), Y.ravel()])[0]
coord_list = list(coord)

mytree = spatial.KDTree(combined_X_Y_arrays)
dist, indexes = mytree.query(coord_list)
...
2

Would it be too slow to wrap a method from scipy.interpolate in a UserExpression?

3

Hi - Do you mean something like this? I don’t think speed would be an issue here. Is this safe to use for a convergence study? I’m getting nans when I call errornorm in the code below. Is the UserExpression evaluated at the quadrature points when I call errornorm?

from dolfin import *
import numpy as np
import matplotlib.pyplot as plt
from scipy import spatial
import scipy.interpolate

x = np.linspace(-5, 5, 100)
y = np.linspace(-10, 10, 100)
X, Y = np.meshgrid(x, y)
# dummy func, note my function is not explicit in x, y
Z = np.sin(X**2 + Y**2) / (X**2 + Y**2)
# fig, ax = plt.subplots()
# plt.contourf(X, Y, Z)
# plt.axis('equal')
f = scipy.interpolate.LinearNDInterpolator(list(zip(X.ravel(), Y.ravel())), Z.ravel())

class Solution(UserExpression):
    def __init__(self, degree=2, **kwargs):
        super().__init__(**kwargs)
        
    def eval(self, values, x):
        values[0] = f(x[0], x[1])

    def value_shape(self):
        return ()  

soln = Solution()
for N in [20, 40, 80]: # try for different mesh resolutions
    mesh = RectangleMesh(Point(-5, -10), Point(5, 10), N, N)
    V = FunctionSpace(mesh, 'CG', 1)
    z = project(soln, V)
    z_ = project(Expression(('sin(pow(x[0], 2) + pow(x[1], 2)) / (pow(x[0], 2) + pow(x[1], 2))'), degree=2), V)
    e = errornorm(z_, z, degree_rise=4)
    print(e)  # why nan?
4

The nan comes from a divide-by-zero in z_, since there is a node at \mathbf{x}=\mathbf{0}. If you use a different expression, the error converges until saturating on the resolution of the 100\times 100 grid that is being interpolated from.

5

Of course, thank you! :sweat_smile: