function x = LagrangeFEM_2D(meshfile)
%       meshfile = path to the files (meshfile.node, .ele, .edge) generated
%       by Triangle.

% A function solving the problem 
%           -Div(k Grad u(x)) + q u(x) = f(x) in a given domain
%           u = g_d on part of the boundary (indicator 1)
%           du/dn = g_n on part of the boundary (indicator 2)
%           du/dn = -delta u (i.e. weakly imposing 0 bdry conditions) on
%           part of the boundary (indicator 3) <- this takes care of your
%           problem 4.

%NOTE: Again, while this function is more or less general (as you can
%      freely change k,q,f,g_d,g_n and the domain), it is written to
%      provide a simple solution to your 2D programing assignment, and
%      as such is not flawless FE code!
%      In particular, while I do assemble the global (sparse) matrix to
%      keep the code clean, I give an example (in class, and on my web
%      page), of how we can EASILY modify this to store this matrix in a
%      somewhat sparse format (using a lot less memory). 

% Coefficient definitions (the functions k,q and f are defined
%  at the end of the file.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  [V,El,Ed] = readmesh(meshfile);
  
  total_dof = size(V,1);

  % Populating the matrices and the right hand side.
  A1=zeros(total_dof,total_dof);
  A0=zeros(total_dof,total_dof); 
  B =zeros(total_dof,1);
  for i=1:size(El,1)
      % compute the affine transfromation
      EVC = V(El(i,:),1:2)'; % <-- element vertex coordinates 
      b = zeros(2,1); T = zeros(2,2); 
      b =  EVC(:,1); T(:,1) = EVC(:,2) - b; T(:,2) = EVC(:,3) - b;  
      Tinv = inv(T);
      % the transformation K^ -> K (reference to element) is now x = T*X+b.
      % the transformation K -> K^ is then X = Tinv*(x-b)
      
      [wval,tval] = quadrature(T,b);
      quadrature_length = size(wval,2);

      for q_ind = 1:quadrature_length
           % note that the "quadrature" function is mapping
           % reference quadrature points to points in the 
           % element, while here we map those points back to
           % evaluate the basis functions. This is obviously
           % redundant, and I only wrote it this way for clarity.
           x = tval(:,q_ind); X = Tinv*(x - b);
           J = abs(det(T));
           
          for b1_ind=1:3
              phi_1val = basis_value(b1_ind,X);
              phi_1grad = Tinv'*basis_grad(b1_ind,X);

              for b2_ind=1:3
                  phi_2val = basis_value(b2_ind,X);
                  phi_2grad = Tinv'*basis_grad(b2_ind,X);
                  
                  % contribution to stiffness matrix
                  A1(El(i,b1_ind),El(i,b2_ind)) = A1(El(i,b1_ind),El(i,b2_ind)) + wval(q_ind)*k(x)*(phi_1grad'*phi_2grad)*J;
                  
                  % and to mass matrix;
                  A0(El(i,b1_ind),El(i,b2_ind)) = A0(El(i,b1_ind),El(i,b2_ind)) + wval(q_ind)*q(x)*(phi_1val*phi_2val)*J;
              end
              
              B(El(i,b1_ind)) = B(El(i,b1_ind)) + wval(q_ind) * f(x) * phi_1val*J;
          end
      end
  end
  
  % the next part is a cycle over the edges to calculate the 
  %   contributions of the neumann / mixed boundary conditions
  % for simplicity I use the Simpson's rule directly (avoiding
  %   transformation to a reference edge, although implementing
  %   such a transformation is trivial - see above). 
  for i=1:size(Ed,1)
      % the points
      l = V(Ed(i,1),1:2)'; r = V(Ed(i,2),1:2)'; m = (l+r)/2;
      % length
      e = norm(r-l);
      % the term integral (over part of the boundary) g_n*v (part of the
      % RHS functional)
      if (Ed(i,3)==2)
          %RHS contributions
          B(Ed(i,1)) = B(Ed(i,1)) +e/6*(g_n(l)+2*g_n(m));
          B(Ed(i,2)) = B(Ed(i,2)) +e/6*(g_n(r)+2*g_n(m));
      end
      % the term integral (over part of the boundary) delta*u*v (part of
      % the bilinear form)
      if (Ed(i,3)==3)
          %Matrix Contributions
          A0(Ed(i,1),Ed(i,1)) = A0(Ed(i,1),Ed(i,1)) + e/6*(delta(l)+delta(m));
          A0(Ed(i,1),Ed(i,2)) = A0(Ed(i,1),Ed(i,2)) + e/6*(delta(m));
          A0(Ed(i,2),Ed(i,1)) = A0(Ed(i,2),Ed(i,1)) + e/6*(delta(m));
          A0(Ed(i,2),Ed(i,2)) = A0(Ed(i,2),Ed(i,2)) + e/6*(delta(r)+delta(m));
      end
  end
  
  % The global matrix is now:
  A=A0+A1;
  
  % Imposing the boundary conditions.
  for i=1:total_dof
      if (V(i,3)==1)
          for j=1:total_dof
              A(i,j) = 0;
              B(j) = B(j) - A(j,i)*g_d(V(i,1:2)');
              A(j,i) = 0;
          end
          A(i,i) = 1;
          B(i) = g_d(V(i,1:2)');
      end
  end
  
  % Solving,
  x=A\B

  % and plotting the solution
  plot_surf(V,El,x);
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  function y = basis_value(b,X)
    % the basis functions are 1-x-y, x, y
    values = [1-X(1)-X(2), X(1), X(2)];
    y = values(b);
  end

  function y = basis_grad(b,X)
    % the basis gradients are [-1,-1],[1,0],[0,1]
    values = [-1 -1; 1, 0; 0, 1];
    y = values(b,:)';
  end

  function [w,t] = quadrature(T,b)
      w = [1/6,1/6,1/6];
      t = zeros(2,3);
      t(:,1) = T*[0.5;0] + b; t(:,2) = T*[0.5;0.5] + b; t(:,3) = T*[0;0.5] + b;  
  end
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  % The variable coefficents of the problem we're solving.
  function y = k(x)
      y = 1;
  end
  
  function y = q(x)
      y = 1;
  end
  
  function y = f(x)
    y = x(1)*x(2);
  end
  
  function y = g_d(x)
    y=0;
  end

  function y = g_n(x)
    y=1;
  end

  function y = delta(x)
    y=1;
  end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

end