Documentation

Boundary Conditions by Writing Functions

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

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 n·(cu)+qu=g on ∂Ω.

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.

 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 = 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 Systems of PDEs) are

hu=rn·(cu)+qu=g+hμ.

The notation n·(cu) means the N-by-1 matrix with (i,1)-component

j=1N(cos(α)ci,j,1,1x+cos(α)ci,j,1,2y+sin(α)ci,j,2,1x+sin(α)ci,j,2,2y)uj,

where the outward normal vector of the boundary n=(cos(α),sin(α)). 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μ 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 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=(011000110)

    and

    g=(1+x201+y2)

  • On the circular segments (numbers 5–8),

    q=(01+x22+y20001+x41+y40)

    and

    g=(cos(πx)0tanh(x+y))

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

Related Examples

Was this topic helpful?