Added matrix calculation

include/matrix_element_calculator.hpp

+#ifndef MATRIX_ELEMENT_CALCULATOR_HPP_INCLUDED
+#define MATRIX_ELEMENT_CALCULATOR_HPP_INCLUDED
+
+template <value_t, params_t, vars_t, drum_t>
+class MatrixElementCalculator{
+        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 element at (index1, index2)
+                virtual value_t operator()(const size_t index1, const size_t index2, const std::vector<drum_t>& drums) const noexcept = 0
+};
+
+#endif

include/system.hpp

@@ -3,7 +3,7 @@
 #include <vector>
 #include <iostream>
 
-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>
+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)
@@ -50,11 +50,17 @@ class System{
                         return grabber_.save();
                 }
 
+                std::vector<value_t> get_matrix(const matele_t& mec){
+                        calculate_matrix();
+                        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);
@@ -64,6 +70,32 @@ class System{
                         }
                 }
 
+                //Calculates dynamical matrix for the currently pushed values
+                void calculate_matrix(const matelecalc_t& mec){
+                        //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_[i][i] = mec(i, drums_);
+                                //couplings
+                                if(n0_index != drums_.size()-1) [[likely]] //actually a neighbour
+                                        matrix_[i][n0_index] = mec(i, n0_index, drums_);
+                                if(n1_index != drums_.size()-1) [[likely]] //actually a neighbour
+                                        matrix_[i][n1_index] = mec(i, n1_index, drums_);
+                                if(n2_index != drums_.size()-1) [[likely]] //actually a neighbour
+                                        matrix_[i][n2_index] = mec(i, n2_index, drums_);
+                        }
+                }
+
                 std::vector<drum_t> drums_; //The last drum is NOT TO BE TOUCHED!
                 sysparams_t sysparams_;
                 grabber_t grabber_;
@@ -71,6 +103,8 @@ class System{
                 force_t force_;
                 value_t dt_, t_end_;
                 value_t time_;
+
+                std::vector<value_t> matrix_; //Matrix of this problem
 };
 
 #endif