Added jacobian stuff, and tasmet_variable, updated triplets

src/sys/CMakeLists.txt

@@ -0,0 +1,10 @@
+add_library(sys
+  # enginesystem.cpp
+  globalconf.cpp
+  tasmet_variable.cpp
+  jaccol.cpp
+  jacobian.cpp
+  jacrow.cpp
+  # tasystem.cpp
+  triplets.cpp
+  )

src/sys/globalconf.cpp

@@ -3,59 +3,64 @@
 #include "tasmet_exception.h"
 #include "tasmet_io.h"
 
-GlobalConf::GlobalConf(us Nf,d freq)
-
+GlobalConf::GlobalConf(us Nf,d freq):
+    _Nf(Nf)
 {
 
-    set(Nf,freq);
     TRACE(10,"GlobalConf constructor done");
+
+    if(Nf>=constants::maxNf)
+        throw TaSMETError("Too large number of frequencies given");
+
+
+    us Ns=this->Ns();
+    
+    // Reinitialize all operators
+    _iDFT=zeros<dmat>(Ns,Ns);
+    _fDFT=zeros<dmat>(Ns,Ns);
+
+    _fDFT.row(0).fill(1.0/double(Ns));
+
+    for(us i=1;i<=_Nf;i++){
+        for(us j=0; j<Ns;j++){
+            //Row i+1 (cosine components)
+            _fDFT(2*i-1,j)=2.0*cos(2.0*number_pi*double(i)*double(j)/double(Ns))/double(Ns);
+            //Row i (sine components)
+            _fDFT(2*i,j)=-2.0*sin(2.0*number_pi*double(i)*double(j)/double(Ns))/double(Ns);
+        }
+    }
+
+    _iDFT.col(0).fill(1.0);	// Steady part
+    for(us k=0;k<Ns;k++){
+        for (us n=1;n<=_Nf;n++){
+            _iDFT(k,2*n-1)=cos(2.0*number_pi*double(n)*double(k)/Ns);
+            _iDFT(k,2*n)=-sin(2.0*number_pi*double(n)*double(k)/Ns);
+        }
+    }
+
+    setfreq(freq);
+
 }
 void GlobalConf::show() const {
     cout << "------- Global configuration ------ \n";        
     cout << "------- Number of harmonics to solve for: "<< _Nf <<"\n";
     cout << "------- Fundamental frequency           : " << _omg/2/number_pi << " Hz\n";           
 }
