Documentation

### This is machine translation

Mouseover text to see original. Click the button below to return to the English version of the page.

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, see `generateMesh` for mesh generation and `solvepde` for PDE solution.

## Syntax

```[u,p,e,t] = adaptmesh(g,b,c,a,f) [u,p,e,t] = adaptmesh(g,b,c,a,f,'PropertyName',PropertyValue) ```

## Description

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

`$-\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 2-D 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 name-value pairs.

PropertyValueDefault Description
`Maxt`

positive integer

`inf`

Maximum number of new triangles

`Ngen`

positive integer

`10`

Maximum number of triangle generations

`Mesh`

`p1`, `e1`, `t1`

`initmesh`

Initial mesh

`Tripick`

MATLAB® function

`pdeadworst`

Triangle selection method

`Par`

numeric

`0.5`

Function parameter

`Rmethod`

```'longest' | 'regular'```

`'longest'`

Triangle refinement method

`Nonlin`

```'on' | 'off'```

`'off'`

Use nonlinear solver

`Toln`

numeric

`1e-4`

Nonlinear tolerance

`Init``u0``0`

Nonlinear initial value

`Jac````'fixed | 'lumped' | 'full'````'fixed'`

Nonlinear Jacobian calculation

`norm`

numeric``` | inf```

`inf`

Nonlinear residual norm

`MesherVersion`

`'R2013a'` | `'preR2013a'`

`'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 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 `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`.

## Examples

collapse all

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.

```[u,p,e,t]=adaptmesh('cirsg','cirsb',1,0,0,'maxt',500,... 'tripick','pdeadworst','ngen',inf); ```
```Number of triangles: 197 Number of triangles: 201 Number of triangles: 216 Number of triangles: 233 Number of triangles: 254 Number of triangles: 265 Number of triangles: 313 Number of triangles: 344 Number of triangles: 417 Number of triangles: 475 Number of triangles: 629 Maximum number of triangles obtained. ```
```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 ```
`size(t,2)`
```ans = 629 ```

The maximum absolute error is 0.0028, with 629 triangles.

`pdemesh(p,e,t)` Test how many refinements you need with a uniform triangle net.

```[p,e,t]=initmesh('cirsg'); [p,e,t]=refinemesh('cirsg',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.0121 ```
`size(t,2)`
```ans = 788 ```
```[p,e,t]=refinemesh('cirsg',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.0078 ```
`size(t,2)`
```ans = 3152 ```
`pdemesh(p,e,t)` Uniform refinement with 3152 triangles produces an error of 0.0078. This error is over three times as large as that produced by the adaptive method (0.0028) with many fewer triangles (629). Typically, a problem with regular solution has an $O\left({h}^{2}\right)$ error. However, this solution is singular since $u\approx {r}^{1/3}$ at the origin.

## Diagnostics

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`

## Algorithms

### Mesh Refinement for Improving Solution Accuracy

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

### Error Estimate for 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, 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 $‖\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 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)

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

### Mesh Refinement Functions

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.

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

 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.

##### Support Get trial now