Documentation

## Boundary Conditions by Writing Functions

### Note

THIS PAGE DESCRIBES THE LEGACY WORKFLOW. New features might not be compatible with the legacy workflow. For the corresponding step in the recommended workflow, see Specify Boundary Conditions.

### About Boundary Conditions by Writing Functions

This section shows how to express boundary conditions for 2-D geometry using the legacy function syntax. However, the recommended way to express boundary conditions is to use Specify Boundary Conditions.

To use this legacy syntax, write the functions using the templates in Boundary Conditions for Scalar PDE or Boundary Conditions for PDE Systems.

### Boundary Conditions for Scalar PDE

For a scalar PDE, some boundary segments can have Dirichlet conditions, and some boundary segments can have generalized Neumann conditions.

Dirichlet boundary conditions are

hu = r,

where h and r can be functions of x, y, the solution u, the edge segment index, and, for parabolic and hyperbolic equations, time.

Generalized Neumann boundary conditions are $\stackrel{\to }{n}\text{\hspace{0.17em}}·\text{\hspace{0.17em}}\left(c\nabla u\right)+qu=g$ on ∂Ω.

$\stackrel{\to }{n}$ is the outward unit normal. g and q are functions defined on ∂Ω, and can be functions of x, y, the solution u, the edge segment index, and, for parabolic and hyperbolic equations, time.

To write a function file, say `pdebound.m`, use the following syntax:

`[qmatrix,gmatrix,hmatrix,rmatrix] = pdebound(p,e,u,time)`

Your function returns matrices `qmatrix`, `gmatrix`, `hmatrix`, and `rmatrix`, based on these inputs:

• `p` — Points in the mesh (Mesh Data)

• `e` — Finite element edges in the mesh, a subset of all the edges (Mesh Data)

• `u` — Solution of the PDE

• `time` — Time, for parabolic or hyperbolic PDE only

If your boundary conditions do not depend on `u` or `time`, those inputs are `[]`. If your boundary conditions do depend on `u` or `time`, then when `u` or `time` are `NaN`, ensure that the outputs such as `qmatrix` consist of matrices of `NaN` of the correct size. This signals to solvers, such as `parabolic`, to use a time-dependent or solution-dependent algorithm.

Before specifying boundary conditions, you need to know the boundary labels. See Identify Boundary Labels.

The PDE solver, such as `assempde` or `adaptmesh`, passes a matrix `p` of points and `e` of edges. `e` has seven rows and `ne` columns, where you do not necessarily know in advance the size `ne`.

• `p` is a 2-by-Np matrix, where `p(1,k)` is the x-coordinate of point `k`, and `p(2,k)` is the y-coordinate of point `k`.

• `e` is a 7-by-`ne` matrix, where

• `e(1,k)` is the index of the first point of edge `k`.

• `e(2,k)` is the index of the second point of edge `k`.

• `e(5,k)` is the label of the geometry edge of edge `k` (see Identify Boundary Labels).

`e` contains an entry for every finite element edge that lies on an exterior boundary.

Use the following template for your boundary file.

```function [qmatrix,gmatrix,hmatrix,rmatrix] = pdebound(p,e,u,time) ne = size(e,2); % number of edges qmatrix = zeros(1,ne); gmatrix = qmatrix; hmatrix = zeros(1,2*ne); rmatrix = hmatrix; for k = 1:ne x1 = p(1,e(1,k)); % x at first point in segment x2 = p(1,e(2,k)); % x at second point in segment xm = (x1 + x2)/2; % x at segment midpoint y1 = p(2,e(1,k)); % y at first point in segment y2 = p(2,e(2,k)); % y at second point in segment ym = (y1 + y2)/2; % y at segment midpoint switch e(5,k) case {some_edge_labels} % Fill in hmatrix,rmatrix or qmatrix,gmatrix case {another_list_of_edge_labels} % Fill in hmatrix,rmatrix or qmatrix,gmatrix otherwise % Fill in hmatrix,rmatrix or qmatrix,gmatrix end end```

For each column `k` in `e`, entry `k` of `rmatrix` is the value of `rmatrix` at the first point in the edge, and entry `ne` + `k` is the value at the second point in the edge. For example, if r = x2 + y4, then write these lines:

