%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function finite_element_example
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Solve Laplace equation using finite element method.  Based on
%% Burden and Faires, Chapter 12.4.  Simplified to include only
%% specified value b.c rather than a combination of value and flux.
%% Eli, APM202, 201302

%% INPUT: 
%% (A) grid specification in specify_grid: vertices location, which
%%     vertices are in what triangle element, and which vertices
%%     are on the boundary.
%% (B) boundary conditions in set_boundary_conditions: value of
%%     function specified at each bondary vertex; see also
%%     the function analytic_soln.
%% (C) The PDE to be solved is specified to be Laplace equation via
%%     the specification of the matrix equation in
%%     calculate_matrix_eqn.  May be changed there to other
%%     elliptic equations. 

clear all; clf;
set(0,'defaulttextfontsize',12)
set(0,'defaultaxesfontsize',12)

[V,T,bndry_vertices,p]=specify_grid;
plot_elements_and_vertices(V,T,bndry_vertices,p);
[a,b,c]=calculate_basis_functions(p);
gamma=set_boundary_conditions(bndry_vertices,V,p);
[A,rhs]=calculate_matrix_eqn(a,b,c,T,p,gamma);
gamma=solve_matrix_eqn_and_plot(A,rhs,V,gamma,p);

fclose(p.fid);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [V,T,bndry_vertices,p]=specify_grid
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% output: g: geometry of grid; p: parameters; see below for others.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% <<<< start of user input for grid specification >>>> %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% specify node/vertex locations:  MUST specify the interior
%% vertices first!
V=[0.2 0.2; %1
   0.4 0.2; %2
   0.3 0.1; %3
   0.5 0.15; %4
   0.3 0.3; %5
   0.05 0.35; %6
   0.0 0.2; %7
   0.075 0.075; %8
   0.2 0.0; %9
   0.4 0.0; %10
   0.55 0.05; %11
   0.45 0.35; %12
   0.2 0.4; %13
   0.6 0.2]; %14

% Specify which three vertices compose triangle it in T(it,:).  The
% three basis functions \phi_{1,2,3} on triangle element i are setup
% such that \phi_1 is equal to 1 on the first vertex listed in T(i,:)
% and zero on the other two vertices; similarly, \phi_2 is 1 on the
% second vertex, and \phi_3 is 1 on the third.
T=[1 3 9; %triangle element #1
   2 3 10; %2
   1 6 7; %3
   1 2 3; %4
   2 4 10; %5
   4 14 11; %6
   1 7 8; %7
   1 8 9; %8
   3 9 10; %9
   4 10 11; %10
   2 4 12; %11
   4 14 12; %12
   1 6 13; %13
   1 2 5; %14
   1 13 5; %15
   2 12 5; %16
   12 13 5]; %17

%% specify which vertices lie on the boundary: these MUST be
%% specified AFTER the interior vertices in V:
bndry_vertices=[6:14];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% <<<< end of user input for grid specification >>>> %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

p.Ntriangles=length(T(:,1));
p.Nvertices=length(V(:,1));
p.Nbndry_vertices=length(bndry_vertices);
p.Ninterior_vertices=p.Nvertices-p.Nbndry_vertices;

%% find list of interior vertices:
interior_vertices=1:p.Ninterior_vertices;

% Find which triangles contain each vertex (max 10 trianlges per
% vertex), store it in p.Triangles_with_vrtx(iv,:); also set
% p.nTriangles_with_vrtx(iv) to be the number of triangles which contain
% this vertex:
p.nTriangles_with_vrtx=zeros(p.Nvertices,1); 
p.Triangles_with_vrtx=zeros(p.Nvertices,10); 
for it=1:p.Ntriangles
  for jv=1:3
    nvertex=T(it,jv);
    ntriangles_in_vertex=p.nTriangles_with_vrtx(nvertex);
    p.nTriangles_with_vrtx(nvertex)=p.nTriangles_with_vrtx(nvertex)+1;
    p.Triangles_with_vrtx(nvertex,ntriangles_in_vertex+1)=it;
  end
