Accelerating the pace of engineering and science

# Documentation

### Improving Solution Accuracy Using Mesh Refinement

Partial Differential Equation Toolbox™ software has a function for global, uniform mesh refinement. 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.

### Error Estimate for the FEM Solution

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 uh holds in the L2-norm $‖\cdot ‖$:

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

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

${D}_{h}\left(v\right)={\left(\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 Ei, the set of all interior edges of the triangulation. The coefficients α and β are independent of the triangulation. This bound is turned into an elementwise error indicator function E(K) for element K by summing the contributions from its edges.

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

–∇ · (cu) + au = f

is

$E\left(K\right)=\alpha {‖h\left(f-au\right)‖}_{K}+\beta {\left(\frac{1}{2}\sum _{\tau \in \partial K}{h}_{\tau }^{2}{\left({n}_{\tau }·\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 L2 norm is computed over the element K. This error indicator is computed by the pdejmps function.

### Mesh Refinement Functions

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.

### Mesh Refinement Termination Criteria

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.