adaptmesh

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 (c \nabla u) + a u = f$

for (x,y) ∊ Ω, or the elliptic system PDE problem

$- \nabla \cdot (c \otimes \nabla u) + a u = 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.

b describes the boundary conditions of the PDE problem.

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.

Property	Value	Default	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

Adaptive Mesh Generation and Mesh Refinement

Open Live Script

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:

size(t,2)

ans = 699

Plot the mesh.

pdemesh(p,e,t)

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

size(t,2)

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

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.

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. [1], [2]. For Poisson's equation –Δu = f on Ω, the following error estimate for the FEM-solution u_h holds in the L₂-norm $‖ \cdot ‖$ :

$‖ \nabla (u - u_{h}) ‖ \leq α ‖ h f ‖ + β D_{h} (u_{h})$

where h = h(x) is the local mesh size, and

$D_{h} (v) = {(\sum_{τ \in E_{1}} h_{τ}^{2} {[\frac{\partial v}{\partial n_{τ}}]}^{2})}^{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)

$E (K) = α {‖ h (f - a u) ‖}_{K} + β {(\frac{1}{2} \sum_{τ \in \partial K} h_{τ}^{2} {(n_{τ} \cdot c \nabla u_{h})}^{2})}^{1 / 2}$

where $n_{τ}$ is the unit normal of the edge τ and the braced term is the jump in flux across the element edge. The L₂ 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 [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 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.

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.

Documentation