diff --git a/projects/braidingTightBinding/include/coupler_braid.hpp b/projects/braidingTightBinding/include/coupler_braid.hpp
index 18e909d193a6b311b27dce5e7204a8870678b0f3..7c192de4840f9155efe94c42afa4fc60da25a945 100644
--- a/projects/braidingTightBinding/include/coupler_braid.hpp
+++ b/projects/braidingTightBinding/include/coupler_braid.hpp
@@ -6,55 +6,57 @@
template <typename value_t, typename drum_t, typename params_t>
class CouplerBraid: public Coupler<value_t, drum_t>{
- public:
- //constructors
- CouplerBraid(const std::vector<Vortex<value_t> >& vortices): vortices_(vortices) {}
- ~CouplerBraid() = default;
-
- void precompute(const value_t t_end, const value_t dt, const std::vector<drum_t>& drum_vec) noexcept final override
- {
- //we need the drum parameters later on, so extract it
- params_.reserve(drum_vec.size());
- for(auto d: drum_vec)
- params_.push_back(d.get_parameters());
-
- //we also need the neighbour vectors nicely indexed
- for(auto d: drum_vec){
- if(d.get_parameters().sublattice == 'A'){
- n_vecs_.push_back(d.get_parameters().s0);
- n_vecs_.push_back(d.get_parameters().s1);
- n_vecs_.push_back(d.get_parameters().s2);
- break;
- }
- }
- }
- void step(value_t dt) noexcept final
- {
- //move vortices or alpha here
- }
- value_t operator()(const size_t drum_index, const size_t neighbour_index) const noexcept final override
- {
- Vec2<value_t> v = params_[drum_index].position;
- Vec2<value_t> s = n_vecs_[neighbour_index];
-
- //the formula is only valid on the 'A' sublattice
- if(params_[drum_index].sublattice == 'B')
- v = v - s;
-
- std::complex<value_t> ret_val = 1.; //vortices are multiplicative
- for(auto tex: vortices_){
- ret_val *= tex.distortion(v);
- }
- //add kekule
- ret_val *= vortices_.begin()->kekule(v, s);
-
- return 1. + std::real<value_t>(ret_val); //return with added offset, as in eliska's code
- }
-
- private:
- std::vector<Vortex<value_t> > vortices_;
- std::vector<params_t> params_;
- std::vector<Vec2<value_t> > n_vecs_;
+public:
+ //constructors
+ CouplerBraid(const std::vector<Vortex<value_t> >& vortices): vortices_(vortices) {}
+ ~CouplerBraid() = default;
+
+ void precompute(const value_t t_end, const value_t dt, const std::vector<drum_t>& drum_vec) noexcept final override
+ {
+ //we need the drum parameters later on, so extract it
+ params_.reserve(drum_vec.size());
+ for(auto d: drum_vec){
+ params_.push_back(d.get_parameters());
+ }
+
+ //we also need the neighbour vectors nicely indexed
+ for(auto d: drum_vec){
+ if(d.get_parameters().sublattice == 'A'){
+ n_vecs_.push_back(d.get_parameters().s0);
+ n_vecs_.push_back(d.get_parameters().s1);
+ n_vecs_.push_back(d.get_parameters().s2);
+ break;
+ }
+ }
+ }
+ void step(value_t dt) noexcept final
+ {
+ //move vortices or alpha here
+ }
+ value_t operator()(const size_t drum_index, const size_t neighbour_index) const noexcept final override
+ {
+ Vec2<value_t> v = params_[drum_index].position;
+ Vec2<value_t> s = n_vecs_[neighbour_index];
+
+ //the formula is only valid on the 'A' sublattice
+ if(params_[drum_index].sublattice == 'B'){
+ v = v - s;
+ }
+
+ std::complex<value_t> ret_val = 1.; //vortices are multiplicative
+ for(auto tex: vortices_){
+ ret_val *= tex.distortion(v);
+ }
+ //add kekule
+ ret_val *= vortices_.begin()->kekule(v, s);
+
+ return 1. + std::real<value_t>(ret_val); //return with added offset, as in eliska's code
+ }
+
+private:
+ std::vector<Vortex<value_t> > vortices_;
+ std::vector<params_t> params_;
+ std::vector<Vec2<value_t> > n_vecs_;
};
#endif
diff --git a/projects/braidingTightBinding/include/coupler_const.hpp b/projects/braidingTightBinding/include/coupler_const.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..8dfea90f2ecd13a015430f9fd8506423ff568d66
--- /dev/null
+++ b/projects/braidingTightBinding/include/coupler_const.hpp
@@ -0,0 +1,30 @@
+#ifndef COUPLER_CONST_HPP_INCLUDED
+#define COUPLER_CONST_HPP_INCLUDED
+#include <vector>
+#include <coupler.hpp>
+
+template <typename value_t, typename drum_t, typename params_t>
+class CouplerConst: public Coupler<value_t, drum_t>{
+public:
+ //constructors
+ CouplerConst(value_t value): value_(value) {}
+ ~CouplerConst() = default;
+
+ void precompute(const value_t t_end, const value_t dt, const std::vector<drum_t>& drum_vec) noexcept final override
+ {
+ //nothing to do
+ }
+ void step(value_t dt) noexcept final
+ {
+ //nothing to do
+ }
+ value_t operator()(const size_t drum_index, const size_t neighbour_index) const noexcept final override
+ {
+ return value_;
+ }
+
+private:
+ value_t value_;
+};
+
+#endif
diff --git a/projects/braidingTightBinding/include/driver_braid.hpp b/projects/braidingTightBinding/include/driver_braid.hpp
index e36a2adf4367a7224cb64f14767b5e9f33a87c3d..7ae77aaf6b02723feb32a6cc49664524397e9460 100644
--- a/projects/braidingTightBinding/include/driver_braid.hpp
+++ b/projects/braidingTightBinding/include/driver_braid.hpp
@@ -5,23 +5,23 @@
template <typename value_t, typename drum_t>
class DriverBraid: public Driver<value_t, drum_t>{
- public:
- //constructors
- DriverBraid() = default;
- ~DriverBraid() = default;
+public:
+ //constructors
+ DriverBraid() = default;
+ ~DriverBraid() = default;
- void precompute(const value_t t_end, const value_t dt, const std::vector<drum_t>& drum_vec) noexcept final override
- {
+ void precompute(const value_t t_end, const value_t dt, const std::vector<drum_t>& drum_vec) noexcept final override
+ {
- }
- void step(value_t dt) noexcept final
- {
+ }
+ void step(value_t dt) noexcept final
+ {
- }
- value_t operator()(const size_t drum_index) const noexcept final override
- {
- return 0;
- }
+ }
+ value_t operator()(const size_t drum_index) const noexcept final override
+ {
+ return 0;
+ }
};
#endif
diff --git a/projects/braidingTightBinding/include/drum_parameters_braid.hpp b/projects/braidingTightBinding/include/drum_parameters_braid.hpp
index 7211663064acbbe15b6d87708e93c596bbedc989..de4c4b8a3b81bc796511f9f10d4f8c7b98db97d4 100644
--- a/projects/braidingTightBinding/include/drum_parameters_braid.hpp
+++ b/projects/braidingTightBinding/include/drum_parameters_braid.hpp
@@ -7,58 +7,59 @@
//TODO: Create an interface. It would be a trivial one, though.
template <typename value_t>
class DrumParametersBraid{
- public:
- /* Arguments:
- * k0: diagonal matrix element
- * k1: coupling matrix element prefactor 1
- * k2: coupling matrix element prefactor 2
- * c: damping coefficient
- * position: position of the drum in real space
- * sublattice: 'A'/'a' or 'B'/'b', identifies sublattice drum is part of
- */
- DrumParametersBraid(const value_t k0, const value_t k1, const value_t k2, const value_t c, const Vec2<value_t>& position, char sublattice) noexcept
- : k0(k0), k1(k1), k2(k2), c(c), position(position), sublattice(sublattice)
- {
- generate_sublattice(sublattice, value_t(1.));
- }
- DrumParametersBraid() = delete; //no default initialization allowed
- DrumParametersBraid(const DrumParametersBraid&) = default;
- DrumParametersBraid& operator=(const DrumParametersBraid&) = default;
- ~DrumParametersBraid() = default;
+public:
+ /* Arguments:
+ * k0: diagonal matrix element
+ * k1: coupling matrix element prefactor 1
+ * k2: coupling matrix element prefactor 2
+ * c: damping coefficient
+ * position: position of the drum in real space
+ * sublattice: 'A'/'a' or 'B'/'b', identifies sublattice drum is part of
+ */
+ DrumParametersBraid(const value_t k0, const value_t k1, const value_t k2, const value_t c, const Vec2<value_t>& position, char sublattice) noexcept
+ : k0(k0), k1(k1), k2(k2), c(c), position(position), sublattice(sublattice)
+ {
+ generate_sublattice(sublattice, value_t(1.));
+ }
+ DrumParametersBraid() = delete; //no default initialization allowed
+ DrumParametersBraid(const DrumParametersBraid&) = default;
+ DrumParametersBraid& operator=(const DrumParametersBraid&) = default;
+ ~DrumParametersBraid() = default;
- private:
- int generate_sublattice(char sublattice, value_t lc) noexcept{
- if(sublattice == 'A' or sublattice == 'a'){
- s0 = Vec2<value_t> (0, -lc);
- s1 = Vec2<value_t> (-std::sqrt(3)*lc/2., lc/2.);
- s2 = Vec2<value_t> (std::sqrt(3)*lc/2., lc/2.);
- return 0;
- }
- else if(sublattice == 'B' or sublattice == 'b'){
- s0 = -1*(Vec2<value_t> (0, -lc));
- s1 = -1*(Vec2<value_t> (-std::sqrt(3)*lc/2., lc/2.));
- s2 = -1*(Vec2<value_t> (std::sqrt(3)*lc/2., lc/2.));
- return 0;
- }
- else
- return -1; //invalid sublattice identifier
- }
+private:
+ int generate_sublattice(char sublattice, value_t lc) noexcept{
+ if(sublattice == 'A' or sublattice == 'a'){
+ s0 = Vec2<value_t> (0, -lc);
+ s1 = Vec2<value_t> (-std::sqrt(3)*lc/2., lc/2.);
+ s2 = Vec2<value_t> (std::sqrt(3)*lc/2., lc/2.);
+ return 0;
+ }
+ else if(sublattice == 'B' or sublattice == 'b'){
+ s0 = -1*(Vec2<value_t> (0, -lc));
+ s1 = -1*(Vec2<value_t> (-std::sqrt(3)*lc/2., lc/2.));
+ s2 = -1*(Vec2<value_t> (std::sqrt(3)*lc/2., lc/2.));
+ return 0;
+ }
+ else{
+ return -1; //invalid sublattice identifier
+ }
+ }
- public:
- value_t k0; //onsite prefactor
- value_t k1; //coupling prefactor 1
- value_t k2; //coupling prefactor 2
- value_t c; //damping coefficient
- char sublattice; //sublattice the drum is part of (governs electrode geometry)
- Vec2<value_t> s0, s1, s2; //vectors connecting to neighbours, ordered according to drum.hpp
- Vec2<value_t> position; //position of the drum in real space
+public:
+ value_t k0; //onsite prefactor
+ value_t k1; //coupling prefactor 1
+ value_t k2; //coupling prefactor 2
+ value_t c; //damping coefficient
+ char sublattice; //sublattice the drum is part of (governs electrode geometry)
+ Vec2<value_t> s0, s1, s2; //vectors connecting to neighbours, ordered according to drum.hpp
+ Vec2<value_t> position; //position of the drum in real space
};
template<typename value_t>
std::ostream& operator<<(std::ostream& os, const DrumParametersBraid<value_t>& dp)
{
- return os << "Sublattice: " << dp.sublattice << ", k0: " << dp.k0 << ", k1: " << dp.k1 << ", k2: " << dp.k2 << std::endl << "s0: " << dp.s0 << ", s1: " << dp.s1 << ", s2: " << dp.s2;
+ return os << "Sublattice: " << dp.sublattice << ", k0: " << dp.k0 << ", k1: " << dp.k1 << ", k2: " << dp.k2 << std::endl << "s0: " << dp.s0 << ", s1: " << dp.s1 << ", s2: " << dp.s2;
}
diff --git a/projects/braidingTightBinding/include/drum_variables_braid.hpp b/projects/braidingTightBinding/include/drum_variables_braid.hpp
index d8a7170ef82ca03abeed3fa3d8bbfc8032146e73..00f4f9e1ba5f709c5ac807b5601ff17f5dfb2274 100644
--- a/projects/braidingTightBinding/include/drum_variables_braid.hpp
+++ b/projects/braidingTightBinding/include/drum_variables_braid.hpp
@@ -8,16 +8,16 @@
//TODO: Generate an interface. It would be trivial, though.
template <typename value_t>
class DrumVariablesBraid{
- public:
- //constructors
- DrumVariablesBraid() = default;
- DrumVariablesBraid& operator=(const DrumVariablesBraid&) = default;
- ~DrumVariablesBraid() = default;
+public:
+ //constructors
+ DrumVariablesBraid() = default;
+ DrumVariablesBraid& operator=(const DrumVariablesBraid&) = default;
+ ~DrumVariablesBraid() = default;
- value_t t0, t1, t2; //t + dt, i.e. baseline plus vortices along bond 0,1,2
- value_t V; //Coupling to central electrode
- value_t x, xdot; //Current Elongation (position) and velocity
- value_t x_temp, xdot_temp; //used as calculation arguments in rk steps
+ value_t t0, t1, t2; //t + dt, i.e. baseline plus vortices along bond 0,1,2
+ value_t V; //Coupling to central electrode
+ value_t x, xdot; //Current Elongation (position) and velocity
+ value_t x_temp, xdot_temp; //used as calculation arguments in rk steps
};
#endif
diff --git a/projects/braidingTightBinding/include/force_braid.hpp b/projects/braidingTightBinding/include/force_braid.hpp
index 62748139d9021ddc31906ff8fcc7273a3e4c9133..da25331aa2b63559fa434adbbd166cb11a86648e 100644
--- a/projects/braidingTightBinding/include/force_braid.hpp
+++ b/projects/braidingTightBinding/include/force_braid.hpp
@@ -9,35 +9,35 @@
template <typename value_t, typename params_t, typename vars_t, typename buffer_t>
class ForceBraid: public Force<value_t, params_t, vars_t, buffer_t>{
- public:
- ForceBraid() = default;
- ~ForceBraid() = default;
+public:
+ ForceBraid() = default;
+ ~ForceBraid() = default;
- 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 final override
- {
- //fetch data
- //note that no copying is done here, these are all references
- //TODO: implement drive
- const params_t& drum_params = drum.get_parameters(); //parameters of drum
- const vars_t& drum_vars = drum.get_variables(); //variables of drum
- const vars_t& n0_vars = neighbour0.get_variables(); //variables of neighbour 0
- const vars_t& n1_vars = neighbour1.get_variables(); //variables of neighbour 1
- const vars_t& n2_vars = neighbour2.get_variables(); //variables of neighbour 2
+ 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 final override
+ {
+ //fetch data
+ //note that no copying is done here, these are all references
+ //TODO: implement drive
+ const params_t& drum_params = drum.get_parameters(); //parameters of drum
+ const vars_t& drum_vars = drum.get_variables(); //variables of drum
+ const vars_t& n0_vars = neighbour0.get_variables(); //variables of neighbour 0
+ const vars_t& n1_vars = neighbour1.get_variables(); //variables of neighbour 1
+ const vars_t& n2_vars = neighbour2.get_variables(); //variables of neighbour 2
- value_t part1 = - drum_params.k0 * drum_vars.x_temp;
- value_t part2 = - drum_params.c * drum_vars.xdot_temp;
- value_t part3 = 0.;
- part3 -= drum_params.k1 + drum_params.k2 * n0_vars.t0 * n0_vars.x_temp;
- part3 -= drum_params.k1 + drum_params.k2 * n1_vars.t1 * n1_vars.x_temp;
- part3 -= drum_params.k1 + drum_params.k2 * n2_vars.t2 * n2_vars.x_temp;
+ value_t part1 = - drum_params.k0 * drum_vars.x_temp;
+ value_t part2 = - drum_params.c * drum_vars.xdot_temp;
+ value_t part3 = 0.;
+ part3 -= drum_params.k1 + drum_params.k2 * n0_vars.t0 * n0_vars.x_temp;
+ part3 -= drum_params.k1 + drum_params.k2 * n1_vars.t1 * n1_vars.x_temp;
+ part3 -= drum_params.k1 + drum_params.k2 * n2_vars.t2 * n2_vars.x_temp;
- return part1 + part2 + part3;
- }
+ return part1 + part2 + part3;
+ }
};
diff --git a/projects/braidingTightBinding/include/grabber.hpp b/projects/braidingTightBinding/include/grabber.hpp
index 3b5e95c604ed23b8496d555fec2cc30069f7e46f..1b4af91d186208f45306d66876a10bdb63f39b05 100644
--- a/projects/braidingTightBinding/include/grabber.hpp
+++ b/projects/braidingTightBinding/include/grabber.hpp
@@ -7,134 +7,137 @@
template<typename value_t>
struct DataStorage{
- DataStorage(const size_t num_drums)
- {
- x.reserve(num_drums);
- xdot.reserve(num_drums);
- t0.reserve(num_drums);
- t1.reserve(num_drums);
- t2.reserve(num_drums);
- drive.reserve(num_drums);
- }
-
- ~DataStorage() = default;
+ DataStorage(const size_t num_drums)
+ {
+ x.reserve(num_drums);
+ xdot.reserve(num_drums);
+ t0.reserve(num_drums);
+ t1.reserve(num_drums);
+ t2.reserve(num_drums);
+ drive.reserve(num_drums);
+ }
+
+ ~DataStorage() = default;
- value_t time;
- std::vector<value_t> x;
- std::vector<value_t> xdot;
- std::vector<value_t> t0;
- std::vector<value_t> t1;
- std::vector<value_t> t2;
- std::vector<value_t> drive;
+ value_t time;
+ std::vector<value_t> x;
+ std::vector<value_t> xdot;
+ std::vector<value_t> t0;
+ std::vector<value_t> t1;
+ std::vector<value_t> t2;
+ std::vector<value_t> drive;
};
template <typename value_t, typename drum_t>
class Grabber{
- public:
- Grabber(const size_t grab_every, const std::string params_file, const std::string adjacency_file, const std::string dynamic_file)
- : grab_every_(grab_every), grab_cnt_(0), par_f_(params_file), adj_f_(adjacency_file), dyn_f_(dynamic_file) {}
- Grabber() = delete; //no default initialization
- Grabber(const Grabber&) = default;
- ~Grabber() = default;
+public:
+ Grabber(const size_t grab_every, const std::string params_file, const std::string adjacency_file, const std::string dynamic_file)
+ : grab_every_(grab_every), grab_cnt_(0), par_f_(params_file), adj_f_(adjacency_file), dyn_f_(dynamic_file) {}
+ Grabber() = delete; //no default initialization
+ Grabber(const Grabber&) = default;
+ ~Grabber() = default;
- //call in system constructor
- void init(const value_t t_end, const value_t dt, const std::vector<drum_t>& drums, const std::vector<std::vector<int> >& adjacency) noexcept
- {
- drums_ = drums;
- adjacency_ = adjacency;
- t_end_ = t_end;
- dt_ = dt;
- }
+ //call in system constructor
+ void init(const value_t t_end, const value_t dt, const std::vector<drum_t>& drums, const std::vector<std::vector<int> >& adjacency) noexcept
+ {
+ drums_ = drums;
+ adjacency_ = adjacency;
+ t_end_ = t_end;
+ dt_ = dt;
+ }
- //call after each simulation step, checks itself if it should save
- bool grab(const std::vector<drum_t>& drums, const value_t time)
- {
- if(grab_cnt_++ % grab_every_ == 0){
- data_.push_back(DataStorage<value_t>(drums_.size()));
- data_.back().time = time;
- for(auto d: drums){
- data_.back().x.push_back(d.get_variables().x);
- data_.back().xdot.push_back(d.get_variables().xdot);
- data_.back().t0.push_back(d.get_variables().t0);
- data_.back().t1.push_back(d.get_variables().t1);
- data_.back().t2.push_back(d.get_variables().t2);
- data_.back().drive.push_back(d.get_variables().V);
- }
- return true;
- }
- else
- return false;
- }
+ //call after each simulation step, checks itself if it should save
+ bool grab(const std::vector<drum_t>& drums, const value_t time)
+ {
+ if(grab_cnt_++ % grab_every_ == 0){
+ data_.push_back(DataStorage<value_t>(drums_.size()));
+ data_.back().time = time;
+ for(auto d: drums){
+ data_.back().x.push_back(d.get_variables().x);
+ data_.back().xdot.push_back(d.get_variables().xdot);
+ data_.back().t0.push_back(d.get_variables().t0);
+ data_.back().t1.push_back(d.get_variables().t1);
+ data_.back().t2.push_back(d.get_variables().t2);
+ data_.back().drive.push_back(d.get_variables().V);
+ }
+ return true;
+ }
+ else{
+ return false;
+ }
+ }
- //prints the data out. call at end of it all.
- bool save(){
- //TODO: Check for file open
- //TODO: Overload operator<< for DrumParameters and DataStorage.
- //print parameters
- std::fstream par (par_f_, par.out);
- par << "# The last drum is to be ignored\n";
- par << "# k0 \t k1 \t k2 \t c \t pos_x \t pos_y \t sublattice\n";
- for(auto d: drums_){
- par << d.get_parameters().k0 << " \t";
- par << d.get_parameters().k1 << " \t";
- par << d.get_parameters().k2 << " \t";
- par << d.get_parameters().c << " \t";
- par << d.get_parameters().position.x() << " \t";
- par << d.get_parameters().position.y() << " \t";
- par << d.get_parameters().sublattice << "\n";
- }
- par.close();
+ //prints the data out. call at end of it all.
+ bool save(){
+ //TODO: Check for file open
+ //TODO: Overload operator<< for DrumParameters and DataStorage.
+ //print parameters
+ std::fstream par (par_f_, par.out);
+ par << "# The last drum is to be ignored\n";
+ par << "# k0 \t k1 \t k2 \t c \t pos_x \t pos_y \t sublattice\n";
+ for(auto d: drums_){
+ par << d.get_parameters().k0 << " \t";
+ par << d.get_parameters().k1 << " \t";
+ par << d.get_parameters().k2 << " \t";
+ par << d.get_parameters().c << " \t";
+ par << d.get_parameters().position.x() << " \t";
+ par << d.get_parameters().position.y() << " \t";
+ par << d.get_parameters().sublattice << "\n";
+ }
+ par.close();
- //print adjacency vectors
- std::fstream adj (adj_f_, adj.out);
- adj << "# <-1> means <no neighbour here>\n";
- adj << "# The last drum is to be ignored\n";
- adj << "# n0 \t n1 \t n2\n";
- for(auto v: adjacency_){
- for(auto i: v){
- if(i == drums_.size()-1)
- adj << -1 << " \t"; //no neighbour present
- else
- adj << i << " \t"; //neighbour
- }
- adj << "\n";
- }
- adj.close();
+ //print adjacency vectors
+ std::fstream adj (adj_f_, adj.out);
+ adj << "# <-1> means <no neighbour here>\n";
+ adj << "# The last drum is to be ignored\n";
+ adj << "# n0 \t n1 \t n2\n";
+ for(auto v: adjacency_){
+ for(auto i: v){
+ if(i == drums_.size()-1){
+ adj << -1 << " \t"; //no neighbour present
+ }
+ else{
+ adj << i << " \t"; //neighbour
+ }
+ }
+ adj << "\n";
+ }
+ adj.close();
- std::fstream dyn (dyn_f_, dyn.out);
- dyn << "# If adjacency of a bond is -1 in " << adj_f_;
- dyn << ", the coupling is to be ignored\n";
- dyn << "# The last drum is to be ignored\n";
- dyn << "# Number of drums : " << drums_.size() << "\n";
- dyn << "# t \t drum_index \t x \t xdot \t t0 \t t1 \t t2 \t drive\n";
- for(auto entry: data_){
- for(size_t i = 0; i < entry.x.size(); ++i){
- dyn << entry.time << " \t";
- dyn << i << " \t";
- dyn << entry.x[i] << " \t";
- dyn << entry.xdot[i] << " \t";
- dyn << entry.t0[i] << " \t";
- dyn << entry.t1[i] << " \t";
- dyn << entry.t2[i] << " \t";
- dyn << entry.drive[i] << "\n";
- }
- }
- dyn.close();
+ std::fstream dyn (dyn_f_, dyn.out);
+ dyn << "# If adjacency of a bond is -1 in " << adj_f_;
+ dyn << ", the coupling is to be ignored\n";
+ dyn << "# The last drum is to be ignored\n";
+ dyn << "# Number of drums : " << drums_.size() << "\n";
+ dyn << "# t \t drum_index \t x \t xdot \t t0 \t t1 \t t2 \t drive\n";
+ for(auto entry: data_){
+ for(size_t i = 0; i < entry.x.size(); ++i){
+ dyn << entry.time << " \t";
+ dyn << i << " \t";
+ dyn << entry.x[i] << " \t";
+ dyn << entry.xdot[i] << " \t";
+ dyn << entry.t0[i] << " \t";
+ dyn << entry.t1[i] << " \t";
+ dyn << entry.t2[i] << " \t";
+ dyn << entry.drive[i] << "\n";
+ }
+ }
+ dyn.close();
- return true;
- }
+ return true;
+ }
- private:
- size_t grab_every_; //grab data every grab_every_ steps
- size_t grab_cnt_; //grab count
- value_t t_end_, dt_; //maximal simulation time and timestep
- std::vector<std::vector<int> > adjacency_; //adjacency vector
- std::vector<drum_t> drums_; //drums in initial configuration, grab params from here
- std::list<DataStorage<value_t> > data_; //here goes the dynamical data
- //filenames
- std::string par_f_; //parameters file
- std::string adj_f_; //adjacency file
- std::string dyn_f_; //dynamics file (this is the gold)
+private:
+ size_t grab_every_; //grab data every grab_every_ steps
+ size_t grab_cnt_; //grab count
+ value_t t_end_, dt_; //maximal simulation time and timestep
+ std::vector<std::vector<int> > adjacency_; //adjacency vector
+ std::vector<drum_t> drums_; //drums in initial configuration, grab params from here
+ std::list<DataStorage<value_t> > data_; //here goes the dynamical data
+ //filenames
+ std::string par_f_; //parameters file
+ std::string adj_f_; //adjacency file
+ std::string dyn_f_; //dynamics file (this is the gold)
};
#endif
diff --git a/projects/braidingTightBinding/include/matrix_element_calculator_braid.hpp b/projects/braidingTightBinding/include/matrix_element_calculator_braid.hpp
index e09455181531bc4510ea629624333421502f8c20..3060c0c572588ae7ba1e920f2921d8f64c10ac92 100644
--- a/projects/braidingTightBinding/include/matrix_element_calculator_braid.hpp
+++ b/projects/braidingTightBinding/include/matrix_element_calculator_braid.hpp
@@ -4,32 +4,32 @@
template <typename value_t, typename params_t, typename vars_t, typename drum_t>
class MatrixElementCalculatorBraid: public MatrixElementCalculator<value_t, params_t, vars_t, drum_t>{
- public:
- MatrixElementCalculatorBraid() = default;
- ~MatrixElementCalculatorBraid() = default;
+public:
+ MatrixElementCalculatorBraid() = default;
+ ~MatrixElementCalculatorBraid() = default;
- //diagonal element at (index, index)
- value_t operator()(const size_t index, const std::vector<drum_t>& drums) const noexcept final override
- {
- return drums[index].get_parameters().k0;
- }
+ //diagonal element at (index, index)
+ value_t operator()(const size_t index, const std::vector<drum_t>& drums) const noexcept final override
+ {
+ return drums[index].get_parameters().k0;
+ }
- //coupling elements at (index1, index2)
- //for neighbour 0
- value_t operator()(const size_t index1, const size_t index2, const std::vector<drum_t>& drums) const noexcept final override
- {
- return drums[index1].get_parameters().k1 + drums[index1].get_parameters().k2 * drums[index2].get_variables().t0;
- }
- //for neighbour 1
- value_t operator()(const size_t index1, const size_t index2, const std::vector<drum_t>& drums, const int) const noexcept final override
- {
- return drums[index1].get_parameters().k1 + drums[index1].get_parameters().k2 * drums[index2].get_variables().t1;
- }
- //for neighbour 2
- value_t operator()(const size_t index1, const size_t index2, const std::vector<drum_t>& drums, const int, const int) const noexcept final override
- {
- return drums[index1].get_parameters().k1 + drums[index1].get_parameters().k2 * drums[index2].get_variables().t2;
- }
+ //coupling elements at (index1, index2)
+ //for neighbour 0
+ value_t operator()(const size_t index1, const size_t index2, const std::vector<drum_t>& drums) const noexcept final override
+ {
+ return drums[index1].get_parameters().k1 + drums[index1].get_parameters().k2 * drums[index2].get_variables().t0;
+ }
+ //for neighbour 1
+ value_t operator()(const size_t index1, const size_t index2, const std::vector<drum_t>& drums, const int) const noexcept final override
+ {
+ return drums[index1].get_parameters().k1 + drums[index1].get_parameters().k2 * drums[index2].get_variables().t1;
+ }
+ //for neighbour 2
+ value_t operator()(const size_t index1, const size_t index2, const std::vector<drum_t>& drums, const int, const int) const noexcept final override
+ {
+ return drums[index1].get_parameters().k1 + drums[index1].get_parameters().k2 * drums[index2].get_variables().t2;
+ }
};
#endif
diff --git a/projects/braidingTightBinding/include/rbcomb_generator_braid.hpp b/projects/braidingTightBinding/include/rbcomb_generator_braid.hpp
index 3420ea29516eb754e1fd2aa51a02dfa4b2d94be7..1b0c33192e8237bddd8dfce01072c4eb80cdb5d4 100644
--- a/projects/braidingTightBinding/include/rbcomb_generator_braid.hpp
+++ b/projects/braidingTightBinding/include/rbcomb_generator_braid.hpp
@@ -7,381 +7,394 @@
template <typename value_t, typename params_t, typename vars_t, typename sbuffer_t>
class RbcombGeneratorBraid: public LatticeGenerator<value_t, params_t, vars_t, sbuffer_t> {
- public:
- RbcombGeneratorBraid(size_t x_extent, size_t y_extent) noexcept : x_extent_(x_extent), y_extent_(y_extent) {}
- RbcombGeneratorBraid() = delete;
- ~RbcombGeneratorBraid() = default;
-
- //the params_t argument shall contain the position of the first drum that is placed.
- //the last drum (return.first.back()) signifies inexistent neighbours, it is not part of the lattice.
- //It shall be of sublattice 'A'.
- std::pair<std::vector<Drum<value_t, params_t, vars_t, sbuffer_t> >, std::vector<std::vector<int> > >
- operator()(const params_t& drum_params) noexcept final override
- {
- using drum_t = Drum<value_t, params_t, vars_t, sbuffer_t>;
- std::vector<drum_t> current_line, next_line;
- current_line.reserve(x_extent_);
- next_line.reserve(x_extent_);
-
- //calculate vectors to next drums
- Vec2<value_t> vec_right = drum_params.s2 - drum_params.s1; //drum to the right in same row
- Vec2<value_t> vec_down = drum_params.s0; //drum below
- Vec2<value_t> vec_upleft = drum_params.s1; //drum at top left
- Vec2<value_t> vec_upright = drum_params.s2; //drum at top right
-
- //prepare memory
- drums_.reserve(x_extent_);
- adjacency_vector_.reserve(x_extent_);
- //push first drum
- drums_.push_back(drum_t(drum_params));
- adjacency_vector_.push_back(std::vector<int>(3,-1)); //initialized with 'no neighbour here'
-
- //initialize the first line to bootstrap the algorithm
- for(size_t i = 1; i < x_extent_; ++i){
- params_t new_params (drum_params);
- new_params.position = drums_.back().get_parameters().position + vec_right;
- drums_.push_back(drum_t(new_params));
- adjacency_vector_.push_back(std::vector<int>(3,-1));
- }
-
- //add y_extent_ inequivalent hexagon rows
- for(size_t i = 0; i < y_extent_; ++i){
- generate_layer(vec_down, vec_upleft, vec_upright);
- }
-
- //terminate by adding another B-row
- terminate_lattice(vec_upleft, vec_upright);
-
- //push a drum to the end that has no neighbours and signifies the inexistent neighbour
- drums_.push_back(drum_t(drum_params));
- adjacency_vector_.push_back(std::vector<int>(3,-1));
- //make no-neighbours point here
- for(size_t i = 0; i < adjacency_vector_.size()-1; ++i){
- for(size_t j = 0; j < 3; ++j){
- if(adjacency_vector_[i][j] == -1){
- adjacency_vector_[i][j] = drums_.size()-1;
- }
- }
- }
-
- /* Deploy without consistency check
- if(check_consistency() == false)
- std::cout << "INCONSISTENT LATTICE DETECTED!\n";
- */
-
- return std::pair<std::vector<drum_t>, std::vector<std::vector<int> > > (drums_, adjacency_vector_);
- }
-
- private:
- //check if there exists a drum in drum_vec that is close to drum_pos.
- //returns the index of the first close drum or -1.
- int find_close(
- const Vec2<value_t>& drum_pos,
- const std::vector<Drum<value_t, params_t, vars_t, sbuffer_t> >& drum_vec) const noexcept
- {
- for(int i = 0; i < drum_vec.size(); ++i){
- auto test_drum = drum_vec[i];
- if(test_drum.get_parameters().position.r_wrt(drum_pos) < 1.e-8)
- return i;
- }
- return -1; //no close drum found
- }
-
- //checks sanity of generated lattice
- //returns true if lattice passed the test
- bool check_consistency() noexcept
- {
- //check for duplicate drums
- for(size_t i = 0; i < drums_.size(); ++i){
- if(i != find_close(drums_[i].get_parameters().position, drums_))
- return false;
- }
- return true;
- }
-
- //arguments are the neighbour vectors down, upleft, upright.
- //adds a new layer of drums to drums_, filling in the adjency_vector_ as well.
- //
- // 4: | x_extent_ drums
- // 3: \/ x_extent_ + 1 drums
- // 2: | x_extent_ + 1 drums
- // 1: \/ x_extent_ drums
- void generate_layer(const Vec2<value_t>& v_d, const Vec2<value_t>& v_ul, const Vec2<value_t>& v_ur) noexcept
- {
- //shorten syntax
- using drum_t = Drum<value_t, params_t, vars_t, sbuffer_t>;
-
- //prepare vector to hold new row
- std::vector<drum_t> new_row;
- new_row.reserve(x_extent_);
-
- //prepare parameter templates for both sublattices
- params_t params_A(drums_[0].get_parameters()); //first drum will always be of sublattice 'A'
- params_t params_B(params_A.k0, params_A.k1, params_A.k2, params_A.c, params_A.position, 'B'); //generate sublattice 'B' type of parameters
-
- //row 1
- //prepare iterator
- auto drum_it = drums_.end();
- drum_it -= static_cast<int>(x_extent_);
-
- //first drum is special
- Vec2<value_t> new_pos_first = drum_it->get_parameters().position + v_ur;
- params_t new_params_first (params_B);
- new_params_first.position = new_pos_first;
- new_row.push_back(drum_t(new_params_first));
- int index_1 = drums_.size() + new_row.size()-1;
- adjacency_vector_.push_back(std::vector<int>(3,-1));
- adjacency_vector_[std::distance(drums_.begin(), drum_it)][2] = drums_.size();
- adjacency_vector_[drums_.size()][2] = std::distance(drums_.begin(), drum_it);
- ++drum_it;
-
- //now all other drums
- while(drum_it != drums_.end()){
- Vec2<value_t> new_pos1 = drum_it->get_parameters().position + v_ul;
- Vec2<value_t> new_pos2 = drum_it->get_parameters().position + v_ur;
-
- //see if these drums already exist and add them to the list
- int index_1 = find_close(new_pos1, new_row);
- if(index_1 == -1){ //not in there yet, so push it
- params_t new_params (params_B);
- new_params.position = new_pos1;
- new_row.push_back(drum_t(new_params));
- index_1 = drums_.size()+new_row.size()-1;
- //make adjacency vector longer
- adjacency_vector_.push_back(std::vector<int>(3,-1));
- }
- else
- index_1 += drums_.size();
- int index_2 = find_close(new_pos2, new_row);
- if(index_2 == -1){ //not in there yet, so push it
- params_t new_params (params_B);
- new_params.position = new_pos2;
- new_row.push_back(drum_t(new_params));
- index_2 = drums_.size()+new_row.size()-1;
- //make adjacency vector longer
- adjacency_vector_.push_back(std::vector<int>(3,-1));
- }
- else
- index_2 += drums_.size();
-
- //add new bonds to the adjacency vector
- int old_index = std::distance(drums_.begin(), drum_it);
- adjacency_vector_[old_index][1] = index_1;
- adjacency_vector_[old_index][2] = index_2;
- adjacency_vector_[index_1][1] = old_index;
- adjacency_vector_[index_2][2] = old_index;
-
- ++drum_it;
- }
-
- //push new drums
- for(auto d: new_row)
- drums_.push_back(d);
-
- //clear new_row
- new_row.clear();
-
- //row 2
- //prepare iterator
- drum_it = drums_.end();
- drum_it -= static_cast<int>(x_extent_);
-
- //first drum is special, we add one to the left
- new_params_first = params_A;
- new_params_first.position = drum_it->get_parameters().position - v_d + v_ul - v_ur;
- new_row.push_back(drum_t(new_params_first)); //this drum has no neighbours yet.
- adjacency_vector_.push_back(std::vector<int>(3,-1));
- //we don't increment the iterator, as we reuse this drum for straight up addition
-
- while(drum_it != drums_.end()-1){
- //These drums are guaranteed to be unique, hence inexistent.
- params_t new_params (params_A);
- new_params.position = drum_it->get_parameters().position - v_d;
- new_row.push_back(drum_t(new_params));
- int index = drums_.size() + new_row.size() - 1;
- adjacency_vector_.push_back(std::vector<int>(3,-1));
-
- int old_index = std::distance(drums_.begin(), drum_it);
- adjacency_vector_[old_index][0] = index;
- adjacency_vector_[index][0] = old_index;
-
- ++drum_it;
- }
-
- //push new drums
- for(auto d: new_row)
- drums_.push_back(d);
-
- //clear new_row
- new_row.clear();
-
- //row 3
- drum_it = drums_.end();
- drum_it -= static_cast<int>(x_extent_); //this row has one drum more
-
- //first drum is special (only goes to ur)
- params_t new_params (params_B);
- new_params.position = drum_it->get_parameters().position + v_ur;
- new_row.push_back(drum_t(new_params));
- int indexx = drums_.size();
- adjacency_vector_.push_back(std::vector<int>(3,-1));
- int old_indexx = std::distance(drums_.begin(), drum_it); //avoid naming conflicts
- adjacency_vector_[old_indexx][2] = indexx;
- adjacency_vector_[indexx][2] = old_indexx;
- ++drum_it;
-
- //now the unspecial ones
- while(drum_it != drums_.end()){
- Vec2<value_t> new_pos1 = drum_it->get_parameters().position + v_ul;
- Vec2<value_t> new_pos2 = drum_it->get_parameters().position + v_ur;
-
- //see if these drums already exist and add them to the list
- int index_1 = find_close(new_pos1, new_row);
- if(index_1 == -1){ //not in there yet, so push it
- params_t new_params (params_B);
- new_params.position = new_pos1;
- new_row.push_back(drum_t(new_params));
- index_1 = drums_.size()+new_row.size()-1;
- //make adjacency vector longer
- adjacency_vector_.push_back(std::vector<int>(3,-1));
- }
- else
- index_1 += drums_.size();
- int index_2 = find_close(new_pos2, new_row);
- if(index_2 == -1){ //not in there yet, so push it
- params_t new_params (params_B);
- new_params.position = new_pos2;
- new_row.push_back(drum_t(new_params));
- index_2 = drums_.size()+new_row.size()-1;
- //make adjacency vector longer
- adjacency_vector_.push_back(std::vector<int>(3,-1));
- }
- else
- index_2 += drums_.size();
-
- //add new bonds to the adjacency vector
- int old_index = std::distance(drums_.begin(), drum_it);
- adjacency_vector_[old_index][1] = index_1;
- adjacency_vector_[old_index][2] = index_2;
- adjacency_vector_[index_1][1] = old_index;
- adjacency_vector_[index_2][2] = old_index;
-
- ++drum_it;
- }
-
- //push new drums
- for(auto d: new_row)
- drums_.push_back(d);
-
- //clear new_row
- new_row.clear();
-
-
-
- //row 4
- //prepare iterator
- drum_it = drums_.end();
- drum_it -= static_cast<int>(x_extent_);
-
- while(drum_it != drums_.end()){
- //These drums are guaranteed to be unique, hence inexistent.
- params_t new_params (params_A);
- new_params.position = drum_it->get_parameters().position - v_d;
- new_row.push_back(drum_t(new_params));
- int index = drums_.size() + new_row.size() - 1;
- adjacency_vector_.push_back(std::vector<int>(3,-1));
-
- int old_index = std::distance(drums_.begin(), drum_it);
- adjacency_vector_[old_index][0] = index;
- adjacency_vector_[index][0] = old_index;
-
- ++drum_it;
- }
-
- //push new drums
- for(auto d: new_row)
- drums_.push_back(d);
-
- //clear new_row
- new_row.clear();
- }
-
- void terminate_lattice(const Vec2<value_t>& v_ul, const Vec2<value_t>& v_ur) noexcept {
- //shorten syntax
- using drum_t = Drum<value_t, params_t, vars_t, sbuffer_t>;
-
- //prepare vector to hold new row
- std::vector<drum_t> new_row;
- new_row.reserve(x_extent_);
-
- //prepare parameter templates for both sublattices
- params_t params_A(drums_[0].get_parameters()); //first drum will always be of sublattice 'A'
- params_t params_B(params_A.k0, params_A.k1, params_A.k2, params_A.c, params_A.position, 'B'); //generate sublattice 'B' type of parameters
-
- //row 1
- //prepare iterator
- auto drum_it = drums_.end();
- drum_it -= static_cast<int>(x_extent_);
-
- //first drum is special
- Vec2<value_t> new_pos_first = drum_it->get_parameters().position + v_ur;
- params_t new_params_first (params_B);
- new_params_first.position = new_pos_first;
- new_row.push_back(drum_t(new_params_first));
- int index_1 = drums_.size() + new_row.size()-1;
- adjacency_vector_.push_back(std::vector<int>(3,-1));
- adjacency_vector_[std::distance(drums_.begin(), drum_it)][2] = drums_.size();
- adjacency_vector_[drums_.size()][2] = std::distance(drums_.begin(), drum_it);
- ++drum_it;
-
- //now all other drums
- while(drum_it != drums_.end()){
- Vec2<value_t> new_pos1 = drum_it->get_parameters().position + v_ul;
- Vec2<value_t> new_pos2 = drum_it->get_parameters().position + v_ur;
-
- //see if these drums already exist and add them to the list
- int index_1 = find_close(new_pos1, new_row);
- if(index_1 == -1){ //not in there yet, so push it
- params_t new_params (params_B);
- new_params.position = new_pos1;
- new_row.push_back(drum_t(new_params));
- index_1 = drums_.size()+new_row.size()-1;
- //make adjacency vector longer
- adjacency_vector_.push_back(std::vector<int>(3,-1));
- }
- else
- index_1 += drums_.size();
- int index_2 = find_close(new_pos2, new_row);
- if(index_2 == -1){ //not in there yet, so push it
- params_t new_params (params_B);
- new_params.position = new_pos2;
- new_row.push_back(drum_t(new_params));
- index_2 = drums_.size()+new_row.size()-1;
- //make adjacency vector longer
- adjacency_vector_.push_back(std::vector<int>(3,-1));
- }
- else
- index_2 += drums_.size();
-
- //add new bonds to the adjacency vector
- int old_index = std::distance(drums_.begin(), drum_it);
- adjacency_vector_[old_index][1] = index_1;
- adjacency_vector_[old_index][2] = index_2;
- adjacency_vector_[index_1][1] = old_index;
- adjacency_vector_[index_2][2] = old_index;
-
- ++drum_it;
- }
-
- //push new drums
- for(auto d: new_row)
- drums_.push_back(d);
-
- }
-
- size_t x_extent_, y_extent_; //Number of equivalent hexagon rows in x and y direction
- std::vector<Drum<value_t, params_t, vars_t, sbuffer_t> > drums_;
- std::vector<std::vector<int> > adjacency_vector_;
+public:
+ RbcombGeneratorBraid(size_t x_extent, size_t y_extent) noexcept : x_extent_(x_extent), y_extent_(y_extent) {}
+ RbcombGeneratorBraid() = delete;
+ ~RbcombGeneratorBraid() = default;
+
+ //the params_t argument shall contain the position of the first drum that is placed.
+ //the last drum (return.first.back()) signifies inexistent neighbours, it is not part of the lattice.
+ //It shall be of sublattice 'A'.
+ std::pair<std::vector<Drum<value_t, params_t, vars_t, sbuffer_t> >, std::vector<std::vector<int> > >
+ operator()(const params_t& drum_params) noexcept final override
+ {
+ using drum_t = Drum<value_t, params_t, vars_t, sbuffer_t>;
+ std::vector<drum_t> current_line, next_line;
+ current_line.reserve(x_extent_);
+ next_line.reserve(x_extent_);
+
+ //calculate vectors to next drums
+ Vec2<value_t> vec_right = drum_params.s2 - drum_params.s1; //drum to the right in same row
+ Vec2<value_t> vec_down = drum_params.s0; //drum below
+ Vec2<value_t> vec_upleft = drum_params.s1; //drum at top left
+ Vec2<value_t> vec_upright = drum_params.s2; //drum at top right
+
+ //prepare memory
+ drums_.reserve(x_extent_);
+ adjacency_vector_.reserve(x_extent_);
+ //push first drum
+ drums_.push_back(drum_t(drum_params));
+ adjacency_vector_.push_back(std::vector<int>(3,-1)); //initialized with 'no neighbour here'
+
+ //initialize the first line to bootstrap the algorithm
+ for(size_t i = 1; i < x_extent_; ++i){
+ params_t new_params (drum_params);
+ new_params.position = drums_.back().get_parameters().position + vec_right;
+ drums_.push_back(drum_t(new_params));
+ adjacency_vector_.push_back(std::vector<int>(3,-1));
+ }
+
+ //add y_extent_ inequivalent hexagon rows
+ for(size_t i = 0; i < y_extent_; ++i){
+ generate_layer(vec_down, vec_upleft, vec_upright);
+ }
+
+ //terminate by adding another B-row
+ terminate_lattice(vec_upleft, vec_upright);
+
+ //push a drum to the end that has no neighbours and signifies the inexistent neighbour
+ drums_.push_back(drum_t(drum_params));
+ adjacency_vector_.push_back(std::vector<int>(3,-1));
+ //make no-neighbours point here
+ for(size_t i = 0; i < adjacency_vector_.size()-1; ++i){
+ for(size_t j = 0; j < 3; ++j){
+ if(adjacency_vector_[i][j] == -1){
+ adjacency_vector_[i][j] = drums_.size()-1;
+ }
+ }
+ }
+
+ /* Deploy without consistency check
+ if(check_consistency() == false)
+ std::cout << "INCONSISTENT LATTICE DETECTED!\n";
+ */
+
+ return std::pair<std::vector<drum_t>, std::vector<std::vector<int> > > (drums_, adjacency_vector_);
+ }
+
+private:
+ //check if there exists a drum in drum_vec that is close to drum_pos.
+ //returns the index of the first close drum or -1.
+ int find_close(
+ const Vec2<value_t>& drum_pos,
+ const std::vector<Drum<value_t, params_t, vars_t, sbuffer_t> >& drum_vec
+ ) const noexcept
+ {
+ for(int i = 0; i < drum_vec.size(); ++i){
+ auto test_drum = drum_vec[i];
+ if(test_drum.get_parameters().position.r_wrt(drum_pos) < 1.e-8){
+ return i;
+ }
+ }
+ return -1; //no close drum found
+ }
+
+ //checks sanity of generated lattice
+ //returns true if lattice passed the test
+ bool check_consistency() noexcept
+ {
+ //check for duplicate drums
+ for(size_t i = 0; i < drums_.size(); ++i){
+ if(i != find_close(drums_[i].get_parameters().position, drums_)){
+ return false;
+ }
+ }
+ return true;
+ }
+
+ //arguments are the neighbour vectors down, upleft, upright.
+ //adds a new layer of drums to drums_, filling in the adjency_vector_ as well.
+ //
+ // 4: | x_extent_ drums
+ // 3: \/ x_extent_ + 1 drums
+ // 2: | x_extent_ + 1 drums
+ // 1: \/ x_extent_ drums
+ void generate_layer(const Vec2<value_t>& v_d, const Vec2<value_t>& v_ul, const Vec2<value_t>& v_ur) noexcept
+ {
+ //shorten syntax
+ using drum_t = Drum<value_t, params_t, vars_t, sbuffer_t>;
+
+ //prepare vector to hold new row
+ std::vector<drum_t> new_row;
+ new_row.reserve(x_extent_);
+
+ //prepare parameter templates for both sublattices
+ params_t params_A(drums_[0].get_parameters()); //first drum will always be of sublattice 'A'
+ params_t params_B(params_A.k0, params_A.k1, params_A.k2, params_A.c, params_A.position, 'B'); //generate sublattice 'B' type of parameters
+
+ //row 1
+ //prepare iterator
+ auto drum_it = drums_.end();
+ drum_it -= static_cast<int>(x_extent_);
+
+ //first drum is special
+ Vec2<value_t> new_pos_first = drum_it->get_parameters().position + v_ur;
+ params_t new_params_first (params_B);
+ new_params_first.position = new_pos_first;
+ new_row.push_back(drum_t(new_params_first));
+ int index_1 = drums_.size() + new_row.size()-1;
+ adjacency_vector_.push_back(std::vector<int>(3,-1));
+ adjacency_vector_[std::distance(drums_.begin(), drum_it)][2] = drums_.size();
+ adjacency_vector_[drums_.size()][2] = std::distance(drums_.begin(), drum_it);
+ ++drum_it;
+
+ //now all other drums
+ while(drum_it != drums_.end()){
+ Vec2<value_t> new_pos1 = drum_it->get_parameters().position + v_ul;
+ Vec2<value_t> new_pos2 = drum_it->get_parameters().position + v_ur;
+
+ //see if these drums already exist and add them to the list
+ int index_1 = find_close(new_pos1, new_row);
+ if(index_1 == -1){ //not in there yet, so push it
+ params_t new_params (params_B);
+ new_params.position = new_pos1;
+ new_row.push_back(drum_t(new_params));
+ index_1 = drums_.size()+new_row.size()-1;
+ //make adjacency vector longer
+ adjacency_vector_.push_back(std::vector<int>(3,-1));
+ }
+ else{
+ index_1 += drums_.size();
+ }
+ int index_2 = find_close(new_pos2, new_row);
+ if(index_2 == -1){ //not in there yet, so push it
+ params_t new_params (params_B);
+ new_params.position = new_pos2;
+ new_row.push_back(drum_t(new_params));
+ index_2 = drums_.size()+new_row.size()-1;
+ //make adjacency vector longer
+ adjacency_vector_.push_back(std::vector<int>(3,-1));
+ }
+ else{
+ index_2 += drums_.size();
+ }
+ //add new bonds to the adjacency vector
+ int old_index = std::distance(drums_.begin(), drum_it);
+ adjacency_vector_[old_index][1] = index_1;
+ adjacency_vector_[old_index][2] = index_2;
+ adjacency_vector_[index_1][1] = old_index;
+ adjacency_vector_[index_2][2] = old_index;
+
+ ++drum_it;
+ }
+
+ //push new drums
+ for(auto d: new_row){
+ drums_.push_back(d);
+ }
+
+ //clear new_row
+ new_row.clear();
+
+ //row 2
+ //prepare iterator
+ drum_it = drums_.end();
+ drum_it -= static_cast<int>(x_extent_);
+
+ //first drum is special, we add one to the left
+ new_params_first = params_A;
+ new_params_first.position = drum_it->get_parameters().position - v_d + v_ul - v_ur;
+ new_row.push_back(drum_t(new_params_first)); //this drum has no neighbours yet.
+ adjacency_vector_.push_back(std::vector<int>(3,-1));
+ //we don't increment the iterator, as we reuse this drum for straight up addition
+
+ while(drum_it != drums_.end()-1){
+ //These drums are guaranteed to be unique, hence inexistent.
+ params_t new_params (params_A);
+ new_params.position = drum_it->get_parameters().position - v_d;
+ new_row.push_back(drum_t(new_params));
+ int index = drums_.size() + new_row.size() - 1;
+ adjacency_vector_.push_back(std::vector<int>(3,-1));
+
+ int old_index = std::distance(drums_.begin(), drum_it);
+ adjacency_vector_[old_index][0] = index;
+ adjacency_vector_[index][0] = old_index;
+
+ ++drum_it;
+ }
+
+ //push new drums
+ for(auto d: new_row){
+ drums_.push_back(d);
+ }
+
+ //clear new_row
+ new_row.clear();
+
+ //row 3
+ drum_it = drums_.end();
+ drum_it -= static_cast<int>(x_extent_); //this row has one drum more
+
+ //first drum is special (only goes to ur)
+ params_t new_params (params_B);
+ new_params.position = drum_it->get_parameters().position + v_ur;
+ new_row.push_back(drum_t(new_params));
+ int indexx = drums_.size();
+ adjacency_vector_.push_back(std::vector<int>(3,-1));
+ int old_indexx = std::distance(drums_.begin(), drum_it); //avoid naming conflicts
+ adjacency_vector_[old_indexx][2] = indexx;
+ adjacency_vector_[indexx][2] = old_indexx;
+ ++drum_it;
+
+ //now the unspecial ones
+ while(drum_it != drums_.end()){
+ Vec2<value_t> new_pos1 = drum_it->get_parameters().position + v_ul;
+ Vec2<value_t> new_pos2 = drum_it->get_parameters().position + v_ur;
+
+ //see if these drums already exist and add them to the list
+ int index_1 = find_close(new_pos1, new_row);
+ if(index_1 == -1){ //not in there yet, so push it
+ params_t new_params (params_B);
+ new_params.position = new_pos1;
+ new_row.push_back(drum_t(new_params));
+ index_1 = drums_.size()+new_row.size()-1;
+ //make adjacency vector longer
+ adjacency_vector_.push_back(std::vector<int>(3,-1));
+ }
+ else{
+ index_1 += drums_.size();
+ }
+ int index_2 = find_close(new_pos2, new_row);
+ if(index_2 == -1){ //not in there yet, so push it
+ params_t new_params (params_B);
+ new_params.position = new_pos2;
+ new_row.push_back(drum_t(new_params));
+ index_2 = drums_.size()+new_row.size()-1;
+ //make adjacency vector longer
+ adjacency_vector_.push_back(std::vector<int>(3,-1));
+ }
+ else{
+ index_2 += drums_.size();
+ }
+
+ //add new bonds to the adjacency vector
+ int old_index = std::distance(drums_.begin(), drum_it);
+ adjacency_vector_[old_index][1] = index_1;
+ adjacency_vector_[old_index][2] = index_2;
+ adjacency_vector_[index_1][1] = old_index;
+ adjacency_vector_[index_2][2] = old_index;
+
+ ++drum_it;
+ }
+
+ //push new drums
+ for(auto d: new_row){
+ drums_.push_back(d);
+ }
+
+ //clear new_row
+ new_row.clear();
+
+
+
+ //row 4
+ //prepare iterator
+ drum_it = drums_.end();
+ drum_it -= static_cast<int>(x_extent_);
+
+ while(drum_it != drums_.end()){
+ //These drums are guaranteed to be unique, hence inexistent.
+ params_t new_params (params_A);
+ new_params.position = drum_it->get_parameters().position - v_d;
+ new_row.push_back(drum_t(new_params));
+ int index = drums_.size() + new_row.size() - 1;
+ adjacency_vector_.push_back(std::vector<int>(3,-1));
+
+ int old_index = std::distance(drums_.begin(), drum_it);
+ adjacency_vector_[old_index][0] = index;
+ adjacency_vector_[index][0] = old_index;
+
+ ++drum_it;
+ }
+
+ //push new drums
+ for(auto d: new_row){
+ drums_.push_back(d);
+ }
+
+ //clear new_row
+ new_row.clear();
+ }
+
+ void terminate_lattice(const Vec2<value_t>& v_ul, const Vec2<value_t>& v_ur) noexcept {
+ //shorten syntax
+ using drum_t = Drum<value_t, params_t, vars_t, sbuffer_t>;
+
+ //prepare vector to hold new row
+ std::vector<drum_t> new_row;
+ new_row.reserve(x_extent_);
+
+ //prepare parameter templates for both sublattices
+ params_t params_A(drums_[0].get_parameters()); //first drum will always be of sublattice 'A'
+ params_t params_B(params_A.k0, params_A.k1, params_A.k2, params_A.c, params_A.position, 'B'); //generate sublattice 'B' type of parameters
+
+ //row 1
+ //prepare iterator
+ auto drum_it = drums_.end();
+ drum_it -= static_cast<int>(x_extent_);
+
+ //first drum is special
+ Vec2<value_t> new_pos_first = drum_it->get_parameters().position + v_ur;
+ params_t new_params_first (params_B);
+ new_params_first.position = new_pos_first;
+ new_row.push_back(drum_t(new_params_first));
+ int index_1 = drums_.size() + new_row.size()-1;
+ adjacency_vector_.push_back(std::vector<int>(3,-1));
+ adjacency_vector_[std::distance(drums_.begin(), drum_it)][2] = drums_.size();
+ adjacency_vector_[drums_.size()][2] = std::distance(drums_.begin(), drum_it);
+ ++drum_it;
+
+ //now all other drums
+ while(drum_it != drums_.end()){
+ Vec2<value_t> new_pos1 = drum_it->get_parameters().position + v_ul;
+ Vec2<value_t> new_pos2 = drum_it->get_parameters().position + v_ur;
+
+ //see if these drums already exist and add them to the list
+ int index_1 = find_close(new_pos1, new_row);
+ if(index_1 == -1){ //not in there yet, so push it
+ params_t new_params (params_B);
+ new_params.position = new_pos1;
+ new_row.push_back(drum_t(new_params));
+ index_1 = drums_.size()+new_row.size()-1;
+ //make adjacency vector longer
+ adjacency_vector_.push_back(std::vector<int>(3,-1));
+ }
+ else{
+ index_1 += drums_.size();
+ }
+ int index_2 = find_close(new_pos2, new_row);
+ if(index_2 == -1){ //not in there yet, so push it
+ params_t new_params (params_B);
+ new_params.position = new_pos2;
+ new_row.push_back(drum_t(new_params));
+ index_2 = drums_.size()+new_row.size()-1;
+ //make adjacency vector longer
+ adjacency_vector_.push_back(std::vector<int>(3,-1));
+ }
+ else{
+ index_2 += drums_.size();
+ }
+
+ //add new bonds to the adjacency vector
+ int old_index = std::distance(drums_.begin(), drum_it);
+ adjacency_vector_[old_index][1] = index_1;
+ adjacency_vector_[old_index][2] = index_2;
+ adjacency_vector_[index_1][1] = old_index;
+ adjacency_vector_[index_2][2] = old_index;
+
+ ++drum_it;
+ }
+
+ //push new drums
+ for(auto d: new_row){
+ drums_.push_back(d);
+ }
+
+ }
+
+ size_t x_extent_, y_extent_; //Number of equivalent hexagon rows in x and y direction
+ std::vector<Drum<value_t, params_t, vars_t, sbuffer_t> > drums_;
+ std::vector<std::vector<int> > adjacency_vector_;
};
#endif
diff --git a/projects/braidingTightBinding/include/vortex.hpp b/projects/braidingTightBinding/include/vortex.hpp
index eb391e17fc69cc337945712c8684f45950aeace7..41233ba70dfb4eebf007cf2dde7c451255bdb4fe 100644
--- a/projects/braidingTightBinding/include/vortex.hpp
+++ b/projects/braidingTightBinding/include/vortex.hpp
@@ -6,43 +6,43 @@
template <typename value_t>
class Vortex{
- public:
- //constructors
- Vortex(const Vec2<value_t> position, const value_t l0, const value_t alpha, const value_t delta, const value_t a) noexcept: position_(position), l0_(l0), alpha_(alpha), delta_(delta), K_(value_t(4.*3.141592653589793/(3.*std::sqrt(3)*a)), value_t(0.)) {}
- Vortex(const Vec2<value_t> position, const Vortex& o) noexcept: position_(position), l0_(o.l0_), alpha_(o.alpha_), delta_(o.delta_), K_(o.K_) {}
- Vortex() = default;
- Vortex(const Vortex&) = default;
- Vortex& operator=(const Vortex&) = default;
- ~Vortex() = default;
+public:
+ //constructors
+ Vortex(const Vec2<value_t> position, const value_t l0, const value_t alpha, const value_t delta, const value_t a) noexcept: position_(position), l0_(l0), alpha_(alpha), delta_(delta), K_(value_t(4.*3.141592653589793/(3.*std::sqrt(3)*a)), value_t(0.)) {}
+ Vortex(const Vec2<value_t> position, const Vortex& o) noexcept: position_(position), l0_(o.l0_), alpha_(o.alpha_), delta_(o.delta_), K_(o.K_) {}
+ Vortex() = default;
+ Vortex(const Vortex&) = default;
+ Vortex& operator=(const Vortex&) = default;
+ ~Vortex() = default;
- //access
- value_t l0() const noexcept { return l0_; }
- value_t alpha() const noexcept { return alpha_; }
- value_t delta() const noexcept { return delta_; }
- Vec2<value_t>& position() noexcept { return position_; }
- const Vec2<value_t>& position() const noexcept { return position_; }
+ //access
+ value_t l0() const noexcept { return l0_; }
+ value_t alpha() const noexcept { return alpha_; }
+ value_t delta() const noexcept { return delta_; }
+ Vec2<value_t>& position() noexcept { return position_; }
+ const Vec2<value_t>& position() const noexcept { return position_; }
- //modifiers
- void set_l0(value_t l0) noexcept { l0_ = l0; }
- void set_alpha(value_t alpha) noexcept { alpha_ = alpha; }
- void set_delta(value_t delta) noexcept { delta_ = delta; }
- void set_position(const Vec2<value_t>& position) noexcept { position_ = position; }
- void move_by(const Vec2<value_t>& translation) noexcept { position_ += translation; }
- void move_to(const Vec2<value_t>& position) noexcept { position_ = position; }
+ //modifiers
+ void set_l0(value_t l0) noexcept { l0_ = l0; }
+ void set_alpha(value_t alpha) noexcept { alpha_ = alpha; }
+ void set_delta(value_t delta) noexcept { delta_ = delta; }
+ void set_position(const Vec2<value_t>& position) noexcept { position_ = position; }
+ void move_by(const Vec2<value_t>& translation) noexcept { position_ += translation; }
+ void move_to(const Vec2<value_t>& position) noexcept { position_ = position; }
- //functional
- //there is no divide_by_zero check for efficiency reasons. l0_ can not be zero upon call.
- std::complex<value_t> distortion(const Vec2<value_t>& v) const
- {
- return std::polar<value_t>(delta_*std::tanh(v.r_wrt(position_)/l0_), alpha_ - v.phi_wrt(position_));
- }
- std::complex<value_t> kekule(const Vec2<value_t>& v, const Vec2<value_t>& s) const{
- return std::polar<value_t>(1., K_*(s + 2.*(v - position_)));
- }
- private:
- Vec2<value_t> position_;
- const Vec2<value_t> K_;
- value_t l0_, alpha_, delta_;
+ //functional
+ //there is no divide_by_zero check for efficiency reasons. l0_ can not be zero upon call.
+ std::complex<value_t> distortion(const Vec2<value_t>& v) const
+ {
+ return std::polar<value_t>(delta_*std::tanh(v.r_wrt(position_)/l0_), alpha_ - v.phi_wrt(position_));
+ }
+ std::complex<value_t> kekule(const Vec2<value_t>& v, const Vec2<value_t>& s) const{
+ return std::polar<value_t>(1., K_*(s + 2.*(v - position_)));
+ }
+private:
+ Vec2<value_t> position_;
+ const Vec2<value_t> K_;
+ value_t l0_, alpha_, delta_;
};
#endif