```rmatrix(k) = x1^2 + y1^4; rmatrix(k+ne) = x2^2 + y2^4;```

The syntax for `hmatrix` is identical: entry `k` of `hmatrix` is the value of `r` at the first point in the edge, and entry `k` + `ne` is the value at the second point in the edge.

For each column `k` in `e`, entry `k` of `qmatrix` is the value of `qmatrix` at the midpoint in the edge. For example, if q = x2 + y4, then write these lines:

`qmatrix(k) = xm^2 + ym^4;`

The syntax for `gmatrix` is identical: entry `k` of `gmatrix` is the value of `gmatrix` at the midpoint in the edge.

If the coefficients depend on the solution `u`, use the element `u(e(1,k))` as the solution value at the first point of edge `k`, and `u(e(2,k))` as the solution value at the second point of edge `k`.

For example, consider the following geometry, a rectangle with a circular hole.

```% Rectangle is code 3, 4 sides, % followed by x-coordinates and then y-coordinates R1 = [3,4,-1,1,1,-1,-.4,-.4,.4,.4]'; % Circle is code 1, center (.5,0), radius .2 C1 = [1,.5,0,.2]'; % Pad C1 with zeros to enable concatenation with R1 C1 = [C1;zeros(length(R1)-length(C1),1)]; geom = [R1,C1]; % Names for the two geometric objects ns = (char('R1','C1'))'; % Set formula sf = 'R1-C1'; % Create geometry gd = decsg(geom,sf,ns); % View geometry pdegplot(gd,'EdgeLabels','on') xlim([-1.1 1.1]) axis equal```

Suppose the boundary conditions on the outer boundary (segments 1 through 4) are Dirichlet, with the value u(x,y) = t(x – y), where t is time. Suppose the circular boundary (segments 5 through 8) has a generalized Neumann condition, with q = 1 and g = x2 + y2.

Write the following boundary file to represent the boundary conditions:

```function [qmatrix,gmatrix,hmatrix,rmatrix] = pdebound(p,e,u,time) ne = size(e,2); % number of edges qmatrix = zeros(1,ne); gmatrix = qmatrix; hmatrix = zeros(1,2*ne); rmatrix = hmatrix; for k = 1:ne x1 = p(1,e(1,k)); % x at first point in segment x2 = p(1,e(2,k)); % x at second point in segment xm = (x1 + x2)/2; % x at segment midpoint y1 = p(2,e(1,k)); % y at first point in segment y2 = p(2,e(2,k)); % y at second point in segment ym = (y1 + y2)/2; % y at segment midpoint switch e(5,k) case {1,2,3,4} % rectangle boundaries hmatrix(k) = 1; hmatrix(k+ne) = 1; rmatrix(k) = time*(x1 - y1); rmatrix(k+ne) = time*(x2 - y2); otherwise % same as case {5,6,7,8}, circle boundaries qmatrix(k) = 1; gmatrix(k) = xm^2 + ym^2; end end```

### Boundary Conditions for PDE Systems

The general mixed-boundary conditions for PDE systems of N equations (see Equations You Can Solve Using Legacy Functions) are

`$\begin{array}{l}hu=r\\ n\text{\hspace{0.17em}}·\text{\hspace{0.17em}}\left(c\otimes \nabla u\right)+qu=g+{h}^{\prime }\mu \end{array}$`

The notation $n\text{\hspace{0.17em}}·\text{\hspace{0.17em}}\left(c\otimes \nabla u\right)$ means the N-by-1 matrix with (i,1)-component

`$\sum _{j=1}^{N}\left(\mathrm{cos}\left(\alpha \right){c}_{i,j,1,1}\frac{\partial }{\partial x}+\mathrm{cos}\left(\alpha \right){c}_{i,j,1,2}\frac{\partial }{\partial y}+\mathrm{sin}\left(\alpha \right){c}_{i,j,2,1}\frac{\partial }{\partial x}+\mathrm{sin}\left(\alpha \right){c}_{i,j,2,2}\frac{\partial }{\partial y}\right)\text{\hspace{0.17em}}{u}_{j}$`

