function [K,b] = assemble_el(xinit,e,IDD,IDD1,IDD2,Ishape,shapeDisplacements)

% [K,b] = assemble_el(xinit,e,IDD,IDD1,IDD2,Ishape,shapeDisplacements)
%          assembles a stiffness matrix and a right-hand side vector
%          both of which arise from a variational formulation
%          of the 2d linear elasticity that determines a shape-to-mesh
%          mapping for optimal shape design.
%          The variational formulation is as follows:
%           Find uH from Vh : a(uH,v) = f(v) for all v in Vh,
%          where
%           a(v,u) = \int_{Omega} lambda(x)*div(uH(x))*div(v(x)) dx +
%                    + \int_{Omega} mu(x) * \sum_{k,l=1}^3 eps_{k,l}(uH(x))*
%                                                          eps_{k,l}(v(x)) dx,
%           f(v) = - \int_{Omega} lambda(x)*div(U(x))*div(v(x)) dx -
%                  - \int_{Omega} mu(x)* sum_{k,l=1}^3 eps_{k,l}(U(x))*
%                                                      eps_{k,l}(v(x)) dx,

%sprintf('%s','assemble_el')

global K b xinit_last KSolved;

if max(size(xinit_last)~=size(xinit))==1
    KSolved = 0;
elseif max(max(xinit_last~=xinit))==1
    KSolved = 0;
end


% number of nodes (number of degrees of freedom divided by 3)
n = max(size(xinit));

% number of elements
ne = max(size(e));

if KSolved==0
    K = sparse(2*n,2*n);
end

b = zeros(2*n,1);

% Lame's coefficients
E = 100e9;
nu = 0.25;
lambda = (nu*E)/((1+nu)*(1-2*nu));
mu = E/(1+nu);


for i=1:ne
    ei = e(i,:);
    x1 = xinit(:,ei(1));
    x2 = xinit(:,ei(2));
    x3 = xinit(:,ei(3));
    
    Ri = [x2-x1 x3-x1];
    ri = x1;
    
    massPointi = Ri*[1/3;1/3] + ri;
    
    [res1,shapeIdx1] = in(Ishape,ei(1));
    [res2,shapeIdx2] = in(Ishape,ei(2));
    [res3,shapeIdx3] = in(Ishape,ei(3));
    
    U1 = zeros(2,1);
    U2 = zeros(2,1);
    U3 = zeros(2,1);
    
    if res1==1
        U1 = shapeDisplacements(:,shapeIdx1);
    end
    
    if res2==1
        U2 = shapeDisplacements(:,shapeIdx2);
    end
    
    if res3==1
        U3 = shapeDisplacements(:,shapeIdx3);
    end
    
    [divUi,epsUi] = diveps(x1,x2,x3,U1,U2,U3);
    
    if KSolved==0
        Aei = aloc_el(lambda,mu,xinit,ei);
    end

    bei = bloc_el(lambda,mu,divUi,epsUi,xinit,ei);
        
    for i1=1:3
        for j1=1:2
          if KSolved==0
            for i2=1:3
                for j2=1:2
                    if (i1==i2 & j1==j2) | (~in(IDD,ei(i1)) & ~in(IDD,ei(i2)) & ~(j1==1 & in(IDD1,ei(i1))) & ~(j1==2 & in(IDD2,ei(i1))) & ~(j2==1 & in(IDD1,ei(i2))) & ~(j2==2 & in(IDD2,ei(i2))))
                        index1 = (ei(i1)-1)*2 + j1;
                        lokindex1 = (i1-1)*2 + j1;
                        index2 = (ei(i2)-1)*2 + j2;
                        lokindex2 = (i2-1)*2 + j2;

                        K(index1,index2) = K(index1,index2) + Aei(lokindex1,lokindex2);
                    end
                end
            end
          end
               
            if ~in(IDD,ei(i1)) & ~(j1==1 & in(IDD1,ei(i1))) & ~(j1==2 & in(IDD2,ei(i1)))
                index1 = (ei(i1)-1)*2 + j1;
                lokindex1 = (i1-1)*2 + j1;
                
                b(index1) = b(index1) + bei(lokindex1);
            end
            
        end
    end
end

KSolved = 1;
