| bcMatrix(N, in_N, qN, gN, in_N0, in_D, hD, rD, in_D0... |
function bcm = bcMatrix(N, in_N, qN, gN, in_N0, in_D, hD, rD, in_D0...
,in_M, qM, gM, hM, rM, in_nb)
% ====== Providing boundary conditions matrix bcm =========
% ======= (system case, e.g. for 2D plane-strain problem, 2 PDEs for 2 dependent variables) =====
%
% A crutial step for solving boundary value problems (BVP) is to define the boundary conditions (BCs).
% In a recent submission,
% "Solve boundary value probelm by command lines," an simple example is
% given for the scalar case (1 PDE for 1 dependent variable)
% subject to Neumann BC only (sub-function "bcNeumann" in that submission).
% I now try to show how to construct the BC
% matrix for system case (Dimession N, i.e. N PDEs for N dependent variables).
% Dirichlet, Neumann, and mixed BCs are considered.
% A more complex example showing how to use this function will be provided
% soon
%
% ===== By Hongxue Cai (h-cai@northwestern.edu) =======
%
%
% ---- Input: ---------->
% N: dimnession of the PDE system (e.g. N = 1 for scalar case, N = 2
% for 2D plane-strain problem)
% in_N: index of non-zero Neumann boundary segments (e.g. if the # 2,
% 6, and 7 geometrical elements are subject to non-zero Neumann BC,
% then in_N = [2 6 7], or [2 7 6] ...;
% in_N = 0 if no non-zero Neumman bc
% qN: q for non-zero Neumann boundary segments, cell array of
% size N^2 x length(in_N) (Note: cell array is convinient as
% q are text exprssions of different length)
% e.g. for 2D plane-strain problem, N = 2, qN{1, 1:length(in_N)} = q11,
% and if in_N = [2 6 7], then qN{1, 1} is q11 for boundary #2,
% qN{1, 2} is q11 for boundary #6; qN{1, 3} is q11 for boundary #7;
% Similarly,
% qN{2, 1:length(in_N)}= q21; qN{3, 1:length(in_N)} = q12; and qN{4, 1:length(in_N)} = q22.
% qN = 0 if no non-zero Neumman bc
% gN: g for non-zero Neumann boundary segments, cell array of
% size N x length(in_N)
% gN = 0 if no non-zero Neumman bc
%
% in_N0: index of zero Neumann boundary segments
% (g1 = g2 = ... = gN = 0)
%
% in_D: index of non-zero Dirichlet boundary segments
% hD: h for non-zero Dirichlet boundary segments, cell array
% rD: r for non-zero Dirichlet boundary segments, cell array
%
% in_D0: index of zero Dirichlet boundary conditions (r1 = r2 = ... = rN = 0)
%
% in_M: index of mixed boundary segments
% qM: q for mixed boundary segments, cell array
% gM: g for mixed boundary segments, cell array
% hM: h for mixed boundary segments, cell array
% rM: r for mixed boundary segments, cell array
%Note: in the case of mixed boundary condition, 1 Dirichlet bc is assumed
%(M =1)
%
% in_nb: index for non-boundary segements, thoses inside the
% geometrical domain who seperate the sub-domains of the geometry
% (an example showing how to use this function and how to use
% different material properties in different sub-domains will be
% given soon
% if non-cell input for q, g, h, r, convert to cell
if isnumeric(qN)
qN = num2cell(qN);
end
if isnumeric(gN)
gN = num2cell(gN);
end
if isnumeric(hD)
hD = num2cell(hD);
end
if isnumeric(rD)
rD = num2cell(rD);
end
if isnumeric(qM)
qM = num2cell(qM);
end
if isnumeric(gM)
gM = num2cell(gM);
end
if isnumeric(hM)
hM = num2cell(hM);
end
if isnumeric(rM)
rM = num2cell(rM);
end
% length of q and g for non-zero Neumann boundaries
for k = 1:length(in_N)
for j = 1:N^2
LqN(j, k) = length(qN{j, k});
end
for j = 1:N
LgN(j, k) = length(gN{j, k});
end
Lsumqg(k) = sum(LqN(:, k)) + sum(LgN(:, k));
end
LN = max(Lsumqg); %maximum rows of text expression
% length of h and r for non-zero Dirichlet boundaries
for k = 1:length(in_D)
for j = 1:N^2
LhD(j, k) = length(hD{j, k});
end
for j = 1:N
LrD(j, k) = length(rD{j, k});
end
Lsumhr(k) = sum(LhD(:, k)) + sum(LrD(:, k));
end
LD = max(Lsumhr); %maximum rows of text expression
% length of q, g, h, and r for mixed boundaries
if in_M ~= 0
for k = 1:length(in_M)
for j = 1:N^2
LqM(j, k) = length(qM{j, k});
% LhM(j, k) = length(hM{j, k});
end
for j = 1:N
LgM(j, k) = length(gM{j, k});
LhM(j, k) = length(hM{j, k});
end
LrM(1, k) = length(rM{1, k}); %for mixed bc, 1 Dirichlet bc assumed
Lsumqghr(k) = sum(LqM(:, k)) + sum(LgM(:, k)) + sum(LhM(:, k)) + sum(LrM(:, k));
end
LM = max(Lsumqghr); %maximum rows of text expression
else
LM = 0;
end
L_N = length(in_N);
L_N0 = length(in_N0);
L_D = length(in_D);
L_D0 = length(in_D0);
L_M = length(in_M);
Lnb = length(in_nb);
%number of column of the bc matrix
coLnbcm = L_N + L_N0 + L_D + L_D0 + L_M + Lnb;
% row # before text expression begins
row0_N = 3 + N^2 + N + - 1; % non-zero Neumann boundaries
M = N;
row0_D = 3 + M - 1; % non-zero Dirichlet boundaries
M = 1; % 1 Dirichlet bc is assumed
row0_M = 3 + N^2 + N + M*N + M - 1; % mixed boundaries
% number of rows of the bc matrix
row_bcm = max([LN+row0_N, LD+row0_D, LM+row0_M]);
bcm = zeros(row_bcm, coLnbcm);
bcm(1,:) = N; % row#1 --> Dimension of the system
% ------------------ Non-zero Neumman boundaries ---------------------->
if in_N ~= 0
for ij = 1:L_N
bcm(2, in_N(ij)) = 0; % M = 0 -> non Dirichlet bc
bcm(3:3+N^2-1,in_N(ij)) = LqN(:, ij); %Length for q11, q21, ..., qNN
bcm(3+N^2:3+N^2+N-1,in_N(ij)) = LgN(ij); %Length for g1, g2, ..., gN
%Note: as M = 0, no need for Lh, Lr, h, r
%
%q for non-zero Nmann boundary segments
nRowN = 3 + N^2 + N;
for j = 1:N^2
bcm(nRowN:nRowN+LqN(j, ij)-1, in_N(ij)) = qN{j,ij};
nRowN = nRowN + LqN(j, ij);
end
%g for non-zero Nmann boundary segments
for j = 1:N
bcm(nRowN:nRowN+LgN(j, ij)-1, in_N(ij)) = gN{j,ij};
nRowN = nRowN + LgN(j, ij);
end
end
end
% <------------------ end of Non-zero Neumman boundaries --------
%
%
% ------------------ Zero Neumman boundaries ---------------------->
if in_N0 ~= 0
bcm(2, in_N0) = 0; % M = 0 -> non Dirichlet bc
bcm(3:3+N^2-1,in_N0) = 1; %Length for q11, q21, ..., qNN
bcm(3+N^2:3+N^2+N-1,in_N0) = 1; %Length for g1, g2, ..., gN
%Note: as M = 0, no need for Lh, Lr, h, r
%
%q
nRowN = 3 + N^2 + N;
bcm(nRowN:nRowN+N^2-1, in_N0) = '0';
nRowN = nRowN + N^2;
%g for non-zero Nmann boundary segments
bcm(nRowN:nRowN+N-1, in_N0) = '0';
nRowN = nRowN + N;
end
% <------------------ Zero Neumman boundaries --------
%
%
%
% ------------------- Non-zero Dirichlet boundaries --------------->
if in_D ~= 0
M = N;
for ij = 1:L_D
bcm(2, in_D(ij)) = M; %M = N -> Dirichlet bc only
bcm(3:3+N^2-1, in_D(ij)) = 1; %Length for q11, q21, ..., qNN
bcm(3+N^2:3+N^2+N-1, in_D(ij)) = 1; %Length for g1, g2, ..., gN
bcm(3+N^2+N:3+N^2+N+M*N-1, in_D(ij)) = LhD(:,ij); %Length for h11, h21, ..., hNN
bcm(3+N^2+N+M*N:3+N^2+N+M*N+M-1, in_D(ij)) = LrD(:,ij); %Length for r1, r2, ..., rN
%
%q = 0
nRowD = 3+N^2+N+M*N+M;
for j = 1:N^2
bcm(nRowD:nRowD+N^2-1,in_D(ij)) = '0';
end
nRowD = nRowD + N^2;
%g = 0
for j = 1:N
bcm(nRowD:nRowD+N-1,in_D(ij)) = '0';
end
nRowD = nRowD + N;
%h
for j = 1:M^2
bcm(nRowD:nRowD+LhD(j, ij)-1, in_D(ij)) = hD{j,ij};
nRowD = nRowD + LhD(j, ij);
end
%r
for j = 1:M
bcm(nRowD:nRowD+LrD(j, ij)-1, in_D(ij)) = rD{j,ij};
nRowD = nRowD + LrD(j, ij);
end
end
end
% <------------------- end of Non-zero Dirichlet boundaries
%
%
% ------- Zero Dirichlet boundaries (rigid walls) ----->
if in_D0 ~= 0
M = N;
bcm(2, in_D0) = M; %M = N -> Dirichlet bc only
bcm(3:3+N^2-1, in_D0) = 1; %Length for q11, q21, ..., qNN
bcm(3+N^2:3+N^2+N-1, in_D0) = 1; %Length for g1, g2, ..., gN
bcm(3+N^2+N:3+N^2+N+M*N-1, in_D0) = 1; %Length for h11, h21, ..., hNN
bcm(3+N^2+N+M*N:3+N^2+N+M*N+M-1, in_D0) = 1; %Length for r1, r2, ..., rN
%
%qij = 0
nRowD = 3+N^2+N+M*N+M;
for j = 1:N^2
bcm(nRowD:nRowD+N^2-1,in_D0) = '0';
end
nRowD = nRowD + N^2;
%g = 0
for j = 1:N
bcm(nRowD:nRowD+N-1,in_D0) = '0';
end
nRowD = nRowD + N;
%hii = 1, hij = 0
bcm(nRowD:nRowD+M^2-1, in_D0) = '0';
bcm(nRowD:M+1:nRowD+M^2-1, in_D0) = '1';
nRowD = nRowD + M^2;
%r = 0
bcm(nRowD:nRowD+M-1, in_D0) = '0';
end
% <------- end of Zero-Dirichlet boundaries (rigid walls) -----
%
%
% ------------------ Mixed boundaries -------------------- ----->
if in_M ~= 0
M =1;
for ij = 1:L_M
bcm(2, in_M(ij)) = M; %M = 1 -> 1 Dirichlet bc
bcm(3:3+N^2-1, in_M(ij)) = LqM(:,ij);; %Length for q11, q21, ..., qNN
bcm(3+N^2:3+N^2+N-1, in_M(ij)) = LgM(:,ij);; %Length for g1, g2, ..., gN
bcm(3+N^2+N:3+N^2+N+M*N-1, in_M(ij)) = LhM(:,ij); %Length for h11, h12, ..., h1N
bcm(3+N^2+N+M*N:3+N^2+N+M*N+M-1, in_M(ij)) = LrM(:,ij); %Length for r1, r2, ..., rN
%
%q
nRowM = 3+N^2+N+M*N+M;
for j = 1:N^2
bcm(nRowM:nRowM+LqM(j, ij)-1, in_M(ij)) = qM{j,ij};
nRowM = nRowM + LqM(j, ij);
end
%g
for j = 1:N
bcm(nRowM:nRowM+LgM(j, ij)-1, in_M(ij)) = gM{j,ij};
nRowM = nRowM + LgM(j, ij);
end
%h
for j = 1:N %1 Dirichlet BC assumed
bcm(nRowM:nRowM+LhM(j, ij)-1, in_M(ij)) = hM{j,ij};
nRowM = nRowM + LhM(j, ij);
end
%r
j = 1; %1, not N, 1 Dirichlet BC assumed
bcm(nRowM:nRowM+LrM(j, ij)-1, in_M(ij)) = rM{j,ij};
nRowM = nRowM + LrM(j, ij);
end
end
% <------------------ end of mixed boundaries --------
%
%
% ---------- non-boundary segments, sub-domain segments ----->
if in_nb ~= 0
for ij = 1:Lnb
bcm(1,in_nb(ij)) = 0; %system number N = 0 for border segments
bcm(2,in_nb(ij)) = N; %Dirichlet boundary number M = N
end
end
% <---------- end of non-boundary segments, sub-domain segments ---- |
|
