Main Content

Create adaptive 2-D mesh and solve PDE

**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.**

`[`

generates an adaptive `u`

,`p`

,`e`

,`t`

] = adaptmesh(`g`

,`b`

,`c`

,`a`

,`f`

)`[p,e,t]`

mesh and returns the solution
`u`

for an elliptic 2-D 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.

Upon termination, the function issues one of these messages:

`Adaption completed`

. (This means that the`Tripick`

function returned zero triangles to refine.)`Maximum number of triangles obtained`

.`Maximum number of refinement passes obtained`

.

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 the 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 maximal 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

Find the number of triangles.

size(t,2)

ans = 699

Plot the mesh.

pdemesh(p,e,t)

Test how many refinements you need with a uniform triangle mesh.

```
[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

Find the number of triangles in this case.

size(t,2)

ans = 816

Refine the mesh one more time. The maximal absolute error for uniform meshing is still greater than for adaptive meshing.

```
[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

Find the number of triangles in this case.

size(t,2)

ans = 3264

Plot the mesh.

pdemesh(p,e,t)

Uniform refinement with more triangles produces a larger error. Typically, a problem with regular solution has an $$O({h}^{2})$$ error. However, this solution is singular since $$u\approx {r}^{1/3}$$ at the origin.

`g`

— Geometry descriptiondecomposed geometry matrix | geometry function | handle to geometry function

Geometry description, specified as a decomposed geometry matrix, a geometry
function, or a handle to the geometry function. For details about a decomposed geometry
matrix, see `decsg`

. For details about a geometry function,
see Parametrized Function for 2-D Geometry Creation.

A geometry function must return the same result for the same input arguments in every function call. Thus, it must not contain functions and expressions designed to return a variety of results, such as random number generators.

**Data Types: **`double`

| `char`

| `string`

| `function_handle`

`b`

— Boundary conditionsboundary matrix | boundary file

Boundary conditions, specified as a boundary matrix or boundary file. Pass a boundary file as a function handle or as a file name. Typically, you export a boundary matrix from the PDE Modeler app.

**Example: **`b = 'circleb1'`

, `b = "circleb1"`

, or
`b = @circleb1`

**Data Types: **`double`

| `char`

| `string`

| `function_handle`

`c`

— PDE coefficientscalar | matrix | character vector | character array | string scalar | string vector | coefficient function

PDE coefficient, specified as a scalar, matrix, character vector, character array,
string scalar, string vector, or coefficient function. `c`

represents
the *c* coefficient in the scalar PDE

$$-\nabla \cdot \left(c\nabla u\right)+au=f$$

or in the system of PDEs

$$-\nabla \cdot \left(c\otimes \nabla u\right)+au=f$$

The coefficients `c`

, `a`

, and
`f`

can depend on the solution `u`

if you use the
nonlinear solver by setting the value of `'Nonlin'`

to
`'on'`

. The coefficients cannot be functions of the time
`t`

.

**Example: **`'cosh(x+y.^2)'`

**Data Types: **`double`

| `char`

| `string`

| `function_handle`

`a`

— PDE coefficientscalar | matrix | character vector | character array | string scalar | string vector | coefficient function

PDE coefficient, specified as a scalar, matrix, character vector, character array,
string scalar, string vector, or coefficient function. `a`

represents
the *a* coefficient in the scalar PDE

$$-\nabla \cdot \left(c\nabla u\right)+au=f$$

or in the system of PDEs

$$-\nabla \cdot \left(c\otimes \nabla u\right)+au=f$$

The coefficients `c`

, `a`

, and
`f`

can depend on the solution `u`

if you use the
nonlinear solver by setting the value of `'Nonlin'`

to
`'on'`

. The coefficients cannot be functions of the time
`t`

.

**Example: **`2*eye(3)`

**Data Types: **`double`

| `char`

| `string`

| `function_handle`

`f`

— PDE coefficientscalar | matrix | character vector | character array | string scalar | string vector | coefficient function

PDE coefficient, specified as a scalar, matrix, character vector, character array,
string scalar, string vector, or coefficient function. `f`

represents
the *f* coefficient in the scalar PDE

$$-\nabla \cdot \left(c\nabla u\right)+au=f$$

or in the system of PDEs

$$-\nabla \cdot \left(c\otimes \nabla u\right)+au=f$$

The coefficients `c`

, `a`

, and
`f`

can depend on the solution `u`

if you use the
nonlinear solver by setting the value of `'Nonlin'`

to
`'on'`

. The coefficients cannot be function of the time
`t`

.

**Example: **`char('sin(x)';'cos(y)';'tan(z)')`

**Data Types: **`double`

| `char`

| `string`

| `function_handle`

Specify optional
comma-separated pairs of `Name,Value`

arguments. `Name`

is
the argument name and `Value`

is the corresponding value.
`Name`

must appear inside quotes. You can specify several name and value
pair arguments in any order as
`Name1,Value1,...,NameN,ValueN`

.

```
[u,p,e,t] =
adaptmesh(g,'cirsb',1,0,0,'Maxt',500,'Tripick','pdeadworst','Ngen',Inf)
```

`Maxt`

— Maximum number of new triangles`Inf`

(default) | positive integerMaximum number of new triangles, specified as the comma-separated pair consisting
of `'Maxt'`

and a positive integer.

**Data Types: **`double`

`Ngen`

— Maximum number of triangle generations10 (default) | positive integer

Maximum number of triangle generations, specified as the comma-separated pair
consisting of `'Ngen'`

and a positive integer.

**Data Types: **`double`

`Mesh`

— Initial meshmesh generated by

`initmesh`

(default) | `[p,e,t]`

meshInitial mesh, specified as the comma-separated pair consisting of
`'Mesh'`

and a mesh specified by `[p,e,t]`

triples. By default, the function uses the mesh generated by the
`initmesh`

function.

**Data Types: **`double`

`Tripick`

— Triangle selection methodindices of triangles returned by

`pdeadworst`

(default) | MATLABTriangle selection method, specified as the comma-separated pair consisting of
`'Tripick'`

and a MATLAB function. By default, the function uses the indices of triangles
returned by the `pdeadworst`

function.

Given the error estimate computed by the function `pdejmps`

,
the triangle selection method identifies 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.`Par`

is the optional argument of`adaptmesh`

.

The matrices `cc`

, `aa`

, `ff`

,
and `errf`

all have *N*t columns, where
*N*t 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`

identifies triangles where`errf`

exceeds a fraction of the worst value. The default fraction value is 0.5. You can change it by specifying the`Par`

argument value when calling`adaptmesh`

.`pdeadgsc`

selects triangles using a relative tolerance criterion.

**Data Types: **`double`

`Par`

— Function parameter for triangle selection method0.5 (default) | number

Function parameter for the triangle selection method, specified as the
comma-separated pair consisting of `'Par'`

and a number between 0 and
1. The `pdeadworst`

triangle selection method uses it as its
`wlevel`

argument. The `pdeadgsc`

triangle
selection method uses it as its `tol`

argument.

**Data Types: **`double`

`Rmethod`

— Triangle refinement method`'longest'`

(default) | `'regular'`

Triangle refinement method, specified as the comma-separated pair consisting of
`'Rmethod'`

and either `'longest'`

or
`'regular'`

. For details on the refinement method, see `refinemesh`

.

**Data Types: **`char`

| `string`

`Nonlin`

— Toggle to use a nonlinear solver`'off'`

(default) | `'on'`

Toggle to use a nonlinear solver, specified as the comma-separated pair consisting
of `'Nonlin'`

and `'on'`

or
`'off'`

.

**Data Types: **`char`

| `string`

`Toln`

— Nonlinear tolerance`1e-4`

(default) | positive numberNonlinear tolerance, specified as the comma-separated pair consisting of
`'Toln'`

and a positive number. The function passes this argument
to `pdenonlin`

, which iterates until the magnitude
of the residual is less than `Toln`

.

**Data Types: **`double`

`Init`

— Nonlinear initial value0 (default) | scalar | vector of characters | vector of numbers

Nonlinear initial value, specified as the comma-separated pair consisting of
`'Init'`

and a scalar, a vector of characters, or a vector of
numbers. The function passes this argument to `pdenonlin`

, which uses it as an initial guess in its
`'U0'`

argument.

**Data Types: **`double`

`Jac`

— Nonlinear Jacobian calculation`'fixed'`

(default) | `'lumped'`

| `'full'`

Nonlinear Jacobian calculation, specified as the comma-separated pair consisting
of `'Jac'`

and either `'fixed'`

,
`'lumped'`

, or `'full'`

. The function passes this
argument to `pdenonlin`

, which uses it as an initial guess
in its `'Jacobian'`

argument.

**Data Types: **`char`

| `string`

`Norm`

— Nonlinear solver residual norm`Inf`

(default) | p value for LNonlinear solver residual norm, specified as the comma-separated pair consisting
of `'Norm'`

and `p`

value for
L^{p} norm. `p`

can be any positive real
value, `Inf`

, or `-Inf`

. The `p`

norm of a vector `v`

is `sum(abs(v)^p)^(1/p)`

. The
function passes this argument to `pdenonlin`

, which uses it as an initial guess in its
`'Norm'`

argument.

**Data Types: **`double`

| `char`

| `string`

`MesherVersion`

— Algorithm for generating initial mesh`'preR2013a'`

(default) | `'R2013a'`

Algorithm for generating initial mesh, specified as the comma-separated pair
consisting of `'MesherVersion'`

and either
`'R2013a'`

or `'preR2013a'`

. The
`'R2013a'`

algorithm runs faster, and can triangulate more
geometries than the `'preR2013a'`

algorithm. Both algorithms use
Delaunay triangulation.

**Data Types: **`char`

| `string`

`u`

— PDE solutionvector

PDE solution, returned as a vector.

If the PDE is scalar, meaning that it 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.

`p`

— Mesh points2-by-

`Np`

matrixMesh points, returned as a 2-by-`Np`

matrix. `Np`

is the number of points (nodes) in the mesh. Column `k`

of
`p`

consists of the *x*-coordinate of point
`k`

in `p(1,k)`

and the
*y*-coordinate of point `k`

in
`p(2,k)`

. For details, see Mesh Data as [p,e,t] Triples.

`e`

— Mesh edges7-by-

`Ne`

matrixMesh edges, returned as a 7-by-`Ne`

matrix. `Ne`

is the number of edges in the mesh. An edge is a pair of points in `p`

containing a boundary between subdomains, or containing an outer boundary. For details,
see Mesh Data as [p,e,t] Triples.

`t`

— Mesh elements4-by-

`Nt`

matrixMesh elements, returned as a 4-by-`Nt`

matrix.
`Nt`

is the number of triangles in the mesh.

The `t(i,k)`

, with `i`

ranging from 1 through
`end-1`

, contain indices to the corner points of element
`k`

. For details, see Mesh Data as [p,e,t] Triples. The last row,
`t(end,k)`

, contains the subdomain numbers of the elements.

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 mid-sides, 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 L-shaped membrane. In such cases, it is more efficient to refine the mesh
selectively. The selection that is based on estimates of errors in the computed solutions is
called *adaptive mesh refinement*.

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 stops34 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 (computes
an estimate of the element error contribution), the mesh refiner (selects and subdivides
elements), and the termination criteria.

The adaptation is a feedback process. It is easily applied to a larger range of problems
than those for which its design was tailored. Estimates, selection criteria, and so on must
be optimal for giving the most accurate solution at fixed cost or lowest computational
effort for a given accuracy. Such results have been proven only for model problems, but
generally, the equidistribution heuristic has been found nearly 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.

The error indicator function used in the software is an elementwise estimate of the
contribution, based on [1] and [2]. 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 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

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(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 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 mesh refinement 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 in [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.

There are two selection criteria. 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
specified 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 might never vanish, regardless of element size, making equidistribution
impossible.

[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.

You have a modified version of this example. Do you want to open this example with your edits?

You clicked a link that corresponds to this MATLAB command:

Run the command by entering it in the MATLAB Command Window. Web browsers do not support MATLAB commands.

Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .

Select web siteYou can also select a web site from the following list:

Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.

- América Latina (Español)
- Canada (English)
- United States (English)

- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)

- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)