An Introduction to MMPDElab

TL;DR

This paper introduces MMPDElab, a MATLAB package designed for adaptive mesh movement and finite element solutions of second-order PDEs across various dimensions, utilizing simplicial meshes.

Contribution

It presents MMPDElab, a new MATLAB tool for adaptive mesh movement and finite element solutions of PDEs using simplicial meshes in multiple dimensions.

Findings

Supports adaptive mesh movement in 1D, 2D, and 3D.

Enables finite element solutions of second-order PDEs.

Utilizes simplicial meshes for flexibility.

Abstract

This article presents an introduction to MMPDElab, a package written in MATLAB for adaptive mesh movement and adaptive moving mesh P1 finite element solution of second-order partial different equations having continuous solutions in one, two, and three spatial dimensions. MMPDElab uses simplicial meshes.

Equations

$$du={\left[(\nabla u^{(1)}{)}^T,...,(\nabla u^{(npde)}{)}^T\right]}_{N_v\times (d\ast npde)},$$

$$du(i,:)=\left[\frac{\partial u^{(1)}}{\partial x_{1}},...,\frac{\partial u^{(1)}}{\partial x_{d}},...,\frac{\partial u^{(npde)}}{\partial x_{1}},...,\frac{\partial u^{(npde)}}{\partial x_{d}}\right](\mbox{\boldmath$x$}_{i}),\quad i=1,...,N_{v}.$$

$$M(i,:)=\left[\mathbb{M}_{11},...,\mathbb{M}_{d1},...,\mathbb{M}_{1d},...,\mathbb{M}_{dd}\right](\mbox{\boldmath$x$}_{i}),\quad i=1,...,N_{v}.$$

$$\displaystyle\sum_{i=1}^{npde}\int_{\Omega}F_{i}(\nabla u,u,u_{t},\nabla v^{(i)},v^{(i)},\mbox{\boldmath$x$},t)d\mbox{\boldmath$x$}+\sum_{i=1}^{npde}\int_{\Gamma_{N}^{(i)}}G_{i}(\nabla u,u,v^{(i)},\mbox{\boldmath$x$},t)ds=0,$$

$$\forall v^{(i)}\in V^{(i)},i=1,...,npde,0<t\le T$$

$$\displaystyle R_{i}(u,\mbox{\boldmath$x$},t)=0,\qquad\text{ on }\Gamma_{D}^{(i)},\quad i=1,...,npde$$

$$u_t=ϵu_{xx}-u{u}_x,x\in \Omega \equiv (0,1),t\in (0,1]$$

$$u(t,x)=u_{exact}(t,x),x\text{ on }\partial \Omega ,t\in (0,1]$$

$$u(0,x)=u_{exact}(0,x),x\in \Omega$$

$$u_{exact}(t,x)=\frac{0.1e^{\frac{-x+0.5-4.95t}{20ϵ}}+0.5e^{\frac{-x+0.5-0.75t}{4ϵ}}+e^{\frac{-x+0.375}{2ϵ}}}{e^{\frac{-x+0.5-4.95t}{20ϵ}}+e^{\frac{-x+0.5-0.75t}{4ϵ}}+e^{\frac{-x+0.375}{2ϵ}}}.$$

$${\int }_{\Omega }(u_{t}v+ϵu_x{v}_x+u{u}_{x}v)dx=0,\forall v\in V\equiv H_0^1(\Omega ).$$

$$u_t=u_{xx}+u_{yy}+(13{\pi }^2-1)u_{exact}(t,x,y),(x,y)\in \Omega \equiv (0,1)\times (0,1),t\in (0,1]$$

