function [imax jmax bc nbn K] = geobc(fn);
fid = fopen(fn,'r');
bc = [];
xx = fgetl(fid);
ij = str2num(xx); % first line
imax = ij(1); jmax = ij(2);
K = imax*jmax; % total number of nodes

xx = fgetl(fid);
while xx ~= -1
    bc = [bc; str2num(xx)];
    xx = fgetl(fid);
end
[nbn J] = size(bc);
if J == 3 % bc in the form of (i, j, v)
    bc = [(bc(:,1)-1)*jmax+bc(:,2) bc(:,3)]; % bc changed into (k, v)
end
return;