%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                                                                      %
%%                                                                                      %
%% Syfte:    Löser, med cG(1)-metoden, differentialekvationen                           %
%%                                                                                      %
%%              -div(a(x,y)*grad(u(x,y))) + c(x,y)*u(x,y) = f(x,y),    i Omega,         %
%%                                                                                      %
%%           med Robinrandvillkor                                                       %
%%                                                                                      %
%%              -a*Du/Dn = gamma*(u - g_D) + g_N, på randen, Gamma, till Omega,         %
%%                                                                                      %
%%           och "plottar" lösningen.                                                   %
%%                                                                                      %
%%           (Du/Dn betecknar skalärprodukten mellan grad(u) och den                    %
%%            utåtriktade enhetsnormalen, n, d.v.s. riktningsderivatan                  %
%%            av u längs riktningen vinkelrätt ut ur Omega.)                            %
%%                                                                                      %
%%                                                                                      %
%% Variabler:                                                                           %  
%%                                                                                      %
%%           p: "Point matrix"         (Matris med dimension 2 x nnodes)                %
%%           e: "Edge matrix"          (Matris med dimension 7 x nedges)                %
%%           t: "Triangle matrix"      (Matris med dimension 4 x ntri)                  %
%%                                                                                      %
%%              p, e, och t, som innehåller information om trianguleringen              %
%%                        av Omega, exporteras från GUI:n 'pdetool'.                    %
%%                                                                                      %
%%              Viktigast att känna till:                                               %
%%                                                                                      %
%%              'p' innehåller alla nodkoordinater:                                     %
%%              p(1,i) är x-koordinaten för nod nummer i.                               %
%%              p(2,i) är y-koordinaten för nod nummer i.                               %
%%              (1 <= i <= nnodes)                                                      %
%%                                                                                      %
%%              'e' innehåller information om de triangelsidor som ligger               %
%%                  på ett randsegment:                                                 %
%%              e(1,l) är nodnumret för sida l:s "startpunkt".                          %
%%              e(2,l) är nodnumret för sida l:s "slutpunkt".                           %
%%              e(5,l) är numret på det randsegment varpå sida l ligger.                %
%%              (1 <= l <= nedges)                                                      %
%%                                                                                      %
%%              't' ger information om vilka tre noder som är hörnnoder i               %
%%                  en given triangel:                                                  %
%%              t(1,n) är nodnumret för "första" hörnet i triangel nummer n.            %
%%              t(2,n) är nodnumret för "andra" hörnet i triangel nummer n.             %
%%              t(3,n) är nodnumret för "tredje" hörnet i triangel nummer n.            %
%%              1 -> 2 -> 3 moturs riktning.                                            %
%%              (1 <= n <= ntri)                                                        %
%%                                                                                      %
%%           A:  styvhetsmatris        (Matris med dimension nnodes x nnodes)           %
%%          Mc:  massmatris            (Matris med dimension nnodes x nnodes)           %
%%           R:  randmatris            (Matris med dimension nnodes x nnodes)           %
%%                                                                                      %
%%           b:  lastvektor            (Vektor med dimension nnodes x 1)                %
%%          rv:  randvektor            (Vektor med dimension nnodes x 1)                %
%%                                                                                      %
%%           U:  lösning               (Vektor med dimension nnodes x 1)                %
%%                                                                                      %
%%                                                                                      %
%% Funktioner:             a_2D.m  definierar funktionen       a(x,y),   i Omega        %
%%                         c_2D.m  definierar funktionen       c(x,y),   i Omega        %
%%                         f_2D.m  definierar funktionen       f(x,y),   i Omega        %
%%                                                                                      %
%%                     gamma_2D.m  definierar funktionen   gamma(x,y),  på Gamma        %
%%                       g_D_2D.m  definierar funktionen     g_D(x,y),  på Gamma        %
%%                       g_N_2D.m  definierar funktionen     g_N(x,y),  på Gamma        %
%%                                                                                      %
%%    MyFirst2DPoissonAssembler.m  beräknar A, Mc, R, b, och rv                         %
%%                                                                                      %
%%                                                                                      %
%%                                                                                      %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                                                                      %
%% Exempel:  Du vill lösa ekvationen                                                    %
%%                                                                                      %
%%             -div(grad(u(x,y))) = 1, i Omega,                                         %
%%                                                                                      %
%%           där Omega är rektangeln:  -1 < x < 1, -0.8 < y < 0.8,                      %
%%                                                                                      %
%%           med randvillkor                                                            %
%%                                                                                      %
%%              u = 0.84 - (x^2)/4,                      -1 < x < 1, y = -0.8,          %
%%              u = 0.75 - (y^2)/4,                      x = 1, -0.8 < y < 0.8,         %
%%             -Du/Dy = -Du/Dn = 10/21*(u + (x^2)/4),    -1 < x < 1, y = 0.8,           %
%%              Du/Dx = -Du/Dn = 2/3*(u + (y^2)/4),      x = -1, -0.8 < y < 0.8.        %
%%                                                                                      %
%%           Börja med att rita och triangulera Omega i GUI:n 'pdetool',                %
%%           varefter du exporterar matriserna p, e, och t.                             %
%%                                                                                      %
%%           Definiera därefter funktionerna 'a', 'c' och 'f':                          %
%%                                                                                      %
%%             function diffusion = a_2D(x,y,delomradesnummer)                          %
%%             diffusion = 1;                                                           %
%%                                                                                      %
%%             function reaktion = c_2D(x,y,delomradesnummer)                           %
%%             reaktion = 0;                                                            %
%%                                                                                      %
%%             function source = f_2D(x,y,delomradesnummer)                             %
%%             source = 1;                                                              %
%%                                                                                      %
%%           Notera att vi i detta exempel endast har ett delområde: hela Omega.        %
%%           (I vissa tillämpningar, t.ex. då Omega består av flera material,           %
%%           kan varje material utgöra ett delområde och materialparametern a           %
%%           definieras av olika uttryck i olika delområden.) För att erhålla           %
%%           information om vilka delområden Omega är uppdelat i, samt deras nummer,    %
%%           väljer du, inifrån 'pdetool', "Boundary Mode" i menyn "Boundary",          %
%%           och markerar därefter "Show Subdomain Labels" i samma meny.                %
%%                                                                                      %
%%           Notera i 'pdetool' också vilket randsegmentnummer som hör                  %
%%           till respektive randvillkor. Denna information får du genom                %
%%           att välja "Boundary Mode" i menyn "Boundary", och därefter                 %
%%           markera "Show Edge Labels" i samma meny. Antag att de är                   %
%%           numrerade 1, 2, 3, 4, i den ordning som de skrivits ovan.                  %
%%                                                                                      %
%%           Definiera slutligen "randvillkorsfunktionerna" 'gamma', 'g_D', och 'g_N':  %
%%                                                                                      %
%%             function permeabilitet = gamma_2D(x,y,rand_nr)                           %
%%                                                                                      %
%%             if rand_nr == 1                                                          %
%%               permeabilitet = 10^7; % gamma(x,y) på randsegment 1.                   %
%%             elseif rand_nr == 2                                                      %
%%               permeabilitet = 10^7; % gamma(x,y) på randsegment 2.                   %
%%             elseif rand_nr == 3                                                      %
%%               permeabilitet = 10/21; % gamma(x,y) på randsegment 3.                  %
%%             elseif rand_nr == 4                                                      %
%%               permeabilitet = 2/3; % gamma(x,y) på randsegment 4.                    %
%%             end                                                                      %
%%                                                                                      %
%%                                                                                      %
%%             function omgivning = g_D_2D(x,y,rand_nr)                                 %
%%                                                                                      %
%%             if rand_nr == 1                                                          %
%%               omgivning = 0.84 - (x^2)/4; % g_D(x,y) på randsegment 1.               %
%%             elseif rand_nr == 2                                                      %
%%               omgivning = 0.75 - (y^2)/4; % g_D(x,y) på randsegment 2.               %
%%             elseif rand_nr == 3                                                      %
%%               omgivning = -(x^2)/4; % g_D(x,y) på randsegment 3.                     %
%%             elseif rand_nr == 4                                                      %
%%               omgivning = -(y^2)/4; % g_D(x,y) på randsegment 4.                     %
%%             end                                                                      %
%%                                                                                      %
%%                                                                                      %
%%             function flux = g_N_2D(x,y,rand_nr)                                      %
%%                                                                                      %
%%             if rand_nr == 1                                                          %
%%               flux = 0; % g_N(x,y) på randsegment 1.                                 %
%%             elseif rand_nr == 2                                                      %
%%               flux = 0; % g_N(x,y) på randsegment 2.                                 %
%%             elseif rand_nr == 3                                                      %
%%               flux = 0; % g_N(x,y) på randsegment 3.                                 %
%%             elseif rand_nr == 4                                                      %
%%               flux = 0; % g_N(x,y) på randsegment 4.                                 %
%%             end                                                                      %
%%                                                                                      %
%%                                                                                      %
%%           Du kan nu starta programmet med kommandot:                                 %
%%                                                                                      %
%%             >> MyFirst2DPoissonSolver 		 		 		%
%%                                                                                      %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                    %
% Assemblera matriser och vektorer.  %
%                                    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


