Adaptive 2-D 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,
generateMesh for mesh generation
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
perform adaptive mesh generation and PDE solution for elliptic
problems with 2-D geometry. Optional arguments are given as
property name/property value pairs.
The function produces a solution u to the elliptic scalar PDE problem
for (x,y) ∊ Ω, or the elliptic system PDE problem
with the problem geometry and boundary conditions given by
b. The mesh is described by the
The solution u is represented as the solution vector
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
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
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
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
b describes the boundary conditions of the PDE problem.
The adapted triangular mesh of the PDE problem is given by the mesh data
details on the mesh data representation, see Mesh Data as [p,e,t] Triples.
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 name-value pairs.
Maximum number of new triangles
Maximum number of triangle generations
Triangle selection method
Triangle refinement method
Use nonlinear solver
Nonlinear initial value
Nonlinear Jacobian calculation
Nonlinear residual norm
Algorithm for generating initial mesh
Par is passed to the
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
the input mesh data. This triangular mesh is used as starting mesh
for the adaptive algorithm. For details on the mesh data representation,
initmesh. If no initial mesh
is provided, the result of a call to
no options is used as the initial mesh.
The triangle selection method,
Tripick, is a user-definable 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
t represent the current
generation of triangles;
ff are the current coefficients for the PDE problem, expanded to the
u is the current solution;
the computed error estimate; and
par, the function parameter, is given to
adaptmesh as optional argument. The matrices
errf all have
Nt columns, where Nt is the current number of
triangles. The numbers of rows in
ff are exactly the same as the input arguments
errf has one row for each
equation in the system. The two standard triangle selection methods are
selects triangles where
errf exceeds a fraction (the default is 0.5) of the
worst value, and
pdeadgsc selects triangles using a relative tolerance
The refinement method is
regular. For details
on the refinement method, see
MesherVersion property chooses the algorithm for mesh generation. The
'R2013a' algorithm runs faster and can triangulate more geometries than
'preR2013a' algorithm. Both algorithms use Delaunay
The adaptive algorithm can also solve nonlinear PDE problems.
For nonlinear PDE problems, the
must be set to
on. The nonlinear tolerance
nonlinear initial value
u0, nonlinear Jacobian
Jac, and nonlinear residual norm
passed to the nonlinear solver
Solve the Laplace equation over a circle sector, with Dirichlet boundary conditions u = cos(2/3atan2( y , x )) along the arc and u = 0 along the straight lines, and compare the resulting solution to the exact solution. Set options so that
adaptmesh refines the triangles using the worst error criterion until it obtains a mesh with at least 500 triangles.
c45 = cos(pi/4); L1 = [2 -c45 0 c45 0 1 0 0 0 0]'; L2 = [2 -c45 0 -c45 0 1 0 0 0 0]'; C1 = [1 -c45 c45 -c45 -c45 1 0 0 0 1]'; C2 = [1 c45 c45 -c45 c45 1 0 0 0 1]'; C3 = [1 c45 -c45 c45 c45 1 0 0 0 1]'; g = [L1 L2 C1 C2 C3]; [u,p,e,t] = adaptmesh(g,'cirsb',1,0,0,'maxt',500,... 'tripick','pdeadworst','ngen',inf);
Number of triangles: 204 Number of triangles: 208 Number of triangles: 217 Number of triangles: 230 Number of triangles: 265 Number of triangles: 274 Number of triangles: 332 Number of triangles: 347 Number of triangles: 460 Number of triangles: 477 Number of triangles: 699 Maximum number of triangles obtained.
Find the maximum absolute error.
x = p(1,:); y = p(2,:); exact = ((x.^2 + y.^2).^(1/3).*cos(2/3*atan2(y,x)))'; max(abs(u - exact))
ans = 0.0028
The number of triangles in this case is:
ans = 699
Plot the mesh.
Test how many refinements you need with a uniform triangle net.
[p,e,t] = initmesh(g); [p,e,t] = refinemesh(g,p,e,t); u = assempde('cirsb',p,e,t,1,0,0); x = p(1,:); y = p(2,:); exact = ((x.^2 + y.^2).^(1/3).*cos(2/3*atan2(y,x)))'; max(abs(u - exact))
ans = 0.0116
ans = 816
[p,e,t] = refinemesh(g,p,e,t); u = assempde('cirsb',p,e,t,1,0,0); x = p(1,:); y = p(2,:); exact = ((x.^2+y.^2).^(1/3).*cos(2/3*atan2(y,x)))'; max(abs(u - exact))
ans = 0.0075
ans = 3264
Plot the mesh.
Uniform refinement with more triangles produces a larger error. Typically, a problem with regular solution has an error. However, this solution is singular since at the origin.
Upon termination, one of the following messages is displayed:
Adaption completed (This means
Tripick function returned zero triangles
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 2-D 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
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 L-shaped 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.
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
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
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. , . For Poisson's equation –Δu = f on Ω, the following error estimate for the FEM-solution uh holds in the L2-norm :
where h = h(x) is the local mesh size, and
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 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 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)|
where is the unit normal of the edge τ and the braced term is
the jump in flux across the element edge. The L2
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 ill-conditioning, such problems 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 is used by Rosenberg-Stenger , 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, 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.
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.
 Johnson, C. Numerical Solution of Partial Differential Equations by the Finite Element Method. Lund, Sweden: Studentlitteratur, 1987.
 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.
 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.