Difference between revisions of "Metnum03-Farhan Aditya Wibowo"

From ccitonlinewiki
Jump to: navigation, search
(Tugas Besar)
(Tugas Besar)
Line 147: Line 147:
 
berikut adalah coding yang dipakai
 
berikut adalah coding yang dipakai
  
 +
model Trusses_3D_Tugas_Besar_Safety 
 
//define initial variable
 
//define initial variable
 
parameter Integer Points=size(P,1); //Number of Points
 
parameter Integer Points=size(P,1); //Number of Points
 
parameter Integer Trusses=size(C,1); //Number of Trusses
 
parameter Integer Trusses=size(C,1); //Number of Trusses
parameter Real Yield=215e6; //Yield Strength (Pa)
+
parameter Real Yield= (nilai yield) ; //Yield Strength Material
parameter Real Area=0.000224;  //Area L Profile (Dimension=0.03, Thickness=0,004) (m2)
+
parameter Real Area= (nilai area) ;  //Luas Siku
parameter Real Elas=193e9;    //Elasticity SS 304  (Pa)
+
parameter Real Elas= (nilai elastisitas) ;    //Elasticity Material
 
 
 
//define connection
 
//define connection
parameter Integer C[:,2]=[1,5;  
+
parameter Integer C[Trusses,2]=[1,5; //vertical 1st floor
                                 2,6;
+
                                 2,6; //vertical 1st floor
                                 3,7;
+
                                 3,7; //vertical 1st floor
                                 4,8;
+
                                 4,8; //vertical 1st floor
                                 5,6;  //1st floor
+
                                 5,6;  //horizontal 1st floor
                                 6,7;  //1st floor
+
                                 6,7;  //horizontal 1st floor
                                 7,8;  //1st floor
+
                                 7,8;  //horizontal 1st floor
                                 5,8;  //1st floor
+
                                 5,8;  //horizontal 1st floor
                                 5,9;
+
                                 5,9; //vertical 2st floor
                                 6,10;
+
                                 6,10; //vertical 2st floor
                                 7,11;
+
                                 7,11; //vertical 2st floor
                                 8,12;
+
                                 8,12; //vertical 2st floor
                                 9,10; //2nd floor
+
                                 9,10; //horizontal 2st floor
                                 10,11;//2nd floor  
+
                                 10,11; //horizontal 2st floor
                                 11,12;//2nd floor
+
                                 11,12; //horizontal 2st floor
                                 9,12; //2nd floor
+
                                 9,12; //horizontal 2st floor
                                 9,13;
+
                                 9,13; //vertical 3st floor
                                 10,14;
+
                                 10,14; //vertical 3st floor
                                 11,15;
+
                                 11,15; //vertical 3st floor
                                 12,16;
+
                                 12,16; //vertical 3st floor
                                 13,14;//3rd floor
+
                                 13,14; //horizontal 3st floor
                                 14,15;//3rd floor
+
                                 14,15; //horizontal 3st floor
                                 15,16;//3rd floor
+
                                 15,16; //horizontal 3st floor
                                 13,16];//3rd floor
+
                                 13,16]; //horizontal 3st floor
                                                             
+
                                                           
 
//define coordinates (please put orderly)
 
