This function assumes that the unknown u, which represents the discrete
solution the continuous second-order 1D PDE operator
L U = f,
with continuous boundary condition
a0 U + b0 dU/dn = g,
are given at the 1D cell centers and vertices. Furthermore, all discrete
calculations are performed at the 1D cell centers and vertices.
The function receives as input quantities associated to the discrete
analog of the continuous problem given by the squared linear system
A u = b
where A is the discrete analog of L and b is the discrete analog of g,
both constructed by the user without boundary conditions.
The function output is the modified square linear system
A u = b
where both A and b include boundary condition information.
The boundary condition is always one of the following forms:
For Dirichlet set: a0 not equal zero and b0 = 0.
For Neumann set : a0 = 0 and b0 not equal zero.
For Robin set : both a0 and b0 not equal zero.
For Periodic set : both a0 = 0 and b0 = 0.
For periodic bc, it is assumed that not only u but also du/dn are the same
in both extremes of the domain since a second-order PDE is assumed.
The code assumes the following assertions:
assert(k >= 2, 'k >= 2');
assert(mod(k, 2) == 0, 'k % 2 = 0');
assert(m >= 2*k+1, ['m >= ' num2str(2*k+1) ' for k = ' num2str(k)]);
Parameters:
output
A : Linear operator with boundary conditions added
b : Right hand side with boundary conditions added
input
A : Linear operator without boundary conditions added
b : Right hand side without boundary conditions added
k : Order of accuracy
m : Number of cells
dx : Step size
dc : a0 (2x1 vector for left and right vertices, resp.)
nc : b0 (2x1 vector for left and right vertices, resp.)
v : g (2x1 vector for left and right vertices, resp.)
----------------------------------------------------------------------------
SPDX-License-Identifier: GPL-3.0-or-later
© 2008-2024 San Diego State University Research Foundation (SDSURF).
See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details.
----------------------------------------------------------------------------
0001 function [A, b] = addBC1D(A, b, k, m, dx, dc, nc, v) 0002 % This function assumes that the unknown u, which represents the discrete 0003 % solution the continuous second-order 1D PDE operator 0004 % L U = f, 0005 % with continuous boundary condition 0006 % a0 U + b0 dU/dn = g, 0007 % are given at the 1D cell centers and vertices. Furthermore, all discrete 0008 % calculations are performed at the 1D cell centers and vertices. 0009 % 0010 % The function receives as input quantities associated to the discrete 0011 % analog of the continuous problem given by the squared linear system 0012 % A u = b 0013 % where A is the discrete analog of L and b is the discrete analog of g, 0014 % both constructed by the user without boundary conditions. 0015 % The function output is the modified square linear system 0016 % A u = b 0017 % where both A and b include boundary condition information. 0018 % 0019 % The boundary condition is always one of the following forms: 0020 % 0021 % For Dirichlet set: a0 not equal zero and b0 = 0. 0022 % For Neumann set : a0 = 0 and b0 not equal zero. 0023 % For Robin set : both a0 and b0 not equal zero. 0024 % For Periodic set : both a0 = 0 and b0 = 0. 0025 % 0026 % For periodic bc, it is assumed that not only u but also du/dn are the same 0027 % in both extremes of the domain since a second-order PDE is assumed. 0028 % 0029 % The code assumes the following assertions: 0030 % assert(k >= 2, 'k >= 2'); 0031 % assert(mod(k, 2) == 0, 'k % 2 = 0'); 0032 % assert(m >= 2*k+1, ['m >= ' num2str(2*k+1) ' for k = ' num2str(k)]); 0033 % 0034 % Parameters: 0035 % output 0036 % A : Linear operator with boundary conditions added 0037 % b : Right hand side with boundary conditions added 0038 % 0039 % input 0040 % A : Linear operator without boundary conditions added 0041 % b : Right hand side without boundary conditions added 0042 % k : Order of accuracy 0043 % m : Number of cells 0044 % dx : Step size 0045 % dc : a0 (2x1 vector for left and right vertices, resp.) 0046 % nc : b0 (2x1 vector for left and right vertices, resp.) 0047 % v : g (2x1 vector for left and right vertices, resp.) 0048 % ---------------------------------------------------------------------------- 0049 % SPDX-License-Identifier: GPL-3.0-or-later 0050 % © 2008-2024 San Diego State University Research Foundation (SDSURF). 0051 % See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details. 0052 % ---------------------------------------------------------------------------- 0053 0054 % verify bc sizes and square linear system 0055 assert(all(size(dc) == [2 1]), 'dc is a 2x1 vector'); 0056 assert(all(size(nc) == [2 1]), 'nc is a 2x1 vector'); 0057 assert(all(size(v) == [2 1]), 'v is a 2x1 vector'); 0058 assert(size(A,1) == size(A,2), 'A is a square matrix'); 0059 assert(size(A,2) == numel(b), 'b size = A columns'); 0060 0061 % remove first and last rows of A 0062 vec = sparse(2,1); 0063 vec(1) = 1; 0064 vec(2) = size(A,1); 0065 0066 [rows,cols,s] = find(A(vec,:)); 0067 A = A - sparse(vec(rows), cols, s, size(A,1), size(A,2)); 0068 0069 % remove first and last coefficients of right-hand-side vector b 0070 b(vec) = 0; 0071 0072 [Abcl,Abcr] = addBC1Dlhs(k, m, dx, dc, nc); 0073 A = A + Abcl + Abcr; 0074 b = addBC1Drhs(b, dc, nc, v, vec); 0075 end