function Q=my_quad2genf(f,p,t,R1,h)
Nnodes=size(p,2);
Ntri=size(t,2);

% Initierar värdet av integralen innan 
% summeringen över trianglar börjat

Q=0;

%Loopar över antalet trianglar

for n=1:Ntri
    
    % Definierar noderna i,j och k till triangel
    % Kn som hörnen med nummer t(1,n), t(2,n) resp.
    % t(3,n)
    
    i=t(1,n);
    j=t(2,n);
    k=t(3,n);
    
    % Definierar koordinaterna (x,y) för de tre noderna
    % till triangel Kn
    
    xnodi=p(1,i);
    xnodj=p(1,j);
    xnodk=p(1,k);
    
    ynodi=p(2,i);
    ynodj=p(2,j);
    ynodk=p(2,k);
    
    % Beräknar koordinaten för mittpunkten på Kn
    
    xmitt=(xnodi+xnodj+xnodk)/3;
    ymitt=(ynodi+ynodj+ynodk)/3;
    
    % Beräknar arean A av triangel Kn med noder i,j,k
    % som beloppet av kryssprodukten av två vektorer 
    % som spänner upp Kn (se Linjär Algebra)
    
    A=abs((xnodj-xnodi)*(ynodk-ynodi)-(xnodk-xnodi)*(ynodj-ynodi))/2;
    
    % Beräknar integralen med mittpunktskvadratur
    
%    Q=Q+feval(f,xmitt,ymitt)*A;
     Q=Q+f(xmitt,ymitt)*A;
    
end

% En plot av funktionsytan på cirkelskiva med radie R1

r=0:h:R1;
theta=0:h:2*pi+.1;
[R,Theta]=meshgrid(r,theta);
X1=R.*cos(Theta);
X2=R.*sin(Theta);
%surf(X1,X2,feval(f,X1,X2))
surf(X1,X2,f(X1,X2))
