How do I implement a custom operator involving derivatives?

Question

Hi,

So I am trying to reimplement a simulation of piezoelectric discs in FEniCS (the group I'm with already have a working version in MATLAB), and I ran into a problem I can't seem to find a way to implement.

So the weak formulation of the governing equations can be stated in cartesian and cylindrical coordinates (and probably others as well, but that is not really relevant), and the cartesian one is quite straightforward; one just ends up with a bunch of matrices that needs to be added and/or multiplied with the test and trial functions, so as far as I can see that would be easy in FEniCS as well.

The problem arises when I want to use the cylindrical coordinate system (it is more optimized for the problem because one dimension becomes irrelevant, and this is what the already working MATLAB code uses so comparisons are easier), because I need to define an operator that looks like this:

$$
[L_u] =
\begin{bmatrix}
\frac{\delta}{\delta r} & 0 \\
\frac{1}{r} & 0 \\
0 & \frac{\delta}{\delta z} \\
\frac{\delta}{\delta z} & \frac{\delta}{\delta r} \\
\end{bmatrix}
$$

How would I define something like this in FEniCS? I assume that what I need is something similar to how i.e. the grad function is defined but i can't seem to find anything like that in the examples (I might be looking in all the wrong places). I also don't really see how I would acquire the value for r so I can divide by it (see second line in the operator) since the variable r doesn't really pop up in my equations other than for this coordinate transform.

Thanks in advance, and apologies for the (probably) dumb question ^^,

EDIT:

I tried to implement it like this, but this only gives me an UFLException: Cannot create a tensor by joining subexpressions with different shapes.

def Lu(x):
    r = Expression("x[0]")

    ddr = x.dx(0)
    ddz = x.dx(1)

    return as_matrix([[ddr, 0],
                      [1/r, 0],
                      [0, ddz],
                      [ddz, ddr]])

chris_richardson · Answer 1 · Jun 18, 2016

You can define r as r=Expression("x[0]") i.e. the x-coordinate. You may find it easier to multiply your equation through by r to avoid singularities at r=0.

answered Jun 18, 2016 by chris_richardson FEniCS Expert (31,740 points)

Thank you, that partially answers my question, but how do I define the operator?

commented Jun 18, 2016 by hagen FEniCS Novice (120 points)

How do I implement a custom operator involving derivatives?

1 Answer