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.

## Syntax

``[u,p,e,t] = adaptmesh(g,b,c,a,f)``
``[u,p,e,t] = adaptmesh(g,b,c,a,f,Name,Value)``

## Description

example

````[u,p,e,t] = adaptmesh(g,b,c,a,f)` generates an adaptive `[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`.```

example

````[u,p,e,t] = adaptmesh(g,b,c,a,f,Name,Value)` performs adaptive mesh generation and PDE solution for elliptic 2-D PDE problems using one or more `Name,Value` pair arguments.```

## 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 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\left({h}^{2}\right)$ error. However, this solution is singular since $u\approx {r}^{1/3}$ at the origin.

## Input Arguments

collapse all

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`

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`

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`

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`

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`

### Name-Value Arguments

Specify optional pairs of arguments as `Name1=Value1,...,NameN=ValueN`, where `Name` is the argument name and `Value` is the corresponding value. Name-value arguments must appear after other arguments, but the order of the pairs does not matter.

Before R2021a, use commas to separate each name and value, and enclose `Name` in quotes.

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

Maximum number of new triangles, specified as the comma-separated pair consisting of `'Maxt'` and a positive integer.

Data Types: `double`

Maximum number of triangle generations, specified as the comma-separated pair consisting of `'Ngen'` and a positive integer.

Data Types: `double`

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

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

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`

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`

Toggle to use a nonlinear solver, specified as the comma-separated pair consisting of `'Nonlin'` and `'on'` or `'off'`.

Data Types: `char` | `string`

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

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`

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`

Nonlinear solver residual norm, specified as the comma-separated pair consisting of `'Norm'` and `p` value for Lp 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`

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`

## Output Arguments

collapse all

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.

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

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

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

## Algorithms

collapse all

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

### Error Estimate for FEM Solution

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

### 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 might never vanish, regardless of element size, making equidistribution impossible.

## References

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

## Version History

Introduced before R2006a