diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 03d3806..9e601a8 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -7,7 +7,7 @@ include_directories(
   ${ARMADILLO_INCLUDE_DIRS}
   .
   protobuf
-  # duct
+  duct
   # duct/cell
   # duct/connectors
   # duct/eq
@@ -29,8 +29,10 @@ add_library(tasmet_src
   tasmet_tracer.cpp
   tasmet_exception.cpp
   tasmet_assert.cpp
+  tasmet_pyeval.cpp
 
   duct/grid.cpp
+  duct/geom.cpp
 
   funcs/bessel.cpp
   funcs/cbessj.cpp
diff --git a/src/duct/duct.h b/src/duct/duct.h
new file mode 100644
index 0000000..3e3d3dd
--- /dev/null
+++ b/src/duct/duct.h
@@ -0,0 +1,34 @@
+// duct.h
+//
+// Author: J.A. de Jong
+//
+// Description:
+//
+//////////////////////////////////////////////////////////////////////
+#pragma once
+#ifndef DUCT_H
+#define DUCT_H
+#include "segment.h"
+#include "duct.pb.h"
+
+class Equation;
+class Drag;
+class Heat;
+
+class Duct {
+    Duct(const string& name,const Geom& geom);
+public:
+
+    const Geom& geom() const;
+
+
+    // Parsing a protobuf to generate a NEW Duct
+    static Duct* parseProto(const pb::Duct&);
+
+
+    // Solving
+
+};
+
+#endif // DUCT_H
+//////////////////////////////////////////////////////////////////////
diff --git a/src/duct/geom.cpp b/src/duct/geom.cpp
new file mode 100644
index 0000000..fd00588
--- /dev/null
+++ b/src/duct/geom.cpp
@@ -0,0 +1,101 @@
+// geom.cpp
+//
+// last-edit-by: J.A. de Jong
+//
+// Description:
+//
+//////////////////////////////////////////////////////////////////////
+#include "geom.h"
+#include "grid.h"
+#include "duct.pb.h"
+#include "tasmet_pyeval.h"
+#include <memory>
+#include "tasmet_exception.h"
+
+Geom::Geom(const pb::Duct& duct) {
+
+    // Step one: create the grid
+    std::unique_ptr<Grid> grid;
+    switch (duct.gridtype()) {
+    case pb::Linear:
+        grid = std::unique_ptr<Grid>(new LinearGrid(duct.ngp(),
+                                                    duct.length()));
+        break;
+    case pb::BlGrid:
+            grid = std::unique_ptr<Grid>(new BlGrid(duct.length(),
+                                                    duct.dxb(),
+                                                    duct.dxmid()));
+        break;
+    default:
+        return;
+        break;
+    }
+
+    // Store x
+    x = grid->getx();
+
+    EvaluateFun Sfun(duct.area());
+    Sfun.addGlobalDef("L",duct.length());
+    S = Sfun(x);
+
+    for(auto& Si:S) {
+        if(Si <=0){
+            throw TaSMETError("Invalid cross-sectional area function");
+        }
+    }
+
+
+    EvaluateFun phifun(duct.phi());
+    Sfun.addGlobalDef("L",duct.length());
+    phi = Sfun(x);
+
+    EvaluateFun rhfun(duct.rh());
+    Sfun.addGlobalDef("L",duct.length());
+    rh = Sfun(x);
+
+    cshape = duct.cshape();
+}
+string Geom::shapeString() const {
+    return pb::Cshape_Name(cshape);
+}
+
+bool Geom::isPrismatic() const {
+    bool ok = true;
+    d S0 = S(0);
+    for(auto& Si : S) {
+        ok&= Si==S0;
+
+    }
+    return ok;
+}
+d Geom::getFluidVolume() const {
+    d Vf=0;
+    us i=0;
+
+    Vf+=0.5*(x(i+1)-x(i))*Sf(i); // First 'cell'
+
+    for(i=1;i<ngp()-1;i++){    // Middle cells
+        Vf+=0.5*(x(i)-x(i-1))*Sf(i);
+        Vf+=0.5*(x(i+1)-x(i))*Sf(i);
+    }
+
+    Vf+=0.5*(x(i)-x(i-1))*Sf(i); // Last 'cell'
+    return Vf;
+  }
+  d Geom::getSolidVolume() const {
+    d Vs=0;
+    int i=0;
+
+    Vs+=0.5*(x(i+1)-x(i))*Ss(i); // First 'cell'
+
+    for(int i=1;i<ngp()-1;i++){    // Middle cells
+        Vs+=0.5*(x(i)-x(i-1))*Ss(i);
+        Vs+=0.5*(x(i+1)-x(i))*Ss(i);
+    }
+
+    Vs+=0.5*(x(i)-x(i-1))*Ss(i); // Last 'cell'
+    return Vs;
+  }
+
+
+//////////////////////////////////////////////////////////////////////
diff --git a/src/duct/geom.h b/src/duct/geom.h
new file mode 100644
index 0000000..24ebe7e
--- /dev/null
+++ b/src/duct/geom.h
@@ -0,0 +1,44 @@
+// geom.h
+//
+// Author: J.A. de Jong
+//
+// Description:
+//
+//////////////////////////////////////////////////////////////////////
+#pragma once
+#ifndef GEOM_H
+#define GEOM_H
+
+#include "tasmet_constants.h"
+#include "tasmet_types.h"
+#include "duct.pb.h"
+
+class Geom{
+
+public:
+    // Parse geometry from duct
+    Geom(const pb::Duct& duct);
+
+    // Public members
+    vd x,S,phi,rh;
+    pb::Cshape cshape;
+
+    // Shape keyword
+    string shapeString() const;
+
+    bool isPrismatic() const;
+
+    us ngp() const {return x.size();}
+    d L() const {return x(x.size()-1);}
+    d Sleft() const {return S(0);}
+    d Sright() const {return S(ngp()-1);}
+
+    d Sf(us i) const {return S(i)*phi(i);}       // Fluid-occupied cross-sectional area of cell-wall
+    d Ss(us i) const {return (1-phi(i))*S(i);}       // Solid-occupied cross-sectional area
+
+    d getFluidVolume() const;
+    d getSolidVolume() const;
+
+};                            /* class Geom */
+
+#endif /* _GEOM_H_ */
diff --git a/src/duct/grid.cpp b/src/duct/grid.cpp
index 1ff7a5d..16c5fe7 100644
--- a/src/duct/grid.cpp
+++ b/src/duct/grid.cpp
@@ -20,15 +20,15 @@ LinearGrid::LinearGrid(us ngp,d L):
         throw TaSMETError(error);
     }
     if(L<=0){
-        throw TaSMETError("Illegal length chosen.");
+        throw TaSMETError("LinearGrid: illegal length chosen.");
     }
 }
 BlGrid::BlGrid(d L,d dxb,d dxmid):
     L(L),dxb(dxb),dxmid(dxmid)
 {
     TRACE(15,"BlGrid::BlGrid()");
-    if(L<=0 || dxb<=0 || dxmid < dxb){
-        throw TaSMETError("Illegal length chosen.");
+    if(L<=0 || dxb<=0 || dxmid <= dxb){
+        throw TaSMETError("BlGrid: illegal length chosen.");
     }
     
 }
@@ -45,6 +45,5 @@ vd BlGrid::getx() const {
     vd xi=linspace(0,1,N);
     
     return L*0.5*(1+tanh(delta*(xi-0.5))/tanh(delta/2));
-
     
 }
diff --git a/src/duct/grid.h b/src/duct/grid.h
index ce2369e..bba4d7b 100644
--- a/src/duct/grid.h
+++ b/src/duct/grid.h
@@ -14,7 +14,7 @@ public:
 
 class LinearGrid: public Grid{
     us ngp;
-    d L;			// Length of the Duct
+    d L; // Length of the Duct
 public:
     LinearGrid(us ngp,d L);
     LinearGrid(const LinearGrid& g):LinearGrid(g.ngp,g.L){}
@@ -26,7 +26,7 @@ public:
 // dxb: boundary layer grid spacing (minimal grid spacing)
 // xb: boundary layer tickness
 // dxmid: spacing in the middle part
-  
+
 
 class BlGrid:public Grid{
     d L,dxb,dxmid;
diff --git a/src/gui/CMakeLists.txt b/src/gui/CMakeLists.txt
index af70fe2..d833a74 100644
--- a/src/gui/CMakeLists.txt
+++ b/src/gui/CMakeLists.txt
@@ -1,6 +1,7 @@
 #================
 # The TaSMET GUI
 #================
+
 add_library(tasmet_gui
   mainwindow.cpp
   add_duct_dialog.cpp
diff --git a/src/gui/add_duct_dialog.cpp b/src/gui/add_duct_dialog.cpp
index 0ba9459..701b45f 100644
--- a/src/gui/add_duct_dialog.cpp
+++ b/src/gui/add_duct_dialog.cpp
@@ -7,7 +7,7 @@
 //////////////////////////////////////////////////////////////////////
 #include "ui_add_duct_dialog.h"
 #include "add_duct_dialog.h"
-#include <PythonQt.h>
+
 #include <QSignalBlocker>
 #include <qcustomplot.h>
 #include "tasmet_constants.h"
@@ -18,6 +18,7 @@
 #include "duct/grid.h"
 #include <QVector>
 #include "tasmet_qt.h"
+#include "geom.h"
 
 DECLARE_ENUM(PreviewShow,CSArea,Porosity,HydraulicRadius,SolidTemperatureFunction)
 
@@ -148,59 +149,6 @@ void AddDuctDialog::set(const pb::Duct& duct) {
     _dialog->stempfunc->setText(QString::fromStdString(duct.stempfunc()));
 }
 
-class PyEvaluate {
-    PythonQtObjectPtr _context;
-    PythonQt* _pyqt;
-public:
-    PyEvaluate(const std::string& fun_return) {
-
-        _pyqt = PythonQt::self();
-        if(_pyqt->hadError()) {
-            TRACE(15,"Previous error in script");
-            _pyqt->handleError();
-            _pyqt->clearError();
-        }
-
-        _context = _pyqt->getMainModule();
-        
-
-        std::stringstream script;
-        // script << "print(\"hoi\")\n";
-        script << "def myfun(x,L):\n    return ";
-        script << fun_return << "\n";
-
-        _context.evalScript(QString::fromStdString(script.str()));
-
-        if(_pyqt->hadError()) {
-            _pyqt->clearError();
-            throw TaSMETError("Script error");
-        }
-
-    }
-    vd operator()(const d L,const vd& x) {
-
-        vd y(x.size());
-
-        QVariant res;
-
-        for(us i=0;i<x.size();i++) {
-
-            QVariantList args;
-            args << x(i) << L;
-
-            res = _context.call("myfun",args);
-            y(i) = res.toDouble();
-            
-            if(_pyqt->hadError()) {
-                _pyqt->clearError();
-                throw TaSMETError("Script error");
-            }
-        }
-        return y;
-    }
-};
-
-
 
 void AddDuctDialog::changed(){
     TRACE(15,"AddDuctDialog::changed()");
@@ -226,79 +174,58 @@ void AddDuctDialog::changed(){
     // Empty the graph
     _plot->graph(0)->setData(QVector<d>(),QVector<d>());
 
-    std::unique_ptr<Grid> grid;
-
     switch (_duct.gridtype()) {
     case pb::Linear:
-        try {
-            grid = std::unique_ptr<Grid>(new LinearGrid(_duct.ngp(),_duct.length()));
-        }
-        catch(...) {
-            return;
-        }
         _dialog->ngp->setEnabled(true);        
         _dialog->dxb->setEnabled(false);
         _dialog->dxmid->setEnabled(false);
         break;
     case pb::BlGrid:
-        
-        try {
-            grid = std::unique_ptr<Grid>(new BlGrid(_duct.length(),
-                                                    _duct.dxb(),
-                                                    _duct.dxmid()));
-        }
-        catch(...) {
-            return;
-        }
-        _duct.clear_ngp();
         _dialog->ngp->setEnabled(false);        
         _dialog->dxb->setEnabled(true);
         _dialog->dxmid->setEnabled(true);
         break;
     default:
+        tasmet_assert(false,"Invalid grid type");
         return;
         break;
     }
 
-    vd x = grid->getx();
-    vd y;
-
+    
     PreviewShow pshow = (PreviewShow) _dialog->previewshow->currentIndex();
     
-    std::string pyeval_return_type;
+    std::unique_ptr<Geom> geom;
+    try {
+        geom = std::unique_ptr<Geom>(new Geom(_duct));
+    }
+    catch(TaSMETError& e) {
+        return;
+    }
+    vd x = geom->x;
+    vd y;
 
     switch (pshow) {
     case CSArea:
-        pyeval_return_type = _duct.area();
         _plot->yAxis->setLabel("S [m^2]");
+        y = geom->S;
         break;
     case Porosity:
-        pyeval_return_type = _duct.phi();
         _plot->yAxis->setLabel("phi [-]");
+        y = geom->phi;
         break;
     case HydraulicRadius:
-        pyeval_return_type = _duct.rh();
-        _plot->yAxis->setLabel("rh [m]");        
+        _plot->yAxis->setLabel("rh [m]");
+        y = geom->rh;
         break;
     case SolidTemperatureFunction:
-        pyeval_return_type = _duct.stempfunc();
-        _plot->yAxis->setLabel("Ts [K]");        
+        _plot->yAxis->setLabel("Ts [K]");
+        y = geom->rh;
         break;
     default:
         tasmet_assert(false,"Unhandled PreviewShow case");
         break;
     }
 
-    try {
-        PyEvaluate pyeval(pyeval_return_type);
-        y = pyeval(_duct.length(),x);
-
-    }
-    catch(TaSMETError& e) {
-        
-        return;
-    }
-
     QVector<d> qx = from_arma(x);
     QVector<d> qy = from_arma(y);
 
diff --git a/src/gui/add_duct_dialog.ui b/src/gui/add_duct_dialog.ui
index ce4e8d7..6a1df7d 100644
--- a/src/gui/add_duct_dialog.ui
+++ b/src/gui/add_duct_dialog.ui
@@ -7,7 +7,7 @@
     <x>0</x>
     <y>0</y>
     <width>841</width>
-    <height>876</height>
+    <height>654</height>
    </rect>
   </property>
   <property name="windowTitle">
diff --git a/src/protobuf/duct.proto b/src/protobuf/duct.proto
index 11e1045..c521e76 100644
--- a/src/protobuf/duct.proto
+++ b/src/protobuf/duct.proto
@@ -11,7 +11,8 @@ enum GridType {
 }
 enum HeatTransferModel {
      Isentropic = 0;
-     LaminarHeatTransfer = 1;
+     NoHeatTransfer = 1;     
+     LaminarHeatTransfer = 2;
 }
 enum DragModel {
      Inviscid = 0;
diff --git a/src/sys/segment.h b/src/sys/segment.h
index 22a72f1..cb95920 100644
--- a/src/sys/segment.h
+++ b/src/sys/segment.h
@@ -23,15 +23,12 @@ class TaSystem;
 
 class Segment{
 
-    us _id; // Unique number for each segment in a
-    // TaSystem
-
     // User identifier
     std::string _name;
 
 protected:
-    Segment(us id,const std::string& name): _id(id),_name(name) {}
-    Segment(const Segment& o): Segment(o._id,o._name){}
+    Segment(const std::string& name): _name(name) {}
+    Segment(const Segment& o): Segment(o._name){}
     Segment& operator=(const Segment&)=delete;
 public:
     virtual ~Segment(){}
@@ -40,9 +37,7 @@ public:
     
     // Get and set name
     const std::string& getName() const{return _name;} // This one is just the name
-    us getid() const{return _id;} // This one is just the name
     void setName(const std::string& name){ _name = name; } // This one is just the name
-    void setid(const us id){ _id=id;} // Set ID
     
     // Tell a TaSystem whether this Segment arbitrates Mass or
     // not. The special return value of -1 tells it does not. If it
diff --git a/src/sys/tasystem.cpp b/src/sys/tasystem.cpp
index b6d9ef2..2156a69 100644
--- a/src/sys/tasystem.cpp
+++ b/src/sys/tasystem.cpp
@@ -68,17 +68,17 @@ int TaSystem::getArbitrateMassEq(const vus& neqs) const {
     }
     return arbitrateMassEq;
 }
-TaSystem& TaSystem::operator+=(const Segment& s){
+TaSystem& TaSystem::add(const us id,const Segment& s){
 
-    TRACE(24,"TaSystem::operator+=(Seg)");
+    // 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[s.getid()]=s.copy();
+    // 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) {
diff --git a/src/sys/tasystem.h b/src/sys/tasystem.h
index bc75875..f20a737 100644
--- a/src/sys/tasystem.h
+++ b/src/sys/tasystem.h
@@ -48,7 +48,7 @@ public:
 
     us nSegments() const {return _segs.size();}
 
-    TaSystem& operator+=(const Segment& s);	// Add a segment to the
+    TaSystem& add(const us id,const Segment& s);	// Add a segment to the
     // system. It creates a copy
 
     dmat showJac();
diff --git a/src/tasmet_pyeval.cpp b/src/tasmet_pyeval.cpp
new file mode 100644
index 0000000..eecb449
--- /dev/null
+++ b/src/tasmet_pyeval.cpp
@@ -0,0 +1,100 @@
+// tasmet_pyeval.cpp
+//
+// last-edit-by: J.A. de Jong 
+// 
+// Description:
+//
+//////////////////////////////////////////////////////////////////////
+
+#include "tasmet_pyeval.h"
+#include <PythonQt.h>
+#include "tasmet_io.h"
+#include "tasmet_exception.h"
+#include <QString>
+#include <QVariant>
+#include <QVariantList>
+#include <sstream>
+
+inline void getError(PythonQt* pyqt) {
+    if(pyqt->hadError()) {
+        pyqt->handleError();
+        pyqt->clearError();
+
+        throw TaSMETError("Script error");
+    }
+}
+
+
+EvaluateFun::EvaluateFun(const string& fun_return) {
+
+    _pyqt = PythonQt::self();
+
+    try {
+        getError(_pyqt);
+    }
+    catch(TaSMETError& e){
+        cerr << "Uncleared error in script" << endl;
+    }
+
+    PythonQtObjectPtr context = _pyqt->getMainModule();
+
+    std::stringstream script;
+
+    script << "def myfun(x):\n    return ";
+    script << fun_return << "\n";
+
+    context.evalScript(QString::fromStdString(script.str()));
+
+    // See if we have an error. If so, throw
+    getError(_pyqt);
+
+}
+EvaluateFun::~EvaluateFun() {
+    try {
+        getError(_pyqt);
+    }
+    catch(...) {
+
+    }
+    PythonQtObjectPtr context = _pyqt->getMainModule();
+    string delfun =  "del myfun(x):\n ";
+
+    context.evalScript(QString::fromStdString(delfun));
+
+}
+void EvaluateFun::addGlobalDef(const string& name,const d value) {
+
+    PythonQtObjectPtr context = _pyqt->getMainModule();
+
+    std::stringstream script;
+    script << name << " = " ;
+    script << value << "\n";
+        
+    context.evalScript(QString::fromStdString(script.str()));        
+
+    getError(_pyqt);
+
+}
+vd EvaluateFun::operator()(const vd& x) {
+
+    PythonQtObjectPtr context = _pyqt->getMainModule();
+
+    vd y(x.size());
+
+    QVariant res;
+
+    for(us i=0;i<x.size();i++) {
+
+        QVariantList args;
+        args << x(i);
+
+        res = context.call("myfun",args);
+        y(i) = res.toDouble();
+            
+        getError(_pyqt);
+    }
+    return y;
+}
+
+
+//////////////////////////////////////////////////////////////////////
diff --git a/src/tasmet_pyeval.h b/src/tasmet_pyeval.h
new file mode 100644
index 0000000..9ede6fb
--- /dev/null
+++ b/src/tasmet_pyeval.h
@@ -0,0 +1,37 @@
+// tasmet_pyeval.h
+//
+// Author: J.A. de Jong 
+//
+// Description: This file implements a helper class to evaluate simple
+// 1D single-parameter python math functions
+//////////////////////////////////////////////////////////////////////
+#pragma once
+#ifndef TASMET_PYEVAL_H
+#define TASMET_PYEVAL_H
+#include "tasmet_types.h"
+
+class PythonQt;
+
+// Helper class to evaluate simple 1D math function evaluations in Python
+class EvaluateFun {
+
+    PythonQt* _pyqt;
+public:
+    EvaluateFun(const string& fun_return);
+    
+    // Add a global definition to the namespace
+    void addGlobalDef(const string& name,
+                      const d value);
+    
+    // Evaluate a single function at multiple points and return
+    // the result for each point
+    vd operator()(const vd& points);
+
+    // Used to cleanup the python namespace
+    ~EvaluateFun();
+};
+
+
+
+#endif // TASMET_PYEVAL_H
+//////////////////////////////////////////////////////////////////////