diff --git a/include/tug/Boundary.hpp b/include/tug/Boundary.hpp
index 2881d6a..10b9460 100644
--- a/include/tug/Boundary.hpp
+++ b/include/tug/Boundary.hpp
@@ -10,9 +10,6 @@
 #include "Grid.hpp"
 #include <cstddef>
 
-using namespace std;
-using namespace Eigen;
-
 /**
  * @brief Enum defining the two implemented boundary conditions.
  *
@@ -156,7 +153,7 @@ public:
    * @return vector<BoundaryElement> Contains the boundary conditions as
    * BoundaryElement objects.
    */
-  const vector<BoundaryElement> getBoundarySide(BC_SIDE side);
+  const std::vector<BoundaryElement> getBoundarySide(BC_SIDE side);
 
   /**
    * @brief Get thes Boundary Side Values as a vector. Value is -1 in case some
@@ -165,7 +162,7 @@ public:
    * @param side Boundary side for which the values are to be returned.
    * @return VectorXd Vector with values as doubles.
    */
-  VectorXd getBoundarySideValues(BC_SIDE side);
+  Eigen::VectorXd getBoundarySideValues(BC_SIDE side);
 
   /**
    * @brief Returns the boundary condition of a specified element on a given
@@ -207,7 +204,7 @@ public:
 private:
   Grid grid; // Boundary is directly dependent on the dimensions of a predefined
 
-  vector<vector<BoundaryElement>>
+  std::vector<std::vector<BoundaryElement>>
       boundaries; // Vector with Boundary Element information
 };
 
diff --git a/include/tug/Grid.hpp b/include/tug/Grid.hpp
index d40e281..1319677 100644
--- a/include/tug/Grid.hpp
+++ b/include/tug/Grid.hpp
@@ -11,8 +11,6 @@
 #include <Eigen/Core>
 #include <Eigen/Sparse>
 
-using namespace Eigen;
-
 class Grid {
 public:
   /**
@@ -45,7 +43,7 @@ public:
    * must have correct dimensions as defined in row and col. (Or length, in 1D
    * case).
    */
-  void setConcentrations(MatrixXd concentrations);
+  void setConcentrations(Eigen::MatrixXd concentrations);
 
   /**
    * @brief Gets the concentrations matrix for a Grid.
@@ -53,7 +51,7 @@ public:
    * @return MatrixXd An Eigen3 matrix holding the concentrations and having the
    *                  same dimensions as the grid.
    */
-  const MatrixXd getConcentrations();
+  const Eigen::MatrixXd getConcentrations();
 
   /**
    * @brief Set the alpha coefficients of a 1D-Grid. Grid must be one
@@ -62,7 +60,7 @@ public:
    * @param alpha An Eigen3 MatrixXd with 1 row holding the alpha coefficients.
    * Matrix columns must have same size as length of grid.
    */
-  void setAlpha(MatrixXd alpha);
+  void setAlpha(Eigen::MatrixXd alpha);
 
   /**
    * @brief Set the alpha coefficients of a 2D-Grid. Grid must be two
@@ -73,7 +71,7 @@ public:
    * @param alphaY An Eigen3 MatrixXd holding the alpha coefficients in
    * y-direction. Matrix must be of same size as the grid.
    */
-  void setAlpha(MatrixXd alphaX, MatrixXd alphaY);
+  void setAlpha(Eigen::MatrixXd alphaX, Eigen::MatrixXd alphaY);
 
   /**
    * @brief Gets the matrix of alpha coefficients of a 1D-Grid. Grid must be one
@@ -81,7 +79,7 @@ public:
    *
    * @return MatrixXd A matrix with 1 row holding the alpha coefficients.
    */
-  const MatrixXd getAlpha();
+  const Eigen::MatrixXd getAlpha();
 
   /**
    * @brief Gets the matrix of alpha coefficients in x-direction of a 2D-Grid.
@@ -89,7 +87,7 @@ public:
    *
    * @return MatrixXd A matrix holding the alpha coefficients in x-direction.
    */
-  const MatrixXd getAlphaX();
+  const Eigen::MatrixXd getAlphaX();
 
   /**
    * @brief Gets the matrix of alpha coefficients in y-direction of a 2D-Grid.
@@ -97,7 +95,7 @@ public:
    *
    * @return MatrixXd A matrix holding the alpha coefficients in y-direction.
    */
