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[0] = eps/pi/(x[0]**2 + x[1]**2 + eps**2)
C = FiniteElement("Lagrange", mesh.ufl_cell(), degree = 2)
C = FunctionSpace(mesh, C)
zerotop_concentration = DirichletBC(C, 0, "on_boundary && x[1] > 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[0] = eps/pi/(x[0]**2 + x[1]**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[0] = sourceExpression
Of course this does not work because value[0] is defined directly as x[0], and not as Expression('x[0]', 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[0] = self.sourceExpression(x)
def value_shape(self):
return ()
# Test:
sourceExpression = Expression('x[0]', 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[0]**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[0] = (1/(2 * pi)) * (param / ( (x[0] - 0.5)**2 + (x[1] - 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[0] = (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] - 0.5)**2 + (x[1] - 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!