Partial Differential Equation Toolbox™ software has a function for global, uniform mesh refinement for 2-D geometry. It divides each triangle into four similar triangles by creating new corners at the midsides, adjusting for curved boundaries. You can assess the accuracy of the numerical solution by comparing results from a sequence of successively refined meshes. If the solution is smooth enough, more accurate results may be obtained by extrapolation.

The solutions of equations often have geometric features like
localized strong gradients. An example of engineering importance in
elasticity is the stress concentration occurring at reentrant corners
such as the MATLAB^{®} L-shaped membrane. Then it is more economical
to refine the mesh selectively, i.e., only where it is needed. When
the selection is based on estimates of errors in the computed solutions,
a posteriori estimates, we speak of *adaptive mesh refinement*.
See `adaptmesh`

for an example of
the computational savings where global refinement needs more than
6000 elements to compete with an adaptively refined mesh of 500 elements.

The adaptive refinement generates a sequence of solutions on
successively finer meshes, at each stage selecting and refining those
elements that are judged to contribute most to the error. The process
is terminated when the maximum number of elements is exceeded, when
each triangle contributes less than a preset tolerance, or when an
iteration limit is reached. You can provide an initial mesh, or let `adaptmesh`

call `initmesh`

automatically. You also choose selection
and termination criteria parameters. The three components of the algorithm
are the error indicator function, which computes an estimate of the
element error contribution, the mesh refiner, which selects and subdivides
elements, and the termination criteria.

The adaptation is a feedback process. As such, it is easily
applied to a larger range of problems than those for which its design
was tailored. You want estimates, selection criteria, etc., to be
optimal in the sense of giving the most accurate solution at fixed
cost or lowest computational effort for a given accuracy. Such results
have been proved only for model problems, but generally, the equidistribution
heuristic has been found near optimal. Element sizes should be chosen
such that each element contributes the same to the error. The theory
of adaptive schemes makes use of *a priori* bounds
for solutions in terms of the source function *f*.
For nonelliptic problems such a bound may not exist, while the refinement
scheme is still well defined and has been found to work well.

The error indicator function used in the software is an elementwise
estimate of the contribution, based on the work of C. Johnson et al. [5], [6]. For Poisson's
equation –Δ*u* = *f* on
Ω, the following error estimate for the FEM-solution *u _{h}* holds
in the

$$\Vert \nabla (u-{u}_{h})\Vert \le \alpha \Vert hf\Vert +\beta {D}_{h}({u}_{h}),$$

where *h* = *h*(*x*) is
the local mesh size, and

$${D}_{h}(v)={\left({\displaystyle \sum _{\tau \in {E}_{1}}{h}_{\tau}^{2}{\left[\frac{\partial v}{\partial {n}_{\tau}}\right]}^{2}}\right)}^{1/2}.$$

The braced quantity is the jump in normal derivative of *v* across
edge *τ*, *h*_{τ} is
the length of edge *τ*, and the sum runs over *E _{i}*,
the set of all interior edges of the triangulation. The coefficients

The general form of the error indicator function for the elliptic equation

–∇ · (*c*∇*u*)
+ *au* = *f*

is

$$E\left(K\right)=\alpha {\Vert h\left(f-au\right)\Vert}_{K}+\beta {\left(\frac{1}{2}{\displaystyle \sum _{\tau \in \partial K}{h}_{\tau}^{2}{\left({n}_{\tau}\xb7\text{\hspace{0.17em}}c\nabla {u}_{h}\right)}^{2}}\right)}^{1/2},$$

where $${n}_{\tau}$$ is the unit normal of edge *τ* and
the braced term is the jump in flux across the element edge. The *L*_{2} norm
is computed over the element *K*. This error indicator
is computed by the `pdejmps`

function.

Partial Differential Equation Toolbox software is geared to elliptic problems. For reasons of accuracy and ill-conditioning, they require the elements not to deviate too much from being equilateral. Thus, even at essentially one-dimensional solution features, such as boundary layers, the refinement technique must guarantee reasonably shaped triangles.

When an element is refined, new nodes appear on its midsides,
and if the neighbor triangle is not refined in a similar way, it is
said to have *hanging nodes*. The final triangulation
must have no hanging nodes, and they are removed by splitting neighbor
triangles. To avoid further deterioration of triangle quality in successive
generations, the "longest edge bisection" scheme Rosenberg-Stenger [8] is used, in which the
longest side of a triangle is always split, whenever any of the sides
have hanging nodes. This guarantees that no angle is ever smaller
than half the smallest angle of the original triangulation.

Two selection criteria can be used. One, `pdeadworst`

,
refines all elements with value of the error indicator larger than
half the worst of any element. The other, `pdeadgsc`

,
refines all elements with an indicator value exceeding a user-defined
dimensionless tolerance. The comparison with the tolerance is properly
scaled with respect to domain and solution size, etc.

For smooth solutions, error equidistribution can be achieved
by the `pdeadgsc`

selection if the maximum number
of elements is large enough. The `pdeadworst`

adaptation
only terminates when the maximum number of elements has been exceeded
or when the iteration limit is reached. This mode is natural when
the solution exhibits singularities. The error indicator of the elements
next to the singularity may never vanish, regardless of element size,
and equidistribution is too much to hope for.

Was this topic helpful?