//define coordinates (please put orderly)
parameter Real P[:,6]=[0.3,-0.375,0,1,1,1;     //1
+
parameter Real P[Points,3]=[ 0   ,0 ,0,1,1,1; //1
                            -0.3,-0.375,0,1,1,1;   //2
+
                          0.75,0 ,0,1,1,1; //2
                            -0.3,0.375,0,1,1,1;     //3
+
                          0.75,0.6,0,1,1,1; //3
                            0.3,0.375,0,1,1,1;     //4
+
                          0   ,0.6,0,1,1,1; //4
                           
+
                          0   ,0 ,0.3,0,0,0; //5
                            0.3,-0.375,0.6,0,0,0;   //5
+
                          0.75,0 ,0.3,0,0,0;  //6
                            -0.3,-0.375,0.6,0,0,0;  //6
+
                          0.75,0.6,0.3,0,0,0; //7
                            -0.3,0.375,0.6,0,0,0;   //7
+
                          0   ,0.6,0.3,0,0,0; //8
                            0.3,0.375,0.6,0,0,0;   //8
+
                          0   ,0 ,1.05,0,0,0; //9
                           
+
                          0.75,0 ,1.05,0,0,0;  //10   
                            0.3,-0.375,1.2,0,0,0;   //9
+
                          0.75,0.6,1.05,0,0,0; //11
                            -0.3,-0.375,1.2,0,0,0;  //10   
+
                          0   ,0.6,1.05,0,0,0; //12
                            -0.3,0.375,1.2,0,0,0;   //11
+
                          0   ,0 ,1.8,0,0,0; //13
                            0.3,0.375,1.2,0,0,0;   //12
+
                          0.75,0 ,1.8,0,0,0;  //14
                           
+
                          0.75,0.6,1.8,0,0,0; //15
                            0.3,-0.375,1.8,0,0,0;   //13
+
                          0   ,0.6,1.8,0,0,0]; //16
                            -0.3,-0.375,1.8,0,0,0;  //14
 
                            -0.3,0.375,1.8,0,0,0;   //15
 
                            0.3,0.375,1.8,0,0,0];   //16
 
 
                              
 
                              
 
//define external force (please put orderly)
 
//define external force (please put orderly)
Line 217: Line 214:
 
                             0,0,-1000,  
 
                             0,0,-1000,  
 
                             0,0,-1000,  
 
                             0,0,-1000,  
                             0,0,-500};  
+
                             0,0,-500};
  
 +
//define external force (please put orderly)
 +
parameter Real F[Points*3]={0,0,0,
 +
                          0,0,0,
 +
                          0,0,0,
 +
                          0,0,0,
 +
                          0,0,0,
 +
                          0,0,0,
 +
                          0,0,0,
 +
                          0,0,0,
 +
                          0,0,0,
 +
                          0,0,0,
 +
                          0,0,0,
 +
                          0,0,0,
 +
                          0,0,-1000,
 +
                          0,0,-500,
 +
                          0,0,-500,
 +
                          0,0,-1000};
 
//solution
 
//solution
 
Real displacement[N], reaction[N];
 
Real displacement[N], reaction[N];
 
Real check[3];
 
Real check[3];
 
 
Real stress1[Trusses];
 
Real stress1[Trusses];
 
Real safety[Trusses];
 
Real safety[Trusses];
 
Real dis[3];
 
Real dis[3];
 
Real Str[3];
 
Real Str[3];
 
 
protected
 
protected
parameter Integer N=3*Points;
+
parameter Integer N=3*Points;
Real q1[3], q2[3], g[N,N], G[N,N], G_star[N,N], id[N,N]=identity(N), cx, cy, cz, L, X[3,3];
+
Real q1[3], q2[3], g[N,N], G[N,N], G_star[N,N], id[N,N]=identity(N), cx, cy, cz, L, X[3,3];
Real err=10e-10, ers=10e-6;
+
Real err=10e-15, ers=10e-8;
 
 
 
algorithm
 
algorithm
 
//Creating Global Matrix
 
//Creating Global Matrix
Line 238: Line 249:
 
for i in 1:Trusses loop
 
for i in 1:Trusses loop
 
  for j in 1:3 loop
 
  for j in 1:3 loop
  q1[j]:=P[C[i,1],j];
+
  q1[j]:=P[C[i,1],j];
  q2[j]:=P[C[i,2],j];
+
  q2[j]:=P[C[i,2],j];
  end for;
+
  end for;      
     
+
  //Solving Matrix
    //Solving Matrix
+
  L:=Modelica.Math.Vectors.length(q2-q1);
    L:=Modelica.Math.Vectors.length(q2-q1);
+
  cx:=(q2[1]-q1[1])/L;
    cx:=(q2[1]-q1[1])/L;
+
  cy:=(q2[2]-q1[2])/L;
    cy:=(q2[2]-q1[2])/L;
+
  cz:=(q2[3]-q1[3])/L;  
    cz:=(q2[3]-q1[3])/L;  
+
  X:=(Area*Elas/L)*[cx^2,cx*cy,cx*cz;
    X:=(Area*Elas/L)*[cx^2,cx*cy,cx*cz;
+
                    cy*cx,cy^2,cy*cz;
                      cy*cx,cy^2,cy*cz;
+
                    cz*cx,cz*cy,cz^2];
                      cz*cx,cz*cy,cz^2];
 
 
 
     //Transforming to global matrix
 
     //Transforming to global matrix
    g:=zeros(N,N);  
