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

week 10 group 1

parent 4634d308
build
*~
\ No newline at end of file
# Assignment
This week we look into how we can use existing code to make our code faster/better/simpler/etc. We actually have been doing this all semester long for the direct-stiffness-method example, where we used the external library armadillo. Today, we will have a more detailed look into how to include an external library. Please follow these steps:
1. refresh your memory about the code by looking into the `Model.hpp` and `Model.cpp` of week 9 (previous week, not current week). Think about how armadillo is used and if you see potential for improvement. Make a slide about it.
2. now consider the code provided in week 10 (current week). Specifically look into `Model.hpp`, `Model.cpp`, and `MatrixArma.hpp` and `VectorArma.hpp`. What is different? Study the `MatrixArma.hpp` and `VectorArma.hpp` files and each of the functions implemented for these classes. What could the purpose be for this change? What are the advantages and disadvantages? Summarize new approach on slides and provide thoughts about it.
3. in the `Model` class: where is armadillo still occuring? What needs to be done to remove it from the `Model` class? summarize on slides.
4. an additional file is provided as a starting point for solving this issue. Have a look and modify it to start removing the last bits of armadillo in the `Model` class. Note that it needs to be added to the `CMakeList.txt` file to be compiled. Explain your approach on slides
5. For your discussion: Given the objective of having armadillo included in a limited part of the code, was the current approach successful? Do you see additional potential for improvement? What are they? Summarize on slides.
6. Based on what you have learned today and over the past weeks, what classes would you implement for this code? What would be the members of these class? Summarize on slides.
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()
# ------------------------------------------------------
# test
option(TEHPC_TEST "tests" OFF)
if(${TEHPC_TEST})
enable_testing()
add_subdirectory(tests)
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
VectorArma.hpp
MatrixArma.hpp
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})
#ifndef MATRIXARMAHEADER
#define MATRIXARMAHEADER
#include <armadillo>
#include <iostream>
#include "VectorArma.hpp"
template <typename T>
class MatrixArma {
public:
// constructor
MatrixArma() {}
MatrixArma(int i, int j) : mat(i,j) {}
MatrixArma(const char * c) : mat(c) {}
MatrixArma(const std::string & string) : mat(string) {}
MatrixArma(const std::initializer_list< std::initializer_list<T> > & ilist) : mat(ilist) {}
MatrixArma(const arma::Mat<T> & m) : mat(m) {} // for operators
// destructor
~MatrixArma() {}
// FUNCTIONS
public:
void zeros() { this->mat.zeros(); }
MatrixArma<T> t() { return MatrixArma<T>(this->mat.t()); }
void swap(MatrixArma<T> & other) { this->mat.swap(other.mat); }
arma::Row<T> as_row() { return this->mat.as_row(); }
arma::Col<T> as_col() { return this->mat.as_col(); }
// OPERATORS
public:
T& operator()(int row, int col) { return this->mat(row,col); }
T operator()(int row, int col) const { return this->mat(row,col); }
friend std::ostream & operator << (std::ostream &out, const MatrixArma<T> &m) { out << m.mat; return out; }
MatrixArma<T> operator*(const MatrixArma<T> & m1) { return MatrixArma<T>(this->mat * m1.mat); }
MatrixArma<T> operator+(const MatrixArma<T> & m1) { return MatrixArma<T>(this->mat + m1.mat); }
VectorArma<T> operator*(const VectorArma<T> & v1) { return VectorArma<T>(this->mat * v1.storage()); }
arma::vec operator*(const arma::vec & other) { return (this->mat * other); }
// ACCESSORS
public:
arma::Mat<T> & storage() { return this->mat; }
// MEMBERS
protected:
arma::Mat<T> mat;
};
typedef MatrixArma<double> MatrixArmaD;
typedef MatrixArma<int> MatrixArmaI;
#endif
#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<MatrixArmaD>(&nodes, fin);
read_line<MatrixArmaI>(&elements, fin);
read_line<VectorArmaD>(&E, fin);
read_line<VectorArmaD>(&A, fin);
read_line<MatrixArmaI>(&pli, fin);
read_line<MatrixArmaD>(&pl, fin);
read_line<MatrixArmaI>(&pdi, fin);
read_line<MatrixArmaD>(&pd, fin);
}
void Model::local_stiffness_element(MatrixArmaD & k_e,
int el) {
double dx = nodes(elements(el,1)-1,0) - nodes(elements(el,0)-1,0);
double dy = nodes(elements(el,1)-1,1) - nodes(elements(el,0)-1,1);
double L = sqrt(dx*dx + dy*dy);
double cx = dx / L; // Sin of element angle with respect to global coordinates
double cy = dy / L; // Cosin of element angle with respect to global coordinates
double k = (A(el)*E(el))/L;
// local K
MatrixArmaD k_e_local( {{ k, 0, -k, 0},
{ 0, 0, 0, 0},
{-k, 0, k, 0},
{ 0, 0, 0, 0}});
// Rotation matrix from local angle to global
MatrixArmaD R ({{ cx, cy, 0, 0},
{-cy, cx, 0, 0},
{ 0, 0, cx, cy},
{ 0, 0, -cy, cx}});
MatrixArmaD 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
}
//---------------------------------------------------------------------------------------
// ASSEMBLY FUNCTION
// Assemble stiffness matrix from input struct containing all parsed input values from input file.
// -----------------------------------------------------------------------------------
void Model::assemble() {
MatrixArmaI 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 / )
MatrixArmaD* p_k_f = & K_f;
MatrixArmaD K_f(Nnp*n,Nnp*n);
K_f.zeros();
// Assemble floating K
for (int i=0;i<Nel;++i) {
MatrixArmaD k_e;
this->local_stiffness_element(k_e, i);
MatrixArmaD K_e(Nnp*n,Nnp*n);
K_e.zeros();
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 (){
MatrixArmaD* p_k = & K;
VectorArmaD* p_f = & f;
MatrixArmaD* p_pd = & pd;
MatrixArmaD K (K_f); // Start from floating structure
VectorArmaD f (pl.as_row().as_col()); // Initialize load
VectorArmaU pdl = VectorArmaU(find(pdi.as_row())); // contains 1 for prescribed, 0 for free (displacement dofs)
VectorArmaD pd ((*p_pd).as_row().as_col()); // prescribed displacements.
for (int i=0;(unsigned)i< pdl.size();++i){
f(pdl(i)) = pd(pdl(i));
for (int j=0; j<Nnp*n; ++j)
K(pdl(i),j) = 0.;
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 (){
VectorArmaD* p_r = & r;
VectorArmaD r = K_f * d;
VectorArmaU pdl = VectorArmaU(find(pli.as_row()));
for (int i=0;(unsigned)i<pdl.size();++i){
r(pdl(i)) = 0;
}
(*p_r).swap(r);
}
void Model::solve()
{
*(& d) = VectorArmaD(arma::solve (K.storage(),f.storage())); // 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 "MatrixArma.hpp"
#include "VectorArma.hpp"
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);
protected:
void local_stiffness_element(MatrixArmaD & k_e, int el);
private:
int Nel,Nen,Nnp,n; // Number of elements, nodes-per-element, nodes, dofs per node
MatrixArmaD nodes; // Coordinates of all nodes
MatrixArmaI elements; // Node index pair from every element
VectorArmaD E, A; // Young modulus and cross section of every element
MatrixArmaI pli, pdi; // constrained dofs matrix (1 -> constrained / 0 -> free)
MatrixArmaD pl, pd; // load and displacement boundary conditions for every dof
// Struct to store all vectors and matrices of the 2d truss stiffness model.
MatrixArmaD K,K_f; // Final stiffness matrix + unconstrained stiffness matrix (floating).
VectorArmaD f,d,r; // Load + displacement + reaction.
};
#endif
#ifndef SOLVERARMAHEADER
#define SOLVERARMAHEADER
#include <armadillo>
class SolverArma {
public:
// constructor
SolverArma(Model & model) : model(model) {}
// FUNCTIONS
public:
void solve();
// MEMBERS
protected:
Model & model; // reference to model
};
#endif
#ifndef VECTORARMAHEADER
#define VECTORARMAHEADER
#include <armadillo>
template <typename T>
class VectorArma {
public:
// constructor
VectorArma() {}
VectorArma(const arma::Col<T> & col) : vec(col) {}
//VectorArma(const arma::Row<T> & row) : vec(row) {}
// FUNCTIONS
public:
unsigned int size() { return this->vec.n_elem; }
void swap(VectorArma<T> & other) { this->vec.swap(other.vec); }
// OPERATORS
public:
T& operator()(int i) { return this->vec(i); }
T operator()(int i) const { return this->vec(i); }
friend std::ostream & operator << (std::ostream &out, const VectorArma<T> &v) { out << v.vec; return out; }
//VectorArma<T> & operator=(const arma::Col<T> & r) { return VectorArma<T>(r); }
//VectorArma<T> & operator=(const arma::uvec & r) { return VectorArma<T>(r); }
// ACCESSORS
public:
arma::Col<T> & storage() { return this->vec; }
const arma::Col<T> & storage() const { return this->vec; }
// MEMBERS
protected:
arma::Col<T> vec;
};
typedef VectorArma<double> VectorArmaD;
typedef VectorArma<arma::uword> VectorArmaU;
#endif
# define a function to implement test
function(tehpc_add_test TEST_NAME)
add_simulation(${TEST_NAME})
add_test(${TEST_NAME} ${TEST_NAME})
endfunction()
# dummy test
#tehpc_add_test(dummy.test)
#include <iostream>
int main()
{
std::cout << "start dummy test" << std::endl;
/*
if (some-condition) {
std::cerr << "test failed" << std::endl;
return 1; // failure
}
*/
// all tests passed: success
std::cout << "end dummy test" << std::endl;
return 0; // success
}
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