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

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

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 $$\overrightarrow{n}\text{\hspace{0.17em}}\xb7\text{\hspace{0.17em}}\left(c\nabla u\right)+qu=g$$ on ∂Ω.

$$\overrightarrow{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-*N*matrix, where_{p}`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* = *x*^{2} + *y*^{4},
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* = *x*^{2} + *y*^{4},
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.

Code for generating the figure

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* = *x*^{2} + *y*^{2}.

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

The general mixed-boundary conditions for PDE systems of *N* equations
(see Systems of PDEs) are

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

The notation $$n\text{\hspace{0.17em}}\xb7\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}(\alpha ){c}_{i,j,1,1}\frac{\partial}{\partial x}+\mathrm{cos}(\alpha ){c}_{i,j,1,2}\frac{\partial}{\partial y}+\mathrm{sin}(\alpha ){c}_{i,j,2,1}\frac{\partial}{\partial x}+\mathrm{sin}(\alpha ){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}(\alpha ),\mathrm{sin}(\alpha )\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-*N*matrix, where_{p}`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 Systems of PDEs. 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.

Code for generating the figure

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}(\pi x)\\ 0\\ \mathrm{tanh}(x+y)\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

Was this topic helpful?