+
  g:=zeros(N,N);  
    for m,n in 1:3 loop
+
  for m,n in 1:3 loop
      g[3*(C[i,1]-1)+m,3*(C[i,1]-1)+n]:=X[m,n];
+
    g[3*(C[i,1]-1)+m,3*(C[i,1]-1)+n]:=X[m,n];
      g[3*(C[i,2]-1)+m,3*(C[i,2]-1)+n]:=X[m,n];
+
    g[3*(C[i,2]-1)+m,3*(C[i,2]-1)+n]:=X[m,n];
      g[3*(C[i,2]-1)+m,3*(C[i,1]-1)+n]:=-X[m,n];
+
    g[3*(C[i,2]-1)+m,3*(C[i,1]-1)+n]:=-X[m,n];
      g[3*(C[i,1]-1)+m,3*(C[i,2]-1)+n]:=-X[m,n];
+
    g[3*(C[i,1]-1)+m,3*(C[i,2]-1)+n]:=-X[m,n];
    end for;
+
  end for;  
 
 
  G_star:=G+g;
 
  G_star:=G+g;
 
  G:=G_star;
 
  G:=G_star;
 
end for;
 
end for;
 
 
//Implementing boundary
 
//Implementing boundary
 
for x in 1:Points loop
 
for x in 1:Points loop
  if P[x,4] <> 0 then
+
if P[x,4] <> 0 then
    for a in 1:Points*3 loop
+
  for a in 1:Points*3 loop
      G[(x*3)-2,a]:=0;
+
    G[(x*3)-2,a]:=0;
      G[(x*3)-2,(x*3)-2]:=1;
+
    G[(x*3)-2,(x*3)-2]:=1;
    end for;
+
  end for;
  end if;
+
end if;
  if P[x,5] <> 0 then
+
if P[x,5] <> 0 then
    for a in 1:Points*3 loop
+
  for a in 1:Points*3 loop
      G[(x*3)-1,a]:=0;
+
    G[(x*3)-1,a]:=0;
      G[(x*3)-1,(x*3)-1]:=1;
+
    G[(x*3)-1,(x*3)-1]:=1;
    end for;
+
  end for;
  end if;
+
end if;
  if P[x,6] <> 0 then
+
if P[x,6] <> 0 then
    for a in 1:Points*3 loop
+
  for a in 1:Points*3 loop
      G[x*3,a]:=0;
+
    G[x*3,a]:=0;
      G[x*3,x*3]:=1;
+
    G[x*3,x*3]:=1;
    end for;
+
  end for;
  end if;
+
end if;
 
end for;
 
end for;
 
 
//Solving displacement
 
//Solving displacement
 
displacement:=Modelica.Math.Matrices.solve(G,F);
 
displacement:=Modelica.Math.Matrices.solve(G,F);
 
 
//Solving reaction
 
//Solving reaction
 
reaction:=(G_star*displacement)-F;
 
reaction:=(G_star*displacement)-F;
 
 
//Eliminating float error
 
//Eliminating float error
 
for i in 1:N loop
 
for i in 1:N loop
Line 297: Line 301:
 
  displacement[i]:=if abs(displacement[i])<=err then 0 else displacement[i];
 
  displacement[i]:=if abs(displacement[i])<=err then 0 else displacement[i];
 
end for;
 
end for;
 
 
//Checking Force
 
//Checking Force
 
check[1]:=sum({reaction[i] for i in (1:3:(N-2))})+sum({F[i] for i in (1:3:(N-2))});
 
check[1]:=sum({reaction[i] for i in (1:3:(N-2))})+sum({F[i] for i in (1:3:(N-2))});
 
check[2]:=sum({reaction[i] for i in (2:3:(N-1))})+sum({F[i] for i in (2:3:(N-1))});
 
check[2]:=sum({reaction[i] for i in (2:3:(N-1))})+sum({F[i] for i in (2:3:(N-1))});
check[3]:=sum({reaction[i] for i in (3:3:N)})+sum({F[i] for i in (3:3:N)});
+
check[3]:=sum({reaction[i] for i in (3:3:N)})+sum({F[i] for i in (3:3:N)});
 
 
 
