To receive notifications about scheduled maintenance, please subscribe to the mailing-list gitlab-operations@sympa.ethz.ch. You can subscribe to the mailing-list at https://sympa.ethz.ch

Commit 51ceece1 authored by dkammer's avatar dkammer
Browse files

add week 6

parent a03b7738
2 # Nel (Number of elements)
2 # Nen (Number of nodes per element)
3 # Nnp (Number of nodes)
2 # n (Number of dofs / degrees of freedom per node)
1 0; 0 0; 1 1 # nodes (node coordinates (x,y))
1 3; 2 3 # elements (Nodes connecting each element)
1 1 # E (Young Modulus of all elements)
1 1 # A (Cross-section area of all elements)
0 0; 0 0; 1 1 # pli (Prescribed Load Index)
0 0; 0 0; 10 0 # pl (Prescribed Load)
1 1; 1 1; 0 0 # pdi (Prescribed Displacement Index)
0 0; 0 0; 0 0 # pd (Prescribed Displacement)
#### Stiffness Method
##### C++ code
###### 1. Compile and Link source files.
In order to compile the source file, just type:
`g++ stiff_cpp_arma.cc -o stiff_cpp_arma.x -std=c++11 -O2 -larmadillo`
If you encounter any linking-related problems, you will probably have to manually add the absolute path to armadillo's include/ and lib/ directories to ${CPLUS_INCLUDE_PATH} and ${LD_LIBRARY_PATH}/${LIBRARY_PATH} respectively.
If armadillo was compiled manually, you need to include the path, e.g.
`g++ stiff_cpp_arma.cc -o stiff_cpp_arma -std=c++11 -O2 -I /path/to/unpacked/armadillo-710.3.0/include -DARMA_DONT_USE_WRAPPER -lopenblas`
###### 2. Execute
Make sure the input file is in the same directory as the executable and just type `./stiff_cpp_arma.x`
#include <iostream> // print to console (cout)
#include <cmath> // math operations
#include <fstream> // file streaming for input file
#include <armadillo> // armadillo library
/* Armadillo syntax:
* arma::mat = arma::Mat<double> -> matrix of doubles.
* arma::imat = arma::Mat<sword> -> matrix of signed long integers (64-bit architecture)*.
* arma::vec = arma::Col<double> -> column vector of doubles.
* arma::ivec = arma::Col<sword> -> column vector of signed long integers.
* arma::uvec = arma::Col<uword> -> column vector of unsigned long integers.
* arma::rowvec = arma::Row<double> -> row vector of doubles.
* * In a 64-bit architecture, a "word" is 64-bit long, hence the typenames "sword" and "uword".
*/
#define LEN 256 // buffer length for input reading
// Struct to store all input parameters from file ( extracted by read_input() function)
struct input_t {
int Nel,Nen,Nnp,n; // Number of elements, nodes-per-element, nodes, dofs per node
arma::mat nodes; // Coordinates of all nodes
arma::imat elements; // Node index pair from every element
arma::vec E, A; // Young modulus and cross section of every element
arma::imat pli, pdi; // constrained dofs matrix (1 -> constrained / 0 -> free)
arma::mat pl, pd; // load and displacement boundary conditions for every dof
};
// Struct to store all vectors and matrices of the 2d truss stiffness model.
struct model_t {
arma::mat K,K_f; // Final stiffness matrix + unconstrained stiffness matrix (floating).
arma::vec f,d,r; // Load + displacement + reaction.
};
// FUNCTION DECLARATIONS
template <typename T>
void read_line(T*, std::ifstream&);
void read_line(int*, std::ifstream&);
input_t read_input(const char* );
void assemble (input_t*, model_t*);
void apply_BC (input_t*, model_t*);
void compute_reaction (input_t*, model_t*);
int main()
{
input_t inp = read_input("2dtruss.inp"); // Read input file
model_t model;
assemble(&inp, &model); // Get unconstrained stiffness matrix + load vector
apply_BC(&inp, &model);
model.d = arma::solve (model.K,model.f); // Solve linear system and get displacement vector
compute_reaction(&inp,&model); // Calculate reaction
std::cout << "K_f:\n" << model.K_f << "\n";
std::cout << "K:\n" << model.K << "\n";
std::cout << "d:\n" << model.d << "\n";
std::cout << "r:\n" << model.r << "\n";
return 0;
}
/* ------------------------------------- /
* ----- FUNCTION DEFINITIONS ---------- /
* -------------------------------------*/
// templated function to read any armadillo object from line string.
// (armadillo constructor from string -> ok for small matrices but not efficient).
template <typename T>
void read_line(T * el, std::ifstream &fin){
char line[LEN]; // Buffer for storing a line from the file
fin.getline(line,LEN,'#'); // Load line from file until '#' character is reached
*el = T(line); // Convert 'line' string to T type (mat, imat, vec, etc)
fin.getline(line,LEN); // Move the file stream pointer to the next line.
}
// read integer from line.
void read_line(int * el, std::ifstream &fin){
char line[LEN];
fin.getline(line,LEN,'#');
*el = atoi(line);
fin.getline(line,LEN);
}
// INPUT FUNCTION : Reads input from file. Takes file name and returns input struct.
// ---------------------------------------------------------------------------------
input_t read_input(const char* file){
std::ifstream fin(file); // Loads file stream
input_t inp; // Creates an input struct to store the data
read_line(&inp.Nel, fin);
read_line(&inp.Nen, fin);
read_line(&inp.Nnp, fin);
read_line(&inp.n, fin);
read_line<arma::mat>(&inp.nodes, fin);
read_line<arma::imat>(&inp.elements, fin);
read_line<arma::vec>(&inp.E, fin);
read_line<arma::vec>(&inp.A, fin);
read_line<arma::imat>(&inp.pli, fin);
read_line<arma::mat>(&inp.pl, fin);
read_line<arma::imat>(&inp.pdi, fin);
read_line<arma::mat>(&inp.pd, fin);
return inp;
}
//---------------------------------------------------------------------------------------
// ASSEMBLY FUNCTION
// Assemble stiffness matrix from input struct containing all parsed input values from input file.
// -----------------------------------------------------------------------------------
void assemble (input_t * inp, model_t * model){
int Nel = inp->Nel;
int Nen = inp->Nen;
int Nnp = inp->Nnp;
int n = inp->n;
arma::vec L(Nel),cx(Nel),cy(Nel); // Store length values for every element
arma::vec A = inp->A; // Cross-section of every element
arma::vec E = inp->E; // Young modulus of every element
arma::mat nodes = inp->nodes; // Nodes matrix
arma::imat elements = inp->elements; // Elements matrix
for (int i=0;i<Nel;++i){
double dx = nodes(elements(i,1)-1,0) - nodes(elements(i,0)-1,0);
double dy = nodes(elements(i,1)-1,1) - nodes(elements(i,0)-1,1);
L(i) = sqrt(dx*dx + dy*dy);
cx(i) = dx / L(i); // Sinus of element angle with respect to global coordinates
cy(i) = dy/ L(i); // Cosinus of element angle with respect to global coordinates
}
arma::imat LM(Nen*n,Nel); // Gather arma::matrix --> Converts local indeces of dofs to global
for (int i=0;i<Nel;++i){
for (int j=0;j<Nen;++j){
for (int k=0;k<n;++k){
LM(2*j+k,i) = k + n * elements(i,j) - 2;
}
}
}
// Floating stiffness matrix (No boundary conditions / )
arma::mat K_f(Nnp*n,Nnp*n,arma::fill::zeros);
arma::mat k_e, rhs; // Defined below
// Assemble floating K
for (int i=0;i<Nel;++i){
arma::mat K_e(Nnp*n,Nnp*n,arma::fill::zeros);
double k = (A(i)*E(i))/L(i);
arma::mat k_e_local( {{k,0,-k,0},{0,0,0,0},{-k,0,k,0},{0,0,0,0}}); // local K for element i
arma::mat R ({{cx(i) ,cy(i) ,0 ,0},{-cy(i) ,cx(i) ,0 ,0},{0 ,0 ,cx(i) ,cy(i)},{0,0,-cy(i),cx(i)}}); // Rotation matrix for element i from local angle to global
rhs = k_e_local* R; // Store right-hand-side
k_e = R.t()*rhs; // Rotated local K of element i (still local indeces)
// R.t() -> Transpose of R. For orthogonal matrices: R^-1 = R^T
for (int j=0;j<Nen*n;++j){
for (int k=0;k<Nen*n;++k){
K_e(LM(j,i),LM(k,i)) = k_e(j,k); // Re-arrange local k_e indices into global K
}
}
K_f = K_f + K_e; // Add the rotated and globally adjusted dofs of element i to global K
}
/* Swap exchanges the content the two objects with each other.
* We could also write `model->K = K` which would copy the matrix named `K`
* into the member `K` of the `model` object, however this would copy the matrix.
* Since we do not need the matrix any longer, we can do this.
*/
model->K_f.swap(K_f);
}
// applies boundary conditions on displacements
void apply_BC (input_t * inp, model_t * model){
int Nnp = inp->Nnp;
int n = inp->n;
arma::mat K (model->K_f); // Start from floating structure
arma::vec f (inp->pl.as_row().as_col()); // Initialize load
arma::uvec pdl = find(inp->pdi.as_row()); // contains 1 for prescribed, 0 for free (displacement dofs)
arma::vec pd (inp->pd.as_row().as_col()); // prescribed displacements.
for (int i=0;i< pdl.n_elem;++i){
f(pdl(i)) = pd(pdl(i));
K.row(pdl(i)) = arma::rowvec (Nnp*n,arma::fill::zeros);
K(pdl(i),pdl(i)) = 1;
}
model->K.swap(K);
model->f.swap(f);
}
// REACTION FUNCTION : Computes reaction to load. Takes floating stiffness matrix
// + displacements + input parameters (prescriptions) and returns the reaction vector.
void compute_reaction (input_t * inp, model_t * model){
arma::vec r = model->K_f * model->d;
arma::uvec pdl = find(inp->pli.as_row());
for (int i=0;i<pdl.n_elem;++i){
r(pdl(i)) = 0;
}
model->r.swap(r);
}
clear;clc;
%% INPUT
% Length parameter
L = 1;
% Young's modulus parameter
E = 1;
% Load parameter
P = 10;
% Cross-section area parameter
A = 1;
% nodes(i,1) and nodes(i,2) are, respectively, the x and y coordinates of node i with respect to the global Cartesian coordiante system
nodes = [L 0; 0 0; L L];
% element(i,1) and elements(i,2) are, repectively, the first and second nodes of element i
elements = [1 3; 2 3];
% E_array(i) is Young's modulus of element i
E_array = [E; E];
% A_array(i) is the cross-section area of element i
A_array = [A; A];
% EA_array(i) is the axial rigidity of element i
EA_array = E_array.*A_array; % .* element-wise multiplication
% PRESCRIBED DISPLACEMENT AND PRESCRIBED LOADS ARRAYS
% Prescribed Load Index
% PLI(i,1) = 1 means fu (concentrated load in the x direction) at node i is prescribed
% PLI(i,1) = 0 means fu (concentrated load in the x direction) at node i is free
% PLI(i,2) = 1 means fv (concentrated load in the x direction) at node i is prescribed
% PLI(i,2) = 0 means fv (concentrated load in the x direction) at node i is free
PLI = [0 0; 0 0; 1 1];
% Prescribed Load values
% PD gives the value of the prescribed concentrated loads that correspond to PDI
PL = [0 0; 0 0; 10 0];
% Prescribed Displacement Index
% PDI(i,1) = 1 means u (displacement in the x direction) at node i is prescribed
% PDI(i,1) = 0 means u (displacement in the x direction) at node i is free
% PDI(i,2) = 1 means v (displacement in the x direction) at node i is prescribed
% PDI(i,2) = 0 means v (displacement in the x direction) at node i is free
PDI = [1 1; 1 1; 0 0];
% Prescribed Displacement values
% PD gives the value of the prescribed displacments that correspond to PDI
PD = [0 0; 0 0; 0 0];
%% PREPROCESSING
% Number of elements and number of element nodes
[Nel,Nen]=size(elements);
% Number of nodal points and number of nodal degrees of freedom
[Nnp,n]=size(nodes);
% Element length and direction cosines vectors
L_array = zeros(Nel,1);
cx_array = zeros(Nel,1); % cos of angle between element (from node 1 to node 2) to global x axis
cy_array = zeros(Nel,1); % cos of angle between element (from node 1 to node 2) to global y axis
for i=1:Nel
dx = nodes(elements(i,2),1) - nodes(elements(i,1),1) ; % Difference between x values of node 2 and node 1 of element i
dy = nodes(elements(i,2),2) - nodes(elements(i,1),2) ; % Difference between y values of node 2 and node 1 of element i
L_array(i) = (dx^2+dy^2)^0.5; % Pythogoras
cx_array(i) = dx/L_array(i);
cy_array(i) = dy/L_array(i);
end
% Vectorize PLI array
PLIV = zeros(Nnp*n,1);
for i =1:Nnp
PLIV(n*i) = PLI(i,2);
PLIV(n*i-1) = PLI(i,1);
end
% Vectorize PL array by turning each row to a column and stacking the
% columns. Can also done with reshape(PLV.',[],1)
PLV = zeros(Nnp*n,1);
for i =1:Nnp
PLV(n*i) = PL(i,2);
PLV(n*i-1) = PL(i,1);
end
% Vectorize PDI array
PDIV = zeros(Nnp*n,1);
for i =1:Nnp
PDIV(n*i) = PDI(i,2);
PDIV(n*i-1) = PDI(i,1);
end
% Vectorize PD array
PDV = zeros(Nnp*n,1);
for i =1:Nnp
PDV(n*i) = PD(i,2);
PDV(n*i-1) = PD(i,1);
end
% Relate local and global (Cartesian) dofs
% LM(k,i) is the global DOF corresponding to the k'th local DOF of element
% i (k ranges from 1 to Nen*n - for a 2d truss that is 4)
LM=zeros(n*Nen,Nel);
for el=1:Nel
for node=1:Nen
LM(node*n-(n-1):node*n,el) = elements(el,node)*n-(n-1):elements(el,node)*n;
end
end
%% ASSEMBLY OF GLOBAL EQUATIONS
% Initialize global stiffness matrix K assuming no prescribed displacements
% (floating structure)
K_float = zeros(Nnp*n);
% Calculate element contributions to K and assemble K
for el=1:Nel
K_el = zeros(Nnp*n);
% Element parameters
EA = EA_array(el);
L = L_array(el);
cx = cx_array(el);
cy = cy_array(el);
% Element stiffness matrix in local coordinates (aligned with element axis)
k_el_local = EA/L* [1 0 -1 0; 0 0 0 0; -1 0 1 0; 0 0 0 0];
% Rotation matrix from element directions to global directions (parallel to global Cartesian axes)
R = [cx cy 0 0; -cy cx 0 0; 0 0 cx cy; 0 0 -cy cx];
% Element stiffness matrix in local Cartesian coordinates (aligned with global Cartesian axes)
k_el = R^-1 * k_el_local * R;
% Expand k_el to the dimensions of K
for i=1:Nen*n
for j=1:Nen*n
K_el(LM(i,el),LM(j,el)) = k_el(i,j); % K_el = expanded element stiffness matrix
end
end
% Assemble K by superposing the K_el
K_float = K_float + K_el;
end
fprintf("K_f: \n");
for i=1:Nnp*n
for j=1:Nnp*n
fprintf("%4.2f\t",K_float(i,j));
end
fprintf("\n");
end
% Account for essential boundary conditions by modifying K_float to K and the RHS vector f
K = K_float; % Start with the "floating structure"
f = PLV; % Initially, f is the external load vector
% Find the indices of DOF's where displacement is prescribed
PDL = find(PDIV==1);
% For those DOF's change the equations to explicitly read the essential
% boundary conditions 1 * d(i) = f(i) = d_prescribed(i)
for i=1:length(PDL)
% Set the i'th RHS entry (i being the DOF index of a presribed displacement
% to the value of the prescribed displacement (zero if there is a support)
f(PDL(i)) = PDV(PDL(i));
% Nullify row (equation) i
K(PDL(i),:) = 0;
% Set to one the (i,i) coefficient
K(PDL(i),PDL(i)) = 1;
% In the end the i'th line reads 1 * d(i) = f(i) = d_prescribed(i)
end
fprintf("K: \n");
for i=1:Nnp*n
for j=1:Nnp*n
fprintf("%4.2f\t",K(i,j));
end
fprintf("\n");
end
%% SOLUTION OF GLOBAL EQUATIONS FOR UNKNOWN DISPLACEMNTS
d = K\f;
%% POST PROCESSING
fprintf("d: \n");
for j=1:Nnp*n
fprintf("%4.2f\t",d(j));
end
fprintf("\n");
% Calculation of reactions
r = K_float*d;
% Set to NaN (not a number) entries of r where external load is prescribed
% (and therefore there are no reactions by definition)
PLL = find(PLIV==1);
for i=1:length(PLL)
r(PLL(i)) = NaN;
end
fprintf("r: \n");
for j=1:Nnp*n
fprintf("%4.2f\t ",r(j));
end
fprintf("\n");
# Assignment
In this project, you will apply your knowledge of classes in C++ on simulation code. Please follow the steps below and create a few slides summarizing your approach and discussion.
The code needed for this project is in the CPP folder. The MATLAB folder is there for your reference. Further, you find documentation on the simulation method in the `doc` folder.
This project builds on the previous project of last week. The overall objective is to modify the model struct (not the input) into a class and use it in the main function.
Go into the `CPP` directory.
1. read the `README.md`
2. read the c++ code and gain SOME understanding of what to expect from the code.
3. compile and execute the code. Does it do what you expected?
4. create an additional model.hpp and model.cpp file
5. implement a class model in these files that is an advanced version of the struct.
6. extend the class to include as many of the functions as member functions as makes sense and you have time to do
7. modify the main function to use an object of type model.
8. compile and execute, and check that it reproduces the same results as before.
9. note what else you could do to improve the model class
Your slides should summarize all of your steps above.
\ No newline at end of file
#### DIRECT STIFFNESS METHOD
The Direct Stiffness Method is the simplest finite-element method. It is used in structural analysis to compute the reactions of a load or force acting on a structure based on the material stiffness properties of all its constitutive elements.
In the end, it only consists in solving a linear system of the form K * d = f, where d is a vector containing the displacements of all degrees of freedom involved, f is athe loading vector and K the stiffness matrix of the structure.
##### SYSTEM
Now, the system that will be solved in this example is a simple 2D truss consisting of just 2 elements and 2 supports, with a fixed horizontal loading on the elbow connecting both elements, as illustrated below:
![2dtruss](2dtruss.png)
##### ALGORITHM
In order to calculate the reactions, we will first assemble the floating stiffness matrix of the structure, apply boundary conditions from the supports, solve the resulting linear system for node displacements and compute the reactions directly as r = K<sub>f</sub> * d, K<sub>f</sub> being the floating or unconstrained stiffness matrix.
###### ASSEMBLY
The most tedious task among all others is the assembly of the floating stiffness matrix K<sub>f</sub>. For this, we need to:
1. Compute the local stiffness matrix of every element k<sub>e</sub> in its local coordinates (x<sup>e</sup>,y<sup>e</sup>).
2. Transform k<sub>e</sub> to global coordinates (x,y) through a rotation R.
k<sub>e</sub> (x,y) = R<sup>T</sup> k<sub>e</sub> R.
3. Assemble local matrices of size (2 x dim, 2 x dim) into the global matrix of size (nodes x dim) by replacing local indeces with global indeces (Gather).
###### BOUNDARY CONDITIONS
Once we have the floating structure, we apply the boundary conditions set by the supports and prescribed loading by setting:
1. All entries in the stiffness matrix that relate to constrained degrees of freedom to zero.
2. Force / loading vector to prescribed values.
###### SOLUTION
For solving the linear system, we use a predefined solver from the armadillo library (linear algebra). However, you could also implement you own solver. A typical linear solver would follow these steps:
1. LU factorization/decomposition i.e. Cholesky :
Factorization into lower and upper triangular matrices L and U.
K = L * U.
2. Forward-substitution algorithm to solve for y = U * d:
L * (y) = f.
3. Backward-substitution algortihm to solve for d:
U * d = y.
###### REACTIONS
Reactions are computed by simply calculating the force from the obtained displacements on the floating matrix:
r = K <sub>f</sub> * d.
#### Lennard-Jones MD code
##### C ++ code
###### 1. Compile
g++ ljmd.cc -o ljmd.x
###### 2. Execute
./ljmd.x
###### 3. INPUT FILES
- ljmc.cc (source file)
- argon\_108.inp (simulation data)
- argon\_108.rest (restart file / initial conditions)
###### 4. OUTPUT FILES
- argon\_108.xyz (positions in (x,y,z) coordinates of every atom in the system for every 'nfi'-th time-step)
- argon\_108.dat ( system data for every 'nfi'-th time-step)
- current time-step
- temperature
- kinetic energy
- potential energy
- total energy
108 # natoms
39.948 # mass in AMU
0.2379 # epsilon in kcal/mol
3.405 # sigma in angstrom
8.5 # rcut in angstrom
17.1580 # box length (in angstrom)
10000 # nr MD steps
5.0 # MD time step (in fs)
100 # output print frequency
6.67103294321331 -10.6146871435653 12.6336939877734
1.06574058650169 -3.33432278188177 -2.59038677851747
-1.78412295775301 -16.5259458407765 4.61680014503288
9.81030288841332 -5.82405963687712 3.65645974729484
15.3247211142637 -4.78053306974291 8.26672732196501
2.11253838645732 -4.72955953073222 6.75738627120688
-4.56209227856753 -3.20329426260617 5.55255473461665
5.86663918205102 -10.4139411857802 1.86416360433577
5.21527405147765 -28.9277726235844 6.19270143781154
13.7051046236876 -4.97120493339350 1.94629699552930
4.19613239880890 -14.9790822892036 16.5333880558744
16.7521502712298 -10.0804241645332 14.5028498925954
12.6350229329042 -1.40635636262803 2.38627163152721
-1.08439966554446 2.41557722385393 8.37791671838088
8.96146834305310 -19.0272834036820 4.05806985477701
11.1459036450055 3.18552372924961 6.723294858788574E-002
16.0690989728094 -1.03744620980343 18.5655733607861
18.6984479366938 -17.6491643729990 3.60417710784249
-3.13161415889863 -6.75371567663538 5.00611042020118
17.4072659061416 -7.88472086048433 7.56603710600785
19.6499870844957 -10.2813545329547 3.47269993677700
9.14358801056816 10.5769423320094 11.0434815971320
12.3852978120607 -11.0158733347720 20.0337606688683
6.18023614761844 -4.34746306558901 5.42729289210231
0.436403345315857 -1.36655907890507 7.37333479343561
-0.920905510346089 -3.30961460032849 21.5616988921137
3.31118739595332 -20.7760460532019 20.3172712920931
0.507105497066932 -21.6567949737129 0.933840061914709
14.4213709093942 -20.3431287076219 15.4309530884964
5.56191131389967 -7.25872501357231 12.4568959982054
3.06146523952825 -14.4087233421704 6.69325558343884
-1.37943237515652 -10.1053635146924 4.64502788848255