# 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? Mar 4, 2020 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... Mar 5, 2020 at 3:59
• Take a look here: dealii.org/developer/doxygen/deal.II/step_47.html Mar 5, 2020 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 :-) Mar 6, 2020 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, 2020 at 14:45

## 3 Answers

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 Mar 6, 2020 at 20:20
• @lightxbulb, would you mind adding an answer related to the main points in the linked paper? Mar 10, 2020 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). Mar 10, 2020 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. May 14, 2020 at 10:57
• @knl I finished with the implementation of the Morley element using Kirby's paper. See my answer to my own question below if you're interested. I believe that my derivation is correct, even though I still get nonsensical results from my current implementation. Jun 3, 2020 at 23:50

This text contains a brief discussion on the Morley element:

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

• Can you elaborate on the points touched in the paper? Jun 4, 2020 at 19:17
• "The requested URL was not found on this server." Apr 29 at 9:29

I have implemented Morley elements for the biharmonic equation according to Kirby's paper:

arxiv.org/abs/1706.09017

I was specifically interested in the 2D problem:

$$\Delta^2 u(x) = 0, x \in \Omega$$ $$u(x) = g(x), \frac{\partial u}{\partial n}(x) = \frac{\partial g}{\partial n}(x), x \in \partial \Omega$$

Disclaimer: I am not certain my implementation is working entirely correctly, though the theory should be fine.

In Kirby's paper, a transformation matrix is derived $$V = M^T$$ (equation $$(39)$$), that allows one to transform the reference nodal basis functions $$\hat{\Psi}$$ to the physical element's nodal basis functions $$\Psi$$, through: $$\Psi = M(\hat{\Psi} \circ F)$$. The matrix $$M$$ was computable in cases where the 6x6 matrix inversion for the simpler method in @knl's answer was not. For the above problem, only $$\int \Delta\Psi_i\Delta\Psi_j$$ is required to be computed per element.

For a canonical reference element I used a triangle with vertices $$v_0 = (0,0), v_1 = (1,0), v_2 = (0,1)$$, edge midpoints: $$m_0 = (0.5,0.5), m_1 = (0, 0.5), m_2 = (0.5, 0)$$ and corresponding normals at those: $$n_0 = (\frac{1}{\sqrt{2}}, \frac{1}{\sqrt{2}}), n_1 = (1,0), n_2 = (0,-1)$$. I write every basis function as: $$\hat{\Psi}_i = a_0 + a_1x + a_2y + a_3x^2 + a_4xy+a_5y^2$$. To compute the nodal basis functions I have used $$\hat{N}(\hat{\Psi}) = I$$, where: $$\hat{N} = (\delta_{v_0}, \delta_{v_1}, \delta_{v_2}, \delta^{n_0}_{m_0}, \delta^{n_1}_{m_1}, \delta^{n_2}_{m_2})$$,

such that $$\delta_{v}(f) = f(v)$$ and $$\delta^{n}_{v}(f) = n^T\nabla f(v)$$.

The system: $$\hat{N}(\hat{\Psi}) = I$$ with the above reference nodes yields functions:

$$\hat{\Psi}_0 = 1 -x -y + 0x^2 + 2xy + 0y^2$$ $$\hat{\Psi}_1 = 0 +\frac{x}{2} +\frac{y}{2} + \frac{x^2}{2} - xy - \frac{y^2}{2}$$ $$\hat{\Psi}_2 = 0 +\frac{x}{2} +\frac{y}{2} - \frac{x^2}{2} - xy + \frac{y^2}{2}$$ $$\hat{\Psi}_3 = 0 -\frac{x}{\sqrt{2}} -\frac{y}{\sqrt{2}} + \frac{x^2}{\sqrt{2}} + \frac{2xy}{\sqrt{2}} + \frac{y^2}{\sqrt{2}}$$ $$\hat{\Psi}_4 = 0 +x +0y -x^2 + 0xy + 0y^2$$ $$\hat{\Psi}_5 = 0 +0x -y + 0x^2 + 0xy + y^2$$