for i in 1:3 loop
 
for i in 1:3 loop
  check[i] := if abs(check[i])<=ers then 0 else check[i];
+
check[i] := if abs(check[i])<=ers then 0 else check[i];
 
end for;
 
end for;
 
 
//Calculating stress in each truss
 
//Calculating stress in each truss
 
for i in 1:Trusses loop
 
for i in 1:Trusses loop
for j in 1:3 loop
+
for j in 1:3 loop
  q1[j]:=P[C[i,1],j];
+
  q1[j]:=P[C[i,1],j];
  q2[j]:=P[C[i,2],j];
+
  q2[j]:=P[C[i,2],j];
  dis[j]:=abs(displacement[3*(C[i,1]-1)+j]-displacement[3*(C[i,2]-1)+j]);
+
  dis[j]:=abs(displacement[3*(C[i,1]-1)+j]-displacement[3*(C[i,2]-1)+j]);
end for;
+
end for;      
     
+
  //Solving Matrix
    //Solving Matrix
+
  L:=Modelica.Math.Vectors.length(q2-q1);
    L:=Modelica.Math.Vectors.length(q2-q1);
+
  cx:=(q2[1]-q1[1])/L;
    cx:=(q2[1]-q1[1])/L;
+
  cy:=(q2[2]-q1[2])/L;
    cy:=(q2[2]-q1[2])/L;
+
  cz:=(q2[3]-q1[3])/L;  
    cz:=(q2[3]-q1[3])/L;  
+
  X:=(Elas/L)*[cx^2,cx*cy,cx*cz;
    X:=(Elas/L)*[cx^2,cx*cy,cx*cz;
+
                cy*cx,cy^2,cy*cz;
                cy*cx,cy^2,cy*cz;
+
                cz*cx,cz*cy,cz^2];  
                cz*cx,cz*cy,cz^2];
+
  Str:=(X*dis);
   
+
  stress1[i]:=Modelica.Math.Vectors.length(Str);
    Str:=(X*dis);
 
    stress1[i]:=Modelica.Math.Vectors.length(Str);
 
 
end for;
 
end for;
 
 
//Safety factor
 
//Safety factor
 
for i in 1:Trusses loop
 
for i in 1:Trusses loop
  if stress1[i]>0 then
+
if stress1[i]>0 then
    safety[i]:=Yield/stress1[i];
+
  safety[i]:=Yield/stress1[i];
  else
+
else
    safety[i]:=0;
+
  safety[i]:=0;
  end if;  
+
end if;  
 
end for;
 
end for;
 
 
end Trusses_3D_Tugas_Besar_Safety;
 
end Trusses_3D_Tugas_Besar_Safety;

Revision as of 13:06, 4 January 2021

بِسْمِ اللَّهِ الرَّحْمَنِ الرَّحِيم

dengan ini saya mengisi halaman ini untuk memenuhi mata kuliah metode numerik

Bio Data

Nama  : Farhan Aditya Wibowo

NPM  : 1706024665

Program Studi : S1 Teknik Mesin Parallel

Pertemuan 1 (19 November 2020)

Pada kesempatan pertemuan sebelumnya kelas kami diberikan fasilitas untuk berdiskusi untuk membahas hal apa saja yang telah didapat dari kelas metode numerik dan pengaplikasian dari yang telah didapat, dari situ saya meraskan bahwa pada kehidupan sehari-hari kitapun harus ada suatu target seperti layaknya kita mengerjakan metode numerik dengan target tugas selesai dan mengerti apa yang di beri tahu oleh orang yang lebih berpengalaman dalam bidang tersebut. Lalu selanjutnya kami membahas tentang konsep tak hingga yang sebenarnya bukan bilangan, melainkan hanya sebuah konsep dibuat manusia karena manusia memiliki batasan dan hanya maha penciptalah yang tahu akan tak hingga tersebut. Pada sesi selanjutnya kami diberi informasi bahwa pemberian materi akan berupa perangkat lunak Open Modelica. Open Modelica merupakan penyimulasi sistem dengan data input kode untuk melakukan suatu penelitian terhadap sistem tersebut.

