function B=BMatrix (nx,ny,nz,varargin)

% Macierz B(:,:) pomnozona przez wektor n = nx*ny*nz wartosci funkcji 
% na siatce 3D zwraca wektor (nx+2)*(ny+2)*(nz+2) wartosci na 
% wiekszej siatce. Dolozone sa wartosci brzegowe rowne zero (gdy
% na wejsciu sa tylko 3 argumenty nx, ny, nz) lub wartosci brzegowe
% bedace powtorzeniem sasiadow - gdy sa wiecej niz 3 argumenty, np. 
% (nx, ny, nz, 'D').

if length(varargin)==0 
    zero_bd = 1;
else
    zero_bd = 0;
end;

n = nx*ny*nz; % tyle jest punktow siatki
M = (nx+2)*(ny+2)*(nz+2); % tyle jest punktow powiekszonej siatki

B = zeros(M,n); 

if zero_bd
    for i=1:M % petla po wezlach wiekszej siatki
        [xk,yk,zk] = wspolrzedne(nx+2,ny+2,nz+2,i); % wspolrzedne punktu w wiekszej siatce
        % Dla punktu na brzegu - zostawiamy zero
        if (xk==1)||(yk==1)||(zk==1)||(xk==(nx+2))||(yk==(ny+2))||(zk==(nz+2))
        else
            nr = nr_punktu(nx,ny,xk-1,yk-1,zk-1); % numer odp. punktu w mniejszej siatce
            B(i,nr) = 1;
        end;
    end;
else
    for i=1:M % petla po wezlach wiekszej siatki
        [xk,yk,zk] = wspolrzedne(nx+2,ny+2,nz+2,i); % wspolrzedne punktu w wiekszej siatce
        % Dla punktu na brzegu - zostawiamy zero
        if (xk==1)||(yk==1)||(zk==1)||(xk==(nx+2))||(yk==(ny+2))||(zk==(nz+2))
            if (xk==1)||(xk==(nx+2))
                if yk==1
                    yk=2;
                end;
                if yk==(ny+2)
                    yk=ny+1;
                end;
                if zk==1
                    zk=2;
                end;
                if zk==(nz+2)
                    zk=nz+1;
                end;
                if (xk==1)
                    nr = nr_punktu(nx,ny,1,yk-1,zk-1); % numer odp. punktu w mniejszej siatce
                else
                    nr = nr_punktu(nx,ny,nx,yk-1,zk-1); 
                end;
                B(i,nr) = 1;
            elseif (yk==1)||(yk==(ny+2))
                if zk==1
                    zk=2;
                end;
                if zk==(nz+2)
                    zk=nz+1;
                end;
                if (yk==1)
                    nr = nr_punktu(nx,ny,xk-1,1,zk-1); % numer odp. punktu w mniejszej siatce
                else
                    nr = nr_punktu(nx,ny,xk-1,ny,zk-1); 
                end;
                B(i,nr) = 1;
            else
                if (zk==1)
                    nr = nr_punktu(nx,ny,xk-1,yk-1,1); % numer odp. punktu w mniejszej siatce
                else
                    nr = nr_punktu(nx,ny,xk-1,yk-1,nz); 
                end;
                B(i,nr) = 1;
            end;
        else
            nr = nr_punktu(nx,ny,xk-1,yk-1,zk-1); % numer odp. punktu w mniejszej siatce
            B(i,nr) = 1;
        end;
    end;
end;