-void GlobalConf::set(us Nf,d freq){
+void GlobalConf::setomg(d omg){
     TRACE(15,"GlobalConf::set(_Nf,freq)");
-    d omg = 2*number_pi*freq;
-    //ctor
+
+
+    _DDTfd=zeros<dmat>(Ns(),Ns());
+
     // Sanity checks
     if(omg<constants::minomg && omg>constants::maxomg)
         throw TaSMETError("Illegal frequency given");
-    if(Nf>=constants::maxNf)
-        throw TaSMETError("Too large number of frequencies given");
 
-    this->_Nf=Nf;
     this->_omg=omg;
 
-    us Ns=this->Ns();
-    // Reinitialize all operators
-    iDFT_=zeros<dmat>(Ns,Ns);
-    fDFT_=zeros<dmat>(Ns,Ns);
-
-    DDTfd_=zeros<dmat>(Ns,Ns);
-    fDFT_.row(0).fill(1.0/double(Ns));
-
     for(us i=1;i<=_Nf;i++){
-        for(us j=0; j<Ns;j++){
-            //Row i+1 (cosine components)
-            fDFT_(2*i-1,j)=2.0*cos(2.0*number_pi*double(i)*double(j)/double(Ns))/double(Ns);
-            //Row i (sine components)
-            fDFT_(2*i,j)=-2.0*sin(2.0*number_pi*double(i)*double(j)/double(Ns))/double(Ns);
-        }
-    }
-
-    iDFT_.col(0).fill(1.0);	// Steady part
-    for(us k=0;k<Ns;k++){
-        for (us n=1;n<=_Nf;n++){
-            iDFT_(k,2*n-1)=cos(2.0*number_pi*double(n)*double(k)/Ns);
-            iDFT_(k,2*n)=-sin(2.0*number_pi*double(n)*double(k)/Ns);
-        }
-    }
-
-    for(us i=1;i<=_Nf;i++){
-        DDTfd_(2*i-1,2*i  )=-double(i)*_omg;
-        DDTfd_(2*i  ,2*i-1)=double(i)*_omg;
+        _DDTfd(2*i-1,2*i  )=-double(i)*_omg;
+        _DDTfd(2*i  ,2*i-1)=double(i)*_omg;
     }
 
 }

src/sys/globalconf.h

@@ -10,14 +10,19 @@
 #ifndef _GLOBALCONF_H_
 #define _GLOBALCONF_H_
 
+#include <memory>
 #include "tasmet_types.h"
 #include "tasmet_tracer.h"
 
+class GlobalConf;
+
+typedef std::shared_ptr<GlobalConf> gc_ptr;
+
 class GlobalConf{
     d _omg;		// The "base" frequency in rad/s
     us _Nf;			// Number of frequencies to solve for
 
-    dmat fDFT_,iDFT_,DDTfd_;
+    dmat _fDFT,_iDFT,_DDTfd;
     // d Wfo_=0;			// First order 'upwind' factor. If
     // Wfo=-1, interpolation is done from
     // the left side. If Wfo=0,
@@ -26,6 +31,10 @@ class GlobalConf{
     // the right side
 public:
     GlobalConf(us Nf,d freq);
+
+    const dmat& iDFT = _iDFT; //inverse discrete Fourier transform matrix
+    const dmat& fDFT =  _fDFT; //forward discrete Fourier transform matrix
+    const dmat& DDTfd = _DDTfd; //Derivative in frequency domain
     
     us Nf() const {return _Nf;}
     us Ns() const {return 2*_Nf+1;}    
@@ -35,19 +44,11 @@ public:
     d getfreq() const {return _omg/2/number_pi;}
     // d meshPeclet(const Gas& gas,d dx,d u) const {return u*dx*gas.rho0()*gas().cp(T0())/gas().kappa(T0());}
 
-    void set(us Nf,d freq);	// Set data for new frequency and
-    // number of samples
-
-    void setNf(us Nf){set(Nf,getfreq());}  
-    void setfreq(d freq){set(Nf(),freq);}
-    void setomg(d omg)  {set(Nf(),omg/(2*number_pi));}
-
-    const dmat& iDFT() const {return iDFT_;} //inverse discrete Fourier transform matrix
-    const dmat& fDFT() const {return fDFT_;} //forward discrete Fourier transform matrix
-    const dmat& DDTfd() const {return DDTfd_;}//Derivative in frequency domain
+    void setfreq(d freq){setomg(2*number_pi*freq);}
+    void setomg(d omg);
 
     void show() const;
-
+    
 }; /* Class GlobalConf */
 
 #endif /* _GLOBALCONF_H_ */

src/sys/jaccol.cpp

@@ -0,0 +1,70 @@
+// jaccol.cpp
+//
+// last-edit-by: J.A. de Jong 
+// 
+// Description:
+//
+//////////////////////////////////////////////////////////////////////
+
+#include "jaccol.h"
+#include "tasmet_variable.h"
+
+
+JacCol::JacCol(const Variable& thevar):
+    coldof_(thevar.getDofNr()),
+    data_(thevar.getGc().Ns(),thevar.getGc().Ns(),fillwith::zeros)
+{  }
+JacCol::JacCol(us coldof,const GlobalConf& gc):
+    coldof_(coldof),
+    data_(gc.Ns(),gc.Ns(),fillwith::zeros)
+{  }
+JacCol::JacCol(us coldof,const dmat& data):
+    coldof_(coldof),
+    data_(data)
+{  }
+
+JacCol::JacCol(const Variable& thevar,const dmat& data):
+    coldof_(thevar.getDofNr()),
+    data_(data)
+{  }
+// JacCol::JacCol(JacCol&& j):
+//   tobeadded(std::move(j.tobeadded)),
+//   coldof_(j.coldof_),
+//   data_(std::move(j.data_))
+// {
+//   TRACE(45,"JacCol::JacCol(JacCol&&)");
+// }
+JacCol::JacCol(const JacCol& j):
+    tobeadded(j.tobeadded),
+    coldof_(j.coldof_),
+    data_(j.data_)
+{
+    TRACE(5,"JacCol::JacCol(const JacCol& j)");
+}
+JacCol& JacCol::operator=(const JacCol& j){
+    TRACE(10,"JacCol::operator=(const JacCol& j)");
+    data_=j.data_;
+    coldof_=j.coldof_;
+    tobeadded=j.tobeadded;
+    return *this;
+}
+// JacCol& JacCol::operator=(JacCol&& j){
+//   TRACE(45,"JacCol::operator=(const JacCol& j)");
+//   data_=std::move(j.data_);
+//   coldof_=j.coldof_;
+//   tobeadded=j.tobeadded;
+//   return *this;
+// }
+  
+JacCol& JacCol::operator+=(const dmat& data){
+    data_+=data;
+    return *this;
+}
+JacCol& JacCol::operator*=(const d& val){
+    data_*=val;
+    return *this;
+}
+
+
+
+//////////////////////////////////////////////////////////////////////

src/sys/jaccol.h

@@ -0,0 +1,47 @@
+// jaccol.h
+//
+// Author: J.A. de Jong 
+//
+// Description:
+// A column block of the Jacobian matrix
+//////////////////////////////////////////////////////////////////////
+#pragma once
+#ifndef JACCOL_H
+#define JACCOL_H
+#include "tasmet_types.h"
+
+class Variable;
+class GlobalConf;
+
+
+class JacCol{                // Column block of Jacobian
+    bool tobeadded=true;
+    us coldof_;                  // First dof of column
+    dmat data_;                  // data
+public:
+    JacCol(const Variable&);
+    JacCol(us coldof,const GlobalConf&);
+    JacCol(us coldof,const dmat&);    
+    JacCol(const Variable&,const dmat&);
+    // JacCol(JacCol&&);           // Move data
+    JacCol(const JacCol&);
+    void prePostMultiply(const dmat& pre,const dmat& post){ data_=pre*data_*post;}
+    JacCol& operator=(const JacCol&);
+    // JacCol& operator=(JacCol&&);
+
+    // Negate all terms
+    JacCol operator-();
+    bool isToAdd() const {return tobeadded;}
+    void setToAdd(bool set){tobeadded=set;} // If this is set to
+    // false, this column will not be added to the row.
+    JacCol& operator+=(const dmat& add);
+    JacCol& operator*=(const double&);
+    dmat& data(){return data_;}
+    const dmat& const_data() const {return data_;}
+    const us& getColDof() const{return coldof_;}
+    void show() const;
+};
+  
+
+#endif // JACCOL_H
+//////////////////////////////////////////////////////////////////////

src/sys/jacobian.cpp

@@ -0,0 +1,41 @@
+#include "jacobian.h"
+#include "triplets.h"
+#include "tasmet_tracer.h"
+
+void Jacobian::operator+=(const Jacobian& other){
+    TRACE(2,"Jacobian::append()");
+    jacrows.insert(jacrows.end(),other.jacrows.begin(),other.jacrows.end());
+}
+void Jacobian::operator+=(const JacRow& other){
+    TRACE(2,"Jacobian::append()");
+    jacrows.push_back(other);
+}
+Jacobian::operator TripletList() const{
+    TRACE(18,"Jacobian::operator Tripletlist()");
+    int insertrow,insertcol;
+    TripletList res(ndofs_);
+    const dmat& typicaldatacel=jacrows.at(0).jaccols.at(0).const_data();
+    us size=typicaldatacel.n_rows;
+
+    us i,j;
+
+    for(const JacRow& row: jacrows) {
+        insertrow=row.getRowDof();
+        for(const JacCol& col: row.jaccols){
+            insertcol=col.getColDof();
+            if(insertcol>=0){
+                const dmat& data=col.const_data();
+                for(i=0;i<size;i++){
+                    for(j=0;j<size;j++){
+                        if(data(i,j)!=0)
+                            res.addTriplet(Triplet(i+insertrow,j+insertcol,data(i,j)));
+                    }
+                }
+            } // insertcol>0
+        }   // for loop over cols
+    }     // for loop over rows
+    // TRACE(10,"SFSG");
+    return res;
+}
+  
+

src/sys/jacobian.h

@@ -0,0 +1,30 @@
+// jacobian.h
+//
+// Author: J.A. de Jong 
+//
+// Description:
+//
+//////////////////////////////////////////////////////////////////////
+#pragma once
+#ifndef JACOBIAN_H
+#define JACOBIAN_H
+
+#include "jacrow.h"
+
+class TripletList;
+
+class Jacobian{
+    us ndofs_;
+public:
+    Jacobian(us ndofs):ndofs_(ndofs){}
+    vector<JacRow> jacrows;
+    void operator+=(const Jacobian&);
+    void operator+=(const JacRow&);
+    operator TripletList() const;
+};
+
+  
+
+
+#endif // JACOBIAN_H
+//////////////////////////////////////////////////////////////////////

src/sys/jacrow.cpp

@@ -0,0 +1,58 @@
+#include "jacobian.h"
+#include "tasmet_variable.h"
+#include "globalconf.h"
+
+
+
+JacRow::JacRow(const JacCol& j):
+    JacRow(-1,1)
+{
+    (*this)+=j;
+}
+JacRow::JacRow(us rowdof,const JacCol& j):
+    JacRow(rowdof,1)
+{
+    (*this)+=j;
+}
+// JacRow& JacRow::operator+=(JacCol&& j){
+//   TRACE(45,"JacRow::operator+=(JacCol&& j)");
+//   jaccols.emplace_back(std::move(j));
+// }
+JacRow& JacRow::operator+=(const JacCol& j){
+    TRACE(10,"JacRow::operator+=(const JacCol& j)");
+    if(j.isToAdd())
+        jaccols.emplace_back(j);
+    return *this;
+}
+    
+// JacRow& JacRow::addCol(const JacCol& jaccol){
+//   if(jaccol.isToAdd())
+//     jaccols.push_back(jaccol);
+//   return *this;
+// }  
+JacRow& JacRow::operator*=(const d& val){
+    TRACE(2,"Jacobian::operator*=()");
+    for(JacCol& col: jaccols)
+        col.data()*=val;
+    return *this;
+}
+JacRow& JacRow::operator+=(const JacRow& other){
+    TRACE(2,"Jacobian::operator*=()");
+    this->jaccols.reserve(this->jaccols.capacity()+other.jaccols.size()-this->jaccols.size());
+    for(const JacCol& col :other.jaccols)
+        (*this)+=col;
+    return *this;
+}
+JacRow JacRow::operator-() const{
+    TRACE(15,"JacRow::operator-()");
+    JacRow result(*this);
+    for (JacCol& jaccol : result.jaccols) 
+        jaccol*=-1;
+
+    return result;
+}
+void JacRow::prePostMultiply(const dmat& pre,const dmat& post) {
+    for(JacCol& jaccol: jaccols)
+        jaccol.prePostMultiply(pre,post);
+}
+

src/sys/jacrow.h

@@ -0,0 +1,44 @@
+// jacrow.h
+//
+// Author: J.A. de Jong 
+//
+// Description:
+//
+//////////////////////////////////////////////////////////////////////
+#pragma once
+#ifndef JACROW_H
+#define JACROW_H
+
+#include "jaccol.h"
+
+class Variable;
+
+  
+class JacRow{                 // Row in Jacobian matrix
+    int rowdof_=-1;            // Number of first row, default is
+    // invalid state
+public:
+    // Negate all terms
+    JacRow operator-() const;
+
+    vector<JacCol> jaccols;     // Column blocks
+    JacRow(const JacCol&);
+    JacRow(us rowdof,const JacCol&); // Initialization with only one column
+    JacRow(int rowdofnr,us cols=2): rowdof_(rowdofnr){ jaccols.reserve(cols);}
+    // void addCol(const JacCol& jaccol);
+    // JacRow& operator+=(JacCol&&);
+    JacRow& operator+=(const JacCol&);
+    JacRow& operator+=(const JacRow& jacrow);
+    JacRow& operator*=(const d& val); // Multiply all terms with constant value
+    void prePostMultiply(const dmat& pre,const dmat& post);
+    const int& getRowDof() const {return rowdof_;}
+    void setRowDof(us dofnr){rowdof_=dofnr;}
+    void show() const;
+};
+
+  
+
+
+
+#endif // JACROW_H
+//////////////////////////////////////////////////////////////////////

src/sys/tasmet_variable.cpp

@@ -0,0 +1,225 @@
+// var.cpp
+//
+// last-edit-by: J.A. de Jong 
+// 
+// Description:
+//
+//////////////////////////////////////////////////////////////////////
+#include "tasmet_variable.h"
+#include "tasmet_exception.h"
+#define Ns (_gc->Ns())
+#define Nf (_gc->Nf())
+#define fDFT (_gc->fDFT)
+#define iDFT (_gc->iDFT)
+#define DDTfd (_gc->DDTfd)
+
+//******************************************************************************** Operators
+  
+Variable operator*(const double& scalar,const Variable& other){ // Pre-multiplication with scalar
+    TRACE(0,"Variable::operator*(scalar,Variable)");
+    return Variable(other._gc,scalar*other._adata);
+}
+Variable Variable::operator/(const Variable& Variable2) const {
+    TRACE(0,"Variable::operator/()");
+    return Variable(this->_gc,this->_tdata/Variable2._tdata,false);
+}
+Variable Variable::operator/(d val) const {
+    TRACE(0,"Variable::operator/(double)");
+    return Variable(this->_gc,this->tdata()/val,false);
+}
+Variable Variable::operator-(const Variable& other) const{
+    TRACE(0,"Variable::operator-(Variable)");
+    return Variable(this->_gc,_adata-other._adata);
+}
+Variable Variable::operator*(d scalar) const {
+    TRACE(0,"Variable::operator*(scalar)");	// Post-multiplication with scalar
+    return Variable(this->_gc,scalar*_adata);
+}
+Variable Variable::operator+(const Variable& other) const{
+    TRACE(0,"Variable::operator+()");
+    return Variable(this->_gc,_adata+other._adata);
+}
+Variable Variable::operator*(const Variable& Variable2) const { // Multiply two
+    TRACE(0,"Variable::operator*(const Variable& Variable2) const");
+    return Variable(this->_gc,_tdata%Variable2._tdata,false);
+}
+//***************************************** The Variable class
+Variable::Variable(const gc_ptr& gc,double initval): _gc(gc)
+{
+    TRACE(0,"Variable::Variable(us ndofs, double initval)");
+
+    _tdata=initval*ones(Ns);
+    _adata=zeros(Ns);
+    _adata(0)=initval;
+}
+Variable::Variable(const gc_ptr& gc,const vd& data,bool adata): _gc(gc)
+{ // Create a tasystem and fill it with time data.
+    TRACE(0,"Variable::Variable(gc,_tdata)");
+    if(adata){
+        this->_adata=data;
+        _tdata=iDFT*_adata;
+    }
+    else{
+        this->_tdata=data;
+        _adata=fDFT*_tdata;
+    }
+}
+Variable::Variable(const gc_ptr& gc,const vc& data)
+    :Variable(gc)
+{ // Create a tasystem and fill it with time data.
+    if(data.size()!=Nf+1)
+        throw TaSMETError("Wrong size of amplitude vector given. Does the"
+                      "vector size correspond to Ns?");
+    setadata(data);
+}
+void Variable::setGc(const gc_ptr& gc){
+    TRACE(10,"Variable::setGc()");
+    this->_gc = gc;
+    updateNf();
+}
+void Variable::resetHarmonics(){
+    TRACE(10,"Variable::resetHarmonics()");
+    if(Nf>1){
+        for(us i=3;i<Ns;i++)
+            _adata(i)=0;
+    }
+    _tdata=iDFT*_adata;
+}
+Variable::Variable(const Variable& other){
+    TRACE(0,"Variable::Variable(const Variable& other)");
+    this->_gc=other._gc;
+    this->_tdata=other._tdata;
+    this->_adata=other._adata;
+}
+Variable& Variable::operator=(const Variable& other){
+    // THIS WOULD COUPLE TO THE WRONG GLOBALCONF when setRes is used
+    // between ducts!!!!!
+    if(this!=&other){
+        if(!this->_gc) {
+            this->_gc=other._gc;
+        }
+        this->_tdata=other._tdata;
+        this->_adata=other._adata;
+        updateNf();
+    }
+    return *this;
+}
+void Variable::updateNf(){
+    TRACE(0,"Variable::updateNf()");
+    us oldsize=_adata.size();
+    assert(oldsize==_tdata.size());
+    assert(_gc);
+    us newsize=Ns;
+    if(oldsize!=newsize){
+        _adata.resize(newsize);
+        _tdata=iDFT*_adata;
+    }
+}
+// Get methods (which require implementation)
+d Variable::operator()(us i) const {//Extract result at specific frequency
+    TRACE(-2,"Variable::operator()("<<i<<"), Ns: "<< Ns);
+    if(i>=Ns)
+        throw TaSMETError("Invalid frequency number!");
+    TRACE(-1,"_adata: "<<_adata);
+    return _adata(i);
+}  
+vc Variable::getcRes() const
+{
+    TRACE(-2,"Variable::getcRes()");
+    //	TRACE(0,"_adata:" << _adata);
+    vc cadata(Nf+1);
+    cadata(0)=_adata(0);
+    for(us i=1;i<Nf+1;i++)
+        cadata(i)=_adata(2*i-1)+I*_adata(2*i); //The minus is very important
+    //	TRACE(0,"Resulting cadata:" << cadata);
+    return cadata;
+}
+// Set methods
+void Variable::setadata(us freqnr,d val) { //Set result for specific frequency zero,real one, -imag one, etc
+    TRACE(-2,"Variable::setadata("<<freqnr<<","<<val<<")");
+    if(freqnr>=Ns)
+        throw TaSMETError("Tried to change a value beyond length of the vector");
+    _adata[freqnr]=val;
+    _tdata=iDFT*_adata;
+}
+void Variable::setadata(const vc& res)
+{
+    TRACE(0,"Variable::setadata(const vc& res)");
+    assert(res.size()==Nf+1);
+    _adata(0)=res(0).real();
+    for(us i=1;i<Nf+1;i++){
+        _adata(2*i-1)=res(i).real();
+        _adata(2*i)=res(i).imag();
+    }
+    _tdata=iDFT*_adata;
+}
+d min(us a,us b){return a<=b?a:b;}
+void Variable::setadata(const vd& val) {
+    TRACE(0,"Variable::setadata(const vd& val)");
+    _adata.zeros();
+    us minsize=min(val.size(),_adata.size());
+    for(us i=0;i<minsize;i++)
+        _adata(i)=val(i);
+    _tdata=iDFT*_adata;
+}
+void Variable::settdata(double val) {
+    TRACE(0,"Variable::settdata(double val)");
+    for (auto& td: _tdata) {
+        td=val;
+    }
+    _adata=fDFT*_tdata;
+}
+void Variable::settdata(const vd& val) {
+    TRACE(0,"Variable::settdata(vd& val)");
+    assert(val.size()==_tdata.size());
+    _tdata=val;
+    _adata=fDFT*_tdata;
+}
+
+vd Variable::timeResponse(us nperiod,us ninst) const {
+    TRACE(15,"vd Variable::timeResponse()");
+    vd t=timeResponseTime(nperiod,ninst);
+    vc cres=getcRes();
+    vd res(t.size(),fillwith::zeros);
+    c omg=_gc->getomg();
+    for(us i=0;i<Nf+1;i++)
+        res+=real(cres(i)*exp(((d) i)*omg*I*t));
+    return res;
+}
+vd Variable::timeResponseTime(us nperiod,us ninst) const{
+    d T=1/_gc->getfreq();
+    return linspace(0,nperiod*T,ninst);
+}
+dmat Variable::freqMultiplyMat() const{
+    TRACE(0,"Variable::freqMultiplyMat()");
+    dmat result(Ns,Ns,fillwith::zeros);
+    result(0,0)=_adata(0);
+    if(Nf>0){
+        for(us j=1;j<Nf+1;j++){
+            result(2*j-1,2*j-1)= _adata(2*j-1);
+            result(2*j-1,2*j  )=-_adata(2*j  ); // Yes only one
+            // minus sign
+            result(2*j  ,2*j-1)= _adata(2*j);      
+            result(2*j  ,2*j  )= _adata(2*j-1);      
+        }	// end for loop
+    }	// if(Nf>0)
+    return result;
+}
+
+//Get a tasystem which is the time derivative of the current one
+Variable Variable::ddt() const {
+
+    vd newadata=DDTfd*_adata;
+
+    return Variable(_gc,newadata);
+}
+// The product
+
+//***************************************** End of the Variable class
+
+
+
+
+
+
+

src/sys/tasmet_variable.h

@@ -0,0 +1,115 @@
+// tasmet_variable.h
+//
+// Author: J.A. de Jong 
+//
+// Description:
+//
+//////////////////////////////////////////////////////////////////////
+#pragma once
+#ifndef TASMET_VARIABLE_H
+#define TASMET_VARIABLE_H
+
+#include "globalconf.h"
+#include <tasmet_types.h>
+#include <tasmet_assert.h>
+
+class Variable;			// Forward declaration
+Variable operator*(const double&,const Variable&);
+// Variable operator+(const d&,const Variable&);
+
+class Variable {
+
+    friend Variable operator*(const double&,const Variable&);
+
+    int _dofnr=-1;
+    vd _tdata,_adata;
+
+    gc_ptr _gc;
+
+public:
+    void setDofNr(us Dofnr){_dofnr=Dofnr;}
+    int getDofNr() const{return _dofnr;}
+
+    Variable(const gc_ptr&,double); // Initialize with one time-average value
+    Variable(const Variable& o);
+    Variable(const gc_ptr&);	// Initialize with zeros
+
+    // Initialize with amplitudedata. With _tdata if adata is set to false
+    Variable(const gc_ptr&,const vd& data,bool adata=true);
+
+    // Assign with frequency domain data
+    Variable(const gc_ptr&,const vc& data);
+    
+    Variable& operator=(const Variable&);
+      
+    ~Variable(){}
+
+    void setGc(const gc_ptr&);
+    const GlobalConf& getGc() const {return *_gc;}
+    void updateNf();
+    
+    // Get methods
+    d operator()(us i) const; // Extract amplitude data result at
+                              // specific frequency
+
+    operator vd() const {return _adata;}
+
+    // Extract data
+    const vd& tdata() const  {return _tdata; } //Get time data
+    const vd& adata() const  {return _adata; } //Get time data                                                 //vector
+
+    // Obtain a time response vector
+    vd timeResponse(us nperiod=2,us ninst=100) const;
+    // Obtain the corresponding time vector
+    vd timeResponseTime(us nperiod=2,us ninst=100) const;
+
+    dmat diagt() const {return diagmat(_tdata);}
+    dmat diag() const {return diagmat(_adata);}    
+
+    //Set/reset methods
+    void resetHarmonics();
+    void updateNDofs(us new_ndofs);
+
+    void setadata(const vd& values); //Set amplitude data vector to these values
+    void settdata(double value); //Set time data to specific value for all time
+    void settdata(const vd& values);
+    void setadata(us freq,double val); //Set result vector at specific frequency
+    void setadata(const vc& values); //Set result vector to these values,
+
+    us size() const {return _adata.size();}
+    // Specific methods to the result using time domain data
+    // Operations ********************
+
+    Variable ddt() const;       // Time-derivative of this variable
+    
+    Variable operator/(const Variable& Variable2) const; // Time-domain division operator
+    Variable operator/(d value) const; // Time-domain or frequency domain division operator
+
+    // Multiply two Variables in time domain
+    Variable operator*(const Variable& tasystem) const;
+
+    // Multiply a tasystem with a scalar. This operation is possible
+    // for both frequency and time domain data
+    Variable operator*(d scalar) const;
+
+    // add two Variableiables
+    Variable operator+(const Variable& other) const;
+    //Subtract two Variableiables
+    Variable operator-(const Variable& Variable2) const;
+    // with Note multiplication is defined outside of the class
+
+    // If we need to multiply two numbers in frequency domain, this
+    // corresponds to a matrix-vector multiplication (cosines and
+    // sines) are mixed up due to the complex numbers. This product
+    // can be obtained by getting the matrix-Variable of the first
+    // Variable. The following function will give the effective matrix
+    dmat freqMultiplyMat() const;
+
+    // Get the result vector in complex form
+    vc getcRes() const;
+    
+};
+
+
+#endif /* TASMET_VARIABLE_H_ */
+

src/sys/triplets.cpp

@@ -3,6 +3,9 @@
 #include "tasmet_io.h"
 #include "tasmet_exception.h"
 
+TripletList::TripletList(us ndofs): _ndofs(ndofs) {
+    TRACE(15,"TripletList::TripletList()");
+}
 TripletList::operator sdmat() const {
       
     TRACE(15,"TripletList::operator sdmat()");
@@ -60,11 +63,12 @@ void TripletList::show()  const {
         cout << "Row: " << t.row << " , column: " << t.col << " , value: " << t.value << "\n";
     }
 }
-void TripletList::multiplyTriplets(const d& factor){
+TripletList& TripletList::operator*=(d factor){
     TRACE(15,"multiplyTriplets()");
     for(auto tr: triplets){
         tr.value*= factor;
     }
+    return *this;
 }
 
 void TripletList::reserveExtraDofs(us n){
@@ -77,8 +81,12 @@ void TripletList::shiftTriplets(int nrows,int ncols){
     // shift the position of the values in a matrix. nrows and ncols
     // can be negative numbers.
     TRACE(15,"shiftTriplets()");
-    TRACE(100,"EXTRA CHECKS HERE!");
     for(auto tr: triplets){
+        #if TASMET_DEBUG == 1
+        if(tr.col+ncols >= _ndofs || (tr.row+nrows >= _ndofs)) {
+            FATAL("Out of bounds shift");
+        }
+        #endif
         tr.col+=ncols;
         tr.row+=nrows;
     }

src/sys/triplets.h

@@ -36,7 +36,7 @@ public:
     // Make one row zero
     void zeroOutRow(us rownr);
 
-    void multiplyTriplets(const d& multiplicationfactor);
+    TripletList& operator*=(d multiplicationfactor);
 
     // Add to capacity
     void reserveExtraDofs(us n);