Post-processing: integrate on boundary product of functions

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

    n = FacetNormal(mesh)
    flux = -k*dot(grad(u), n)*ds
    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

1 Answer

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  
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  
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:

For better performance, you can also put more complicated C++ code in the string that is passed to the standard Expression constructor, but I have less experience with that (and whatever caveats might be involved).  
written 5 months ago by David Kamensky  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »