function [q,g,h,r]=geomb(p,e,u,time)
%Boundary condition data
fprintf('Function geomb being evaluated \n');
b=[
 1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   %N
 1   1   1   0   0   0   0   0   1   1   1   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   %M
 1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   %Q
 1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   %G
 1   1   1   48  48  48  48  48  1   1   1   48  48  48  48  48  48  48  48  48  48  48  48  48  48  48  48  48  48  48  %H
 1   1   1   48  48  48  48  48  1   1   1   48  48  48  48  48  48  48  48  48  48  48  48  48  48  48  48  48  48  48  %R
 48  48  48  0   0   0   0   0   48  48  48  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   %Q
 48  48  48  0   0   0   0   0   48  48  48  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   %G
 49  49  49  0   0   0   0   0   49  49  49  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   %H
 48  48  48  0   0   0   0   0   48  48  48  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   %R 
];

if any(size(u))
  [q,g,h,r]=pdeexpdi(p,e,u,time,b);
else
  [q,g,h,r]=pdeexpdi(p,e,time,b);
end
