[firedrake] Interpolating data to Function

[Lawrence Mitchell]

> On 30 Jul 2014, at 21:21, Tuomas Karna <tuomas.karna at gmail.com> wrote:
> 
> Hi Lawrence,
> 
>> On 07/30/2014 02:13 AM, Lawrence Mitchell wrote:
>>> On 30/07/14 02:56, Tuomas Karna wrote:
>>> Hi Firedrake-people,
>>> 
>>> I'm working with simple 2D ocean model and I'd like to load bathymetric
>>> data to a Function, interpolating to the correct coordinates of the
>>> nodes. Currently I get the x/y coordinates from specialized x/y function
>>> and then just write the data directly to Function.dat.data array. So
>>> something like this:
>>> 
>>> P1 = FunctionSpace(mesh, "CG", 1)
>>> bath = Function(P1,name='bathymetry')
>>> x_func = Function(P1).interpolate(Expression('x[0]'))
>>> y_func = Function(P1).interpolate(Expression('x[1]'))
>>> def interpolateBath(x,y):
>>>     # interpolate here
>>>     return 3.0
>>> bath.dat.data[:] = interpolateBath(x_func.dat.data, y_func.dat.data)
>>> 
>>> Is there a cleaner way of doing this?
>> I think there is unfortunately currently not.
> OK, just checking if I was missing something obvious. 

Note that we have a branch in the works that allows you to define a python class to be used in function.interpolate (this is the same as DOLFIN's python expressions). So this may well end up doing what you want. You'll need an updated version of pyop2 and the python_expressions firedrake branch if you want to try it. 

>>> I have a similar question for the boundaries. I have external data
>>> defined on a line which I'd like to interpolate to the nodes of a
>>> certain boundary. I could use a Function for this too, if I knew the
>>> indices of the boundary nodes. I've hacked something together using
>>> mesh.exterior_facets and FunctionSpace.exterior_facet_boundary_node_map,
>>> just wondering if there is an easier way?
>> If your mesh has a marker for that boundary, (e.g. a boundary id in
>> Gmsh), then you can get the indices of the boundary nodes (by which I
>> presume you mean the vertices) with:
>> 
>> mesh = Mesh("input.msh")
>> 
>> V = FunctionSpace(mesh, 'CG', 1)
>> 
>> bc = DirichletBC(V, 0.0, integer_id_of_boundary)
>> 
>> boundary_nodes = bc.nodes
>> 
>> 
>> under the hood, this likely does something very similar to what you did
>> with exterior_facets and the exterior_facet_boundary_node_map.
> Thanks, I have boundary markers so bc.nodes works fine. In parallel however it also returns indices that exceed the size of x_func.dat.data, maybe the halo nodes are included? Is there an easy way to get the resident (non-halo) nodes only?

You can do one of two things:

1. Restrict the bc nodes using numpy.take. You can find out the number of local degrees of freedom in a function in f.dof_dset.size

2. Just access the day with the halo values in place:

x_func.dat.data_with_halos

I think the second option is probably simpler

Lawrence