% Solve the problem -\Delta u + c*u = f with Dirichlet boundary conditions.
% The domain is the disk centered at the origin of radius 1.
% Its triangulation is read from the end of this file.

function main()

   c=0; % a parameter in the equation, see above

   white=0.99*[1, 1, 1];
   blue = [0, 129, 205]/256;
   green = [0, 200,  70]/256;
   
% load the triangulation from the end of the file.
   dummy_arg=0;
   P=get_points(dummy_arg);
   T=get_triangles(dummy_arg);

% find the number of points and the number of triangles
   [Np, k]=size(P);
   [Nt, k]=size(T);
   Nb = 30; % first Nb points in P are on the boundary

% plot the triangulation
   lw = 1.4;
   figure(1); clf; hold on; axis equal; axis off;
   for l=1:Nt
      i=T(l, 1); j=T(l, 2); k=T(l, 3);
      
% plot the three edges in each triangle
      plot([P(i, 1), P(j, 1)], [P(i, 2), P(j, 2)], 'linewidth', lw, 'color', blue)
      plot([P(i, 1), P(k, 1)], [P(i, 2), P(k, 2)], 'linewidth', lw, 'color', blue)
      plot([P(j, 1), P(k, 1)], [P(j, 2), P(k, 2)], 'linewidth', lw, 'color', blue)
   end

% a hack to deal with bounding box issues
   s=1.05;
   plot(-s, -s, '*', 'color', white); plot(s, s, '*', 'color', white); 

% save as eps and svg (needs the plot2svg function)
   saveas(gcf, 'triangulation.eps', 'psc2') 
%  plot2svg('Finite_element_triangulation.svg');
   
% right-hand side
%  f=inline ('5-x.^2-y.^2', 'x', 'y');
   f=inline ('4+0*x.^2+0*y.^2', 'x', 'y');

% the values of f at the nodes in the triangulation
   F = f(P(:, 1), P(:, 2));
   RHS = 0*F; RHS = RHS((Nb+1):Np); % will solve A*U=RHS

% an empty sparse matrix of size Np by Np
   A = sparse(Np-Nb, Np-Nb);
   
% iterate through triangles 
   for l=1:Nt

%     fill in the matrix
      i=T(l, 1); j=T(l, 2); k=T(l, 3);

      [pii, pij, pik, pjj, pjk, pkk, gii, gij, gik, gjj, gjk, gkk] = calc_elems (i, j, k, P(i, :), P(j, :), P(k, :));

%     One has to consider the cases when a given vertex is on the boundary, or not.
%     If yes, it can't be included in the matrix
      if i > Nb
	 A(i-Nb, i-Nb) = A(i-Nb, i-Nb) + gii + c*pii;
      end

      if i > Nb & j > Nb
	 A(i-Nb, j-Nb) = A(i-Nb, j-Nb) + gij + c*pij;
	 A(j-Nb, i-Nb) = A(i-Nb, j-Nb);
      end


      if i > Nb & k > Nb
	 A(i-Nb, k-Nb) = A(i-Nb, k-Nb) + gik + c*pik;
	 A(k-Nb, i-Nb) = A(i-Nb, k-Nb);
      end

      if j > Nb
	 A(j-Nb, j-Nb) = A(j-Nb, j-Nb) + gjj + c*pjj;
      end

      if j > Nb & k > Nb
	 A(j-Nb, k-Nb) = A(j-Nb, k-Nb) + gjk + c*pjk;
	 A(k-Nb, j-Nb) = A(j-Nb, k-Nb);
      end

      if k > Nb
	 A(k-Nb, k-Nb) = A(k-Nb, k-Nb) + gkk + c*pkk;
      end

%     add the appropriate contributions to the right-hand terms
      if i > Nb 
	 RHS(i-Nb) = RHS(i-Nb) + F(i)*pii + F(j)*pij + F(k)*pik;
      end

      if j > Nb
	 RHS(j-Nb) = RHS(j-Nb) + F(i)*pij + F(j)*pjj + F(k)*pjk;
      end

      if k > Nb
	 RHS(k-Nb) = RHS(k-Nb) + F(i)*pik + F(j)*pjk + F(k)*pkk;
      end
      
   end

%  plot the sparse matrix (more exactly, its sign, then it is easier to see the pattern)
   figure(6); imagesc(sign(abs(A))); colormap (1-gray); axis ij; axis equal; axis off;

%  save as eps and svg (needs the plot2svg function)
   saveas(gcf, 'sparse_dirichlet.eps', 'psc2')