where the outward normal vector of the boundary $n=\left(\mathrm{cos}\left(\alpha \right),\mathrm{sin}\left(\alpha \right)\right)$. For each edge segment there are M Dirichlet conditions and the h-matrix is M-by-N, M ≥ 0. The generalized Neumann condition contains a source ${h}^{\prime }\mu$ where the solver computes Lagrange multipliers µ such that the Dirichlet conditions are satisfied.

To write a function file, say `pdebound.m`, use the following syntax:

`[qmatrix,gmatrix,hmatrix,rmatrix] = pdebound(p,e,u,time)`

Your function returns matrices `qmatrix`, `gmatrix`, `hmatrix`, and `rmatrix`, based on these inputs:

• `p` — Points in the mesh (Mesh Data)

• `e` — Finite element edges in the mesh, a subset of all the edges (Mesh Data)

• `u` — Solution of the PDE

• `time` — Time, for parabolic or hyperbolic PDE only

If your boundary conditions do not depend on `u` or `time`, those inputs are `[]`. If your boundary conditions do depend on `u` or `time`, then when `u` or `time` are `NaN`, ensure that the outputs such as `qmatrix` consist of matrices of `NaN` of the correct size. This signals to solvers, such as `parabolic`, to use a time-dependent or solution-dependent algorithm.

Before specifying boundary conditions, you need to know the boundary labels. See Identify Boundary Labels.

A PDE solver, such as `assempde` or `adaptmesh`, passes a matrix `p` of points and `e` of edges. `e` has seven rows and `ne` columns, where you do not necessarily know in advance the size `ne`.

• `p` is a 2-by-Np matrix, where `p(1,k)` is the x-coordinate of point `k`, and `p(2,k)` is the y-coordinate of point `k`.

• `e` is a 7-by-`ne` matrix, where

• `e(1,k)` is the index of the first point of edge `k`.

• `e(2,k)` is the index of the second point of edge `k`.

• `e(5,k)` is the label of the geometry edge of edge `k` (see Identify Boundary Labels).

`e` contains an entry for every finite element edge that lies on an exterior boundary.

Let N be the dimension of the system of PDEs; see Equations You Can Solve Using Legacy Functions. Use the following template for your boundary file.

```function [qmatrix,gmatrix,hmatrix,rmatrix] = pdebound(p,e,u,time) N = 3; % Set N = the number of equations ne = size(e,2); % number of edges qmatrix = zeros(N^2,ne); gmatrix = zeros(N,ne); hmatrix = zeros(N^2,2*ne); rmatrix = zeros(N,2*ne); for k = 1:ne x1 = p(1,e(1,k)); % x at first point in segment x2 = p(1,e(2,k)); % x at second point in segment xm = (x1 + x2)/2; % x at segment midpoint y1 = p(2,e(1,k)); % y at first point in segment y2 = p(2,e(2,k)); % y at second point in segment ym = (y1 + y2)/2; % y at segment midpoint switch e(5,k) case {some_edge_labels} % Fill in hmatrix,rmatrix or qmatrix,gmatrix case {another_list_of_edge_labels} % Fill in hmatrix,rmatrix or qmatrix,gmatrix otherwise % Fill in hmatrix,rmatrix or qmatrix,gmatrix end end```

For the boundary file, you represent the matrix h for each edge segment as a vector, taking the matrix column-wise, as `hmatrix(:)`. Column `k` of `hmatrix` corresponds to the matrix at the first edge point `e(1,k)`, and column `k` + `ne` corresponds to the matrix at the second edge point `e(2,k)`.

Similarly, you represent each vector r for an edge as a column in the matrix `rmatrix`. Column `k` corresponds to the vector at the first edge point `e(1,k)`, and column `k` + `ne` corresponds to the vector at the second edge point `e(2,k)`.

Represent the entries for the matrix q for each edge segment as a vector, `qmatrix(:)`, similar to the matrix `hmatrix(:)`. Similarly, represent g for each edge segment is a column vector in the matrix `gmatrix`. Unlike h and r, which have two columns for each segment, q and g have just one column for each segment, which is the value of the function at the midpoint of the edge segment.

