Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • engelerp/rbcomb-simulation
1 result
Show changes
Hide whitespace changes
Inline Side-by-side
Showing
with 422 additions and 424 deletions
projects/braidingTightBinding/lib/driver.hpp
View file @ 5f273894
......@@ -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
projects/braidingTightBinding/lib/drum.hpp
View file @ 5f273894
......@@ -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
projects/braidingTightBinding/lib/force.hpp
View file @ 5f273894
......@@ -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
projects/braidingTightBinding/lib/lattice_generator.hpp
View file @ 5f273894
......@@ -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
projects/braidingTightBinding/lib/matrix_element_calculator.hpp
View file @ 5f273894
......@@ -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
projects/braidingTightBinding/lib/rk4_buffer.hpp
View file @ 5f273894
......@@ -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
};
......
projects/braidingTightBinding/lib/rk4_stepper.hpp
View file @ 5f273894
#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
projects/braidingTightBinding/lib/system.hpp
View file @ 5f273894
......@@ -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
projects/braidingTightBinding/lib/system_parameters.hpp
View file @ 5f273894
......@@ -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
projects/braidingTightBinding/lib/vec2.hpp
View file @ 5f273894
......@@ -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