PRILOG ZA STUDENTE ZA IZRADU SEMINARSKOG RADA 5.

a)

Konvektivno-difuziona jednačina (stacionarna) glasi:

U prilogu je program koji rješava ovu diferencijalnu jednačinu, na 2D domenu dimenzija Lx x

Ly.Strujno polje ⃗⃗ je poznato i definisano postavkom: vx=C∙x, vy= - C∙y. Shema Dekartove
(kartezijanske strukturisane) numeričke mreže i oznaka data je na slici Sl.1. Jednačina je
diskretizovana metodom konačnih (kontrolisanih) zapremina (CV, FV- metod). Dat je Matlab program
koji je napisan za ovu svrhu, u skladu sa oznakama i postavkom. Program sadrži komentare
koji objašnjavaju svaku cjelinu programa.

Predmetni asistent
Milan Šekularac

Sl.1. 2D domen i shema Dekartove strukturisane mreže konačnih zapremina sa oznakama

%% ######################################################################
%% PREDMET KME, VJEZBE.
%% Milan Sekularac, predmetni asistent
%% 11.CAS, STACIONARNA TRANSPORTNA JEDNACINA SA KONVEKCIJOM I DIFUZIJOM
%***********************************************************************
% PROGRAM RJESAVA JEDNACINU STACIONARNOG KONVEKTIVNO-DIFUZIONOG
% TRANSPORTA STRUJOM FLUIDA, SKALARNE FIZICKE VELICINE "fi".
%
% d/dx(ro*u*fi)+d/dy(ro*v*fi)=d/dx(L*dfi/dx)+d/dy(L*dfi/dy)+S (1)
%
% ili kao: div( ro*U*fi ) = div( L*grad(fi) ) + S
%
%***********************************************************************
% JEDNACINA JE DISKRETIZOVANA METODOM KONACNIH ZAPREMINA, DOMEN JE
% PRAVOUGAONI, DIMENZIJA Lx x Ly. ULAZNA GRANICA JE SJEVERNA (N),
% IZLAZNA JE ISTOCNA (E). GRANICNI USLOVI SU DIRICHLET-OVI ZA N,W,S
% GRANICE I NOJMANOV ZA IZLAZNU (E) GRANICU: grad(fi) x ne = 0.
%***********************************************************************
% STRUJNO POLJE JE POZNATO UNAPRIJED I ZADATO KAO: V(u,v) = f(x,y)
% PROJEKCIJE BRZINE SU: u=x, v=-y,
% (*Postavka je uzeta iz primjera u knjizi J.Ferziger, M.Peric
% "Computational methods for fluid dynamics").
%#######################################################################

close all % zatvori otvorene grafike ako ih ima