-  const MatrixXd getAlphaY();
+  const Eigen::MatrixXd getAlphaY();
 
   /**
    * @brief Gets the dimensions of the grid.
@@ -166,16 +164,16 @@ public:
   double getDeltaRow();
 
 private:
-  int col;                 // number of grid columns
-  int row;                 // number of grid rows
-  int dim;                 // 1D or 2D
-  double domainCol;        // number of domain columns
-  double domainRow;        // number of domain rows
-  double deltaCol;         // delta in x-direction (between columns)
-  double deltaRow;         // delta in y-direction (between rows)
-  MatrixXd concentrations; // Matrix holding grid concentrations
-  MatrixXd alphaX;         // Matrix holding alpha coefficients in x-direction
-  MatrixXd alphaY;         // Matrix holding alpha coefficients in y-direction
+  int col;                        // number of grid columns
+  int row;                        // number of grid rows
+  int dim;                        // 1D or 2D
+  double domainCol;               // number of domain columns
+  double domainRow;               // number of domain rows
+  double deltaCol;                // delta in x-direction (between columns)
+  double deltaRow;                // delta in y-direction (between rows)
+  Eigen::MatrixXd concentrations; // Matrix holding grid concentrations
+  Eigen::MatrixXd alphaX; // Matrix holding alpha coefficients in x-direction
+  Eigen::MatrixXd alphaY; // Matrix holding alpha coefficients in y-direction
 };
 
 #endif // GRID_H_
diff --git a/include/tug/Simulation.hpp b/include/tug/Simulation.hpp
index 60e87ac..571ef79 100644
--- a/include/tug/Simulation.hpp
+++ b/include/tug/Simulation.hpp
@@ -12,8 +12,6 @@
 #include "Boundary.hpp"
 #include "Grid.hpp"
 
-using namespace std;
-
 /**
  * @brief Enum defining the two implemented solution approaches.
  *
@@ -193,7 +191,7 @@ public:
    *
    * @return string Filename with configured simulation parameters.
    */
-  string createCSVfile();
+  std::string createCSVfile();
 
   /**
    * @brief Writes the currently calculated concentration values of the grid
@@ -202,7 +200,7 @@ public:
    * @param filename Name of the file to which the concentration values are
    *                 to be written.
    */
