% script file mean_tris.m
% Exercise 1

clc; clear all; close all;

n = 50; % number of points along one boundary


x = 2*rand(ceil(0.25*n^2),1)-1;
y = 2*rand(ceil(0.25*n^2),1)-1;

% add border points
% unit circle:
phi = linspace(0,2*pi,n)';
% delete end point not to be twice
phi(end) = [];
bd_x = cos(phi);
bd_y = sin(phi);

n = n-1; % last point removed

% boundary value:
ud = 0.2*sin(5*phi);



% domain: out of [-1,1]^2
% leave only those which are inside of the circle
OUT = find(x.^2+y.^2>=1);
x(OUT) = [];
y(OUT) = [];

Nodes.x = [x;bd_x];
Nodes.y = [y;bd_y];

N = length(Nodes.x); % total number of nodes

TRI = delaunay(Nodes.x,Nodes.y);

figure(1);
clf
P=trimesh(TRI,Nodes.x,Nodes.y);
set(P,'color','k');
hold on
%plot(x,y,'b*')
axis equal;

% the rhs: zero at all inner points,
% ud values at the boundary points
RHS = zeros(N,1);
RHS(N-n+1:end) = [ud];

A = sysmat_min_surf(N,n,TRI,Nodes);

U = A\RHS;

figure(2)
trisurf(TRI,Nodes.x,Nodes.y,U);
xlabel('x')
ylabel('y')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Exercise 2
% now do the ordering:
P = symrcm(A);
A1 = A(P,P);

% rhs reordered:
RHS1 = RHS(P);

U1 = A1\RHS1; % U1 = U(P) by construction

% visualise U1 using a reordered TRIs
% find the inverse permutation vector
Pinv = zeros(N,1);
for ctr=1:N
  Pinv(ctr) = find(P==ctr);
end

% create the new triangulation array
trin = size(TRI,1); % number of triangles
TRI1 = zeros(trin,3);
for ctr=1:trin
  TRI1(ctr,1) = Pinv(TRI(ctr,1));
  TRI1(ctr,2) = Pinv(TRI(ctr,2));
  TRI1(ctr,3) = Pinv(TRI(ctr,3));    
end

% the new nodes
Nodes1.x = Nodes.x(P);
Nodes1.y = Nodes.y(P);

figure(3)
trisurf(TRI1,Nodes1.x,Nodes1.y,U1);
xlabel('x')
ylabel('y')