clear all % obrisi sve iz radne memorije
clc % obrisi sve sa ekrana komandnog prozora
%########################################################################
#
%% FIZICKI PARAMETRI I POSTAVKA _____
Lx = 3.;% % dimenzija domena u x-pravcu [m] Ly| |
Ly = 1; % % dimenzija domena u y-pravcu [m] |_Lx__|
% STRUJNO POLJE JE DATO POSTAVKOM: u = Cuv*x, v = -Cuv*y.
Cuv = 0.05; % Cuv je faktor skaliranja intenziteta projekcija brzine
%########################################################################
#
%% NUMERICKI PODACI
Ni = 30.; %9.; % broj vert. linija po kojima leze centri CV-a
Nj = 30.; %8.; % broj hor. linija po kojima leze centri CV-a
% NAPOMENA: KOMENTARI (U VEZI S MREZOM/INDEKSIMA) SE ODNOSE NA MREZU
% KOJA JE DATA UZ ZADATAK KAO ILUSTRACIJA: Ni x Nj = 9 x 7
%% REZULTUJUCI NUMERICKI KORACI NA MREZI KONACNIH ZAPREMINA
deltax = Lx/(Ni-2); % razmak izmedju 2 stran. susjednih CV-a u x pravcu
deltay = Ly/(Nj-2); % razmak izmedju 2 stran. susjednih CV-a u y pravcu
%% RASTOJANJA
dely = deltay; % razmak izmedju 2 centralna cvora susjednih CV-a
delx = deltax; % razmak izmedju 2 centralna cvora susjednih CV-a
%% Numericki parametri za iterativnu proceduru:
omega = 0.85; % relaksacioni Cuv kod SOR [/]
epsilon = 10.^-8.; % kriterijum za zaustavljanje iteriranja [/]
MaxIter = 5000; % Maximalni broj iteracija !? Prati uslov "Rez < epsilon" !!!
%########################################################################
#
%% FIZICKE OSOBINE SUSPTANCE (FLUIDA)
% KOEFICIJENT DIFUZIJE (GAMA):
L = 1; % koeficijent difuzije Npr.[W/mK]
% SPECIFICNA TOPLOTA:
c = 1.; % specificna toplota materijala [J/kg K]
ro = 1.2; % gustina materijala [kg / m^3]
ro = ro*c; % s obzirom na oznacavanje pri diskretizaciji!
%########################################################################
#
%% IZVORNI CLAN S
% Neka je izvor kontstantan za svaku CV-u i jednak Sc
Sc = 0.0; % izvorni clan - konstantni dio % Npr.[W/m^3]
%IZVORNI CLAN JE U OPSTEM SLUCAJU: S := Sc+Sp*(fi)p
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
~~~~~~
% (*FUSNOTA O IZVORNOM CLANU:
% DA JE JEDNACINA U ZADATKU ENERGIJSKA JEDNACINA ZA UNUTRASNJU ENERGIJU
% (ODN. TEMPERATURU) BIO BI PRISUTAN IZVORNI CLAN U FORMI FUNKCIJE "VISKOZNE
% DISIPACIJE" JEDNAK:
% mi*{ (du/dy+dv/dx)^2+2[(du/dx)^2+(dv/dy)^2]-2/3(du/dx+dv/dy)^2 }
% GDJE JE mi DINAMICKA VISKOZNOST FLUIDA. SVODI SE NA 4*mi. OVDJE JE ZANEMAREN.
% ON BI FIGURISAO U IZVORNOM CLANU, I TO KAO KONSTANTNI DIO Sc, BUDUCI DA
% OVAJ IZVORNI CLAN NE ZAVISI OD ZAVISNO PROMJENJIVE ("fi") VEC OD STRUJNOG
% POLJA ("u", "v")).
%########################################################################
#

%########################################################################
#
%% GRANICNI USLOVI
% Za ivice geometrije (N,W,S,E): DIRIHLEOVI => poznate vrijednosti fi na
% tim granicama
% Napomena: Vidi sliku domena (numericke mreze) i oznacavanja.
fi = zeros(1,Ni*Nj); % VEKTOR RJESENJA, ALOKACIJA I POCETNI USLOV!!!

GUvarijanta = 1; % =1 => varijanta (1) iz skripte, odn. (2)