Pertemuan 2 (16 November 2020)

Membuat suatu program untuk menyelesaikan persamaan aljabar simultan, berikut merupakan program yang saya buat: Screen Shot 2020-11-29 at 18.18.27.png

Dengan metode gauss elimination dan function solve yang disediakan oleh software modelica, saya mendapatkan hasil seperti berikut :

Tugasmetnumbowo2.png

Pertemuan 3 ( 23 November 2020)

Pada pertemuan ini kami membuktikan studi dari buku Numerical Methods for Engineers 7th Edition

Tugasmetnumbowo3.png

dari soal diketahui matriks [K][X]=[W]

komponen [K] dan [W] diketahui, lalu mencari komponen [X]

Pada Open Modelica saya menggunakan metode Gauss Elimination karena pada cara manual kita seharusnya mencari variable dari Matrik K dan W.

Tugasmetnumbowo4.png

dengan begitu ini hasil dari yang saya kerjakan

Tugasmetnumbowo5.png

Pertemuan 3 example 2.1(23 November 2020)

Di pertemuan hari ini, Pak Dai menjelaskan tentang aplikasi metode numerik untuk permasalahan-permasalahan teknik. Salah satunya adalah permasalahan sistem pegas-massa.

Setelah kelas, Pak Dai memberikan tugas untuk menyelesaikan soal berikut:

Tugasmetnumbowo16.png

Untuk menyelesaikan soal ini perlu dilakukan pengelompokan menjadi node dan elemen seperti pada tabel berikut:

Tugasmetnumbowo17.png

lalu perlu dilakukan perhitungan nilai kekakuan pada elemen. Untuk elemen 1,3,5, dan 6 nilai kekakuannya adalah 4,22 x 10^5 lb/in. sedangkan untuk elemen 2, dan 4 nilai kekakuannya adalah 2,98 x 10^5 lb/in.

setelah itu perlu dilakukan analisis kekakuan pada tiap elemen dalam matriks koordinat global, kemudian dijumlahkan untuk mendapatkan K global. berikut adalah hasil penjumlahan dari nilai kekakuan tiap elemen:

Tugasmetnumbowo18.png

disederhanakan menjadi

Tugasmetnumbowo9.png

setelah mendapat matriks kekakuan, diterapkan kondisi batas dan beban. Untuk node 1 dan 3 adalah fixed, maka U1X=0, U1Y=0, U3X=0, U3Y=0 dan beban diberikan pada node 4 dan 5 sebesar F4Y= -500lb dan F5Y= -500lb

dengan menggunakan Hukum Hooke yaitu F=k.x, maka didapatkan persamaan defleksi menjadi

Tugasmetnumbowo10.png

Persamaan diatas kemudian diselesaikan dengan OpenModelica, berikut saya lampirkan coding yang saya buat

Tugasmetnumbowo11.png

setelah melakukan pengecekan dan simulasi, saya melakukan plotting. berikut hasil plotting simulasi tersebut.

Tugasmetnumbowo12.png

hasil dapat dilihat pada panel di bagian kanan bawah gambar.


Untuk mendapatkan gaya reaksi, persamaannya adalah:

{R}=[K].{U}-{F}

diterapkan matriks-matriks yang sudah diketahui, didapat:

Tugasmetnumbowo13.png

Persamaan diatas dapat diselesaikan oleh openmodelica. berikut adalah coding yang saya buat

Tugasmetnumbowo14.png

hasil simulasi tersebut di plot sebagai berikut:

Tugasmetnumbowo15.png

nilai R tiap node dapat dilihat di bagian kanan bawah gambar

Pertemuan 4 (30 November 2020)

Kuismetnumbowo1.jpg

Kuismetnumbowo2.jpg

Pertemuan 14 Desember 2020

Aplikasi Metode Numerik dalam Kasus Optimasi

Tugas Besar

Objektif

Mencari cost yang paling optimal pada problem di bawah ini

Tugas Besar Metnum Geometri Jos.jpg

Constraint

- Spesifikasi L (Panjang) dan geometri rangka truss

- Gaya beban terhadap struktur (1000 N dan 2000 N)

Metode Pengerjaan

- Mencari data material

- Mencari beberapa parameter seperti stress dan defleksi

