%% fem_50d applies the finite element method to Laplace's equation.
function fem_50d
% input  a prefix, e.g. A7_
% files read  A7_coord.dat   one x y pair per line, for nodes in order  1, 2, ...
%             A7_elem3.dat   three node numbers per line, triangles, any order
%                            function stima3 is applied to these cooddinates
%             A7_elem4.dat   four node numbers per line, quadralaterals
%                            function stima4 is applied to these coordinates
%             A7_dirich.dat  two node numbers per line, dirichlet boundary
%                            function u_d is applied to these coordinates
%             A7_neum.dat    two node numbers per line, neumann bounday
%                            function g is applied to these coordinates
%
%  Discussion:
%    FEM_50d is a set of MATLAB routines to apply the finite
%    element method to solving Laplace's equation in an arbitrary
%    region, using about 50 lines of MATLAB code.
%    
%    FEM_50d is partly a demonstration, to show how little it
%    takes to implement the finite element method (at least using
%    every possible MATLAB shortcut.)  The user supplies datafiles
%    that specify the geometry of the region and its arrangement
%    into triangular and quadrilateral elements, and the location
%    and type of the boundary conditions, which can be any mixture
%    of Neumann and dirichlet.
%    
%    The unknown state variable U(x,y) is assumed to satisfy
%    Laplace's equation:
%      -Uxx(x,y) - Uyy(x,y) = F(x,y) in Omega
%    with dirichlet boundary conditions
%      U(x,y) = U_D(x,y) on Gamma_D
%    and Neumann boundary conditions on the outward normal derivative:
%      Un(x,y) = G(x,y) on Gamma_N
%    If Gamma designates the boundary of the region Omega,
%    then we presume that
%      Gamma = Gamma_D + Gamma_N
%    but the user is free to determine which boundary conditions to
%    apply.  Note, however, that the problem will generally be singular
%    unless at least one dirichlet boundary condition is specified.
%    
%    The code uses piecewise linear basis functions for triangular elements,
%    and piecewise isoparametric bilinear basis functions for quadrilateral
%    elements.
%    
%    The user is required to supply a number of data files and MATLAB
%    functions that specify the location of nodes, the grouping of nodes
%    into elements, the location and value of boundary conditions, and 
%    the right hand side function in Laplace's equation.  Note that the
%    fact that the geometry is completely up to the user means that
%    just about any two dimensional region can be handled, with arbitrary
%    shape, including holes and islands.
%    
%  Modified:
%    29 March 2004
%    23 February 2008 JSS
%     3 March 2008 JSS
%  Reference:
%    Jochen Alberty, Carsten Carstensen, Stefan Funken,
%    Remarks Around 50 Lines of MATLAB:
%    Short Finite Element Implementation,
%    Numerical Algorithms,
%    Volume 20, pages 117-137, 1999.
%
  clear
  format compact
  prename = input('enter prefix to file names: ', 's')
  diary 'fem_50d.out'
  disp([prename 'fem_50d.out'])
  debug = input('enter debug level:0=none, 1=input and matrices, 2=assembly ')
%
%  Read the nodal coordinatesinate data file.
%  ordered list of  x  y  pairs, coordinates for nodes
%  nodes must be sequentially numbered 1, 2, 3, ...
  coordinates = load(strcat(prename,'coord.dat'));
  if debug>0
    disp([prename 'coord.dat ='])
	disp(coordinates)
  end % debug
%
%  Read the triangular element data file.
%  three integer node numbers per line
  elements3 = load(strcat(prename,'elem3.dat'));
  if debug>0
    disp([prename 'elem3.dat ='])
	disp(elements3)
  end % debug
%
%  Read the quadrilateral element data file.
%  four integer node numbers per line
  elements4 = load(strcat(prename,'elem4.dat'));
  if debug>0
    disp([prename 'elem4.dat ='])
	disp(elements4)
  end % debug
%
%  Read the dirichlet boundary edge data file.
%  two integer node numbers per line, function u_d sets values
%  there must be at least one pair
%  other boundary edges are set in  neumann
  dirichlet = load(strcat(prename,'dirich.dat'));
  if debug>0
    disp([prename 'dirich.dat ='])
	disp(dirichlet)
  end % debug
%
%  Read the Neumann boundary edge data file.
%  two integer node numbers per line, function g  sets values
  neumann = load(strcat(prename,'neum.dat'));
  if debug>0
    disp([prename 'neum.dat ='])
	disp(neumann)
  end % debug

  A = sparse ( size(coordinates,1), size(coordinates,1) );
  b = sparse ( size(coordinates,1), 1 );

%
%  Assembly from triangles.
%
  for j = 1 : size(elements3,1)
    VV = coordinates(elements3(j,:),:);
    MM = stima3(coordinates(elements3(j,:),:));
	if debug>1
	  disp(['MM = stima3(VV) at j=' num2str(j) ' nodes(' ...
	        num2str(elements3(j,1)) ',' num2str(elements3(j,2)) ',' ...
			num2str(elements3(j,3)) ')'])
	  disp('VV=')
	  disp(VV)
	  disp('MM=')
	  disp(MM)
	end % debug
    A(elements3(j,:),elements3(j,:)) = A(elements3(j,:),elements3(j,:)) + MM;
  end
  if debug>0
    disp('assembly from triangles A=')
    disp(A)
  end % debug
