### Custom Expression Syntax

225

views

0

I want to add a scalar valued function fun(x,y) into the left hand side of my variational formulation:

$\int_{\Omega}fun\left(x,y\right)\cdot grad\left(u\left(x,y\right)\right)\cdot grad\left(v\left(x,y\right)\right)d\left(x,y\right)$∫

This function will later on evaluate a random field which is saved in a numpy Matrix A. My idea is to write a Python function fun(x,y) which implements this evaluation. I understood that this requires calling fun(x,y) in a custom expression, which is then added to the variational formulation.The problem I am facing now is that I am not able to figure out the correct syntax for this.

This should then look similar to this:

$\int_{\Omega}fun\left(x,y\right)\cdot grad\left(u\left(x,y\right)\right)\cdot grad\left(v\left(x,y\right)\right)d\left(x,y\right)$∫

_{Ω}`ƒ``u``n`(`x`,`y`)·`g``r``a``d`(`u`(`x`,`y`))·`g``r``a``d`(`v`(`x`,`y`))`d`(`x`,`y`) .This function will later on evaluate a random field which is saved in a numpy Matrix A. My idea is to write a Python function fun(x,y) which implements this evaluation. I understood that this requires calling fun(x,y) in a custom expression, which is then added to the variational formulation.The problem I am facing now is that I am not able to figure out the correct syntax for this.

This should then look similar to this:

```
from dolfin import *
import numpy
mesh = UnitSquareMesh(8, 8)
A = numpy.identity(8)
def fun(x,y):
return A[floor(8*x),ceil(8*y)]
V = FunctionSpace(mesh, 'Lagrange', 1)
u = TrialFunction(V)
v = TestFunction(V)
class MyExpression(Expression):
def eval(self, x, y):
val[0] = fun(x,y)
x = Expression('x[0]', degree = 1)
y = Expression('x[1]', degree = 1)
my_expr = MyExpression(degree = 1)
a = inner(my_expr(x,y)*nabla_grad(u), nabla_grad(v))*dx
L = Constant(1.0)*v*dx
u = Function(V)
solve(a == L, u, bc)
```

However this gives me the error that it can not evaluate x and y in my_expr(x,y).

If I do not define x and y as expressions beforehand I get the error that x and y are not defined (neither are x[0] and x[1]).

Do you know how I could change the code in order to evaluate the custom expression at the x and y coordinates in the variational formulation?

My idea would be to somehow give the finite element space to MyExpression and get x and y from the current element.

Thank you.

Community: FEniCS Project

### 1 Answer

3

how about this?

```
class MyExpression(Expression):
def eval(self, val, x):
val[0] = A[floor(8*x[0]),ceil(8*x[1])]
```

**edit:**

i can do better, here's working code:

```
from fenics import *
import numpy as np
mesh = UnitSquareMesh(8, 8)
A = np.identity(8)
V = FunctionSpace(mesh, 'Lagrange', 1)
u = TrialFunction(V)
v = TestFunction(V)
class MyExpression(Expression):
def eval(self, val, x):
idx = int(np.floor(8*x[0])) - 1
idy = int(np.ceil(8*x[1])) - 1
val[0] = A[idx, idy]
my_expr = MyExpression(element=V._ufl_element, degree = 1)
a = inner(my_expr*nabla_grad(u), nabla_grad(v))*dx
L = Constant(1.0)*v*dx
u = Function(V)
solve(a == L, u)
```

1

Great that works. Thanks a lot. I will continue working from this.

written
9 months ago by
Christoph Strößner

This sound great and easier than my idea.

Do you also know how I call my_expr in the varational formulation?

`a = inner(my_expr()*nabla_grad(u), nabla_grad(v))*dx`

This does not work, as my_expr requires at least one input argument.

written
9 months ago by
Christoph Strößner

see updated answer above!

written
9 months ago by
pf4d

Thanks a lot. I will continue working from here.

written
9 months ago by
Christoph Strößner

Please login to add an answer/comment or follow this question.