%% Heat equation with circular transport effects on a circular plate
% The general transport equation has the following generic form:
% \frac{\partial u}{\partial t} = \vec v \cdot \mathbf{\nabla} \, u \, - \,  \mathbf{\nabla} \cdot \Gamma \,\mathbf{\nabla} \phi + S
% Mi a következő egyenletet szeretnénk megoldani:
% \frac{\partial u}{\partial t} = \vec v \cdot \mathbf{\nabla} \, u + k \Delta u , ahol \vec v = \frac{1}{r + \epsilon}\pmatrix{ \sin(\phi) \cr -\cos(\phi) }, (r,\phi) a polárkoordinátás megadása (x,y) pontnak. Tehát \vec vegy örvénylő mozgást kölcsönöz az egyenletnek és az örvénylés erősödik a központ felé.
%
% A Matlab PDE Toolbox végeselemű megoldói a következő formájú PDE egyenletekkel tudnak mit kezdeni
% m(t,x,u,\nabla u) \frac{\partial^2 u}{\partial t^2} + d(t,x,u,\nabla u) \frac{\partial u}{\partial t} - \nabla(c(t,x,u,\nabla u) \nabla u) + a(t,x,u,\nabla u)\, u = f(t,x,u,\nabla u)
% A peremfeltételek pedig:
% Dirichlet: h(t,x,u) \, u = r(t,x,u)
% Neumann: \vec n \cdot c(t,x,u,\nabla u) \cdot \nabla u + q(t,x,u) \, u = g(t,x,u), ahol \vec n a perem x pontjáből kifelé mutató a peremre merőleges normálvektor
% Az m, d, c, a, f együtthatókra vannak kisebb-nagyobb megkötések, lásd: Specify m, d, a coefficients, Specify f coefficient, Specify c coefficient. Egy együttható-függvények egy alkalmas megválasztásával a kívánt PDE alakot érhetjük el:
% \underbrace{m(t,x,u,\nabla u) \frac{\partial^2 u}{\partial t^2}}_{0} + \underbrace{d(t,x,u,\nabla u) \frac{\partial u}{\partial t}}_{\frac{\partial u}{\partial t}} - \underbrace{\nabla(c(t,x,u,\nabla u) \nabla u)}_{k \Delta u} + \underbrace{a(t,x,u,\nabla u)\, u}_{0} = \underbrace{f(t,x,u,\nabla u)}_{\vec v \cdot \mathbf{\nabla} \, u}
% Tehát:
% \left\{\begin{array}{l}
% m(t,x,u,\nabla u) = 0
% \\
% d(t,x,u,\nabla u) = 1
% \\
% c(t,x,u,\nabla u) = k
% \\
% a(t,x,u,\nabla u) = 0
% \\
% f(t,x,u,a,\nabla u) = \vec \nabla \cdot \vec v \,u
% \end{array}\right.
%
% Mindenekelőtt létre kell hoznunk (a pdeModeller GUI segítségével) egy tartományt, amelyen megoldjuk a PDE-t. A pdecirc segítségével előre legyártunk két kört az x,y síkon. A pdeModeller ezt a két kört már látni fogja. A pdeModeller-ben érdemes a koordináta tengelyek határait megfelelően kihúzni: Options --> Axes Limits ... --> Check "Auto". Ezt követően mentsük el a geometriai modelljét a PDE-nek a "gd", "sf", "ns" változókba: Draw --> Export Geometry Description, ...
% Végül mentsük a változókat: save(geomFileName,"gd","sf","ns")

geomFileName = "pdeGeom_C1_C2_r5.mat";

%{
    pdecirc(0,0,5)
    pdecirc(2,2,1)
    pdeModeler

    save(geomFileName,"gd","sf","ns")
%}

if ~exist(geomFileName,'file')
    return
else
    load(geomFileName,"gd","sf","ns")
end

%
% A már korábban elmentett geometriára vonatkozó információk (gd, sf, ns) segítsgével létrehozzuk a geometrikai modellt a decsg függvény meghívásával.

pdeGeom = decsg(gd,sf,ns);
pdegplot(pdeGeom,"EdgeLabels","on","FaceLabels","on")
axis equal

%
%% PDE modell
% Létrehozunk egy PDE modellt egyetlen ismeretnel függvénnyel (lehetne akár több is: u_1(t,x), \dots, u_n(t,x). Ilyen esetben a PDE rendszerről beszélünk).

numberOfPDE = 1;
pdeModel = createpde(numberOfPDE);

% A geometriai modellt hozzáadjuk a PDE modellhez:

geometryFromEdges(pdeModel,pdeGeom);

%
%% Együttható függvények deklarálása
% Az f(t,x,u,\nabla u) függvény a Matlab fájl végén van deklarálva.

k = 0.005;

specifyCoefficients(pdeModel, ...
    "m",0,...
    "d",1,...
    "c",k,...
    "a",0,...
    "f",@fcoef);

%
%% Peremfeltételek megadása

applyBoundaryCondition(pdeModel,"neumann", ...
    "Edge",1:4, ...
    "g",0, ...
    "q",0);

%
%% A végeselem módszerhez felosztjuk a teret

msh = generateMesh(pdeModel,'Hmax',0.5);
figure;
pdemesh(pdeModel);
axis equal

%
%% Kezdeti feltétel megadása

setInitialConditions(pdeModel,0); % @(r) sin(r.x + r.y);
setInitialConditions(pdeModel,5,"Face",2);

%
%% PDE megoldása

tlist = logspace(-2,2,100);
results = solvepde(pdeModel,tlist);

sol = results.NodalSolution;

%
%% Megjelenítés

fig = figure(123);
% fig.Visible = 'on';
delete(fig.Children);
ax = axes(fig);

Tri = msh.Elements(1:3,:)';
for i = 1:numel(tlist)
    Sf = trisurf(Tri,msh.Nodes(1,:)',msh.Nodes(2,:)',sol(:,i));
    zlim([-1.1,5.1])
    clim([0,1])
    drawnow
end

%
% Az f(t,x,u,\nabla u) függvény segítségével implementáljuk a transzport jelenséget.
% f(t,x,u,\nabla u) = \vec v \nabla u = \frac{1}{r + \epsilon} \pmatrix{ \sin(\phi) \cr -\cos(\phi) }^\top \, \pmatrix{ \frac{\partial u}{\partial x} \cr \frac{\partial u}{\partial y} } = \frac{1}{r + \epsilon} \Big( \sin(\phi) \frac{\partial u}{\partial x} - \cos(\phi) \frac{\partial u}{\partial y} \Big)

function f = fcoef(location,state)
    x = location.x;
    y = location.y;
    r = sqrt(x.^2 + y.^2);
    phi = atan2(y,x);
    f = ( sin(phi).*state.ux - cos(phi).*state.uy ) ./ (r + 0.1);
end

function f = fcoef_befele(location,state)
    x = location.x;
    y = location.y;
    r = sqrt(x.^2 + y.^2);
    phi = atan2(y,x);
    
    v_rotX = sin(phi)./(r+0.1);
    v_rotY = -cos(phi)./(r+0.1);
    v_beX = x;
    v_beY = y;

    lambda = 0.95;
    vX = lambda*v_rotX + (1-lambda)*v_beX;
    vY = lambda*v_rotY + (1-lambda)*v_beY;

    f = vX.*state.ux + vY.*state.uy;
end