%
%  Assembly from quadralaterals.
%
  for j = 1 : size(elements4,1)
    VV = coordinates(elements3(j,:),:);
    MM = stima4(coordinates(elements4(j,:),:));
	if debug>1
	  disp(['MM = stima3(VV) at j=' num2str(j) ' nodes(' ...
	        num2str(elements4(j,1)) ',' num2str(elements4(j,2)) ',' ...
			num2str(elements4(j,3)) ',' num2str(elements4(j,4)) ')'])
	  disp('VV=')
	  disp(VV)
	  disp('MM=')
	  disp(MM)
	end % debug
    A(elements4(j,:),elements4(j,:)) = A(elements4(j,:),elements4(j,:)) + MM;
  end
  if debug>0
    disp('assembly from triangles and rectangles A=')
    disp(A)
  end % debug

%
%  Volume Forces from triangles.
%
  for j = 1 : size(elements3,1)
    b(elements3(j,:)) = b(elements3(j,:)) ...
      + det( [1,1,1; coordinates(elements3(j,:),:)'] ) * ...
      f(sum(coordinates(elements3(j,:),:))/3)/6;
  end
  if debug>0
    disp('forces from triangles b=')
    disp(b)
  end % debug
%
%  Volume Forces from quadralaterals.
%
  for j = 1 : size(elements4,1)
    b(elements4(j,:)) = b(elements4(j,:)) ...
      + det([1,1,1; coordinates(elements4(j,1:3),:)'] ) * ...
      f(sum(coordinates(elements4(j,:),:))/4)/4;
  end
  if debug>0
    disp('forces from triangles and rectangles b=')
    disp(b)
  end % debug

%
%  Neumann conditions.
%
  for j = 1 : size(neumann,1)
    GG = norm(coordinates(neumann(j,1),:) - coordinates(neumann(j,2),:)) * ...
         g(sum(coordinates(neumann(j,:),:))/2)/2;
    if debug>2
	  disp(['Neumann at j=' num2str(j)])
      disp(GG)
	end % debug
    b(neumann(j,:)) = b(neumann(j,:)) + GG
  end
  if debug>0
    disp('using g() add neumann conditions b=')
    disp(b)
  end % debug

%
%  Determine which nodes are associated with dirichlet conditions.
%  Assign the corresponding entries of U, and adjust the right hand side.
%
  u = sparse ( size(coordinates,1), 1 );
  BoundNodes = unique ( dirichlet );
  if debug>1
    disp('BoundNodes=')
    disp(BoundNodes)
  end % debug
  u(BoundNodes) = u_d ( coordinates(BoundNodes,:) );
  if debug>0
    disp('using u_d() add dirichlet conditions u=')
    disp(u)
  end % debug

  b = b - A * u;
  if debug>0
    disp('apply u to b, dirichlet conditions b=')
    disp(b)
  end % debug

%
%  Compute the solution by solving A * u = b 
%  for the, un bound, remaining unknown values of u.
%
  FreeNodes = setdiff ( 1:size(coordinates,1), BoundNodes );

  u(FreeNodes) = A(FreeNodes,FreeNodes) \ b(FreeNodes);
  if debug>0
    disp('solve for A * u = b  u=')
    disp(u)
  end % debug
  
%
%  Graphic representation.
%
  show ( elements3, elements4, coordinates, full(u), dirichlet );
%
%  check solution
%
  for i=1:size(coordinates,1)
    disp([num2str(i) '  x= ' num2str(coordinates(i,1)) '  y= '  num2str(coordinates(i,2)) ...
	      ' analytic solution= ' num2str(uana(coordinates(i,1),coordinates(i,2)))])
  end
  diary off
  return % logical end of function  fem_50d
  

%% uana computes boundary values and analytic solution for checking
%  Discussion:
%    This is generally unknown yet is used here for testing
%  Parameters:
%    Input, real pair  xx,yy are x,y coordinates
%    Output, value of solution at  x,y
%
  function ana = uana(xx, yy)
    ana= 1.0+(2.0/11.0)*xx*xx+(3.0/11.0)*yy*yy;
  end % uana


  %% f evaluates the right hand side of Laplace's equation.
  %  Discussion:
  %    This routine must be changed by the user to reflect a particular problem.
  %  Parameters:
  %    Input, real U(N,M), contains the M-dimensional coordinates of N points.
  %    Output, VALUE(N), contains the value of the right hand side of Laplace's
  %    equation at each of the points.
  function valuef = f ( uf )
    valuef = ones ( size ( uf, 1 ), 1 );
    valuef = valuef.*(-10/11.0);
  end % f


  %% g evaluates the outward normal values assigned at Neumann boundary conditions.
  %  Discussion:
  %    This routine must be changed by the user to reflect a particular problem.
  %  Parameters:
  %    Input, real U(N,M), contains the M-dimensional coordinates of N points.
  %    Output, VALUE(N), contains the value of outward normal at each point
  %    where a Neumann boundary condition is applied.
  function valueg = g ( ug )
    valueg = zeros ( size ( ug, 1 ), 1 );
  end % g


  %% u_d evaluates the dirichlet boundary conditions.
  %  Discussion:
  %    The user must supply the appropriate routine for a given problem
  %  Parameters:
  %    Input, real U(N,M), contains the M-dimensional coordinates of N points.
  %    Output, VALUE(N), contains the value of the dirichlet boundary
  %    condition at each point.
  function valued = u_d ( ud )
    valued = zeros ( size ( ud, 1 ), 1 );
    for kk=1:size(ud,1)
      % U(x,y) = 1 + 2/11 x^2 + 3/11 y^2 solution values on boundary
      valued(kk)=uana(ud(kk,1),ud(kk,2));
    end
  end % u_d


  %% STIMA3 determines the local stiffness matrix for a triangular element.
  %  Discussion:
  %
  %    Although this routine is intended for 2D usage, the same formulas
  %    work for tetrahedral elements in 3D.  The spatial dimension intended
  %    is determined implicitly, from the spatial dimension of the vertices.
  %  Parameters:
  %
  %    Input, real VERTICES(1:(D+1),1:D), contains the D-dimensional 
  %    coordinates of the vertices.
  %
  %    Output, real M(1:(D+1),1:(D+1)), the local stiffness matrix 
  %    for the element.
  function M = stima3 ( vertices )
    d = size ( vertices, 2 );
    D_eta = [ ones(1,d+1); vertices' ] \ [ zeros(1,d); eye(d) ];
    M = det ( [ ones(1,d+1); vertices' ] ) * D_eta * D_eta' / prod ( 1:d );
  end % stima3
  

  %% STIMA4 determines the local stiffness matrix for a quadrilateral element.
  %  Parameters:
  %    Input, real VERTICES(1:4,1:2), contains the coordinates of the vertices.
  %    Output, real M(1:4,1:4), the local stiffness matrix for the element.%
  function M = stima4 ( vertices )
    D_Phi = [ vertices(2,:) - vertices(1,:); vertices(4,:) - vertices(1,:) ]';
    B = inv ( D_Phi' * D_Phi );
    C1 =  [  2, -2; -2,  2 ] * B(1,1) ...
        + [  3,  0;  0, -3 ] * B(1,2) ...
        + [  2,  1;  1,  2 ] * B(2,2);
    C2 =  [ -1,  1;  1, -1 ] * B(1,1) ...
        + [ -3,  0;  0,  3 ] * B(1,2) ...
        + [ -1, -2; -2, -1 ] * B(2,2);
    M = det ( D_Phi ) * [ C1 C2; C2 C1 ] / 6;
  end % stima4


  %% SHOW displays the solution of the finite element computation.
  %  Parameters:
  %    Input, integer elements3(N3,3), the nodes that make up each triangle.
  %    Input, integer elements4(N4,4), the nodes that make up each quadrilateral.
  %    Input, real coordinates(N,1:2), the coordinates of each node.
  %    Input, real U(N), the finite element coefficients which represent the solution.
  %    There is one coefficient associated with each node.
  function show ( elements3, elements4, coordinates, us, dirichlet )
    figure(1)
    hold off
  %
  %  Display the information associated with triangular elements.
    trisurf ( elements3, coordinates(:,1), coordinates(:,2), us' );
  %
  %  Retain the previous image, and overlay the information associated
  %  with quadrilateral elements.
    hold on
    trisurf ( elements4, coordinates(:,1), coordinates(:,2), us' );
  %
  %  Define the initial viewing angle.
  %
    view ( -67.5, 30 );  
    title ( 'Solution to the Problem' )
    hold off
    figure(2)
    for kk=1:size(elements3,1)
      plot([coordinates(elements3(kk,1),1)   ...
            coordinates(elements3(kk,2),1)   ...
            coordinates(elements3(kk,3),1)   ...
            coordinates(elements3(kk,1),1)], ...
           [coordinates(elements3(kk,1),2)   ...
            coordinates(elements3(kk,2),2)   ...
            coordinates(elements3(kk,3),2)   ...
            coordinates(elements3(kk,1),2)],':r')
      hold on
    end
    for kk=1:size(dirichlet,1)
      plot([coordinates(dirichlet(kk,1),1)   ...
            coordinates(dirichlet(kk,2),1)], ...
           [coordinates(dirichlet(kk,1),2)   ...
            coordinates(dirichlet(kk,2),2)],'-b')
      hold on
    end
    title('triangularization')
    xlabel('X')
    ylabel('Y')
    axis tight
    axis square
    grid off
    hold on
    xb=coordinates(BoundNodes,1);
	yb=coordinates(BoundNodes,2);
	plot(xb,yb,'.b')
	hold on
    xb=coordinates(FreeNodes,1);
	yb=coordinates(FreeNodes,2);
	plot(xb,yb,'.r')
	hold off
  end % show
end % fem_50d