end

%% find the locations of the three vertices for triangle (i),
%% [p.x1(i),p.y1(i)], [p.x2(i),p.y2(i)], [p.x3(i),p.y3(i)]:
for i=1:p.Ntriangles
  p.x1(i)=V(T(i,1),1);p.x2(i)=V(T(i,2),1);p.x3(i)=V(T(i,3),1);
  p.y1(i)=V(T(i,1),2);p.y2(i)=V(T(i,2),2);p.y3(i)=V(T(i,3),2);
end

%% calculate center ("incenter") of the triangles for plotting
%% purposes:
C=zeros(p.Ntriangles,2);
for i=1:p.Ntriangles
  % the length of the sides opposite to the three vertices:
  L1=sqrt((p.x2(i)-p.x3(i))^2+(p.y2(i)-p.y3(i))^2);
  L2=sqrt((p.x1(i)-p.x3(i))^2+(p.y1(i)-p.y3(i))^2);
  L3=sqrt((p.x2(i)-p.x1(i))^2+(p.y2(i)-p.y1(i))^2);
  % the center:
  p.cntr(i,:)=(L1*[p.x1(i),p.y1(i)]+L2*[p.x2(i),p.y2(i)]+L3*[p.x3(i),p.y3(i)]) ...
         /(L1+L2+L3);
end

%% calculate area of each element:
for i=1:p.Ntriangles
  p.area(i)=abs(p.x1(i)*(p.y2(i)-p.y3(i))+p.x2(i)*(p.y3(i)-p.y1(i))+p.x3(i)*(p.y1(i)-p.y2(i)))/2;
end

%% print grid info diagnostics:
p.fid=fopen('Output/finite_element_info.txt','w');
fprintf(p.fid,'** which vertices compose each triangle:\n');
for i=1:p.Ntriangles; fprintf(p.fid,'T(%d,:)=%d, %d, %d;\n',i,T(i,:)); end
fprintf(p.fid,'** area of each triangle:\n');
for i=1:p.Ntriangles; fprintf(p.fid,'p.area(%d)=%g;\n',i,p.area(i)); end
fprintf(p.fid,'** which triangles contain each vertex:\n');
for iv=1:p.Nvertices; 
  ntriangles_per_vrtx=p.nTriangles_with_vrtx(iv);
  fprintf(p.fid,'p.Triangles_with_vrtx(%d,:)=',iv);
  fprintf(p.fid,'%d,',p.Triangles_with_vrtx(iv,1:ntriangles_per_vrtx));
  fprintf(p.fid,'\n');
end
fprintf(p.fid,'** coordinates of each vertex:\n');
for i=1:p.Nvertices; fprintf(p.fid,'vertex coordinates: V(%d,:)=%g, %g;\n',i,V(i,:)); end
fprintf(p.fid,'** the boundary vertices are:\n');
fprintf(p.fid,'bndry_vertices(:)=');fprintf(p.fid,'%d,',bndry_vertices(:));fprintf(p.fid,'\n');
fprintf(p.fid,'** the interior vertices are:\n');
fprintf(p.fid,'interior_vertices(:)=');fprintf(p.fid,'%d,',interior_vertices(:));fprintf(p.fid,'\n');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function plot_elements_and_vertices(V,T,bndry_vertices,p)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Plot triangles and mark vertices:
figure(1); clf; subplot(1,2,1);
for i=1:p.Ntriangles
  for j=1:3
    hold on
    % plot triangle borders:
    plot([V(T(i,1),1),V(T(i,2),1)] ...
         ,[V(T(i,1),2),V(T(i,2),2)],'r-');
    plot([V(T(i,2),1),V(T(i,3),1)] ...
         ,[V(T(i,2),2),V(T(i,3),2)],'r-');
    plot([V(T(i,3),1),V(T(i,1),1)] ...
         ,[V(T(i,3),2),V(T(i,1),2)],'r-');
  end