For example, consider the following geometry, a rectangle with a circular hole.

```% Rectangle is code 3, 4 sides, % followed by x-coordinates and then y-coordinates R1 = [3,4,-1,1,1,-1,-.4,-.4,.4,.4]'; % Circle is code 1, center (.5,0), radius .2 C1 = [1,.5,0,.2]'; % Pad C1 with zeros to enable concatenation with R1 C1 = [C1;zeros(length(R1)-length(C1),1)]; geom = [R1,C1]; % Names for the two geometric objects ns = (char('R1','C1'))'; % Set formula sf = 'R1-C1'; % Create geometry gd = decsg(geom,sf,ns); % View geometry pdegplot(gd,'EdgeLabels','on') xlim([-1.1 1.1]) axis equal```

Suppose N = 3. Suppose the boundary conditions are mixed. There is `M` = 1 Dirichlet condition:

• The first component of u = 0 on the rectangular segments (numbers 1–4). So h(1,1) = 1 and r(1) = 0 for those segments.

• The second components of u = 0 on the circular segments (numbers 5–8). So h(2,2) = 1 and r(2) = 0 for those segments.

• On the rectangular segments (numbers 1–4),

`$q=\left(\begin{array}{ccc}0& 1& 1\\ 0& 0& 0\\ 1& 1& 0\end{array}\right)$`

and

`$g=\left(\begin{array}{c}1+{x}^{2}\\ 0\\ 1+{y}^{2}\end{array}\right)$`
• On the circular segments (numbers 5–8),

`$q=\left(\begin{array}{ccc}0& 1+{x}^{2}& 2+{y}^{2}\\ 0& 0& 0\\ 1+{x}^{4}& 1+{y}^{4}& 0\end{array}\right)$`

and

`$g=\left(\begin{array}{c}\mathrm{cos}\left(\pi x\right)\\ 0\\ \mathrm{tanh}\left(x+y\right)\end{array}\right)$`

Write the following boundary file to represent the boundary conditions:

```function [qmatrix,gmatrix,hmatrix,rmatrix] = pdebound(p,e,u,time) N = 3; ne = size(e,2); % number of edges qmatrix = zeros(N^2,ne); gmatrix = zeros(N,ne); hmatrix = zeros(N^2,2*ne); rmatrix = zeros(N,2*ne); for k = 1:ne x1 = p(1,e(1,k)); % x at first point in segment x2 = p(1,e(2,k)); % x at second point in segment xm = (x1 + x2)/2; % x at segment midpoint y1 = p(2,e(1,k)); % y at first point in segment y2 = p(2,e(2,k)); % y at second point in segment ym = (y1 + y2)/2; % y at segment midpoint switch e(5,k) case {1,2,3,4} hk = zeros(N); hk(1,1) = 1; hk = hk(:); hmatrix(:,k) = hk; hmatrix(:,k+ne) = hk; rk = zeros(N,1); % Not strictly necessary rmatrix(:,k) = rk; % These are already 0 rmatrix(:,k+ne) = rk; qk = zeros(N); qk(1,2) = 1; qk(1,3) = 1; qk(3,1) = 1; qk(3,2) = 1; qk = qk(:); qmatrix(:,k) = qk; gk = zeros(N,1); gk(1) = 1+xm^2; gk(3) = 1+ym^2; gmatrix(:,k) = gk; case {5,6,7,8} hk = zeros(N); hk(2,2) = 1; hk = hk(:); hmatrix(:,k) = hk; hmatrix(:,k+ne) = hk; rk = zeros(N,1); % Not strictly necessary rmatrix(:,k) = rk; % These are already 0 rmatrix(:,k+ne) = rk; qk = zeros(N); qk(1,2) = 1+xm^2; qk(1,3) = 2+ym^2; qk(3,1) = 1+xm^4; qk(3,2) = 1+ym^4; qk = qk(:); qmatrix(:,k) = qk; gk = zeros(N,1); gk(1) = cos(pi*xm); gk(3) = tanh(xm*ym); gmatrix(:,k) = gk; end end```

Get trial now