# bcMatrix, boundary condition Matrix

### Hongxue Cai (view profile)

15 Feb 2006 (Updated )

Provide boundary condition matrix for 2D mechncal finite element modeling

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  ----```