Let 
denote a triangulation of 
and let
be the space of 
piecewise linear polynomials
associated with .
PLTMG usually represents such a piecewise polynomial using
the standard nodal basis;
a function can then be specified by giving its values at
the vertices.
Let 
be the subspace of whose elements are
zero at the knots of lying on 
(the portion of the boundary with Dirichlet boundary conditions),
and which satisfy the continuity condtions on 
.
Similarly,
let 
be the affine space of whose elements
satisfy the Dirichlet boundary conditions at knots of
lying on and the continuity
conditions on .
The discrete equations solved by PLTMG
are: find 
, and 
such that
where
If the Jacobian is not self-adjoint, some upwinding terms based on the Scharfetter-Gummel discretization scheme [7] are added.
When the continuation process is used, we use a normalization equation of the form:
The scalar 
is the steplength.
PLTMG uses two different
normalization equations [10, 29].
Most frequently, we use:
Here 
is a parameter selected by PLTMG; by choosing
and properly, it is possible to achieve
target values in either 
or 
.
The vector 
is the current solution point
and 
the current unit tangent
vector. The
scalar 
is defined formally using the chain rule
for differentiation:
The other normalization used in PLTMG is based on the pseudo-arclength method, characterized by:
The normalization (
) is used when
() is not well defined (e.g.,
if 
). For reasonable
choices of , these will be isolated points on the
solution manifold, such as symmetry breaking bifurcation points.
In such instances,
() is used on a temporary basis until ()
is defined again.
There are only seven relatively short subroutines
(ELEASM, ELENBC, ELEDBC, ELEUN, ELEUNP, ELEMTX, and ELEBDI)
that have knowledge of the discretization
procedure and the form of the differential equation
()-() and
the functional in ().
ELEASM is used for assembling stiffness matrices,
right hand sides, and the contribution to from P1XY
for a single element. ELEMTX assembles inexpensive element stiffness matrices
used in computing the ordering of the vertices for the hierarchical
basis iteration.
ELENBC computes contributions to
the natural boundary conditions from a single
element edge.
ELEDBC evaluates Dirichlet boundary conditions
and initial guesses along a single element edge.
ELEUN computes an integral on a single interior
element edge that is required for the a posteriori error
estimates. ELEUNP does a similar calculation for linked edges.
ELEBDI computes the contribution to from P2XY
for a single edge.
These routines are written to accept general quadrature rules. The details of the rule (weights and abscissas) are read from tables (that are initialized in DATA statements) making it easy to change the formulas. Currently ELEASM uses a 3-point rule and ELEMTX a 1-point rule. ELENBC, ELEUN, ELEUNP and ELEBDI use 2-point-Gauss rules. The organization within these routines favors clarity at the expense of efficiency, making it reasonably simple to change some aspects of the problem class and discretization.
Global stiffness matrices are stored in the sparse matrix format described in [18]. All nonzeros are stored in a linear array; first the diagonal entries are stored, followed by the upper triangular entries, stored column by column. If the matrix is nonsymmetric, this is followed by the lower triangular entries, stored row by row. Symmetric and nonsymmetric storage is governed by the parameter ISPD; ISPD=1 designates the symmetric storage option, while ISPD=0 indicates nonsymmetric storage.
