kill(all)$

maxima : fvm : non regular grid ( non-orthogonal,grid =triangles),diffusion flux output ,version 1.0
peter.vlasschaert@gmail.com,04/02/2023

load packages :
--------------------

load(lrats)$

find w vector

construction : with geogebra 6 classic

graphics
------------

s1:Centroid triangle = G ,H
s2: intersect(s1,|AC|) =  q
w1,w2 = unit vectors (along |AC| ,                    along |GH| )
w         = vector ,perpendicular  on the line segment (=|AC|)
-------------------------------------------------------------------------------------

computation of the flux:
==================

---------------------------------------------------------------------------------------------------------

G = p , reference triangle
H = c1,neighbourhood triangle
w1 = unit vector ,direction:  A to C: positive :  e[η], counter clock wise.
w2 = unit vector ,direction : G to H : positive : e[z]
w   = vector         , perpendicular to  e[η]
-----------------------------------------------------------------------------------------
AA[f] = w , on side of atriangle ( ≠ unit vector) ,e[η] perpendicular |AC|
-----------------------------------------------------------------------------------------
cartesian grid              : unit vector ,e[x],e[y]
non-orthogonal grid  : unit vector , e[η],e[z]
------------------------------------------------------------------------------------------
length : |AC|  = Δη
length : |GH|  = Δz
------------------------------------------------------------------------------------------

-------------------------------------------------------------------------------------------------------

Volume infinitesimale
------------------------------

→ first step : centroid

ref : triangle : p

→ centroid : p ( figure : G)

p1a:[x[p],y[p]];

neighborhood: triangle : c1 ( call adjacent triangle on to ref)

→ centroid : c1 (figure : H)

p1b:[x[c1],y[c1]];

→  second step : define unit vectors :e[η],e[z]

p2a:e[η] =((x[A]-x[C])*e[x]+(y[A]-y[C])*e[y])/Δη;

p2b:e[z] =((x[c1]-x[p])*e[x]+(y[c1]-y[p])*e[y])/Δz;

→ third step : find A[f], perpendicular e[η]

rem : A[f]= function(e[η]),vector

p3:AA[f] =AA[fx]*e[x]+AA[fy]*e[y];

1e) unit vector : relations,used scalair product or dot product

p3a:[e[x]*e[x]=1,e[y]*e[y]=1,e[x]*e[y]=0];

2e) vector :A[fx],A[fy]

p3b:AA[f]=(y[A]-y[C])*e[x]-(x[A]-x[C])*e[y];

3e) test A[f] perpendicular e[η] = scalair product = 0

p3c:ratsimp(p3b*p2a);

p3c1:(lratsubst(p3a,p3c));

Diffusion flux : out the face , |AC|
--------------------------------------------

depends(Φ,[x,y]);

→ first step :gradient : cartesian coordinates (2D)

p4: ∇(Φ) = diff(Φ,x)*e[x]+diff(Φ,y)*e[y];

p4a:expand(D*p4);

→ second step : integration over CV 

simp :false$

 > conversion : CV  to CS

p4aa:'integrate(∇.'(D*∇(Φ)),V) ='integrate('(Γ*'∇(Φ)).'unit(n),A);

 > discrete : over all faces ,sum operator

p4bb:lhs(p4aa)='∑*'(( B1)[f].(B2)[f]);

simp:true$

f ={w,e,s,n} ,faces(=CS) of CV
--------------------------------------
B1:D*∇(Φ)       ⇒ vector , D = diffusion constant
B2:(vec(dA))   ⇒ vector
B1[f] . B2[f]    ⇒ dot product,see : definition with infix operator
--------------------------------------
vec(dA) = unit(n)*dA
-------------------------------------

→ third step : for face (=f) between G( = p) and H(= c1)  triangle
use :p3b

p4c:p4a*p3;

dot product : use p3a

p4c1:expand(p4c);

p4c2:lratsubst(p3a,p4c1);

→ fourth step : diif(Φ,x),diif(Φ,y) & diif(Φ,η),diif(Φ,z)

depends(Φ,[η,z]);

p5a:diff(Φ);

p5a1:part(p5a,1,1) = 'diff(Φ,x)*'diff(x,η)+ 'diff(Φ,y)*'diff(y,η);
p5a2:part(p5a,2,1) = 'diff(Φ,x)*'diff(x,z)+ 'diff(Φ,y)*'diff(y,z);

p5a3:[p5a1,p5a2];

p5a33:[lhs(p5a1),lhs(p5a2)];

use : p4

p5a4:[ratcoeff(rhs(p4),e[x]),ratcoeff(rhs(p4),e[y])];

