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