%  plot2svg('Finite_element_sparse_matrix.svg');

%  calculate U, then add zeros for the points on the boundary
   U_calc = A\RHS;
   U_calc = [0*(1:Nb) U_calc']';

%  exact solution
   u = inline('1-x.^2-y.^2', 'x', 'y');
   U_exact = u(P(:, 1), P(:, 2));

%  plot the computed solution
   figure(2); clf; 
   plot_solution(T, P, U_calc, lw, green)

%  save as eps and svg (needs the plot2svg function)
   saveas(gcf, 'computed_dirichlet.eps', 'psc2')
%  plot2svg('Finite_element_solution.svg');
   
%   plot the exact solution
%   figure(3); clf; 
%   plot_solution(T, P, U_exact)
%   title('Exact Dirichlet solution')
%   saveas(gcf, 'exact_dirichlet.eps', 'psc2')
   
   disp(sprintf('Error between exact solution and computed solution is %0.9g', max(abs(U_exact-U_calc))))
   
% given the l-th triangle, calculate the integrals
% pij = int_K(lambda_i*lambda_j) and  gij = int_K(grad lambda_i dot grad lambda_j)
function [pii, pij, pik, pjj, pjk, pkk, gii, gij, gik, gjj, gjk, gkk] = calc_elems (i, j, k, P1, P2, P3)

%  extract the x and y coordinates of the points
   x1 = P1(1); y1 = P1(2); 
   x2 = P2(1); y2 = P2(2); 
   x3 = P3(1); y3 = P3(2);

%  triangle area
   AT = abs((x2-x1)*(y3-y1)-(x3-x1)*(y2-y1))/2;

%  the integrals of lambda_i*lambda_j
   pii = AT/6;  pij = AT/12; pik = AT/12;
                pjj = AT/6;  pjk = AT/12;
		             pkk = AT/6;

			     
%  the integrals of grad lambda_i dot grad lambda_j    			     
   gii = ((x2-x3)*(x2-x3)+(y2-y3)*(y2-y3))/(4*AT);
   gjj = ((x1-x3)*(x1-x3)+(y1-y3)*(y1-y3))/(4*AT);
   gkk = ((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2))/(4*AT);

   gij = -((x3-x1)*(x3-x2)+(y3-y1)*(y3-y2))/(4*AT);
   gik = -((x2-x1)*(x2-x3)+(y2-y1)*(y2-y3))/(4*AT);
   gjk = -((x1-x2)*(x1-x3)+(y1-y2)*(y1-y3))/(4*AT);


function plot_solution (T, P, U, lw, color, fs)

   [Np, k]=size(P);
   [Nt, k]=size(T);

% create a path in 3D that will trace the outline of the solution
   X=zeros(5*Nt, 1);
   Y=zeros(5*Nt, 1);
   Z=zeros(5*Nt, 1);

   
   for l=1:Nt
      i=T(l, 1); j=T(l, 2); k=T(l, 3);

      m = 5*l-4;
      X(m)   = P(i, 1); Y(m)   = P(i, 2); Z(m)   = U(i);
      X(m+1) = P(j, 1); Y(m+1) = P(j, 2); Z(m+1) = U(j);
      X(m+2) = P(k, 1); Y(m+2) = P(k, 2); Z(m+2) = U(k);
      X(m+3) = P(i, 1); Y(m+3) = P(i, 2); Z(m+3) = U(i);
      X(m+4) = NaN;     Y(m+4) = NaN; Z(m+4) = NaN;
   end

%  plot the solution
   plot3(X, Y, Z,   'linewidth', lw, 'color', color)
   set(gca, 'fontsize', 15)

   

function P=get_points(dummy_arg)

   P=[1 0
0.97799999999999998 0.20799999999999999
0.91400000000000003 0.40699999999999997
0.80900000000000005 0.58799999999999997
0.66900000000000004 0.74299999999999999
0.5 0.86599999999999999
0.309 0.95099999999999996
0.105 0.995
-0.105 0.995
-0.309 0.95099999999999996
-0.5 0.86599999999999999
-0.66900000000000004 0.74299999999999999
-0.80900000000000005 0.58799999999999997
-0.91400000000000003 0.40699999999999997
-0.97799999999999998 0.20799999999999999
-1 5.6700000000000004e-16
-0.97799999999999998 -0.20799999999999999
-0.91400000000000003 -0.40699999999999997
-0.80900000000000005 -0.58799999999999997
-0.66900000000000004 -0.74299999999999999
-0.5 -0.86599999999999999
-0.309 -0.95099999999999996
-0.105 -0.995
0.105 -0.995
0.309 -0.95099999999999996
0.5 -0.86599999999999999
0.66900000000000004 -0.74299999999999999
0.80900000000000005 -0.58799999999999997
0.91400000000000003 -0.40699999999999997
0.97799999999999998 -0.20799999999999999
2.6291902682773483e-17 -0.00065384615384583118
-0.29530436727611131 -0.40714990300538867
0.20468717553848803 -0.45950882973942581
0.5000339709144086 -0.052282439231331621
0.2044416534667938 0.45895712720185405
-0.37389783786652392 0.33573030516976354
-0.56520453690073769 -0.059175479864500904
0.55140920581395014 0.3176130751959379
-0.13138695201728712 0.62243041389833154
-0.1444609125190702 -0.68304604895205301
0.51896266786078527 -0.46675660322909635
-0.67677962425194449 0.22091656257348974
-0.61723328415727241 -0.35579831401388728
0.15396059620050356 -0.7270900369296075
0.73917109220890242 0.07757578859901848
0.70634528325694623 -0.23042511622333953
0.15390386456341718 0.72682700843038894
0.43656290869075642 0.60123684202225869
-0.43638925515227728 0.60099824488402309
-0.44644144213994863 -0.6148097863548887
-0.78986314125549695 -0.082937447632792677
0.48306021940714738 -0.66512339089274719
-0.71711649590963977 0.4137416136492385
0.7196608156574994 0.41521760024330068
-0.17239398182683685 0.81255391574260738
0.83205471084583948 -0.087400017493309945
0.3403176445451454 0.76427847186026787
-0.55164668399041139 0.34772453168234385
-0.58119856273351911 0.5140347159468972
-4.3977950621443274e-18 -0.83738246960490337
0.72531999275454517 -0.41850054828302335
-0.34213203269662834 -0.76835550876536474
-0.62483523734501689 -0.56238344018259578
0.56491731978950166 0.48814678592215621
0.38917253432269977 0.4180653756110998
0.23717545979432902 0.16912185462803336
-0.02596228168575496 0.28616935365733154
0.3938288333314971 0.24434831196145793
0.5426706476687696 0.13031124265795693
0.38005096874664096 0.07200907885348734
0.24522905811277632 -0.17819696266523521
0.22329825233727255 -0.0049527105560499022
-0.018125285606249857 -0.28382039345137922
-0.25412870038086932 -0.1265506189486573
-0.22001244253386174 0.12401796918826143
0.70594582936521699 0.24504049798687536
-0.28923932941907254 0.68786790036032985
-0.27673116940573905 0.49861546727758305
-0.19580411520445754 0.33991543187901008
-0.06480669781759768 0.45864775798987351
0.039811209497243899 0.59660089503236968
-0.012452577623367307 0.7574818502043279
0.36004239628137363 -0.5404259782287143
0.3198317338457855 -0.70312416491673124
0.2617529370159869 0.31866372248564073
0.11364902800796572 0.34001459209079515
-0.80612573707151369 -0.26251531242500947
-0.84358781498716018 0.088619865046719065
-0.44904385214632619 0.45802230787668152
0.35347968445219013 -0.076483406156372766
0.43529169673866364 -0.27201412098929839
0.33966312039894614 -0.38995864108359496
-0.44160050993872146 0.12473708976956265
0.22042251667383625 -0.60109068660701248
0.01773990267719644 -0.55195139956951267
0.57741967278428696 0.63172121266253345
0.83616296331029605 0.18384745602550973
-0.42194773652236522 0.74460017792571787
-0.15191200145175998 -0.50793571011947636
-0.28468886832575607 -0.60129916940148143
-0.15211412155416426 -0.35582508357904519
0.091328532274082594 0.85396134404515744
0.25268219760086502 -0.82436443072803334
-0.56522176794370282 0.65147861937265739
-0.17960003622540516 -0.84596380431778762
0.64262958611925292 -0.57845575520448644
-0.82977670051153962 0.27012165242582176
-0.75613161550480157 -0.43637469407737101
0.84682787874323973 -0.27560544844003693
0.66460522715024917 -0.070226896389436425
0.27135187110185849 0.61028308052190072
-0.5291687129171514 -0.72847571124389077
-0.46589500975993081 -0.44182174455360729
-0.44717423342713147 -0.23275286852251509
-0.91658877313709719 -0.096341120235654235
0.071238170548373975 0.15038100559175938
-0.07537518167909911 0.14811510831848701
0.14411045609387388 -0.30720135657166003
0.046251548820153726 -0.41167351754462878
0.080240926729569728 -0.1479533766093776
-0.090696855217268765 -0.13701177299039563
-0.14601986575596898 -0.0018723983808613188
-0.65587042533541751 -0.19614678101972682
-0.68222829844326438 0.056477001652015607
0.5820334614227094 -0.32431245274666742
0.56939947784766054 -0.17974699975143729
-0.33054333836607347 0.20546019299713006
-0.54252415463210835 0.21357850028719924
-0.397822414790862 -0.038177513011169929
-0.0015128273166917471 -0.6940216317375888
0.12814058587287039 -0.86580177941165415
-0.30513569295646459 -0.26238649633237193];

function T = get_triangles (dummy_arg)
   T=[51 87 123
3 54 76
24 60 23
20 63 19
21 112 20
36 89 58
105 22 23
37 123 114
108 18 19
20 112 63
44 94 95
25 103 131
30 109 29
115 88 16
62 50 112
88 15 16
14 15 107
53 13 14
27 106 52
47 81 111
97 1 2
17 115 16
13 53 59
118 73 119
104 98 11
75 93 129
52 84 26
69 110 45
27 52 26
34 90 91
29 109 61
28 106 27
56 30 1
24 131 60
33 94 83
104 11 12
79 80 78
38 68 69
95 40 130
1 45 56
14 107 53
48 6 57
7 8 102
82 102 9
68 65 85
9 10 55
6 96 5
96 6 48
87 51 115
65 48 111
101 74 132
4 54 3
76 97 3
117 116 67
38 65 68
100 50 62
124 93 128
90 70 72
126 91 125
29 61 28
82 9 55
47 82 81
6 7 57
48 65 64
77 49 78
49 89 78
21 22 62
50 100 113
17 18 87
37 93 124
83 52 41
52 83 84
58 53 42
53 58 59
5 96 4
76 54 38
11 98 10
39 81 82
110 56 45
109 56 46
7 47 57
111 57 47
128 58 42
36 78 89
49 104 59
59 104 13
130 44 95
105 60 40
125 61 46
106 61 41
105 62 22
100 62 40
113 63 50
108 63 43
38 54 64
64 54 4
85 86 66
38 64 65
92 33 83
68 66 70
80 39 78
74 121 122
85 65 35
68 70 69
110 69 34
76 69 45
72 71 90
34 69 70
72 70 66
118 120 73
116 72 66
120 72 31
101 121 74
100 32 113
117 31 116
114 43 113
75 117 79
93 75 127
69 76 38
97 76 45
39 55 77
77 55 10
79 78 36
39 77 78
79 36 127
79 67 80
86 80 67
35 111 81
39 80 81
81 80 86
102 82 47
39 82 55
92 41 91
83 94 84
44 103 84
84 103 26
86 85 35
66 68 85
86 35 81
86 67 116
108 87 18
123 87 43
124 51 37
107 88 42
49 59 89
58 89 59
91 90 71
70 90 34
71 118 92
91 41 125
41 92 83
71 92 91
36 58 128
129 93 37
94 33 95
44 84 94
119 33 118
119 99 95
48 64 96
4 96 64
1 97 45
3 97 2
49 77 98
10 98 77
40 95 99
119 101 99
40 99 100
32 100 99
101 73 121
32 99 101
7 102 47
9 102 8
44 60 131
26 103 25
98 104 49
13 104 12
60 105 23
62 105 40
61 106 28
52 106 41
88 107 15
53 107 42
63 108 19
87 108 43
56 109 30
61 109 46
126 110 34
56 110 46
65 111 35
57 111 48
62 112 21
63 112 50
114 32 132
63 113 43
32 114 113
129 114 74
87 115 17
88 115 51
72 116 31
86 116 66
79 117 67
117 75 122
92 118 33
120 118 71
101 119 73
33 119 95
72 120 71
121 120 31
120 121 73
122 121 31
117 122 31
122 75 129
114 123 43
51 123 37
42 88 124
51 124 88
61 125 41
46 110 126
91 126 34
46 126 125
79 127 75
128 127 36
124 128 42
128 93 127
114 129 37
129 74 122
40 60 130
44 130 60
25 131 24
44 131 103
101 132 32
74 114 132];