 # How to define source term function

Hi,

I need to use a parametric source function (an approximative delta) and I tried to do the following, where u is a given and defined function already.

``````class nascentDelta(Expression):
def __init__(self, eps):
self.eps = eps
def eval(self, values, x):
eps = self.eps
values = eps/pi/(x**2 + x**2 + eps**2)

C = FiniteElement("Lagrange", mesh.ufl_cell(), degree = 2)
C = FunctionSpace(mesh, C)
zerotop_concentration = DirichletBC(C, 0, "on_boundary && x > 1 - DOLFIN_EPS")
bcc = [zerotop_concentration]
c = TrialFunction(C)
d = TestFunction(C)

a = inner(u, grad(c)) * d * dx + (c.dx(0) * d.dx(0) + c.dx(1) * d.dx(1)) * dx - inner(c.dx(1), d.dx(1)) * ds(1)
delta = nascentDelta(0.01)
L = inner(delta, d) * dx

c = Function(C)
solve(a == L, c, bcc)
``````

But this gives me the following recursion error:

File “/usr/local/lib/python3.6/dist-packages/dolfin/function/expression.py”, line 432, in getattr
return self._parameters[name]
[Previous line repeated 327 more times]
RecursionError: maximum recursion depth exceeded while calling a Python object

This is based on the idea detailed at https://fenicsproject.org/qa/7941/applying-dirac-function-using-fenics-pointsource-function/?show=7941#q7941

Could anyone help me and tell what I am doing wrong? Thank you so much.

Custom Python definitions of expressions should be subclasses of `UserExpression`, not `Expression` in the current version of FEniCS (2019.1). See, e.g., the following example:

``````from dolfin import *

# Minimal example reproducing error from post:

#class e(Expression):
#    def __init__(self, eps):
#        self.eps = eps
#eInstance = e(1.0)

# In version 2019.1, you want to subclass UserExpression:
class e(UserExpression):
def __init__(self,eps,**kwargs):
# Call superclass constructor with keyword arguments to properly
# set up the instance:
super().__init__(**kwargs)
# Perform custom setup tasks for the subclass after that:
self.eps = eps

def eval(self, values, x):
eps = self.eps
values = eps/pi/(x**2 + x**2 + eps**2)

def value_shape(self):
return ()

# Test:

# NOTE: UserExpressions are always interpolated in finite element spaces during
# assembly of forms.  You can control which space by passing a degree
# keyword argument.  The default value, if nothing is passed, is degree=2.
# For the given regularized delta function, it will always be an approximation,
# since the formula is non-polynomial.  (An alternative would be to define
# it directly in UFL, using `x = SpatialCoordinate(mesh)`.)
eInstance = e(1.0,degree=4)

mesh = UnitSquareMesh(1,1)
print(assemble(eInstance*dx(domain=mesh)))
``````
1 Like

Would it be somehow possible in the above example to receive an Expression argument in `e`, so it would be

`sourceExpression = Expression('0.0', degree = 2)`

and in the source term function:

`def __init__(self,eps, sourceExpression, **kwargs):`
`self.sourceExpression = sourceExpression`

and in the eval definition I would like to use this argument as

`values = sourceExpression`

Of course this does not work because `value` is defined directly as `x`, and not as `Expression('x', degree = 2)`.

Is there a way to set its value to an argument instead of a fix definition?

Thank you!

You could pass `x` through to `sourceExpression`, as follows:

``````from dolfin import *

class e(UserExpression):
def __init__(self,eps,sourceExpression,**kwargs):
super().__init__(**kwargs)
self.eps = eps
self.sourceExpression = sourceExpression
def eval(self, values, x):
values = self.sourceExpression(x)
def value_shape(self):
return ()

# Test:
sourceExpression = Expression('x', degree=2)
eInstance = e(1.0, sourceExpression, degree=2)
mesh = UnitSquareMesh(1,1)
print(assemble(eInstance*dx(domain=mesh)))
``````

However, in the above example, it’s not clear what the advantage would be over just using `sourceExpression` directly.

