stiffness_matrix, a FreeFem++ code which investigates the computation of the stiffness 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 stiffness matrix K as
K(i,j) = Integral ( x in Omega ) Phi'(I)(X) Phi'(J)(X) dx
We note that it must be the case that K is a symmetric positive definite
matrix, in part because, for any vector C, the expression { C' * K * C )
can be rewritten as
C' K 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, and its derivative is 0 everywhere, so
C' K C = 1' * K * 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 ) 0 * Phi(J)(X) dx
= 0.
The following FreeFem++ commands create the stiffness matrix for a given problem:
mesh Th = square ( 3, 3 );
fespace Vh ( Th, P1 );
Vh uh;
Vh vh;
varf bid ( uh, vh ) = int2d ( Th ) ( dx(uh) * dx(vh) + dy(uh) * dy(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.0E+30 -0.5 0.0 0.0 -0.5 0.0 0.0 0.0 ...
2: -0.5 1.0E+30 -0.5 0.0 0.0 -1.0 0.0 0.0 ...
3: 0.0 -0.5 1.0E+30 -0.5 0.0 0.0 -1.0 0.0 ...
4: 0.0 0.0 -0.5 1.0E+30 0.0 0.0 0.0 -0.5 ...
5: -0.5 0.0 0.0 0.0 1.0E+30 -1.0 0.0 0.0 ...
6: 0.0 -1.0 0.0 0.0 -1.0 4.0 -1.0 0.0 ...
The entries 1.0E+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 ) ( dx(uh) * dx(vh) + dy(uh) * dy(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: 1.0 -0.5 0.0 0.0 -0.5 0.0 0.0 0.0 ...
2: -0.5 2.0 -0.5 0.0 0.0 -1.0 0.0 0.0 ...
3: 0.0 -0.5 2.0 -0.5 0.0 0.0 -1.0 0.0 ...
4: 0.0 0.0 -0.5 1.0 0.0 0.0 0.0 -0.5 ...
5: -0.5 0.0 0.0 0.0 2.0 -1.0 0.0 0.0 ...
6: 0.0 -1.0 0.0 0.0 -1.0 4.0 -1.0 0.0 ...
To me, this suggests only that the internal mechanism for computing the 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 stiffness matrix, and we thought it would be a good thing to write it all down as best we could.
The computer code and data files described and made available on this web page are distributed under the MIT license
