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

173

views

0

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

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

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:

https://www.allanswered.com/post/qvzpa/#rezmw

and references cited therein.

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:

https://www.allanswered.com/post/klzvn/#jrzer

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).

https://www.allanswered.com/post/klzvn/#jrzer

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.

x = SpatialCoordinate(mesh)raisesUFLException: 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?