Adaptive 2D mesh generation and PDE solution
This page describes the legacy workflow. New features might
not be compatible with the legacy workflow. In the recommended workflow,
see generateMesh
for mesh generation
and solvepde
for PDE solution.
[u,p,e,t] = adaptmesh(g,b,c,a,f)
[u,p,e,t]
= adaptmesh(g,b,c,a,f,'PropertyName',PropertyValue)
[u,p,e,t] = adaptmesh(g,b,c,a,f)
and
[u,p,e,t]
= adaptmesh(g,b,c,a,f,'PropertyName',PropertyValue)
perform adaptive mesh generation and PDE solution for elliptic
problems with 2D geometry. Optional arguments are given as
property name/property value pairs.
The function produces a solution u to the elliptic scalar PDE problem
$$\nabla \cdot \left(c\nabla u\right)+au=f$$
for (x,y) ∊ Ω, or the elliptic system PDE problem
$$\nabla \cdot \left(c\otimes \nabla u\right)+au=f$$
with the problem geometry and boundary conditions given by g
and
b
. The mesh is described by the p
,
e
, and t
matrices.
The solution u is represented as the solution vector u
.
If the PDE is scalar, meaning that is has only one equation, then u
is a
column vector representing the solution u at each node in the mesh. If the
PDE is a system of N > 1 equations, then u
is a
column vector with N*Np
elements, where Np
is the number
of nodes in the mesh. The first Np
elements of u
represent the solution of equation 1, the next Np
elements represent the
solution of equation 2, and so on.
The algorithm works by solving a sequence of PDE problems using refined triangular meshes. The
first triangular mesh generation is provided as an optional argument to
adaptmesh
or obtained by a call to initmesh
without
options. Subsequent generations of triangular meshes are obtained by solving the PDE problem,
computing an error estimate, selecting a set of triangles based on the error estimate, and
then refining the triangles. The solution to the PDE problem is then recomputed. The loop
continues until the triangle selection method selects no further triangles, until the maximum
number of triangles is attained, or until the maximum number of triangle generations is
reached.
g
describes the geometry of the PDE problem. g
can be a
decomposed geometry matrix, the name of a geometry file, or a function handle to a geometry
file. For details, see Three Ways to Create 2D Geometry.
b
describes the boundary conditions of the PDE problem. For details, see
Boundary Conditions by Writing Functions.
The adapted triangular mesh of the PDE problem is given by the mesh data p
,
e
, and t
. For
details on the mesh data representation, see Mesh Data as [p,e,t] Triples.
The coefficients c
, a
, and f
of the
PDE problem can be given in a wide variety of ways. In the context of
adaptmesh
, the coefficients can depend on u
if the
nonlinear solver is enabled using the property nonlin
. The coefficients
cannot depend on t
, the time.
adaptmesh
accepts the following property namevalue pairs.
Property  Value  Default  Description 

Maxt  positive integer  inf  Maximum number of new triangles 
Ngen  positive integer  10  Maximum number of triangle generations 
Mesh 
 initmesh
 Initial mesh 
Tripick  MATLAB^{®} function  pdeadworst  Triangle selection method 
Par  numeric  0.5  Function parameter 
Rmethod 
 'longest'  Triangle refinement method 
Nonlin 
 'off'  Use nonlinear solver 
Toln  numeric  1e4  Nonlinear tolerance 
Init  u0  0  Nonlinear initial value 
Jac  'fixed  'lumped' 
'full'  'fixed'  Nonlinear Jacobian calculation 
norm  numeric  inf  Nonlinear residual norm 

 'preR2013a'  Algorithm for generating initial mesh 
