kill(all)$

maxima: finite element (=discrete), version 1.0,part 2
peter.vlasschaert@gmail.com,04/01/2021

load packages,option
=================

load(lrats)$ ratprint:false $

physical system :ki (= physical constant)

see:
maxima: finite element (=discrete), version 1.0,part 1
peter.vlasschaert@gmail.com,03/01/2021

def1:

numel(ll):=[ll,ll+1];

def2:stiffness matrix

K:genmatrix(v, 3+1, 3+1);

def3:elemental matrix

elmat(h):=matrix([k[h],-k[h]],[-k[h],k[h]]);

def4: zero matrix

ss:makelist(makelist(v[i2,j2]:0,i2,1,3+1),j2,1,3+1);

def5: elemental matrix into def2

ellmatt(ei,mm):=block([ik,jk,ii,jj],
s1:elmat(ei),
s2:makelist(numel(jj),jj,1,mm),    
s3:makelist(jj,jj,0,mm-1),/*correction 2by2 matrix*/         
for ik:s2[ei][1] thru s2[ei][2] do(
    for jk:s2[ei][1] thru s2[ei][2] do(            
        v[ik,jk]:elmat(ei)[ik-s3[ei],jk-s3[ei]]        
    )  
  )   
)$  

def6:combine all elemental matrix

sum(ellmatt(ii,4-1),ii,1,4-1)$

ev(K,nouns);
x:k[1]$/* physical matrix:ll =1*/
y:k[2]$/* physical matrix:ll =2*/
v[2,2]:v[2,2]+x$
v[3,3]:v[3,3]+y$
K:ev(K,nouns);
KK:K;

solution : figure 1
==============

matsize(AA):=[length(AA),length(transpose(AA))];

z1:matsize(K);

step 1 : system of equations "physical model"
-------------------------------------------------------------

lhs: system

z2:makelist(u[i],i,1,z1[1]);
z2a:transpose(z2);

z2b:K.z2a;

→ rem

z2b[1][1];

rhs: system

z2c:transpose(makelist(F[i],i,1,z1[1]));

→ rem

z2c[1][1];

lhs=rhs

z2d:makelist(z2b[i][1]=z2c[i][1],i,1,z1[1]);

z2d1:reverse(z2d);

step 2: constraint
------------------------

u[1]:0;
F[1];    /* reaction force    */
F[2]:0;  /* no internal force */
F[3]:0;  /* no internal force */
F[4]:P;  /* external force    */

step 3:apply constraint on system
---------------------------------------------

z3a:ev(z2d,nouns);

z3a1:sort(listofvars(z3a));

z3a2:solve(z3a,[u[2],u[3],u[4]]);

step 4:BC ( = Boundary Condition , u[1]=0)
---------------------------------------------------------
see: website
new: 20/12/2020:maxima :apply general BC on a stiff system ( 'finite element'),version 1.0
download  (Maxima,code)       :maximageneralbcappliedtomodelfiniteelmementversion1
download  (Maxima,png)         :maximageneralbcappliedtomodelfiniteelmementversion1outp

→ u[1] : delete first row & delete first column (use :submatrix 'maxima')

z3b:submatrix(1,KK,1);

z3b1:transpose([0,0,P]);
z3b2:transpose([u[2],u[3],u[4]]);

z3b3:z3b.z3b2;

z3b3a:z3b.z3b2 =z3b1;

lhs=rhs

z3b4:makelist(z3b3[i][1]=z3b1[i][1],i,1,z1[1]-1);

step 5: solve:u[i],i=2,3,4
---------------------------------
because : u[1]=0 , fixed 

z3c:flatten(solve(z3b4,[u[2],u[3],u[4]]));

step 6: postprocessing
---------------------------------
use :z3c into z3a

z3c1:lratsubst(z3c,z3a);

rem :reaction force :F[1] = -P @ the wall ( = inside the wall)