[firedrake] Interpolating data to Function

[Tuomas Karna]

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

Cheers,

Tuomas

> Cheers,
>
> Lawrence