# Assembling Element Matrices in Finite Element Method Problem

5 views (last 30 days)
Anilcan Taner on 28 Oct 2019
Commented: Fabio Freschi on 29 Oct 2019
I should start with mentioning I am die hard newbie who has term project about a hydraulic problem which I have to solve with Matlab.
So far I've calculted KK{i} 3x3 matrices but now I have to put them in a element matrices to proceed forward with my solution.
The matrix should look like this with the help of KK{i} matrices;
( The KK{1} and KK{2} in the solution manual )
( As you can see they have the placements in overall matrix that we're planning to achieve. The 6 1 2 that I've crossed comes from the ElemX, ElemY and ElemZ)
( For those helpful people here's the interpolation scheme)
The Overall Matrix that we're planning to achieve looks like;
I've found some helpful formulations but I am unable to fit them to my code.
So sorry if it is a lame or too chaotic to understand..
Don't mind the Turkish explainations in the code.
echo off all
clc
nodeX = [0,5,9.17,12,0,5,9.17,12,0,5,8,9.17,12]; % noktalarin x koordinatlari
nodeY = [8,8,8,8,4,4,5.5,5.5,0,0,0,2.83,4]; % noktalarin y koordinatlari
elemX = [5,6,6,7,7,7,5,9,10,11,6,7,7]; % i icin olusan ucgenlerin global nokta karsiliklari
elemY = [1,1,2,2,3,4,6,6,6,6,7,12,8]; % j icin
elemZ = [6,2,7,3,4,8,9,10,11,12,12,13,13]; % k icin
varTable = [];
fprintf("e xi\txj\txk\tyi\tyj\tyk\tck\tci\tcj\tbi\tbj\tbk\tA\n"); % tablo basligini formatla
for i = 1:length(elemX)
% sirasiyla elemanlari okutuyoruz
currElemX = elemX(i);
currElemY = elemY(i);
currElemZ = elemZ(i);
% i ye karsilik gelen global noktalarin koordinat karsiliklari
% bulunuyor
xi = nodeX(elemX(i));
xj = nodeX(elemY(i));
xk = nodeX(elemZ(i));
yi = nodeY(elemX(i));
yj = nodeY(elemY(i));
yk = nodeY(elemZ(i));
% gerekli hesaplama yapiliyor
ck = xj - xi;
ci = xk - xj;
cj = xi - xk;
bi = yj - yk;
bj = yk - yi;
bk = yi - yj;
% alan hesabi
A = abs(0.5 * (xi * yj + xj * yk + xk * yi - xi * yk - xj * yi - xk * yj));
%dongunun disinda kullanmak icin bulunan degerleri sakla
varTable = [varTable ; ck,ci,cj,bi,bj,bk,A];
% ekrana uygun formatla yazdiriliyor
fprintf("%d %.2f\t%.2f\t%.2f\t%.2f\t%.2f\t%.2f\t%.2f\t%.2f\t%.2f\t%.2f\t%.2f\t%.2f\t%.4f\t\n",i,xi,xj,xk,yi,yj,yk,ck,ci,cj,bi,bj,bk,A);
end
Kcal= array2table(varTable, 'Variablenames',{'ck','ci','cj','bi','bj','bk','A'}) %K 3x3 matris hesabi icin tablo
%tablodaki degerleri disari cekiyorum sutunlari ayri olarak
%A='Tablename'.'variablename' komutu ile
%n sembolu sondaki sadece karisiklik onlemek amacli
cint=Kcal.ci;
cjnt=Kcal.cj;
cknt=Kcal.ck;
bint=Kcal.bi;
bjnt=Kcal.bj;
bknt=Kcal.bk;
Ant=Kcal.A;
%Transpoze ediyoruz satira cevirmek icin sutun matrisini
bin=bint.';
bjn=bjnt.';
bkn=bknt.';
cin=cint.';
cjn=cjnt.';
ckn=cknt.';
An=Ant.';
%Yazdiracagimiz icin basliklar \t ile tab bosluk biraktirliyor \n ile yeni
%satir yazdiriliyor.
fprintf("e\t K11\tK12\tK13\tK21\tK22\tK23\tK31\tK32\tK33\n"); % tablo basligini formatla
%Degiskenler icin bos matris
varTable2 = [];
for i=1:height(Kcal) % Kcal sutun olarak dizildigi icin height kullandik.
currcin=cin(i);
currcjn=cjn(i);
currckn=ckn(i);
currbin=bin(i);
currbjn=bjn(i);
currbkn=bkn(i);
currAn=An(i);
cinn=cin(i);
cjnn=cjn(i);
cknn=ckn(i);
binn=bin(i);
bjnn=bjn(i);
bknn=bkn(i);
Ann=An(i);
% K Row x Column olarak isimlendirir isek, 3x3 matrisin degerleri
K11=(1/(4 * Ann)) * ( binn^2 + cinn^2);
K12=(1/(4 * Ann)) * ((binn * bjnn) + (cinn * cjnn));
K13=(1/(4 * Ann)) * ((binn * bknn) + (cinn * cknn));
K21=(1/(4 * Ann)) * ((binn * bjnn) + (cinn * cjnn));
K22=(1/(4 * Ann)) * ((bjnn)^2 + (cjnn)^2);
K23=(1/(4 * Ann)) * ((bjnn * bknn) + (cjnn * cknn));
K31=(1/(4 * Ann)) * ((binn * bknn) + (cinn * cknn));
K32=(1/(4 * Ann)) * ((bjnn * bknn) + (cjnn * cknn));
K33=(1/(4 * Ann)) * ((bknn)^2 + (cknn)^2);
varTable2 = [varTable2 ; K11, K12,K13, K21, K22, K23, K31, K32, K33]; %Degiskenleri tabloya aliyoruz
fprintf("%d\t%.4f\t%.4f\t%.4f\t%.4f\t%.4f\t%.4f\t%.4f\t%.4f\t%.4f\t\n",i,K11,K12,K13,K21,K22,K23,K31,K32,K33);
end
Kcal2= array2table(varTable2, 'Variablenames',{'K11','K12','K13','K21','K22','K23','K31', 'K32', 'K33'}) %K 3x3 matris hesabi icin tablo2
%K matrislerini olusturmak icin verilerimiz hazir.
%Siradaki islem bu verileri 3x3 Matrislere yerlestirmek.
% 1.SIRA: K1 K2 K3 VERILERINI AYRI AYRI CEKMEK
for i=1:9
%GEREKSIZ%
varK = [1:1:i]; %1 den i ye kadar rakamlardan olusan satir matrisi yazdirdik
%GEREKSIZ%
K{i}=varTable2(i,:); %1x9 matrisleri K{1} K{2} olarak ayirir row matrisi olarak
end
for i=1:9
KK{i}=reshape(K{i},[3,3]); %% KK{i} adiyla 3x3 matrislerimizi elde ettik.
%%%%%%% cell olarak kaydetmistik bunlari double yapmamiz lazim
%%%% cell ve double arasi fark nedir???
end

Fabio Freschi on 28 Oct 2019
Edited: Fabio Freschi on 28 Oct 2019
If I have understood correctly, you want to store the local 3x3 matrix in the global matrix. To do so you can (I would say must) use the command sparse to create the sparse matrix. sparse accetps as inputs the row, col indices vectors and the values. So the idea is
Inside the loop over elements
1) create the local 3x3 matrix
2) create the indices
3) add them in the global vectors
After the loop, create the global matrix. I can't go through your code, comments are hard to understand. But I can use a code snippet
% initialization (it is possible to correctly preallcoate these vectors, but make the code more complex to understand)
iRowGlo = [];
jColGlo = [];
KGlo = [];
for i = 1:NumOfElems
KLoc = rand(3); % dummy values
% you should also have here the indices of nodes, you have your own technique to get the values
idx = [5 1 6]; % your first picture
% now create the row/col indices
[jColLoc,iRowLoc] = meshgrid(idx,idx); % this makes all the combinations
% load the indices and values in the global vectors
iRowGlo = [iRowGlo; iRowLoc(:)]; % colon makes the vector a column
jColGlo = [jColGlo; jColLoc(:)];
KGlo = [KGlo; KLoc(:)];
end
% now use sparse to assemble the complete matrix
K = sparse(iRowGlo,jColGlo,Kglo,NumOfNodes,NumOfNodes);
note that
1) the matrix is sparse, that is the correct way to store finite element matices
2) for repeated values of iRowGlo jColGlo, sparse sums the values, as erequired by the finite element assembly
3) you don't need to store the local matrix in a cell array, just reuse the variables at each iteration

Anilcan Taner on 28 Oct 2019
It really opened a way for me I will try to implement it to my problem. Thanks
Fabio Freschi on 29 Oct 2019
Let me know if you need more help