end

%% plot vertices:
for i=1:p.Nvertices
  plot(V(i,1),V(i,2),'r*','markersize',12)
end

%% highlight boundary vertices:
for i=1:length(bndry_vertices)
  plot(V(bndry_vertices(i),1),V(bndry_vertices(i),2),'bo','markersize',12)
end

%% write triangle number in each of the triangles:
for i=1:p.Ntriangles
  text(p.cntr(i,1),p.cntr(i,2),sprintf('T_{%d}',i));
end

%% mark vertices with their numbers:
for i=1:p.Nvertices
  offset=0.007;
  text(V(i,1)+offset,V(i,2)+offset,sprintf('E_{%d}',i));
end
xlim([0 0.7]);
ylim([0 0.5]);
box on
xlabel('x');
ylabel('y');
title('Triangle elements; bndary (o) & interior (*) vertices')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [a,b,c]=calculate_basis_functions(p)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% calculate three basis functions for each triangle element.  The first
% basis function on triangle it is given, for example, by
% a(it,1)+b(it,1)*x+c(it,1)*y.  See explanatory comment before the
% definition of T(i,j) regarding the order of the basis functions.
a=zeros(p.Ntriangles,3); b=a; c=a;
for i=1:p.Ntriangles
  A=[1 p.x1(i) p.y1(i); 1 p.x2(i) p.y2(i); 1 p.x3(i) p.y3(i)];
  TMP=A\[1 0 0]';
  a(i,1)=TMP(1); b(i,1)=TMP(2); c(i,1)=TMP(3);
  TMP=A\[0 1 0]';
  a(i,2)=TMP(1); b(i,2)=TMP(2); c(i,2)=TMP(3);
  TMP=A\[0 0 1]';
  a(i,3)=TMP(1); b(i,3)=TMP(2); c(i,3)=TMP(3);
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function gamma=set_boundary_conditions(bndry_vertices,V,p)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Specify b.c here, by setting the gamma coeffs corresponding to
%% specified boundary values (as opposed to boundary conditions
%% that involve derivatives that are not included in this example).

gamma=zeros(p.Nvertices,1);
for iv=1:p.Nbndry_vertices
  x=V(bndry_vertices(iv),1);
  y=V(bndry_vertices(iv),2);
  gamma(p.Ninterior_vertices+iv)=analytic_soln(x,y);
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function ib=vrtx_basis_func(iv,it,T)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% given vertex number iv and triangle element number it, find which
%% basis function is nonzero.  Result is in the range [1,2,3].  See
%% comment near definition of T.
ib=find(T(it,:)==iv);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function joint_triangles=find_joint_triangles(iv,jv,p)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% return list of triangles which include both vertices iv and jv.
list1=p.Triangles_with_vrtx(iv,1:p.nTriangles_with_vrtx(iv));
list2=p.Triangles_with_vrtx(jv,1:p.nTriangles_with_vrtx(jv));
joint_triangles=zeros(10,1);
njoint_triangles=0;
for i=1:length(list1)
  indx=find(list1(i)==list2);
  if not(isempty(indx)) && list1(i)~=0
    njoint_triangles=njoint_triangles+1;
    joint_triangles(njoint_triangles)=list1(i);
  end
end
joint_triangles=joint_triangles(1:njoint_triangles);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [A,rhs]=calculate_matrix_eqn(a,b,c,T,p,gamma)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% calculate the A=alpha_{ij} and rhs=beta_i for the matrix equation

