Not switching to Qwt. We use QCustomplot

CMakeLists.txt

@@ -14,6 +14,7 @@ add_definitions(-DTRACER=1)
 
 add_definitions(-DTASMET_FLOAT=64)
 add_definitions(-DTASMET_DEBUG=1)
+add_definitions(-DARMA_DONT_USE_WRAPPER)
 #====================================================
 # Compiler settings ************************************************
 #====================================================
@@ -48,12 +49,14 @@ set(CMAKE_AUTOMOC ON)
 set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O2 -march=native -mtune=native")
 # set(CMAKE_CLANG "${CMAKE_GCC} -march=native -mtune=native -fopenmp")
 
+
 # To disable bound checking on std::vector, and to disable assertions
 # add_definitions(-DNDEBUG)
 
 # To increase speed on Armadillo
 # add_definitions(-DARMA_NO_DEBUG)
 
+
 # Disable traces
 # add_definitions(-DTRACER=0)
 
@@ -128,5 +131,5 @@ add_subdirectory(testing)
 
 # add_executable(tasmet src/main.cpp src/main.moc)
 add_executable(tasmet src/main.cpp)
-target_link_libraries(tasmet tasmet_gui tasmet_src messages PythonQt Qt5::Widgets)
+target_link_libraries(tasmet tasmet_gui tasmet_src messages PythonQt Qt5::Widgets openblas)
 

src/duct/duct.cpp

@@ -7,33 +7,95 @@
 //////////////////////////////////////////////////////////////////////
 #include "duct.h"
 #include "tasystem.h"
+#include "tasmet_assert.h"
+#include "tasmet_pyeval.h"
 
-Duct::Duct(const pb::Duct& duct):
+Duct::Duct(const us id,const pb::Duct& duct):
-    Segment(duct.name()),
+    Segment(id,duct.name()),
     Geom(duct)
 {
 
-    
+    const char* invTsfun = "Invalid solid-temperature prescribing function";
+    EvaluateFun Tsfun(duct.stempfunc(),invTsfun);
+    Tsfun.addGlobalDef("L",duct.length());
+    _Tsprescribed = Tsfun(x);
+    if(min(_Tsprescribed) < constants::min_T0 ||
+       max(_Tsprescribed) > constants::max_T0) {
+        throw TaSMETError(invTsfun);
+    }
+
+
+    switch (duct.htmodel()) {
+    case pb::Isentropic: {
+        break;
+    }
+    default:
+        tasmet_assert(false,"Invalid heat transfer model");
+    break;
+    }
 
 
 }
-void Duct::updateSolution(const TaSystem&,const vd&) const {
-
-}
-vd Duct::getSolution(const TaSystem& sys,vd& sol,const us insert_start) {
-
-    
-
-}
-
-
 Duct::Duct(const Duct& other):
     Segment(other),
-    Geom(other)
+    Geom(other),
+    _Tsprescribed(other._Tsprescribed)
 {
+    // Do something with the equations here
+}
+Duct* Duct::copy() const {
+    return new Duct(*this);
+}
+Duct::~Duct() {
+    // for(Equation* eq: _eqs){
+    //     delete eq;
+    // }
+}
+void Duct::residual(const TaSystem& sys,arma::subview_col<d> && residual) const {
 
+    const arma::subview_col<d> sol = sys.getSolution(_id);
+
+    VARTRACE(15,sol(0));
+    
 
 }
