diff --git a/projects/braidingTightBinding/lib/coupler.hpp b/projects/braidingTightBinding/lib/coupler.hpp
index ed35420a649e80033fd3202682ac4c36dc97156e..1a50b6d0c4ef973bc558d4e5b130b39cfa3b8576 100644
--- a/projects/braidingTightBinding/lib/coupler.hpp
+++ b/projects/braidingTightBinding/lib/coupler.hpp
@@ -4,14 +4,14 @@
template <typename value_t, typename drum_t>
class Coupler{
- public:
- //constructors
- Coupler() = default;
- ~Coupler() = default;
+public:
+ //constructors
+ Coupler() = default;
+ ~Coupler() = default;
- virtual void precompute(const value_t t_end, const value_t dt, const std::vector<drum_t>& drum_vec) noexcept = 0;
- virtual void step(value_t dt) noexcept = 0;
- virtual value_t operator()(const size_t drum_index, const size_t neighbour_index) const noexcept = 0;
+ virtual void precompute(const value_t t_end, const value_t dt, const std::vector<drum_t>& drum_vec) noexcept = 0;
+ virtual void step(value_t dt) noexcept = 0;
+ virtual value_t operator()(const size_t drum_index, const size_t neighbour_index) const noexcept = 0;
};
#endif
diff --git a/projects/braidingTightBinding/lib/diagonalizer.hpp b/projects/braidingTightBinding/lib/diagonalizer.hpp
index 7d0ab0c04630534f0789e6e635855e409621ef54..75e4364b7deda801c8d3ebfc0c4fd86365079f84 100644
--- a/projects/braidingTightBinding/lib/diagonalizer.hpp
+++ b/projects/braidingTightBinding/lib/diagonalizer.hpp
@@ -18,74 +18,78 @@ extern "C" void dsyev_(
//TODO: Rule of 7 conformity (nontrivial heap resource)
class Diagonalizer{
- public:
- Diagonalizer() = default;
- ~Diagonalizer() = default;
-
- std::vector<double> ev(const std::vector<double>& matrix, const size_t N){
- //prepare memory
- matrix_ = matrix;
- if(eigenvalues_.size() != N){
- eigenvalues_.clear();
- eigenvalues_.reserve(N - eigenvalues_.capacity());
- }
- for(size_t i = 0; i < N; ++i)
- eigenvalues_.push_back(0.);
-
- //prepare workspace
- N_ = N;
- prepare_workspace('N');
-
- dsyev_('N', 'L', N_, matrix_.data(), N_, eigenvalues_.data(), work_, lwork_, info_);
-
- if(info_ != 0)
- throw("Diagonalization failed!");
-
- //clean up
- delete[] work_;
-
- return eigenvalues_;
- }
-
- std::pair<std::vector<double>, std::vector<double> > evv(const std::vector<double>& matrix, const size_t N){
- //prepare memory
- matrix_ = matrix;
- if(eigenvalues_.size() != N){
- eigenvalues_.clear();
- eigenvalues_.reserve(N - eigenvalues_.capacity());
- }
- for(size_t i = 0; i < N; ++i)
- eigenvalues_.push_back(0.);
-
- //prepare workspace
- N_ = N;
- prepare_workspace('V');
-
- dsyev_('V', 'L', N_, matrix_.data(), N_, eigenvalues_.data(), work_, lwork_, info_);
-
- if(info_ != 0)
- throw("Diagonalization failed!");
-
- //clean up
- delete[] work_;
-
- return std::pair<std::vector<double>, std::vector<double> > (eigenvalues_, matrix_);
- }
-
- private:
- void prepare_workspace(char type){
- dsyev_(type, 'L', N_, matrix_.data(), N_, eigenvalues_.data(), &dwork_, -1, info_);
- lwork_ = static_cast<int>(dwork_);
- work_ = new double[lwork_];
- }
-
- int info_;
- double dwork_;
- int lwork_;
- int N_; //linear matrix dimension
- std::vector<double> matrix_; //the matrix
- std::vector<double> eigenvalues_; //eigenvalues
- double* work_;
+public:
+ Diagonalizer() = default;
+ ~Diagonalizer() = default;
+
+ std::vector<double> ev(const std::vector<double>& matrix, const size_t N){
+ //prepare memory
+ matrix_ = matrix;
+ if(eigenvalues_.size() != N){
+ eigenvalues_.clear();
+ eigenvalues_.reserve(N - eigenvalues_.capacity());
+ }
+ for(size_t i = 0; i < N; ++i){
+ eigenvalues_.push_back(0.);
+ }
+
+ //prepare workspace
+ N_ = N;
+ prepare_workspace('N');
+
+ dsyev_('N', 'L', N_, matrix_.data(), N_, eigenvalues_.data(), work_, lwork_, info_);
+
+ if(info_ != 0){
+ throw("Diagonalization failed!");
+ }
+
+ //clean up
+ delete[] work_;
+
+ return eigenvalues_;
+ }
+
+ std::pair<std::vector<double>, std::vector<double> > evv(const std::vector<double>& matrix, const size_t N){
+ //prepare memory
+ matrix_ = matrix;
+ if(eigenvalues_.size() != N){
+ eigenvalues_.clear();
+ eigenvalues_.reserve(N - eigenvalues_.capacity());
+ }
+ for(size_t i = 0; i < N; ++i){
+ eigenvalues_.push_back(0.);
+ }
+
+ //prepare workspace
+ N_ = N;
+ prepare_workspace('V');
+
+ dsyev_('V', 'L', N_, matrix_.data(), N_, eigenvalues_.data(), work_, lwork_, info_);
+
+ if(info_ != 0){
+ throw("Diagonalization failed!");
+ }
+
+ //clean up
+ delete[] work_;
+
+ return std::pair<std::vector<double>, std::vector<double> > (eigenvalues_, matrix_);
+ }
+
+private:
+ void prepare_workspace(char type){
+ dsyev_(type, 'L', N_, matrix_.data(), N_, eigenvalues_.data(), &dwork_, -1, info_);
+ lwork_ = static_cast<int>(dwork_);
+ work_ = new double[lwork_];
+ }
+
+ int info_;
+ double dwork_;
+ int lwork_;
+ int N_; //linear matrix dimension
+ std::vector<double> matrix_; //the matrix
+ std::vector<double> eigenvalues_; //eigenvalues
+ double* work_;
};
#endif
diff --git a/projects/braidingTightBinding/lib/driver.hpp b/projects/braidingTightBinding/lib/driver.hpp
index b1d4217d7f0257a3b189e187e2427c346690fa82..4b820b677727ad3bbb1d75d6ea5a4f1102e7457c 100644
--- a/projects/braidingTightBinding/lib/driver.hpp
+++ b/projects/braidingTightBinding/lib/driver.hpp
@@ -5,14 +5,14 @@
template <typename value_t, typename drum_t>
class Driver{
- public:
- //constructors
- Driver() = default;
- ~Driver() = default;
+public:
+ //constructors
+ Driver() = default;
+ ~Driver() = default;
- virtual void precompute(const value_t t_end, const value_t dt, const std::vector<drum_t>& drum_vec) noexcept = 0;
- virtual void step(value_t dt) noexcept = 0;
- virtual value_t operator()(const size_t drum_index) const noexcept = 0;
+ virtual void precompute(const value_t t_end, const value_t dt, const std::vector<drum_t>& drum_vec) noexcept = 0;
+ virtual void step(value_t dt) noexcept = 0;
+ virtual value_t operator()(const size_t drum_index) const noexcept = 0;
};
#endif
diff --git a/projects/braidingTightBinding/lib/drum.hpp b/projects/braidingTightBinding/lib/drum.hpp
index ee5f3f7252e624f2c9849bbfac659a284f458ca1..ffa629229e59457d794aaa1650af959221608191 100644
--- a/projects/braidingTightBinding/lib/drum.hpp
+++ b/projects/braidingTightBinding/lib/drum.hpp
@@ -10,67 +10,67 @@
template <typename value_t, typename params_t, typename vars_t, typename sbuffer_t>
class Drum{
- public:
- //'structors
- Drum(const params_t& dp) noexcept: params_(dp) {}
- Drum() = delete; //default construction not allowed
- Drum(const Drum&) = default;
- Drum& operator=(const Drum&) = default;
- ~Drum() = default;
+public:
+ //'structors
+ Drum(const params_t& dp) noexcept: params_(dp) {}
+ Drum() = delete; //default construction not allowed
+ Drum(const Drum&) = default;
+ Drum& operator=(const Drum&) = default;
+ ~Drum() = default;
- //access
- params_t& get_parameters() noexcept
- {
- return params_;
- }
+ //access
+ params_t& get_parameters() noexcept
+ {
+ return params_;
+ }
- const params_t& get_parameters() const noexcept
- {
- return params_;
- }
+ const params_t& get_parameters() const noexcept
+ {
+ return params_;
+ }
- vars_t& get_variables() noexcept
- {
- return vars_;
- }
+ vars_t& get_variables() noexcept
+ {
+ return vars_;
+ }
- const vars_t& get_variables() const noexcept
- {
- return vars_;
- }
+ const vars_t& get_variables() const noexcept
+ {
+ return vars_;
+ }
- sbuffer_t& get_sbuffer() noexcept
- {
- return sbuffer_;
- }
+ sbuffer_t& get_sbuffer() noexcept
+ {
+ return sbuffer_;
+ }
- const sbuffer_t& get_sbuffer() const noexcept
- {
- return sbuffer_;
- }
+ const sbuffer_t& get_sbuffer() const noexcept
+ {
+ return sbuffer_;
+ }
- //modifiers
- void set_coupling_0(value_t t0)
- {
- vars_.t0 = t0;
- }
- void set_coupling_1(value_t t1)
- {
- vars_.t1 = t1;
- }
- void set_coupling_2(value_t t2)
- {
- vars_.t2 = t2;
- }
- void set_drive(value_t V)
- {
- vars_.V = V;
- }
+ //modifiers
+ void set_coupling_0(value_t t0)
+ {
+ vars_.t0 = t0;
+ }
+ void set_coupling_1(value_t t1)
+ {
+ vars_.t1 = t1;
+ }
+ void set_coupling_2(value_t t2)
+ {
+ vars_.t2 = t2;
+ }
+ void set_drive(value_t V)
+ {
+ vars_.V = V;
+ }
- private:
- params_t params_; // constant drum parameters
- vars_t vars_; // drum variables
- sbuffer_t sbuffer_; //stepper buffer
+private:
+ params_t params_; // constant drum parameters
+ vars_t vars_; // drum variables
+ sbuffer_t sbuffer_; //stepper buffer
};
#endif
diff --git a/projects/braidingTightBinding/lib/force.hpp b/projects/braidingTightBinding/lib/force.hpp
index 094fffcfa4cbfb930090cd77860839ea80089668..eff0bcd459702feb8d5dd2bf696c889b0324e5f4 100644
--- a/projects/braidingTightBinding/lib/force.hpp
+++ b/projects/braidingTightBinding/lib/force.hpp
@@ -5,18 +5,18 @@
//Force evaluation interface
template <typename value_t, typename params_t, typename vars_t, typename buffer_t>
class Force{
- public:
- //constructors
- Force() = default;
- ~Force() = default;
+public:
+ //constructors
+ Force() = default;
+ ~Force() = default;
- //virtual functional
- virtual value_t operator()(
- const Drum<value_t, params_t, vars_t, buffer_t>& drum, //Drum we're calculating the force on
- const Drum<value_t, params_t, vars_t, buffer_t>& neighbour0, //Neighbour 0
- const Drum<value_t, params_t, vars_t, buffer_t>& neighbour1, //Neighbour 1
- const Drum<value_t, params_t, vars_t, buffer_t>& neighbour2, //Neighbour 2
- const value_t time) const noexcept = 0;
+ //virtual functional
+ virtual value_t operator()(
+ const Drum<value_t, params_t, vars_t, buffer_t>& drum, //Drum we're calculating the force on
+ const Drum<value_t, params_t, vars_t, buffer_t>& neighbour0, //Neighbour 0
+ const Drum<value_t, params_t, vars_t, buffer_t>& neighbour1, //Neighbour 1
+ const Drum<value_t, params_t, vars_t, buffer_t>& neighbour2, //Neighbour 2
+ const value_t time) const noexcept = 0;
};
#endif
diff --git a/projects/braidingTightBinding/lib/lattice_generator.hpp b/projects/braidingTightBinding/lib/lattice_generator.hpp
index 20819df1be678f9073606d752c15345dc8ec47be..d3b6253bdafc805cf513b91d189b26e53d0e7d2a 100644
--- a/projects/braidingTightBinding/lib/lattice_generator.hpp
+++ b/projects/braidingTightBinding/lib/lattice_generator.hpp
@@ -6,15 +6,15 @@
template <typename value_t, typename params_t, typename vars_t, typename sbuffer_t>
class LatticeGenerator{
- public:
- LatticeGenerator() = default;
- ~LatticeGenerator() = default;
+public:
+ LatticeGenerator() = default;
+ ~LatticeGenerator() = default;
- //returns a std::pair of a std::vector of drums in the system and an adjacency vector of size_t.
- //the last drum (.back()) is the drum seen when no neighbour is present.
- //TODO: later, overload with a vector of drum parameters to set each drum's properties individually
- virtual std::pair<std::vector<Drum<value_t, params_t, vars_t, sbuffer_t> >, std::vector<std::vector<int> > >
- operator()(const params_t&) noexcept = 0;
+ //returns a std::pair of a std::vector of drums in the system and an adjacency vector of size_t.
+ //the last drum (.back()) is the drum seen when no neighbour is present.
+ //TODO: later, overload with a vector of drum parameters to set each drum's properties individually
+ virtual std::pair<std::vector<Drum<value_t, params_t, vars_t, sbuffer_t> >, std::vector<std::vector<int> > >
+ operator()(const params_t&) noexcept = 0;
};
#endif
diff --git a/projects/braidingTightBinding/lib/matrix_element_calculator.hpp b/projects/braidingTightBinding/lib/matrix_element_calculator.hpp
index 05aaa4bcafbeb4874275c853de2893fbf4aadbb9..7f3c53ca76b2abe8643c35a1ef2f2825aaee0c5c 100644
--- a/projects/braidingTightBinding/lib/matrix_element_calculator.hpp
+++ b/projects/braidingTightBinding/lib/matrix_element_calculator.hpp
@@ -3,19 +3,19 @@
template <typename value_t, typename params_t, typename vars_t, typename drum_t>
class MatrixElementCalculator{
- public:
- MatrixElementCalculator() = default;
- ~MatrixElementCalculator() = default;
+public:
+ MatrixElementCalculator() = default;
+ ~MatrixElementCalculator() = default;
- //diagonal element at (index, index)
- virtual value_t operator()(const size_t index, const std::vector<drum_t>& drums) const noexcept = 0;
- //coupling elements at (index1, index2)
- //for neighbour 0
- virtual value_t operator()(const size_t index1, const size_t index2, const std::vector<drum_t>& drums) const noexcept = 0;
- //for neighbour 1
- virtual value_t operator()(const size_t index1, const size_t index2, const std::vector<drum_t>& drums, const int) const noexcept = 0;
- //for neighbour 2
- virtual value_t operator()(const size_t index1, const size_t index2, const std::vector<drum_t>& drums, const int, const int) const noexcept = 0;
+ //diagonal element at (index, index)
+ virtual value_t operator()(const size_t index, const std::vector<drum_t>& drums) const noexcept = 0;
+ //coupling elements at (index1, index2)
+ //for neighbour 0
+ virtual value_t operator()(const size_t index1, const size_t index2, const std::vector<drum_t>& drums) const noexcept = 0;
+ //for neighbour 1
+ virtual value_t operator()(const size_t index1, const size_t index2, const std::vector<drum_t>& drums, const int) const noexcept = 0;
+ //for neighbour 2
+ virtual value_t operator()(const size_t index1, const size_t index2, const std::vector<drum_t>& drums, const int, const int) const noexcept = 0;
};
#endif
diff --git a/projects/braidingTightBinding/lib/rk4_buffer.hpp b/projects/braidingTightBinding/lib/rk4_buffer.hpp
index 5e466227f195f00f9c191801332b5fa115a93b11..3706980bd3706a326090629196cff4a4b7db8a28 100644
--- a/projects/braidingTightBinding/lib/rk4_buffer.hpp
+++ b/projects/braidingTightBinding/lib/rk4_buffer.hpp
@@ -3,13 +3,13 @@
template <typename value_t>
class RK4Buffer{
- public:
- RK4Buffer() = default;
- ~RK4Buffer() = default;
+public:
+ RK4Buffer() = default;
+ ~RK4Buffer() = default;
- public:
- value_t k1, k2, k3, k4; //for velocity
- value_t l1, l2, l3, l4; //for position
+public:
+ value_t k1, k2, k3, k4; //for velocity
+ value_t l1, l2, l3, l4; //for position
};
diff --git a/projects/braidingTightBinding/lib/rk4_stepper.hpp b/projects/braidingTightBinding/lib/rk4_stepper.hpp
index 1951a33ea67a7c0e880961fcd0ce56696c86a31d..61b6b8207d112f73a43d888c26096e1067bae6f9 100644
--- a/projects/braidingTightBinding/lib/rk4_stepper.hpp
+++ b/projects/braidingTightBinding/lib/rk4_stepper.hpp
@@ -1,95 +1,93 @@
#ifndef RK4_STEPPER_HPP_INCLUDED
#define RK4_STEPPER_HPP_INCLUDED
-template <typename value_t, typename params_t, typename vars_t, typename buffer_t,
+template <typename value_t, typename params_t, typename vars_t, typename buffer_t,
typename force_t>
class Rk4Stepper{
- public:
- Rk4Stepper() = default;
- ~Rk4Stepper() = default;
+public:
+ Rk4Stepper() = default;
+ ~Rk4Stepper() = default;
- //time signifies the state of the system, i.e. the time when the complete step was started.
- //operations performed at t
- void step_1(const force_t& force,
- std::vector<Drum<value_t, params_t, vars_t, buffer_t> >& drums,
- const std::vector<std::vector<int> >& adj_vec,
- const value_t dt,
- const value_t time) const noexcept
- {
- //step 1: calculate k1 and l1 for all drums
- for(size_t i = 0; i < drums.size()-1; ++i){
- //TODO: Empty neighbours are now signified by adj_vec[i][j] == -1. We could also push a 0-drum to the end and point to that.
- drums[i].get_sbuffer().k1 = dt*(force(drums[i], drums[adj_vec[i][0]], drums[adj_vec[i][1]], drums[adj_vec[i][2]], time));
- drums[i].get_sbuffer().l1 = dt*drums[i].get_variables().xdot;
- }
- //update x_temp, xdot_temp
- for(size_t i = 0; i < drums.size()-1; ++i){
- drums[i].get_variables().xdot_temp = drums[i].get_variables().xdot + drums[i].get_sbuffer().k1/2.;
- drums[i].get_variables().x_temp = drums[i].get_variables().x + drums[i].get_sbuffer().l1/2.;
- }
- //now it's time for the owner to update the coupler and driver to t+dt/2.
+ //time signifies the state of the system, i.e. the time when the complete step was started.
+ //operations performed at t
+ void step_1(const force_t& force,
+ std::vector<Drum<value_t, params_t, vars_t, buffer_t> >& drums,
+ const std::vector<std::vector<int> >& adj_vec,
+ const value_t dt,
+ const value_t time) const noexcept
+ {
+ //step 1: calculate k1 and l1 for all drums
+ for(size_t i = 0; i < drums.size()-1; ++i){
+ //TODO: Empty neighbours are now signified by adj_vec[i][j] == -1. We could also push a 0-drum to the end and point to that.
+ drums[i].get_sbuffer().k1 = dt*(force(drums[i], drums[adj_vec[i][0]], drums[adj_vec[i][1]], drums[adj_vec[i][2]], time));
+ drums[i].get_sbuffer().l1 = dt*drums[i].get_variables().xdot;
+ }
+ //update x_temp, xdot_temp
+ for(size_t i = 0; i < drums.size()-1; ++i){
+ drums[i].get_variables().xdot_temp = drums[i].get_variables().xdot + drums[i].get_sbuffer().k1/2.;
+ drums[i].get_variables().x_temp = drums[i].get_variables().x + drums[i].get_sbuffer().l1/2.;
+ }
+ //now it's time for the owner to update the coupler and driver to t+dt/2.
+ }
- }
+ //operations performed at t+dt/2
+ void step_2(const force_t& force,
+ std::vector<Drum<value_t, params_t, vars_t, buffer_t> >& drums,
+ const std::vector<std::vector<int> >& adj_vec,
+ const value_t dt,
+ const value_t time) const noexcept
+ {
+ //step 2: calculate k2 and l2 for all drums
+ for(size_t i = 0; i < drums.size()-1; ++i){
+ //TODO: Empty neighbours are now signified by adj_vec[i][j] == -1. We could also push a 0-drum to the end and point to that.
+ drums[i].get_sbuffer().k2 = dt*(force(drums[i], drums[adj_vec[i][0]], drums[adj_vec[i][1]], drums[adj_vec[i][2]], time));
+ drums[i].get_sbuffer().l2 = dt*drums[i].get_variables().xdot_temp;
+ }
+ //update x_temp, xdot_temp
+ for(size_t i = 0; i < drums.size()-1; ++i){
+ drums[i].get_variables().xdot_temp = drums[i].get_variables().xdot + drums[i].get_sbuffer().k2/2.;
+ drums[i].get_variables().x_temp = drums[i].get_variables().x + drums[i].get_sbuffer().l2/2.;
+ }
- //operations performed at t+dt/2
- void step_2(const force_t& force,
- std::vector<Drum<value_t, params_t, vars_t, buffer_t> >& drums,
- const std::vector<std::vector<int> >& adj_vec,
- const value_t dt,
- const value_t time) const noexcept
- {
- //step 2: calculate k2 and l2 for all drums
- for(size_t i = 0; i < drums.size()-1; ++i){
- //TODO: Empty neighbours are now signified by adj_vec[i][j] == -1. We could also push a 0-drum to the end and point to that.
- drums[i].get_sbuffer().k2 = dt*(force(drums[i], drums[adj_vec[i][0]], drums[adj_vec[i][1]], drums[adj_vec[i][2]], time));
- drums[i].get_sbuffer().l2 = dt*drums[i].get_variables().xdot_temp;
- }
- //update x_temp, xdot_temp
- for(size_t i = 0; i < drums.size()-1; ++i){
- drums[i].get_variables().xdot_temp = drums[i].get_variables().xdot + drums[i].get_sbuffer().k2/2.;
- drums[i].get_variables().x_temp = drums[i].get_variables().x + drums[i].get_sbuffer().l2/2.;
- }
+ //step 3: calculate k3 and l3 for all drums
+ for(size_t i = 0; i < drums.size()-1; ++i){
+ //TODO: Empty neighbours are now signified by adj_vec[i][j] == -1. We could also push a 0-drum to the end and point to that.
+ drums[i].get_sbuffer().k3 = dt*(force(drums[i], drums[adj_vec[i][0]], drums[adj_vec[i][1]], drums[adj_vec[i][2]], time));
+ drums[i].get_sbuffer().l3 = dt*drums[i].get_variables().xdot_temp;
+ }
+ //update x_temp, xdot_temp
+ for(size_t i = 0; i < drums.size()-1; ++i){
+ drums[i].get_variables().xdot_temp = drums[i].get_variables().xdot + drums[i].get_sbuffer().k3;
+ drums[i].get_variables().x_temp = drums[i].get_variables().x + drums[i].get_sbuffer().l3;
+ }
+ //now it's time for the owner to update the coupler and driver to t+dt/2.
+ }
- //step 3: calculate k3 and l3 for all drums
- for(size_t i = 0; i < drums.size()-1; ++i){
- //TODO: Empty neighbours are now signified by adj_vec[i][j] == -1. We could also push a 0-drum to the end and point to that.
- drums[i].get_sbuffer().k3 = dt*(force(drums[i], drums[adj_vec[i][0]], drums[adj_vec[i][1]], drums[adj_vec[i][2]], time));
- drums[i].get_sbuffer().l3 = dt*drums[i].get_variables().xdot_temp;
- }
- //update x_temp, xdot_temp
- for(size_t i = 0; i < drums.size()-1; ++i){
- drums[i].get_variables().xdot_temp = drums[i].get_variables().xdot + drums[i].get_sbuffer().k3;
- drums[i].get_variables().x_temp = drums[i].get_variables().x + drums[i].get_sbuffer().l3;
- }
- //now it's time for the owner to update the coupler and driver to t+dt/2.
- }
-
-
- //operations performed at t+dt
- void step_3(const force_t& force,
- std::vector<Drum<value_t, params_t, vars_t, buffer_t> >& drums,
- const std::vector<std::vector<int> >& adj_vec,
- const value_t dt,
- const value_t time) const noexcept
- {
- //step 4: calculate k4 and l4 for all drums
- for(size_t i = 0; i < drums.size()-1; ++i){
- //TODO: Empty neighbours are now signified by adj_vec[i][j] == -1. We could also push a 0-drum to the end and point to that.
- drums[i].get_sbuffer().k4 = dt*(force(drums[i], drums[adj_vec[i][0]], drums[adj_vec[i][1]], drums[adj_vec[i][2]], time));
- drums[i].get_sbuffer().l4 = dt*drums[i].get_variables().xdot_temp;
- }
- //update x, xdot, x_temp, xdot_temp to t+dt
- for(size_t i = 0; i < drums.size()-1; ++i){
- drums[i].get_variables().xdot_temp = drums[i].get_variables().xdot + drums[i].get_sbuffer().k4;
- drums[i].get_variables().x_temp = drums[i].get_variables().x + drums[i].get_sbuffer().l4;
- drums[i].get_variables().xdot = drums[i].get_variables().xdot + drums[i].get_sbuffer().k4;
- drums[i].get_variables().x = drums[i].get_variables().x + drums[i].get_sbuffer().l4;
- }
- }
+ //operations performed at t+dt
+ void step_3(const force_t& force,
+ std::vector<Drum<value_t, params_t, vars_t, buffer_t> >& drums,
+ const std::vector<std::vector<int> >& adj_vec,
+ const value_t dt,
+ const value_t time) const noexcept
+ {
+ //step 4: calculate k4 and l4 for all drums
+ for(size_t i = 0; i < drums.size()-1; ++i){
+ //TODO: Empty neighbours are now signified by adj_vec[i][j] == -1. We could also push a 0-drum to the end and point to that.
+ drums[i].get_sbuffer().k4 = dt*(force(drums[i], drums[adj_vec[i][0]], drums[adj_vec[i][1]], drums[adj_vec[i][2]], time));
+ drums[i].get_sbuffer().l4 = dt*drums[i].get_variables().xdot_temp;
+ }
+ //update x, xdot, x_temp, xdot_temp to t+dt
+ for(size_t i = 0; i < drums.size()-1; ++i){
+ drums[i].get_variables().xdot_temp = drums[i].get_variables().xdot + drums[i].get_sbuffer().k4;
+ drums[i].get_variables().x_temp = drums[i].get_variables().x + drums[i].get_sbuffer().l4;
+ drums[i].get_variables().xdot = drums[i].get_variables().xdot + drums[i].get_sbuffer().k4;
+ drums[i].get_variables().x = drums[i].get_variables().x + drums[i].get_sbuffer().l4;
+ }
+ }
};
#endif
diff --git a/projects/braidingTightBinding/lib/system.hpp b/projects/braidingTightBinding/lib/system.hpp
index 37855c12fa61b561f5846b77ecca5f3b1a7f452b..53788862242ed9273b91ce94a959d7d9a85e7fd7 100644
--- a/projects/braidingTightBinding/lib/system.hpp
+++ b/projects/braidingTightBinding/lib/system.hpp
@@ -5,111 +5,111 @@
template<typename value_t, typename drum_t, typename grabber_t, typename sysparams_t, typename force_t, typename coupler_t, typename driver_t, typename stepper_t, typename matelecalc_t>
class System{
- public:
- System(const value_t t_end, const value_t dt, const std::vector<drum_t>& drums, const stepper_t stepper, const force_t force, const sysparams_t sysparams, const grabber_t grabber)
- : drums_(drums), stepper_(stepper), force_(force), sysparams_(sysparams), grabber_(grabber), t_end_(t_end), dt_(dt), time_(0.)
- {
- sysparams_.coupler.precompute(t_end, dt, drums);
- sysparams_.driver.precompute(t_end, dt, drums);
- grabber_.init(t_end, dt, drums, sysparams.adjacency_vector);
-
- push_dc(); //push the initial values
-
- }
- ~System() = default;
-
- void simulate(){
- while(time_ <= t_end_){
- grabber_.grab(drums_, time_);
- step();
- }
- }
-
- void step(){
- push_dc();
- stepper_.step_1(force_, drums_, sysparams_.adjacency_vector, dt_, time_);
- step_dc(dt_/2.);
- push_dc();
- stepper_.step_2(force_, drums_, sysparams_.adjacency_vector, dt_, time_);
- step_dc(dt_/2.);
- push_dc();
- stepper_.step_3(force_, drums_, sysparams_.adjacency_vector, dt_, time_);
-
- time_ += dt_;
- }
-
- void reset_time(){
- time_ = value_t(0.);
- }
-
- void set_step(value_t dt){
- dt_ = dt;
- }
-
- bool save(){
- return grabber_.save();
- }
-
- std::vector<value_t> get_matrix(const matelecalc_t& mec){
- calculate_matrix(mec);
- return matrix_;
- }
-
- private:
- void step_dc(const value_t dt){
- sysparams_.coupler.step(dt);
- sysparams_.driver.step(dt);
- }
-
- void push_dc(){
- for(size_t i = 0; i < drums_.size()-1; ++i){
- drums_[i].get_variables().V = sysparams_.driver(i);
- drums_[i].get_variables().t0 = sysparams_.coupler(i,0);
- drums_[i].get_variables().t1 = sysparams_.coupler(i,1);
- drums_[i].get_variables().t2 = sysparams_.coupler(i,2);
- }
- }
-
- //Calculates dynamical matrix for the currently pushed values
- void calculate_matrix(const matelecalc_t& mec){
- //for convenience
- size_t N = drums_.size()-1;
- //clear the matrix
- matrix_.clear();
- //prepare memory
- matrix_.reserve((drums_.size()-1)*(drums_.size()-1) - matrix_.capacity());
- //initialize
- for(size_t i = 0; i < (drums_.size()-1) * (drums_.size()-1); ++i)
- matrix_.push_back(value_t(0.));
-
- for(size_t i = 0; i < drums_.size()-1; ++i){
- size_t n0_index = sysparams_.adjacency_vector[i][0];
- size_t n1_index = sysparams_.adjacency_vector[i][1];
- size_t n2_index = sysparams_.adjacency_vector[i][2];
- //diagonal
- matrix_[N*i + i] = mec(i, drums_);
- //couplings
- //neighbour 0
- if(n0_index != drums_.size()-1) [[likely]] //actually a neighbour
- matrix_[N*i + n0_index] = mec(i, n0_index, drums_);
- //neighbour 1
- if(n1_index != drums_.size()-1) [[likely]] //actually a neighbour
- matrix_[N*i + n1_index] = mec(i, n1_index, drums_, 1);
- //neighbour 2
- if(n2_index != drums_.size()-1) [[likely]] //actually a neighbour
- matrix_[N*i + n2_index] = mec(i, n2_index, drums_, 2, 2);
- }
- }
-
- std::vector<drum_t> drums_; //The last drum is NOT TO BE TOUCHED!
- sysparams_t sysparams_;
- grabber_t grabber_;
- stepper_t stepper_;
- force_t force_;
- value_t dt_, t_end_;
- value_t time_;
-
- std::vector<value_t> matrix_; //Matrix of this problem
+public:
+ System(const value_t t_end, const value_t dt, const std::vector<drum_t>& drums, const stepper_t stepper, const force_t force, const sysparams_t sysparams, const grabber_t grabber)
+ : drums_(drums), stepper_(stepper), force_(force), sysparams_(sysparams), grabber_(grabber), t_end_(t_end), dt_(dt), time_(0.)
+ {
+ sysparams_.coupler.precompute(t_end, dt, drums);
+ sysparams_.driver.precompute(t_end, dt, drums);
+ grabber_.init(t_end, dt, drums, sysparams.adjacency_vector);
+
+ push_dc(); //push the initial values
+
+ }
+ ~System() = default;
+
+ void simulate(){
+ while(time_ <= t_end_){
+ grabber_.grab(drums_, time_);
+ step();
+ }
+ }
+
+ void step(){
+ push_dc();
+ stepper_.step_1(force_, drums_, sysparams_.adjacency_vector, dt_, time_);
+ step_dc(dt_/2.);
+ push_dc();
+ stepper_.step_2(force_, drums_, sysparams_.adjacency_vector, dt_, time_);
+ step_dc(dt_/2.);
+ push_dc();
+ stepper_.step_3(force_, drums_, sysparams_.adjacency_vector, dt_, time_);
+
+ time_ += dt_;
+ }
+
+ void reset_time(){
+ time_ = value_t(0.);
+ }
+
+ void set_step(value_t dt){
+ dt_ = dt;
+ }
+
+ bool save(){
+ return grabber_.save();
+ }
+
+ std::vector<value_t> get_matrix(const matelecalc_t& mec){
+ calculate_matrix(mec);
+ return matrix_;
+ }
+
+private:
+ void step_dc(const value_t dt){
+ sysparams_.coupler.step(dt);
+ sysparams_.driver.step(dt);
+ }
+
+ void push_dc(){
+ for(size_t i = 0; i < drums_.size()-1; ++i){
+ drums_[i].get_variables().V = sysparams_.driver(i);
+ drums_[i].get_variables().t0 = sysparams_.coupler(i,0);
+ drums_[i].get_variables().t1 = sysparams_.coupler(i,1);
+ drums_[i].get_variables().t2 = sysparams_.coupler(i,2);
+ }
+ }
+
+ //Calculates dynamical matrix for the currently pushed values
+ void calculate_matrix(const matelecalc_t& mec){
+ //for convenience
+ size_t N = drums_.size()-1;
+ //clear the matrix
+ matrix_.clear();
+ //prepare memory
+ matrix_.reserve((drums_.size()-1)*(drums_.size()-1) - matrix_.capacity());
+ //initialize
+ for(size_t i = 0; i < (drums_.size()-1) * (drums_.size()-1); ++i)
+ matrix_.push_back(value_t(0.));
+
+ for(size_t i = 0; i < drums_.size()-1; ++i){
+ size_t n0_index = sysparams_.adjacency_vector[i][0];
+ size_t n1_index = sysparams_.adjacency_vector[i][1];
+ size_t n2_index = sysparams_.adjacency_vector[i][2];
+ //diagonal
+ matrix_[N*i + i] = mec(i, drums_);
+ //couplings
+ //neighbour 0
+ if(n0_index != drums_.size()-1) [[likely]] //actually a neighbour
+ matrix_[N*i + n0_index] = mec(i, n0_index, drums_);
+ //neighbour 1
+ if(n1_index != drums_.size()-1) [[likely]] //actually a neighbour
+ matrix_[N*i + n1_index] = mec(i, n1_index, drums_, 1);
+ //neighbour 2
+ if(n2_index != drums_.size()-1) [[likely]] //actually a neighbour
+ matrix_[N*i + n2_index] = mec(i, n2_index, drums_, 2, 2);
+ }
+ }
+
+ std::vector<drum_t> drums_; //The last drum is NOT TO BE TOUCHED!
+ sysparams_t sysparams_;
+ grabber_t grabber_;
+ stepper_t stepper_;
+ force_t force_;
+ value_t dt_, t_end_;
+ value_t time_;
+
+ std::vector<value_t> matrix_; //Matrix of this problem
};
#endif
diff --git a/projects/braidingTightBinding/lib/system_parameters.hpp b/projects/braidingTightBinding/lib/system_parameters.hpp
index 637d7792be595c26bc2c0eeb8a89652b5dbc19be..8c904899061f5d8c4bd2eefafe48451913406198 100644
--- a/projects/braidingTightBinding/lib/system_parameters.hpp
+++ b/projects/braidingTightBinding/lib/system_parameters.hpp
@@ -4,15 +4,15 @@
template <typename coupler_t, typename driver_t>
class SystemParameters{
- public:
- //constructors
- SystemParameters(coupler_t coupler, driver_t driver, std::vector< std::vector<int> > adjacency_vector) noexcept: coupler(coupler), driver(driver), adjacency_vector(adjacency_vector) {}
+public:
+ //constructors
+ SystemParameters(coupler_t coupler, driver_t driver, std::vector< std::vector<int> > adjacency_vector) noexcept: coupler(coupler), driver(driver), adjacency_vector(adjacency_vector) {}
-
- public:
- coupler_t coupler;
- driver_t driver;
- std::vector< std::vector<int> > adjacency_vector;
+
+public:
+ coupler_t coupler;
+ driver_t driver;
+ std::vector< std::vector<int> > adjacency_vector;
};
#endif
diff --git a/projects/braidingTightBinding/lib/vec2.hpp b/projects/braidingTightBinding/lib/vec2.hpp
index a2bb4c3aabc8ba0cace7bc557ab5d14c4afd1287..9f9d5e0f5733bcbb07cf7615201f15d504707d45 100644
--- a/projects/braidingTightBinding/lib/vec2.hpp
+++ b/projects/braidingTightBinding/lib/vec2.hpp
@@ -5,146 +5,146 @@
template <typename value_t>
class Vec2{
- public:
- //constructors
- Vec2(const value_t x, const value_t y) noexcept: x_(x), y_(y) {}
- Vec2(const Vec2<value_t>& other) = default;
- Vec2& operator=(const Vec2&) = default;
- Vec2() = default;
- ~Vec2() = default;
-
-
- //access
- value_t x() const noexcept { return x_; }
- value_t y() const noexcept { return y_; }
- Vec2 normalized() const { return *this/this->norm(); }
-
- //TODO: Precalculate these and flag when recomputation is needed
- value_t r() const noexcept { return this->norm(); }
- value_t phi() const
- {
- value_t angle (0);
- angle = std::atan2(y_, x_);
- if(std::isnan(angle)){
- throw("ATAN2 NAN OCCURRED");
- return 0.;
- }
- return angle;
- }
-
- //function members
- value_t r_wrt(const Vec2& origin) const noexcept { return (*this - origin).r(); }
- value_t phi_wrt(const Vec2& origin) const { return (*this - origin).phi(); }
-
- value_t norm() const noexcept { return std::sqrt(x_*x_ + y_*y_); }
- value_t norm_sq() const noexcept { return x_*x_ + y_* y_; }
-
- //modifiers
- Vec2& rotate(const Vec2& center, const value_t degrees)
- {
- *this -= center;
- value_t temp = x_;
- const value_t radians = degrees*3.141592653589793/180.;
-
- x_ = std::cos(radians)*x_ - std::sin(radians)*y_;
- y_ = std::sin(radians)*temp + std::cos(radians)*y_;
-
- *this += center;
-
- return *this;
- }
-
- Vec2& normalize()
- {
- value_t length = norm();
- if(length == value_t(0))
- throw("NORMALIZE ZERO VEC2 ATTEMPTED");
- x_ /= length;
- y_ /= length;
- return *this;
- }
-
- //overloading
- Vec2& operator+=(const Vec2& rhs) noexcept
- {
- x_ += rhs.x_;
- y_ += rhs.y_;
- return *this;
- }
- Vec2& operator-=(const Vec2& rhs) noexcept
- {
- x_ -= rhs.x_;
- y_ -= rhs.y_;
- return *this;
- }
- Vec2& operator*=(const value_t s) noexcept
- {
- x_ *= s;
- y_ *= s;
- return *this;
- }
- Vec2& operator/=(const value_t s)
- {
- if(x_ == value_t(0) and y_ == value_t(0))
- return *this;
- else if(s == value_t(0))
- throw("DIVISION BY ZERO ATTEMPTED");
- x_ /= s;
- y_ /= s;
- return *this;
- }
-
- value_t& operator[](std::size_t i) noexcept
- {
- if(i == 0)
- return x_;
- else
- return y_;
- }
- value_t operator[](std::size_t i) const noexcept
- {
- if(i == 0)
- return x_;
- else
- return y_;
- }
-
- //non-ADL hidden friends
- friend Vec2 operator+(const Vec2& lhs, Vec2 rhs) noexcept
- {
- rhs += lhs;
- return rhs;
- }
- friend Vec2 operator-(Vec2 lhs, const Vec2& rhs) noexcept
- {
- lhs -= rhs;
- return lhs;
- }
- friend Vec2 operator*(const value_t s, Vec2 v) noexcept
- {
- v *= s;
- return v;
- }
- friend Vec2 operator*(Vec2 v, const value_t s) noexcept
- {
- v *= s;
- return v;
- }
- friend Vec2 operator/(Vec2 v, const value_t s)
- {
- v /= s;
- return v;
- }
- friend value_t operator*(const Vec2& v, const Vec2& w) noexcept { return v.x_*w.x_ + v.y_*w.y_; }//dot product
-
- private:
- value_t x_, y_;
+public:
+ //constructors
+ Vec2(const value_t x, const value_t y) noexcept: x_(x), y_(y) {}
+ Vec2(const Vec2<value_t>& other) = default;
+ Vec2& operator=(const Vec2&) = default;
+ Vec2() = default;
+ ~Vec2() = default;
+
+
+ //access
+ value_t x() const noexcept { return x_; }
+ value_t y() const noexcept { return y_; }
+ Vec2 normalized() const { return *this/this->norm(); }
+
+ //TODO: Precalculate these and flag when recomputation is needed
+ value_t r() const noexcept { return this->norm(); }
+ value_t phi() const
+ {
+ value_t angle (0);
+ angle = std::atan2(y_, x_);
+ if(std::isnan(angle)){
+ throw("ATAN2 NAN OCCURRED");
+ return 0.;
+ }
+ return angle;
+ }
+
+ //function members
+ value_t r_wrt(const Vec2& origin) const noexcept { return (*this - origin).r(); }
+ value_t phi_wrt(const Vec2& origin) const { return (*this - origin).phi(); }
+
+ value_t norm() const noexcept { return std::sqrt(x_*x_ + y_*y_); }
+ value_t norm_sq() const noexcept { return x_*x_ + y_* y_; }
+
+ //modifiers
+ Vec2& rotate(const Vec2& center, const value_t degrees)
+ {
+ *this -= center;
+ value_t temp = x_;
+ const value_t radians = degrees*3.141592653589793/180.;
+
+ x_ = std::cos(radians)*x_ - std::sin(radians)*y_;
+ y_ = std::sin(radians)*temp + std::cos(radians)*y_;
+
+ *this += center;
+
+ return *this;
+ }
+
+ Vec2& normalize()
+ {
+ value_t length = norm();
+ if(length == value_t(0))
+ throw("NORMALIZE ZERO VEC2 ATTEMPTED");
+ x_ /= length;
+ y_ /= length;
+ return *this;
+ }
+
+ //overloading
+ Vec2& operator+=(const Vec2& rhs) noexcept
+ {
+ x_ += rhs.x_;
+ y_ += rhs.y_;
+ return *this;
+ }
+ Vec2& operator-=(const Vec2& rhs) noexcept
+ {
+ x_ -= rhs.x_;
+ y_ -= rhs.y_;
+ return *this;
+ }
+ Vec2& operator*=(const value_t s) noexcept
+ {
+ x_ *= s;
+ y_ *= s;
+ return *this;
+ }
+ Vec2& operator/=(const value_t s)
+ {
+ if(x_ == value_t(0) and y_ == value_t(0))
+ return *this;
+ else if(s == value_t(0))
+ throw("DIVISION BY ZERO ATTEMPTED");
+ x_ /= s;
+ y_ /= s;
+ return *this;
+ }
+
+ value_t& operator[](std::size_t i) noexcept
+ {
+ if(i == 0)
+ return x_;
+ else
+ return y_;
+ }
+ value_t operator[](std::size_t i) const noexcept
+ {
+ if(i == 0)
+ return x_;
+ else
+ return y_;
+ }
+
+ //non-ADL hidden friends
+ friend Vec2 operator+(const Vec2& lhs, Vec2 rhs) noexcept
+ {
+ rhs += lhs;
+ return rhs;
+ }
+ friend Vec2 operator-(Vec2 lhs, const Vec2& rhs) noexcept
+ {
+ lhs -= rhs;
+ return lhs;
+ }
+ friend Vec2 operator*(const value_t s, Vec2 v) noexcept
+ {
+ v *= s;
+ return v;
+ }
+ friend Vec2 operator*(Vec2 v, const value_t s) noexcept
+ {
+ v *= s;
+ return v;
+ }
+ friend Vec2 operator/(Vec2 v, const value_t s)
+ {
+ v /= s;
+ return v;
+ }
+ friend value_t operator*(const Vec2& v, const Vec2& w) noexcept { return v.x_*w.x_ + v.y_*w.y_; }//dot product
+
+private:
+ value_t x_, y_;
};
template <typename value_t>
std::ostream& operator<<(std::ostream& os, const Vec2<value_t>& v)
{
- return os << "(" << v.x() << ", " << v.y() << ")";
+ return os << "(" << v.x() << ", " << v.y() << ")";
}
#endif