function [optmesh,bndmesh,vert,L,M] = convomesh(n,a1,b1,a2,b2,lev)

% computes an optimal polynomial mesh on a planar C^2 convex domain,
% closure of the interior of a level curve of a bivariate convex function 
% f defined below by the user; the convex domain is considered as a 
% starlike domain with respect to its barycenter 

% by Federico Piazzon and Marco Vianello, University of Padova
% December 2014


% INPUT: 

% n = polynomial degree 

% [a1,b1] x [a2,b2] = rectangle containing the level curve of f 

% lev = level of f 


% OUTPUT: 

% optmesh = 2-column array of the mesh points

% bndmesh = boundary mesh

% vert = polygon approximating the boundary

% L = approximate length of the boundary

% M = approximate maximal curvature of the boundary



% FUNCTION BODY

NN=10; tol=10^(-4); err=tol+1; alpha=0.5;

% boundary length by polygonal approximation
Lold=0;
while (err > tol) 
NN=2*NN;
x=linspace(a1,b1,NN);
stepx=(b1-a1)/(NN-1);
y=linspace(a2,b2,NN);
stepy=(b2-a2)/(NN-1);
[u,v]=meshgrid(x,y);
z=f(u,v);
levc=contourc(x,y,z,[lev lev]);
pts=[levc(1,2:end);levc(2,2:end)]; 
k=convhull(pts(1,:),pts(2,:));
vert=[pts(1,k);pts(2,k)]';
diff=vert(2:end,:)-vert(1:end-1,:);
delta=sqrt(diff(:,1).^2+diff(:,2).^2);
Lnew=sum(delta);
err=abs(Lnew-Lold);
Lold=Lnew;
end;

% L is the approximate length
L=Lnew;

% computing curvatures at the polygon vertices by the formula 
% curvature = (fxx fy^2 - 2 fxy fx fy + fyy fx^2)/(fx^2+fy^2)^(3/2)

ff=f(vert(:,1),vert(:,2));

% 0=-step,1=there,2=+step

ff21=f(vert(:,1)+stepx,vert(:,2));
ff01=f(vert(:,1)-stepx,vert(:,2));
fx=(ff21-ff01)/(2*stepx);
fxx=(ff21-2*ff+ff01)/stepx^2;

ff22=f(vert(:,1)+stepx,vert(:,2)+stepy);
ff02=f(vert(:,1)-stepx,vert(:,2)+stepy);
fx2=(ff22-ff02)/(2*stepx);
ff20=f(vert(:,1)+stepx,vert(:,2)-stepy);
ff00=f(vert(:,1)-stepx,vert(:,2)-stepy);
fx0=(ff20-ff00)/(2*stepx);
fxy=(fx2-fx0)/(2*stepy);

ff12=f(vert(:,1),vert(:,2)+stepy);
ff10=f(vert(:,1),vert(:,2)-stepy);
fy=(ff12-ff10)/(2*stepy);
fyy=(ff12-2*ff+ff10)/stepy^2;

curv=abs(fxx.*fy.^2-2*fxy.*fx.*fy+fyy.*fx.^2)./(fx.^2+fy.^2).^1.5;

% computing the boundary mesh 
M=max(curv);

N=ceil(M*n*L/(2*alpha)); 
deltas=L/N;go=true;j=1;k=1;
while (go==true) 
partL=0;
while (partL<deltas & go==true)
partL=partL+delta(j);
j=j+1;if j>length(delta) go=false; end;
end;
if partL>deltas 
j=j-1;
end;
bndmesh(k,:)=vert(j,:);k=k+1;
end;

% radial Chebyshev-Lobatto points
ind=(0:2*n);
xi=cos(ind*pi/(2*n));

% computing the optimal mesh and plotting 

% computing the barycenter
c=sum(vert)/length(vert(:,1));
c1=[c(1)*ones(length(vert),1) c(2)*ones(length(vert),1)];
c2=[c(1)*ones(length(bndmesh),1) c(2)*ones(length(bndmesh),1)];

optmesh=[ ];
for i=1:2*n+1
% scaling the boundary points 
vv=(1+xi(i))*(vert-c1)/2;
mm=(1+xi(i))*(bndmesh-c2)/2;
optmesh=[optmesh;mm];
end;

% eliminating center repetitions
optmesh=singlepts(optmesh,12);

plot(optmesh(:,1),optmesh(:,2),'.','MarkerSize',8);
axis equal;

end


function fval=f(u,v);

% % ellipse
% a=1.2;b=0.5;
% fval=u.^2/a^2+v.^2/b^2;

% quartic
a=4;b=1;
fval=a*u.^4+b*v.^4;

end


