boundarymatrix

 b = boundarymatrix([arglist]) by Pruefert
 defines a 1D boundary matrix 
 usuage:
 same BC on both boundary:
 b = boundarymatrix            homogenious Neumann BCs  
 b = boundarymatrix(g)         Neumann BCs
 b = boundarymatrix(q,g)       Robin BCs 
 b = boundarymatrix([],[],h,r) Dirichlet BC

 different BC: e.g. Dirichlet on the left bound, Robin on the right
 boundary
 b = boundarymatrix([],[],h,r,q,g,[],[])
 q,g,h,r are STRINGS like '1' 'sin(s)' etc... but NO inlines
(c) Uwe.Pruefert@tu-berlin.de

0001 function b = boundarymatrix(q,g,h,r,varargin)
0002 % b = boundarymatrix([arglist]) by Pruefert
0003 % defines a 1D boundary matrix
0004 % usuage:
0005 % same BC on both boundary:
0006 % b = boundarymatrix            homogenious Neumann BCs
0007 % b = boundarymatrix(g)         Neumann BCs
0008 % b = boundarymatrix(q,g)       Robin BCs
0009 % b = boundarymatrix([],[],h,r) Dirichlet BC
0010 %
0011 % different BC: e.g. Dirichlet on the left bound, Robin on the right
0012 % boundary
0013 % b = boundarymatrix([],[],h,r,q,g,[],[])
0014 % q,g,h,r are STRINGS like '1' 'sin(s)' etc... but NO inlines
0015 %(c) Uwe.Pruefert@tu-berlin.de
0016 
0017 
0018 switch nargin                  % handle the "minimalistic syntax" option
0019     case {0}
0020         q = '0';
0021         g = '0';
0022         h = [];
0023         r = [];
0024     case {1}  
0025         g =  q;
0026         q = '0';
0027         h = [];
0028         r = [];
0029     case {2}
0030         h = [];
0031         r = [];
0032     otherwise
0033         %
0034 end    
0035 optargs = length(varargin);
0036 nobdrs = max(1,1+optargs/4);    % number of definition blocks, if only one,
0037                                 % we assume the same BCs on both boundaries
0038 
0039 b = cell(nobdrs);               % initialize b as cell array, later we
0040                                 % will overwrite it
0041 if ~mod(optargs,4)==0,
0042     error('argument list has wrong lenght)')
0043 end
0044    
0045 if isempty(h)&&isempty(r),
0046    b1 = 1;                      % only Neumann BCs
0047    b1(2,1) = 0;
0048    b1(3,1) = length(q);
0049    b1(4,1) = length(g);
0050    b1=[b1;double(q)';double(g)'];
0051 elseif isempty(q)&&isempty(g)   % Dirichlet BCs, setting the Neumann part
0052    b1 = 1;                      % to zero
0053    b1(2,1) = 1;
0054         b1(3,1) = 1;
0055         b1(4,1) = 1;
0056         b1(5,1) = length(h);
0057         b1(6,1) = length(r);
0058         b1(7,1) = 48;
0059         b1(8,1) = 48;
0060         b1 = [b1;double(h)';double(r)'];    
0061 end
0062 b{1} = b1;
0063 
0064 % more than one block of BCs in the parameter list
0065 l = 1;
0066 
0067 for k = 2:nobdrs,
0068     q = varargin{l}; l = l+1;
0069     g = varargin{l}; l = l+1;
0070     h = varargin{l}; l = l+1;
0071     r = varargin{l}; l = l+1;
0072     if isempty(h)&&isempty(r),
0073         b1 = 1; %
0074         b1(2,1) = 0;
0075         b1(3,1) = length(q);
0076         b1(4,1) = length(g);
0077         b1 = [b1;double(q)';double(g)']; 
0078         b{k} = b1;
0079     elseif isempty(q)&&isempty(g) %
0080         b1 = 1;
0081         b1(2,1) = 1;
0082         b1(3,1) = 1;
0083         b1(4,1) = 1;
0084         b1(5,1) = length(h);
0085         b1(6,1) = length(r);
0086         b1(7,1) = 48;
0087         b1(8,1) = 48;
0088         b1 = [b1;double(h)';double(r)'];
0089         b{k} = b1;
0090     end
0091 end
0092 
0093                                 % contruct the boundary matrix, we have to
0094                                 % care about the right dimension of ALL
0095                                 % rows when assembling  the matrix
0096                                 
0097 if nobdrs==1,                   % the case of one BCs for all boundaries
0098     bb=[b{1},b{1}];
0099 else
0100     maxlengthb = 0;
0101     for l = 1:k,                % compute the max lenght of all entries of b
0102         if length(b{l})>maxlengthb,
0103             maxlengthb = length(b{l});
0104         end
0105     end
0106     bb = zeros(maxlengthb,nobdrs);  % initalize the matrix
0107     for l = 1:k,
0108         bb(1:length(b{l}),l)=b{l};
0109     end
0110 end
0111 b=bb;                           % return b
0112 end