How to implement a varying coefficient in the bilinear form?

variational formulation
1

This question might sound silly, but I have searched the documentation up and down just to find out, that it is very outdated and that there seems to be no explanation and/or solution at all.
I use FEniCS on Ubuntu WSL on Windows 10 and installed it via the Ubuntu-PPA.

dolfin-version
2019.1.0
fenics-version
2019.1.0

Let’s solve the following problem:
Find u\in V, such that
\begin{align*}\int_\Omega \kappa(x, y)\nabla u\cdot\nabla v\,\mathrm{d}\Omega & = 0 & \forall v\in V\\ u(x, y)&= y &\text{on }\Gamma_D\end{align*}
with \kappa being some scalar function depending on x and y.
For simplicity, let’s assume \Omega is the unit square, \Gamma_D is the left side of the domain and \kappa is a piecewise constant:
\kappa(x, y) := \begin{cases}1 & x \leq \frac 1 2 \\ 10 & \text{otherwise}\end{cases}

The idea is, to have an easy variational formulation, such as

a = kappa*inner(grad(u), grad(v))*dx

I think, this should be possible, right?
The implementation should be straightforward, such as this one here Defining subdomains for different materials or this one here Generating meshes with subdomains.

However, the code in these tutorials subclasses the class Expression:

class K(Expression):
    def __init__(self, materials, k_0, k_1, **kwargs):
        self.materials = materials
        self.k_0 = k_0
        self.k_1 = k_1

    def eval_cell(self, values, x, cell):
        if self.materials[cell.index] == 0:
            values[0] = self.k_0
        else:
            values[0] = self.k_1

kappa = K(materials, 1, 10, degree=0)

Trying to do that leads to an error:

/usr/lib/python3/dist-packages/dolfin/function/expression.py in __getattr__(self, name)
    430             return self._parameters
    431         else:
--> 432             return self._parameters[name]
    433 
    434     def __setattr__(self, name, value):

RecursionError: maximum recursion depth exceeded

The reason is apparently, that Expression should not be used anymore as a base class, but we should use UserExpression instead.
This change dates back to 2017, but the documentation was never updated.

So, let’s try to do it with a UserExpression instead.
Here is an MWE - or at least a small example that should actually work (it doesn’t, though):

from fenics import *

mesh = UnitSquareMesh(8, 8)
V = FunctionSpace(mesh, 'P', 1)

class Omega1(SubDomain):
    def inside(self, x, on_boundary):
        return x[0] >= 0.5 - DOLFIN_EPS

materials = MeshFunction('size_t', mesh, mesh.geometric_dimension())
materials.set_all(0)
subdomain1 = Omega1()
subdomain1.mark(materials, 1)

class DCBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 0, DOLFIN_EPS)

left = DCBoundary()

class K(UserExpression):
    def __init__(self, **kwargs):
        self.materials = kwargs['materials']
        self.k0 = kwargs['k0']
        self.k1 = kwargs['k1']
    def eval_cell(self, values, x, cell):
        if self.materials[cell.index] == 0:
            values[0] = self.k0
        else:
            values[0] = self.k1
        
kappa = K(materials=materials, k0=1, k1=10)

u = TrialFunction(V)
v = TestFunction(V)

uD = Expression('x[1]', degree=1)
bc = DirichletBC(V, uD, left)

a = kappa*inner(grad(u), grad(v))*dx
L = Constant(0)*v*dx

u_ = Function(V)
solve(a == L, u_, bc)

The problem is the multiplication with kappa in the definition of the bilinear form.
It results in the error

AttributeError                            Traceback (most recent call last)
<ipython-input-40-b1cfe9512dba> in <module>()
----> 1 a = kappa*inner(grad(u), grad(v))*dx
      2 L = Constant(0)*v*dx

/usr/lib/python3/dist-packages/ufl/exproperators.py in _mul(self, o)
    191         return NotImplemented
    192     o = as_ufl(o)
--> 193     return _mult(self, o)
    194 
    195 

/usr/lib/python3/dist-packages/ufl/exproperators.py in _mult(a, b)
    122     # - matrix-matrix (A*B, M*grad(u)) => A . B
    123     # - matrix-vector (A*v) => A . v
--> 124     s1, s2 = a.ufl_shape, b.ufl_shape
    125     r1, r2 = len(s1), len(s2)
    126 

/usr/lib/python3/dist-packages/ufl/coefficient.py in ufl_shape(self)
     71     def ufl_shape(self):
     72         "Return the associated UFL shape."
---> 73         return self._ufl_shape
     74 
     75     def ufl_function_space(self):

AttributeError: 'K' object has no attribute '_ufl_shape'

Apparently, I have to extend the UserExpression somehow to be able to take part in the multiplication.
Or do I need something completely different?
Some function object or some UFL object or whatsoever?
I searched a lot and did not really find something in the documentation available.

It does work, when defining kappa as an Expression:

kappa = Expression('x[0] <= 0.5 + DOLFIN_EPS ? 1 : 10', degree=0)

