Now that we have a feeling for how the finite element method works on a theoretical level, let's now work up the practical machinery for handling 2D problems. We note that the extension to 3D is straightforward in terms of the theory and algorithms. The difficulties when extending to 3D (the so-called curse of dimensionality) have to do with the size of the resulting system and in creating the finite element mesh.
Let's start with an arbitrary 2D planar domain 
and break it up
into discrete elements to form a finite dimensional subspace, 
.
For 2D planar domains we have the choice of representing our function (if
we assume for now that we will use linear elements) as either triangles,
or as quadralaterals,
For now, let's restrict our development to triangles, knowing that it is easy to modify our formulae for quadralaterals. We take out a specific triangle from our finite dimensional subspace and apply the previous formulations for the three vertices,
or in matrix form,
We can now solve for the coefficients:
where
If we now define,
we can express our coefficients as
with
were we can permute the indicies 
to find 
Now we can write an expression for a general point in our triangle in terms
of the values at the vertices,
or, most succintly,
with 
, the solution value at node 
, and 
is called
the local shape function or basis function. This can be
expressed in a variety of ways in the literature (depending, usually, on
whether you are reading engineering or mathematical treatments of finite
element analysis):
From our definition of 
we can see that 
at node 
and
at nodes 
. In general this can be expressed as,
where 
is the Kroneker delta.
Our expression for the shape functions is, in fact, a transformation from
cartesian coordinates to a natural set of coordinates, often called area coordinates or barycentric coordinates. If we look at the
transformation for 
for a point within the element with coordinates
, we note that,
which is just a ratio of the areas of the triangles 
and 
.
Thus we have a transformation from a point inside an element in terms of
its barycentric coordinates to the global cartesian system,
Now that we have a suitable set of basis functions, we can find the finite element approximation to our 2D problem. Recall that our orginal problem can be formulated as,
where
and
The finite element approximation to the original boundary value problem is
which has the equivalent form
where
which can be expressed by the matrix and vector elements,
and
Fortunately, the above quantities are easy to evaluate for linear
triangles. As a matter of fact, there are closed form solutions for the
matrix elements 
:
Therefore,
and for the right hand side we have, assuming constant sources,
Thus the local, element matrix looks like
which have the compact forms,
and
Now we must add up all the contributions from each element into a global matrix and global vector.
where 
is equal to the total number of elements in the discretized
solution domain and represents the node numbers (vertices). Let's say
we want to add the contributions from two elements which have global
numbers 
for element one and 
for element two. The local element
matrices look like
and
where the 
denotes symmetry. The compenents of represent the
connectivities in the global matrix and are the contributions from the
) pair. Thus the two element system would have the
following global matrix in terms of the local contributions:
This process continues until we have added in all our element contributions to form the global system. Now we have only to take into account of our boundary conditions and then solve our system.