- Mencari density atau yield untuk dimensi material yang dibutuhkan

- Memilih material yang cocok

Asumsi

- Variasi Stiffness terikat dengan variabel area. Memvariasikan Elastisitas tergolong sulit karena setiap material memiliki range yang tidak teratur dan dalam satu material yang sejenis (struktur biaya tetap) tidak terjadi perubahan nilai elastisitas yang berbanding lurus dengan perubahan biaya.

- Beban akan terdistribusi hanya pada point penghubung (karena bersifat truss)

Data Material

material yang ditinjau adalah ss400, ss316, ss304

Screen Shot 2021-01-04 at 12.52.07.png

Mencari Parameter Stress dan Defleksi

mencari stress dan defleksi menggunakan Modeling Trusses

berikut adalah coding yang dipakai

model Trusses_3D_Tugas_Besar_Safety //define initial variable parameter Integer Points=size(P,1); //Number of Points parameter Integer Trusses=size(C,1); //Number of Trusses parameter Real Yield= (nilai yield) ; //Yield Strength Material parameter Real Area= (nilai area) ; //Luas Siku parameter Real Elas= (nilai elastisitas) ; //Elasticity Material //define connection parameter Integer C[Trusses,2]=[1,5; //vertical 1st floor

                               2,6;  //vertical 1st floor
                               3,7;  //vertical 1st floor
                               4,8;  //vertical 1st floor
                               5,6;  //horizontal 1st floor
                               6,7;  //horizontal 1st floor
                               7,8;  //horizontal 1st floor
                               5,8;  //horizontal 1st floor
                               5,9;  //vertical 2st floor
                               6,10;  //vertical 2st floor
                               7,11;  //vertical 2st floor
                               8,12;  //vertical 2st floor
                               9,10;  //horizontal 2st floor
                               10,11;  //horizontal 2st floor
                               11,12;  //horizontal 2st floor
                               9,12;  //horizontal 2st floor
                               9,13;  //vertical 3st floor
                               10,14;  //vertical 3st floor
                               11,15;  //vertical 3st floor
                               12,16;  //vertical 3st floor
                               13,14;  //horizontal 3st floor
                               14,15;  //horizontal 3st floor
                               15,16;  //horizontal 3st floor
                               13,16];  //horizontal 3st floor
                                                            

//define coordinates (please put orderly) parameter Real P[Points,3]=[ 0 ,0 ,0,1,1,1; //1

                          0.75,0  ,0,1,1,1;  //2
                          0.75,0.6,0,1,1,1;  //3
                          0   ,0.6,0,1,1,1;  //4
                          0   ,0  ,0.3,0,0,0;  //5
                          0.75,0  ,0.3,0,0,0;  //6
                          0.75,0.6,0.3,0,0,0;  //7
                          0   ,0.6,0.3,0,0,0;  //8
                          0   ,0  ,1.05,0,0,0;  //9
                          0.75,0  ,1.05,0,0,0;  //10  
                          0.75,0.6,1.05,0,0,0;  //11
                          0   ,0.6,1.05,0,0,0;  //12
                          0   ,0  ,1.8,0,0,0;  //13
                          0.75,0  ,1.8,0,0,0;  //14
                          0.75,0.6,1.8,0,0,0;  //15
                          0   ,0.6,1.8,0,0,0];  //16
                           

//define external force (please put orderly) parameter Real F[Points*3]={0,0,0,

                           0,0,0, 
                           0,0,0, 
                           0,0,0, 
                           0,0,0, 
                           0,0,0, 
                           0,0,0, 
                           0,0,0, 
                           0,0,0, 
                           0,0,0, 
                           0,0,0, 
                           0,0,0, 
                           0,0,-500, 
                           0,0,-1000, 
                           0,0,-1000, 
                           0,0,-500};

//define external force (please put orderly) parameter Real F[Points*3]={0,0,0,

                          0,0,0, 
                          0,0,0, 
                          0,0,0, 
                          0,0,0, 
                          0,0,0, 
                          0,0,0, 
                          0,0,0, 
                          0,0,0, 
                          0,0,0, 
                          0,0,0, 
                          0,0,0, 
                          0,0,-1000, 
                          0,0,-500, 
                          0,0,-500, 
                          0,0,-1000}; 

