Point evaluation on immersed manifold

[edg-13]

I am trying to evaluate points on a 2D surface in a 3D space. As I understand it, I cannot do this using the Function.at() capability.

I have been able to extract an array of points using Function.vector() but I require the coordinates of the points that the data corresponds to. Is there a definitive way to find a mapping between coordinate and function value in this case?

Note that I have successfully read in the .vtu files with a mapping between coordinate and function value with third party software, but ideally I need the point evaluation during the simulation.

Any suggestions would be great!

[dham]

This is supposed to be fixed if you use the karai branch of Firedrake and use the VertexOnlyMesh approach to point evaluation rather than .at().

[edg-13]

Following from this, is there a way to set a Function in a similar manner? For example, Function.set_at(coord, new_value)..

[dham]

I don't understand what you are asking for. Functions in finite element spaces don't have values that you can set in a pointwise manner. They are the weighted sum of basis functions. What is it you are attempting to achieve?

[edg-13]

I am trying to get a source term which is 0 on most of the mesh but only non zero on some nodes. The mesh has a complicated surface so I am not certain whether it is feasible to define a function using spatial coordinates.

Edit: Upon further thought I would be happy to approximate this delta-like distribution as a sharp Gaussian which is expressible as a sum of basis functions. Thanks for the clarification.

[dham]

Point forcing is the adjoint (transpose) operation to interpolation. You can get this in Firedrake at the moment using a slightly lower level interface to solve your PDE. Here's an example of Poisson forced by two point sources:

from firedrake import *
import matplotlib.pyplot as plt

m = UnitSquareMesh(10, 10)

V = FunctionSpace(m, "CG", 1)
u = TrialFunction(V)
v = TestFunction(V)
sol = Function(V)

vertices = VertexOnlyMesh(m, [[0.3, 0.3], [0.4, 0.4]])
V_p = FunctionSpace(vertices, "DG", 0)
point_vals = Function(V_p)
point_vals.interpolate(Constant(1.0))

bc = DirichletBC(V, 0.0, "on_boundary")

poisson_lhs = assemble(inner(grad(u), grad(v))*dx, bcs=bc)
poisson_rhs = Interpolator(v, V_p).interpolate(point_vals, transpose=True)

solve(poisson_lhs, sol, poisson_rhs)
fig, axes = plt.subplots()
tripcolor(sol, axes=axes)
axes.legend()
plt.show()