Respectively with the following partial derivatives:

$$\partial_{\hat{x}\hat{x}}\hat{\Psi}_0 = 0,\, \partial_{\hat{x}\hat{y}}\hat{\Psi}_0 = 2,\, \partial_{\hat{y}\hat{y}}\hat{\Psi}_0 = 0$$ $$\partial_{\hat{x}\hat{x}}\hat{\Psi}_1 = 1,\, \partial_{\hat{x}\hat{y}}\hat{\Psi}_1 = -1,\, \partial_{\hat{y}\hat{y}}\hat{\Psi}_1 = -1$$ $$\partial_{\hat{x}\hat{x}}\hat{\Psi}_2 = -1,\, \partial_{\hat{x}\hat{y}}\hat{\Psi}_2 = -1,\, \partial_{\hat{y}\hat{y}}\hat{\Psi}_2 = 1$$ $$\partial_{\hat{x}\hat{x}}\hat{\Psi}_3 = \frac{2}{\sqrt{2}},\, \partial_{\hat{x}\hat{y}}\hat{\Psi}_3 = \frac{2}{\sqrt{2}},\, \partial_{\hat{y}\hat{y}}\hat{\Psi}_3 = \frac{2}{\sqrt{2}}$$ $$\partial_{\hat{x}\hat{x}}\hat{\Psi}_4 = -2,\, \partial_{\hat{x}\hat{y}}\hat{\Psi}_4 = 0,\, \partial_{\hat{y}\hat{y}}\hat{\Psi}_4 = 0$$ $$\partial_{\hat{x}\hat{x}}\hat{\Psi}_5 = 0,\, \partial_{\hat{x}\hat{y}}\hat{\Psi}_5 = 0,\, \partial_{\hat{y}\hat{y}}\hat{\Psi}_5 = 2$$

Now using $$\Psi = MF(\hat{\Psi} \circ F)$$ (where $$F$$ is the transformation from the physical element to the reference element) I get:

$$\Delta \Psi_i = \sum_{j=0}^{5}M_{i,j}\Bigg(\left(\left(\frac{\partial \hat{x}}{\partial x}\right)^2 + \left(\frac{\partial \hat{x}}{\partial y}\right)^2\right)\partial_{\hat{x}\hat{x}}\hat{\Psi}_j + \\ 2\left(\frac{\partial \hat{x}}{\partial x}\frac{\partial \hat{y}}{\partial x} + \frac{\partial \hat{x}}{\partial y}\frac{\partial \hat{y}}{\partial y}\right)\partial_{\hat{x}\hat{y}}\hat{\Psi}_j + \\ \left(\left(\frac{\partial \hat{y}}{\partial x}\right)^2 + \left(\frac{\partial \hat{y}}{\partial y}\right)^2\right)\partial_{\hat{y}\hat{y}}\hat{\Psi}_j\Bigg),$$

where $$(\hat{x}, \hat{y}) = F(x, y)$$. All the terms in the above expression are constants so $$\Delta \Psi_i$$ can be evaluated for each $$i$$, and then the element matrix $$K_{i,j} = \int \Delta \Psi_i \Delta \Psi_j = \frac{|det(J_{F^{-1}})|}{2} \Delta \Psi_i \Delta \Psi_j$$ can be computed. The form of $$M$$ is explicitly given in the cited paper.

An important implementation detail is matching the reference element correctly to the physical element wrt the normals. In Kirby's paper this is also discussed, as there the tangents are taken to have their tail at a vertex with a lower global index, and have their head at a vertex with a higher global index. The corresponding normals are: $$n_x = t_y, n_y = -t_x$$, where $$(t_x,t_y)$$ is a tangent. This implies that the triangle vertices are always sorted in ascending order wrt their global indices, which yields a consistent orientation of the normals.