function [X,Y,well]=makegrid(well,aroundAll,aroundEach,dmax,dmin,factor)
% [X,Y,well]=makegrid(well,aroundAll,aroundEach,dmax,dmin,factor)
% TO simplified more radical grid generation
% There is an overall grid with uniform cell size covering all wells and
% a rectangle around the well set (distance aroundAll)
% Around each well there is a fine  grid starting with cell size dmin
% increasing to the cell size of the overall grid by a given growth factor
% these meshes are merged and columns and rows smaller than dmin are
% set to the average of the two closeby grid lines, eliminating one.
% TO 091117 091124

fprintf('Generating computation FDM grid for this model ....\n');

NWEL=length(well);

if aroundAll<aroundEach
    error('aroundAll<aroundEach');
end
if dmax<=dmin
    error('dmax<=dmin');
end
if factor<1
    eror('factor<1');
end

%% round well coordinates to a spacing equal to dmin
for iW=1:NWEL
    well(iW).x=round(well(iW).x/dmin)*dmin;
    well(iW).y=round(well(iW).y/dmin)*dmin;
end

%% initial coarse overall grid
X=(min([well.x])-aroundAll):dmax:(max([well.x])+aroundAll);
Y=(min([well.y])-aroundAll):dmax:(max([well.y])+aroundAll);

%% subgrid to be placed around each well
d=0.5*dmin; dd=dmin;
while sum(d)<aroundEach && dd<dmax
    dd=factor*dd;
    d=[d d(end)+dd];
end
d=cumsum(d);
subgrid=[-d(end:-1:1)  d];

%% merging subgrids with overall grid
for iW=1:NWEL
    X=[X well(iW).x+subgrid];
    Y=[Y well(iW).y+subgrid];
end

X=unique(X); X=cleangrid(X,dmin,'COL');
Y=unique(Y); Y=cleangrid(Y,dmin,'ROW'); Y=Y(end:-1:1);  % let Y always run backward

%% find grid locations of wells

for iW=1:NWEL
    well(iW).ix=find(well(iW).x>X,1,'last');
    well(iW).iy=find(well(iW).y<Y,1,'last');
end 

fprintf('.... Grid construction finished.\n');

function U=cleangrid(U,dmin,dir)
% U=cleangrid(U,dmin,dir)
% remove columns smaller than dmin

N=length(U); dU=abs(diff(U));

i=1;
while i<length(dU);
    if i<length(dU) && dU(i)<dmin
        dU(i+1)=dU(i)+dU(i+1);
        dU(i)=[];
    else
        i=i+1;
    end
end
U=U(1)+cumsum(dU);

fprintf('N%s reduced from %d to %d by removing %ss smaller than %g from initial grid.\n',dir,N-1,length(U)-1,lower(dir),dmin);