alpha=zeros(p.Nvertices,p.Nvertices);
for iv=1:p.Nvertices
  for jv=1:p.Nvertices
    %% integrate over triangles which include both vertices iv,jv: 
    joint_triangles=find_joint_triangles(iv,jv,p);
    if not(isempty(joint_triangles))
      for itt=1:length(joint_triangles)
        it=joint_triangles(itt);
        %% indices for base functions that are not zero on triangle
        %% element it and vertices iv,jv:
        iv_b=vrtx_basis_func(iv,it,T);
        jv_b=vrtx_basis_func(jv,it,T);
        alpha(iv,jv)=alpha(iv,jv)+p.area(it) ...
                 *(b(it,iv_b)*b(it,jv_b)+c(it,iv_b)*c(it,jv_b));
      end
    end
  end
end
A=alpha(1:p.Ninterior_vertices,1:p.Ninterior_vertices);

rhs=zeros(p.Ninterior_vertices,1);
for i=1:p.Ninterior_vertices
  for k=p.Ninterior_vertices+1:p.Nvertices
    rhs(i)=rhs(i)-alpha(i,k)*gamma(k);
  end
end

%% print the matrix equation to info file:
fprintf(p.fid,'** the A matrix:\n');
for i=1:p.Ninterior_vertices; 
  fprintf(p.fid,'%9.2g,',A(i,:)); 
  fprintf(p.fid,'\n'); 
end

% fprintf(p.fid,'** the alpha_ij matrix:\n');
% for i=1:p.Nvertices; 
%   fprintf(p.fid,'%9.2g,',alpha(i,:)); 
%   fprintf(p.fid,'\n'); 
% end

fprintf(p.fid,'** the rhs of the matrix equation:\n');
fprintf(p.fid,'%9.2g,',rhs(:)); fprintf(p.fid,'\n'); 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function gamma=solve_matrix_eqn_and_plot(A,rhs,V,gamma,p);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

gamma(1:p.Ninterior_vertices)=A\rhs;

%% print solution and error:
fprintf(p.fid,'** the solution at all vertices:\n');
fprintf(p.fid,'%9.2g,',gamma(:)); fprintf(p.fid,'\n'); 

fprintf(p.fid,'** the solution minus analytic solution at all vertices:\n');
my_error=gamma*0;
for iv=1:p.Nvertices
  my_error(iv)=gamma(iv)-analytic_soln(V(iv,1),V(iv,2));
end
fprintf(p.fid,'%9.2g,',my_error(:)); fprintf(p.fid,'\n'); 

%% contour finite element and exact solutions:
figure(1); subplot(1,2,2);
% interpolate to a grid:
[xq,yq]=meshgrid(min(V(:,1)):0.001:max(V(:,1)),min(V(:,2)):0.001:max(V(:,2)));
finite_element_soln=griddata(V(:,1),V(:,2),gamma,xq,yq);
exact_soln=analytic_soln(xq,yq);
exact_soln(isnan(finite_element_soln))=NaN;
[C,h]=contourf(xq,yq,finite_element_soln,'edgecolor','none');
colorbarf(C,h);
v=get(h,'LevelList');
hold on
[C,h]=contour(xq,yq,exact_soln,v,'k');
clabel(C,h);
xlim([0 0.7]);
ylim([0 0.5]);
box on
xlabel('x')
ylabel('y')
title('exact solution: black; finite element: color')

%% save as pdf:
set(gcf, 'PaperUnits', 'inches'); set(gcf, 'PaperSize', [10 4]);
set(gcf, 'PaperPosition', [0 0 10 4]); % [left, bottom, width, height];
saveas(gcf,sprintf('Figure-finite_element_example.pdf'));



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function f=analytic_soln(x,y);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% A solution to Laplace eqn (that is, a harmonic function) on
%% which the b.c specified here are based.  If the solution is 
%% linear in x,y, it should be exactly reproduced by the finite
%% element solution.  Otherwise errors are to be expected, which
%% should degrease with element size. 

%f=2+5.*x+7.*y; % linear

f=2+5.*x+7.*y+400.*x.*y; % slightly nonlinear

%f=exp(15.*x).*sin(15.*y); % any harmonic function may be used here.