if GUvarijanta == 1
%% KONSTANTNE VRIJEDNOSTI NA GRANICAMA GEOMETRIJE (VARIJANTA 1 U SKRIPTI)
fis=100.; % ,,south" granica (fi, i=1,9,17,...,65)
fiw=100.; % ,,west" granica
fie=0.; % ,,east" granica
fin=0.; % ,,north" granica
% ZADAVANJE:
% ,,west" granica
for i=1:1:Nj;
fi(i)=fiw; % 1,2,3,...,8. Varijanta (1)
end
% ,,east" granica
for i=Ni*Nj-(Nj-1):1:Ni*Nj
%Varijanta (1) gran.uslova:
fi(i)=fie; % 65,66,67,...,72;
end
% ,,south" & ,,north" granica
k = 0.0; % za jednokratnu upotrebu
for i=1+Nj:Nj:Nj*(Ni-1)+1
% Varijanta (1) gran.uslova:
fi(i)=fis; % 9,17,...,57
fi(i+Nj-1)=fin; % 16,24,...,64
end
elseif GUvarijanta == 2
%% ZA VARIJANTU (2) GRANICNIH USLOVA, OTKRITI DONJI SEGMENT PROGRAMA
%% PROMJENJIVI GRANICNI USLOVI (VARIJANTA 2 U SKRIPTI)
% %................... OBJASNJENJA .....................................
%% % GRANICNI USLOVI (PROFIL fi(x) ili fi(y) DUZ NEKE GRANICE):
%% Npr: W granica: parabola sa maksimumom u Ly/2 jednakim 200
%% i minimumom u y=0 i y=Ly jednakim 100.
%% S i N granica: linearni profil od 100 na x=0 do 20 na x=Lx
%% E granica: constantni profil = 20.
%% Opisane veze:
%% fiw = -4*(cc-aa)/Ly^2. * yw^2.0 + cc;
%% fis = aa + (bb-aa)/Lx*xs;
%% fin = fis;
%% fie = bb;
aa = 120; bb = 30; cc = 200;
% y koordinate cvorova na W i E granici (1.nacin):
yw(1)= 0 ;
yw(2:1:Nj-1) = [deltay/2.:deltay:Ly-deltay/2];
yw(Nj) = Ly;
% 2.nacin:
% yw = cat(2,[0.],[deltay/2.:deltay:Ly-deltay/2],[Ly]);

% y koordinate cvorova na E stranici su:

ye=yw;

% x koordinate cvorova na S i N granici:

xs=[deltax/2.0:deltax:Lx-deltax/2.0];
xn = xs;
end
%########################################################################
#

%########################################################################
#
%% DEFINISANJE PROJEKCIJA BRZINE u,v
% U CENTRIMA STRANICA KONTROLISANIH ZAPREMINA: e,w,n,s
% ZA TE KOORDINATE CE BITI SRACUNATE BRZINE u I v PO DATOJ VEZI IZ POSTAVKE

% x koordinata za centre "e" i "w" stranica kontrolisanih zapremina CV:

for xrednibroj=2:1:Ni-1
k=(xrednibroj-2)* Nj;
for i=Nj+k+2:1:2*Nj-1+k % 9,10,...,13
xe(i)=deltax+(xrednibroj-2)*delx; % 16,17,...,20
xw(i)=(xrednibroj-2)*delx; % ...........
end % 44,45,...,48
end

% y koordinata za centre "n" i "s" stranica kontrolisanih zapremina CV:

for yrednibroj=2:1:Nj-1
k=(yrednibroj-2)* Nj;
for i=Nj+k+2:Nj:Nj+2+(Ni-3)*Nj+k
yn(i)=deltay+(yrednibroj-2)*dely;
ys(i)=(yrednibroj-2)*dely;
end
end
%########################################################################
#

%########################################################################
#
%% Difuzioni fluksevi
% To su "fluksevi" nastali diskretizacijom difuzionog clana.
% Diferencijalna Shema je CDS (centralna):
% npr: De=Le*Ae/delxe = Le*(deltay*1)/delxe,
% gdje je An povrsina sjeverne stranice na kartezijanskoj mrezi,
% delxe je rastojanje od cvorova (tezista konacnih zapremina) P i E

%_____Definisanje Dn ,,north"_________
for k=0:Nj:(Ni-3)*Nj
for i=2+Nj+k:1:2*Nj-2+k
Dn(i)=L*deltax/dely;
end
end
% specificni
for i=2*Nj-1:Nj:Nj*(Ni-1)-1
Dn(i)=2*L*deltax/dely;
end
%_____Definisanje Ds ,,south"_________
for k=0:Nj:(Ni-3)*Nj
for i=3+Nj+k:1:3+Nj+k+Nj-4
Ds(i)=L*deltax/dely;
end
end
% specificni
for i=2+Nj:Nj:Nj*(Ni-1)-Nj+2
Ds(i)=2*L*deltax/dely;
end
%____Definisanje Dw ,,west"___________
for k=0:Nj:(Ni-4)*Nj
for i=2*Nj+2+k:1:3*Nj-1+k
Dw(i)=L*deltay/delx;
end
end
% specificni (8,9,10,11)
for i=Nj+2:1:2*Nj-1
Dw(i)=2*L*deltay/delx;
end
%_____Definisanje De ,,east"
for k=0:Nj:(Ni-4)*Nj
for i=2+Nj+k:1:2*Nj-1+k
 De(i)=L*deltay/delx;
