# Is there a way I can diagnose errors in my FEM formulation from principles?

I'm working through an FEM calculation by hand to verify I understand the algorithm before I write a bunch of code based on incorrect assumptions. However, when I iterate through the calculations, the values diverge. Can I diagnose my issue by inspecting the coefficient matrix or problem setup for concordance with some known properties?

## Setup

### Geometry

The shape is a trapezoid divided into three triangular elements, for which I have approximated the value (potential) within with a linear function: For the sake of an example problem, I assigned the following coordinates and values to the global indices:

| Node | Coordinate | Value (V_i) |
|------+------------+-------------|
|    1 | 1,1        | 100         |
|    2 | 0,0        | 100         |
|    3 | 2,1        | *           |
|    4 | 1,0        | *           |
|    5 | 3,0        | *           |


I computed the area for each triangle $i$ as

$$2A_i = (x^{(i)}_1 y^{(i)}_2 - x^{(i)}_2 y^{(i)}_1) + (x^{(i)}_3 y^{(i)}_1 - x^{(i)}_1 y^{(i)}_3) + (x^{(i)}_2 y^{(i)}_3 - x^{(i)}_3 y^{(i)}_2)$$

where the indices are the local index labels inside each triangle.

I notice that if I place the origin in the upper left, the counter-clockwise labeling no longer produces positive $A$ for triangles. However, as I understand it, that's only a convention and shouldn't affect the calculation.

| Triangle | 2A   |
|----------|------|
| 1        | -1   |
| 2        | -1   |
| 3        | -2   |


### Coefficient Matrix

The coefficient matrix for this problem is given by the symmetric form

\begin{equation} C = \begin{bmatrix} C_{11}^{(1)} + C_{11}^{(2)} && C_{13}^{(1)} && C_{12}^{(2)} && C_{12}^{(1)} + C_{13}^{(2)} && 0 \\ && C_{33}^{(1)} && 0 && C_{23}^{(1)} && 0 \\ && && C_{22}^{(2)} + C_{11}^{(3)} && C_{23}^{(2)} + C_{13}^{(3)} && C_{12}^{(3)} \\ && && && C_{22}^{(1)} + C_{33}^{(2)} + C_{33}^{(3)} && C_{32}^{(3)} \\ && && && && C_{22}^{(3)} \end{bmatrix} \end{equation}

where the local terms are

\begin{align} C_{11}^{(e)} &= \frac{1}{4A} \left[ (y_2 - y_3)^2 + (x_3 - x_2)^2 \right] \\ C_{12}^{(e)} &= \frac{1}{4A} \left[ (y_2 - y_3)(y_3 - y_1) + (x_3 - x_2)(x_1 - x_3) \right] \\ C_{13}^{(e)} &= \frac{1}{4A} \left[ (y_2 - y_3)(y_1 - y_2) + (x_3 - x_2)(x_2 - x_1) \right] \\ C_{21}^{(e)} &= C_{12}^{(e)} \\ C_{22}^{(e)} &= \frac{1}{4A} \left[ (y_3 - y_1)^2 + (x_1 - x_3)^2 \right] \\ C_{23}^{(e)} &= \frac{1}{4A} \left[ (y_3 - y_1)(y_1 - y_2) + (x_1 - x_3)(x_2 - x_1) \right] \\ C_{31}^{(e)} &= C_{13}^{(e)} \\ C_{32}^{(e)} &= C_{23}^{(e)} \\ C_{33}^{(e)} &= \frac{1}{4A} \left[ (y_1 - y_2)^2 + (x_2 - x_1)^2 \right] \end{align}

which I have computed to be (symmetric terms not shown):

\begin{equation} \begin{bmatrix} -0.75 & 0.25 & -0.25 & 0.5 & 0 \\ & -0.25 & 0 & 0.5 & 0 \\ && -0.5 & 0.25 & 0.25 \\ &&& -0.5 & 0.25 \\ &&&& -0.25 \end{bmatrix} \end{equation}

## Solution (Attempts)

According to the text, if I want to iteratively solve for $V_i$, I just have to repeat the following calculation until it converges, setting the unknown values to $0$ initially:

\begin{equation} V_k = - \frac{1}{C_{kk}} \sum^{n}_{i = 1, i \neq k} V_i C_{ki} \end{equation}

### Attempt 1

However, when I do this, the values diverge:

| Iteration |          V3 |          V4 |          V5 |
|-----------+-------------+-------------+-------------|
|         0 |           0 |           0 |           0 |
|         1 |         -50 |         200 |           0 |
|         2 |          50 |         125 |         150 |
|         3 |        87.5 |         300 |         175 |
|         4 |       187.5 |      331.25 |       387.5 |


The equations have some of the right properties (non-adjacent nodes don't directly contribute to the value of a node in the iterative equations, the matrix is symmetric), but even going back and correcting for arithmetic errors, the divergent behavior remains the same.

### Attempt 2

On the other hand, starting with something like

| Node | Coordinate | Value (V_i) |
|------+------------+-------------|
|    1 | 1,1        | 100         |
|    2 | 0,0        | *           |
|    3 | 2,1        | 100         |
|    4 | 1,0        | 0           |
|    5 | 3,0        | *           |


results in immediate convergence.

Is there something I can look for instead of repeating calculations?

• It would be a good idea to test a problem where you have an analytic solution already (e.g. a Laplace equation on a rectangle with some simple boundary conditions). – Ian Jan 30 '17 at 4:12
• I'm not sure I understand the context, but if you are trying to solve the linear system corresponding to the matrix (partially?) shown above "Solution (attempts)", a Jacobi iteration is not guaranteed to converge because the matrix is not diagonally dominant. Thus the problem seems to have little to do with the finite element method but rather with your approach to solving linear systems. – hardmath Jan 30 '17 at 4:14
• @hardmath That's disappointing, especially considering that the author didn't go into the convergence behavior of the iterative algorithm when presenting it. Regarding the context, is there anything I should add to the question to make it more understandable? – bright-star Jan 30 '17 at 4:16
• It is possible that the reason the matrix is not diagonally dominant originates with the matrix assembly steps preceding your attempt to solve the linear system. – hardmath Jan 30 '17 at 4:17
• Yes you won't get Jacobi to work here, but maybe Conjugate-Gradient will. – mathreadler Jan 30 '17 at 5:12

@hardmath pointed out that Jacobi iteration (indicated by the text's author) requires [weak] diagonal dominance. I checked my arithmetic for the coupling matrices, and found that the coefficient matrix was in fact diagonally dominant. After computing with the new matrix and same initial conditions, the solution converged.

| Triangle | 2A   |
|----------|------|
| 1        | -2   |
| 1        | -2   |
| 1        | -4   |


\begin{equation} C = \begin{bmatrix} -1.5 & 0 & 0.5 & 1 & 0 \\ & -0.5 & 0 & 0.5 & 0 \\ && -1 & 0 & 0.5 \\ &&& -1.75 & -0.5 \\ &&&& -0.5 \end{bmatrix} \end{equation}

| Iteration |     V3 |    V4 |     V5 |
|-----------+--------+-------+--------|
|         0 |      0 |     0 |      0 |
|         1 |     50 |   150 |      0 |
|         2 |     50 |   150 |   -100 |
|         3 |      0 |   200 |    -50 |
|         4 |     25 |   175 |   -100 |
|         5 |    -50 |   225 |    -75 |
|         6 |   12.5 | 112.5 | -137.5 |
|         7 | -18.75 | 81.25 |    -50 |
|         8 |     25 |   125 |    -50 |
|         9 |     25 |   125 |    -50 |