function e=model_incompress(m,lambda)
% returns dispacement at the point in the middle of left side
% test for plane elasticity

%m=2;
n1=2*m; n2=3*m; h1=0.5/m;h2=0.5/m;
% lambda=100;
mu=2;
rpoint=[1.5;0.5];


jmf=model_plane(n1,n2,h1,h2,lambda,mu);
% index of node returned
rnode=find(abs(jmf.xyz(1,:)-rpoint(1))+abs(jmf.xyz(2,:)-rpoint(2))<5*eps)
rdof=jmf.id(2,rnode);

K=assemble(jmf.Kloc,jmf.ien,jmf.id);

K=0.5*(K+K');

plot_eig=0;
if plot_eig,
disp('global unconstrained matrix should have 3 zero eigenvalues')
[v,d]=eig(full(K));
e=diag(d);
e(1:5)

plot_elem(jmf.xyz,jmf.ien,'k')
scale=0.2*(h1*n1+h2*n2);
hold on, plot_disp(jmf,scale*v(:,1),'b')
hold on, plot_disp(jmf,scale*v(:,2),'r')
hold on, plot_disp(jmf,scale*v(:,3),'g')
title('model plane and eigenvectors for 3 smallest eigenvalues')
hold off
end


% mesh coordinates = x1 x2 instead of x y
% x=jmf.xyz;

% plot sides (visible in debug mode)

plot_sides=0;
if plot_sides,
plot_elem(jmf.xyz,jmf.ien,'k')
hold on,plot_elem(jmf.xyz,jmf.side{1},'r')
hold on,plot_elem(jmf.xyz,jmf.side{2},'g')
hold on,plot_elem(jmf.xyz,jmf.side{3},'m')
hold on,plot_elem(jmf.xyz,jmf.side{4},'b')
end

% get list, number, and coordinates of nodes on sides
for i=1:4
    s{i}=unique(jmf.side{i});
    o{i}=ones(size(s{i}));
    xx{i}=jmf.xyz(1,s{i})';
    yy{i}=jmf.xyz(2,s{i})';
end

% constraint equations
cnodes=[s{1};s{1};s{3};s{3}];
cdofs=[1*o{1};2*o{1};1*o{3};2*o{3}];
cvals=0.5*[0*o{1};0*o{1};1*o{3};1*xx{3}];
C=point_constraint(jmf,cnodes,cdofs);
% loads

load=zeros(size(K,1),1);

% apply constraints
[Kc,fc]=constrain(K,load,C,cvals);

% solve
u = Kc\fc;

%plot_disp(jmf,0*u,'k')
plot_disp(jmf,u,'b')

%energy=0.5*u'*K*u-u'*load;
e=u(rdof);

return