%function [u,utrue,errL2,errH1] = fem2d()
%
%  program to solve the equation
%
%     -Laplacian u = f   on  [0,1] X [0,1]
%
%  with boundary conditions
%
%           u(x,0) = u(0,y) = 0
%     u_y(x,1) = 0,  u_x(1,y) = -y(1-.5y)
%
%  using the finite element method with piecewise linear
%  approximation functions.   Uses function f
%  to evaluate the right-hand side.
%
%  User is asked for the number of subintervals n.
%  Routine plots the true and computed solutions and their
%  difference.  Returns the computed solution u, the true
%  solution utrue, the L2-norm of the error errL2, and the
%  H1-seminorm of the error errH1.
%
%  Stiffness matrix is stored as a dense matrix and solved
%  using MATLAB operator "\".
%

n = input(' Enter number of subintervals in each direction: ');
h = 1/n;
nsq = n*n;

%  Loop over squares.  Assemble global stiffness matrix A
%  and right-hand side vector b.

A = zeros(nsq,nsq);  b=zeros(nsq,1);

for j=1:n,                           % Outer loop over y
  yjm1 = (j-1)*h;
  yj = j*h;
  for i=1:n,                         % Inner loop over x
    xim1 = (i-1)*h;
    xi = i*h;

      % Handle two triangles within each square

    index_upperright = (j-1)*n + i;          % Indices of nodes in
    index_upperleft = index_upperright - 1;  % each triangle
    index_lowerright = index_upperright - n;
    index_lowerleft = index_lowerright - 1; 
    area = .5*h^2;                           % Area of each triangle

    if i > 1 & j > 1,           % Lower left node is a free node
      k = index_lowerleft;
      A(k,k) = A(k,k) + area*(1/h)^2;   % Contribution from right triangle
      A(k,k+1) = A(k,k+1) - area*(1/h)^2; 
      A(k,k) = A(k,k) + area*(1/h)^2;   % Contribution from left triangle
      A(k,k+n) = A(k,k+n) - area*(1/h)^2;
    end;
    if j > 1,                   % Lower right node is a free node
      k = index_lowerright;
      A(k,k) = A(k,k) + area*((1/h)^2 + (1/h)^2);
      if i > 1, A(k,k-1) = A(k,k-1) - area*(1/h)^2; end;
      A(k,k+n) = A(k,k+n) - area*(1/h)^2;
    end;
    if i > 1,                   % Upper left node is a free node
      k = index_upperleft;
      A(k,k) = A(k,k) + area*((1/h)^2 + (1/h)^2);
      if j > 1, A(k,k-n) = A(k,k-n) - area*(1/h)^2; end;
      A(k,k+1) = A(k,k+1) - area*(1/h)^2;
    end;
                                % Upper right node is always a free node
    k = index_upperright;
    A(k,k) = A(k,k) + area*(1/h)^2;      % Contribution from right triangle
    if j > 1, A(k,k-n) = A(k,k-n) - area*(1/h)^2; end;
    A(k,k) = A(k,k) + area*(1/h)^2;      % Contribution from left triangle
    if i > 1, A(k,k-1) = A(k,k-1) - area*(1/h)^2; end;

%         Right-hand side vector.  Use second order quadrature formula
%         involving function values at midpoints of each side.

      midptr1 = [.5*(xi+xim1),yjm1];    % Midpoints of sides of right triangle
      midptr2 = [xi,.5*(yj+yjm1)];
      midptr3 = [.5*(xi+xim1),.5*(yj+yjm1)];
      midptl1 = [xim1,.5*(yj+yjm1)];    % Midpoints of sides of left triangle
      midptl2 = [.5*(xi+xim1),yj];
      midptl3 = midptr3;
      fr1 = f(midptr1);                 % Evaluate f at midpoints.
      fr2 = f(midptr2);
      fr3 = f(midptr3);
      fl1 = f(midptl1);
      fl2 = f(midptl2);
      fl3 = f(midptl3);
            
%          Translate triangle to have lower left corner at (0,0).
      midptr1 = midptr1 - [xim1,yjm1];
      midptr2 = midptr2 - [xim1,yjm1];
      midptr3 = midptr3 - [xim1,yjm1];
      midptl1 = midptl1 - [xim1,yjm1];
      midptl2 = midptl2 - [xim1,yjm1];
      midptl3 = midptl3 - [xim1,yjm1];

    if i > 1 & j > 1,           % Lower left node is a free node
      k = index_lowerleft;
      phir1 = 1 - midptr1(1)/h;         % Contribution from right triangle
      phir2 = 1 - midptr2(1)/h;
      phir3 = 1 - midptr3(1)/h;
      b(k) = b(k) + (area/3)*(fr1*phir1 + fr2*phir2 + fr3*phir3);
      phil1 = 1 - midptl1(2)/h;         % Contribution from left triangle
      phil2 = 1 - midptl2(2)/h;
      phil3 = 1 - midptl3(2)/h;
      b(k) = b(k) + (area/3)*(fl1*phil1 + fl2*phil2 + fl3*phil3);
    end;
    if j > 1,                   % Lower right node is a free node
      k = index_lowerright;
      phir1 = (1/h)*midptr1(1) - (1/h)*midptr1(2);
      phir2 = (1/h)*midptr2(1) - (1/h)*midptr2(2);
      phir3 = (1/h)*midptr3(1) - (1/h)*midptr3(2);
      b(k) = b(k) + (area/3)*(fr1*phir1 + fr2*phir2 + fr3*phir3);
    end;
    if i > 1,                   % Upper left node is a free node
      k = index_upperleft;
      phil1 = -(1/h)*midptl1(1) + (1/h)*midptl1(2);
      phil2 = -(1/h)*midptl2(1) + (1/h)*midptl2(2);
      phil3 = -(1/h)*midptl3(1) + (1/h)*midptl3(2);
      b(k) = b(k) + (area/3)*(fl1*phil1 + fl2*phil2 + fl3*phil3);
    end;
                                % Upper right node is always a free node
    k = index_upperright;
    phir1 = (1/h)*midptr1(2);            % Contribution from right triangle
    phir2 = (1/h)*midptr2(2);
    phir3 = (1/h)*midptr3(2);
    b(k) = b(k) + (area/3)*(fr1*phir1 + fr2*phir2 + fr3*phir3);
    phil1 = (1/h)*midptl1(1);            % Contribution from left triangle
    phil2 = (1/h)*midptl2(1);
    phil3 = (1/h)*midptl3(1);
    b(k) = b(k) + (area/3)*(fl1*phil1 + fl2*phil2 + fl3*phil3);

