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

week 7: group 1 code

parent 3034f819
cmake_minimum_required(VERSION 3.1)
project(tehpc)
enable_language(CXX)
# flags for compilation (options)
set(CMAKE_CXX_FLAGS "-Wall -Wextra -pedantic -std=c++11 " CACHE STRING "" FORCE)
set(CMAKE_CXX_FLAGS_RELEASE "-O3 -DNDEBUG " CACHE STRING "" FORCE )
# alternative to add_library (but with option)
set(BUILD_SHARED_LIBS ON CACHE BOOL "Build shared libraries.")
mark_as_advanced(BUILD_SHARED_LIBS)
# build type
set(CMAKE_BUILD_TYPE "${CMAKE_BUILD_TYPE}" CACHE STRING
"Build options: None Debug Release RelWithDebInfo MinSizeRel."
FORCE )
# include directories
set(TEHPC_INCLUDE_DIRS
"${CMAKE_CURRENT_BINARY_DIR}/src"
"${CMAKE_CURRENT_SOURCE_DIR}/src"
)
include_directories(${TEHPC_INCLUDE_DIRS})
find_package(Armadillo REQUIRED)
include_directories(${ARMADILLO_INCLUDE_DIRS})
# source directory
add_subdirectory(src)
# ------------------------------------------------------
# examples
# define add_simulation function
function(add_simulation SIM_EXE)
add_executable(${SIM_EXE} ${SIM_EXE}.cpp)
target_link_libraries(${SIM_EXE} tehpc)
endfunction()
option(TEHPC_EXAMPLES "examples" OFF)
if(TEHPC_EXAMPLES)
add_subdirectory(examples)
endif()
This diff is collapsed.
add_simulation(basic_example)
configure_file(basic_example_1.inp basic_example_1.inp COPYONLY)
configure_file(basic_example_1.run basic_example_1.run COPYONLY)
#include <iostream> // print to console (cout)
//#include <cmath> // math operations
//#include <fstream> // file streaming for input file
#include "Model.hpp"
int main(int argc, char *argv[])
{
std::string fname;
std::string odir;
if(argc<3) {
std::cerr << "Not enough arguments:"
<< " ./basic_example <sim-name> <output-directory>" << std::endl;
return 0;
}
else {
fname = argv[1];
odir = argv[2];
}
Model my_model;
std::cout << "read input" << std::endl;
my_model.read_input(fname+".inp");
my_model.assemble(); // Get unconstrained stiffness matrix + load vector
my_model.apply_BC();
std::cout << "solve" << std::endl;
my_model.solve(); // Solve linear system and get displacement vector
my_model.compute_reaction(); // Calculate reaction
my_model.print_results(odir,fname); // print results
return 0;
}
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)
odir=~/dev/output
name=basic_example_1
./basic_example $name $odir 2>&1 | tee $name.progress
cp $name.inp $odir/$name.inp
cp $name.run $odir/$name.run
cp $name.progress $odir/$name.progress
set(TEHPC_HEADER
Model.hpp
#Mesh.hpp
)
set(TEHPC_SRC
Model.cpp
#Mesh.cpp
)
set(TEHPC_DEPEND_FILES ${TEHPC_SRC} ${TEHPC_HEADER} PARENT_SCOPE)
add_library(tehpc ${TEHPC_SRC})
# if linking to other libraries, e.g. MPI
target_link_libraries(tehpc ${ARMADILLO_LIBRARIES})
#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 std::string & 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(const std::string & path,
const std::string & fname)
{
std::string path_to_file = path + "/" + fname + ".out";
std::ofstream ofile(path_to_file, std::ios::out);
ofile << "K_f:\n" << K_f << "\n";
ofile << "K:\n" << K << "\n";
ofile << "d:\n" << d << "\n";
ofile << "r:\n" << r << "\n";
ofile.close();
}
#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 std::string&);
void assemble ();
void apply_BC ();
void solve ();
void compute_reaction ();
void print_results(const std::string & path,
const std::string & fname);
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
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