Par
is passed to the Tripick
function,
which is described later. Normally it is used as tolerance of how
well the solution fits the equation.
No more than Ngen
successive refinements
are attempted. Refinement is also stopped when the number of triangles
in the mesh exceeds Maxt
.
p1
, e1
, and t1
are
the input mesh data. This triangular mesh is used as starting mesh
for the adaptive algorithm. For details on the mesh data representation,
see initmesh
. If no initial mesh
is provided, the result of a call to initmesh
with
no options is used as the initial mesh.
The triangle selection method, Tripick
, is a userdefinable triangle
selection method. Given the error estimate computed by the function
pdejmps
, the triangle selection method selects the triangles to be
refined in the next triangle generation. The function is called using the arguments
p
, t
, cc
, aa
,
ff
, u
, errf
, and
par
. p
and t
represent the current
generation of triangles; cc
, aa
, and
ff
are the current coefficients for the PDE problem, expanded to the
triangle midpoints; u
is the current solution; errf
is
the computed error estimate; and par
, the function parameter, is given to
adaptmesh
as optional argument. The matrices cc
,
aa
, ff
, and errf
all have
Nt columns, where Nt is the current number of
triangles. The numbers of rows in cc
, aa
, and
ff
are exactly the same as the input arguments c
,
a
, and f
. errf
has one row for each
equation in the system. The two standard triangle selection methods are
pdeadworst
and pdeadgsc
. pdeadworst
selects triangles where errf
exceeds a fraction (the default is 0.5) of the
worst value, and pdeadgsc
selects triangles using a relative tolerance
criterion.
The refinement method is longest
or regular
. For details
on the refinement method, see refinemesh
.
The MesherVersion
property chooses the algorithm for mesh generation. The
'R2013a'
algorithm runs faster and can triangulate more geometries than
the 'preR2013a'
algorithm. Both algorithms use Delaunay
triangulation.
The adaptive algorithm can also solve nonlinear PDE problems.
For nonlinear PDE problems, the Nonlin
parameter
must be set to on
. The nonlinear tolerance Toln
,
nonlinear initial value u0
, nonlinear Jacobian
calculation Jac
, and nonlinear residual norm Norm
are
passed to the nonlinear solver pdenonlin
.
Upon termination, one of the following messages is displayed:
Adaption completed
(This means
that the Tripick
function returned zero triangles
to refine.)
Maximum number of triangles obtained
Maximum number of refinement passes obtained
Partial Differential Equation Toolbox™ provides the refinemesh
function for global, uniform mesh
refinement for 2D geometries. 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 can be obtained by
extrapolation.
The solutions of equations often have geometric features such as localized strong
gradients. An example of engineering importance in elasticity is the stress concentration
occurring at reentrant corners, such as the MATLAB Lshaped membrane. In such cases, it is more economical to refine the mesh
selectively, that is, only where it is needed. The selection that is based on estimates of
errors in the computed solutions is called 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, and so on 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 must be chosen so 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 bounds might not exist, while the refinement scheme is still well defined and works 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. [1], [2]. For Poisson's equation –Δu = f on Ω, the following error estimate for the FEMsolution u_{h} holds in the L_{2}norm $$\Vert \cdot \Vert $$:
$$\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 the edge τ, h_{τ} is the length of the edge τ, and the sum runs over E_{i}, 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 the element K by summing the contributions from its edges.
The general form of the error indicator function for the elliptic equation
–∇ · (c∇u) + au = f  (1) 
is
$$E\left(K\right)=\alpha {\Vert h\left(fau\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 the edge τ and the braced term is
the jump in flux across the element edge. The L_{2}
norm is computed over the element K. The pdejmps
function computes this error indicator.
Partial Differential Equation Toolbox software is geared to elliptic problems. For reasons of accuracy and illconditioning, such problems require the elements not to deviate too much from being equilateral. Thus, even at essentially onedimensional 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 is used by RosenbergStenger [3], 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
userdefined dimensionless tolerance. The comparison with the tolerance is properly scaled
with respect to domain, solution size, and so on.
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.
[1] Johnson, C. Numerical Solution of Partial Differential Equations by the Finite Element Method. Lund, Sweden: Studentlitteratur, 1987.
[2] Johnson, C., and K. Eriksson. Adaptive Finite Element Methods for Parabolic Problems I: A Linear Model Problem. SIAM J. Numer. Anal, 28, (1991), pp. 43–77.
[3] Rosenberg, I.G., and F. Stenger, A lower bound on the angles of triangles constructed by bisecting the longest side. Mathematics of Computation. Vol 29, Number 10, 1975, pp 390–395.