But this works only, when the definition of the subdomain is something easy, that can be represented as a formula.
In my work, however, the mesh stems from gmsh and the subdomains are tagged using the physical region mechanism of gmsh.
So i really do need something that works with a MeshFunction.

In fact, my mesh and its subregions are defined as follows:

from dolfin import *
infile_mesh = '/path/to/mesh.xdmf'
infile_boundary = '/path/to/boundary.xdmf'

mesh = Mesh()
mvc = MeshValueCollection("size_t", mesh)

with XDMFFile(infile_mesh) as infile:
    infile.read(mesh)
    infile.read(mvc, "Subdomain")
subdomains = cpp.mesh.MeshFunctionSizet(mesh, mvc)

mvc = MeshValueCollection("size_t", mesh)
with XDMFFile(infile_boundary) as infile:
    infile.read(mvc, "Boundary")
boundaries = cpp.mesh.MeshFunctionSizet(mesh, mvc)

The variable subdomains contains markers for the different regions of the mesh and boundaries contains markers for the different boundaries.

Now, how do I define a varying coefficient, where the value of the coefficient depends on these markers, and where I can use that thing (whatever it is) as factor in the bilinear form?

Best regards
Wolfgang

P.S. I do know, that for a piecewise constant, I can define a custom measure

dx = Measure('dx', domain=mesh, subdomain_data=subdomains)

and then use dx(k) in the definition of the bilinear form.
But this does not work, if kappa is something complicated, because then I have to define a custom measure for every single element, which is definitely not the way to go.

P.P.S. I did read the manuals.
But this is not explained anywhere.
The available documentation on such problems either is totally outdated and uses Expression as a base class, or it uses explicit expressions with easy formulas.
Also, I could not find a single problem in dolfin-demos, with a varying coefficient.
I did a

grep -e '^a = ' --include='*.ufl' -r -- ./dolfin-demos

to inspect all available definitions of bilinear forms.
But - if I’m not mistaken - none of them uses a spatially varying coefficient.

1 Like
2

Hi,
first, try to define the shape of the expression, i.e., use

class K(UserExpression):
def __init__(self, **kwargs):
    self.materials = kwargs['materials']
    self.k0 = kwargs['k0']
    self.k1 = kwargs['k1']
def eval_cell(self, values, x, cell):
    if self.materials[cell.index] == 0:
        values[0] = self.k0
    else:
        values[0] = self.k1
def value_shape(self):
    return (1,)

If this does not work or change anything, maybe try interpolating / projecting the expression to an appropriate function space and then use the obtained function in your weak form.

3

Unfortunately this gives the same AttributeError as above.

Apparently, the following works:

class K(UserExpression):
    def __init__(self, **kwargs):
        super().__init__(kwargs)
        self.materials = kwargs['materials']
        self.k0 = kwargs['k0']
        self.k1 = kwargs['k1']
    def eval_cell(self, values, x, cell):
        if self.materials[cell.index] == 0:
            values[0] = self.k0
        else:
            values[0] = self.k1
    def value_shape(self):
        return ()
            
kappa = K(materials=materials, k0=1, k1=20, element=V.ufl_element())
# The following works, too
# kappa = K(materials=materials, k0=1, k1=20, degree=0)

I have to return an empty list in value_shape, as otherwise I get the error message:

Can only integrate scalar expressions. The integrand is a tensor expression with value shape (1,) and free indices with labels ().

However, there is one caveat left:
If I change one of the parameters k0 or k1 and call solve, the underlying linear system still contains the old values.

kappa = K(materials=materials, k0=1, k1=10, element=V.ufl_element())
solve(a == L, u_, bcs) # solves for old values of k0 and k1

If I want to have different values, I have to reassign the variable of the bilinear form again:

kappa = K(materials=materials, k0=1, k1=10, element=V.ufl_element())
a = kappa*inner(grad(u), grad(v))*dx
solve(a == L, u_, bcs) # solves for new values of k0 and k1

This is probably intended behavior, as the solver routines depend on a only, but a has no knowledge anymore, that it depends on kappa.
But I leave this here as a reference for other people with the same problem, so they are aware of this.

Anyway: Is this the recommended way to implement this?
Are there any subtleties one has to take into account?

4

UserExpression assumes that, if you override the constructor, you will call its superclass constructor. The following modification runs:

class K(UserExpression):
    def __init__(self, materials, k0, k1, **kwargs):
        
        # Call superclass constructor with keyword arguments to properly
        # set up the instance.  It will get confused, though, if you
        # include your own keyword arguments in kwargs.  Since the body
        # of the function assumes these arguments will be passed, there
        # doesn't seem to be a good reason not to make them positional
        # arguments.
        super().__init__(**kwargs)
        
        self.materials = materials #kwargs['materials']
        self.k0 = k0 #kwargs['k0']
        self.k1 = k1 #kwargs['k1']
    def eval_cell(self, values, x, cell):
        if self.materials[cell.index] == 0:
            values[0] = self.k0
        else:
            values[0] = self.k1
    def value_shape(self):
        return ()
        
#kappa = K(materials=materials, k0=1, k1=10)
kappa = K(materials,1,10)
1 Like