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 3034f819 authored by dkammer's avatar dkammer
Browse files

add week 6.5

parent 51ceece1
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)
#include <iostream> // print to console (cout)
#include <cmath> // math operations
#include "Model.hpp"
#define LEN 256 // buffer length for input reading
// 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 Model::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 Model::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.
// ---------------------------------------------------------------------------------
void Model::read_input(const char* file){
std::ifstream fin(file); // Loads file stream
read_line(&Nel, fin);
read_line(&Nen, fin);
read_line(&Nnp, fin);
read_line(&n, fin);
read_line<arma::mat>(&nodes, fin);
read_line<arma::imat>(&elements, fin);
read_line<arma::vec>(&E, fin);
read_line<arma::vec>(&A, fin);
read_line<arma::imat>(&pli, fin);
read_line<arma::mat>(&pl, fin);
read_line<arma::imat>(&pdi, fin);
read_line<arma::mat>(&pd, fin);
}
//---------------------------------------------------------------------------------------
// ASSEMBLY FUNCTION
// Assemble stiffness matrix from input struct containing all parsed input values from input file.
// -----------------------------------------------------------------------------------
void Model::assemble (){
arma::vec L(Nel),cx(Nel),cy(Nel); // Store length values for every element
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* p_k_f = & K_f;
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.
*/
(*p_k_f).swap(K_f);
}
// applies boundary conditions on displacements
void Model::apply_BC (){
arma::mat* p_k = & K;
arma::mat* p_f = & f;
arma::mat* p_pd = & pd;
arma::mat K (K_f); // Start from floating structure
arma::vec f (pl.as_row().as_col()); // Initialize load
arma::uvec pdl = find(pdi.as_row()); // contains 1 for prescribed, 0 for free (displacement dofs)
arma::vec pd ((*p_pd).as_row().as_col()); // prescribed displacements.
for (int i=0;(unsigned)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;
}
(*p_k).swap(K);
(*p_f).swap(f);
}
// REACTION FUNCTION : Computes reaction to load. Takes floating stiffness matrix
// + displacements + input parameters (prescriptions) and returns the reaction vector.
void Model::compute_reaction (){
arma::vec* p_r = & r;
arma::vec r = K_f * d;
arma::uvec pdl = find(pli.as_row());
for (int i=0;(unsigned)i<pdl.n_elem;++i){
r(pdl(i)) = 0;
}
(*p_r).swap(r);
}
void Model::solve()
{
*(& d) = arma::solve (K,f); // Solve linear system and get displacement vector
}
void Model::print_results()
{
std::cout << "K_f:\n" << K_f << "\n";
std::cout << "K:\n" << K << "\n";
std::cout << "d:\n" << d << "\n";
std::cout << "r:\n" << r << "\n";
}
#ifndef MODELHEADERDEF
#define MODELHEADERDEF
#include <fstream> // file streaming for input file
#include <armadillo> // armadillo library
class Model
{
public:
template <typename T>
void read_line(T*, std::ifstream&);
void read_line(int*, std::ifstream&);
void read_input(const char* );
void assemble ();
void apply_BC ();
void solve ();
void compute_reaction ();
void print_results();
private:
// Struct to store all input parameters from file ( extracted by read_input() function)
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.
arma::mat K,K_f; // Final stiffness matrix + unconstrained stiffness matrix (floating).
arma::vec f,d,r; // Load + displacement + reaction.
};
#endif
\ No newline at end of file
#### 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 "Model.hpp"
int main()
{
Model my_model;
my_model.read_input("2dtruss.inp");
my_model.assemble(); // Get unconstrained stiffness matrix + load vector
my_model.apply_BC();
my_model.solve(); // Solve linear system and get displacement vector
my_model.compute_reaction(); // Calculate reaction
my_model.print_results(); // print results
return 0;
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment