function lambda = BiharmonicEig(p, t, e, num)
[S, M]=assemble(p,t); N=length(p);
%%%-------Find boundary nodes----------------------------------------------
bdnodeE = unique([e(1,:),e(2,:)]);
Inode = setdiff(linspace(1,N,N), bdnodeE);
A = S(Inode, :); 
B = M(Inode, Inode); 
% ----------------------- clamped plate eigenvalues -----------------------
disp('clamped plate eigenvalues');
[V,D]=eigs(A*inv(M)*A', B, num,'sm'); 
lambda = diag(D);

function [S, M] = assemble(p, t)
[weight,point]=quad_3;  % get 3 point quadrature rule
nq=length(weight);
yloc=zeros(3,nq);
gyloc=zeros(2,3,nq);
for r=1:nq
  [yloc(:,r),gyloc(:,:,r)]=phiRef(point(:,r));
end
nt = length(t); nv=length(p);
S=sparse(nv,nv); M=sparse(nv,nv);
for it=1:nt
    indices=t(1:3,it)';
    % The coordinates of the vertices of triangle 'it'
    v=p(:,t(1:3,it));    
    B=[v(:,2)-v(:,1),v(:,3)-v(:,1)]; detB=abs(det(B));
	for r=1:nq
        gphi(:,:,r)=(inv(B))'*gyloc(:,:,r);
	end
    % Stiffness matrix
    Sloc=zeros(3,3);
    for r=1:nq
        Sloc=Sloc+((gphi(:,:,r)'*gphi(:,:,r)))*weight(r);
    end
    Sloc=Sloc*detB/2;
    S(indices, indices) = S(indices, indices)+Sloc;
    % Mass matrix
    Mloc=zeros(3,3);
    for r=1:nq
        Mloc=Mloc+((yloc(:,r)*yloc(:,r)'))*weight(r);
    end
    Mloc=Mloc*detB/2;
    M(indices, indices) = M(indices, indices)+Mloc;
end

function [weight,point]=quad_3()
% 3 point quadrature on the reference triangle (2nd order exactly)
w1=1/3; weight=[w1,w1,w1];
point(:,1)=[0;1/2]; point(:,2)=[1/2;0]; point(:,3)=[1/2;1/2];

function [y,grady]=phiRef(xhat)
% Linear Basis Functions on the reference triangle
y=zeros(3,1); grady=zeros(2,3);
y(1) = 1 - xhat(1) - xhat(2); y(2) = xhat(1); y(3) = xhat(2);
grady(:,1) = [-1; -1]; grady(:,2) = [1; 0]; grady(:,3) = [0; 1];