$$\{\begin{array}{ll}u(t,x,y)=0,& (x,y)\text{ on }x=0\text{ and }y=0,t\in (0,1]\\ \frac{\partial u}{\partial x}=2\pi e^{-t}\sin (3\pi y),& (x,y)\text{ on }x=1,t\in (0,1]\\ \frac{\partial u}{\partial y}=-3\pi e^{-t}\sin (2\pi x),& (x,y)\text{ on }y=1,t\in (0,1]\end{array}$$

$$u(0,x,y)=u_{exact}(0,x,y),(x,y)\in \Omega .$$

$$u_{exact}(t,x,y)=e^{-t}\sin (2\pi x)\sin (3\pi y).$$

$${\int }_{\Omega }\left[(u_{t}v+u_x{v}_x+u_y{v}_y)-(13{\pi }^2-1)v{u}_{exact}(t,x,y)\right]dxdy$$

$$+{\int }_0^1\left(-2\pi e^{-t}\sin (3\pi y)\right)v(1,y)dy+{\int }_0^1\left(3\pi e^{-t}\sin (2\pi x)\right)v(x,1)dx=0,\forall v\in V$$

$$\{\begin{array}{ll}{\theta }_t={\theta }_{xx}+{\theta }_{yy}+\frac{{\beta }^2}{2Le}Y{e}^{-\frac{\beta (1-\theta )}{(1-\alpha (1-\theta ))}},& (x,y)\in \Omega ,t\in (0,60]\\ Y_t=\frac{1}{Le}(Y_{xx}+Y_{yy})-\frac{{\beta }^2}{2Le}Y{e}^{-\frac{\beta (1-\theta )}{(1-\alpha (1-\theta ))}},& (x,y)\in \Omega ,t\in (0,60]\end{array}$$

$$\{\begin{array}{ll}\theta =1,Y=0,& \text{ on bfMark = 2}\\ \frac{\partial \theta }{\partial n}=0,\frac{\partial Y}{\partial n}=0,& \text{ on bfMark = 1}\\ \frac{\partial \theta }{\partial n}=-k\theta ,\frac{\partial Y}{\partial n}=0,& \text{ on bfMark = 3}\end{array}$$

$$\{\begin{array}{ll}\theta =1,Y=0,& \text{ for }x\le 7.5\\ \theta =e^{7.5-x},Y=1-e^{Le(7.5-x)},& \text{ for }x>7.5\end{array}$$

$${\int }_{\Omega }\left({\theta }_t{v}^{(1)}+{\theta }_x{v}_x^{(1)}+{\theta }_y{v}_y^{(1)}-\frac{{\beta }^2}{2Le}Y{v}^{(1)}{e}^{-\frac{\beta (1-\theta )}{(1-\alpha (1-\theta ))}}\right)dxdy+{\int }_{{\Gamma }_{\text{bfMark = 3}}}k\theta v^{(1)}ds$$

$$+{\int }_{\Omega }\left(Y_t{v}^{(2)}+\frac{1}{Le}{Y}_x{v}_x^{(2)}+\frac{1}{Le}{Y}_y{v}_y^{(2)}+\frac{{\beta }^2}{2Le}Y{v}^{(2)}{e}^{-\frac{\beta (1-\theta )}{(1-\alpha (1-\theta ))}}\right)dxdy$$

$$.=0,\forall v^{(1)},v^{(2)}\in V$$

$$-(u_{xx}+u_{yy}+u_{zz})=f,(x,y,z)\in \Omega \equiv (0,1)\times (0,1)\times (0,1)$$

$$\{\begin{array}{ll}\frac{\partial u}{\partial x}=2\pi \sin (3\pi y)\sin (\pi z),& (x,y,z)\text{ on }{\Gamma }_N\\ u=u_{exact}(x,y,z),& (x,y,z)\text{ on }{\Gamma }_D\end{array}$$

$$u_{exact}(x,y,z)=\sin (2\pi x)\sin (3\pi y)\sin (\pi z).$$

$${\int }_{\Omega }(u_x{v}_x+u_y{v}_y+u_z{v}_z-fv)dxdydz+{\int }_{{\Gamma }_N}(-2\pi \sin (3\pi y)\sin (\pi z))v(1,y,z)dydz=0,\forall v\in V$$

Full text

An Introduction to MMPDElab

Weizhang Huang Department of Mathematics, the University of Kansas, Lawrence, KS 66045, U.S.A. (whuang@ku.edu).

1 Introduction

MMPDElab is a package written in MATLAB111MATLAB*®* is a trademark of The MathWorks, Inc., Natick, MA 01760. for adaptive mesh movement and adaptive moving mesh P1 finite element solution of second-order partial different equations (PDEs) having continuous solutions in one, two, and three spatial dimensions. It uses simplicial meshes, i.e., line segments in one dimension, triangles in two dimensions, and tetrahedra in three dimensions. The adaptive mesh movement is based on the new implementation [7, 8] of the moving mesh partial differential equation (MMPDE) method [1, 2, 10, 11, 12, 13, 14]. The mesh equation is integrated using either ode45 (an explicit MATLAB ODE solver) or ode15s (an implicit MATLAB ODE solver) while physical PDEs are discretized in space using P1 conforming finite elements on moving meshes and integrated in time with the fifth-order Radau IIA method (an implicit Runge-Kutta method) with a two-step error estimator [5] for time step selection. More information on the moving mesh P1 finite element method can be found from recent applications such as those found in [3, 9, 15, 16, 17, 20, 21].

The source code of MMPDElab can be downloaded at

• •

https://whuang.ku.edu/MMPDElab/mmpdelabv1.html

• •

https://github.com/weizhanghuang/MMPDElab

The functions in MMPDElab can be grouped into three categories:

• •

Matrix operations (with names in the form Matrix_xxx)

• •

Mesh movement (with names in the form MovMesh_xxx)

• •

Moving mesh P1 finite element solution (with names in the form MovFEM_xxx)

The functions in the first category Matrix_xxx perform vectorized computation of basic matrix operations such as multiplication, inversion, and finding transposes and determinants for arrays of matrices of small size (typically $3\times 3$ or smaller). These operations are used by functions in the other two categories which will be explained in the subsequent sections.

We now introduce notation whose understanding is crucial to the use of the package. A simplicial mesh or a simplicial triangulation, $\mathcal{T}_{h}$, of $N$ elements and $N_{v}$ vertices in $d$-dimensions ($d=1$, 2, or 3) is represented in MATLAB by the matrices $X$ and $tri$, where $X$ is a matrix of size $N_{v}\times d$ containing the coordinates of the vertices and $tri$ is a matrix of size $N\times(d+1)$ listing the connectivity of the mesh. More specifically, $X(i,:)$ gives the coordinates of the $i$th vertex \mbox{\boldmathx}_{i} while $tri(j,:)$ contains the global IDs of the vertices of the $j$th element. In MMPDElab, $npde$ components of the physical solution at the vertices are given by the matrix $u$ of size $N_{v}\times npde$, i.e., $u(i,:)$ contains the values of $u$ at the $i$th vertex. Its derivatives with respect to the physical coordinate \mbox{\boldmathx}=(x_{1},x_{2},...,x_{d})^{T} are saved in the form

[TABLE]

where $u^{(k)}$ ($k=1,...,npde$) is the $k$th component of $u$ and $\nabla$ is the gradient operator. Thus,

[TABLE]

The metric tensor or the monitor function, $\mathbb{M}$, is calculated at the vertices and saved in the form

[TABLE]

That is, $M$ has the size $N_{v}\times(d\ast d)$, with each row containing the entries of a matrix of size $d\times d$. $M$ is a good example of an array of matrices of small size. It is emphasized that when a moving mesh function is called, the mesh connectivity is kept fixed while the location of the vertices varies. The user can decide whether or not to change the connectivity in between the calls.

To conclude this section, I am deeply thankful for many colleagues and former graduate students for their invaluable discussion and comments. I am particularly grateful to Dr. Lennard Kamenski who was involved in the project at the early stage.

MMPDElab is a package written in MATLAB for adaptive mesh movement and adaptive moving mesh P1 finite element solution of partial different equations having continuous solutions.

Copyright (C) 2019  Weizhang Huang (whuang@ku.edu)

MMPDElab is free software: you can redistribute it and/or modify
it under the terms of the GNU Affero General Public License as
published by the Free Software Foundation, either version 3 of the
License, or (at your option) any later version.

MMPDElab is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
See the GNU Affero General Public License at
<https://www.gnu.org/licenses/>.

2 Adaptive mesh movement

The adaptive mesh movement can be carried out by calling MovMesh() (based on the $\xi$-formulation of the MMPDE moving mesh method [7, 8]), MovMesh_XM() (based on the $x$-formulation of the MMPDE moving mesh method), or MovMesh_X() (based on the $x$-formulation of the MMPDE moving mesh method with the metric tensor $\mathbb{M}=I$, i.e., without mesh adaptation). The corresponding MMPDE is defined as a gradient flow equation of the meshing functional developed in [6] based on mesh equidistribution and alignment (with its parameters being chosen as $p=1.5$ and $\theta=1/3$). The headers of these functions read as

[Xnew,Ih,Kmin] = MovMesh(tspan,Xi_ref,X,M,tau,tri,tri_bf,nodes_fixed, ... mmpde_int_method,dt0,abstol)

[Xnew,Ih,Kmin] = MovMesh_XM(tspan,X,M,tau,tri,tri_bf,nodes_fixed, ... mmpde_int_method,dt0,abstol,Xi_ref)

[Xnew,Ih,Kmin] = MovMesh_X(tspan,X,tau,tri,tri_bf,nodes_fixed, ... mmpde_int_method,dt0,abstol,Xi_ref)

These functions integrate the corresponding moving mesh equations over a time period specified by tspan. All of the meshes, X (the current mesh), Xnew (the new mesh), and Xi_ref (the reference computational mesh), are assumed to have the same number of vertices and elements and the same connectivity (specified by tri). The input and output variables are explained in the following.

• •

tspan is a vector specifying the time interval for integration.

• •

X, of size $N_{v}\times d$, contains the coordinates of vertices of the current mesh.

• •

Xi_ref, of size $N_{v}\times d$, contains the coordinates of vertices of the reference computational mesh. This mesh, typically chosen as the initial physical mesh, is a mandatory input for MovMesh() but is optional for MovMesh_XM() and MovMesh_X(). In the latter case, when Xi_ref is not supplied, the uniformity of the new physical mesh measured in the metric $\mathbb{M}$ is made with reference to an equilateral simplex.

• •

M, of size $N_{v}\times(d\ast d)$, contains the values of the metric tensor $\mathbb{M}$ at the vertices of X. More specifically, M(i,1:d$\ast$d) gives the metric tensor at the $i$th vertex, i.e., [\mathbb{M}_{11},...,\mathbb{M}_{d1},...,\mathbb{M}_{1d},...,\mathbb{M}_{dd}](\mbox{\boldmathx}_{i}).

• •

tau is the positive parameter used for adjusting the time scale of mesh movement.

• •

tri, of size $N\times(d+1)$, lists the connectivity for all meshes.

• •

tri_bf, of size $N_{bf}\times d$, specifies the boundary facets for all meshes, with each row containing the IDs of the vertices of a facet on the boundary. A boundary facet consists of a point in 1D, a line segment (with two vertices) in 2D, or a triangle (with three vertices) in 3D. tri_bf can be computed using the Matlab function freeBoundary in 2D and 3D.

• •

nodes_fixed is a vector containing the IDs of the vertices, such as corners, which are not allowed to move during the mesh movement.

• •

mmpde_int_method is an optional input variable, specifying that either ode15s (implicit) or ode45 (explicit) is used to integrate the moving mesh equation. The default is ode15s.

• •

dt0 is an optional input variable specifying the initial time stepsize that is used in the time integration of the mesh equation. The default is dt0 = (tspan(end)-tspan(1))/10.

• •

abstol is an optional input variable specifying the absolute tolerance used for time step selection in the time integration of the mesh equation. The default is asbstol = 1e-6 for ode15s and 1e-8 for ode45.

• •

Xnew, of size $N_{v}\times d$, contains the coordinates of vertices of the new mesh.

• •

Ih is an optional output variable giving the value of the meshing functional at the new mesh.

• •

Kmin is an optional output variable giving the minimal element volume.

In addition to MovMesh(), MovMesh_XM(), and MovMesh_X(), the following functions can also be used by the user.

[X,tri] = MovMesh_circle2tri(jmax) This function creates a triangular mesh (X, tri) for the unit circle. 2. 2.

[X,tri] = MovMesh_cube2tet(x,y,z) This function creates a tetrahedral mesh (X, tri) from the cuboid mesh specified by x, y, and z for a cuboid domain. Each subcuboid is divided into 6 tetrahedra. 3. 3.

V = MovMesh_freeBoundary_faceNormal(X,tri,tri_bf) This function computes the unit outward normals for the boundary facets. V has the size of $N_{bf}\times d$. 4. 4.

V = MovMesh_freeBoundary_vertexNormal(X,tri,tri_bf) This function computes the unit outward normals for the boundary vertices. V has the size of $N_{v}\times d$, with the normals for the interior vertices being set to be $[1,...,1]^{T}/\sqrt{d}$. 5. 5.

[Grad,Hessian] = MovMesh_GradHessianRecovery(u,X,tri,tri_bf) This function computes the gradient and Hessian of function $u$ at the vertices using centroid-vortex-centroid-vertex volume-weighted average. 6. 6.

Grad = MovMesh_GradKRecovery(u,X,tri,tri_bf) This function computes the gradient of function $u$ on the elements. 7. 7.

Grad = MovMesh_GradRecovery(u,X,tri,tri_bf) This function computes the gradient of function $u$ at the vertices using volume averaging. 8. 8.

fnew = MovMesh_LinInterp(f,X,QP,tri,tri_bf,useDelaunayTri) This function performs linear interpolation of f (defined on X) at query points QP using triangulation or Delaunay triangulation. useDelaunayTri is a logical variable with value true or false. 9. 9.

[X,tri] = MovMesh_MeshMerge(X1,tri1,X2,tri2) This function merges two non-overlapping meshes (X1,tri1) and (X2,tri2) which may or may not have common boundary segments. 10. 10.

[Qgeo,Qeq,Qali] = MovMesh_MeshQualMeasure(X,tri,M,Linf_norm,Xi_ref) This function computes the geometric, equidistribution, and alignment measures (in maximum norm or $L^{2}$ norm in $\xi$) for the mesh (X, tri) according to the metric tensor. Here, both Linf_norm and Xi_ref are optional input variables. 11. 11.

[Qmax,Ql2] = MovMesh_MeshQualMeasure2(X,tri,M,Xi_ref) This function computes the maximum and $L^{2}$ norm of the mesh quality measure based on a single condition combining both equidistribution and alignment. Xi_ref is an optional input variable. 12. 12.

[X,tri] = MovMesh_MeshRemoveNodes(X1,tri1,ID) This function removes the nodes listed in ID from the existing mesh (X1,tri1). 13. 13.

[XF,TriF,TriF_parent] = MovMesh_MeshUniformRefine(X,Tri,Level) This function uniformly refines a simplicial mesh (Level) times or levels. On each level, an element is refined into $2^{d}$ elements. 14. 14.

M = MovMesh_metric_arclength(u,X,tri,tri_bf) This function computes the arclength metric tensor. 15. 15.

MC = MovMesh_metric_F2C(M,Tri,Tri_parent,TriC) This function computes the metric tensor on a coarse mesh from the metric tensor defined on a fine mesh. 16. 16.

M = MovMesh_metric_intersection(M1,M2) This function computes the intersection of two symmetric and positive definite matrices. When M1 and M2 are diagonal, i.e., $\verb|M1|=\text{diag}(\alpha_{1},...,\alpha_{d})$ and $\verb|M2|=\text{diag}(\beta_{1},...,\beta_{d})$, then $\verb|M|=\text{diag}(\max(\alpha_{1},\beta_{1}),...,\max(\alpha_{d},\beta_{d}))$. The intersection of two general symmetric and positive definite matrices is defined similarly through simultaneous diagonalization. 17. 17.

M = MovMesh_metric_iso(u,X,tri,tri_bf,alpha,m) This function computes the isotropic metric tensor based on the $L^{2}$ norm or the $H^{1}$ seminorm of linear interpolation error ($m=0$ or $m=1$). 18. 18.

MM = MovMesh_metric_smoothing(M,ncycles,X,tri) This function smooths the metric tensor ncycles times by local averaging. 19. 19.

M = MovMesh_metric(u,X,tri,tri_bf,alpha,m) This function computes the metric tensor based on the $L^{2}$ norm or the $H^{1}$ seminorm of linear interpolation error ($m=0$ or $m=1$). 20. 20.

[X,tri] = MovMesh_rect2tri(x,y,job) This function creates a triangular mesh (X, tri) from the rectangular mesh specified by x and y for a rectangular domain. Each rectangle is divided into 2 (for job = 2 or 3) or 4 (for job = 1) triangles. 21. 21.

M1 = Matrix_ceil(M,beta) This function puts a ceiling on the eigenvalues of symmetric and positive definite matrix M such that $\lambda_{max}(M1)\leq\beta$.

Examples using these functions include ex1d_1.m, ex2d_1.m, ex2d_2_Lshape.m, ex2d_3_hole.m, ex2d_4_horseshoe.m, and ex3d_1.m in the subdirectory ./examples.

Troubleshooting. Occasionally one may see an error message like

Error using triangulation The coordinates of the input points must be finite values; Inf and NaN are not permitted.

Error in MovMesh>MovMesh_rhs (line 296) TR = triangulation(tri2,XI2);

when calling MovMesh(), MovMesh_XM(), or MovMesh_X(). Typically this is caused by a stability issue when integrating the MMPDE, and using a smaller initial time step dt0 (e.g., dt0 = 1e-6) may solve the problem.

The generation of initial meshes. This package includes a few functions for generating initial meshes for simple domains, such as MovMesh_circle2tri(), MovMesh_cube2tet(), and

MovMesh_rect2tri(). For complex domains, one can use MATLAB function delaunayTriangulation() or other automatic mesh generators such as DistMesh (in MATLAB) [18] and TetGen (in C++) [19].

3 Adaptive mesh movement P1 finite element solution of PDEs

This package aims to solve the system of PDEs in the weak form: find $u=[u^{(1)},...,u^{(npde)}]\in H^{1}(\Omega)\otimes\cdots\otimes H^{1}(\Omega)$ such that

[TABLE]

subject to the Dirichlet boundary conditions

[TABLE]

where for $i=1,...,npde$, $\Gamma_{D}^{(i)}$ and $\Gamma_{N}^{(i)}$ are the boundary segments corresponding to the Dirichlet and Neumann boundary conditions for $u^{(i)}$, respectively, $\Gamma_{D}^{(i)}\cup\Gamma_{N}^{(i)}=\partial\Omega$, and $V^{(i)}=\{w\in H^{1}(\Omega),\;w=0\text{ on }\Gamma_{D}^{(i)}\}$. The headers of MovFEM() (Initial-Boundary-Value-Problem solver) and MovFEM_bvp() (Boundary-Value-Problem solver) read as

[Unew,dt0,dt1] = MovFEM(t,dt,U,X,Xdot,tri,tri_bf,pdedef, ... fixed_step,relTol,absTol,direct_ls,ControlWeights)

Unew = MovFEM_bvp(U,X,tri,tri_bf,pdedef,nonlinearsolver,MaxIter,Tol)

MovFEM() integrates the system of PDEs on a moving mesh over a time step. Its input and output variables are explained in the following.

• •

t is the current time.

• •

dt is the intended time stepsize for integrating the physical PDEs.

• •

U, of size $N_{v}\times npde$, is the current solution.

• •

X, of size $N_{v}\times d$, contains the coordinates of vertices of the current mesh.

• •

Xdot, of size $N_{v}\times d$, is the nodal mesh velocity.

• •

tri, of size $N\times(d+1)$, lists the connectivity for all meshes.

• •

tri_bf, of size $N_{bf}\times d$, specifies the boundary facets for all meshes.

• •

pdedef is a structure used to define the PDE system in the weak form. It has 5 fields.

1. (i)

pdedef.bfMark, of size $N_{bf}\times 1$, is used to mark the boundary segments (boundary facets). This marking is passed to the definitions of boundary integrals and Dirichlet boundary conditions. 2. (ii)

pdedef.bftype, of size $N_{bf}\times npde$, specifies the types of boundary condition on boundary facets whose numbering is based on tri_bf. pdedef.bftype = 0 for Neumann BCs and pdedef.bftype = 1 for Dirichlet BCs. For example, pdedef.bftype(3,2) = 1 means that variable $u^{(2)}$ has a Dirichlet BC on the 3rd boundary facet while pdedef.bftype(2,1) = 0 specifies that variable $u^{(1)}$ has a Neumann BC on the 2nd boundary facet. 3. (iii)

F = pdedef.volumeInt(du, u, ut, dv, v, x, t, i) This function is used to define $F_{i}$ in the weak form (1), where $v$ and $dv$ are the test function $v^{(i)}$ and its gradient. 4. (iv)

G = pdedef.boundaryInt(du, u, v, x, t, i, bfMark) This function is used to define $G_{i}$ in the weak form (1), where $v$ is the test function $v^{(i)}$. 5. (v)

R = pdedef.dirichletRes(u, x, t, i, bfMark) This function is used to define $R_{i}$ in (2).

• •

fixed_step is an optional input logical variable, indicating whether or not a fixed step is used in time integration. The default is false.

• •

relTol and absTol are optional input variables for the relative and absolute tolerances used for time step selection. The defaults are relTol = 1e-4 and absTol = 1e-6.

• •

direct_ls is an optional input logical variable, indicating whether or not the direct sparse matrix solver is used for solving linear algebraic systems. When direct_ls = false, the BiConjugate Gradients Stabilized Method bicgstab is used. The default is true.

• •

ControlWeights is an optional input variable which is nonnegative vector of size $(N_{v}\ast npde)\times 1$ used to define the weights of the components of the solution for the error estimation used in time step selection.

• •

Unew, of size $N_{v}\times npde$, is the new solution at time t + dt0.

• •

dt0 is the time stepsize actually used to integrate the physical PDEs.

• •

dt1 is the time step size predicted for the next step.

The input and output variables for MovFEM_bvp() are similar to those of MovFEM(). The same weak form (1) and (2) is used for both IBVPs and BVPs. In the latter case, $t$ is a parameter that is not used. Here we list the variables used only in the BVP solver.

• •

nonlinearsolver is an optional input variable for the method used for solving nonlinear algebraic systems, with the choices being newtons and fsolve. The defacult is fsolve.

• •

MaxIter is an optional input variable for the maximum number of iterations allowed for the solution of nonlinear algebraic systems. The default is MaxIter = 300.

• •

Tol is an optional input variable for the tolerance used in the solution of nonlinear algebraic systems. The default is Tol = 1e-6.

The following two functions are used to compute the error when the exact solution is available in the form U = uexact(t,x,varargin).

err = MovFEM_Error_P1L2(uexact,t,X,U,tri,tri_bf,varargin) This function computes the $L^{2}$ norm of the error in P1 finite element approximation. 2. 2.

err = MovFEM_Error_P1Linf(uexact,t,X,U,tri,tri_bf,varargin) This function computes the $L^{\infty}$ norm of the error in P1 finite element approximation.

In the following we give several examples to explain how to define the PDE system through pdedef. More examples can be found in the subdirectory ./examples. A typical flow chart for those examples is shown in Fig. 1.

3.1 Burgers’ equation in 1D

This example, implemented in ex1d_burgers.m, is the IBVP of Burgers’ equation in 1D,

[TABLE]

subject to the Dirichlet boundary condition

[TABLE]

and the initial condition

[TABLE]

where $\epsilon=10^{-3}$ and

[TABLE]

The weak formulation of this example reads as

[TABLE]

The definition of this example in the code is given as

% define PDE system and BCs

% all bcs are dirichlet so no need for marking boundary segments pdedef.bfMark = ones(Nbf,1); pdedef.bftype = ones(Nbf,npde);

pdedef.volumeInt = @pdedef_volumeInt; pdedef.boundaryInt = @pdedef_boundaryInt; pdedef.dirichletRes = @pdedef_dirichletRes;

... ...

function F = pdedef_volumeInt(du, u, ut, dv, v, x, t, ipde) global epsilon; F = ut(:,1).v(:) + epsilondu(:,1).*dv(:,1) + u(:,1).*du(:,1).*v(:);

function G = pdedef_boundaryInt(du, u, v, x, t, ipde, bfMark) G = zeros(size(x,1),1);

function Res = pdedef_dirichletRes(u, x, t, ipde, bfMark) Res = u - uexact(t, x);

3.2 The heat equation in 2D

This example, implemented in ex2d_heat.m, is the IBVP for the heat equation in 2D,

[TABLE]

subject to the boundary conditions

[TABLE]

and the initial condition

[TABLE]

This problem has the exact solution

[TABLE]

The weak formulation reads as

[TABLE]

where $V=\{w\in H^{1}(\Omega),w=0\text{ on }x=0\text{ and }y=0\}$. The definition of this example in the code is given as

% define PDE system and BCs

% mark boundary segments pdedef.bfMark = ones(Nbf,1); % for y = 0 (b1) Xbfm = (X(tri_bf(:,1),:)+X(tri_bf(:,2),:))*0.5; pdedef.bfMark(Xbfm(:,1)<1e-8) = 4; % for x = 0 (b4) pdedef.bfMark(Xbfm(:,1)>1-1e-8) = 2; % for x = 1 (b2) pdedef.bfMark(Xbfm(:,2)>1-1e-8) = 3; % for y = 1 (b3)

% define boundary types pdedef.bftype = ones(Nbf,npde); % for neumann bcs: pdedef.bftype(pdedef.bfMark==2|pdedef.bfMark==3,npde) = 0;

pdedef.volumeInt = @pdedef_volumeInt; pdedef.boundaryInt = @pdedef_boundaryInt; pdedef.dirichletRes = @pdedef_dirichletRes; ... ...

function F = pdedef_volumeInt(du, u, ut, dv, v, x, t, ipde) F = (13pipi-1)*uexact(t,x); F = ut(:,1).*v(:)+du(:,1).*dv(:,1)+du(:,2).*dv(:,2)-F.*v(:);

function G = pdedef_boundaryInt(du, u, v, x, t, ipde, bfMark) G = zeros(size(x,1),1); ID = find(bfMark==2); G(ID) = -2piexp(-t)sin(3pix(ID,2)).v(ID); ID = find(bfMark==3); G(ID) = 3piexp(-t)sin(2pi*x(ID,1)).*v(ID);

function Res = pdedef_dirichletRes(u, x, t, ipde, bfMark) Res = zeros(size(x,1),1); ID = find(bfMark==1|bfMark==4); Res(ID) = u(ID,1)-0.0;

3.3 A combustion model in 2D

This example, implemented in ex2d_combustion.m, is the IBVP for a combustion model (a system of two PDEs) in 2D (see [4]),

[TABLE]

subject to the boundary conditions

[TABLE]

and the initial condition

[TABLE]

where $\Omega$ is shown and the boundary segments are marked as in Fig. 2 and $Le=1$, $\alpha=0.8$, $\beta=10$, and $k=0.1$. The analytical expression of the exact solution is not available. The weak formulation reads as

[TABLE]

where $V=\{v\in H^{1}(\Omega),v=0\text{ on bfMark = 2}\}$. The definition of this example in the code is given as

% define PDE system and BCs

pdedef.bfMark = ones(Nbf,1); Xbfm = (X(tri_bf(:,1),:)+X(tri_bf(:,2),:))*0.5; pdedef.bfMark(Xbfm(:,1) < 1e-8) = 2; pdedef.bfMark(abs(Xbfm(:,1)-15) < 1e-8) = 3; pdedef.bfMark(abs(Xbfm(:,1)-30) < 1e-8) = 3; pdedef.bfMark((abs(Xbfm(:,2)-4) < 1e-8) & ... (Xbfm(:,1) > 15 & Xbfm(:,1) < 30)) = 3; pdedef.bfMark((abs(Xbfm(:,2)-12) < 1e-8) & ... (Xbfm(:,1) > 15 & Xbfm(:,1) < 30)) = 3;

pdedef.bftype = ones(Nbf,npde); pdedef.bftype(pdedef.bfMark==1|pdedef.bfMark==3,:) = 0; pdedef.volumeInt = @pdedef_volumeInt; pdedef.boundaryInt = @pdedef_boundaryInt; pdedef.dirichletRes = @pdedef_dirichletRes;

... ...

function F = pdedef_volumeInt(du, u, ut, dv, v, x, t, ipde)

beta = 10;
alpha = 0.8;
Le = 1;
w = beta^2/(2*Le)*u(:,2).*exp(-beta*(1-u(:,1))./(1-alpha*(1-u(:,1))));
if (ipde==1)
   F = ut(:,1).*v+du(:,1).*dv(:,1)+du(:,2).*dv(:,2) - w.*v;
else
   F = ut(:,2).*v+(du(:,3).*dv(:,1)+du(:,4).*dv(:,2))/Le + w.*v;
end

function G = pdedef_boundaryInt(du, u, v, x, t, ipde, bfMark)

k = 0.1; G = zeros(size(x,1),1); if ipde==1 ID = find(bfMark==3); G(ID) = k*u(ID,1).*v(ID); end

function Res = pdedef_dirichletRes(u, x, t, ipde, bfMark)

Res = zeros(size(x,1),1); ID = find(bfMark==2); if (ipde==1) Res(ID) = u(ID,1)-1; else Res(ID) = u(ID,2)-0; end

3.4 Poisson’s equation in 3D

This example, implemented in ex3d_poisson.m, is the BVP for Poisson’s equation in 3D,

[TABLE]

subject to the boundary conditions

[TABLE]

where $\Gamma_{N}=\{x=1,0<y<1,0<z<1\}$, $\Gamma_{D}=\partial\Omega\setminus\Gamma_{N}$, and $f$ is chosen such that the exact solution of this example is

[TABLE]

The weak formulation of this example reads as

[TABLE]

where $V=\{w\in H^{1}(\Omega),w=0\text{ on }\Gamma_{D}\}$. The definition of this example in the code is given as

% define PDE system and BCs

pdedef.bfMark = ones(Nbf,1); Xbfm = (X(tri_bf(:,1),:)+X(tri_bf(:,2),:)+X(tri_bf(:,3),:))/3; pdedef.bfMark(Xbfm(:,1)>1-1e-8) = 2; % for x=1 pdedef.bftype = ones(Nbf,npde); pdedef.bftype(pdedef.bfMark==2,npde) = 0; % neumann bc for x=1

pdedef.volumeInt = @pdedef_volumeInt; pdedef.boundaryInt = @pdedef_boundaryInt; pdedef.dirichletRes = @pdedef_dirichletRes;

... ...

function F = pdedef_volumeInt(du, u, ut, dv, v, x, t, ipde) F = 14pi^2sin(2pix(:,1)).sin(3pi*x(:,2)).sin(pix(:,3)); F = du(:,1).*dv(:,1)+du(:,2).*dv(:,2)+du(:,3).*dv(:,3)-F.*v(:);

function G = pdedef_boundaryInt(du, u, v, x, t, ipde, bfMark) G = zeros(size(x,1),1); ID = find(bfMark==2); G(ID) = -2pisin(3pix(ID,2)).sin(pix(ID,3)).*v(ID);

function Res = pdedef_dirichletRes(u, x, t, ipde, bfMark) Res = u(:,1) - uexact(t,x);