function theta = adjointAssemble_el(xinit,e,IDD,IDD1,IDD2,Ishape,shapeDisplacements,beta)

% theta = adjointAssemble_el(xinit,e,IDD,IDD1,IDD2,Ishape,shapeDisplacements,beta)
%          applies a vector beta to gradient of a right-hand side vector
%          which arises 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,


theta = zeros(size(shapeDisplacements,1)*size(shapeDisplacements,2),1);

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

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

% 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
    
    [graddivUi,gradepsUi] = graddiveps(x1,x2,x3,U1,U2,U3);

    % size(gradbei) = [6,6]
    gradbei = gradbloc_el(lambda,mu,graddivUi,gradepsUi,xinit,ei);
        
    for i1=1:3
        for j1=1:2
               
            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;
                                
                if res1==1
                    theta(2*shapeIdx1-1) = theta(2*shapeIdx1-1) + gradbei(lokindex1,1)*beta(index1);
                    theta(2*shapeIdx1) = theta(2*shapeIdx1) + gradbei(lokindex1,2)*beta(index1);
                end
                if res2==1
                    theta(2*shapeIdx2-1) = theta(2*shapeIdx2-1) + gradbei(lokindex1,3)*beta(index1);
                    theta(2*shapeIdx2) = theta(2*shapeIdx2) + gradbei(lokindex1,4)*beta(index1);
                end
                if res3==1
                    theta(2*shapeIdx3-1) = theta(2*shapeIdx3-1) + gradbei(lokindex1,5)*beta(index1);
                    theta(2*shapeIdx3) = theta(2*shapeIdx3) + gradbei(lokindex1,6)*beta(index1);
                end
              
%                b(index1) = b(index1) + bei(lokindex1);
            end
            
        end
    end
end
