Specifying the source term in poissons equation from another solution

I have the following set of boundary value problems over the same domain:

  1. Solution U, to Laplaces equation (+ boundary conditions)
  2. I then want a solution to Poissons equation with a source term f=grad(U).grad(U) (+ boundary conditions

My question concerns the specification of f. In solving U I have the function space

V = fem.functionspace(domain, ("Lagrange", 1))

I think the gradient and f should be calculated via

V_grad = fem.functionspace(domain, ("Lagrange", 1, (domain.geometry.dim,)))
grad_U = ufl.grad(U)
f=ufl.inner(grad_U,grad_U)

What I am not overly sure about is the form/type that f has.

Well I am at it - how do I get the array of points that goes with f?

f is a symbolic expression in your code (ufl expression). When you call dolfinx.fem.form or assemble_vector this is turned in to an evaluation of the symbolic expression at the quadrature points of your rhs integral and accumulated into the global vector, doing exactly what you set out to do.

You do not need the array to solve the problem. However if you would like to inspect f, you cal use dolfinx.fem.Expression:

f_compiled = dolfinx.fem.Expression(f, points_in reference_element)
Values = f_compiled.eval(mesh, np.arange(num_cells_local, dtype=np.int32))

As discussed in

What I ended up doing

n=ufl.FacetNormal(macor_mesh)
ds = ufl.Measure("ds", domain=macor_mesh, subdomain_data=facet_marker)
test=np.array([fem.assemble_scalar(fem.form(ufl.inner(grad_u, n) * ds(k))) for k in range(1,7)])
print(test.sum())
test

In this problem I have multiple boundarys (electrodes inside an enclosure - with one set to 1 and all others including the outer boundary set to zero) and solving laplaces equation. Integrals over the boundaries represent power into or out of the boundary (or total charge for emag). So the total has to be zero. It tends to zero but only if the number of elements becomes prohibitively large. I also wonder what interpolation means in this context. Pretty sure it means evaluating a function at a set of points and not interpolating tabulated data at other coordinate points - at least that’s what the tutorials would suggest.

Just to add - Integration measures — FEniCS Tutorial @ Sorbonne gives an explicit section on doing boundary integrals but doesn’t actually evaluate anything - so when I say there are no good tutorials this is exactly what I mean

Please note that this is an open-source code, and the work done for tutorials is done by people that likes to contribute to the community.

There are demos on boundary integrals in:

and

As evaluating a boundary integral doesn’t really differ from evaluating a cell integral (and summing the contributions into a scalar), I did not include another line of code doing that in:

For this problem, you should really look at:

as it exactly talks about how one should evaluate fluxes for finite element problems.