No BSD License  

bcMatrix, boundary condition Matrix

by

 

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

Contact us