Working till generation of the Jacobian. We are going to use has tables stuff for Jacobian matrix to increase speed

src/duct/duct.cpp

@@ -231,7 +231,6 @@ vd Duct::initialSolution() const {
 
 us Duct::getNEqs() const {
     TRACE(15,"Duct::getNEqs()");    
-    us Ns = sys.Ns();
 
     // The number of equations per gridpoint. We have: continuity,
     // momentum, energy, and state
@@ -240,15 +239,15 @@ us Duct::getNEqs() const {
     // When we have to solve a solid heat balance
     number_eqs+= (_ductpb.stempmodel() == pb::HeatBalance ? 1 : 0);
 
-    us neqs = Ns*number_eqs*(ngp()-1);
+    us neqs = number_eqs*(ngp()-1);
 
     // For the last gridpoint, we also have an equation of state
-    neqs += Ns;
+    neqs += 1;
 
     // We also have an extra equation for isentropic. For the energy
     // transport equation, this would result in a boundary condition
     if(_ductpb.htmodel() == pb::Isentropic) {
-        neqs += Ns;
+        neqs += 1;
     }
 
     VARTRACE(15,neqs);
@@ -256,14 +255,13 @@ us Duct::getNEqs() const {
 }
 us Duct::getNDofs() const {
     TRACE(15,"Duct::getNDofs()");        
-    us Ns = sys.Ns();  
  
     // rho,u,T,p
     us nvars_per_gp = 4;
     // Ts maybe
     nvars_per_gp += (_ductpb.stempmodel() == pb::HeatBalance ? 1 : 0);
 
-    return Ns*nvars_per_gp*ngp();
+    return nvars_per_gp*ngp();
 }
 d Duct::getMass() const {
 

src/duct/duct.h

@@ -47,7 +47,7 @@ public:
      *
      * @return The index of the degree of freedom
      */
-    Dof getDof(int varnr,int gp) const;
+    PosId getDof(int varnr,int gp) const;
 
     /// Postprocessing functions
     vd rhox(int t) const { return getvarx(constants::rho,t);}

src/ductbc/adiabaticwall.cpp

@@ -43,12 +43,12 @@ void AdiabaticWall::residualJac(arma::subview_col<d>&& residual
                              ) const {
 
     TRACE(15,"AdiabaticWall::residual()");
-    const pb::Duct& dpb = getDuct(sys).getDuctPb();
+    const Duct& duct = getDuct();
+    const pb::Duct& dpb = duct.getDuctPb();
     us Ns = sys.Ns();
 
-    const Duct& duct = getDuct(sys);
     if(_side == pb::left) {
-        residual.subvec(0,Ns-1) = duct.ut(sys,0);
+        residual.subvec(0,Ns-1) = duct.ut(0);
         if(dpb.htmodel() != pb::Isentropic) {
             // TODO: Put the boundary condition of zero heat flux here
             // residual.subvec(Ns,2*Ns-1) = 
@@ -56,7 +56,7 @@ void AdiabaticWall::residualJac(arma::subview_col<d>&& residual
         }
     }
     else {
-        residual.subvec(0,Ns-1) = duct.ut(sys,-1);
+        residual.subvec(0,Ns-1) = duct.ut(-1);
         if(dpb.htmodel() != pb::Isentropic) {
             // TODO: Put the boundary condition of zero heat flux here
             // residual.subvec(Ns,2*Ns-1) = duct.Tt(sys,-1) - _T;
@@ -72,12 +72,12 @@ us AdiabaticWall::getNEqs() const {
     // dT/dx = 0 --> if htmodel is not Isentropic
     // dTs/dx = 0 => if stempmodel is not Prescribed
 
-    const pb::Duct& dpb = getDuct(sys).getDuctPb();
+    const pb::Duct& dpb = getDuct().getDuctPb();
 
     bool has_solideq = dpb.stempmodel() != pb::Prescribed;
 
-    us neqs = sys.Ns()*(has_solideq ? 2: 1);
-    if(dpb.htmodel() != pb::Isentropic) neqs+= sys.Ns();
+    us neqs = (has_solideq ? 2: 1);
+    if(dpb.htmodel() != pb::Isentropic) neqs+= 1;
 
     VARTRACE(15,neqs);
     return neqs;

src/ductbc/pressurebc.cpp

@@ -110,8 +110,8 @@ us PressureBc::getNEqs() const {
 
     bool has_solideq = dpb.stempmodel() != pb::Prescribed;
 
-    us neqs = sys.Ns()*(has_solideq ? 2: 1);
-    if(dpb.htmodel() != pb::Isentropic) neqs+= sys.Ns();
+    us neqs = (has_solideq ? 2: 1);
+    if(dpb.htmodel() != pb::Isentropic) neqs+= 1;
 
     VARTRACE(15,neqs);
     return neqs;

src/sys/jaccol.cpp

@@ -7,63 +7,6 @@
 //////////////////////////////////////////////////////////////////////
 
 #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

@@ -10,36 +10,39 @@
 #define JACCOL_H
 #include "tasmet_types.h"
 
-class Variable;
-class GlobalConf;
-
+/**
+ * A PosId is a specifier for the position in the global Jacobian
+ * matrix. Both considering the column as well as the row
+ * 
+ */
+struct PosId {
+    int segno;
+    int segoffset;
+};
 
 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;
+public:
+
+    // For more ease of use, we make the members public.
+
+    PosId id;                  /**< The Id of this DOF */
+    dmat data;                 /**< The matrix block */
+
+    JacCol(PosId id): id(id){}
+    JacCol(PosId id,const dmat& data):id(id),data(data){}
+
+    JacCol& operator=(const JacCol& o) { id = o.id; data=o.data; return *this; }
+    JacCol(const JacCol& o): id(o.id),data(o.data){}
+
+    void prePostMultiply(const dmat& pre,const dmat& post){ data=pre*data*post;}
+
+    // Some operations on the data
+    JacCol operator-() { return JacCol(id,-data);}
+
+    JacCol& operator+=(const dmat& add) {data+=add; return *this;}
+    JacCol& operator*=(const double d) {data*=d; return *this;}
+
 };
   
 

src/sys/jacobian.cpp

@@ -10,32 +10,34 @@ 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;
+// Jacobian::operator TripletList() const{
 
-    us i,j;
+//     TRACE(18,"Jacobian::operator Tripletlist()");
 
-    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;
-}
+//     int insertrow,insertcol;
+//     TripletList res(ndofs_);
+//     const dmat& typicaldatacel=jacrows.at(0).jaccols.at(0).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.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

@@ -14,12 +14,14 @@
 class TripletList;
 
 class Jacobian{
-    us ndofs_;
+    us _ndofs;
 public:
-    Jacobian(us ndofs):ndofs_(ndofs){}
+    Jacobian(us ndofs):_ndofs(ndofs){}
     vector<JacRow> jacrows;
+
     void operator+=(const Jacobian&);
     void operator+=(const JacRow&);
+
     operator TripletList() const;
 };
 

src/sys/jacrow.cpp

@@ -3,48 +3,30 @@
 #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);
+    jaccols.push_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){
+JacRow& JacRow::operator*=(const d val){
     TRACE(2,"Jacobian::operator*=()");
     for(JacCol& col: jaccols)
-        col.data()*=val;
+        col*=val;
     return *this;
 }
 JacRow& JacRow::operator+=(const JacRow& other){
-    TRACE(2,"Jacobian::operator*=()");
+    // Reserve some space
     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;

src/sys/jacrow.h

@@ -8,32 +8,43 @@
 #pragma once
 #ifndef JACROW_H
 #define JACROW_H
-
+#include <vector>
 #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
+
+    PosId id;                   /**< Specifies the row for this
+                                   JacRow */
+    std::vector<JacCol> jaccols; /**< The columns in this row */
+
+    /** 
+     * Initialize a Jacobian row with a single column block
+     *
+     * @param id The equation ID
+     * @param col The added column
+     */
+    JacRow(PosId id,const JacCol& col):JacRow(id){jaccols.push_back(col);}
+
+    JacRow(PosId id): id(id){}
+
     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
+    JacRow& operator*=(const d val); // Multiply all terms with constant value
+
+    /** 
+     * Pre- and postmultiply all Column blocks with the given
+     * matrices. Hence the new columns will be col_new = pre*col_old*post
+     *
+     * @param pre : The pre-multiply matrix
+     * @param post: The post-multiply matrix
+     */
     void prePostMultiply(const dmat& pre,const dmat& post);
-    const int& getRowDof() const {return rowdof_;}
-    void setRowDof(us dofnr){rowdof_=dofnr;}
-    void show() const;
+
+
 };