%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%             Risoluzione di un problema di Poisson omogeneo
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all
clc

f = @(x, y) 1;
sol = @(x,y) -(1/4)*(x.^2+y.^2)+1/4; % soluzione esatta
grad = @(x,y) [-(1/2)*x -(1/2)*y];   % gradiente soluzione esatta
h = 5*logspace(-1,-2,10);


for k = 1:length(h)
    
    [ p , t, I, B ] = triangular_mesh('circle', h(k));
    
    
    % Calcolo l'area dei triangoli, dei baricentri e dei punti medi dei lati
    
    for i = 1:length(t(:,1))
        
        % Area
        
        mat=[p(t(i,1),1),p(t(i,1),2),1;
            p(t(i,2),1),p(t(i,2),2),1;
            p(t(i,3),1),p(t(i,3),2),1 ];
        ar(i,1)=(1/2)*abs(det(mat));
        
  
        verts = t(i,:);
        
        % Baricentri
        x_bar(i,1) = mean( p(verts,1) );
        y_bar(i,1) = mean( p(verts,2) );
        
        % Vertici
        x_tr(:,i) = p(verts,1);
        y_tr(:,i) = p(verts,2);
        
        % Punti medi
        x_midedge(i, :) = [mean([x_tr(1,i) x_tr(2,i)]) mean([x_tr(3,i) x_tr(2,i)]) mean([x_tr(1,i) x_tr(3,i)]) ];
        y_midedge(i, :) = [mean([y_tr(1,i) y_tr(2,i)]) mean([y_tr(3,i) y_tr(2,i)])  mean([y_tr(1,i) y_tr(3,i)]) ];
 
        
  
    end
    
 
    
    % Assemblaggio di A e F
    
    A = sparse(length(p), length(p));
    F = sparse(length(t(:,1)),1);
    
    extract = [1 2 3 1 2];
    
    for tt = 1:length(t)
        for iloc = 1:3
            
            iglob = t(tt,iloc);
            F(iglob) = F(iglob) + ar(tt)*f(x_bar(tt), y_bar(tt))/3;
            
            for jloc=1:3
                jglob=t(tt,jloc);
                
                i1 = t(tt,extract(iloc+1));
                i2 = t(tt,extract(iloc+2));
                j1 = t(tt,extract(jloc+1));
                j2 = t(tt,extract(jloc+2));
                
                l1=[p(i2,1)-p(i1,1) p(i2,2)-p(i1,2)];
                l2=[p(j2,1)-p(j1,1) p(j2,2)-p(j1,2)];
                
                A(iglob,jglob) = A(iglob,jglob)+(l1*l2')/(4*(ar(tt)));
            end
        end
    end
     
    A = A(I,I);
    F = F(I);
    
    
    % Risoluzione del sistema lineare
    u_h = zeros(length(p(:,1)),1);
    u_h(I)=A\F;
 
    %----------------- Grafici della soluzione calcolata e dell'interpolata ------------------
    % Grafico della soluzione calcolata
    figure(1)
    subplot(1,2,1)
    trimesh(t,p(:,1),p(:,2),u_h);
    title('Grafico della u_{FEM}')
    
    % Grafico dell'interpolata
    for i = 1:length(p)
        u_ex(i,1) = sol(p(i,1),p(i,2));
    end
    subplot(1,2,2)
    trimesh(t, p(:,1), p(:,2),u_ex);
    title('Grafico dell''interpolata della u_{esatta}')
    
    
    %------------------------------------ Calcolo degli errori ------------------------------------
    % Calcolo dell'errore in norma infinito
    errinf(k) = norm(u_h-u_ex,inf);
    
    
    % Calcolo dell'errore in norma L^2: formula del baricentro
    u_ex_bar = sol(x_bar, y_bar);
    
    for i = 1:length(t(:,1))
        uh_bar(i,1) = (1/3)*(u_h(t(i,1))+u_h(t(i,2))+u_h(t(i,3)));
    end
    
    errL2(k) = (sum((u_ex_bar-uh_bar).^2, 2)'*ar)^(1/2);
    
    
    % Calcolo dell'errore in norma L^2: formula del punto medio
     errL2_medio(k) = 0;
    for i = 1:length(t(:,1))
 
        u1 = u_h(t(i,1));
        u2 = u_h(t(i,2));
        u3 = u_h(t(i,3));
 
        errL2_medio(k) = errL2_medio(k) + ar(i)/3*((sol(x_midedge(i,1),y_midedge(i,1))-(u2+u1)*0.5).^2+(sol(x_midedge(i,2),y_midedge(i,2))-0.5*(u3+u2)).^2 + (sol(x_midedge(i,3),y_midedge(i,3))-(u3+u1)*0.5)^2);
        
    end
      errL2_medio(k) = sqrt(errL2_medio(k));
         
      % Calcolo dell'errore in norma H^1: formula del baricentro
      gradex = grad(x_bar, y_bar);
      
      for j = 1:length(t(:,1))
          
          rot = [0 -1;1 0];
          
          i1 = t(j,1);
          i2 = t(j,2);
          i3 = t(j,3);
          
          l1 = [p(i3,1)-p(i2,1) p(i3,2)-p(i2,2)];
          l2 = [p(i1,1)-p(i3,1) p(i1,2)-p(i3,2)];
          l3 = [p(i2,1)-p(i1,1) p(i2,2)-p(i1,2)];
          
          gradbar(j,:) = (u_h(i1)/(2*ar(j))).*rot*l1' + (u_h(i2)/(2*ar(j))).*rot*l2' + (u_h(i3)/(2*ar(j))).*rot*l3';
          
      end
      
      
       errH1_baric(k)=(((u_ex_bar-uh_bar).^2 +  sum((gradbar-gradex).*(gradbar-gradex),2))'*ar)^(1/2);

      
    % Calcolo dell'errore in norma H^1: fomula del punto medio 
    % (uso ancora gradbar perch� il gradiente di una funzione lineare � costante)
       
    for i = 1:length(t(:,1))
 
        x = grad(x_midedge(i,1), y_midedge(i,1));
        y = grad(x_midedge(i,2), y_midedge(i,2));
        z = grad(x_midedge(i,3), y_midedge(i,3));
        
        xx = (x(1) - gradbar(i,1))^2 + (x(2) - gradbar(i,2))^2;
        yy = (y(1) - gradbar(i,1))^2 + (y(2) - gradbar(i,2))^2;
        zz = (z(1) - gradbar(i,1))^2 + (z(2) - gradbar(i,2))^2;
  
        integrale1(i) = ar(i)*((1/3)*(xx + yy + zz));
    end
    
    
    errH11(k) = ((sum(integrale1) +   errL2_medio(k)^2))^(1/2);
    
     
end

figure(3)
loglog(h,errinf,'-*')
hold on
loglog(h,.1*h.^2,'-')
title('Grafico dell''errore in norma infinito')
legend('Errore in norma infinito', 'h^2')
xlabel('h')
ylabel('Errore')

figure(4)
loglog(h,errL2_medio,'-<')
hold on
loglog(h,h.^2,'-')
title('Grafico dell''errore in norma L^2')
legend('Errore calcolato con la formula del punto medio dei lati','h^2')
xlabel('h')
ylabel('Errore')


figure(5)

loglog(h,errH11,'-<')
hold on
loglog(h,h,'-')
title('Grafico dell''errore in norma H^1')
legend('Errore calcolato con la formula del punto medio dei lati','h')
xlabel('h')
ylabel('Errore')
