kill(all)$

fluid  : divergence of tensor( = tensor),version 1.0,part 4
peter.vlasschaert@gmail.com,26/05/2019

rem : before 'part 3: tensor (= first        order tensor ) = vector '
                        'part 4 :tensor (= second order tensor )                 '

load packages

load(lrats)$

load(linearalgebra)$

kronecker delta : operator
---------------------------------------

1e) : kronecker delta operator: 'maxima'
=================================

p1:kron_delta(e1,e1);

p2a:kron_delta(e2,e1);
p2b:kron_delta(e1,e2);
p2c:kron_delta(x,y);

2e) : my kronecker delta 
=====================

krondelta(s1,s2):=block([],if s1=s2 then 1 else 0);

rem  :    s1 = s2   =>  1
               s1 ≠ s2   =>  0

p2a1:krondelta(e2,e1);
p2b1:krondelta(e1,e2);
p2c1:krondelta(x,y);
p2d1:krondelta(x,x);

divergence of a vector:operator
------------------------------------------------

depends([f],[x,y,z]);

gradf(m) := diff(m,x)*e1+diff(m,y)*e2+diff(m,z)*e3;

p3a:gradf(f);

p3a1:args(p3a);

p3a2:part(p3a1[1],1);
p3a3:part(p3a1[1],2);

tensor representation : 

simp:false;

p3b1: a11.e1.e1+a12.e1.e2+a13.e1.e3;
p3b2: a21.e2.e1+a22.e2.e2+a23.e2.e3;
p3b3: a31.e3.e1+a32.e3.e2+a33.e3.e3;
p3b:flatten(p3b1+p3b2+p3b3);

result :
======
.  :  non commutative operator , see , 'p3b'
see: part3.  (why)

non commutative operator :
=======================

p3d1:args(p3b);

same : representation of  'p3d1' of a tensor.
------------------------------------------------------------------

p3d1a:[e1,e1,e1,e2,e2,e2,e3,e3,e3];

p3d1b:[e1,e2,e3,e1,e2,e3,e1,e2,e3];

p31c:[a11,a12,a13,a21,a22,a23,a31,a32,a33];

p31c1:[a[1,1],a[1,2],a[1,3],a[2,1],a[2,2],a[2,3],a[3,1],a[3,2],a[3,3]];

solve : problem non commutative operator
------------------------------------------------------------------

1e) parts : nabla operator.
======================

p4a1:part(p3a1[1],1);
p4a2:part(p3a1[2],1);
p4a3:part(p3a1[3],1);

simp:true;

2e) use krondelta
===============

by : e3
----------

p5a1:makelist(krondelta(p4a1,p3d1a[i]),i,1,length(p3d1a));

p5a11:makelist(p5a1[i]*p31c1[i],i,1,length(p3d1a));

p5a11a:makelist(p5a1[i]*p3d1b[i],i,1,length(p3d1a));

p5a11b:makelist(p5a11[i]*p5a11a[i],i,1,length(p3d1a));

by : e2
----------

p5a2:makelist(krondelta(p4a2,p3d1a[i]),i,1,length(p3d1a));

p5a22:makelist(p5a2[i]*p31c1[i],i,1,length(p3d1a));

p5a22a:makelist(p5a2[i]*p3d1b[i],i,1,length(p3d1a));

p5a22b:makelist(p5a22[i]*p5a22a[i],i,1,length(p3d1a));

by : e1    
----------

p5a3:makelist(krondelta(p4a3,p3d1a[i]),i,1,length(p3d1a));

p5a33:makelist(p5a3[i]*p31c1[i],i,1,length(p3d1a));

p5a33a:makelist(p5a3[i]*p3d1b[i],i,1,length(p3d1a));

p5a33b:makelist(p5a33[i]*p5a33a[i],i,1,length(p3d1a));

combine
========

p6:map("+",p5a11b,p5a22b,p5a33b);

p6a:[x[1],x[2],x[3]];

e1 :  component
==============

p6a1:ratcoeff(p6,e1);

p6b1:sum(p6a1[i],i,1,length(p3d1a));

p6c1:args(p6b1);

for i:1 thru length(p6c1)
do(
     if i=1 then a[i,1]:'diff(a[i,1],p6a[1]),
     if i=2 then a[i,1]:'diff(a[i,1],p6a[2]) ,
     if i=3 then a[i,1]:'diff(a[i,1],p6a[3]) 
    );    

p6d1:ev(p6b1,simp)*e1;

e2 :  component
==============

p6a2:ratcoeff(p6,e2);

p6b2:sum(p6a2[i],i,1,length(p3d1a));

p6c2:args(p6b2);

for i:1 thru length(p6c1)
do(
     if i=1 then a[i,2]:'diff(a[i,2],p6a[1]),
     if i=2 then a[i,2]:'diff(a[i,2],p6a[2]) ,
     if i=3 then a[i,2]:'diff(a[i,2],p6a[3]) 
    );    

p6d2:ev(p6b2,simp)*e2;

e3 :  component
==============

p6a3:ratcoeff(p6,e3);

p6b3:sum(p6a3[i],i,1,length(p3d1a));

p6c3:args(p6b3);

for i:1 thru length(p6c1)
do(
     if i=1 then a[i,3]:'diff(a[i,3],p6a[1]),
     if i=2 then a[i,3]:'diff(a[i,3],p6a[2]) ,
     if i=3 then a[i,3]:'diff(a[i,3],p6a[3]) 
    );    

p6d3:ev(p6b3,simp)*e3;

final result : div(tensor:second order tensor)
-------------------------------------------------------------------
=>    lower order of operation.

p7:p6d1+p6d2+p6d3;

full index notation  ( = einstein convention)
=====================================
------------------------------------------------------------------
i : summation index
j : free index
------------------------------------------------------------------

p8:'diff(a[i,j],x[i])*e[j];

tensor (=second order)  => tensor ( = first order)  = vector .


result: (website)
==============

---------------------------------------------------------------------------------------------------------------
part 3:div(vector,  'first        order')  = 3 components
part 4:div(tensor, 'second  order')  = 9 components  = (3*3)     components
----------------------------------------------------------------------------------------------------------------
general :
            div(tensor, 'thrid        order') = 27 components = (3*3*3) components
----------------------------------------------------------------------------------------------------------------
example :second order :stress tensor used , navier - stokes equation.
-----------------------------------------------------------------------------------------------------------------