+vd Duct::initialSolution(const TaSystem& sys) const {
+    
+    vd initsol(getNDofs(sys));
+
+    us vars_per_gp = 5;
+    us Ns = sys.Ns();
+
+    const Gas& gas = sys.gas();
+    
+    for(us i=0;i<ngp();i++){
+
+        // Initial density
+        initsol.subvec((i*vars_per_gp+0)*Ns,(i*vars_per_gp+1)*Ns-1) =
+            gas.rho0();
+        
+        // Initial velocity
+        initsol.subvec((i*vars_per_gp+1)*Ns,(i*vars_per_gp+2)*Ns-1) =
+            0;
+
+        // Initial Temperature
+        initsol.subvec((i*vars_per_gp+2)*Ns,(i*vars_per_gp+2)*Ns-1) =
+            _Tsprescribed(i);
+
+        // Initial pressure
+        initsol.subvec((i*vars_per_gp+3)*Ns,(i*vars_per_gp+4)*Ns-1) =
+            gas.p0();
+        
+        // Initial solid temperature
+        initsol.subvec((i*vars_per_gp+4)*Ns,(i*vars_per_gp+5)*Ns-1) =
+            _Tsprescribed(i);
+
+    }
+
+    return initsol;
+
+}
+
 us Duct::getNEqs(const TaSystem& sys) const {
 
     us Ns = sys.Ns();
@@ -41,8 +103,25 @@ us Duct::getNEqs(const TaSystem& sys) const {
 
     return Ns*number_eqs*ngp();
 }
+us Duct::getNDofs(const TaSystem& sys) const {
+    
+    us Ns = sys.Ns();  
 
-Duct* Duct::copy() const {
+    // rho,u,T,p,Ts
-    return new Duct(*this);
+    us nvars_per_gp = 5;
+
+    return Ns*nvars_per_gp*ngp();
 }
+d Duct::getMass(const TaSystem& sys) const {
+
+    return 0;
+}
+void Duct::jac(const TaSystem&,Jacobian&,us dof_start,us eq_start) const {
+
+}
+void Duct::show(const TaSystem&,us verbosity_level) const {
+
+}
+
+
 //////////////////////////////////////////////////////////////////////

src/duct/duct.h

@@ -19,28 +19,31 @@ class TaSystem;
 
 class Duct : public Segment, public Geom {
 
-    Duct(const Duct& other);
-
     // Drag* _drag = nullptr;
     // Heat* _heat = nullptr;
     std::vector<Equation*> _eqs;
+    vd _Tsprescribed;
 
-    std::vector<Variable> _rho;
+protected:
-    std::vector<Variable> _u;
+    Duct(const Duct&);
-    std::vector<Variable> _T;
-    std::vector<Variable> _p;
-    std::vector<Variable> _Ts;
-
 public:
-    Duct(const pb::Duct&);
+    Duct(const us id,const pb::Duct&);
+    ~Duct();
+
+    vd Tsprescribed() const { return _Tsprescribed;}
+
+    // Initialize the solution to something sensible
+    vd initializeSolution(const TaSystem& sys);
+
     virtual Duct* copy() const;
     const Geom& geom() const;
     
     // Solving
-    virtual void residual(const TaSystem&,vd&,const us insertion_start) const;
+    virtual void residual(const TaSystem&,arma::subview_col<d>&& residual) const;
 
-    virtual void updateSolution(const TaSystem&,const vd&);
+    vd initialSolution(const TaSystem&) const;
-    virtual getSolution(const TaSystem&,vd& sol,const us insertion_start) const;
+
+    // virtual void getSolution(const TaSystem&,const us insertion_start) const;
 
     // Return the total number of equations in this segment
     virtual us getNEqs(const TaSystem&) const;

src/duct/geom.cpp

@@ -34,24 +34,32 @@ Geom::Geom(const pb::Duct& duct) {
     // Store x
     x = grid->getx();
 
-    EvaluateFun Sfun(duct.area());
+    const char* invS = "Invalid cross-sectional area function";
+    EvaluateFun Sfun(duct.area(),invS);
     Sfun.addGlobalDef("L",duct.length());
     S = Sfun(x);
 
-    for(auto& Si:S) {
+    if(min(S) <= 0) {
-        if(Si <=0){
+        throw TaSMETError(invS);
-            throw TaSMETError("Invalid cross-sectional area function");
-        }
     }
+    
 
-
+    const char* invpor = "Invalid porosity function";
-    EvaluateFun phifun(duct.phi());
+    EvaluateFun phifun(duct.phi(),invpor);
     Sfun.addGlobalDef("L",duct.length());
     phi = Sfun(x);
 
-    EvaluateFun rhfun(duct.rh());
+    if(min(phi) <= 0 || max(phi) > 1) {
+        throw TaSMETError(invpor);
+    }
+
+    const char* invrh = "Invalid hydraulic radius function";
+    EvaluateFun rhfun(duct.rh(),invrh);
     Sfun.addGlobalDef("L",duct.length());
     rh = Sfun(x);
+    if(min(rh) <= 0 ) {
+        throw TaSMETError(invrh);
+    }
 
     cshape = duct.cshape();
 }

src/duct/grid.cpp

@@ -8,13 +8,14 @@ LinearGrid::LinearGrid(us ngp,d L):
     ngp(ngp),
     L(L)
 {
-    std::stringstream error;
     if(ngp<constants::min_ngp) {
+        std::stringstream error;
         error << "Number of gridpoints is lower than minimum. Minimum is: ";
         error << constants::min_ngp;
         throw TaSMETError(error);
     }
     else if(ngp>constants::max_ngp) {
+        std::stringstream error;
         error << "Maximum number of gridpoints exceeded. Maximum is: ";
         error << constants::max_ngp;
         throw TaSMETError(error);

src/gui/add_duct_dialog.cpp

@@ -7,6 +7,7 @@
 //////////////////////////////////////////////////////////////////////
 #include "ui_add_duct_dialog.h"
 #include "add_duct_dialog.h"
+#include "ui_add_duct_dialog.h"
 
 #include <QSignalBlocker>
 #include <qcustomplot.h>
@@ -15,11 +16,9 @@
 #include "tasmet_tracer.h"
 #include "tasmet_assert.h"
 #include "tasmet_exception.h"
-#include "duct/grid.h"
+#include "duct/duct.h"
-#include <QVector>
 #include "tasmet_qt.h"
-#include "geom.h"
+#include <QMessageBox>
-
 DECLARE_ENUM(PreviewShow,CSArea,Porosity,HydraulicRadius,SolidTemperatureFunction)
 
 
@@ -125,6 +124,20 @@ AddDuctDialog::~AddDuctDialog(){
 }
 void AddDuctDialog::accept(){
     
+    try {
+        std::unique_ptr<Duct> duct (new Duct(0,_duct));
+    }
+    catch(TaSMETError& e) {
+
+ 	QMessageBox msg(QMessageBox::Warning,
+                        "Input parsing error",
+                        e.what());
+
+        msg.exec();
+
+        // Do not finally accept until we do not have any errors
+        return;
+    }
     QDialog::accept();
 }
 void AddDuctDialog::reject(){
@@ -194,14 +207,14 @@ void AddDuctDialog::changed(){
     
     PreviewShow pshow = (PreviewShow) _dialog->previewshow->currentIndex();
     
-    std::unique_ptr<Geom> geom;
+    std::unique_ptr<Duct> duct;
     try {
-        geom = std::unique_ptr<Geom>(new Geom(_duct));
+        duct = std::unique_ptr<Duct>(new Duct(0,_duct));
     }
     catch(TaSMETError& e) {
         return;
     }
-    vd x = geom->x;
+    vd x = duct->x;
     vd y;
 
     switch (pshow) {

src/gui/mainwindow.cpp

@@ -31,11 +31,12 @@ TaSMETMainWindow::TaSMETMainWindow():
     restoreGeometry(settings.value("geometry").toByteArray());
     restoreState(settings.value("windowState").toByteArray());        
         
+
     for(const SystemType& t: SystemType_vec){
         window->systemtype->addItem(SystemTypeToString(t));
     }
-    for(const GasType& t: GasType_vec){
+    for(int gastype = pb::GasType_MIN;gastype<=pb::GasType_MAX;gastype++){
-        window->gastype->addItem(GasTypeToString(t));
+        window->gastype->addItem(QString::fromStdString(GasType_Name((pb::GasType) gastype)));
     }
     for(const SegmentType& t: SegmentType_vec){
         window->segmenttype->addItem(SegmentTypeToString(t));
@@ -81,8 +82,9 @@ void TaSMETMainWindow::on_addsegment_clicked() {
         dialog.set(ductmap[id]);
     }
 
-    int code = dialog.exec();
+    int exitcode = dialog.exec();
-    if(code == QDialog::Accepted) {
+    
+    if(exitcode == QDialog::Accepted) {
 
         VARTRACE(15,dialog.get().name());
         ductmap[id] = dialog.get();

src/gui/solver.ui

@@ -6,80 +6,139 @@
    <rect>
     <x>0</x>
     <y>0</y>
-    <width>747</width>
+    <width>711</width>
-    <height>336</height>
+    <height>494</height>
    </rect>
   </property>
   <property name="windowTitle">
    <string>Dialog</string>
   </property>
-  <widget class="QGroupBox" name="solversettings">
+  <widget class="QWidget" name="">
    <property name="geometry">
     <rect>
      <x>10</x>
      <y>40</y>
-     <width>306</width>
+     <width>651</width>
-     <height>168</height>
+     <height>341</height>
     </rect>
    </property>
-   <property name="title">
+   <layout class="QVBoxLayout" name="verticalLayout_3">
-    <string>Solver Settings</string>
+    <item>
-   </property>
+     <layout class="QHBoxLayout" name="horizontalLayout">
-   <layout class="QGridLayout" name="gridLayout_4">
+      <item>
-    <item row="3" column="0">
+       <widget class="QGroupBox" name="solversettings">
-     <widget class="QLabel" name="label_14">
+        <property name="title">
-      <property name="text">
+         <string>Solver Settings</string>
-       <string>Maximum iterations:</string>
+        </property>
-      </property>
+        <layout class="QGridLayout" name="gridLayout_4">
-     </widget>
+         <item row="3" column="0">
+          <widget class="QLabel" name="label_14">
+           <property name="text">
+            <string>Maximum iterations:</string>
+           </property>
+          </widget>
+         </item>
+         <item row="4" column="0">
+          <widget class="QLabel" name="label_13">
+           <property name="text">
+            <string>Solver type</string>
+           </property>
+          </widget>
+         </item>
+         <item row="4" column="1">
+          <widget class="QComboBox" name="solvertype"/>
+         </item>
+         <item row="2" column="1">
+          <widget class="QLineEdit" name="reltol">
+           <property name="alignment">
+            <set>Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter</set>
+           </property>
+          </widget>
+         </item>
+         <item row="2" column="0">
+          <widget class="QLabel" name="label_12">
+           <property name="text">
+            <string>Solution tolerance:</string>
+           </property>
+          </widget>
+         </item>
+         <item row="1" column="0">
+          <widget class="QLabel" name="label_8">
+           <property name="text">
+            <string>Residual tolerance:</string>
+           </property>
+          </widget>
+         </item>
+         <item row="1" column="1">
+          <widget class="QLineEdit" name="funtol">
+           <property name="alignment">
+            <set>Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter</set>
+           </property>
+          </widget>
+         </item>
+         <item row="3" column="1">
+          <widget class="QLineEdit" name="maxiter">
+           <property name="alignment">
+            <set>Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter</set>
+           </property>
+          </widget>
+         </item>
+        </layout>
+       </widget>
+      </item>
+      <item>
+       <widget class="QGroupBox" name="groupBox">
+        <property name="title">
+         <string>Solve!</string>
+        </property>
+        <layout class="QVBoxLayout" name="verticalLayout">
+         <item>
+          <widget class="QPushButton" name="pushButton">
+           <property name="text">
+            <string>Solve</string>
+           </property>
+          </widget>
+         </item>
+         <item>
+          <widget class="QPushButton" name="pushButton_2">
+           <property name="text">
+            <string>Stop</string>
+           </property>
+          </widget>
+         </item>
+         <item>
+          <widget class="QPushButton" name="pushButton_3">
+           <property name="text">
+            <string>Single iteration</string>
+           </property>
+          </widget>
+         </item>
+        </layout>
+       </widget>
+      </item>
+     </layout>
     </item>
-    <item row="4" column="0">
+    <item>
-     <widget class="QLabel" name="label_13">
+     <widget class="QGroupBox" name="groupBox_2">
-      <property name="text">
+      <property name="title">
-       <string>Solver type</string>
+       <string>Progress</string>
-      </property>
-     </widget>
-    </item>
-    <item row="4" column="1">
-     <widget class="QComboBox" name="solvertype"/>
-    </item>
-    <item row="2" column="1">
-     <widget class="QLineEdit" name="reltol">
-      <property name="alignment">
-       <set>Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter</set>
-      </property>
-     </widget>
-    </item>
-    <item row="2" column="0">
-     <widget class="QLabel" name="label_12">
-      <property name="text">
-       <string>Solution tolerance:</string>
-      </property>
-     </widget>
-    </item>
-    <item row="1" column="0">
-     <widget class="QLabel" name="label_8">
-      <property name="text">
-       <string>Residual tolerance:</string>
-      </property>
-     </widget>
-    </item>
-    <item row="1" column="1">
-     <widget class="QLineEdit" name="funtol">
-      <property name="alignment">
-       <set>Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter</set>
-      </property>
-     </widget>
-    </item>
-    <item row="3" column="1">
-     <widget class="QLineEdit" name="maxiter">
-      <property name="alignment">
-       <set>Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter</set>
       </property>
+      <layout class="QVBoxLayout" name="verticalLayout_2">
+       <item>
+        <widget class="QWidget" name="rel_error" native="true"/>
+       </item>
+       <item>
+        <widget class="QWidget" name="fun_error" native="true"/>
+       </item>
+      </layout>
      </widget>
     </item>
    </layout>
   </widget>
+  <zorder>solversettings</zorder>
+  <zorder>groupBox</zorder>
+  <zorder>groupBox_2</zorder>
+  <zorder>solversettings</zorder>
  </widget>
  <resources/>
  <connections/>

src/material/gas.h

@@ -14,13 +14,20 @@
 #include "tasmet_types.h"
 #include "tasmet_constants.h"
 
+#include "protobuf/gas.pb.h"
+
+using pb::GasType;
+using pb::air;
+using pb::helium;
+using pb::nitrogen;
+
 #define element_wise(varname) \
     vd varname_(T.size());     \
     for(us i=0;i<T.size();i++) \
         varname_ = varname(T(i),p(i));\
     return varname_
 
-DECLARE_ENUM(GasType,air,helium,nitrogen);
+
 
 class Gas{
 public:

src/protobuf/gas.proto

@@ -0,0 +1,9 @@
+
+syntax = "proto2";
+package pb;
+
+enum GasType {
+     air = 0;
+     helium = 1;
+     nitrogen = 2;
+}

src/protobuf/system.proto

@@ -2,15 +2,11 @@ syntax = "proto2";
 package pb;
 
 import "duct.proto";
+import "gas.proto";
 
 enum SystemType {
     TaSystem = 0;
 }
-enum GasType {
-     air = 0;
-     helium = 1;
-     nitrogen = 2;
-}
 message System {
     required SystemType systemtype = 1;
     required uint32 Nf = 2;
@@ -18,7 +14,7 @@ message System {
     required GasType gastype = 4;
     required double p0 = 5;
     required double T0 = 6;
-
     map<uint32,Duct> ducts = 7;
+    repeated double solution = 8;
 
 }

src/solver/solver.cpp

@@ -21,7 +21,7 @@ template<typename system_T,typename result_T>
 Solver<system_T,result_T>::Solver(const system_T& sys){
     _sys = sys.copy();
     if(!_sys)
-        throw TasMETBadAlloc();
+        throw TaSMETBadAlloc();
     _running = false;
 }
 
@@ -51,7 +51,7 @@ void Solver<system_T,result_T>::start(progress_callback* callback,bool wait){
                                                _sys,
                                                callback);
         if(!_solver_thread)
-            throw TasMETBadAlloc();
+            throw TaSMETBadAlloc();
     }
     else {
 

src/sys/segment.h

@@ -23,18 +23,23 @@ class TaSystem;
 
 class Segment{
 
-    // User identifier
+protected:
+    // ID
+    us _id;
+    // Name
     std::string _name;
 
-protected:
+
-    Segment(const std::string& name): _name(name) {}
+    Segment(const us id,const std::string& name): _id(id),_name(name) {}
-    Segment(const Segment& o): Segment(o._name){}
+    Segment(const Segment& o): Segment(o._id,o._name){}
     Segment& operator=(const Segment&)=delete;
 public:
     virtual ~Segment(){}
 
     virtual Segment* copy() const = 0;
     
+    virtual vd initialSolution(const TaSystem&) const = 0;
+
     // Get and set name
     const std::string& getName() const{return _name;} // This one is just the name
     void setName(const std::string& name){ _name = name; } // This one is just the name
@@ -44,13 +49,14 @@ public:
     // does, the derived class should return which equation should be
     // overwritten with the mass arbitration equation.
     virtual int arbitrateMassEq() const {return -1;}
-    virtual void residual(const TaSystem&,vd&,const us insertion_start) const=0;
+    virtual void residual(const TaSystem&,
+                          arma::subview_col<d>&& residual // Here we store the residual
+                          ) const=0;
 
-    virtual void updateSolution(const TaSystem&,const vd&) = 0;
-    virtual getSolution(const TaSystem&,vd& sol,const us insertion_start) const = 0;
 
     // Return the total number of equations in this segment
     virtual us getNEqs(const TaSystem&) const { return 0;}
+    
     // Return the total number of DOFS in this segment
     virtual us getNDofs(const TaSystem&) const { return 0;}
 

src/sys/tasystem.cpp

@@ -12,12 +12,44 @@
 #include "tasmet_assert.h"
 #include "tasmet_exception.h"
 #include "tasmet_constants.h"
+#include "duct.h"
 
-TaSystem::TaSystem(const GlobalConf& gc,const Gas& g):
+TaSystem::TaSystem(const pb::System& sys):
-    GlobalConf(gc)),
+    GlobalConf(sys.nf(),sys.freq())
-    _gas(g.copy())
 {
-    TRACE(14,"TaSystem::TaSystem(gc,gastype)");
+    TRACE(14,"TaSystem::TaSystem()");
+
+    if(sys.systemtype() != pb::TaSystem) {
+        throw TaSMETError("Invalid system type for TaSystem");
+    }
+
+    // Checking parameters T0 and p0 also happens in newGas method
+    _gas = std::unique_ptr<Gas>(Gas::newGas(sys.gastype(),sys.t0(),sys.p0()));
+    if(!_gas) throw  TaSMETBadAlloc();
+
+    // Create all ducts
+    for(const auto& d : sys.ducts()) {
+        try {
+            _segs[d.first] = new Duct(d.first,d.second);
+            if(!_segs[d.first]) throw TaSMETBadAlloc();
+        }
+        catch(TaSMETError e) {
+            // Cleanup the already successfully created Ducts
+            cleanup();
+        }
+    }
+    
+    // Copy solution vector, if valid
+    const auto& sol = sys.solution();
+    us size = sol.size(), i=0;
+    if(size>0) {
+        _solution = vd(size);
+        for(auto& val: sol) {
+            _solution(i) = val;
+            i++;
+        }
+    }
+
 }
 TaSystem::TaSystem(const TaSystem& o):
     GlobalConf(o),                 // Share a ptr to the Global conf
@@ -30,8 +62,8 @@ TaSystem::TaSystem(const TaSystem& o):
 
     for(auto& seg: _segs){
         seg.second = seg.second->copy();
+        if(!seg.second) throw TaSMETBadAlloc();
     }
-
 }
 int TaSystem::getArbitrateMassEq(const vus& neqs) const {
     // Tells the TaSystem which Dof should be overwritten with the
@@ -68,79 +100,21 @@ int TaSystem::getArbitrateMassEq(const vus& neqs) const {
     }
     return arbitrateMassEq;
 }
-TaSystem& TaSystem::add(const us id,const Segment& s){
-
-    // TRACE(24,"TaSystem::add(id,Seg)");
-
-    // if(_segs.find(s.getid())!=_segs.end()){
-    //     std::stringstream error;
-    //     error << "Segment with id " << s.getid() <<
-    //         "already present in the system";
-    //     throw TaSMETError(error);
-    // }
-    _segs[id]=s.copy();
-    return *this;
-}
-void TaSystem::updateSolution(const vd& sol) {
-    TRACE(15,"TaSystem::updateSolution()");
-
-    us firstdof = 0;
-    us ndofs;
-    Segment* seg;
-    for(auto& seg_: _segs) {
-        seg = seg_.second;
-        ndofs = seg->getNDofs(*this);
-
-        firstdof += ndofs;
-
-        seg->updateSolution(*this,sol.subvec(firstdof,firstdof+ndofs-1));
-    }
-}
-vd TaSystem::getSolution() const {
-    us firstdof = 0,i=0;
-    vus ndofs = getNDofs();
-    us total_ndofs = arma::sum(ndofs);
-    vd sol(total_ndofs);
-
-    #ifdef TASMET_DEBUG
-    sol.zeros();
-    #endif
-
-    Segment* seg;
-    for(auto& seg_: _segs) {
-        seg = seg_.second;
-        seg->getSolution(*this,sol,firstdof);
-        firstdof += ndofs(i);
-        i++;
-    }
-    return sol;
-}
 vd TaSystem::residual() const {
     TRACE(15,"TaSystem::residual()");
 
     vus neqs = getNEqs();
     us total_neqs = arma::sum(neqs);
-    vus eqstart(neqs.size());
-
-    us i=0;
-
-    for(us& eqs : neqs) {
-
-        if(i>0) {
-            eqstart(i) = eqstart(i-1) + eqs;
-        }
-        else {
-            eqstart(i) = 0;
-        }
-        i++;                 
-    }
-    assert(i==_segs.size()-1);
 
     if(total_neqs>constants::maxndofs)      {
         throw TaSMETError("Too many DOFS required."
                           " Problem too large.");
     }
 
+    // This vector of indices stores the last equation number + 1 for
+    // each equation set in a Segment
+    vus eqsend = arma::cumsum(neqs);
+
     int arbitrateMassEq = getArbitrateMassEq(neqs);
 
     // This is the mass in the sytem. Only counted when there is an
@@ -148,18 +122,23 @@ vd TaSystem::residual() const {
     d mass=0;
     
     VARTRACE(25,total_neqs);
-    VARTRACE(25,eqstart);
+    VARTRACE(25,eqsend);
 
     vd residual(total_neqs);
 
-    i=0;
+    us i=0;
     const Segment* seg;
     for(auto seg_: _segs) {
 
         seg = seg_.second;
 
-        // Put the residual of the segment in the over-all residual
+        if(i==0) {
-        seg->residual(*this,residual,eqstart(i));
+            // Put the residual of the segment in the over-all residual
+            seg->residual(*this,residual.subvec(0,eqsend(0)-1));
+        }
+        else {
+            seg->residual(*this,residual.subvec(eqsend(i-1),eqsend(i)-1));
+        }
 
         // Count the mass, add it
         if(arbitrateMassEq!=-1) {
@@ -180,7 +159,51 @@ vd TaSystem::residual() const {
 
     return residual;
 }
+vd TaSystem::getSolution() const {
 
+    if(_solution.size() == 0) {
+        // Create the initial solution from the segments
+        // and store it here
+
+        vus ndofs = getNDofs();
+        vus dofend = arma::cumsum(ndofs);
+        us total_dofs = arma::sum(ndofs);
+        vd solution = vd(total_dofs);
+
+        us i=0;
+        const Segment* seg;
+        for(auto& seg_: _segs) {
+            seg = seg_.second;
+
+            if(ndofs(i)>0) {
+                if(i==0) {
+                    solution.subvec(0,ndofs(0)-1) =
+                        seg->initialSolution(*this);
+                }
+                else {
+                    solution.subvec(dofend(i-1),dofend(i)-1) =
+                        seg->initialSolution(*this);
+                }
+                i++;
+            }
+        }
+        return solution;
+    } // if the solution did not yet exist
+
+    return _solution;
+}
+const arma::subview_col<d> TaSystem::getSolution(const us seg_id) const {
+
+    vus ndofs = getNDofs();
+    vus dofsend = arma::cumsum(ndofs);
+
+    if(seg_id == 0) {
+        return _solution.subvec(0,dofsend(0)-1);
+    }
+    else {
+        return _solution.subvec(dofsend(seg_id-1),dofsend(seg_id)-1);
+    }
+}
 vus TaSystem::getNDofs() const  {
     TRACE(0,"TaSystem::getNDofs()");
     vus Ndofs(_segs.size());
@@ -215,7 +238,6 @@ void TaSystem::show(us detailnr){
         }
     } // detailnr>0
 }
-
 TripletList TaSystem::jacTriplets() const {
 
     TRACE(14,"TaSystem::jacobian()");
@@ -286,12 +308,14 @@ dmat TaSystem::showJac(){
 
 TaSystem::~TaSystem() {
     TRACE(25,"~TaSystem()");
+    cleanup();
+}
+void TaSystem::cleanup() {
 
     for(auto& seg: _segs){
         delete seg.second;
     }
-    delete _gas;
+
 }
 
-
 //////////////////////////////////////////////////////////////////////

src/sys/tasystem.h

@@ -9,11 +9,14 @@
 #ifndef _TASYSTEM_H_
 #define _TASYSTEM_H_
 #include "system.h"
-#include "segment.h"
 #include "globalconf.h"
 #include "triplets.h"
 #include "gas.h"
+#include "protobuf/system.pb.h"
 #include <map>
+#include <memory>
+
+class Segment;
 
 // Inherit all global configuration members
 class TaSystem : public GradientNonlinearSystem, public GlobalConf {
@@ -23,17 +26,17 @@ protected:
     
     std::map<us,Segment*> _segs;		
 
-    Gas* _gas = nullptr;
+    std::unique_ptr<Gas> _gas;
+
+    vd _solution;
+    vus _startdof;              // Store the start DOFS for each segment
     
     TaSystem& operator=(const TaSystem& other)=delete;
     TaSystem(const TaSystem& o);
 public:
-    TaSystem(const GlobalConf& gc,const Gas& g);
+    TaSystem(const pb::System&);
 
-    // Set globalconf configuration. Applies updateNf as well.
+    const Gas& gas() const {return *_gas;}
-    void setGc(const GlobalConf& gc);
-    const Gas& getGas() const {return *_gas;}
-    void setGas(const Gas& g);
     
     // Set and get the mass in the system. If the mass is not set
     // before initializing, the mass is computed from the segment's
@@ -46,17 +49,17 @@ public:
 
     us nSegments() const {return _segs.size();}
 
-    TaSystem& add(const us id,const Segment& s);	// Add a segment to the
-    // system. It creates a copy
-
     dmat showJac();
     virtual void show(us detailnr=0);
 
     vd residual() const; 
 
     vd getSolution() const;
+
+    // Obtain the solution vector for the Segment with given id
+    const arma::subview_col<d> getSolution(const us seg_id) const;
     
-    virtual void updateSolution(const vd& resvec);	// Update the solution
+    virtual void updateSolution(const vd& sol) {_solution = sol; }	// Update the solution
 
     // Compute Jacobian matrix. The dampfac value is used in an
     // EngineSystem
@@ -79,7 +82,8 @@ public:
 protected:
     virtual int getArbitrateMassEq(const vus& neqs) const;
     virtual TripletList jacTriplets() const;
-
+private:
+    void cleanup();
 }; // class System
   
 

src/tasmet_constants.h

@@ -10,11 +10,6 @@
 #include "tasmet_enum.h"
 #include "tasmet_types.h"
 
-template<typename T>
-T max(T& t1,T& t2) { return t1>t2?t1:t2;}
-template<typename T>
-T min(T& t1,T& t2) { return t1>t2?t2:t1;}
-
 // Variables and their names
 // Unfortunately to let the code compile with Swig v 2.0, strongly
 // typed enums are not supported. Therefore this is a normal

src/tasmet_exception.cpp

@@ -8,7 +8,14 @@
 
 #include "tasmet_exception.h"
 
-const char*  TasMETBadAlloc::what() const throw() {
+const char* TaSMETError::what() const throw() {
+    return _msg.c_str();
+}
+void TaSMETError::setContext(const std::string& ctx) {
+    std::string oldmsg = _msg;
+    _msg = ctx + ": " + oldmsg;
+}
+const char*  TaSMETBadAlloc::what() const throw() {
 	return "Error: memory allocation failed. "
    		"Please make sure enough memory is available and restart the application";
 }

src/tasmet_exception.h

@@ -14,12 +14,23 @@
 #include <sstream>              // stringstream
 
 class TaSMETError : public std::runtime_error {
+    std::string _msg;
 public:
-    TaSMETError(const std::string& msg = "") : std::runtime_error(msg) {}
+    TaSMETError(const std::string& msg = "") :
-    TaSMETError(const std::stringstream& stream) : std::runtime_error(stream.str()){}
+        std::runtime_error(""),
+        _msg(msg)
+    {}
+    TaSMETError(const std::stringstream& stream) :
+        std::runtime_error(""),
+        _msg(stream.str())
+    {
+        
+    }
+    void setContext(const std::string& ctx);
+    virtual const char* what() const throw ();
 };
 
-class TasMETBadAlloc: public std::bad_alloc {
+class TaSMETBadAlloc: public std::bad_alloc {
 	virtual const char* what() const throw();
 };
 

src/tasmet_pyeval.cpp

@@ -15,22 +15,25 @@
 #include <QVariantList>
 #include <sstream>
 
-inline void getError(PythonQt* pyqt) {
+inline void getError(PythonQt* pyqt,string err_msg) {
     if(pyqt->hadError()) {
         pyqt->handleError();
         pyqt->clearError();
 
-        throw TaSMETError("Script error");
+        throw TaSMETError(err_msg);
     }
 }
 
 
-EvaluateFun::EvaluateFun(const string& fun_return) {
+EvaluateFun::EvaluateFun(const string& fun_return,
+                         const string& err_msg):
+    _err_msg(err_msg)
+{
 
     _pyqt = PythonQt::self();
 
     try {
-        getError(_pyqt);
+        getError(_pyqt,_err_msg);
     }
     catch(TaSMETError& e){
         cerr << "Uncleared error in script" << endl;
@@ -46,12 +49,12 @@ EvaluateFun::EvaluateFun(const string& fun_return) {
     context.evalScript(QString::fromStdString(script.str()));
 
     // See if we have an error. If so, throw
-    getError(_pyqt);
+    getError(_pyqt,_err_msg);
 
 }
 EvaluateFun::~EvaluateFun() {
     try {
-        getError(_pyqt);
+        getError(_pyqt,_err_msg);
     }
     catch(...) {
 
@@ -72,7 +75,7 @@ void EvaluateFun::addGlobalDef(const string& name,const d value) {
         
     context.evalScript(QString::fromStdString(script.str()));        
 
-    getError(_pyqt);
+    getError(_pyqt,_err_msg);
 
 }
 vd EvaluateFun::operator()(const vd& x) {
@@ -91,7 +94,7 @@ vd EvaluateFun::operator()(const vd& x) {
         res = context.call("myfun",args);
         y(i) = res.toDouble();
             
-        getError(_pyqt);
+        getError(_pyqt,_err_msg);
     }
     return y;
 }

src/tasmet_pyeval.h

@@ -16,8 +16,9 @@ class PythonQt;
 class EvaluateFun {
 
     PythonQt* _pyqt;
+    string _err_msg;
 public:
-    EvaluateFun(const string& fun_return);
+    EvaluateFun(const string& fun_return,const string& err_msg = "Script error");
     
     // Add a global definition to the namespace
     void addGlobalDef(const string& name,