Hi Tuomas, On Wednesday, July 30, 2014, Tuomas Karna <tuomas.karna@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.

This has been on the wish list for a while and you prompted us to do something about it today. Expect a merge with a better solution tomorrow our time.

...

...

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?

That's not straightforward, but the other way around is easy. You can just access x_func.dat.data_with_halos . It shouldn't hurt too much to do the duplicate interpolation on the halo nodes. Cheers, David

Cheers,

Tuomas

...

Cheers,

Lawrence

_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk <javascript:;> https://mailman.ic.ac.uk/mailman/listinfo/firedrake

-- Dr David Ham Departments of Mathematics and Computing Imperial College London http://www.imperial.ac.uk/people/david.ham