-  void printConcentrationsCSV(string filename);
+  void printConcentrationsCSV(std::string filename);
 
   /**
    * @brief Method starts the simulation process with the previously set
diff --git a/src/BTCS.cpp b/src/BTCS.cpp
index 0e0e822..3d162ed 100644
--- a/src/BTCS.cpp
+++ b/src/BTCS.cpp
@@ -14,13 +14,9 @@
 #include <tug/Boundary.hpp>
 #include <tug/Grid.hpp>
 
-#define NUM_THREADS_BTCS 10
-
-using namespace Eigen;
-
 // calculates coefficient for left boundary in constant case
-static tuple<double, double>
-calcLeftBoundaryCoeffConstant(MatrixXd &alpha, int rowIndex, double sx) {
+static std::tuple<double, double>
+calcLeftBoundaryCoeffConstant(Eigen::MatrixXd &alpha, int rowIndex, double sx) {
   double centerCoeff;
   double rightCoeff;
 
@@ -33,8 +29,8 @@ calcLeftBoundaryCoeffConstant(MatrixXd &alpha, int rowIndex, double sx) {
 }
 
 // calculates coefficient for left boundary in closed case
-static tuple<double, double>
-calcLeftBoundaryCoeffClosed(MatrixXd &alpha, int rowIndex, double sx) {
+static std::tuple<double, double>
+calcLeftBoundaryCoeffClosed(Eigen::MatrixXd &alpha, int rowIndex, double sx) {
   double centerCoeff;
   double rightCoeff;
 
@@ -46,9 +42,9 @@ calcLeftBoundaryCoeffClosed(MatrixXd &alpha, int rowIndex, double sx) {
 }
 
 // calculates coefficient for right boundary in constant case
-static tuple<double, double> calcRightBoundaryCoeffConstant(MatrixXd &alpha,
-                                                            int rowIndex, int n,
-                                                            double sx) {
+static std::tuple<double, double>
+calcRightBoundaryCoeffConstant(Eigen::MatrixXd &alpha, int rowIndex, int n,
+                               double sx) {
   double leftCoeff;
   double centerCoeff;
 
@@ -62,8 +58,9 @@ static tuple<double, double> calcRightBoundaryCoeffConstant(MatrixXd &alpha,
 }
 
 // calculates coefficient for right boundary in closed case
-static tuple<double, double>
-calcRightBoundaryCoeffClosed(MatrixXd &alpha, int rowIndex, int n, double sx) {
+static std::tuple<double, double>
+calcRightBoundaryCoeffClosed(Eigen::MatrixXd &alpha, int rowIndex, int n,
+                             double sx) {
   double leftCoeff;
   double centerCoeff;
 
@@ -76,15 +73,14 @@ calcRightBoundaryCoeffClosed(MatrixXd &alpha, int rowIndex, int n, double sx) {
 }
 
 // creates coefficient matrix for next time step from alphas in x-direction
-static SparseMatrix<double> createCoeffMatrix(MatrixXd &alpha,
-                                              vector<BoundaryElement> &bcLeft,
-                                              vector<BoundaryElement> &bcRight,
-                                              int numCols, int rowIndex,
-                                              double sx) {
+static Eigen::SparseMatrix<double>
+createCoeffMatrix(Eigen::MatrixXd &alpha, std::vector<BoundaryElement> &bcLeft,
+                  std::vector<BoundaryElement> &bcRight, int numCols,
+                  int rowIndex, double sx) {
 
   // square matrix of column^2 dimension for the coefficients
-  SparseMatrix<double> cm(numCols, numCols);
-  cm.reserve(VectorXi::Constant(numCols, 3));
+  Eigen::SparseMatrix<double> cm(numCols, numCols);
+  cm.reserve(Eigen::VectorXi::Constant(numCols, 3));
 
   // left column
   BC_TYPE type = bcLeft[rowIndex].getType();
@@ -140,8 +136,8 @@ static SparseMatrix<double> createCoeffMatrix(MatrixXd &alpha,
 
 // calculates explicity concentration at top boundary in constant case
 static double calcExplicitConcentrationsTopBoundaryConstant(
-    MatrixXd &concentrations, MatrixXd &alpha, vector<BoundaryElement> &bcTop,
-    int rowIndex, int i, double sy) {
+    Eigen::MatrixXd &concentrations, Eigen::MatrixXd &alpha,
+    std::vector<BoundaryElement> &bcTop, int rowIndex, int i, double sy) {
   double c;
 
   c = sy * calcAlphaIntercell(alpha(rowIndex, i), alpha(rowIndex + 1, i)) *
@@ -156,8 +152,10 @@ static double calcExplicitConcentrationsTopBoundaryConstant(
 }
 
 // calculates explicit concentration at top boundary in closed case
-static double calcExplicitConcentrationsTopBoundaryClosed(
-    MatrixXd &concentrations, MatrixXd &alpha, int rowIndex, int i, double sy) {
+static double
+calcExplicitConcentrationsTopBoundaryClosed(Eigen::MatrixXd &concentrations,
+                                            Eigen::MatrixXd &alpha,
+                                            int rowIndex, int i, double sy) {
   double c;
 
   c = sy * calcAlphaIntercell(alpha(rowIndex, i), alpha(rowIndex + 1, i)) *
@@ -171,8 +169,8 @@ static double calcExplicitConcentrationsTopBoundaryClosed(
 
 // calculates explicit concentration at bottom boundary in constant case
 static double calcExplicitConcentrationsBottomBoundaryConstant(
-    MatrixXd &concentrations, MatrixXd &alpha,
-    vector<BoundaryElement> &bcBottom, int rowIndex, int i, double sy) {
+    Eigen::MatrixXd &concentrations, Eigen::MatrixXd &alpha,
+    std::vector<BoundaryElement> &bcBottom, int rowIndex, int i, double sy) {
   double c;
 
   c = sy * alpha(rowIndex, i) * bcBottom[i].getValue() +
@@ -187,8 +185,10 @@ static double calcExplicitConcentrationsBottomBoundaryConstant(
 }
 
 // calculates explicit concentration at bottom boundary in closed case
-static double calcExplicitConcentrationsBottomBoundaryClosed(
-    MatrixXd &concentrations, MatrixXd &alpha, int rowIndex, int i, double sy) {
+static double
+calcExplicitConcentrationsBottomBoundaryClosed(Eigen::MatrixXd &concentrations,
+                                               Eigen::MatrixXd &alpha,
+                                               int rowIndex, int i, double sy) {
   double c;
 
   c = (1 -
@@ -202,13 +202,14 @@ static double calcExplicitConcentrationsBottomBoundaryClosed(
 
 // creates a solution vector for next time step from the current state of
 // concentrations
-static VectorXd createSolutionVector(
-    MatrixXd &concentrations, MatrixXd &alphaX, MatrixXd &alphaY,
-    vector<BoundaryElement> &bcLeft, vector<BoundaryElement> &bcRight,
-    vector<BoundaryElement> &bcTop, vector<BoundaryElement> &bcBottom,
-    int length, int rowIndex, double sx, double sy) {
+static Eigen::VectorXd createSolutionVector(
+    Eigen::MatrixXd &concentrations, Eigen::MatrixXd &alphaX,
+    Eigen::MatrixXd &alphaY, std::vector<BoundaryElement> &bcLeft,
+    std::vector<BoundaryElement> &bcRight, std::vector<BoundaryElement> &bcTop,
+    std::vector<BoundaryElement> &bcBottom, int length, int rowIndex, double sx,
+    double sy) {
 
-  VectorXd sv(length);
+  Eigen::VectorXd sv(length);
   int numRows = concentrations.rows();
   BC_TYPE type;
 
@@ -283,9 +284,10 @@ static VectorXd createSolutionVector(
 // solver for linear equation system; A corresponds to coefficient matrix,
 // b to the solution vector
 // use of EigenLU solver
-static VectorXd EigenLUAlgorithm(SparseMatrix<double> &A, VectorXd &b) {
+static Eigen::VectorXd EigenLUAlgorithm(Eigen::SparseMatrix<double> &A,
+                                        Eigen::VectorXd &b) {
 
-  SparseLU<SparseMatrix<double>> solver;
+  Eigen::SparseLU<Eigen::SparseMatrix<double>> solver;
   solver.analyzePattern(A);
   solver.factorize(A);
 
@@ -295,7 +297,8 @@ static VectorXd EigenLUAlgorithm(SparseMatrix<double> &A, VectorXd &b) {
 // solver for linear equation system; A corresponds to coefficient matrix,
 // b to the solution vector
 // implementation of Thomas Algorithm
-static VectorXd ThomasAlgorithm(SparseMatrix<double> &A, VectorXd &b) {
+static Eigen::VectorXd ThomasAlgorithm(Eigen::SparseMatrix<double> &A,
+                                       Eigen::VectorXd &b) {
   uint32_t n = b.size();
 
   Eigen::VectorXd a_diag(n);
@@ -337,22 +340,23 @@ static VectorXd ThomasAlgorithm(SparseMatrix<double> &A, VectorXd &b) {
 }
 
 // BTCS solution for 1D grid
-static void BTCS_1D(Grid &grid, Boundary &bc, double timestep,
-                    VectorXd (*solverFunc)(SparseMatrix<double> &A,
-                                           VectorXd &b)) {
+static void
+BTCS_1D(Grid &grid, Boundary &bc, double timestep,
+        Eigen::VectorXd (*solverFunc)(Eigen::SparseMatrix<double> &A,
+                                      Eigen::VectorXd &b)) {
   int length = grid.getLength();
   double sx = timestep / (grid.getDelta() * grid.getDelta());
 
-  VectorXd concentrations_t1(length);
+  Eigen::VectorXd concentrations_t1(length);
 
-  SparseMatrix<double> A;
-  VectorXd b(length);
+  Eigen::SparseMatrix<double> A;
+  Eigen::VectorXd b(length);
 
-  MatrixXd alpha = grid.getAlpha();
-  vector<BoundaryElement> bcLeft = bc.getBoundarySide(BC_SIDE_LEFT);
-  vector<BoundaryElement> bcRight = bc.getBoundarySide(BC_SIDE_RIGHT);
+  Eigen::MatrixXd alpha = grid.getAlpha();
+  std::vector<BoundaryElement> bcLeft = bc.getBoundarySide(BC_SIDE_LEFT);
+  std::vector<BoundaryElement> bcRight = bc.getBoundarySide(BC_SIDE_RIGHT);
 
-  MatrixXd concentrations = grid.getConcentrations();
+  Eigen::MatrixXd concentrations = grid.getConcentrations();
   int rowIndex = 0;
   A = createCoeffMatrix(alpha, bcLeft, bcRight, length, rowIndex,
                         sx); // this is exactly same as in 2D
@@ -376,29 +380,31 @@ static void BTCS_1D(Grid &grid, Boundary &bc, double timestep,
 }
 
 // BTCS solution for 2D grid
-static void BTCS_2D(Grid &grid, Boundary &bc, double timestep,
-                    VectorXd (*solverFunc)(SparseMatrix<double> &A,
-                                           VectorXd &b),
-                    int numThreads) {
+static void
+BTCS_2D(Grid &grid, Boundary &bc, double timestep,
+        Eigen::VectorXd (*solverFunc)(Eigen::SparseMatrix<double> &A,
+                                      Eigen::VectorXd &b),
+        int numThreads) {
   int rowMax = grid.getRow();
   int colMax = grid.getCol();
   double sx = timestep / (2 * grid.getDeltaCol() * grid.getDeltaCol());
   double sy = timestep / (2 * grid.getDeltaRow() * grid.getDeltaRow());
 
-  MatrixXd concentrations_t1 = MatrixXd::Constant(rowMax, colMax, 0);
-  VectorXd row_t1(colMax);
+  Eigen::MatrixXd concentrations_t1 =
+      Eigen::MatrixXd::Constant(rowMax, colMax, 0);
+  Eigen::VectorXd row_t1(colMax);
 
-  SparseMatrix<double> A;
-  VectorXd b;
+  Eigen::SparseMatrix<double> A;
+  Eigen::VectorXd b;
 
-  MatrixXd alphaX = grid.getAlphaX();
-  MatrixXd alphaY = grid.getAlphaY();
-  vector<BoundaryElement> bcLeft = bc.getBoundarySide(BC_SIDE_LEFT);
-  vector<BoundaryElement> bcRight = bc.getBoundarySide(BC_SIDE_RIGHT);
-  vector<BoundaryElement> bcTop = bc.getBoundarySide(BC_SIDE_TOP);
-  vector<BoundaryElement> bcBottom = bc.getBoundarySide(BC_SIDE_BOTTOM);
+  Eigen::MatrixXd alphaX = grid.getAlphaX();
+  Eigen::MatrixXd alphaY = grid.getAlphaY();
+  std::vector<BoundaryElement> bcLeft = bc.getBoundarySide(BC_SIDE_LEFT);
+  std::vector<BoundaryElement> bcRight = bc.getBoundarySide(BC_SIDE_RIGHT);
+  std::vector<BoundaryElement> bcTop = bc.getBoundarySide(BC_SIDE_TOP);
+  std::vector<BoundaryElement> bcBottom = bc.getBoundarySide(BC_SIDE_BOTTOM);
 
-  MatrixXd concentrations = grid.getConcentrations();
+  Eigen::MatrixXd concentrations = grid.getConcentrations();
 
 #pragma omp parallel for num_threads(numThreads) private(A, b, row_t1)
   for (int i = 0; i < rowMax; i++) {
@@ -407,7 +413,7 @@ static void BTCS_2D(Grid &grid, Boundary &bc, double timestep,
     b = createSolutionVector(concentrations, alphaX, alphaY, bcLeft, bcRight,
                              bcTop, bcBottom, colMax, i, sx, sy);
 
-    SparseLU<SparseMatrix<double>> solver;
+    Eigen::SparseLU<Eigen::SparseMatrix<double>> solver;
 
     row_t1 = solverFunc(A, b);
 
diff --git a/src/Boundary.cpp b/src/Boundary.cpp
index ea7c891..083e802 100644
--- a/src/Boundary.cpp
+++ b/src/Boundary.cpp
@@ -5,8 +5,6 @@
 #include <stdexcept>
 #include <tug/Boundary.hpp>
 
-using namespace std;
-
 BoundaryElement::BoundaryElement() {
 
   this->type = BC_TYPE_CLOSED;
@@ -38,22 +36,22 @@ double BoundaryElement::getValue() { return this->value; }
 
 Boundary::Boundary(Grid grid) : grid(grid) {
   if (grid.getDim() == 1) {
-    this->boundaries = vector<vector<BoundaryElement>>(
+    this->boundaries = std::vector<std::vector<BoundaryElement>>(
         2); // in 1D only left and right boundary
 
     this->boundaries[BC_SIDE_LEFT].push_back(BoundaryElement());
     this->boundaries[BC_SIDE_RIGHT].push_back(BoundaryElement());
   } else if (grid.getDim() == 2) {
-    this->boundaries = vector<vector<BoundaryElement>>(4);
+    this->boundaries = std::vector<std::vector<BoundaryElement>>(4);
 
     this->boundaries[BC_SIDE_LEFT] =
-        vector<BoundaryElement>(grid.getRow(), BoundaryElement());
+        std::vector<BoundaryElement>(grid.getRow(), BoundaryElement());
     this->boundaries[BC_SIDE_RIGHT] =
-        vector<BoundaryElement>(grid.getRow(), BoundaryElement());
+        std::vector<BoundaryElement>(grid.getRow(), BoundaryElement());
     this->boundaries[BC_SIDE_TOP] =
-        vector<BoundaryElement>(grid.getCol(), BoundaryElement());
+        std::vector<BoundaryElement>(grid.getCol(), BoundaryElement());
     this->boundaries[BC_SIDE_BOTTOM] =
-        vector<BoundaryElement>(grid.getCol(), BoundaryElement());
+        std::vector<BoundaryElement>(grid.getCol(), BoundaryElement());
   }
 }
 
@@ -72,7 +70,7 @@ void Boundary::setBoundarySideClosed(BC_SIDE side) {
   } else {
     n = grid.getCol();
   }
-  this->boundaries[side] = vector<BoundaryElement>(n, BoundaryElement());
+  this->boundaries[side] = std::vector<BoundaryElement>(n, BoundaryElement());
 }
 
 void Boundary::setBoundarySideConstant(BC_SIDE side, double value) {
@@ -90,7 +88,7 @@ void Boundary::setBoundarySideConstant(BC_SIDE side, double value) {
   } else {
     n = grid.getCol();
   }
-  this->boundaries[side] = vector<BoundaryElement>(n, BoundaryElement(value));
+  this->boundaries[side] = std::vector<BoundaryElement>(n, BoundaryElement(value));
 }
 
 void Boundary::setBoundaryElementClosed(BC_SIDE side, int index) {
@@ -111,7 +109,7 @@ void Boundary::setBoundaryElementConstant(BC_SIDE side, int index,
   this->boundaries[side][index].setValue(value);
 }
 
-const vector<BoundaryElement> Boundary::getBoundarySide(BC_SIDE side) {
+const std::vector<BoundaryElement> Boundary::getBoundarySide(BC_SIDE side) {
   if (grid.getDim() == 1) {
     if ((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)) {
       throw_invalid_argument(
@@ -122,9 +120,9 @@ const vector<BoundaryElement> Boundary::getBoundarySide(BC_SIDE side) {
   return this->boundaries[side];
 }
 
-VectorXd Boundary::getBoundarySideValues(BC_SIDE side) {
+Eigen::VectorXd Boundary::getBoundarySideValues(BC_SIDE side) {
   int length = boundaries[side].size();
-  VectorXd values(length);
+  Eigen::VectorXd values(length);
 
   for (int i = 0; i < length; i++) {
     if (getBoundaryElementType(side, i) == BC_TYPE_CLOSED) {
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 15ce0e8..4bbb32f 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -6,4 +6,4 @@ if(TUG_USE_OPENMP AND OpenMP_CXX_FOUND)
   target_link_libraries(tug OpenMP::OpenMP_CXX)
 endif()
 
-target_include_directories(tug PUBLIC ../include)
+target_include_directories(tug PUBLIC ${PROJECT_SOURCE_DIR}/include)
diff --git a/src/FTCS.cpp b/src/FTCS.cpp
index f6b8780..01eaf8d 100644
--- a/src/FTCS.cpp
+++ b/src/FTCS.cpp
@@ -13,8 +13,6 @@
 #include <omp.h>
 #include <tug/Boundary.hpp>
 
-using namespace std;
-
 // calculates horizontal change on one cell independent of boundary type
 static double calcHorizontalChange(Grid &grid, int &row, int &col) {
 
@@ -222,7 +220,7 @@ static void FTCS_1D(Grid &grid, Boundary &bc, double &timestep) {
   double deltaCol = grid.getDeltaCol();
 
   // matrix for concentrations at time t+1
-  MatrixXd concentrations_t1 = MatrixXd::Constant(1, colMax, 0);
+  Eigen::MatrixXd concentrations_t1 = Eigen::MatrixXd::Constant(1, colMax, 0);
 
   // only one row in 1D case -> row constant at index 0
   int row = 0;
@@ -262,7 +260,7 @@ static void FTCS_2D(Grid &grid, Boundary &bc, double &timestep,
   double deltaCol = grid.getDeltaCol();
 
   // matrix for concentrations at time t+1
-  MatrixXd concentrations_t1 = MatrixXd::Constant(rowMax, colMax, 0);
+  Eigen::MatrixXd concentrations_t1 = Eigen::MatrixXd::Constant(rowMax, colMax, 0);
 
   // inner cells
   // these are independent of the boundary condition type
diff --git a/src/Grid.cpp b/src/Grid.cpp
index c19d935..dad79e9 100644
--- a/src/Grid.cpp
+++ b/src/Grid.cpp
@@ -15,8 +15,8 @@ Grid::Grid(int length) {
   this->deltaCol = double(this->domainCol) / double(this->col); // -> 1
   this->dim = 1;
 
-  this->concentrations = MatrixXd::Constant(1, col, 20);
-  this->alphaX = MatrixXd::Constant(1, col, 1);
+  this->concentrations = Eigen::MatrixXd::Constant(1, col, 20);
+  this->alphaX = Eigen::MatrixXd::Constant(1, col, 1);
 }
 
 Grid::Grid(int row, int col) {
@@ -33,12 +33,12 @@ Grid::Grid(int row, int col) {
   this->deltaCol = double(this->domainCol) / double(this->col); // -> 1
   this->dim = 2;
 
-  this->concentrations = MatrixXd::Constant(row, col, 20);
-  this->alphaX = MatrixXd::Constant(row, col, 1);
-  this->alphaY = MatrixXd::Constant(row, col, 1);
+  this->concentrations = Eigen::MatrixXd::Constant(row, col, 20);
+  this->alphaX = Eigen::MatrixXd::Constant(row, col, 1);
+  this->alphaY = Eigen::MatrixXd::Constant(row, col, 1);
 }
 
-void Grid::setConcentrations(MatrixXd concentrations) {
+void Grid::setConcentrations(Eigen::MatrixXd concentrations) {
   if (concentrations.rows() != this->row ||
       concentrations.cols() != this->col) {
     throw_invalid_argument(
@@ -48,9 +48,9 @@ void Grid::setConcentrations(MatrixXd concentrations) {
   this->concentrations = concentrations;
 }
 
-const MatrixXd Grid::getConcentrations() { return this->concentrations; }
+const Eigen::MatrixXd Grid::getConcentrations() { return this->concentrations; }
 
-void Grid::setAlpha(MatrixXd alpha) {
+void Grid::setAlpha(Eigen::MatrixXd alpha) {
   if (dim != 1) {
     throw_invalid_argument("Grid is not one dimensional, you should probably "
                            "use 2D setter function!");
@@ -63,7 +63,7 @@ void Grid::setAlpha(MatrixXd alpha) {
   this->alphaX = alpha;
 }
 
-void Grid::setAlpha(MatrixXd alphaX, MatrixXd alphaY) {
+void Grid::setAlpha(Eigen::MatrixXd alphaX, Eigen::MatrixXd alphaY) {
   if (dim != 2) {
     throw_invalid_argument("Grid is not two dimensional, you should probably "
                            "use 1D setter function!");
@@ -81,7 +81,7 @@ void Grid::setAlpha(MatrixXd alphaX, MatrixXd alphaY) {
   this->alphaY = alphaY;
 }
 
-const MatrixXd Grid::getAlpha() {
+const Eigen::MatrixXd Grid::getAlpha() {
   if (dim != 1) {
     throw_invalid_argument("Grid is not one dimensional, you should probably "
                            "use either getAlphaX() or getAlphaY()!");
@@ -90,7 +90,7 @@ const MatrixXd Grid::getAlpha() {
   return this->alphaX;
 }
 
-const MatrixXd Grid::getAlphaX() {
+const Eigen::MatrixXd Grid::getAlphaX() {
   if (dim != 2) {
     throw_invalid_argument(
         "Grid is not two dimensional, you should probably use getAlpha()!");
@@ -99,7 +99,7 @@ const MatrixXd Grid::getAlphaX() {
   return this->alphaX;
 }
 
-const MatrixXd Grid::getAlphaY() {
+const Eigen::MatrixXd Grid::getAlphaY() {
   if (dim != 2) {
     throw_invalid_argument(
         "Grid is not two dimensional, you should probably use getAlpha()!");
diff --git a/src/Simulation.cpp b/src/Simulation.cpp
index 0d7357b..4fa0eea 100644
--- a/src/Simulation.cpp
+++ b/src/Simulation.cpp
@@ -63,7 +63,7 @@ void Simulation::setTimestep(double timestep) {
     double deltaColSquare = grid.getDeltaCol() * grid.getDeltaCol();
     double minDeltaSquare;
     double maxAlphaX, maxAlphaY, maxAlpha;
-    string dim;
+    std::string dim;
     if (grid.getDim() == 2) {
       dim = "2D";
 
@@ -91,31 +91,31 @@ void Simulation::setTimestep(double timestep) {
     // not work in some cases double CFL_Wiki = 1 / (4 * maxAlpha *
     // ((1/deltaRowSquare) + (1/deltaColSquare)));
 
-    string approachPrefix =
+    std::string approachPrefix =
         (approach == 0) ? "FTCS" : ((approach == 1) ? "BTCS" : "CRNI");
-    cout << approachPrefix << "_" << dim << " :: CFL condition: " << cfl
-         << endl;
-    cout << approachPrefix << "_" << dim << " :: required dt=" << timestep
-         << endl;
+    std::cout << approachPrefix << "_" << dim << " :: CFL condition: " << cfl
+         << std::endl;
+    std::cout << approachPrefix << "_" << dim << " :: required dt=" << timestep
+         << std::endl;
 
     if (timestep > cfl) {
 
       this->innerIterations = (int)ceil(timestep / cfl);
       this->timestep = timestep / (double)innerIterations;
 
-      cerr << "Warning :: Timestep was adjusted, because of stability "
+      std::cerr << "Warning :: Timestep was adjusted, because of stability "
               "conditions. Time duration was approximately preserved by "
               "adjusting internal number of iterations."
-           << endl;
-      cout << approachPrefix << "_" << dim << " :: Required "
+           << std::endl;
+      std::cout << approachPrefix << "_" << dim << " :: Required "
            << this->innerIterations
-           << " inner iterations with dt=" << this->timestep << endl;
+           << " inner iterations with dt=" << this->timestep << std::endl;
 
     } else {
 
       this->timestep = timestep;
-      cout << approachPrefix << "_" << dim
-           << " :: No inner iterations required, dt=" << timestep << endl;
+      std::cout << approachPrefix << "_" << dim
+           << " :: No inner iterations required, dt=" << timestep << std::endl;
     }
 
   } else {
@@ -134,10 +134,10 @@ void Simulation::setIterations(int iterations) {
 
 void Simulation::setSolver(SOLVER solver) {
   if (this->approach == FTCS_APPROACH) {
-    cerr << "Warning: Solver was set, but FTCS approach initialized. Setting "
+    std::cerr << "Warning: Solver was set, but FTCS approach initialized. Setting "
             "the solver "
             "is thus without effect."
-         << endl;
+         << std::endl;
   }
 
   this->solver = solver;
@@ -148,9 +148,9 @@ void Simulation::setNumberThreads(int numThreads) {
     this->numThreads = numThreads;
   } else {
     int maxThreadNumber = omp_get_num_procs();
-    string outputMessage =
+    std::string outputMessage =
         "Number of threads exceeds the number of processor cores (" +
-        to_string(maxThreadNumber) + ") or is less than 1.";
+        std::to_string(maxThreadNumber) + ") or is less than 1.";
 
     throw_invalid_argument(outputMessage);
   }
@@ -159,28 +159,28 @@ void Simulation::setNumberThreads(int numThreads) {
 int Simulation::getIterations() { return this->iterations; }
 
 void Simulation::printConcentrationsConsole() {
-  cout << grid.getConcentrations() << endl;
-  cout << endl;
+  std::cout << grid.getConcentrations() << std::endl;
+  std::cout << std::endl;
 }
 
-string Simulation::createCSVfile() {
-  ofstream file;
+std::string Simulation::createCSVfile() {
+  std::ofstream file;
   int appendIdent = 0;
-  string appendIdentString;
+  std::string appendIdentString;
 
   // string approachString = (approach == 0) ? "FTCS" : "BTCS";
-  string approachString =
+  std::string approachString =
       (approach == 0) ? "FTCS" : ((approach == 1) ? "BTCS" : "CRNI");
-  string row = to_string(grid.getRow());
-  string col = to_string(grid.getCol());
-  string numIterations = to_string(iterations);
+  std::string row = std::to_string(grid.getRow());
+  std::string col = std::to_string(grid.getCol());
+  std::string numIterations = std::to_string(iterations);
 
-  string filename =
+  std::string filename =
       approachString + "_" + row + "_" + col + "_" + numIterations + ".csv";
 
-  while (filesystem::exists(filename)) {
+  while (std::filesystem::exists(filename)) {
     appendIdent += 1;
-    appendIdentString = to_string(appendIdent);
+    appendIdentString = std::to_string(appendIdent);
     filename = approachString + "_" + row + "_" + col + "_" + numIterations +
                "-" + appendIdentString + ".csv";
   }
@@ -193,16 +193,16 @@ string Simulation::createCSVfile() {
   // adds lines at the beginning of verbose output csv that represent the
   // boundary conditions and their values -1 in case of closed boundary
   if (csv_output == CSV_OUTPUT_XTREME) {
-    IOFormat one_row(StreamPrecision, DontAlignCols, "", " ");
+    Eigen::IOFormat one_row(Eigen::StreamPrecision, Eigen::DontAlignCols, "", " ");
     file << bc.getBoundarySideValues(BC_SIDE_LEFT).format(one_row)
-         << endl; // boundary left
+         << std::endl; // boundary left
     file << bc.getBoundarySideValues(BC_SIDE_RIGHT).format(one_row)
-         << endl; // boundary right
+         << std::endl; // boundary right
     file << bc.getBoundarySideValues(BC_SIDE_TOP).format(one_row)
-         << endl; // boundary top
+         << std::endl; // boundary top
     file << bc.getBoundarySideValues(BC_SIDE_BOTTOM).format(one_row)
-         << endl; // boundary bottom
-    file << endl << endl;
+         << std::endl; // boundary bottom
+    file << std::endl << std::endl;
   }
 
   file.close();
@@ -210,17 +210,17 @@ string Simulation::createCSVfile() {
   return filename;
 }
 
-void Simulation::printConcentrationsCSV(string filename) {
-  ofstream file;
+void Simulation::printConcentrationsCSV(std::string filename) {
+  std::ofstream file;
 
   file.open(filename, std::ios_base::app);
   if (!file) {
     exit(1);
   }
 
-  IOFormat do_not_align(StreamPrecision, DontAlignCols);
-  file << grid.getConcentrations().format(do_not_align) << endl;
-  file << endl << endl;
+  Eigen::IOFormat do_not_align(Eigen::StreamPrecision, Eigen::DontAlignCols);
+  file << grid.getConcentrations().format(do_not_align) << std::endl;
+  file << std::endl << std::endl;
   file.close();
 }
 
@@ -232,7 +232,7 @@ void Simulation::run() {
     throw_invalid_argument("Number of iterations are not set!");
   }
 
-  string filename;
+  std::string filename;
   if (this->console_output > CONSOLE_OUTPUT_OFF) {
     printConcentrationsConsole();
   }
@@ -295,9 +295,9 @@ void Simulation::run() {
     // TODO this implementation is very inefficient!
     // a separate implementation that sets up a specific tridiagonal matrix for
     // Crank-Nicolson would be better
-    MatrixXd concentrations;
-    MatrixXd concentrationsFTCS;
-    MatrixXd concentrationsResult;
+    Eigen::MatrixXd concentrations;
+    Eigen::MatrixXd concentrationsFTCS;
+    Eigen::MatrixXd concentrationsResult;
     for (int i = 0; i < iterations * innerIterations; i++) {
       if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) {
         printConcentrationsConsole();
@@ -328,10 +328,10 @@ void Simulation::run() {
     printConcentrationsCSV(filename);
   }
   if (this->time_measure > TIME_MEASURE_OFF) {
-    string approachString =
+    std::string approachString =
         (approach == 0) ? "FTCS" : ((approach == 1) ? "BTCS" : "CRNI");
-    string dimString = (grid.getDim() == 1) ? "-1D" : "-2D";
-    cout << approachString << dimString << ":: run() finished in "
-         << milliseconds.count() << "ms" << endl;
+    std::string dimString = (grid.getDim() == 1) ? "-1D" : "-2D";
+    std::cout << approachString << dimString << ":: run() finished in "
+         << milliseconds.count() << "ms" << std::endl;
   }
 }