function md2=meshadaptation(md,field,scale,epsilon)
%MESHADAPTATION - remesh a model for a given field
%
%   remesh a model for a given field, use the hessian computation to spread the error on the element
%   scale: 1.2 large scale leads to higher refining.
%   epsilon: .5  large epsilon will refine fewer elements
%
%   Usage:
%      md2=meshadaptation(md,field,scale,epsilon)
%
%   Example:
%      md2=meshadaptation(md,md.vel,1.2,0.5);

%some checks
if ~strcmpi(md.type,'2d')
	error('meshadaptation error message: adaptation for 3d meshes not implemented yet')
end
if length(field)~=md.numberofgrids
	error('meshadaptation error message: input field length shoud be numberofgrids')
end

disp(sprintf('      metric computation') )

%load some variables (it is much faster if the variab;es are loaded from md once for all) 
numberofelements=md.numberofelements;
numberofgrids=md.numberofgrids;
index=md.elements;
x=md.x; y=md.y; z=md.z;

%initialization
alpha=zeros(md.numberofelements,3);
beta=zeros(md.numberofelements,3);
gradx=zeros(md.numberofgrids,1);
grady=zeros(md.numberofgrids,1);
metric=zeros(md.numberofelements,1);

%build some usefull variables
field=field(:);
line=index(:);
summation=1/3*ones(3,1);
linesize=3*numberofelements;
x1=x(index(:,1)); x2=x(index(:,2)); x3=x(index(:,3)); y1=y(index(:,1)); y2=y(index(:,2)); y3=y(index(:,3));

%compute nodal functions coefficients N(x,y)=alpha x + beta y + gamma
invdet=1./(x1.*(y2-y3)-x2.*(y1-y3)+x3.*(y1-y2));
alpha=[invdet.*(y2-y3) invdet.*(y3-y1) invdet.*(y1-y2)];
beta=[invdet.*(x3-x2) invdet.*(x1-x3) invdet.*(x2-x1)];

%Compute gradient for each element
grad_elx=sum(field(index).*alpha,2); 
grad_ely=sum(field(index).*beta,2);

%compute the volume of each element
areas=area(md);

%update weights that holds the volume of all the element holding the grid i
weights=sparse(line,ones(linesize,1),repmat(areas,3,1),numberofgrids,1);

%Compute gradient for each grid (average of the elements around)
gradx=sparse(line,ones(linesize,1),repmat(areas.*grad_elx,3,1),numberofgrids,1);
grady=sparse(line,ones(linesize,1),repmat(areas.*grad_ely,3,1),numberofgrids,1);
gradx=gradx./weights;
grady=grady./weights;

%Compute hessian for each element
hessian=[sum(gradx(index).*alpha,2) sum(grady(index).*beta,2) sum(grady(index).*alpha,2)];

%Compute the new metric
hmin=min(sqrt(areas))/scale;
hmax=max(sqrt(areas));

%first, find the eigen values of eah line of H=[hessian(i,1) hessian(i,3); hessian(i,3)  hessian(i,2)]
a=hessian(:,1); b=hessian(:,3); d=hessian(:,2);
lambda1=0.5*((a+d)+sqrt(4*b.^2+(a-d).^2));
lambda2=0.5*((a+d)-sqrt(4*b.^2+(a-d).^2));

%Modify the eigen values to control the shape of the elements
lambda1=min(max(abs(lambda1)/epsilon,1/hmax^2),1/hmin^2);
lambda2=min(max(abs(lambda2)/epsilon,1/hmax^2),1/hmin^2);

%find the corresponding eigen vectors of the initial eigen values
%NOT USED FOR NOW
%div1=sqrt(1+(b./(lambda1-d)).^2);
%div2=sqrt(1+(b./(lambda2-d)).^2);
%vector1=[1./div1 (b./(lambda1-d))./div1];
%vector2=[1./div2 (b./(lambda2-d))./div2];

%compute the metric
%R=[vector1(:,1) vector2(:,1) vector1(:,2) vector2(:,2)];
%det=vector1(:,1).*vector2(:,2)-vector1(:,2).*vector1(:,2);
%Rinv=[vector2(:,2)./det   -vector2(:,1)./det   -vector1(:,2)./det   vector1(:,1)./det];
metric=20*1./sqrt(lambda1.^2+lambda2.^2);

%ok, metric can't go under 1km resolution.
metric(find(metric<10^6))=10^6;

%Remesh with this new metric
disp(sprintf('      initial number of element: %i', md.numberofelements))
md2=meshrefine(md,full(metric));
disp(sprintf('      new number of elements:    %i', md2.numberofelements))
