clear

h=0.16;


% we need interpolation first
%-----------------------------------------------------------


my_d = @(P) ( max(abs( P(:,1)-0.5), abs(P(:,2) - 0.5)) -0.5 );% .* ( sqrt((P(:,1)-x_c).^2 + (P(:,2)-y_c).^2) - rad_c  );
[P,T]=distmesh2d(my_d,@huniform,h,[-1 -1; 2 2;],[0 0; 1 1; 0 1; 1 0]);
%u = @(x,y) x.^4 + y.^4;
u=@(x,y) exp(x.^2 + y.^2);

% P1 is easy, just compute the value of u on the nodes of the mesh
u_vert = u(P(:,1),P(:,2));

% P0 requires more work: in detail, we need to compute baricenter of each triangle and create a matrix where
% each column contains the x and y coordinates of the three vertices of each triangle. For later purposes, we
% also need the value of the three midpoints of the edges
nb_t = size(T,1);
x_bar = zeros(nb_t,1);
y_bar = zeros(nb_t,1);
x_tr = zeros(3,nb_t);
y_tr = zeros(3,nb_t);
x_midedge = zeros(3,nb_t);
y_midedge = zeros(3,nb_t);

% loop over triangles
for t=1:nb_t
    % the baricenter is the mean of the vertices of a triangle. I access the vertices of a triangle as T(t,:)
    the_verts = T(t,:);
    x_bar(t) = mean( P(the_verts,1) );
    y_bar(t) = mean( P(the_verts,2) );
    x_tr(:,t) = P(the_verts,1);
    y_tr(:,t) = P(the_verts,2);
    x_midedge(:,t) = [mean([x_tr(1,t) x_tr(2,t)]); mean([x_tr(1,t) x_tr(3,t)]); mean([x_tr(3,t) x_tr(2,t)])];
    y_midedge(:,t) = [mean([y_tr(1,t) y_tr(2,t)]); mean([y_tr(1,t) y_tr(3,t)]); mean([y_tr(3,t) y_tr(2,t)])];
end

% now evaluate u in the baricenter 
u_bar = u(x_bar,y_bar);






% quadrature over triangles.
%-----------------------------------------------------------

% we loop over triangles and add the contribution of each triangle to the overall integral

% Init the values of the integrals
Qh_vert = 0;
Qh_bar=0;
Qh_midp = 0;


% loop over triangles
for t=1:nb_t
    % recover the vertices of the triangle and the values of u on each vertex, in the baricenter and in the
    % midpoints
    the_verts = T(t,:);
    u_vert_loc = u_vert(the_verts);
    u_bar_loc = u_bar(t);
    
    u_midp_loc = u(x_midedge(:,t),y_midedge(:,t));
    
    % we need the area of the triangle. Here we use the formula with tensor product of vectors: let A = (P2-P1) x (P3-P1) / 2
    % then Area = length(A) = 3rd component of A (the other two are zero because P2-P1 and P3-P1 are in the xy
    % plane)
    % we could also use Heron's formula, or determinant formula, or built-in Matlab functions
    P1 = [P(the_verts(1),:) 0]; % we need 3D vectors to take tensor products
    P2 = [P(the_verts(2),:) 0];
    P3 = [P(the_verts(3),:) 0];
    
    A = cross(P2-P1,P3-P1)/2; % cross is a built-in matlab function that computes tensor product
    Area = A(3);
    % update quadrature value
    
    Qh_vert = Qh_vert + Area/3*sum(u_vert_loc);
    Qh_bar = Qh_bar + Area*u_bar_loc;
    Qh_midp = Qh_midp + Area/3*sum(u_midp_loc);
    
end

disp(Qh_vert)
disp(Qh_bar)
disp(Qh_midp)