end
end
% specificni (26,27,28,29)
for i=(Ni-1)*Nj-Nj+2:1:Nj*(Ni-1)-1
De(i)=2*L*deltay/delx;
end
%########################################################################
#

%########################################################################
#
%% DEFINISANJE KOEFICIJENATA U TRANSPORTNOJ JEDNACINI
% SHEMA ZA KONVEKTIVNE FLUKSEVE JE ,,POWER LAW" SHEMA (S.PATANKAR).
% DEFINISU SE PRVO BRZINE ue, uw, us, un ZA SVAKU KONTROLISANU ZAPREMINU CV
% NA OSNOVU NJH SE RACUNAJU KONVEKTIVNI FLUKSEVI Fe, Fw, Fs, Fn
% npr: Fe=(ro*u*fi)e*deltay.
% Koeficijenti su konstatni pa se sracunavaju jednim prolazom prije GS
% solvera

%% DEFINISANJE BRZINA
% PETLJA PO SVIM CELIJAMA (vidi mrezu)
% (za sve CV-e (bez ivicnih tacaka))

% ALOKACIJA (nije neophodna u MAtlabu):

ue=zeros(1,(Ni-3)*Nj + 2*Nj-1);
uw=zeros(1,(Ni-3)*Nj + 2*Nj-1);
vn=zeros(1,(Ni-3)*Nj + 2*Nj-1);
vs=zeros(1,(Ni-3)*Nj + 2*Nj-1);
for k=0.:Nj:(Ni-3)*Nj
for i=Nj+2+k:1:2*Nj-1+k
% POSTAVKA: u=C*x, v=-C*y
ue(i) = Cuv*xe(i);
uw(i) = Cuv*xw(i);
vs(i) = - Cuv*ys(i);
vn(i) = - Cuv*yn(i);
end
end
%%

%% SRACUNAVANJE KOEFICIJENATA aw, as, an, ae, MATRICE SISTEMA:

% PETLJA PO SVIM CELIJAMA (vidi mrezu Ni=9,Nj=8):
for k=0.:Nj:(Ni-3)*Nj % za sve CV-e (bez ivicnih tacaka)
% Za mrezu Nj=8, Ni=9:
for i=Nj+2+k:1:2*Nj-1+k % 10,11,...,15
% 18,19,...,23
% 58,59,...,63

% ae(i):
%%% %%% POWER LAW shema:
%%% % A(|Pe|)=max(0,(1-0.1*|Pe|))^5)
%%% if (1-0.1*abs(ro*ue(i)*deltay/De(i))).^5 > 0
%%% max1e(i)=(1-0.1*abs(ro*ue(i)*deltay/De(i))).^5;
%%% else
%%% max1e(i)=0;
%%% end
%%% % max(-Fe,0)
%%% if -ro*ue(i)*deltay > 0
%%% max2e(i)=-ro*ue(i)*deltay;
%%% else
%%% max2e(i)=0;
%%% end
%%% ae(i)=De(i)*max1e(i)+max2e(i);

%%% UPWIND
if (-ro*ue(i)*deltay) > 0
maxe(i)= - ro*ue(i)*deltay;
else
maxe(i)=0;
end
ae(i) = De(i) + maxe(i);
 % aw(i):
%%% %%% POWER LAW shema:
%%% % A(|Pw|)=max(0,(1-0.1*|Pw|))^5)
%%% if (1-0.1*abs(ro*uw(i)*deltay/Dw(i))).^5 > 0
%%% max1w(i)=(1-0.1*abs(ro*uw(i)*deltay/Dw(i))).^5;
%%% else
%%% max1w(i)=0;
%%% end
%%% % max(Fw,0)
%%% if ro*uw(i)*deltay > 0
%%% max2w(i)=ro*uw(i)*deltay;
%%% else
%%% max2w(i)=0;
%%% end
%%% aw(i)=Dw(i)*max1w(i)+max2w(i);

