We have written a two dimensional, adaptive, finite element MHD code,
called FEMHD. The code solves the equations
(7), (9), and
(8),
(10)
or (12),
(13).
The code
is written in modular form.
One module generates the mesh and various arrays which specify
its connectivity. The fundamental mesh arrays are the lists
of N vertices 
and an array which lists the
three vertex indices of each of the
triangles.
From these arrays, useful secondary arrays are constructed. These arrays list the neighboring vertices of each vertex, which are connected by a triangle side. Also listed are the number of neighbors, as well as the triangle numbers of the triangles sharing a particular vertex.
Another module contains routines for constructing and and assembling the mass and stiffness matrices, the Poisson bracket tensor, and for inverting sparse symmetric matrices. We invert the stiffness and mass matrices using a sparse matrix Conjugate Gradient algorithm with an Incomplete Cholesky preconditioner.
The discretized equations are solved in a separate module. The time derivatives are represented by second order finite differences. The equations are advanced explicitly, using a leapfrog time differencing, but with the dissipative terms in the equations (7), (8) solved implicitly, by combining them with the mass matrix.