%         Extra boundary term.  Use 2-point Gauss quadrature here.
    if i==n,
      pt1 = .5*(yj+yjm1) - h/sqrt(3);    % Gauss quadrature points
      pt2 = .5*(yj+yjm1) + h/sqrt(3);
      usubn1 = -pt1*(1-.5*pt1);          % u_n at quadrature points
      usubn2 = -pt2*(1-.5*pt2);
      pt1 = pt1 - yjm1;                  % Translate interval
      pt2 = pt2 - yjm1;

      k = index_upperright;              % Upper right basis function
      phi1 = (1/h)*pt1;
      phi2 = (1/h)*pt2;
      b(k) = b(k) + .5*h*(usubn1*phi1 + usubn2*phi2);
      if j > 1,              
        k = index_lowerright;            % Lower right basis function
        phi1 = 1 - (1/h)*pt1;
        phi2 = 1 - (1/h)*pt2;
        b(k) = b(k) + .5*h*(usubn1*phi1 + usubn2*phi2);
      end;
    end;
      
  end;
end;

%  Solve linear system.

u = A\b;

%  Plot true and approximate solutions and error.

u2d = reshape(u,n,n);
u_true = zeros(n,n);
for j=1:n,
  yj = j*h;
  for i=1:n,
    xi = i*h;
    u_true(i,j) = xi*(1-xi)*yj*(1-.5*yj);
  end;
end;

surf(u2d)
xlabel('x')
ylabel('y')
zlabel('ucomp')
pause(5)
figure(2)
surf(u_true)
xlabel('x')
ylabel('y')
zlabel('utrue')
pause(5)
figure(3)
surf(u_true-u2d)
xlabel('x')
ylabel('y')
zlabel('error')
pause(5)

%  Compute (approximate) L2-norm of error.
%  Replace integrals by second order quadrature rule involving
%  function values at the midpoints of sides of triangles.

errL2 = 0;  
for j=1:n,
  yj = j*h;
  yjm1 = (j-1)*h;
  for i=1:n,
    xi = i*h;
    xim1 = (i-1)*h;

    midptr1 = [.5*(xi+xim1),yjm1];    % Midpoints of sides of right triangle
    midptr2 = [xi,.5*(yj+yjm1)];
    midptr3 = [.5*(xi+xim1),.5*(yj+yjm1)];
    midptl1 = [xim1,.5*(yj+yjm1)];    % Midpoints of sides of left triangle
    midptl2 = [.5*(xi+xim1),yj];
    midptl3 = midptr3;
%
    utruer1 = midptr1(1)*(1-midptr1(1))*midptr1(2)*(1-.5*midptr1(2));
    utruer2 = midptr2(1)*(1-midptr2(1))*midptr2(2)*(1-.5*midptr2(2));
    utruer3 = midptr3(1)*(1-midptr3(1))*midptr3(2)*(1-.5*midptr3(2));
    utruel1 = midptl1(1)*(1-midptl1(1))*midptl1(2)*(1-.5*midptl1(2));
    utruel2 = midptl2(1)*(1-midptl2(1))*midptl2(2)*(1-.5*midptl2(2));
    utruel3 = midptl3(1)*(1-midptl3(1))*midptl3(2)*(1-.5*midptl3(2));
%
    if i > 1 & j > 1, 
      ur1 = .5*(u2d(i,j-1)+u2d(i-1,j-1)); 
    elseif i==1 & j > 1,
      ur1 = .5*u2d(i,j-1);
    else
      ur1 = 0;
    end;
    if j > 1, 
      ur2 = .5*(u2d(i,j)+u2d(i,j-1));
    else
      ur2 = .5*u2d(i,j);
    end;
    if i > 1 & j > 1,
      ur3 = .5*(u2d(i,j)+u2d(i-1,j-1));
    else
      ur3 = .5*u2d(i,j);
    end;
    if i > 1 & j > 1,
      ul1 = .5*(u2d(i-1,j)+u2d(i-1,j-1));
    elseif j==1 & i > 1,
      ul1 = .5*u2d(i-1,j);
    else
      ul1 = 0;
    end;
    if i > 1,
      ul2 = .5*(u2d(i,j)+u2d(i-1,j));
    else
      ul2 = .5*u2d(i,j);
    end;
    ul3 = ur3;

    errL2 = errL2 + (area/3)*((utruer1-ur1)^2 + (utruer2-ur2)^2 + ...
                              (utruer3-ur3)^2);
    errL2 = errL2 + (area/3)*((utruel1-ul1)^2 + (utruel2-ul2)^2 + ...
                              (utruel3-ul3)^2);
  end;
end;
errL2 = sqrt(errL2),