Interpolating matrix data to mesh

variational formulation
1

Hi! I have two matrices, say X,Y where i have coordinates of x,y, and i have a matrix Z where i save a value for each coordinate (the matrices have the same dimension). So, in Matlab with this data i create an interpolator (matlab function created with scatteredInterpolant) and evaluate this on my mesh nodes to make the function (or vector in Matlab) i need. How can i do this in Fenics? (I already have the mesh defined)
Thanks for the help!

2

Consider using the package scipy.interpolate to accomplish this. Particularly the functions bisplrep and bisplev are useful. Take the following example for assigning a scalar value of Temperature at different timesteps. Column 0 and 1 of the txt file are the coordinates of the prescribed domain and each column afterwards the corresponding temperature at different timesteps

from dolfin import *
import numpy as np
from scipy.interpolate import bisplrep
from scipy.interpolate import bisplev

data = np.loadtxt("test.txt")
x, y = data[:,0], data[:,1]
T = data[:,2:]

class ExpressionFromScipyFunction(UserExpression):
    def __init__(self, x, y, f, t, **kwargs):
        self._f = f
        self._x = x
        self._y = y
        self._t = t
        UserExpression.__init__(self, **kwargs)
    def eval(self, values, x):
        self.interp = bisplrep(self._x, self._y, self._f[:,self._t])
        values[:] = bisplev(x[0], x[1], self.interp, dx = 0, dy = 0)
    def values_shape(self):
        return (1,)

mesh = RectangleMesh(Point(0,0),Point(0.0255,0.0255),52,52)
V = FunctionSpace(mesh, 'CG', 1)

expression = ExpressionFromScipyFunction(x, y, T, 0, element=V.ufl_element())

expression_out = interpolate(expression, V)
expression._t = 1 # Next Time Step

The above solution works for scalar functions. More work would be needed to expand to prescribed vector values for your coordinates. The mesh used in the example corresponds to the domain the imported data comes from.

1 Like
3

What’s the difference if i project

h = project(expression, V)

instead of

h = interpolate(expression, V)
4

I’ll let the experts answer this: What is the difference between interpolate() and project() ? and Project(x) works, interpolate(x) does not. I am using a UserExpression with values presribed at the dofs so interpolate works well for me.

3 Likes
5

This works, but it takes too long (my data dimension is 250000) i tried with griddata (maybe would work faster i’m not sure) and i get the error:

class ExpressionFromScipyFunction(UserExpression):
    def __init__(self, x, y, f, t, **kwargs):
        self._f = f
        self._x = x
        self._y = y
        self._t = t
        UserExpression.__init__(self, **kwargs)
    def eval(self, values, x):
        values[:] = griddata((self.x,self.y),self._f[:,self._t],(x[0],x[1]))
        #self.interp = bisplrep(self._x, self._y, self._f[:,self._t])
        #values[:] = bisplev(x[0], x[1], self.interp, dx = 0, dy = 0)
    def values_shape(self):
        return (1,)

AttributeError: 'ExpressionFromScipyFunction' object has no attribute 'x'
6

Ivan,

The error you get here is part of the class definition. The variables are self._x and self._y.

values[:] = griddata((self._x,self._y),self._f[:,self._t],(x[0],x[1]))

The problem with griddata is it returns an array of interpolated values at least from my understanding and you would be calling a whole new interpolation for every point in your data. I realize now that what my previous code also suffers from this:

class ExpressionFromScipyFunction(UserExpression):
    def __init__(self, x, y, f, t, **kwargs):
        self._f = f
        self._x = x
        self._y = y
        self._t = t
        UserExpression.__init__(self, **kwargs)
    def eval(self, values, x):
        values[:] = bisplev(x[0], x[1], self.interp, dx = 0, dy = 0)
    def UpdateInterp(self):
        self.interp = bisplrep(self._x, self._y, self._f[:,self._t])
    def values_shape(self):
        return (1,)

Consider trying this, where you call UpdateInterp after instantiating the class.

expression = ExpressionFromScipyFunction(x, y, T, 0, element=V.ufl_element())
expression.UpdateInterp()
expression_out = interpolate(expression, V)
expression._t = 1 # Next Time Step

For me this speeds up the interpolation incredibly since it’s not remaking the interpolation everytime. I can’t speak for griddata but consider trying this approach.

1 Like
7

This solved it! Thanks a lot Sven! I was running the interpolation for 5 days and now it took 3 min.