The notation employed in this section is very powerful, yet very
simple. Rather than use a matrix notation, the convention is used
whereby ANY repeated subscript is assumed to imply a summation over its
appropriate range (i and j over coordinates x and y (and z if
a 3-D formulation is desired), and k and l over 1 to 6 for a 2-D
3-node triangular element (1 to 8 for a 2-D 4-node quadrilateral element
or 1 to 24 for a 3-D 8-node brick element). When a comma appears in a
subscript list, this implies partial differentiation with respect to
the appropriate coordinate. Thus 
has an implied summation
over the coordinates, and represents 
(
if 3-D).
Conservation of linear momentum:
Constitutive equation or flow law:
where
and
and
The differential equation is now:
Doing the variational and substituting for 
:
or
We have the vector identity:
which we can use to expand the first term of the integral equation above, as well as the fact that the first term can be simplified by the following
and we obtain
The second integral in the right hand side is an area integral of a divergence, and so can be converted to a line integral of the normal over the appropriate portion of the boundary. The above equation becomes:
where
is the normal stress on the surface.
The above equations apply for the domain as a whole, but they also
apply for a particular subdomain, or element. Thus we would have a set
of element equations corresponding to the above, where only 
and 
would be changed to 
and 
for element e. We now express the vectors 
and 
as sums of products of shape functions and values of 
and 
or
and 
at each of the k nodes in an element. We have:
and
where
and
and
and
and
and
for 3-node triangular elements. The 
's are bilinear shape
functions defined to be 1 at precisely one node, and 0 at all the
rest (
).
With the above substitution for 
and 
, derivative terms such
as 
and 
become
and
because only the 
's depend on x and y, since the 
's are
values defined only at the nodes. Our element equation becomes
Since 
is arbitrary we can generate the k+1 equations
representable as the matrix equation:
where m and n now run from 1 to k+1.
The matrix K is
where the k+1 line of this matrix stems from the incompressibility condition for ice:
Integrating 
over the area;
and then substituting for 
or
The vector 
is
The vector 
is
The vector 
is
The 
part of this looks like this:
where
