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
https://mailman.ic.ac.uk/mailman/listinfo/firedrake
--