We divide the domain where the problem is defined into a series of small,
non-overlapping elements 
(quadrilaterals in the 2D case, and hexahedral in 3D, the elements
may have linear
or arbitrary curvilinear geometry). Using appropriately defined
degrees of freedom U
(unknowns of the discretized problem), and corresponding shape
functions N, we obtain
an approximate, finite dimensional formulation of the problem:
where 
indicates a finite approximation of the continuum function on the
finite element grid of elements.
The definition of the shape functions and degrees of freedom are given below.
An evaluation of integrals in the definition of B(.,.) and L(.) is
performed element
by element and resulting element matrices (element stiffness matrix 
,
and element load vector 
) are assembled to form a global set of
algebraic equations
The discretization of the domain (grid, i.e. the definition of the
elements, and appropriate shape functions) must result in a continuous
approximation of the solution, which can usually be satisfied if and only if
the two neighboring
elements, sharing a common edge or face, share this common boundary
fully (``1-to-1'' rule). Otherwise, discontinuities (gaps) in the solution
may appear. Computation of local matrices and is always performed
by numerical integration. An explicit evaluation of the integrals is
practically
impossible, even for linearly geometric elements, since for 
and
- versions
of the method, the number of different element ``types'' may exceed any
reasonable value (Limiting the maximum order of approximation
to 8, we may
encounter 32768 different elements in an isotropic 2D problem, and 
in the anisotropic 3D case). We will discuss some problems related to building
element stiffness matrices and right hand side vectors, and to global
element matrices assembly in later sections, since, due to the use of
irregular
meshes, it is not a straightforward adaptation of standard
algorithms.
After solving the resulting set of simultaneous algebraic equations (whether
it is linear or not, is of no interest here), one has to compute various error
estimators/indicators, which will guide the adaptation process, pointing to
elements with local unsatisfactory discretization, so the mesh may be locally
improved and better solution obtained. At some moment in this process, one
either obtains the solution with satisfactory value of global error, or one
exceeds the allowable resources (CPU time, memory, or disk size). It has been
shown that, due to so called ``exponential convergence'', --adaptive
strategy
can produce the best possible solution within given resource limits.