coefmatrix(p5a3,p5a4)*(-1);

p6:flatten(solve(p5a3,p5a4));

fifth step : combine : fourth step & third step.

use :p6 & p4c2

p7:ratsimp(subst(p6,p4c2));

p7a1:num(rhs(p7));
p7a2:denom(rhs(p7));

sixth step : derivatives discrete 
use : p7a2

p8a1:args(p7a2)*[1,-1];

p8a2:sort(flatten(makelist([part(p8a1[i],1),part(p8a1[i],2)],i,1,length(p8a1))));

rem : discrete :'diff(x,z,1) = Δx/Δz,'diff(x,η,1) =Δx/Δη
                       ,'diff(y,z,1)=Δy/Δz,'diff(y,η,1) =Δy/Δη

use: p2a,p2b

p8a1a:p2a;
p8a2b:p2b;

→ x ,with respect Δη,Δz

p8a4:p8a2[1] = ratcoeff(rhs(p8a2b),e[x]);
p8b4:p8a2[2] = ratcoeff(rhs(p8a1a),e[x]);

→ y ,with respect Δη,Δz

p8c4:p8a2[3] = ratcoeff(rhs(p8a2b),e[y]);
p8d4:p8a2[4] = ratcoeff(rhs(p8a1a),e[y]);

final result : discrete derivatives :

p8e5:[p8a4,p8b4,p8c4,p8d4];

seventh step :combine : sixth step & fifth step
use :p7

p9:subst(p8e5,p7);

p9a:denom(rhs(p9));

p9aa:ratsimp(expand(p9a));

use : p3b.p2b

p9a1:expand(subst([p2b,p3b],AA[f]*e[z]/Δη));

use : p3a

p9a2:lratsubst(p3a,p9a1);

result : denom : rhs
--------------------------

p9a3:AA[f]*e[z]/Δη;

check :

rem : expand → ratsimp

is(p9aa=p9a2);

find  here : num construct ?
------------------------------------

p9b:num(rhs(p9));

use p5a33,['diff(Φ,η,1),'diff(Φ,z,1)]

p9c1:ratcoeff(p9b,p5a33[1]);
p9c2:ratcoeff(p9b,p5a33[2]);

first : p9c2,'diff(Φ,z,1)
--------------

use :p3b,p3

x1a:reverse(p3b);

use :rules for dot product

--------------------------------------------------------------------------

x1b1:lratsubst(p3a,p3*e[x]);
x1b2:lratsubst(p3a,p3*e[y]);

x1b3:p3;

x1c:reverse(lhs(x1a) =rhs(x1b3));

x1ca:reverse(ratcoeff(x1c,e[x]));
x1ca1:reverse(x1ca);
x1cb:reverse(ratcoeff(x1c,e[y]));
x1cb1:reverse(x1cb);

--------------------------------------------------------------------------------

use : p9c2 ,final

p9cc2:factor(lratsubst([x1ca,x1cb],p9c2));

D*AA[f].AA[f]/Δη

first : p9c1
--------------
use :x1ca,x1cb

p10a:lratsubst([x1ca1,x1cb1],p9c1);

computation:e[z]*e[η]
use :p1

p10b:D*p2a*p2b;

p10b1:-Δη*p10b;

p10b2:lratsubst(p3a,p10b1);

check :p10a & rhs(p10b2)

p10b3:is(p10a=rhs(p10b2));

final result :use p9cc2 & p10b2 , num : p9a3,use p5a33 "derivatives"
--------------------------------------------------------------------------------------

simp:true$

p11a:lhs(p9) =ratsimp(p9cc2/p9a3*p5a33[1]+lhs(p10b2)/p9a3*p5a33[2]);

rem :x1ca1,x1cb1

p12:p3;

p12a:lratsubst([x1ca1,x1cb1],p12);

p12b:p12a*p12a;

p12c:(lratsubst(p3a,expand(p12b)));

p12c1:lhs(p12c) =factorsum(rhs(p12c));/*distance (Δη)^2*/

p12c2:reverse(lhs(p12c1) =(Δη)^2);

use :p11a

p12d:subst([p12c2],p11a);

rem : cartesian coordinates,"see above"

p12dx:(AA[f])^2 =(AA[fx])^2+(AA[fy])^2;
p12dx1:reverse(p12dx);

final : final : result:
D*AA[f]*∇(Φ)=(-1)*1/AA[f]*1/e[z]*(-D*(AA[f])^2)*('diff(Φ,η,1))+D*e[z]*AA[f]^2*('diff(Φ,z,1))*e[η])

----------------------------------------------------------------------------------------------------------------------------------------

reason : AA[f],vector .