[firedrake] Interpolating data to Function

[David Ham]

Hi Tuomas,

The nicer Python interpolation support is now live. The documentation is
at:
http://www.firedrakeproject.org/firedrake.html#module-firedrake.expression

Hope that's useful.

Regards,

David


On 30 July 2014 21:45, Tuomas Karna <tuomas.karna at gmail.com> wrote:

>
> On 07/30/2014 01:34 PM, Lawrence Mitchell wrote:
> >> 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.
> Sounds good, I'll take a look.
> >
> >>>> 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
> OK, both approaches seem to work fine.
>
> Thanks,
>
> Tuomas
>
> >
> > Lawrence
> >
>
>
> _______________________________________________
> firedrake mailing list
> firedrake at imperial.ac.uk
> 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