[A, Mc, R, b, rv] = MyFirst2DPoissonAssembler(p, e, t);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                             %
% Lös ekvationssystemet.      %
%                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


VL = A + Mc + R;
hl = b + rv;

 U = VL\hl;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                    %
% Plotta lösningen.                  %
%                                    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  


pdesurf(p, t, U)
xlabel('x')
ylabel('y')
title('U(x,y)')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                                                         %
%%           Den exakta lösningen är (kontrollera!)                        %
%%                                                                         %
%%             u(x,y) = 1 - (x^2 + y^2)/4,                                 %
%%                                                                         %
%%           så om du vill kan du beräkna felet:                           %
%%                                                                         %
%%             >> x = p(1,:)'; % Kolonnvektor.                             %
%%             >> y = p(2,:)';                                             %
%%             >> u = 1 - (x.^2 + y.^2)/4; % Blir också kolonnvektor...    %
%%             >> fel = u - U; % ... så dimensionerna stämmer överens!     %
%%                                                                         %
%%           och därefter "plotta" det:                                    %
%%                                                                         %
%%             >> pdesurf(p, t, fel)                                       %
%%                                                                         %
%%           En approximation till max-normen av felet får du med          %
%%                                                                         %
%%             >> max(abs(fel))                                            %
%%                                                                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