//solution Real displacement[N], reaction[N]; Real check[3]; Real stress1[Trusses]; Real safety[Trusses]; Real dis[3]; Real Str[3]; protected

parameter Integer N=3*Points;
Real q1[3], q2[3], g[N,N], G[N,N], G_star[N,N], id[N,N]=identity(N), cx, cy, cz, L, X[3,3];
Real err=10e-15, ers=10e-8;

algorithm //Creating Global Matrix G:=id; for i in 1:Trusses loop

for j in 1:3 loop
 q1[j]:=P[C[i,1],j];
 q2[j]:=P[C[i,2],j];
end for;       
  //Solving Matrix
  L:=Modelica.Math.Vectors.length(q2-q1);
  cx:=(q2[1]-q1[1])/L;
  cy:=(q2[2]-q1[2])/L;
  cz:=(q2[3]-q1[3])/L; 
  X:=(Area*Elas/L)*[cx^2,cx*cy,cx*cz;
                    cy*cx,cy^2,cy*cz;
                    cz*cx,cz*cy,cz^2];
   //Transforming to global matrix
  g:=zeros(N,N); 
  for m,n in 1:3 loop
    g[3*(C[i,1]-1)+m,3*(C[i,1]-1)+n]:=X[m,n];
    g[3*(C[i,2]-1)+m,3*(C[i,2]-1)+n]:=X[m,n];
    g[3*(C[i,2]-1)+m,3*(C[i,1]-1)+n]:=-X[m,n];
    g[3*(C[i,1]-1)+m,3*(C[i,2]-1)+n]:=-X[m,n];
  end for;   
G_star:=G+g;
G:=G_star;

end for; //Implementing boundary for x in 1:Points loop

if P[x,4] <> 0 then
  for a in 1:Points*3 loop
    G[(x*3)-2,a]:=0;
    G[(x*3)-2,(x*3)-2]:=1;
  end for;
end if;
if P[x,5] <> 0 then
  for a in 1:Points*3 loop
    G[(x*3)-1,a]:=0;
    G[(x*3)-1,(x*3)-1]:=1;
  end for;
end if;
if P[x,6] <> 0 then
  for a in 1:Points*3 loop
    G[x*3,a]:=0;
    G[x*3,x*3]:=1;
  end for;
end if;

end for; //Solving displacement displacement:=Modelica.Math.Matrices.solve(G,F); //Solving reaction reaction:=(G_star*displacement)-F; //Eliminating float error for i in 1:N loop

reaction[i]:=if abs(reaction[i])<=err then 0 else reaction[i];
displacement[i]:=if abs(displacement[i])<=err then 0 else displacement[i];

end for; //Checking Force check[1]:=sum({reaction[i] for i in (1:3:(N-2))})+sum({F[i] for i in (1:3:(N-2))}); check[2]:=sum({reaction[i] for i in (2:3:(N-1))})+sum({F[i] for i in (2:3:(N-1))}); check[3]:=sum({reaction[i] for i in (3:3:N)})+sum({F[i] for i in (3:3:N)}); for i in 1:3 loop

check[i] := if abs(check[i])<=ers then 0 else check[i];

end for; //Calculating stress in each truss for i in 1:Trusses loop for j in 1:3 loop

 q1[j]:=P[C[i,1],j];
 q2[j]:=P[C[i,2],j];
 dis[j]:=abs(displacement[3*(C[i,1]-1)+j]-displacement[3*(C[i,2]-1)+j]);

end for;

  //Solving Matrix
  L:=Modelica.Math.Vectors.length(q2-q1);
  cx:=(q2[1]-q1[1])/L;
  cy:=(q2[2]-q1[2])/L;
  cz:=(q2[3]-q1[3])/L; 
  X:=(Elas/L)*[cx^2,cx*cy,cx*cz;
               cy*cx,cy^2,cy*cz;
               cz*cx,cz*cy,cz^2];    
  Str:=(X*dis);
  stress1[i]:=Modelica.Math.Vectors.length(Str);

end for; //Safety factor for i in 1:Trusses loop

if stress1[i]>0 then
  safety[i]:=Yield/stress1[i];
else
  safety[i]:=0;
end if; 

end for; end Trusses_3D_Tugas_Besar_Safety;