[firedrake] Interpolating data to Function

[Lawrence Mitchell]

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.
> 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.

Cheers,

Lawrence