1 Like

Thank you! Using

``````sourceExpr = Expression('x**1.5', degree = 10)
``````

gives me

raise RuntimeError(“Unable to compile C++ code with dijitso”)
RuntimeError: Unable to compile C++ code with dijitso

Have you dealt with something like this before? I have tried to use some ideas suggested on other online threads, none of them worked.

The string passed to the `Expression` constructor needs to be C++ code. C++ does not have a `**` operator. You can instead use the `pow` function from the C math library, e.g., a^b = `pow(a,b)`.

1 Like

If possible, I would like to ask for help / advice on using parametric expressions in the following kind of data structures. I am running and testing the solver in different scenarios, so I define ProblemData as a collection of data i.e. boundary conditions, parameters and forcing functions. Then I create instances of such objects and pass these one by one to the solver. For example we have these two problem data:

``````**#problemdata0**
f0 = Expression('0.0', degree = 2)
bcs_list0 = []

class f_c_pd0(UserExpression):

def __init__(self,param,**kwargs):
super().__init__(**kwargs)
self.param = param

def eval(self, values, x):
param = self.param
values = (1/(2 * pi)) * (param / ( (x - 0.5)**2 + (x - 0.5)**2 + param**2 )**(1.5) )

def value_shape(self):
return ()

problem_data0 = ProblemData(f0, bcs_list0, f_c_pd0)

**#problemdata1**
f1 = Expression('1.0', degree = 2)
bcs_list1 = []

class f_c_pd1(UserExpression):

def __init__(self,param,**kwargs):
super().__init__(**kwargs)
self.param = param

def eval(self, values, x):
param = self.param
values = (1/(2 * pi)) * param

def value_shape(self):
return ()

problem_data1 = ProblemData(f1, bcs_list1, f_c_pd1)
``````

Then we have a `solver` function, where the `f_c` forcing term gets its parameter defined:

``````def solver(problemdata):
param = 1.0
source_instance = problemdata.f_c_pd(param)
L = inner(source_instance, testfunction) * dx
``````

The code I have does work, but defining a whole new `class f_c_pd_i(UserExpression)` for every single problemdata seems odd, given that it is more than 10 lines and only one line changes in their definition.
Is there a more efficient way to define a parametric source term in a problemdata object in a way that we can set that parameter in runtime and then create an instance of the source term and use it in the weak form as `L = inner(source_instance, testfunction) * dx` ?

Thank you so much for any ideas!

It seems to me like, at least in the given examples, the standard `Expression` class can be used. Note that you can pass extra parameters as keyword arguments to the `Expression` constructor. For example, you could do something like the following:

``````param = 1.0
source_instance = Expression("(1/(2 * pi)) * param",degree=0,param=param)
``````

and then update `param` later, like

``````source_instance.param = 2.0
``````

The `UserExpression` class is mainly for terms involving complicated custom logic that cannot be implemented through `Expression`.

Further, I prefer to just define spatially-varying terms in UFL where I can get away with it; `Expression` is mainly needed if the expression will be interpolated or used as Dirichlet data. For a source term that is integrated in a variational form, you can do something like

``````x = SpatialCoordinate(mesh)
param = Constant(1.0)
source_instance = (1/(2 * pi)) * (param / ( (x - 0.5)**2 + (x - 0.5)**2 + param**2 )**(1.5) )
``````

which has the advantage of not depending on the choice of `degree` in the `Expression` or `UserExpression`.

Thank you! I tried out your first suggestion, I think it works great.

I appreciate your second scomment, but the `mesh` is unknown (and thus `x = SpatialCoordinate(mesh)` as well) at the point where `problemdata` is initialised, since for each `problemdata` the solver is run for a loop of finer and finer meshes. I could have an update() function that updates `mesh` and `x = SpatialCoordinate(mesh)` within the `problemdata` and redefines `source_instance` with the problemdata-specific definition, but then I more or less get back to the original point where I had a separate class definition for each source (if I understand well?)

Thank you again!