%%% UPWIND:
if (ro*uw(i)*deltay) > 0
maxw(i)=ro*uw(i)*deltay;
else
maxw(i)=0;
end
aw(i) = Dw(i) + maxw(i);

% an(i):
%%% %%% POWER LAW shema:
%%% % A(|Pn|)=max(0,(1-0.1*|Pn|))^5)
%%% if (1-0.1*abs(ro*vn(i)*deltax/Dn(i))).^5 > 0
%%% max1n(i)=(1-0.1*abs(ro*vn(i)*deltax/Dn(i))).^5;
%%% else
%%% max1n(i)=0;
%%% end
%%% % max(-Fn,0)
%%% if -ro*vn(i)*deltax > 0
%%% max2n(i)=-ro*vn(i)*deltax;
%%% else
%%% max2n(i)=0;
%%% end
%%% an(i)=Dn(i)*max1n(i)+max2n(i);

%%% UPWIND:
if (-ro*vn(i)*deltax) > 0
maxn(i)= -ro*vn(i)*deltax;
else
maxn(i)=0;
end
an(i)=Dn(i) + maxn(i);

% as(i):
%%% %%% POWER LAW shema:
%%% % A(|Ps|)=max(0,(1-0.1*|Ps|))^5)
%%% if (1-0.1*abs(ro*vs(i)*deltax/Ds(i))).^5 > 0
%%% max1s(i)=(1-0.1*abs(ro*vs(i)*deltax/Ds(i))).^5;
%%% else
%%% max1s(i)=0;
%%% end
%%% % max(Fs,0)
%%% if ro*vs(i)*deltax > 0
%%% max2s(i)=ro*vs(i)*deltax;
%%% else
%%% max2s(i)=0;
%%% end
%%% as(i)=Ds(i)*max1s(i)+max2s(i);

%%%% UPWIND
if (ro*vs(i)*deltax) > 0
maxs(i)=ro*vs(i)*deltax;
else
maxs(i)=0;
end
 as(i)=Ds(i)+maxs(i);
end
end
% => KOEFICIJENTI aw,as,an,ae SU DEFINISANI.
%########################################################################
#

%########################################################################
#
%% GLAVNI KOEFICIJENT ap I SLOBODNI CLAN SISTEMA JEDNACINA bp
% Iz diskretizacija: Koeficijent ap je u opstem slucaju (vidi vjezbe):
% ap = ae+as+aw+an + ap0 - Sp*deltax*deltay;
% - kod stacionarnih problema nema ap0
% - ako izvorni clan S ne zavisi od zavisno promenjive "fi", onda je
% Sp nula.

% Sracunavanje po svim CV (u proracunskom dijelu domena)

% tj. za sve CV-e u unutrasnjosti (vidi mrezu):
ap = zeros(1,(Ni-1)*Nj-1); % ALOKACIJA
b = zeros(1,(Ni-1)*Nj-1); % ALOKACIJA
for k=0:Nj:Nj*(Ni-3)
for i=2+Nj+k:1:2*Nj-1+k
ap(i) = ae(i)+aw(i)+as(i)+an(i); % ukupno ap (glavni koeficijent)
b(i) = Sc*deltax*deltay; % slobodni clan (za pocetno fi)
end
end
%########################################################################
#

%########################################################################
#
%########################################################################
#
%########################################################################
#
%% RJESAVANJE SISTEMA ALGEBARSKIH JEDNACINA - GAUS SEIDELOVIM ITERATIVNIM
% SOLVEROM SA RELAKSACIJOM
%
% Iterativno resavanje sistema linearnih algebarskih jednacina.
% Algoritam je Gaus-Seidel sa relaksacijom

Rez = epsilon + 1.; % Normirani rezidual, nulta iteracija

rezidual = zeros(1,(Ni-1)*Nj-1); % Alokacija
brojac = 0.; % pomocni brojac iteracije
while Rez > epsilon % Kriterijum za prestanak iteriranja
brojac = brojac + 1. ;

