[firedrake] Interpolating data to Function

[David Ham]

Hi Tuomas,

On Wednesday, July 30, 2014, 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.


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 at 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
-------------- next part --------------
HTML attachment scrubbed and removed