# MASSMATRIX How Does FreeFem++ Define a Mass Matrix?

MASSMATRIX, a FreeFem++ script which investigates the computation of the mass matrix.

Let Omega be a domain, and let X be a set of nodes associated with a mesh of Omega. Let Phi(I)(X) be the I-th basis function, associated with node I. Then the finite element method defines the mass matrix M as

```        M(i,j) = Integral ( x in Omega ) Phi(I)(X) Phi(J)(X) dx
```
We note that it must be the case that M is a symmetric positive definite matrix, in part because, for any vector C, the expression { C' * M * C ) can be rewritten as
```        C' M C = Sum ( I, J ) ( Integral ( x in Omega ) C(I) * Phi(I)(X) Phi(J)(X) * C(J) dx )
= || C*Phi ||^2 >= 0
```
For a Lagrangian basis Phi, if C is a vector of 1's, then the corresponding finite element function is 1 everywhere, so
```        C' M C = 1' * M * 1
= Sum ( I, J ) Integral ( x in Omega ) Phi(I)(X) Phi(J)(X)  dx
= Sum ( J ) Integral ( x in Omega ) Sum ( I ) Phi(I)(x) Phi(J)(X) dx
= Sum ( J ) Integral ( x in Omega ) 1 * Phi(J)(X) dx
= Integral ( x in Omega ) 1 dx
= Area ( Omega ).
```

The following FreeFem++ commands create the mass matrix for a given problem:

```        mesh Th = square ( 3, 3 );
fespace Vh ( Th, P1 );
Vh uh;
Vh vh;
varf bid ( uh, vh ) = int2d ( Th ) ( uh * vh ) + on ( 1, 2, 3, 4, uh = 0.0 );
matrix ad = bid ( Vh, Vh );
```

Imagine our surprise to see the matrix begin with:

```         1         2         3         4         5         6         7         8
1: 1.00E+30  0.46E-02  0.00E+00  0.00E+00  0.46E-02  0.92E-02  0.00E+00  0.00E+00 ...
2: 0.46E-02  1.00E+30  0.46E-02  0.00E+00  0.00E+00  0.92E-02  0.92E-02  0.00E+00 ...
3: 0.00E+00  0.46E-02  1.00E+30  0.46E-02  0.00E+00  0.00E+00  0.92E-02  0.92E-02 ...
```
The entries 1.00E+30 are total surprises. However, this surprise becomes explicable when we realize that the strange entries only occur on diagonal entries, and only on diagonal entries associated with boundary nodes. For a problem with Dirichlet conditions, the basis functions that would ordinarily be generated at the corresponding node are not needed, at least not in the generation of the finite element equations. But somehow, in this example, these basis functions are being generated, because the offdiagonal elements are correct.

To see this, generate the same problem, but with all Neumann conditions:

```        mesh Th = square ( 3, 3 );
fespace Vh ( Th, P1 );
Vh uh;
Vh vh;
varf bin ( uh, vh ) = int2d ( Th ) ( uh * vh );
matrix an = bin ( Vh, Vh );
```

The only change occurs to the diagonal elements, which are suddenly "cured":

```         1         2         3         4         5         6         7         8
1: 0.18E-01  0.46E-02  0.00E+00  0.00E+00  0.46E-02  0.92E-02  0.00E+00  0.00E+00 ...
2: 0.46E-02  0.27E-02  0.46E-02  0.00E+00  0.00E+00  0.92E-02  0.92E-02  0.00E+00 ...
3: 0.00E+00  0.46E-02  0.27E-02  0.46E-02  0.00E+00  0.00E+00  0.92E-02  0.92E-02 ...
```

To me, this suggests only that the internal mechanism for computing the mass matrix or stiffness matrix must have a patch that says that if a node is associated with a boundary condition, then diagonal elements of the matrix should be replace by 1.0E+30. Presumably, when a nonzero Dirichlet value V is used, the corresponding right hand side becomes 1.0E+30 * V.

This seems to me an awkward ad hoc and unmathematical way to deal with boundary conditions, and especially unfortunate since I have not, as yet, found any documentation that reports that this is going on, let alone explains why it is done this way. However, we had a major loss of faith in our coding abilities when we first encountered this anomaly while trying to compute a mass matrix, and we thought it would be a good thing to write it all down as best we could.