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.
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. Cheers, Lawrence