Custom Expression Syntax


225
views
0
9 months ago by
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)$Ωƒ un(x,y)·grad(u(x,y))·grad(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
9 months ago by
pf4d  
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.

Similar posts:
Search »