%% GS GS GS GS GS GS GS GS GS GS GS GS GS GS GS GS GS GS GS GS
% SRCE SOLVERA: RACUNANJE PETLJOM PO SVIM CV (PRORACUNSKIM
% CVOROVIMA "P", VIDI SKICU U SKRIPTI), ZA SVAKU CV (BEZ
% CVOROVA NA GRANICAMA ), PO FORMULI:
% fi^(k+1)=omega*[(b + suma_i(a*(fi)^(k))i]/ap + (1-omega)*fi^(k),
% gdje je omega koef. relaksacije (izmedju 0 i 1):

for k=0:Nj:(Ni-3)*Nj
for i=2+Nj+k:1:2*Nj-1+k
fi(i)=omega*(b(i)+as(i)*fi(i-1)+aw(i)*fi(i-Nj)+an(i)*fi(i+1)+ae(i)*fi(i+Nj))./ap(i) + (1-omega)*fi(i);
end
end

%#####################################################################
% GRANICNI USLOV NA IZLAZNOJ "E" GRANICI JE NOJMANOV: grad(fi)*nE = 0.
% tj. VRIJEDNOST "fi" NA "E" GRANICI JEDNAKA JE VRIJEDNOSTIMA UZ GRANICU.
% Za izlazne (e) cvorove => fi(i)=fi(i-Nj) jer je (dfi/dx)e=0
for i= Nj*(Ni-1)+2:1:Ni*Nj-1 % 66,67,...,71
fi(i)=fi(i-Nj);
end
%#####################################################################
% KRAJ ITERACIONE SHEME (TEKUCE ITERACIJE)
%% GS GS GS GS GS GS GS GS GS GS GS GS GS GS GS GS GS GS GS GS

%############################################################
%RRRRRRRRRRRRRRRRRRRRRRRRRRRRRR
% Svi fi(i) su sada sracunati. Sada se moze odrediti vektor
% Reziduala: u k-toj iteraciji to je vektor ciji su elementi:
% Rp^(k) = b + suma(a_i*fi_i^(k))_susjedi - ap*(fip)^(k), tj.
% to je u (k)-toj iteraciji vektor sa, u principu, ne-nultim
% elementima. Iteriranjem se moduo njegovih elemenata smanjuje
% i tezi nuli.
% Koji element vektora je mjerodavan: (i) srednja vrijednost
% vektora, (ii) maksimalni element vektora, (ii) rms, i sl.
% Ovdje je usvojen maximalni element kao mjerodavan
Rez = 0.0; % Naziv za maximalni element i za normirani rezidual
for k=0:Nj:(Ni-3)*Nj
for i=2+Nj+k:1:2*Nj-1+k
% Ri =
rezidual(i)=abs(b(i)+as(i)*fi(i-1)+aw(i)*fi(i-Nj)+an(i)*fi(i+1)+ae(i)*fi(i+Nj)-ap(i)*fi(i));
%%% % 1. NACIN DA SE ODREDI MAXIMALNI ELEMENT (PO MODULU) U VEKTORU?
%%% if abs(rezidual(i)) > Rez
%%% Rez = abs(rezidual(i)); %Najveci se zove Rez
%%% end
end
end
% 2. NACIN DA SE ODREDI max(abs(Ri))
% KORISTECI MATLAB-ove funkcije:
Rez = max(abs(rezidual));
% VRIJEDNOST max(R) POSLE PRVE ITERACIJE
% BICE SACUVANA KAO Rez_1, RADI NORMIRANJA:
if brojac == 1.
Rez_1 = Rez; % Najveci rezidual, nakon 1.iteracije
end
% NEKA JE NORMIRANI REZIDUAL OZNACEN ISTOM
% OZNAKOM KAO MAKSIMALNI, PREPISIVANJEM
% PREKO:
Rez = Rez/Rez_1;
% MJERODAVNO ZA PREKID ITERIRANJA - KADA Rez < epsilon !!!
%RRRRRRRRRRRRRRRRRRRRRRRRRRRRRR
%############################################################

% GRANICNIK U SMISLU MAX. BROJA ITERACIJA KOJI DOPUSTAMO:

if brojac > MaxIter; % Gornja granica broja iteracija GS solvera
break
end
end % KRAJ GAUS-SEIDEL SOLVERA
%########################################################################
#
%########################################################################
#
%%
#########################################################################
########

%########################################################################
#########
%% CRTANJE GRAFIKA "fi"
% Treba resenje dobijeno u formi vektora "fi", dimenzije (1,Ni*Nj)
% pretvoriti u matricu dimenzija kako je i sama mreza, u skladu sa
% sintaksom funkcija za crtanje funkcije dvije promjenjive. Vidi prve
% casove vjezbi.
%***********************************************************************
% | |
% | |
 for j=1:1:Ni % | (n) |
k=(j-1)*Nj; % |_____________|
for i=1:1:Nj % \ | \ |
FI(i,j)=fi(i+k); % \ | \ |
end % (w) \ | (domen) \ | (e)
end % Ly \ | \|
% \|____________|
% (s) Lx
%***********************************************************************
% (x,y) KOORDINATE TACAKA U KOJIMA JE DEFINISANA VRIJEDNOST "fi"
% To su: koordinate cvorova na ivicama i centara konacnih zapremina
% unutar domena.
x = cat(2,[0. deltax/2.],[3*deltax/2.0:deltax:(Lx-deltax/2.0)], [Lx]);
y = cat(2,[0. deltay/2.],[3*deltay/2.0:deltay:(Ly-deltay/2.0)], [Ly]);
[xp,yp] = meshgrid(x,y); % DOMEN (x,y) TACAKA

%Grafici rjesenja:
figure(1), pcolor(xp,yp,FI), shading faceted, colorbar('vert'),
title(' RJESENJE: \phi(x,y) \newline NA PRESJECIMA LINIJA NALAZE SE PRORACUNSKI
CVOROVI'),
xlabel('0<x<Lx'), ylabel('0<y<Ly'),
axis image

figure(2), surfc(xp,yp,FI), shading interp, material shiny, title('RJESENJE \phi=f(x,y)'),

xlabel('0<x<Lx'),
ylabel('0<y<Ly'),
zlabel('Vrijednost funkcije \phi');
axis square

figure(3), A=contourf(xp,yp,FI); clabel(A); colorbar('vert'),

title('RJESENJE \phi(x,y)'), xlabel('0<x<Lx'), ylabel('0<y<Ly'),
zlabel('Vrijednost \phi');
axis image

figure(4), pcolor(xp,yp,FI), axis square,shading interp,colorbar('vert'),

title('FUNKCIJA \phi U DOMENU'),
xlabel('0<x<Lx'),ylabel('0<y<Ly'),
axis image

figure(5), A = contour(xp,yp,FI);clabel(A), axis square,

title('FUNKCIJA \phi U DOMENU'),
xlabel('0<x<Lx'),ylabel('0<y<Ly'),
axis image
%%
########################################################################
break
% KRAJ PROGRAMA

%% IMA LI RAZLIKE IZMEDJU UPWIND I POWER LAW REZULTATA NA ISTOJ MREZI ?

fiupwind = fi;
% prethodno sacuvano:
load fipowerlaw.mat
fipowerlaw = fi;
RAZLIKA = (fiupwind - fipowerlaw)./fiupwind*100; %% [%]

% | |
for j=1:1:Ni % | (n) |
k=(j-1)*Nj; % |_____________|
for i=1:1:Nj % \ | \ |
FI(i,j)=RAZLIKA(i+k); % \ | \ |
end % (w) \ | (domen) \ | (e)
end % Ly \ | \|
% \|____________|
% (s) Lx

figure(111), pcolor(xp,yp,FI),
shading interp,colorbar('vert'),
title('RAZLIKA U % U DOMENU'),
xlabel('0<x<Lx'),
ylabel('0<y<Ly'),
axis image
%%%% R: Na mrezi 30X30: max. razlike je 0.016118596333645 [%]