# Morley element implementation reference

I am looking for a detailed reference on the implementation of the Morley element for FEM, specifically for the biharmonic equation. By detailed, I mean that it should discuss the problems associated with transforming the reference element, as well as handling different boundary conditions.

• Have you read any reference already? – nicoguaro Mar 4 at 17:53
• Why would you want to use the Morley element (or any other plate element, for that matter)? There are much better techniques today for the biharmonic equation... – Wolfgang Bangerth Mar 5 at 3:59
• Take a look here: dealii.org/developer/doxygen/deal.II/step_47.html – Wolfgang Bangerth Mar 5 at 20:34
• We can't integrate exactly in all but the very simplest cases. And in all of these simplest cases, the quadrature yields the exact integral :-) – Wolfgang Bangerth Mar 6 at 14:10
• I prefer implementing Morley elements directly in the global coordinate system without using a reference element. I simply specify any polynomial basis for $P_2$. For example, $1, x, y, xy, x^2, y^2$ and then perform change of basis in each global element to find the global basis functions. – knl Mar 6 at 14:45

I can write my experiences here because I do not have any book references at hand.

Consider a triangular element with the corner points $$\boldsymbol{x}_i \in \mathbb{R}^2$$, $$i=1,2,3$$. The degrees of freedom for the Morley element are $$F_i(v)=v(\boldsymbol{x}_i), \quad i=1,2,3.$$ and $$F_4(v)=\frac{\partial v}{\partial \boldsymbol{n}}\left(\frac12(\boldsymbol{x}_1 + \boldsymbol{x}_2)\right),\quad F_5(v)=\frac{\partial v}{\partial \boldsymbol{n}}\left(\frac12(\boldsymbol{x}_1 + \boldsymbol{x}_3)\right),\quad F_6(v)=\frac{\partial v}{\partial \boldsymbol{n}}\left(\frac12(\boldsymbol{x}_2 + \boldsymbol{x}_3)\right),$$ where $$\boldsymbol{n}$$ is the outward unit normal.

Now let $$p_i$$, $$i=1,\dots,6$$, denote the following power basis for $$P_2$$: $$p_1=1, ~ p_2=x, ~ p_3=y, ~ p_4=x^2, ~ p_5=xy, ~ p_6 = y^2.$$ For each element you can find the global basis functions corresponding to the above degrees of freedom by inverting the matrix $$V_{ij} = F_i(p_i)$$, $$i,j=1,\dots,6$$. The $$i$$'th row of $$V^{-1}$$ corresponds to the coefficients in the following representation of the global finite element basis: $$\varphi_i(\boldsymbol{x}) = \sum_{j=1}^6 V^{-1}_{ij} p_j(\boldsymbol{x}).$$ This follows from the fact that the degrees of freedom and the basis functions should satisfy $$F_i(\varphi_j) = \delta_{ij}$$.

You can now calculate $$V^{-1}$$ in each element and then evaluate the resulting basis functions at the global quadrature points, or possibly first differentiate the resulting basis functions and then evaluate. In general $$V$$ can have a fairly bad condition number but I have never encountered such issues in pratical computations expect when implementing the Argyris basis and computing with an overrefined mesh.

You can also take a look at the generic Python implementation that I have successfully used in the past for Morley elements.

• Thank you, I'll try to implement it in the following days. I also finally managed to find a good paper on transforming elements (Morley included): arxiv.org/abs/1706.09017 – lightxbulb Mar 6 at 20:20
• @lightxbulb, would you mind adding an answer related to the main points in the linked paper? – nicoguaro Mar 10 at 19:35
• @nicoguaro I am working on that. I want to get a working implementation on several methods, and then I'll try to summarize my findings and I'll accept knl's answer. The biggest advantage from Kirby's paper seems to be that a 6x6 matrix inversion is not required, and one can obviously work directly on the reference element (which also simplifies analytical integration somewhat). – lightxbulb Mar 10 at 20:58
• @knl I tried what you have described, but for my meshing I get non-invertible matrices even with QR decomposition. To avoid this I implemented the transformation described in Kirby's paper, but my conjugate gradient solver still doesn't produce a reasonable result. Is there something specific with how one goes about setting Dirichlet boundary conditions for the Morley element compared to Lagrange elements? In the Lagrange case I was subtracting the relevant column multiplied by the Dirichlet value from the load vector in order to get a symmetric matrix. – lightxbulb May 14 at 10:57

This text contains a brief discussion on the Morley element:

http://www.csc.kth.se/~jjan/transfer/fenics-superparametric/Hansbo.pdf