function [C,K,F] = FEM(eq,t,U)
%FEM Provede prostorovou diskretizaci lineárního problému metodou konečných prvků.
% Vrací matice C,K,F, stacionární problém je aproximován SLR: K*U=F a
% nestacionární převede na počáteční problém pro SODR1, který lze řešit
% theta metodou viz theta.mat. Jestliže je problém lineární, jsou i
% výsledná SLR, SODR1 lineární, jinak je třeba řešit jako nelineární
% problémy.
% Vstupy: eq - strukture p5edstavující data přímé úlohy vedení tepla, t -
% čas, ve kterém vyhodnocuji matice C, K a F a U - řešení nelineární úlohy
% v předchozí iteratici
% Výstupy: matice C, K a vektor F vyhodnocené v t a v U
% Počty uzlů, prvků a hran
if ~exist('t','var')
    t = 0;
end
if ~exist('U','var')
    U = zeros(size(eq.domain.p,2),1);
end
num_pts = size(eq.domain.p,2);
num_edge = max(eq.domain.e(5,:));
num_elem = size(eq.domain.t,2);
% Okrajové podmínky
dir_p = cell(size(eq.dir_edge,2),1); % Body, kde je předepsána Dirichletova OP
for i = 1:size(eq.dir_edge,2)
    dir_p{i} = unique([eq.domain.e(1,eq.domain.e(5,:)==eq.dir_edge(i)),eq.domain.e(2,eq.domain.e(5,:)==eq.dir_edge(i))]);
end
neu_e = cell(size(eq.neum_edge,2),1); % Hrany, kde je předepsána Neumannova OP
for i = 1:size(eq.neum_edge,2)
    neu_e{i} = vertcat(eq.domain.e(1,eq.domain.e(5,:)==eq.neum_edge(i)),eq.domain.e(2,eq.domain.e(5,:)==eq.neum_edge(i)));
end
% Inicializace globálních matic C, K a vektoru F
C = zeros(num_pts,num_pts);
K = zeros(num_pts,num_pts);
F = zeros(num_pts,1);
% Cyklus přes elementy
for i = 1:num_elem
    Te = eq.domain.t(1:3,i); % Prvek
    Pe = eq.domain.p(:,Te); % Souřadnice vrcholů prvku
    Se = [mean(Pe(1,:)),mean(Pe(2,:))]; % Souřadnice těžiště prvku
    de = (Pe(2,3)-Pe(2,1))*(Pe(1,2)-Pe(1,1))-(Pe(1,1)-Pe(1,3))*(Pe(2,1)-Pe(2,2));
    re = (Pe(2,2)-Pe(2,3))*(Pe(2,3)-Pe(2,1))+(Pe(1,3)-Pe(1,2))*(Pe(1,1)-Pe(1,3));
    se = (Pe(2,2)-Pe(2,3))*(Pe(2,1)-Pe(2,2))+(Pe(1,3)-Pe(1,2))*(Pe(1,2)-Pe(1,1));
    te = (Pe(2,3)-Pe(2,1))*(Pe(2,1)-Pe(2,2))+(Pe(1,1)-Pe(1,3))*(Pe(1,2)-Pe(1,1));
    Ke1 = (0.5/abs(de))*eq.p(Se(1),Se(2),t).*[-re-se,re,se;re,-re-te,te;se,te,-se-te];
    Ke2 = (abs(de)/6).*diag([eq.q(Pe(1,1),Pe(2,1),t,U(Te(1))),eq.q(Pe(1,2),Pe(2,2),t,U(Te(2))),eq.q(Pe(1,3),Pe(2,3),t,U(Te(3)))]);
    if ~eq.stationary
        Ce = (abs(de)/6).*diag([eq.c(Pe(1,1),Pe(2,1),t,U(Te(1))),eq.c(Pe(1,2),Pe(2,2),t,U(Te(2))),eq.c(Pe(1,3),Pe(2,3),t,U(Te(3)))]);
    end
    Ke = Ke1+Ke2; % Lokální matice K
    Fe = (abs(de)/6).*[eq.g(Pe(1,1),Pe(2,1),t,U(Te(1)));eq.g(Pe(1,2),Pe(2,2),t,U(Te(2)));eq.g(Pe(1,3),Pe(2,3),t,U(Te(3)))]; % Lokální vektor F
    % Přiřazení do globální matice C, K, F
    for j = 1:3
        jg = eq.domain.t(j,i);
        for k = 1:3
            kg = eq.domain.t(k,i);
            if ~eq.stationary
                C(jg,kg) = C(jg,kg) + Ce(j,k);
            end
            K(jg,kg) = K(jg,kg) + Ke(j,k);
        end
        F(jg) = F(jg) + Fe(j);
    end
end
% Cyklus přes hrany, kde je předepsána Neumannova OP
for i = eq.neum_edge
    alpha = eq.alpha{i};
    beta = eq.beta{i};
    I =find(eq.neum_edge == i);
    for j = 1:size(neu_e{I},2)
        Pe = eq.domain.p(:,neu_e{I}(:,j)); % Souřadnice prvku S
        de = sqrt((Pe(1,2)-Pe(1,1))^2+(Pe(2,2)-Pe(2,1))^2);
        Ks = 0.5*de.*[alpha(Pe(1,1),Pe(2,1),t,U(neu_e{I}(1,j))),0;0,alpha(Pe(1,2),Pe(2,2),t,U(neu_e{I}(2,j)))]; % Lokální matice K
        Fs = 0.5*de.*[beta(Pe(1,1),Pe(2,1),t,U(neu_e{I}(1,j)));beta(Pe(1,2),Pe(2,2),t,U(neu_e{I}(2,j)))]; % Lokální vektor F
        for k = 1:2
            kg = neu_e{I}(k,j);
            for l = 1:2
                lg = neu_e{I}(l,j);
                K(kg,lg) = K(kg,lg) + Ks(k,l);
            end
            F(kg) = F(kg) + Fs(k);
        end
    end
end
% Cyklus přes body, kde je předepsána Dirichletova OP
kappa = 1e10;
for i = 1:size(eq.dir_edge,2)
    for j = 1:size(dir_p{i},2)
        jg = dir_p{i}(j);
        Pe = eq.domain.p(:,jg);
        F(jg) = eq.f{eq.dir_edge(i)}(Pe(1),Pe(2),t)*K(jg,jg)*kappa;
        K(jg,jg) = K(jg,jg)*kappa;
    end
end
end

