Initial condition from a regular Real Data grid

mesh
1

Hi everyone,
Does anyone know how to use real 2D data (on regular lat/lon grid) as initial condition in a mesh like this one?
thanks!!

p.s: To generate that figure I used the following code:

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

data = np.loadtxt('test.csv')
x,y,values = data[:,0], data[:,1], data[:,2]
interpolant = interp2d(x,y, values, kind='linear') 

class MyExpression1(Expression):
    def __init__(self, dado,*args, **kwargs):
        self.dado = dado
    def eval(self, values, x):
        values[:] = self.dado(*x)

        if(values[:]>5):
        	values[:] = 5
        if(values[:] < -2):
        	values[:] = -2

mesh = Mesh("test.xml")
expression = MyExpression1(interpolant, degree=0, domain=mesh)

V = FunctionSpace(mesh, expression.ufl_element())
expression = interpolate(expression,V)
r = plot(expression)
plt.colorbar(r)
plot(mesh)                                                                     
plt.show()

with the test.csv looking like:

6.061990000000000000e+05 -8.039260000000000000e+05 3
-2.800670000000000000e+05 -1.885540000000000000e+06 2
4.814150000000000000e+05 -8.780600000000000000e+05 3
-1.055750000000000000e+06 -3.797740000000000000e+05 -1
3.794040000000000000e+05 -1.073310000000000000e+06 2
-9.699220000000000000e+05 -1.430850000000000000e+06 1
4.433700000000000000e+05 -1.006770000000000000e+06 0

But to use a real 2D data I’m not able to figure out how to regrid regular grid to unstructured mesh.

2

After

a_func = interpolate(expression, V)

Then a_func should be a dolfin.Function in the V FunctionSpace, which is normally what you would want as an initial value, or am I misunderstanding your question?

4

Why can’t you just use interp2d as you showed above? You will have to do it on one core and then save the result and reload it using multiple cores if you are running in parallel, but that shouldn’t be a problem?

5

@tormod, in that case I used coords = mesh.coordinates() to get vertex coordinates, then I generated a random data (for coords) between (-2,5) only to test it, but how can I do it if my input data has not the same meshe’s coordinates ?

6

I’m not sure I understand your question. There is no mesh.coords() anywhere in your code. Also, looking at the code the arguments degree=0, domain=mesh are never used.

To answer your question: FEniCS will query your Expression with the positions x that it needs to perform the interpolation. You do not need to specify these locations yourself, FEniCS knows which points it needs (not the same locations as the vertex coordinates in cases other than CG1 functions). SciPy’s interp2d is able to compute the interpolated value at any point inside the regular grid, so that should be unproblematic as long as the FEniCS mesh is contained inside the regular lat-lon grid