# How to implement a varying coefficient in the bilinear form?

**URL:** <https://fenicsproject.discourse.group/t/how-to-implement-a-varying-coefficient-in-the-bilinear-form/1988>\
**Category:** variational formulation\
**Created:** [November 26, 2019, 12:25pm UTC](https://fenicsproject.discourse.group/t/how-to-implement-a-varying-coefficient-in-the-bilinear-form/1988 "2019-11-26T12:25:19Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![hoelzlw](https://avatars.discourse-cdn.com/v4/letter/h/43a26b/32.png) [@hoelzlw](https://fenicsproject.discourse.group/u/hoelzlw)\
**Post date:** [November 26, 2019, 12:25pm UTC](https://fenicsproject.discourse.group/t/how-to-implement-a-varying-coefficient-in-the-bilinear-form/1988/1 "2019-11-26T12:25:19Z")

</div>

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.

```auto
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

```python
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](https://fenicsproject.org/pub/tutorial/html/._ftut1013.html#ftut:possion:2D:2mat:impl) or this one here [Generating meshes with subdomains](https://fenicsproject.org/pub/tutorial/sphinx1/._ftut1005.html#generating-meshes-with-subdomains).

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

```python
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:

```nohighlight
/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):

```python
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

```nohighlight
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`:

```python
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:

```python
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

```python
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

```auto
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.

---

<div class="post-metadata">

**Author:** ![plugged](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/plugged/32/5115_2.png) [@plugged](https://fenicsproject.discourse.group/u/plugged)\
**Post date:** [November 26, 2019, 12:45pm UTC](https://fenicsproject.discourse.group/t/how-to-implement-a-varying-coefficient-in-the-bilinear-form/1988/2 "2019-11-26T12:45:39Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![hoelzlw](https://avatars.discourse-cdn.com/v4/letter/h/43a26b/32.png) [@hoelzlw](https://fenicsproject.discourse.group/u/hoelzlw)\
**Post date:** [November 26, 2019, 1:24pm UTC](https://fenicsproject.discourse.group/t/how-to-implement-a-varying-coefficient-in-the-bilinear-form/1988/3 "2019-11-26T13:24:42Z")

</div>

Unfortunately this gives the same `AttributeError` as above.

Apparently, the following works:

```python
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:

```nohighlight
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.

```python
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:

```python
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?

---

<div class="post-metadata">

**Author:** ![kamensky](https://avatars.discourse-cdn.com/v4/letter/k/e95f7d/32.png) [@kamensky](https://fenicsproject.discourse.group/u/kamensky)\
**Post date:** [November 26, 2019, 2:42pm UTC](https://fenicsproject.discourse.group/t/how-to-implement-a-varying-coefficient-in-the-bilinear-form/1988/4 "2019-11-26T14:42:57Z")

</div>

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

```python
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)

```