### Post-processing: integrate on boundary product of functions

173
views
0
5 months ago by
From chapter 5 of "solving PDEs in Python", I get that I can compute a flux

n = FacetNormal(mesh)
total_flux = assemble(flux)

In the aforementioned example, k is a constant. How can I extend this to the case where k is a function as such

def k(x):
# big complicated computation
return result
Community: FEniCS Project

5
5 months ago by
If $k$ is strictly a function of position, say, $k(\mathbf{x}) = x_0^2 + x_1^2 + 1$, you could define it like
def k(x):
return x[0]**2 + x[1]**2 + 1.0​
and then use it directly to define the flux, passing the spatial coordinates as an argument:
x = SpatialCoordinate(mesh)
flux = -k(x)*dot(grad(u),n)*ds​

The assembly would then use the same line of code.

An additional remark, though, on flux extraction:  If this flux is on a Dirichlet boundary where you are prescribing $u$, you might look into conservative flux extraction instead of direct evaluation.  See, e.g., my post here:

and references cited therein.

x = SpatialCoordinate(mesh) raises UFLException: Invalid cell <Mesh of topological dimension 2 (triangles) with 2507 vertices and 4847 cells, ordered>. Are there any caveats to using this class, or am I doing something wrong?
written 5 months ago by Isabelle Santos
1
I've never run into that error before; maybe it's something specific to your mesh.  Another alternative to try (again assuming $k$ is a function of position) would be to just define $k$ as an Expression, e.g.,

k = Expression("x[0]*x[0] + x[1]*x[1] + 1.0",degree=2)​
and then use it in the definition of the flux without passing any argument, i.e.,
flux = -k*dot(grad(u),n)*ds​​
written 5 months ago by David Kamensky
Can I use loops and conditionals in an Expression ?
written 5 months ago by Isabelle Santos
2
Yes, you can create a subclass of Expression and override the eval() method to execute arbitrary Python code.  For an example of this, see my answer here: