kill(all)$

maxima:fvm:non-orthogonal mesh ,between two cells { p[0],p[1]},version 1.0
peter.vlasschaert@gmail.com,26/12/2022

ref : maxima : nonuniform mesh diffusion 2D with sourc S, gauss divg theorema(J[q] flux,LHS  wall),version 1.0

load packages
------------------

load(lrats)$

math : model " diffusion + S =0"
-----------------------------------------

gauss divergence theorema
----------------------------------

simp :false$

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

ma2:lhs(ma1)='∑*'(( B1)[f].(B2)[f]);/*sum,with dot product inside*/

B1= Γ*∇(Φ) =Γ*∇Φ →( Γ*∇(Φ))[f]  ,B2 =A → A[f],f =face

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

simple :Γ*∇Φ =Γ*∇(Φ)=Γ*Φ[x]*e[x]+Γ*Φ[y]*e[y]

ma3s:Γ*∇(Φ) =Γ*Φ[x]*e[x]+Γ*Φ[y]*e[y];

dot-rule
-----------

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

coordinate transformation:
-----------------------------------
x → x(ε,η)
y → y(ε,η)
-----------------------------------
p[0] =centroid ,cell volume
p[1] =centroip1aa:Φ[ε] =Φ[x]*x[ε] +Φ[y]*y[ε];d ,cell volume

p1a:'diff(Φ,ε) = 'diff(Φ,x)*'diff(x,ε)+ 'diff(Φ,y)*'diff(y,ε);
p1b:'diff(Φ,η) = 'diff(Φ,x)*'diff(x,η)+ 'diff(Φ,y)*'diff(y,η);

simple notation :{p1a,p1b}

p1aa:Φ[ε] =Φ[x]*x[ε] +Φ[y]*y[ε];
p1bb:Φ[η] =Φ[x]*x[η] +Φ[y]*y[η];

solve : {Φ[x],Φ[y]}

p1cc:flatten(solve([p1aa,p1bb],[Φ[x],Φ[y]]));

computation : unit_vectors ={e[ε],e[η]}
---------------------------------------------------

→ e[ε]: along p[0] = p to p[1]=a
→ e[η]: along A to B

simp:true;

p2a:e[ε] =((x[a]-x[p])*e[x]+(y[a]-y[p])*e[y])/Δε;
p2b:e[η] =((x[B]-x[A])*e[x]+(y[B]-y[A])*e[y])/Δη;

p2c:lratsubst(pd,p2a*p2b);

vector : A[f] perpendicular on A to B
------------------------------------------------

p3a:A[f] =A[x]*e[x]+A[y]*e[y];

suppose :vector C =A[c] =(y[B]-y[A])*e[y]+(x[B]-x[A])*e[x]

p3b:A[c] =(y[B]-y[A])*e[y]+(x[B]-x[A])*e[x];

orthogonal : A[f],A[c]

p3c1:A[f]*A[c] =0;/*condition*/
p3c2:expand(rhs(p3a)*rhs(p3b)) =0;

use : dot rule

p3c3:lratsubst(pd,p3c2);

set :A[x] =y[B]-y[A]

p3c4:A[x] =y[B]-y[A];

→ subst : p3c4 into p3c3

p3c4a:subst([p3c4],p3c3);

→ find : A[y]

p3c4b:solve(p3c4a,A[y])[1];

→ plug in : p3a,{p3c4b,p3c4}

p3d:[p3c4,p3c4b];
p3da:subst(p3d,p3a);

→ use :{ma3s,p3a}

p4a:rhs(ma3s);

p4b:rhs(p3a);

face =f,see above : math

p4c:p4a*p4b;

dot product :  "f"

p4c1:lratsubst(pd,expand(p4c));

use : simple p1cc

p4c2:ratsimp(lratsubst(p1cc,p4c1));

→ num

p4c2a:(args(p4c2));

p4c2b:p4c2a[1];

p4c2b1:collectterms(num(p4c2b),Φ[η],Φ[ε]);
p4c2b2:denom(p4c2b);

first part :x[ε]*y[η]-y[ε]*x[η]
---------------------------------------
see :p2a,p2b

differential :

p5a1:x[ε]=(x[a]-x[p])/Δε;
p5a2:y[ε]=(y[a]-y[p])/Δε;
p5b1:x[η]=(x[B]-x[A])/Δη;
p5b2:y[η]=(y[B]-y[A])/Δη;
p5cx:[p5a1,p5a2,p5b1,p5b2];

into p4c2b2

p5xx:subst(p5cx,p4c2b2);

p5x1:ratsimp(p5xx);

use computation :e[η],A[f],p3da

p6a:A[f]*e[ε] =lratsubst(pd,rhs(p3da))*rhs(p2a);

p6a1:expand(lratsubst(pd,rhs(p2a)*rhs(p3da)));

find : mq

p7b : p5x1=mq*p6a1;

p7b1:solve(p7b,mq);

final :num(p4c2b2)

p8a:(p4c2b2)=subst(p7b1,mq*lhs(p6a));

use : p4c2b1

p9a:p4c2b1;

declare(Γ,constant)$

coef:Φ[ε]

p9a1:ratsimp(ratcoeff(p9a,Φ[ε]));

p9a1a:subst(p5cx,p9a1);

p9a1ax:Γ*(A[f]*A[f])/Δη;

p9a1ay:p9a1=p9a1ax;

coef:Φ[η]

p9a2:ratsimp(ratcoeff(p9a,Φ[η]));

p9a2a:subst(p5cx,p9a2);

use :{p3c4,p3c4b}

p9a2ax:ratexpand(lratsubst([p3c4,p3c4b],p9a2a));

p9a2axx:ratexpand(1/Γ*(p9a2ax));

see :p2c

p9a2ax1:ratexpand(ratexpand(p2c)*Δη);

p10a:is(p9a2axx=rhs(p9a2ax1));

p10b:Γ*lhs(p9a2ax1) *(-1)=p9a2;

p10bb:reverse(p10b);

final result : f=face

→ p4c2b1,p4c2b2

p11a:p4c2b1/p4c2b2;

use :p10bb,p9a1ay,p8a

p11b:[p10bb,p9a1ay,p8a];

p11c:lratsubst(p11b,p11a);

result :(Γ*∇[f](Φ).A[f]),Γ =constant
---------------------------------------------

→ first final result :

p12:Γ*∇[f](Φ).A[f] =(p11c);

rem : see figure 1 

what about :Δη an A[f]

→ step 1:distance : Δη ,x-y coordinates

p13a:(Δη)^2=(y[B]-y[A])^2+(x[B]-x[A])^2;

p13a1:expand(p13a);

→ step 2: definition A[f]

use :p3da

p13b:p3da;

p13c:expand(p13b*p13b);

→ step 3: dot rule

use:pd

p13d:lhs(p13c)=(lratsubst(pd,rhs(p13c)));

p13d1:lhs(p13d)=factorsum(rhs(p13d));

→ step 4:final solution

p14 : [lhs(p13a)=lhs(p13d1),is(rhs(p13a)=rhs(p13d1))];

p14a:[p14[1]];

rem : A[f]^2 =A[f]*Δ[f]=A[f].A[f], last is dot product.

→ second final result : "p14a"

use : p12

p15:subst(p14a,p12);

rem : A[f]^2 =A[f]*Δ[f]=A[f].A[f], last is dot product.

rem : distrib,multthru,expandwrt,xthru,combine,facsum,ratexpand,
         collectterms,collect,factorsum......

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