Mise a jour

Makefile

@@ -2,23 +2,50 @@ QT_FLAGS	=-fPIC `pkg-config --cflags --libs Qt5OpenGL`
 QT_LIB		=`pkg-config --libs Qt5OpenGL` -lGL -lGLU
 QT_MOC		= moc
 GPP		= g++
-FLAGS 		= -g -O3 -Isrc
+FLAGS 		= -g -O3 -Isrc -Wfatal-errors
 
-QT_FILES	= mainwindow input input_physics input_time input_geometry input_initial_state input_view view input_tank input_pump input_pump_tab input_cloud input_clouds_tab
-QT_MOC_FILES	= mainwindow input_physics input input_view input_geometry input_initial_state input_tank input_pump_tab input_pump input_cloud input_clouds_tab
-MATH_FILES	= poly3 algo spline
-SRC_FILES	= physics time geometry kernel initial_state source
+QT_INPUT_FILES	= main_window data input physics time geometry initial_state view tank pump pump_tab cloud clouds_tab
+QT_INPUT_MOC_FILES	= main_window physics input view geometry initial_state tank pump_tab pump cloud clouds_tab
+MATH_INPUT_FILES	= poly3 algo spline
+SRC_INPUT_FILES = physics time geometry initial_state source input_data
 
-QT_OBJS		= $(addprefix obj/qt/,$(addsuffix .o,$(QT_FILES)))
-QT_MOCS		= $(addprefix moc/,$(addsuffix .cpp,$(QT_MOC_FILES)))
-MATH_OBJS	= $(addprefix obj/math/,$(addsuffix .o,$(MATH_FILES)))
-SRC_OBJS	= $(addprefix obj/,$(addsuffix .o,$(SRC_FILES)))
-EXE		= RichardsFastSlow
+QT_OUTPUT_FILES	= main_window kernel output_view output
+QT_OUTPUT_MOC_FILES	= main_window kernel output
+MATH_OUTPUT_FILES	= algo
+SRC_OUTPUT_FILES = physics time geometry initial_state source input_data solution kernel
 
-all : $(EXE)
 
-$(EXE) : src/main.cpp $(MATH_OBJS) $(SRC_OBJS) $(QT_MOCS) $(QT_OBJS)
-	$(GPP) $(FLAGS) $(QT_FLAGS) $^ $(QT_LIB) -o $@
+QT_INPUT_OBJS	= $(addprefix obj/qt/input/,$(addsuffix .o,$(QT_INPUT_FILES)))
+QT_INPUT_MOCS	= $(addprefix moc/input/,$(addsuffix .cpp,$(QT_INPUT_MOC_FILES)))
+MATH_INPUT_OBJS	= $(addprefix obj/math/,$(addsuffix .o,$(MATH_INPUT_FILES)))
+SRC_INPUT_OBJS = $(addprefix obj/,$(addsuffix .o,$(SRC_INPUT_FILES)))
+
+QT_OUTPUT_OBJS	= $(addprefix obj/qt/output/,$(addsuffix .o,$(QT_OUTPUT_FILES)))
+QT_OUTPUT_MOCS	= $(addprefix moc/output/,$(addsuffix .cpp,$(QT_OUTPUT_MOC_FILES)))
+MATH_OUTPUT_OBJS	= $(addprefix obj/math/,$(addsuffix .o,$(MATH_OUTPUT_FILES)))
+SRC_OUTPUT_OBJS = $(addprefix obj/,$(addsuffix .o,$(SRC_OUTPUT_FILES)))
+
+MATH_KERNEL_OBJS	= $(addprefix obj/math/,$(addsuffix .o,$(MATH_OUTPUT_FILES)))
+SRC_KERNEL_OBJS = $(addprefix obj/,$(addsuffix .o,$(SRC_OUTPUT_FILES)))
+
+EXE_INPUT = RichardsFastSlowInput
+EXE_OUTPUT = RichardsFastSlowOutput
+EXE_KERNEL = RichardsFastSlowKernel
+
+all: $(EXE_OUTPUT)
+input: $(EXE_INPUT)
+output: $(EXE_OUTPUT)
+kernel: $(EXE_KERNEL)
+
+$(EXE_INPUT) : src/main_input.cpp $(SRC_INPUT_OBJS) $(MATH_INPUT_OBJS) $(QT_INPUT_MOCS) $(QT_INPUT_OBJS)
+			$(GPP) $(FLAGS) $(QT_FLAGS) $^ $(QT_LIB) -o $@
+
+$(EXE_OUTPUT) : src/main_output.cpp $(SRC_OUTPUT_OBJS) $(MATH_OUTPUT_OBJS) $(QT_OUTPUT_MOCS) $(QT_OUTPUT_OBJS)
+			$(GPP) $(FLAGS) $(QT_FLAGS) $^ $(QT_LIB) -o $@
+
+$(EXE_KERNEL) : src/main_kernel.cpp $(SRC_KERNEL_OBJS) $(MATH_KERNEL_OBJS)
+			$(GPP) $(FLAGS) $^ -o $@
+
 
 obj/qt/%.o : src/qt/%.cpp src/qt/%.hpp
 	$(GPP) $(FLAGS) $(QT_FLAGS) -c $< -o $@
@@ -33,4 +60,4 @@ moc/%.cpp : src/qt/%.hpp
 	$(QT_MOC) $< -o $@
 
 clean:
-	-$(RM) $(EXE) moc/* obj/*.o obj/math/*.o obj/qt/*.o src/*~ src/qt/*~ src/math/*~
+	-$(RM) $(EXE) moc/input/*.cpp moc/output/*.cpp moc/*.cpp obj/*.o obj/math/*.o obj/qt/input/*.o obj/qt/output/*.o obj/qt/*.o src/*~ src/qt/*~ src/math/*~

moc/input.cpp

@@ -1,186 +0,0 @@
-/****************************************************************************
-** Meta object code from reading C++ file 'input.hpp'
-**
-** Created by: The Qt Meta Object Compiler version 67 (Qt 5.15.2)
-**
-** WARNING! All changes made in this file will be lost!
-*****************************************************************************/
-
-#include <memory>
-#include "../src/qt/input.hpp"
-#include <QtCore/qbytearray.h>
-#include <QtCore/qmetatype.h>
-#if !defined(Q_MOC_OUTPUT_REVISION)
-#error "The header file 'input.hpp' doesn't include <QObject>."
-#elif Q_MOC_OUTPUT_REVISION != 67
-#error "This file was generated using the moc from 5.15.2. It"
-#error "cannot be used with the include files from this version of Qt."
-#error "(The moc has changed too much.)"
-#endif
-
-QT_BEGIN_MOC_NAMESPACE
-QT_WARNING_PUSH
-QT_WARNING_DISABLE_DEPRECATED
-struct qt_meta_stringdata_QtInput_t {
-    QByteArrayData data[12];
-    char stringdata0[116];
-};
-#define QT_MOC_LITERAL(idx, ofs, len) \
-    Q_STATIC_BYTE_ARRAY_DATA_HEADER_INITIALIZER_WITH_OFFSET(len, \
-    qptrdiff(offsetof(qt_meta_stringdata_QtInput_t, stringdata0) + ofs \
-        - idx * sizeof(QByteArrayData)) \
-    )
-static const qt_meta_stringdata_QtInput_t qt_meta_stringdata_QtInput = {
-    {
-QT_MOC_LITERAL(0, 0, 7), // "QtInput"
-QT_MOC_LITERAL(1, 8, 10), // "run_signal"
-QT_MOC_LITERAL(2, 19, 0), // ""
-QT_MOC_LITERAL(3, 20, 11), // "exit_signal"
-QT_MOC_LITERAL(4, 32, 3), // "run"
-QT_MOC_LITERAL(5, 36, 4), // "save"
-QT_MOC_LITERAL(6, 41, 6), // "cancel"
-QT_MOC_LITERAL(7, 48, 14), // "changeTabIndex"
-QT_MOC_LITERAL(8, 63, 5), // "index"
-QT_MOC_LITERAL(9, 69, 14), // "updateGeometry"
-QT_MOC_LITERAL(10, 84, 18), // "updateInitialState"
-QT_MOC_LITERAL(11, 103, 12) // "updateSource"
-
-    },
-    "QtInput\0run_signal\0\0exit_signal\0run\0"
-    "save\0cancel\0changeTabIndex\0index\0"
-    "updateGeometry\0updateInitialState\0"
-    "updateSource"
-};
-#undef QT_MOC_LITERAL
-
-static const uint qt_meta_data_QtInput[] = {
-
- // content:
-       8,       // revision
-       0,       // classname
-       0,    0, // classinfo
-       9,   14, // methods
-       0,    0, // properties
-       0,    0, // enums/sets
-       0,    0, // constructors
-       0,       // flags
-       2,       // signalCount
-
- // signals: name, argc, parameters, tag, flags
-       1,    0,   59,    2, 0x06 /* Public */,
-       3,    0,   60,    2, 0x06 /* Public */,
-
- // slots: name, argc, parameters, tag, flags
-       4,    0,   61,    2, 0x08 /* Private */,
-       5,    0,   62,    2, 0x08 /* Private */,
-       6,    0,   63,    2, 0x08 /* Private */,
-       7,    1,   64,    2, 0x08 /* Private */,
-       9,    0,   67,    2, 0x08 /* Private */,
-      10,    0,   68,    2, 0x08 /* Private */,
-      11,    0,   69,    2, 0x08 /* Private */,
-
- // signals: parameters
-    QMetaType::Void,
-    QMetaType::Void,
-
- // slots: parameters
-    QMetaType::Void,
-    QMetaType::Void,
-    QMetaType::Void,
-    QMetaType::Void, QMetaType::Int,    8,
-    QMetaType::Void,
-    QMetaType::Void,
-    QMetaType::Void,
-
-       0        // eod
-};
-
-void QtInput::qt_static_metacall(QObject *_o, QMetaObject::Call _c, int _id, void **_a)
-{
-    if (_c == QMetaObject::InvokeMetaMethod) {
-        auto *_t = static_cast<QtInput *>(_o);
-        (void)_t;
-        switch (_id) {
-        case 0: _t->run_signal(); break;
-        case 1: _t->exit_signal(); break;
-        case 2: _t->run(); break;
-        case 3: _t->save(); break;
-        case 4: _t->cancel(); break;
-        case 5: _t->changeTabIndex((*reinterpret_cast< int(*)>(_a[1]))); break;
-        case 6: _t->updateGeometry(); break;
-        case 7: _t->updateInitialState(); break;
-        case 8: _t->updateSource(); break;
-        default: ;
-        }
-    } else if (_c == QMetaObject::IndexOfMethod) {
-        int *result = reinterpret_cast<int *>(_a[0]);
-        {
-            using _t = void (QtInput::*)();
-            if (*reinterpret_cast<_t *>(_a[1]) == static_cast<_t>(&QtInput::run_signal)) {
-                *result = 0;
-                return;
-            }
-        }
-        {
-            using _t = void (QtInput::*)();
-            if (*reinterpret_cast<_t *>(_a[1]) == static_cast<_t>(&QtInput::exit_signal)) {
-                *result = 1;
-                return;
-            }
-        }
-    }
-}
-
-QT_INIT_METAOBJECT const QMetaObject QtInput::staticMetaObject = { {
-    QMetaObject::SuperData::link<QWidget::staticMetaObject>(),
-    qt_meta_stringdata_QtInput.data,
-    qt_meta_data_QtInput,
-    qt_static_metacall,
-    nullptr,
-    nullptr
-} };
-
-
-const QMetaObject *QtInput::metaObject() const
-{
-    return QObject::d_ptr->metaObject ? QObject::d_ptr->dynamicMetaObject() : &staticMetaObject;
-}
-
-void *QtInput::qt_metacast(const char *_clname)
-{
-    if (!_clname) return nullptr;
-    if (!strcmp(_clname, qt_meta_stringdata_QtInput.stringdata0))
-        return static_cast<void*>(this);
-    return QWidget::qt_metacast(_clname);
-}
-
-int QtInput::qt_metacall(QMetaObject::Call _c, int _id, void **_a)
-{
-    _id = QWidget::qt_metacall(_c, _id, _a);
-    if (_id < 0)
-        return _id;
-    if (_c == QMetaObject::InvokeMetaMethod) {
-        if (_id < 9)
-            qt_static_metacall(this, _c, _id, _a);
-        _id -= 9;
-    } else if (_c == QMetaObject::RegisterMethodArgumentMetaType) {
-        if (_id < 9)
-            *reinterpret_cast<int*>(_a[0]) = -1;
-        _id -= 9;
-    }
-    return _id;
-}
-
-// SIGNAL 0
-void QtInput::run_signal()
-{
-    QMetaObject::activate(this, &staticMetaObject, 0, nullptr);
-}
-
-// SIGNAL 1
-void QtInput::exit_signal()
-{
-    QMetaObject::activate(this, &staticMetaObject, 1, nullptr);
-}
-QT_WARNING_POP
-QT_END_MOC_NAMESPACE

moc/input_geometry.cpp

@@ -1,143 +0,0 @@
-/****************************************************************************
-** Meta object code from reading C++ file 'input_geometry.hpp'
-**
-** Created by: The Qt Meta Object Compiler version 67 (Qt 5.15.2)
-**
-** WARNING! All changes made in this file will be lost!
-*****************************************************************************/
-
-#include <memory>
-#include "../src/qt/input_geometry.hpp"
-#include <QtCore/qbytearray.h>
-#include <QtCore/qmetatype.h>
-#if !defined(Q_MOC_OUTPUT_REVISION)
-#error "The header file 'input_geometry.hpp' doesn't include <QObject>."
-#elif Q_MOC_OUTPUT_REVISION != 67
-#error "This file was generated using the moc from 5.15.2. It"
-#error "cannot be used with the include files from this version of Qt."
-#error "(The moc has changed too much.)"
-#endif
-
-QT_BEGIN_MOC_NAMESPACE
-QT_WARNING_PUSH
-QT_WARNING_DISABLE_DEPRECATED
-struct qt_meta_stringdata_QtInputGeometry_t {
-    QByteArrayData data[4];
-    char stringdata0[53];
-};
-#define QT_MOC_LITERAL(idx, ofs, len) \
-    Q_STATIC_BYTE_ARRAY_DATA_HEADER_INITIALIZER_WITH_OFFSET(len, \
-    qptrdiff(offsetof(qt_meta_stringdata_QtInputGeometry_t, stringdata0) + ofs \
-        - idx * sizeof(QByteArrayData)) \
-    )
-static const qt_meta_stringdata_QtInputGeometry_t qt_meta_stringdata_QtInputGeometry = {
-    {
-QT_MOC_LITERAL(0, 0, 15), // "QtInputGeometry"
-QT_MOC_LITERAL(1, 16, 15), // "geometryChanged"
-QT_MOC_LITERAL(2, 32, 0), // ""
-QT_MOC_LITERAL(3, 33, 19) // "emitGeometryChanged"
-
-    },
-    "QtInputGeometry\0geometryChanged\0\0"
-    "emitGeometryChanged"
-};
-#undef QT_MOC_LITERAL
-
-static const uint qt_meta_data_QtInputGeometry[] = {
-
- // content:
-       8,       // revision
-       0,       // classname
-       0,    0, // classinfo
-       2,   14, // methods
-       0,    0, // properties
-       0,    0, // enums/sets
-       0,    0, // constructors
-       0,       // flags
-       1,       // signalCount
-
- // signals: name, argc, parameters, tag, flags
-       1,    0,   24,    2, 0x06 /* Public */,
-
- // slots: name, argc, parameters, tag, flags
-       3,    0,   25,    2, 0x0a /* Public */,
-
- // signals: parameters
-    QMetaType::Void,
-
- // slots: parameters
-    QMetaType::Void,
-
-       0        // eod
-};
-
-void QtInputGeometry::qt_static_metacall(QObject *_o, QMetaObject::Call _c, int _id, void **_a)
-{
-    if (_c == QMetaObject::InvokeMetaMethod) {
-        auto *_t = static_cast<QtInputGeometry *>(_o);
-        (void)_t;
-        switch (_id) {
-        case 0: _t->geometryChanged(); break;
-        case 1: _t->emitGeometryChanged(); break;
-        default: ;
-        }
-    } else if (_c == QMetaObject::IndexOfMethod) {
-        int *result = reinterpret_cast<int *>(_a[0]);
-        {
-            using _t = void (QtInputGeometry::*)();
-            if (*reinterpret_cast<_t *>(_a[1]) == static_cast<_t>(&QtInputGeometry::geometryChanged)) {
-                *result = 0;
-                return;
-            }
-        }
-    }
-    (void)_a;
-}
-
-QT_INIT_METAOBJECT const QMetaObject QtInputGeometry::staticMetaObject = { {
-    QMetaObject::SuperData::link<QWidget::staticMetaObject>(),
-    qt_meta_stringdata_QtInputGeometry.data,
-    qt_meta_data_QtInputGeometry,
-    qt_static_metacall,
-    nullptr,
-    nullptr
-} };
-
-
-const QMetaObject *QtInputGeometry::metaObject() const
-{
-    return QObject::d_ptr->metaObject ? QObject::d_ptr->dynamicMetaObject() : &staticMetaObject;
-}
-
-void *QtInputGeometry::qt_metacast(const char *_clname)
-{
-    if (!_clname) return nullptr;
-    if (!strcmp(_clname, qt_meta_stringdata_QtInputGeometry.stringdata0))
-        return static_cast<void*>(this);
-    return QWidget::qt_metacast(_clname);
-}
-
-int QtInputGeometry::qt_metacall(QMetaObject::Call _c, int _id, void **_a)
-{
-    _id = QWidget::qt_metacall(_c, _id, _a);
-    if (_id < 0)
-        return _id;
-    if (_c == QMetaObject::InvokeMetaMethod) {
-        if (_id < 2)
-            qt_static_metacall(this, _c, _id, _a);
-        _id -= 2;
-    } else if (_c == QMetaObject::RegisterMethodArgumentMetaType) {
-        if (_id < 2)
-            *reinterpret_cast<int*>(_a[0]) = -1;
-        _id -= 2;
-    }
-    return _id;
-}
-
-// SIGNAL 0
-void QtInputGeometry::geometryChanged()
-{
-    QMetaObject::activate(this, &staticMetaObject, 0, nullptr);
-}
-QT_WARNING_POP
-QT_END_MOC_NAMESPACE

moc/input_physics.cpp

@@ -1,119 +0,0 @@
-/****************************************************************************
-** Meta object code from reading C++ file 'input_physics.hpp'
-**
-** Created by: The Qt Meta Object Compiler version 67 (Qt 5.15.2)
-**
-** WARNING! All changes made in this file will be lost!
-*****************************************************************************/
-
-#include <memory>
-#include "../src/qt/input_physics.hpp"
-#include <QtCore/qbytearray.h>
-#include <QtCore/qmetatype.h>
-#if !defined(Q_MOC_OUTPUT_REVISION)
-#error "The header file 'input_physics.hpp' doesn't include <QObject>."
-#elif Q_MOC_OUTPUT_REVISION != 67
-#error "This file was generated using the moc from 5.15.2. It"
-#error "cannot be used with the include files from this version of Qt."
-#error "(The moc has changed too much.)"
-#endif
-
-QT_BEGIN_MOC_NAMESPACE
-QT_WARNING_PUSH
-QT_WARNING_DISABLE_DEPRECATED
-struct qt_meta_stringdata_QtInputPhysics_t {
-    QByteArrayData data[4];
-    char stringdata0[35];
-};
-#define QT_MOC_LITERAL(idx, ofs, len) \
-    Q_STATIC_BYTE_ARRAY_DATA_HEADER_INITIALIZER_WITH_OFFSET(len, \
-    qptrdiff(offsetof(qt_meta_stringdata_QtInputPhysics_t, stringdata0) + ofs \
-        - idx * sizeof(QByteArrayData)) \
-    )
-static const qt_meta_stringdata_QtInputPhysics_t qt_meta_stringdata_QtInputPhysics = {
-    {
-QT_MOC_LITERAL(0, 0, 14), // "QtInputPhysics"
-QT_MOC_LITERAL(1, 15, 12), // "modelChoosed"
-QT_MOC_LITERAL(2, 28, 0), // ""
-QT_MOC_LITERAL(3, 29, 5) // "index"
-
-    },
-    "QtInputPhysics\0modelChoosed\0\0index"
-};
-#undef QT_MOC_LITERAL
-
-static const uint qt_meta_data_QtInputPhysics[] = {
-
- // content:
-       8,       // revision
-       0,       // classname
-       0,    0, // classinfo
-       1,   14, // methods
-       0,    0, // properties
-       0,    0, // enums/sets
-       0,    0, // constructors
-       0,       // flags
-       0,       // signalCount
-
- // slots: name, argc, parameters, tag, flags
-       1,    1,   19,    2, 0x08 /* Private */,
-
- // slots: parameters
-    QMetaType::Void, QMetaType::Int,    3,
-
-       0        // eod
-};
-
-void QtInputPhysics::qt_static_metacall(QObject *_o, QMetaObject::Call _c, int _id, void **_a)
-{
-    if (_c == QMetaObject::InvokeMetaMethod) {
-        auto *_t = static_cast<QtInputPhysics *>(_o);
-        (void)_t;
-        switch (_id) {
-        case 0: _t->modelChoosed((*reinterpret_cast< int(*)>(_a[1]))); break;
-        default: ;
-        }
-    }
-}
-
-QT_INIT_METAOBJECT const QMetaObject QtInputPhysics::staticMetaObject = { {
-    QMetaObject::SuperData::link<QWidget::staticMetaObject>(),
-    qt_meta_stringdata_QtInputPhysics.data,
-    qt_meta_data_QtInputPhysics,
-    qt_static_metacall,
-    nullptr,
-    nullptr
-} };
-
-
-const QMetaObject *QtInputPhysics::metaObject() const
-{
-    return QObject::d_ptr->metaObject ? QObject::d_ptr->dynamicMetaObject() : &staticMetaObject;
-}
-
-void *QtInputPhysics::qt_metacast(const char *_clname)
-{
-    if (!_clname) return nullptr;
-    if (!strcmp(_clname, qt_meta_stringdata_QtInputPhysics.stringdata0))
-        return static_cast<void*>(this);
-    return QWidget::qt_metacast(_clname);
-}
-
-int QtInputPhysics::qt_metacall(QMetaObject::Call _c, int _id, void **_a)
-{
-    _id = QWidget::qt_metacall(_c, _id, _a);
-    if (_id < 0)
-        return _id;
-    if (_c == QMetaObject::InvokeMetaMethod) {
-        if (_id < 1)
-            qt_static_metacall(this, _c, _id, _a);
-        _id -= 1;
-    } else if (_c == QMetaObject::RegisterMethodArgumentMetaType) {
-        if (_id < 1)
-            *reinterpret_cast<int*>(_a[0]) = -1;
-        _id -= 1;
-    }
-    return _id;
-}
-QT_WARNING_POP
-QT_END_MOC_NAMESPACE

moc/input_view.cpp

@@ -1,133 +0,0 @@
-/****************************************************************************
-** Meta object code from reading C++ file 'input_view.hpp'
-**
-** Created by: The Qt Meta Object Compiler version 67 (Qt 5.15.2)
-**
-** WARNING! All changes made in this file will be lost!
-*****************************************************************************/
-
-#include <memory>
-#include "../src/qt/input_view.hpp"
-#include <QtCore/qbytearray.h>
-#include <QtCore/qmetatype.h>
-#if !defined(Q_MOC_OUTPUT_REVISION)
-#error "The header file 'input_view.hpp' doesn't include <QObject>."
-#elif Q_MOC_OUTPUT_REVISION != 67
-#error "This file was generated using the moc from 5.15.2. It"
-#error "cannot be used with the include files from this version of Qt."
-#error "(The moc has changed too much.)"
-#endif
-
-QT_BEGIN_MOC_NAMESPACE
-QT_WARNING_PUSH
-QT_WARNING_DISABLE_DEPRECATED
-struct qt_meta_stringdata_QtInputView_t {
-    QByteArrayData data[7];
-    char stringdata0[70];
-};
-#define QT_MOC_LITERAL(idx, ofs, len) \
-    Q_STATIC_BYTE_ARRAY_DATA_HEADER_INITIALIZER_WITH_OFFSET(len, \
-    qptrdiff(offsetof(qt_meta_stringdata_QtInputView_t, stringdata0) + ofs \
-        - idx * sizeof(QByteArrayData)) \
-    )
-static const qt_meta_stringdata_QtInputView_t qt_meta_stringdata_QtInputView = {
-    {
-QT_MOC_LITERAL(0, 0, 11), // "QtInputView"
-QT_MOC_LITERAL(1, 12, 14), // "updateGeometry"
-QT_MOC_LITERAL(2, 27, 0), // ""
-QT_MOC_LITERAL(3, 28, 18), // "updateInitialState"
-QT_MOC_LITERAL(4, 47, 12), // "updateSource"
-QT_MOC_LITERAL(5, 60, 7), // "setTime"
-QT_MOC_LITERAL(6, 68, 1) // "v"
-
-    },
-    "QtInputView\0updateGeometry\0\0"
-    "updateInitialState\0updateSource\0setTime\0"
-    "v"
-};
-#undef QT_MOC_LITERAL
-
-static const uint qt_meta_data_QtInputView[] = {
-
- // content:
-       8,       // revision
-       0,       // classname
-       0,    0, // classinfo
-       4,   14, // methods
-       0,    0, // properties
-       0,    0, // enums/sets
-       0,    0, // constructors
-       0,       // flags
-       0,       // signalCount
-
- // slots: name, argc, parameters, tag, flags
-       1,    0,   34,    2, 0x0a /* Public */,
-       3,    0,   35,    2, 0x0a /* Public */,
-       4,    0,   36,    2, 0x0a /* Public */,
-       5,    1,   37,    2, 0x0a /* Public */,
-
- // slots: parameters
-    QMetaType::Void,
-    QMetaType::Void,
-    QMetaType::Void,
-    QMetaType::Void, QMetaType::Int,    6,
-
-       0        // eod
-};
-
-void QtInputView::qt_static_metacall(QObject *_o, QMetaObject::Call _c, int _id, void **_a)
-{
-    if (_c == QMetaObject::InvokeMetaMethod) {
-        auto *_t = static_cast<QtInputView *>(_o);
-        (void)_t;
-        switch (_id) {
-        case 0: _t->updateGeometry(); break;
-        case 1: _t->updateInitialState(); break;
-        case 2: _t->updateSource(); break;
-        case 3: _t->setTime((*reinterpret_cast< int(*)>(_a[1]))); break;
-        default: ;
-        }
-    }
-}
-
-QT_INIT_METAOBJECT const QMetaObject QtInputView::staticMetaObject = { {
-    QMetaObject::SuperData::link<QtView::staticMetaObject>(),
-    qt_meta_stringdata_QtInputView.data,
-    qt_meta_data_QtInputView,
-    qt_static_metacall,
-    nullptr,
-    nullptr
-} };
-
-
-const QMetaObject *QtInputView::metaObject() const
-{
-    return QObject::d_ptr->metaObject ? QObject::d_ptr->dynamicMetaObject() : &staticMetaObject;
-}
-
-void *QtInputView::qt_metacast(const char *_clname)
-{
-    if (!_clname) return nullptr;
-    if (!strcmp(_clname, qt_meta_stringdata_QtInputView.stringdata0))
-        return static_cast<void*>(this);
-    return QtView::qt_metacast(_clname);
-}
-
-int QtInputView::qt_metacall(QMetaObject::Call _c, int _id, void **_a)
-{
-    _id = QtView::qt_metacall(_c, _id, _a);
-    if (_id < 0)
-        return _id;
-    if (_c == QMetaObject::InvokeMetaMethod) {
-        if (_id < 4)
-            qt_static_metacall(this, _c, _id, _a);
-        _id -= 4;
-    } else if (_c == QMetaObject::RegisterMethodArgumentMetaType) {
-        if (_id < 4)
-            *reinterpret_cast<int*>(_a[0]) = -1;
-        _id -= 4;
-    }
-    return _id;
-}
-QT_WARNING_POP
-QT_END_MOC_NAMESPACE

moc/mainwindow.cpp

@@ -1,136 +0,0 @@
-/****************************************************************************
-** Meta object code from reading C++ file 'mainwindow.hpp'
-**
-** Created by: The Qt Meta Object Compiler version 67 (Qt 5.15.2)
-**
-** WARNING! All changes made in this file will be lost!
-*****************************************************************************/
-
-#include <memory>
-#include "../src/qt/mainwindow.hpp"
-#include <QtCore/qbytearray.h>
-#include <QtCore/qmetatype.h>
-#if !defined(Q_MOC_OUTPUT_REVISION)
-#error "The header file 'mainwindow.hpp' doesn't include <QObject>."
-#elif Q_MOC_OUTPUT_REVISION != 67
-#error "This file was generated using the moc from 5.15.2. It"
-#error "cannot be used with the include files from this version of Qt."
-#error "(The moc has changed too much.)"
-#endif
-
-QT_BEGIN_MOC_NAMESPACE
-QT_WARNING_PUSH
-QT_WARNING_DISABLE_DEPRECATED
-struct qt_meta_stringdata_QtMainWindow_t {
-    QByteArrayData data[7];
-    char stringdata0[61];
-};
-#define QT_MOC_LITERAL(idx, ofs, len) \
-    Q_STATIC_BYTE_ARRAY_DATA_HEADER_INITIALIZER_WITH_OFFSET(len, \
-    qptrdiff(offsetof(qt_meta_stringdata_QtMainWindow_t, stringdata0) + ofs \
-        - idx * sizeof(QByteArrayData)) \
-    )
-static const qt_meta_stringdata_QtMainWindow_t qt_meta_stringdata_QtMainWindow = {
-    {
-QT_MOC_LITERAL(0, 0, 12), // "QtMainWindow"
-QT_MOC_LITERAL(1, 13, 9), // "new_input"
-QT_MOC_LITERAL(2, 23, 0), // ""
-QT_MOC_LITERAL(3, 24, 10), // "load_input"
-QT_MOC_LITERAL(4, 35, 4), // "exit"
-QT_MOC_LITERAL(5, 40, 9), // "run_input"
-QT_MOC_LITERAL(6, 50, 10) // "exit_input"
-
-    },
-    "QtMainWindow\0new_input\0\0load_input\0"
-    "exit\0run_input\0exit_input"
-};
-#undef QT_MOC_LITERAL
-
-static const uint qt_meta_data_QtMainWindow[] = {
-
- // content:
-       8,       // revision
-       0,       // classname
-       0,    0, // classinfo
-       5,   14, // methods
-       0,    0, // properties
-       0,    0, // enums/sets
-       0,    0, // constructors
-       0,       // flags
-       0,       // signalCount
-
- // slots: name, argc, parameters, tag, flags
-       1,    0,   39,    2, 0x08 /* Private */,
-       3,    0,   40,    2, 0x08 /* Private */,
-       4,    0,   41,    2, 0x08 /* Private */,
-       5,    0,   42,    2, 0x08 /* Private */,
-       6,    0,   43,    2, 0x08 /* Private */,
-
- // slots: parameters
-    QMetaType::Void,
-    QMetaType::Void,
-    QMetaType::Void,
-    QMetaType::Void,
-    QMetaType::Void,
-
-       0        // eod
-};
-
-void QtMainWindow::qt_static_metacall(QObject *_o, QMetaObject::Call _c, int _id, void **_a)
-{
-    if (_c == QMetaObject::InvokeMetaMethod) {
-        auto *_t = static_cast<QtMainWindow *>(_o);
-        (void)_t;
-        switch (_id) {
-        case 0: _t->new_input(); break;
-        case 1: _t->load_input(); break;
-        case 2: _t->exit(); break;
-        case 3: _t->run_input(); break;
-        case 4: _t->exit_input(); break;
-        default: ;
-        }
-    }
-    (void)_a;
-}
-
-QT_INIT_METAOBJECT const QMetaObject QtMainWindow::staticMetaObject = { {
-    QMetaObject::SuperData::link<QMainWindow::staticMetaObject>(),
-    qt_meta_stringdata_QtMainWindow.data,
-    qt_meta_data_QtMainWindow,
-    qt_static_metacall,
-    nullptr,
-    nullptr
-} };
-
-
-const QMetaObject *QtMainWindow::metaObject() const
-{
-    return QObject::d_ptr->metaObject ? QObject::d_ptr->dynamicMetaObject() : &staticMetaObject;
-}
-
-void *QtMainWindow::qt_metacast(const char *_clname)
-{
-    if (!_clname) return nullptr;
-    if (!strcmp(_clname, qt_meta_stringdata_QtMainWindow.stringdata0))
-        return static_cast<void*>(this);
-    return QMainWindow::qt_metacast(_clname);
-}
-
-int QtMainWindow::qt_metacall(QMetaObject::Call _c, int _id, void **_a)
-{
-    _id = QMainWindow::qt_metacall(_c, _id, _a);
-    if (_id < 0)
-        return _id;
-    if (_c == QMetaObject::InvokeMetaMethod) {
-        if (_id < 5)
-            qt_static_metacall(this, _c, _id, _a);
-        _id -= 5;
-    } else if (_c == QMetaObject::RegisterMethodArgumentMetaType) {
-        if (_id < 5)
-            *reinterpret_cast<int*>(_a[0]) = -1;
-        _id -= 5;
-    }
-    return _id;
-}
-QT_WARNING_POP
-QT_END_MOC_NAMESPACE

moc/view_solution.cpp

@@ -1,118 +0,0 @@
-/****************************************************************************
-** Meta object code from reading C++ file 'view_solution.hpp'
-**
-** Created by: The Qt Meta Object Compiler version 67 (Qt 5.15.2)
-**
-** WARNING! All changes made in this file will be lost!
-*****************************************************************************/
-
-#include <memory>
-#include "../src/qt/view_solution.hpp"
-#include <QtCore/qbytearray.h>
-#include <QtCore/qmetatype.h>
-#if !defined(Q_MOC_OUTPUT_REVISION)
-#error "The header file 'view_solution.hpp' doesn't include <QObject>."
-#elif Q_MOC_OUTPUT_REVISION != 67
-#error "This file was generated using the moc from 5.15.2. It"
-#error "cannot be used with the include files from this version of Qt."
-#error "(The moc has changed too much.)"
-#endif
-
-QT_BEGIN_MOC_NAMESPACE
-QT_WARNING_PUSH
-QT_WARNING_DISABLE_DEPRECATED
-struct qt_meta_stringdata_QtViewSolution_t {
-    QByteArrayData data[3];
-    char stringdata0[28];
-};
-#define QT_MOC_LITERAL(idx, ofs, len) \
-    Q_STATIC_BYTE_ARRAY_DATA_HEADER_INITIALIZER_WITH_OFFSET(len, \
-    qptrdiff(offsetof(qt_meta_stringdata_QtViewSolution_t, stringdata0) + ofs \
-        - idx * sizeof(QByteArrayData)) \
-    )
-static const qt_meta_stringdata_QtViewSolution_t qt_meta_stringdata_QtViewSolution = {
-    {
-QT_MOC_LITERAL(0, 0, 14), // "QtViewSolution"
-QT_MOC_LITERAL(1, 15, 11), // "time_change"
-QT_MOC_LITERAL(2, 27, 0) // ""
-
-    },
-    "QtViewSolution\0time_change\0"
-};
-#undef QT_MOC_LITERAL
-
-static const uint qt_meta_data_QtViewSolution[] = {
-
- // content:
-       8,       // revision
-       0,       // classname
-       0,    0, // classinfo
-       1,   14, // methods
-       0,    0, // properties
-       0,    0, // enums/sets
-       0,    0, // constructors
-       0,       // flags
-       0,       // signalCount
-
- // slots: name, argc, parameters, tag, flags
-       1,    1,   19,    2, 0x08 /* Private */,
-
- // slots: parameters
-    QMetaType::Void, QMetaType::Int,    2,
-
-       0        // eod
-};
-
-void QtViewSolution::qt_static_metacall(QObject *_o, QMetaObject::Call _c, int _id, void **_a)
-{
-    if (_c == QMetaObject::InvokeMetaMethod) {
-        auto *_t = static_cast<QtViewSolution *>(_o);
-        (void)_t;
-        switch (_id) {
-        case 0: _t->time_change((*reinterpret_cast< int(*)>(_a[1]))); break;
-        default: ;
-        }
-    }
-}
-
-QT_INIT_METAOBJECT const QMetaObject QtViewSolution::staticMetaObject = { {
-    QMetaObject::SuperData::link<QWidget::staticMetaObject>(),
-    qt_meta_stringdata_QtViewSolution.data,
-    qt_meta_data_QtViewSolution,
-    qt_static_metacall,
-    nullptr,
-    nullptr
-} };
-
-
-const QMetaObject *QtViewSolution::metaObject() const
-{
-    return QObject::d_ptr->metaObject ? QObject::d_ptr->dynamicMetaObject() : &staticMetaObject;
-}
-
-void *QtViewSolution::qt_metacast(const char *_clname)
-{
-    if (!_clname) return nullptr;
-    if (!strcmp(_clname, qt_meta_stringdata_QtViewSolution.stringdata0))
-        return static_cast<void*>(this);
-    return QWidget::qt_metacast(_clname);
-}
-
-int QtViewSolution::qt_metacall(QMetaObject::Call _c, int _id, void **_a)
-{
-    _id = QWidget::qt_metacall(_c, _id, _a);
-    if (_id < 0)
-        return _id;
-    if (_c == QMetaObject::InvokeMetaMethod) {
-        if (_id < 1)
-            qt_static_metacall(this, _c, _id, _a);
-        _id -= 1;
-    } else if (_c == QMetaObject::RegisterMethodArgumentMetaType) {
-        if (_id < 1)
-            *reinterpret_cast<int*>(_a[0]) = -1;
-        _id -= 1;
-    }
-    return _id;
-}
-QT_WARNING_POP
-QT_END_MOC_NAMESPACE

src/geometry.cpp

@@ -15,6 +15,22 @@ Geometry::Geometry(){
   Z=nullptr;
 }
 
+
+Geometry::~Geometry(){
+  if(hsoil!=nullptr){
+    delete[] hsoil;
+    delete[] dhsoil;
+    delete[] hbot;
+    delete[] dhbot;
+    delete[] dZ;
+    delete[] nZ;
+    for(size_t i=0;i<nX;++i){
+      delete[] Z[i];
+    }
+    delete[] Z;
+  }
+}
+
 void
 Geometry::initZ(double dZ_avg,bool init){
   if(init){
@@ -30,76 +46,49 @@ Geometry::initZ(double dZ_avg,bool init){
     nZ[k]=n+1;
     if(init) Z[k]=new double[n+1];
     for(size_t j=0;j<n+1;++j){
-      Z[k][j]=factor*(hbot[k]+j*dz);
+      Z[k][j]=hbot[k]+j*dz;
     }
   }
 }
 
 void
-Geometry::initDefault(){
-  lX=10;
-  nX=100;
-  nX_max=nmax_Qt;
-  dX=lX/(nX-1);
-  hsoil=new double[nmax_Qt];
-  dhsoil=new double[nmax_Qt];
-  hbot=new double[nmax_Qt];
-  dhbot=new double[nmax_Qt];
-  nZ=new size_t[nmax_Qt];
-  dZ=new double[nmax_Qt];
-  Z=new double*[nmax_Qt];
-  for(size_t i=0;i<nmax_Qt;++i){
-    Z[i]=new double[nmax_Qt];
-  }
+Geometry::save(fstream& file) {
+  file.write((char*)&lX,sizeof(double));
+  file.write((char*)&nX,sizeof(size_t));
+  file.write((char*)&dX,sizeof(double));
+  file.write((char*)hsoil,nX*sizeof(double));
+  file.write((char*)dhsoil,nX*sizeof(double));
+  file.write((char*)hbot,nX*sizeof(double));
+  file.write((char*)dhbot,nX*sizeof(double));
+  file.write((char*)nZ,nX*sizeof(size_t));
+  file.write((char*)dZ,nX*sizeof(double));
   for(size_t i=0;i<nX;++i){
-    hsoil[i]=0.8;
-    hbot[i]=0.2;
-    dhsoil[i]=0;
-    dhbot[i]=0;
+    file.write((char*)Z[i],nZ[i]*sizeof(double));
   }
-  depth=5;
-  nZ_max=100;
-  factor=5/0.6;
-  initZ(depth/nZ_max,false);
 }
 
 void
-Geometry::update(double _lX,size_t _nX,double _depth,size_t _nZ_max,Spline& Ssoil,Spline& Sbot){
-  lX=_lX;
-  depth=_depth;
-  nX=_nX;
-  nZ_max=_nZ_max;
-  dX=lX/(nX-1);
-  double hs_max=-inf,hb_min=inf;
-  double v;
-  for(size_t k=0;k<nX;++k){
-    double x=k*dX/lX;
-    v=hsoil[k]=Ssoil(x);
-    hs_max=max(v,hs_max);
-    dhsoil[k]=Ssoil.derivate(x);
-    v=hbot[k]=Sbot(x);
-    hb_min=min(v,hb_min);
-    dhbot[k]=Sbot.derivate(x);
+Geometry::load(fstream& file,bool init) {
+  file.read((char*)&lX,sizeof(double));
+  file.read((char*)&nX,sizeof(size_t));
+  file.read((char*)&dX,sizeof(double));
+  if(init){
+    hsoil=new double[nX];
+    dhsoil=new double[nX];
+    hbot=new double[nX];
+    dhbot=new double[nX];
+    nZ=new size_t[nX];
+    dZ=new double[nX];
+    Z=new double*[nX];
   }
-  double dZ_avg=(hs_max-hb_min)/_nZ_max;
-  factor=depth/(hs_max-hb_min);
-  initZ(dZ_avg,false);
-}
-
-Geometry::~Geometry(){
-  if(hsoil!=nullptr){
-    delete[] hsoil;
-    delete[] dhsoil;
-    delete[] hbot;
-    delete[] dhbot;
-    delete[] dZ;
-    delete[] nZ;
-    for(size_t i=0;i<nX_max;++i){
-      delete[] Z[i];
-    }
-    delete[] Z;
-
+  file.read((char*)hsoil,nX*sizeof(double));
+  file.read((char*)dhsoil,nX*sizeof(double));
+  file.read((char*)hbot,nX*sizeof(double));
+  file.read((char*)dhbot,nX*sizeof(double));
+  file.read((char*)nZ,nX*sizeof(size_t));
+  file.read((char*)dZ,nX*sizeof(double));
+  for(size_t i=0;i<nX;++i){
+    if(init) Z[i]=new double[nZ[i]];
+    file.read((char*)Z[i],nZ[i]*sizeof(double));
   }
-
-
 }

src/geometry.hpp

@@ -2,6 +2,7 @@
 #define GEOMETRY_HPP
 
 #include <iostream>
+#include <fstream>
 #include <limits>
 
 #include "math/spline.hpp"
@@ -15,33 +16,15 @@ extern double inf;
 //! The Geometry class contains all geometric parameters of the domain.
 class Geometry{
 public:
-
-  static constexpr size_t const nmax_Qt=400;
-  //!
-  Geometry();
-
-  //! Geometry constructor
-  Geometry(double _lX,size_t _nX,size_t _nZ,Func _hsoil,Func _dhsoil,Func _hbot,Func _dhbot);
-
-  ~Geometry();
-
-  void initDefault();
-
   //! Horizontal length of the domain
   double lX;
 
   //! Number of horizontal steps
   size_t nX;
-  size_t nX_max;
 
   //! Horizontal step
   double dX;
 
-  size_t nZ_max;
-
-  double depth;
-
-
   //! Level of the soil depending on X, e.g, hsoil[k]=level of the soil at X=k*dX.
   //! Vector of size nX.
   double* hsoil;
@@ -59,17 +42,23 @@ public:
 
   //! Number of vertical step at a given X, vector of size nX.
   size_t* nZ;
-
   //! Vertical step at a given X, vector of size nX.
   double* dZ;
 
   //! Vertical considered positions at a given X, vector of vectors of size nX. For each k, Z[k] is a vector of size nZ[k]
   double** Z;
 
-  double factor;
+  //!
+  Geometry();
+
+  ~Geometry();
 
   void initZ(double dZ_avg,bool init=true);
 
-  void update(double _lX,size_t _nX,double _depth,size_t _nZ_max,Spline& hbot,Spline& hsoil);
+  void save(fstream& file);
+
+  void load(fstream& file,bool init);
 };
+
+
 #endif

src/horizontal.hpp

@@ -0,0 +1,13 @@
+#ifndef HORIZONTAL_HPP
+#define HORIZONTAL_HPP
+
+
+class Horizontal{
+private:
+public:
+  Horizontal(Geometry& G);
+  bool
+
+};
+
+#endif

src/initial_state.cpp

@@ -10,40 +10,50 @@ Tank::Tank(){
   delta_top=0.05;
 }
 
+void
+Tank::save(fstream& file){
+  file.write((char*)&saturation,sizeof(double));
+  file.write((char*)&left,sizeof(double));
+  file.write((char*)&right,sizeof(double));
+  file.write((char*)&bottom,sizeof(double));
+  file.write((char*)&top,sizeof(double));
+  file.write((char*)&delta_left,sizeof(double));
+  file.write((char*)&delta_right,sizeof(double));
+  file.write((char*)&delta_bottom,sizeof(double));
+  file.write((char*)&delta_top,sizeof(double));
+}
 
-
+void
+Tank::load(fstream& file){
+  file.read((char*)&saturation,sizeof(double));
+  file.read((char*)&left,sizeof(double));
+  file.read((char*)&right,sizeof(double));
+  file.read((char*)&bottom,sizeof(double));
+  file.read((char*)&top,sizeof(double));
+  file.read((char*)&delta_left,sizeof(double));
+  file.read((char*)&delta_right,sizeof(double));
+  file.read((char*)&delta_bottom,sizeof(double));
+  file.read((char*)&delta_top,sizeof(double));
+}
 
 InitialState::InitialState(){
+  geometry=nullptr;
   hsat=nullptr;
   Pinit=nullptr;
 }
 
-void
-InitialState::update(Geometry& _geometry,Spline& _Ssat){
-  Ssat=&_Ssat;
-  geometry=&_geometry;
-  if(hsat!=nullptr){
-    delete[] hsat;
-    for(size_t i=0;i<geometry->nX;++i){
-      delete[] Pinit[i];
-    }
-    delete[] Pinit;
-  }
-  hsat=new double[geometry->nX];
-  Pinit=new double*[geometry->nX];
-
-  for(size_t i=0;i<geometry->nX;++i){
-    size_t nZ=geometry->nZ[i];
-    Pinit[i]=new double[nZ];
-  }
-  updatePressure();
+InitialState::InitialState(Geometry* _geometry){
+  geometry=_geometry;
+  hsat=nullptr;
+  Pinit=nullptr;
 }
+
 Tank*
 InitialState::addTank(){
   Tank* tank=new Tank;
   tanks.push_back(tank);
+  //updatePressure();
   return tank;
-  updatePressure();
 }
 
 void
@@ -55,10 +65,9 @@ InitialState::removeTank(Tank* tank){
       return;
     }
   }
-  updatePressure();
+  //updatePressure();
 }
 
-
 InitialState::~InitialState(){
   if(hsat!=nullptr){
     delete[] hsat;
@@ -66,20 +75,19 @@ InitialState::~InitialState(){
       delete[] Pinit[i];
     }
     delete[] Pinit;
-  }
-  for(auto it=tanks.begin();it!=tanks.end();++it){
-    delete(*it);
+    for(auto it=tanks.begin();it!=tanks.end();++it){
+      delete(*it);
+    }
   }
 }
 
 void
 InitialState::updatePressure(){
-
   for(size_t i=0;i<geometry->nX;++i){
     double x=double(i)/(geometry->nX-1);
-    hsat[i]=(*Ssat)(x)*geometry->factor;
     size_t nZ=geometry->nZ[i];
-    double temp=hsat[i]+Physics::model_data[0]/(Physics::g*Physics::rho); //TODO : replace Physics::model_data[0] by Psat in future
+    double temp=hsat[i]+Psat/(Physics::g*Physics::rho); //TODO : replace Physics::model_data[0] by Psat in future
+    //cout<<i<<" -> "<<nZ<<endl;
     for(size_t j=0;j<nZ;++j){
       double z=geometry->Z[i][j];
       Pinit[i][j]=Physics::rho*Physics::g*(temp-z);
@@ -96,33 +104,63 @@ InitialState::updatePressure(){
         double db=tank->delta_bottom;
 
         double Ps=Physics::s_inv(S);
-      cout<<S<<" "<<Ps<<endl;
-      double Px,Pz;
-      double zn=z/geometry->factor;
-      if(x<=l){
-        Px=(Ps-Physics::Psec)*(x-l)/dl+Ps;
-      }
-      else if(x>=r){
-        Px=(Physics::Psec-Ps)*(x-r)/dr+Ps;
-      }
-      else{
-        Px=Ps;
-      }
-      if(zn<=b){
-        Pz=(Ps-Physics::Psec)*(zn-b)/db+Ps;
-      }
-      else if(zn>=t){
-        Pz=(Physics::Psec-Ps)*(zn-t)/dt+Ps;
+        double Px,Pz;
+        if(x<=l){
+          Px=(Ps-Physics::Psec)*(x-l)/dl+Ps;
+        }
+        else if(x>=r){
+          Px=(Physics::Psec-Ps)*(x-r)/dr+Ps;
+        }
+        else{
+          Px=Ps;
+        }
+        if(z<=b){
+          Pz=(Ps-Physics::Psec)*(z-b)/db+Ps;
+        }
+        else if(z>=t){
+          Pz=(Physics::Psec-Ps)*(z-t)/dt+Ps;
+        }
+        else{
+          Pz=Ps;
+        }
+        Pinit[i][j]=max(max(min(Px,Pz),Pinit[i][j]),Physics::Psec);
       }
-      else{
-        Pz=Ps;
-      }
-      Pinit[i][j]=max(max(min(Px,Pz),Pinit[i][j]),Physics::Psec);
+      //TODO : voir considrer max avec psec
     }
+  }
+}
 
+void
+InitialState::save(fstream& file){
+  file.write((char*)hsat,geometry->nX*sizeof(double));
+  file.write((char*)hov,geometry->nX*sizeof(double));
+  for(size_t i=0;i<geometry->nX;++i){
+    file.write((char*)Pinit[i],geometry->nZ[i]*sizeof(double));
+  }
+  size_t nt=tanks.size();
+  file.write((char*)&nt,sizeof(double));
+  for(auto it=tanks.begin();it!=tanks.end();++it){
+    (*it)->save(file);
+  }
+}
 
-
-      //TODO : voir considrer max avec psec
-    }
+void
+InitialState::load(fstream& file,bool init){
+  if(init){
+    hsat=new double[geometry->nX];
+    hov=new double[geometry->nX];
+    Pinit=new double*[geometry->nX];
+  }
+  file.read((char*)hsat,geometry->nX*sizeof(double));
+  file.read((char*)hov,geometry->nX*sizeof(double));
+  for(size_t i=0;i<geometry->nX;++i){
+    if(init) Pinit[i]=new double[geometry->nZ[i]];
+    file.read((char*)Pinit[i],geometry->nZ[i]*sizeof(double));
+  }
+  size_t nt;
+  file.read((char*)&nt,sizeof(size_t));
+  for(size_t i=0;i<nt;++i){
+    Tank* tank=addTank();
+    tank->load(file);
   }
 }

src/initial_state.hpp

@@ -2,6 +2,7 @@
 #define INITIAL_STATE_HPP
 
 #include <list>
+#include <fstream>
 #include "physics.hpp"
 #include "geometry.hpp"
 #include "math/spline.hpp"
@@ -14,27 +15,31 @@ public:
   double left,right,bottom,top;
   double delta_left,delta_right,delta_bottom,delta_top;
   Tank();
+  void save(fstream& file);
+  void load(fstream& file);
 };
 
-
-
 class InitialState{
 private:
-
-  Spline* Ssat;
   Geometry* geometry;
 public:
+  //! Vector of size nX
   double* hsat;
+
+
+  double* hov;
+
+  //! Vector of vector of size nX. For each k, Pinit[k] is a vector of size geometry->nZ[k]
   double** Pinit;
   list<Tank*> tanks;
   InitialState();
+  InitialState(Geometry* geometry);
   ~InitialState();
-  void update(Geometry& geometry,Spline& Ssat);
   Tank* addTank();
   void removeTank(Tank* tank);
-
   void updatePressure();
+  void save(fstream& file);
+  void load(fstream& file,bool init);
 };
 
-
 #endif

src/input_data.cpp

@@ -0,0 +1,21 @@
+#include "input_data.hpp"
+
+void
+InputData::load(fstream& file,bool init){
+  Physics::load(file);
+  Time::load(file);
+  geometry.load(file,init);
+  initial_state->load(file,init);
+  source.load(file);
+  file.read((char*)&factor,sizeof(double));
+}
+
+void
+InputData::save(fstream& file){
+  Physics::save(file);
+  Time::save(file);
+  geometry.save(file);
+  initial_state->save(file);
+  source.save(file);
+  file.write((char*)&factor,sizeof(double));
+}

src/input_data.hpp

@@ -0,0 +1,37 @@
+#ifndef INPUT_DATA_HPP
+#define INPUT_DATA_HPP
+
+#include "physics.hpp"
+#include "time.hpp"
+#include "geometry.hpp"
+#include "initial_state.hpp"
+#include "source.hpp"
+
+
+class InputData{
+public:
+  Geometry geometry;
+  InitialState* initial_state;
+  Source source;
+  double factor;
+  InputData();
+  ~InputData();
+  void load(fstream& file,bool init=false);
+  void save(fstream& file);
+};
+
+
+
+inline
+InputData::InputData(){
+  initial_state=new InitialState(&geometry);
+}
+
+inline InputData::~InputData(){
+  if(initial_state!=nullptr){
+    delete initial_state;
+  }
+}
+
+
+#endif

src/kernel.cpp

@@ -1,2 +1,491 @@
 #include "kernel.hpp"
 
+Kernel::Kernel(string filename){
+  fstream file;
+  file.open(filename.c_str(),fstream::in|fstream::binary);
+  InputData::load(file,true);
+
+  initSolutions();
+
+  Solution& S0=*solution[0];
+
+  //Init pressure
+  for(size_t x=0;x<geometry.nX;++x){
+    for(size_t z=0;z<geometry.nZ[x];++z){
+      S0.P[x][z]=initial_state->Pinit[x][z];
+    }
+  }
+  //Init hydraulic head
+  for(size_t x=0;x<geometry.nX;++x){
+    S0.hydr[x]=initial_state->hsat[x]+Psat/(Physics::rho*Physics::g);
+  }
+
+  //Init overland level
+  for(size_t x=0;x<geometry.nX;++x){
+    S0.hov[x]=Physics::invPsol(geometry.hsoil[x],S0.P[x][geometry.nZ[x]-1]); //TODO[Chrisophe] à initiliser
+  }
+
+  //Init l
+  for(size_t x=0;x<geometry.nX;++x){
+    double t=initial_state->hsat[x];
+    S0.l[x]=t;
+    S0.hsat[x]=t;
+  }
+
+  //Init Pl
+  for(size_t x=0;x<geometry.nX;++x){
+    S0.Pl[x]=Psat;
+  }
+
+
+  /*indice_x_Richards=new bool[geometry.nX];
+
+  div_w=new double*[geometry.nX];
+  P_temp=new double*[geometry.nX];
+  for(size_t i=0;i<geometry.nX;++i){
+    div_w[i]=new double[geometry.nZ[i]];
+    P_temp[i]=new double[geometry.nZ[i]];
+  }
+  t_sub_A=new double*[geometry.nX];
+  t_sub_B=new double*[geometry.nX];
+  t_sub_C=new double*[geometry.nX];
+  t_diag_A=new double*[geometry.nX];
+  t_diag_B=new double*[geometry.nX];
+  t_diag_C=new double*[geometry.nX];
+  t_sup_A=new double*[geometry.nX];
+  t_sup_B=new double*[geometry.nX];
+  t_sup_C=new double*[geometry.nX];
+  t_F=new double*[geometry.nX];
+  t_G=new double*[geometry.nX];
+  t_S=new double*[geometry.nX];
+  t_H=new double*[geometry.nX];
+  t_R=new double*[geometry.nX];
+  t_I=new double*[geometry.nX];
+  t_BB=new double*[geometry.nX];
+
+  for(size_t i=0;i<geometry.nX;++i){
+    size_t nZ=geometry.nZ[i];
+    t_sub_A[i]=new double[nZ-2];
+    t_sub_B[i]=new double[nZ-2];
+    t_sub_C[i]=new double[nZ-2];
+    t_diag_A[i]=new double[nZ-1];
+    t_diag_B[i]=new double[nZ-1];
+    t_diag_C[i]=new double[nZ-1];
+    t_sup_A[i]=new double[nZ-2];
+    t_sup_B[i]=new double[nZ-2];
+    t_sup_C[i]=new double[nZ-2];
+    t_F[i]=new double[nZ-1];
+    t_G[i]=new double[nZ-1];
+    t_S[i]=new double[nZ-1];
+    t_H[i]=new double[nZ-1];
+    t_R[i]=new double[nZ-1];
+    t_I[i]=new double[nZ-1];
+    t_BB[i]=new double[nZ-1];
+  }*/
+
+  delete initial_state;
+  initial_state=nullptr;
+
+  file.close();
+  step=0;
+  m=1;
+  scheme_error=0;
+}
+
+void
+Kernel::initSolutions(){
+  Solution* S0=new Solution();
+  S0->init(geometry);
+  solution[0]=S0;
+  for(size_t i=1;i<=Time::nT;++i) solution[i]=nullptr;
+  for(size_t i=0;i<Time::nT+2;++i){
+    Solution* S=new Solution();
+    S->init(geometry);
+    pool_solutions.push(S);
+  }
+}
+
+bool
+Kernel::next(){
+  //cout<<"[Kernel::next]"<<endl;
+  if(m>1) m=m/2;
+  spaceSolution();
+  if(scheme_error>0) return false;
+  ++step;
+  //cout<<"[next] done"<<endl;
+  return true;
+}
+
+void
+Kernel::spaceSolution(){
+
+  //cout<<"[Kernel::spaceSolution]"<<endl;
+  size_t n=step+1; //We want solution at time n
+  size_t k=0; //we work on time interval [n-1+k/m,n-1+(k+1)/m]
+  Solution *S=solution[step];
+  while(m<=max_time_subdivisions){
+    //cout<<" k,m = "<<k<<','<<m<<endl;
+    S=Piccard(S);
+    if(S==nullptr){
+      m*=2;
+      k=0;
+      S=solution[step];
+    }
+    else{
+      ++k;
+      if(k==m) break;
+    }
+  }
+  if(m>max_time_subdivisions){
+    cerr<<"Max time subdivision reached"<<endl;
+    scheme_error=1;
+  }
+  //cout<<"[spaceSolution] done"<<endl;
+}
+
+/*Solution*
+Kernel::Piccard(Solution* Sin){
+  //Return a valid solution or nullptr if not converge
+  //cout<<"[Kernel::Piccard]"<<endl;
+  memcpy(Sin.l,Sin.hsat,sizeof(double)*geometry.nX);
+  //Compute Pl with l=hsat of input solution and P=pressure of input_solution
+
+  Solution* Spl=pool_solutions.pop();
+  compute_Pl(*Spl,Sin.hsat,Sin.P);
+  double dt=Time::dT/m; // util ????
+
+  size_t count=0;
+
+  //Specicy to compute Richards on all columns
+  for(size_t i=0;i<geometry.nX;++i) indice_x_Richards[i]=true;
+
+  double error=0;
+  Solution* Snew=allVerticalRichards(Sin,Sin,Spl);
+  if(Snew==nullptr) return false; //allVerti... must compute sol_out.P and sol_out.hsat
+  //cout<<"[Piccard] done"<<endl;
+  //Quitte de facon anticipe
+  return true;
+  compute_l(Snew,error);
+  compute_Pl(Snew,Snew.l,Snew.P);
+  if(!horizontalProblem(Sin,Snew,error)) return false; //Set hydr of Snew voir mettre hydr_in
+  if(!overlandProblem(Sin,Sin,Snew,error)) return false; //Set hov of Snew
+
+  //Comput error
+  double error_prev=error;
+
+  while(error>=tolerence_Piccard and count<max_iterations_Piccard and (abs(error_prev-error)>tolerence_Piccard/100 or error<oscilation_Piccard)){
+    Solution* Sold=Snew;
+    error_prev=error;
+    ++count;
+    error=0;
+    Snew=allVerticalRichards(Sin,Sold,Sold);
+    if(Snew==nullptr) return false;
+    compute_l(Snew,error);
+    compute_Pl(Snew,Snew.l,Snew.P);
+    if(!horizontalProblem(Sold,Snew,error)) return false;
+    if(!overlandProblem(Sin,Sold,Snew,error)) return false;
+  }
+
+  compute_mass(Snew);
+  Snew.nstep_Piccard=count;
+  //cout<<"[Piccard] done"<<endl;
+  if(error>=tolerence_Piccard){
+    return nullptr;
+    //Faire attention à la pile -> rempiler des solutions
+  }
+  return Snew;
+}
+
+
+void
+Kernel::compute_l(Solution& S,double& error) {
+  //Update S.l using S.hsat
+  bool e=0;
+  for(size_t ix=0;ix<geometry.nX;++ix){
+    double a=S.l[ix];
+    double b=max(S.hsat[ix],a);
+    S.l[ix]=b;
+    double t=b-a;
+    e+=t*t;
+  }
+  error+=sqrt(e);
+}
+
+void
+Kernel::compute_Pl(Solution& S,const double* h,double** P){
+  //cout<<"[Kernel::compute_Pl]"<<endl;
+  for(size_t ix=0;ix<geometry.nX;++ix){
+    double* Px=P[ix];
+    if(h[ix]==geometry.hsoil[ix]){
+      S.Pl[ix]=Px[geometry.nZ[ix]-1];
+      //if(S.Pl[ix]!=-2000) cout<<"error "<<ix<<endl;
+    }
+    else{
+      size_t a=(h[ix]-geometry.hbot[ix])/geometry.dZ[ix];
+      double p1=Px[a];
+      double p2=Px[a+1];
+      S.Pl[ix]=p1+(p2-p1)/geometry.dZ[ix]*(h[ix]-geometry.Z[ix][a]);
+      //cout<<ix<<" : "<<S.Pl[ix]<<endl;
+
+    }
+  }
+}
+
+bool
+Kernel::allVerticalRichards(Solution& S_0,Solution& S_s,Solution& S_sp){
+  //cout<<"[Kernel::allVerticalRichards]"<<endl;
+  for(size_t ix=0;ix<geometry.nX;++ix){
+    if(indice_x_Richards[ix]){
+      if(!Richards1dEvolutiveTime(ix,S_0,S_s,S_sp)) return false;
+    }
+  }
+  //cout<<"[Kernel::allVerticalRichards] done"<<endl;
+  return true;
+}
+
+bool
+Kernel::Richards1dEvolutiveTime(size_t ix,const Solution& S_0,const Solution& S_s,Solution& S_sp){
+  //cout<<"[Kernel::Richards1dEvolutiveTime]"<<endl;
+  //cout<<" ix = "<<ix<<endl;
+  //Warning flux_bot and div_w are always zero
+  double flux_bot=bottomFluxRichards(ix,S_s);
+  compute_div_w(ix,S_s);
+  size_t q=1;
+  size_t l=0;
+  bool conv;
+  //cout<<"Time::dT "<<Time::dT<<endl;
+  double dt=Time::dT/m*3600; //in second
+  while(l<q and q<=max_Richards_time_subdivisions){
+    //cout<<" l,q = "<<l<<','<<q<<endl;
+    if(l=0){
+      conv=Richards1d(ix,S_0,S_s,S_sp,dt,flux_bot);
+    }
+    else{
+      conv=Richards1d(ix,S_0,S_sp,S_sp,dt,flux_bot);
+    }
+    if(!conv){
+      q*=2;
+      dt/=2;
+      l=0;
+    }
+    else{
+      ++l;
+    }
+  }
+  //{cout<<"Continue ?"<<endl;char a;cin>>a;}
+
+  if(q>max_Richards_time_subdivisions) return false;
+  return true;
+}
+
+void
+Kernel::compute_mass(Solution& S){
+
+}
+
+bool
+Kernel::horizontalProblem(Solution& S,double& error){
+  //Compute hydr from P, l, Pl, sources, time
+  //return error otherwise
+  return true;
+}
+
+bool
+Kernel::overlandProblem(const Solution& S_in,Solution& S_out,double& error){
+  //Compute hov
+  //return error otherwise
+  return true;
+}
+
+void
+Kernel::saturationInterface(const Solution& S_in,Solution& S_out){
+//  Solution& S_in=*solution[sol_in];
+//  Solution& S_out=*solution[sol_out];
+//  for(size_t ix=0;ix<geometry.nX;++ix){
+//    double* Px_in=S_in.P[ix];
+//    size_t iz=0;
+//    for(;iz<geometry.nZ[ix] and Px_in[iz]>Psat;++iz);
+//    if(iz>0 and Px_in[iz]<=Psat){
+//      S_out.hsat[ix]=(Psat-Px_in[iz-1])*geometry.dZ[ix]/(Px_in[iz]-Px_in[iz-1])+geometry.Z[ix][iz-1];
+//    }
+//    else{
+//      S_out.hsat[ix]=(iz==0)?geometry.hbot[ix]:geometry.hsoil[ix];
+//    }
+//  }
+}
+
+double
+Kernel::bottomFluxRichards(size_t ix,const Solution& S){
+  return 0;
+}
+
+void
+Kernel::compute_div_w(size_t ix,const Solution& S){
+  double* div_w_ix=div_w[ix];
+  for(size_t iz=0;iz<geometry.nZ[ix];++iz){
+    div_w_ix[iz]=0;
+  }
+}
+
+bool
+Kernel::Richards1d(size_t ix,const Solution& S_0,const Solution& S_s,Solution& S_sp,double dt,double flux_bot){
+  //cout<<"[Richards1d]"<<endl;
+  size_t nZ=geometry.nZ[ix];
+  double* P0=S_0.P[ix];
+  double* Ps=S_s.P[ix];
+  double* Psp=S_sp.P[ix];
+  double n_P0=norm2(P0,nZ);
+  if(weightedRichards1d(ix,P0,Ps,Psp,dt,flux_bot,1,n_P0)) return true;
+  return weightedRichards1d(ix,P0,Ps,Psp,dt,flux_bot,0.5,n_P0);
+}
+
+bool
+Kernel::weightedRichards1d(size_t ix,double* P0,double* Ps,double* Psp,double dt,double flux_bot,double r,double n_P0){
+//  cout<<"[weightedRichards1d]"<<endl;
+  //cout<<"r = "<<r<<endl;
+  //cout<<"ix = "<<ix<<endl;
+  size_t nZ=geometry.nZ[ix];
+  double n0=norm2(P0,nZ);
+  double* Pi=Ps;
+  double* Po=P_temp[ix];
+  size_t count=0;
+  double e=+inf;
+  while(e>=tolerence_Richards and count<max_iterations_Richards){
+    ++count;
+    solveSystemRichards(ix,P0,Pi,Po,flux_bot,dt,count);
+    e=error2(Pi,Po,nZ)/n0;
+    //cout<<count<<" : "<<e<<" vs "<<tolerence_Richards<<endl;
+    if(count==1){
+      Pi=Po;
+      Po=Psp;
+    }
+    else{
+      swap(Pi,Po);
+    }
+  }
+  if(e<tolerence_Richards){
+    //cout<<"[converge]"<<endl;
+    memcpy(Psp,Pi,nZ*sizeof(double));
+    return true;
+  }
+  //cout<<"[don't converge]"<<endl;
+  return false;
+}
+
+void
+Kernel::solveSystemRichards(size_t ix,double* P0,double* P,double* Pout,double flux_bot,double dt,size_t count){
+  //TODO treat fpump
+  //cout<<"[solveSystemRichards]"<<endl;
+  //cout<<"P0 = "<<P0<<endl;
+  //cout<<"Pin = "<<P<<endl;
+  //cout<<"Pout = "<<Pout<<endl;
+
+  size_t nZ=geometry.nZ[ix];
+  assert(nZ>=3);
+  double dZ=geometry.dZ[ix];
+
+
+
+  double* sup_A=t_sup_A[ix];
+  double* diag_A=t_diag_A[ix];
+  double* sub_A=t_sub_A[ix];
+  double* sup_B=t_sup_B[ix];
+  double* diag_B=t_diag_B[ix];
+  double* sub_B=t_sub_B[ix];
+  double* sup_C=t_sup_C[ix];
+  double* diag_C=t_diag_C[ix];
+  double* sub_C=t_sub_C[ix];
+  double* F=t_F[ix];
+  double* G=t_G[ix];
+  double* S=t_S[ix];
+  double* H=t_H[ix];
+  double* R=t_R[ix];
+  double* I=t_I[ix];
+  double* BB=t_BB[ix];
+
+  //Compute A
+  diag_A[0]=(Physics::phi*dZ*Physics::ds(P[0]))/(2*dt);
+  for(size_t i=1;i<nZ-1;++i){
+    diag_A[i]=(Physics::phi*dZ*Physics::ds(P[i]))/dt;
+  }
+
+
+  //Compute B
+
+  //TODO : A optimiser boucle par rapport aux appels de Kr
+  double alpha=Physics::k0/(2*Physics::rho*Physics::g*dZ);
+  diag_B[0]=alpha*(Physics::kr(P[0])+Physics::kr(P[1]));
+  sup_B[0]=-diag_B[0];
+  diag_B[nZ-2]=alpha*(Physics::kr(P[nZ-3])+2*Physics::kr(P[nZ-2])+Physics::kr(P[nZ-1]));
+  sub_B[nZ-3]=-alpha*(Physics::kr(P[nZ-3])+Physics::kr(P[nZ-2]));
+  for(size_t i=1;i<nZ-2;++i){
+    diag_B[i]=alpha*(Physics::kr(P[i-1])+2*Physics::kr(P[i])+Physics::kr(P[i+1]));
+    sub_B[i-1]=-alpha*(Physics::kr(P[i-1])+Physics::kr(P[i]));
+    sup_B[i]=-alpha*(Physics::kr(P[i])+Physics::kr(P[i+1]));
+  }
+
+
+  //Compute C
+  double hk=Physics::k0/2;
+  double temp=1/(dZ*Physics::rho*Physics::g);
+  diag_C[0]=-hk*Physics::dkr(P[0])*(temp*(P[1]-P[0])+1);
+  sup_C[0]=-hk*Physics::dkr(P[1])*(temp*(P[1]-P[0])+1);
+  diag_C[nZ-2]=hk*Physics::dkr(P[nZ-3])*temp*(-P[nZ-3]+2*P[nZ-2]-P[nZ-1]);
+  sub_C[nZ-3]=hk*Physics::dkr(P[nZ-3])*(temp*(P[nZ-2]-P[nZ-3])+1);
+  for(size_t i=1;i<nZ-2;++i){
+    diag_C[i]=hk*Physics::dkr(P[i])*temp*(-P[i-1]+2*P[i]-P[i+1]);
+    sub_C[i-1]=hk*Physics::dkr(P[i-1])*(temp*(P[i]-P[i-1])+1);
+    sup_C[i]=-hk*Physics::dkr(P[i+1])*(temp*(P[i+1]-P[i])+1);
+  }
+
+  //Compute G
+  clear(G,nZ-1);
+  G[nZ-2]=-alpha*(Physics::kr(P[nZ-2])+Physics::kr(P[nZ-1]))*P[nZ-1];
+
+  //Compute S
+  S[0]=-hk*(Physics::kr(P[0])+Physics::kr(P[1]));
+  for(size_t i=1;i<nZ-1;++i){
+    S[i]=hk*(Physics::kr(P[i-1])-Physics::kr(P[i+1]));
+  }
+
+  //Compute H
+  clear(H,nZ-1);
+  H[nZ-2]=-hk*Physics::dkr(P[nZ-1])*(temp*(P[nZ-1]-P[nZ-2])+1);
+
+  //Compute R
+  clear(R,nZ-1);
+  R[0]=-flux_bot;
+
+  //Compute I
+  I[0]=(Physics::phi*dZ*(Physics::s(P[0])-Physics::s(P0[0])))/(2*dt);
+  for(size_t i=1;i<nZ-1;++i){
+    I[i]=(Physics::phi*dZ*(Physics::s(P[i])-Physics::s(P0[i])))/dt;
+  }
+
+  //Compute BB
+  clear(BB,nZ-1);
+  //TODO : Add BB computation from fpump
+
+  //Compute F
+  F[0]=div_w[ix][0]*dZ+R[0]-G[0]-I[0]-S[0]-H[0]+(diag_A[0]+diag_C[0])*P[0]+sup_C[0]*P[1];
+  for(size_t i=1;i<nZ-2;++i){
+    F[i]=div_w[ix][i]*dZ+R[i]-G[i]-I[i]-S[i]-H[i]+(diag_A[i]+diag_C[i])*P[i]+sub_C[i-1]*P[i-1]+sup_C[i]*P[i+1];
+  }
+  F[nZ-2]=div_w[ix][nZ-2]*dZ+R[nZ-2]-G[nZ-2]-I[nZ-2]-S[nZ-2]-H[nZ-2]+(diag_A[nZ-2]+diag_C[nZ-2])*P[nZ-2]+sub_C[nZ-3]*P[nZ-3];
+
+// /  if(ix==39)  display("F",F,nZ-1);
+  //TODO Add contribution of BB in F
+
+  for(size_t i=0;i<nZ-2;++i){
+    sub_A[i]=(sub_B[i]+sub_C[i]);
+    diag_A[i]+=(diag_B[i]+diag_C[i]);
+    sup_A[i]=(sup_B[i]+sup_C[i]);
+  }
+  diag_A[nZ-2]+=(diag_B[nZ-2]+diag_C[nZ-2]);
+
+
+
+ Thomas(nZ-1,sub_A,diag_A,sup_A,F,Pout);
+ Pout[nZ-1]=P[nZ-1];
+}
+*/

src/kernel.hpp

@@ -1,20 +1,85 @@
 #ifndef KERNEL_HPP
 #define KERNEL_HPP
-
+#include <cstdio>
+#include <cstring>
+#include <fstream>
+#include <stack>
 #include "physics.hpp"
 #include "geometry.hpp"
 #include "time.hpp"
+#include "input_data.hpp"
+#include "solution.hpp"
+#include "piccard.hpp"
+
+using namespace std;
+
+class Kernel:public InputData{
+private:
+  static constexpr double tolerence_Piccard=1e-6;
+  static constexpr double oscilation_Piccard=1e-4;
+  static constexpr size_t max_iterations_Piccard=50;
+  static constexpr size_t max_time_subdivisions=32;
+  static constexpr size_t max_Richards_time_subdivisions=256;
+  static constexpr size_t max_iterations_Richards=50;
+  static constexpr double tolerence_Richards=1e-10;
+
+  //For which colum we have to compute Richards
+  bool* indice_x_Richards;
+  int scheme_error;
+
 
 
-class Kernel{
+  void initSolutions();
+
+  void spaceSolution();
+
+  /*double bottomFluxRichards(size_t ix,const Solution& S_k);
+  Solution* Piccard(Solution* Sin);
+  bool Richards1dEvolutiveTime(size_t ix,const Solution& S_in,const Solution& S_k,Solution& S_kp);
+  bool allVerticalRichards(Solution& S_in,Solution& S_k,Solution& S_kp);
+  void compute_l(Solution& S,double& error);
+  void compute_Pl(Solution& S,const double* h,double** P);
+  void compute_mass(Solution& S);
+  void compute_div_w(size_t ix,const Solution& S);
+  bool horizontalProblem(Solution& S,double& error);
+  bool overlandProblem(const Solution& S_in,Solution& S_out,double& error);
+  void saturationInterface(const Solution& S_in,Solution& S_out);
+  bool Richards1d(size_t ix,const Solution& S_0,const Solution& S_s,Solution& S_sp,double dt,double flux_bot);
+  bool weightedRichards1d(size_t ix,double* P0,double* Ps,double* Psp,double dt,double flux_bot,double r,double n_P0);
+  void solveSystemRichards(size_t ix,double* P0,double* Pin,double* Pout,double flux_bot,double dt,size_t count);*/
+//  double computeError(const Solution& S_in,const Solution& S_out,size_t sol_temp);
 public:
-  Geometry& geometry;
   size_t step;
-  double*** P;
-  
-  size_t nZ(size_t ix);
-  double Z(size_t ix,size_t iz);
-  void next();
+  //Number of sub time intervals of [step,step+1]
+  size_t m;
+  stack<Solution*> pool_solutions;
+  Solution** solution;
+/*  double** div_w;
+  double** P_temp;
+
+  double** t_sub_A;
+  double** t_sub_B;
+  double** t_sub_C;
+  double** t_diag_A;
+  double** t_diag_B;
+  double** t_diag_C;
+  double** t_sup_A;
+  double** t_sup_B;
+  double** t_sup_C;
+  double** t_F;
+  double** t_G;
+  double** t_S;
+  double** t_H;
+  double** t_R;
+  double** t_I;
+  double** t_BB;*/
+  //size_t nZ(size_t ix);
+  //double Z(size_t ix,size_t iz);
+
+  Kernel(string filename);
+  void init(fstream& file);
+  bool next();
+
 };
 
 inline size_t

src/main.cpp

@@ -6,21 +6,7 @@
 
 using namespace std;
 
-double hsoil(double x){
-  return 3+cos(x);
-}
-
-double dhsoil(double x){
-  return -sin(x);
-}
 
-double hbot(double x){
-  return sin(x);
-}
-
-double dhbot(double x){
-  return cos(x);
-}
 
 int main(int argc,char** argv){
   QApplication app(argc,argv);
@@ -34,5 +20,5 @@ int main(int argc,char** argv){
   window->resize(1280,1024);
   window->show();
   return app.exec();
-  
+
 }

src/main_input.cpp

@@ -0,0 +1,15 @@
+#include <iostream>
+#include <QApplication>
+#include "qt/input/main_window.hpp"
+
+using namespace std;
+
+int main(int argc,char** argv){
+  QApplication app(argc,argv);
+  QtMainWindow* window;
+  window=new QtMainWindow;
+  window->resize(1280,1024);
+  window->show();
+  return app.exec();
+
+}

src/main_kernel.cpp

@@ -0,0 +1,13 @@
+#include <iostream>
+#include "kernel.hpp"
+
+
+using namespace std;
+
+int main(int argc,char** argv){
+  Kernel kernel("inputs/kernel-2.input");
+  for(size_t i=0;i<Time::nT;++i){
+    cout<<i<<" of "<<Time::nT<<endl;
+    kernel.next();
+  }
+}

src/main_output.cpp

@@ -0,0 +1,14 @@
+#include <iostream>
+#include <QApplication>
+#include "qt/output/main_window.hpp"
+
+using namespace std;
+
+int main(int argc,char** argv){
+  QApplication app(argc,argv);
+  QtMainWindow* window;
+  window=new QtMainWindow;
+  window->resize(1600,1000);
+  window->show();
+  return app.exec();
+}

src/math/algo.cpp

@@ -14,3 +14,39 @@ void Thomas(size_t n,double* a,double* b,double* c,double* d,double* x){
     last=x[i-1]=d[i-1]-c[i-1]*last;
   }
 }
+
+double
+norm2(double* u,size_t n){
+  double r=0;
+  for(size_t i=0;i<n;++i){
+    double t=u[i];
+    r+=(t*t);
+  }
+  return sqrt(r);
+}
+
+double
+error2(double* u,double* v,size_t n){
+  double r=0;
+  for(size_t i=0;i<n;++i){
+    double t=u[i]-v[i];
+    r+=(t*t);
+  }
+  return sqrt(r);
+}
+
+void
+clear(double* u,size_t n){
+  for(size_t i=0;i<n;++i){
+    u[i]=0;
+  }
+}
+
+void
+display(string name,double* u,size_t n){
+  cout<<"------| "<<name<<" |------"<<endl;
+  for(size_t i=0;i<n;++i){
+    cout<<i<<" : "<<u[i]<<endl;
+  }
+  cout<<"--------------------------"<<endl;
+}

src/math/algo.hpp

@@ -1,8 +1,12 @@
 #ifndef ALGO_HPP
 #define ALGO_HPP
 
+
+#include <cmath>
 #include <iostream>
+#include <string>
 
+using namespace std;
 /*! Thomas algorithm to solve AX=B where A is a tridiagonal matrix
   \param n size of the matrix A
   \param a (n-1) sub-diagonal coefficients of matrix A
@@ -13,5 +17,8 @@
 */
 
 void Thomas(size_t n,double* a,double* b,double* c,double* d,double* x);
-
+double norm2(double* u,size_t n);
+double error2(double* u,double* v,size_t n);
+void clear(double* u,size_t n);
+void display(string name,double* u,size_t n);
 #endif

src/math/spline.cpp

@@ -47,9 +47,9 @@ void Spline::compute(){
     q[i].b=m[i]/2;
     q[i].a=(m[i+1]-m[i])/(6*h[i]);
     q[i].c=(P[i+1].y-P[i].y)/h[i]-h[i]*(m[i+1]-m[i])/6-(h[i]*m[i])/2;
-    
+
   }
-  
+
   delete[] Asub;
   delete[] Adiag;
   delete[] Asup;
@@ -62,7 +62,9 @@ size_t
 Spline::findIndex(double x) const{
 //Find minimal i s.t. P[i].x>x
   size_t i=0;
-  while(i<n-1 and P[i].x<=x) ++i;
+  while(i<n-1 and P[i].x<=x){
+    ++i;
+  }
   // Here i=n-1 or P[i].x>x
   return i-1;
 }
@@ -70,7 +72,7 @@ Spline::findIndex(double x) const{
 double
 Spline::operator()(double x) const{
   size_t i=findIndex(x);
-  return q[i](x-P[i].x);	     
+  return q[i](x-P[i].x);
 }
 
 double

src/math/spline.hpp

@@ -2,6 +2,7 @@
 #define SPLINE_HPP
 
 #include <iostream>
+#include <cmath>
 #include "math/point.hpp"
 #include "math/poly3.hpp"
 #include "math/algo.hpp"

src/old/kernel.cpp

@@ -0,0 +1,489 @@
+#include "kernel.hpp"
+
+Kernel::Kernel(string filename){
+  fstream file;
+  file.open(filename.c_str(),fstream::in|fstream::binary);
+  InputData::load(file,true);
+
+  initSolutions();
+
+  Solution& S0=*solution[0];
+
+  //Init pressure
+  for(size_t x=0;x<geometry.nX;++x){
+    for(size_t z=0;z<geometry.nZ[x];++z){
+      S0.P[x][z]=initial_state->Pinit[x][z];
+    }
+  }
+  //Init hydraulic head
+  for(size_t x=0;x<geometry.nX;++x){
+    S0.hydr[x]=initial_state->hsat[x]+Psat/(Physics::rho*Physics::g);
+  }
+
+  //Init overland level
+  for(size_t x=0;x<geometry.nX;++x){
+    S0.hov[x]=Physics::invPsol(geometry.hsoil[x],S0.P[x][geometry.nZ[x]-1]); //TODO[Chrisophe] à initiliser
+  }
+
+  //Init l
+  for(size_t x=0;x<geometry.nX;++x){
+    double t=initial_state->hsat[x];
+    S0.l[x]=t;
+    S0.hsat[x]=t;
+  }
+
+  //Init Pl
+  for(size_t x=0;x<geometry.nX;++x){
+    S0.Pl[x]=Psat;
+  }
+
+
+  indice_x_Richards=new bool[geometry.nX];
+
+  div_w=new double*[geometry.nX];
+  P_temp=new double*[geometry.nX];
+  for(size_t i=0;i<geometry.nX;++i){
+    div_w[i]=new double[geometry.nZ[i]];
+    P_temp[i]=new double[geometry.nZ[i]];
+  }
+  t_sub_A=new double*[geometry.nX];
+  t_sub_B=new double*[geometry.nX];
+  t_sub_C=new double*[geometry.nX];
+  t_diag_A=new double*[geometry.nX];
+  t_diag_B=new double*[geometry.nX];
+  t_diag_C=new double*[geometry.nX];
+  t_sup_A=new double*[geometry.nX];
+  t_sup_B=new double*[geometry.nX];
+  t_sup_C=new double*[geometry.nX];
+  t_F=new double*[geometry.nX];
+  t_G=new double*[geometry.nX];
+  t_S=new double*[geometry.nX];
+  t_H=new double*[geometry.nX];
+  t_R=new double*[geometry.nX];
+  t_I=new double*[geometry.nX];
+  t_BB=new double*[geometry.nX];
+
+  for(size_t i=0;i<geometry.nX;++i){
+    size_t nZ=geometry.nZ[i];
+    t_sub_A[i]=new double[nZ-2];
+    t_sub_B[i]=new double[nZ-2];
+    t_sub_C[i]=new double[nZ-2];
+    t_diag_A[i]=new double[nZ-1];
+    t_diag_B[i]=new double[nZ-1];
+    t_diag_C[i]=new double[nZ-1];
+    t_sup_A[i]=new double[nZ-2];
+    t_sup_B[i]=new double[nZ-2];
+    t_sup_C[i]=new double[nZ-2];
+    t_F[i]=new double[nZ-1];
+    t_G[i]=new double[nZ-1];
+    t_S[i]=new double[nZ-1];
+    t_H[i]=new double[nZ-1];
+    t_R[i]=new double[nZ-1];
+    t_I[i]=new double[nZ-1];
+    t_BB[i]=new double[nZ-1];
+  }
+
+  delete initial_state;
+  initial_state=nullptr;
+
+  file.close();
+  step=0;
+  m=1;
+  scheme_error=0;
+}
+
+void
+Kernel::initSolutions(){
+  Solution* S0=new Solution();
+  S0->init(geometry);
+  solution[0]=S0;
+  for(size_t i=1;i<=Time::nT;++i) solution[i]=nullptr;
+  for(size_t i=0;i<Time::nT+2;++i){
+    Solution* S=new Solution();
+    S->init(geometry);
+    pool_solutions.push(S);
+  }
+}
+
+bool
+Kernel::next(){
+  //cout<<"[Kernel::next]"<<endl;
+  if(m>1) m=m/2;
+  spaceSolution();
+  if(scheme_error>0) return false;
+  ++step;
+  //cout<<"[next] done"<<endl;
+  return true;
+}
+
+void
+Kernel::spaceSolution(){
+  //cout<<"[Kernel::spaceSolution]"<<endl;
+  size_t n=step+1; //We want solution at time n
+  size_t k=0; //we work on time interval [n-1+k/m,n-1+(k+1)/m]
+  Solution *S=solution[step];
+  while(m<=max_time_subdivisions){
+    //cout<<" k,m = "<<k<<','<<m<<endl;
+    S=Piccard(S);
+    if(S==nullptr){
+      m*=2;
+      k=0;
+      S=solution[step];
+    }
+    else{
+      ++k;
+      if(k==m) break;
+    }
+  }
+  if(m>max_time_subdivisions){
+    cerr<<"Max time subdivision reached"<<endl;
+    scheme_error=1;
+  }
+  //cout<<"[spaceSolution] done"<<endl;
+}
+
+Solution*
+Kernel::Piccard(Solution* Sin){
+  //Return a valid solution or nullptr if not converge
+  //cout<<"[Kernel::Piccard]"<<endl;
+  memcpy(Sin.l,Sin.hsat,sizeof(double)*geometry.nX);
+  //Compute Pl with l=hsat of input solution and P=pressure of input_solution
+
+  Solution* Spl=pool_solutions.pop();
+  compute_Pl(*Spl,Sin.hsat,Sin.P);
+  double dt=Time::dT/m; // util ????
+
+  size_t count=0;
+
+  //Specicy to compute Richards on all columns
+  for(size_t i=0;i<geometry.nX;++i) indice_x_Richards[i]=true;
+
+  double error=0;
+  Solution* Snew=allVerticalRichards(Sin,Sin,Spl);
+  if(Snew==nullptr) return false; //allVerti... must compute sol_out.P and sol_out.hsat
+  //cout<<"[Piccard] done"<<endl;
+  //Quitte de facon anticipe
+  return true;
+  compute_l(Snew,error);
+  compute_Pl(Snew,Snew.l,Snew.P);
+  if(!horizontalProblem(Sin,Snew,error)) return false; //Set hydr of Snew voir mettre hydr_in
+  if(!overlandProblem(Sin,Sin,Snew,error)) return false; //Set hov of Snew
+
+  //Comput error
+  double error_prev=error;
+
+  while(error>=tolerence_Piccard and count<max_iterations_Piccard and (abs(error_prev-error)>tolerence_Piccard/100 or error<oscilation_Piccard)){
+    Solution* Sold=Snew;
+    error_prev=error;
+    ++count;
+    error=0;
+    Snew=allVerticalRichards(Sin,Sold,Sold);
+    if(Snew==nullptr) return false;
+    compute_l(Snew,error);
+    compute_Pl(Snew,Snew.l,Snew.P);
+    if(!horizontalProblem(Sold,Snew,error)) return false;
+    if(!overlandProblem(Sin,Sold,Snew,error)) return false;
+  }
+
+  compute_mass(Snew);
+  Snew.nstep_Piccard=count;
+  //cout<<"[Piccard] done"<<endl;
+  if(error>=tolerence_Piccard){
+    return nullptr;
+    //Faire attention à la pile -> rempiler des solutions
+  }
+  return Snew;
+}
+
+
+void
+Kernel::compute_l(Solution& S,double& error) {
+  //Update S.l using S.hsat
+  bool e=0;
+  for(size_t ix=0;ix<geometry.nX;++ix){
+    double a=S.l[ix];
+    double b=max(S.hsat[ix],a);
+    S.l[ix]=b;
+    double t=b-a;
+    e+=t*t;
+  }
+  error+=sqrt(e);
+}
+
+void
+Kernel::compute_Pl(Solution& S,const double* h,double** P){
+  //cout<<"[Kernel::compute_Pl]"<<endl;
+  for(size_t ix=0;ix<geometry.nX;++ix){
+    double* Px=P[ix];
+    if(h[ix]==geometry.hsoil[ix]){
+      S.Pl[ix]=Px[geometry.nZ[ix]-1];
+      //if(S.Pl[ix]!=-2000) cout<<"error "<<ix<<endl;
+    }
+    else{
+      size_t a=(h[ix]-geometry.hbot[ix])/geometry.dZ[ix];
+      double p1=Px[a];
+      double p2=Px[a+1];
+      S.Pl[ix]=p1+(p2-p1)/geometry.dZ[ix]*(h[ix]-geometry.Z[ix][a]);
+      //cout<<ix<<" : "<<S.Pl[ix]<<endl;
+
+    }
+  }
+}
+
+bool
+Kernel::allVerticalRichards(Solution& S_0,Solution& S_s,Solution& S_sp){
+  //cout<<"[Kernel::allVerticalRichards]"<<endl;
+  for(size_t ix=0;ix<geometry.nX;++ix){
+    if(indice_x_Richards[ix]){
+      if(!Richards1dEvolutiveTime(ix,S_0,S_s,S_sp)) return false;
+    }
+  }
+  //cout<<"[Kernel::allVerticalRichards] done"<<endl;
+  return true;
+}
+
+bool
+Kernel::Richards1dEvolutiveTime(size_t ix,const Solution& S_0,const Solution& S_s,Solution& S_sp){
+  //cout<<"[Kernel::Richards1dEvolutiveTime]"<<endl;
+  //cout<<" ix = "<<ix<<endl;
+  //Warning flux_bot and div_w are always zero
+  double flux_bot=bottomFluxRichards(ix,S_s);
+  compute_div_w(ix,S_s);
+  size_t q=1;
+  size_t l=0;
+  bool conv;
+  //cout<<"Time::dT "<<Time::dT<<endl;
+  double dt=Time::dT/m*3600; //in second
+  while(l<q and q<=max_Richards_time_subdivisions){
+    //cout<<" l,q = "<<l<<','<<q<<endl;
+    if(l=0){
+      conv=Richards1d(ix,S_0,S_s,S_sp,dt,flux_bot);
+    }
+    else{
+      conv=Richards1d(ix,S_0,S_sp,S_sp,dt,flux_bot);
+    }
+    if(!conv){
+      q*=2;
+      dt/=2;
+      l=0;
+    }
+    else{
+      ++l;
+    }
+  }
+  //{cout<<"Continue ?"<<endl;char a;cin>>a;}
+
+  if(q>max_Richards_time_subdivisions) return false;
+  return true;
+}
+
+void
+Kernel::compute_mass(Solution& S){
+
+}
+
+bool
+Kernel::horizontalProblem(Solution& S,double& error){
+  //Compute hydr from P, l, Pl, sources, time
+  //return error otherwise
+  return true;
+}
+
+bool
+Kernel::overlandProblem(const Solution& S_in,Solution& S_out,double& error){
+  //Compute hov
+  //return error otherwise
+  return true;
+}
+
+void
+Kernel::saturationInterface(const Solution& S_in,Solution& S_out){
+/*  Solution& S_in=*solution[sol_in];
+  Solution& S_out=*solution[sol_out];
+  for(size_t ix=0;ix<geometry.nX;++ix){
+    double* Px_in=S_in.P[ix];
+    size_t iz=0;
+    for(;iz<geometry.nZ[ix] and Px_in[iz]>Psat;++iz);
+    if(iz>0 and Px_in[iz]<=Psat){
+      S_out.hsat[ix]=(Psat-Px_in[iz-1])*geometry.dZ[ix]/(Px_in[iz]-Px_in[iz-1])+geometry.Z[ix][iz-1];
+    }
+    else{
+      S_out.hsat[ix]=(iz==0)?geometry.hbot[ix]:geometry.hsoil[ix];
+    }
+  }*/
+}
+
+double
+Kernel::bottomFluxRichards(size_t ix,const Solution& S){
+  return 0;
+}
+
+void
+Kernel::compute_div_w(size_t ix,const Solution& S){
+  double* div_w_ix=div_w[ix];
+  for(size_t iz=0;iz<geometry.nZ[ix];++iz){
+    div_w_ix[iz]=0;
+  }
+}
+
+bool
+Kernel::Richards1d(size_t ix,const Solution& S_0,const Solution& S_s,Solution& S_sp,double dt,double flux_bot){
+  //cout<<"[Richards1d]"<<endl;
+  size_t nZ=geometry.nZ[ix];
+  double* P0=S_0.P[ix];
+  double* Ps=S_s.P[ix];
+  double* Psp=S_sp.P[ix];
+  double n_P0=norm2(P0,nZ);
+  if(weightedRichards1d(ix,P0,Ps,Psp,dt,flux_bot,1,n_P0)) return true;
+  return weightedRichards1d(ix,P0,Ps,Psp,dt,flux_bot,0.5,n_P0);
+}
+
+bool
+Kernel::weightedRichards1d(size_t ix,double* P0,double* Ps,double* Psp,double dt,double flux_bot,double r,double n_P0){
+//  cout<<"[weightedRichards1d]"<<endl;
+  //cout<<"r = "<<r<<endl;
+  //cout<<"ix = "<<ix<<endl;
+  size_t nZ=geometry.nZ[ix];
+  double n0=norm2(P0,nZ);
+  double* Pi=Ps;
+  double* Po=P_temp[ix];
+  size_t count=0;
+  double e=+inf;
+  while(e>=tolerence_Richards and count<max_iterations_Richards){
+    ++count;
+    solveSystemRichards(ix,P0,Pi,Po,flux_bot,dt,count);
+    e=error2(Pi,Po,nZ)/n0;
+    //cout<<count<<" : "<<e<<" vs "<<tolerence_Richards<<endl;
+    if(count==1){
+      Pi=Po;
+      Po=Psp;
+    }
+    else{
+      swap(Pi,Po);
+    }
+  }
+  if(e<tolerence_Richards){
+    //cout<<"[converge]"<<endl;
+    memcpy(Psp,Pi,nZ*sizeof(double));
+    return true;
+  }
+  //cout<<"[don't converge]"<<endl;
+  return false;
+}
+
+void
+Kernel::solveSystemRichards(size_t ix,double* P0,double* P,double* Pout,double flux_bot,double dt,size_t count){
+  //TODO treat fpump
+  /*cout<<"[solveSystemRichards]"<<endl;
+  cout<<"P0 = "<<P0<<endl;
+  cout<<"Pin = "<<P<<endl;
+  cout<<"Pout = "<<Pout<<endl;*/
+
+  size_t nZ=geometry.nZ[ix];
+  assert(nZ>=3);
+  double dZ=geometry.dZ[ix];
+
+
+
+  double* sup_A=t_sup_A[ix];
+  double* diag_A=t_diag_A[ix];
+  double* sub_A=t_sub_A[ix];
+  double* sup_B=t_sup_B[ix];
+  double* diag_B=t_diag_B[ix];
+  double* sub_B=t_sub_B[ix];
+  double* sup_C=t_sup_C[ix];
+  double* diag_C=t_diag_C[ix];
+  double* sub_C=t_sub_C[ix];
+  double* F=t_F[ix];
+  double* G=t_G[ix];
+  double* S=t_S[ix];
+  double* H=t_H[ix];
+  double* R=t_R[ix];
+  double* I=t_I[ix];
+  double* BB=t_BB[ix];
+
+  //Compute A
+  diag_A[0]=(Physics::phi*dZ*Physics::ds(P[0]))/(2*dt);
+  for(size_t i=1;i<nZ-1;++i){
+    diag_A[i]=(Physics::phi*dZ*Physics::ds(P[i]))/dt;
+  }
+
+
+  //Compute B
+
+  //TODO : A optimiser boucle par rapport aux appels de Kr
+  double alpha=Physics::k0/(2*Physics::rho*Physics::g*dZ);
+  diag_B[0]=alpha*(Physics::kr(P[0])+Physics::kr(P[1]));
+  sup_B[0]=-diag_B[0];
+  diag_B[nZ-2]=alpha*(Physics::kr(P[nZ-3])+2*Physics::kr(P[nZ-2])+Physics::kr(P[nZ-1]));
+  sub_B[nZ-3]=-alpha*(Physics::kr(P[nZ-3])+Physics::kr(P[nZ-2]));
+  for(size_t i=1;i<nZ-2;++i){
+    diag_B[i]=alpha*(Physics::kr(P[i-1])+2*Physics::kr(P[i])+Physics::kr(P[i+1]));
+    sub_B[i-1]=-alpha*(Physics::kr(P[i-1])+Physics::kr(P[i]));
+    sup_B[i]=-alpha*(Physics::kr(P[i])+Physics::kr(P[i+1]));
+  }
+
+
+  //Compute C
+  double hk=Physics::k0/2;
+  double temp=1/(dZ*Physics::rho*Physics::g);
+  diag_C[0]=-hk*Physics::dkr(P[0])*(temp*(P[1]-P[0])+1);
+  sup_C[0]=-hk*Physics::dkr(P[1])*(temp*(P[1]-P[0])+1);
+  diag_C[nZ-2]=hk*Physics::dkr(P[nZ-3])*temp*(-P[nZ-3]+2*P[nZ-2]-P[nZ-1]);
+  sub_C[nZ-3]=hk*Physics::dkr(P[nZ-3])*(temp*(P[nZ-2]-P[nZ-3])+1);
+  for(size_t i=1;i<nZ-2;++i){
+    diag_C[i]=hk*Physics::dkr(P[i])*temp*(-P[i-1]+2*P[i]-P[i+1]);
+    sub_C[i-1]=hk*Physics::dkr(P[i-1])*(temp*(P[i]-P[i-1])+1);
+    sup_C[i]=-hk*Physics::dkr(P[i+1])*(temp*(P[i+1]-P[i])+1);
+  }
+
+  //Compute G
+  clear(G,nZ-1);
+  G[nZ-2]=-alpha*(Physics::kr(P[nZ-2])+Physics::kr(P[nZ-1]))*P[nZ-1];
+
+  //Compute S
+  S[0]=-hk*(Physics::kr(P[0])+Physics::kr(P[1]));
+  for(size_t i=1;i<nZ-1;++i){
+    S[i]=hk*(Physics::kr(P[i-1])-Physics::kr(P[i+1]));
+  }
+
+  //Compute H
+  clear(H,nZ-1);
+  H[nZ-2]=-hk*Physics::dkr(P[nZ-1])*(temp*(P[nZ-1]-P[nZ-2])+1);
+
+  //Compute R
+  clear(R,nZ-1);
+  R[0]=-flux_bot;
+
+  //Compute I
+  I[0]=(Physics::phi*dZ*(Physics::s(P[0])-Physics::s(P0[0])))/(2*dt);
+  for(size_t i=1;i<nZ-1;++i){
+    I[i]=(Physics::phi*dZ*(Physics::s(P[i])-Physics::s(P0[i])))/dt;
+  }
+
+  //Compute BB
+  clear(BB,nZ-1);
+  //TODO : Add BB computation from fpump
+
+  //Compute F
+  F[0]=div_w[ix][0]*dZ+R[0]-G[0]-I[0]-S[0]-H[0]+(diag_A[0]+diag_C[0])*P[0]+sup_C[0]*P[1];
+  for(size_t i=1;i<nZ-2;++i){
+    F[i]=div_w[ix][i]*dZ+R[i]-G[i]-I[i]-S[i]-H[i]+(diag_A[i]+diag_C[i])*P[i]+sub_C[i-1]*P[i-1]+sup_C[i]*P[i+1];
+  }
+  F[nZ-2]=div_w[ix][nZ-2]*dZ+R[nZ-2]-G[nZ-2]-I[nZ-2]-S[nZ-2]-H[nZ-2]+(diag_A[nZ-2]+diag_C[nZ-2])*P[nZ-2]+sub_C[nZ-3]*P[nZ-3];
+
+// /  if(ix==39)  display("F",F,nZ-1);
+  //TODO Add contribution of BB in F
+
+  for(size_t i=0;i<nZ-2;++i){
+    sub_A[i]=(sub_B[i]+sub_C[i]);
+    diag_A[i]+=(diag_B[i]+diag_C[i]);
+    sup_A[i]=(sup_B[i]+sup_C[i]);
+  }
+  diag_A[nZ-2]+=(diag_B[nZ-2]+diag_C[nZ-2]);
+
+
+
+ Thomas(nZ-1,sub_A,diag_A,sup_A,F,Pout);
+ Pout[nZ-1]=P[nZ-1];
+}

src/old/kernel.hpp

@@ -0,0 +1,94 @@
+#ifndef KERNEL_HPP
+#define KERNEL_HPP
+#include <cstdio>
+#include <cstring>
+#include <fstream>
+#include <stack>
+#include "physics.hpp"
+#include "geometry.hpp"
+#include "time.hpp"
+#include "input_data.hpp"
+#include "solution.hpp"
+#include "piccard.hpp"
+
+using namespace std;
+
+class Kernel:public InputData{
+private:
+  static constexpr double tolerence_Piccard=1e-6;
+  static constexpr double oscilation_Piccard=1e-4;
+  static constexpr size_t max_iterations_Piccard=50;
+  static constexpr size_t max_time_subdivisions=32;
+  static constexpr size_t max_Richards_time_subdivisions=256;
+  static constexpr size_t max_iterations_Richards=50;
+  static constexpr double tolerence_Richards=1e-10;
+
+  //For which colum we have to compute Richards
+  bool* indice_x_Richards;
+  int scheme_error;
+
+
+
+  void initSolutions();
+
+  void spaceSolution();
+
+  double bottomFluxRichards(size_t ix,const Solution& S_k);
+  Solution* Piccard(Solution* Sin);
+  bool Richards1dEvolutiveTime(size_t ix,const Solution& S_in,const Solution& S_k,Solution& S_kp);
+  bool allVerticalRichards(Solution& S_in,Solution& S_k,Solution& S_kp);
+  void compute_l(Solution& S,double& error);
+  void compute_Pl(Solution& S,const double* h,double** P);
+  void compute_mass(Solution& S);
+  void compute_div_w(size_t ix,const Solution& S);
+  bool horizontalProblem(Solution& S,double& error);
+  bool overlandProblem(const Solution& S_in,Solution& S_out,double& error);
+  void saturationInterface(const Solution& S_in,Solution& S_out);
+  bool Richards1d(size_t ix,const Solution& S_0,const Solution& S_s,Solution& S_sp,double dt,double flux_bot);
+  bool weightedRichards1d(size_t ix,double* P0,double* Ps,double* Psp,double dt,double flux_bot,double r,double n_P0);
+  void solveSystemRichards(size_t ix,double* P0,double* Pin,double* Pout,double flux_bot,double dt,size_t count);
+//  double computeError(const Solution& S_in,const Solution& S_out,size_t sol_temp);
+public:
+  size_t step;
+  //Number of sub time intervals of [step,step+1]
+  size_t m;
+  Solution* solution;
+  double** div_w;
+  double** P_temp;
+
+  double** t_sub_A;
+  double** t_sub_B;
+  double** t_sub_C;
+  double** t_diag_A;
+  double** t_diag_B;
+  double** t_diag_C;
+  double** t_sup_A;
+  double** t_sup_B;
+  double** t_sup_C;
+  double** t_F;
+  double** t_G;
+  double** t_S;
+  double** t_H;
+  double** t_R;
+  double** t_I;
+  double** t_BB;
+  size_t nZ(size_t ix);
+  double Z(size_t ix,size_t iz);
+
+  Kernel(string filename);
+  void init(fstream& file);
+  bool next();
+
+};
+
+inline size_t
+Kernel::nZ(size_t ix){
+  return geometry.nZ[ix];
+}
+
+inline double
+Kernel::Z(size_t ix,size_t iz){
+  return geometry.Z[ix][iz];
+}
+
+#endif

src/physics.cpp

@@ -3,98 +3,147 @@
 // Brooks and Corey model
 //------------------------
 
-namespace Physics{
-  double g=9.81;
-  double rho=1000;
-  double phi=0.3;
-  double k0=3e-5;
-  double nivrivsat=0.01;
-  double Psec=-10000;
-  double (*s)(double)=&s_BC;
-  double (*ds)(double)=&ds_BC;
-  void (*s_ds)(double,double&,double&)=&s_ds_BC;
-  double (*s_inv)(double)=&s_inv_BC;
-  double (*kr)(double)=&kr_BC;
-  double (*dkr)(double)=&dkr_BC;
-  void (*kr_dkr)(double,double&,double&)=&kr_dkr_BC;
-  double model_data[6]={-2000,0,3,11,0,0};
-
-#define Psat model_data[0]
-#define sres model_data[1]
-#define lambda model_data[2]
-#define alpha model_data[3]
-
-  void setModel(Model model){
-    switch(model){
+double Physics::g;
+double Physics::rho;
+double Physics::phi;
+double Physics::k0;
+double Physics::nivrivsat;
+double Physics::Psec;
+Physics::Model Physics::model;
+double (*Physics::s)(double);
+double (*Physics::ds)(double);
+void (*Physics::s_ds)(double,double&,double&);
+double (*Physics::s_inv)(double);
+double (*Physics::kr)(double);
+double (*Physics::dkr)(double);
+void (*Physics::kr_dkr)(double,double&,double&);
+double Physics::model_data[6];
+
+
+#define sres Physics::model_data[1]
+#define lambda Physics::model_data[2]
+#define alpha Physics::model_data[3]
+
+void Physics::setModel(Model m){
+  model=m;
+  switch(model){
     case BrooksCorey:
-      s=&s_BC;
-      ds=&ds_BC;
-      s_ds=&s_ds_BC;
-      kr=&kr_BC;
-      dkr=&dkr_BC;
-      kr_dkr=&kr_dkr_BC;
-      Psat=-2000;
-      sres=0;
-      lambda=3;
-      alpha=11;
-      break;
+    cout<<"Set model BC"<<endl;
+    s=&s_BC;
+    ds=&ds_BC;
+    s_ds=&s_ds_BC;
+    s_inv=&s_inv_BC;
+    kr=&kr_BC;
+    dkr=&dkr_BC;
+    cout<<"dkr : "<<&dkr<<endl;
+    kr_dkr=&kr_dkr_BC;
+    break;
     default:
-      assert(false);
-    };
-  }
+    assert(false);
+  };
+}
 
-  double
-  s_BC(double P){
-    if(P>=Psat) return 1;
-    return sres+(1-sres)*pow(Psat/P,lambda);
-  }
+double
+Physics::s_BC(double P){
+  if(P>=Psat) return 1;
+  return sres+(1-sres)*pow(Psat/P,lambda);
+}
 
-  double
-  ds_BC(double P){
-    if(P>=Psat) return 0;
-    return ((sres-1)*lambda*pow(Psat/P,lambda))/P;
-  }
+double
+Physics::ds_BC(double P){
+  if(P>=Psat) return 0;
+  return ((sres-1)*lambda*pow(Psat/P,lambda))/P;
+}
 
-  void
-  s_ds_BC(double P,double& v,double& dv){
-    if(P>=Psat){
-      v=1;
-      dv=0;
-    }
-    else{
-      double t=(1-sres)*pow(Psat/P,lambda);
-      v=sres+t;
-      dv=-(lambda*t)/P;
-    }
+void
+Physics::s_ds_BC(double P,double& v,double& dv){
+  if(P>=Psat){
+    v=1;
+    dv=0;
   }
-
-  double
-  s_inv_BC(double S){
-    return pow((S-sres)/(1-sres),-1/lambda)*Psat;
+  else{
+    double t=(1-sres)*pow(Psat/P,lambda);
+    v=sres+t;
+    dv=-(lambda*t)/P;
   }
+}
+
+double
+Physics::s_inv_BC(double S){
+  return pow((S-sres)/(1-sres),-1/lambda)*Psat;
+}
+
+double
+Physics::kr_BC(double P){
+  if(P>=Psat) return 1;
+  return pow(Psat/P,alpha);
+}
+
+double
+Physics::dkr_BC(double P){
+  if(P>=Psat) return 0;
+  return -alpha*pow(Psat/P,alpha)/P;
+}
 
-  double
-  kr_BC(double P){
-    if(P>=Psat) return 1;
-    return pow(Psat/P,alpha);
+void
+Physics::kr_dkr_BC(double P,double& v,double& dv){
+  if(P>=Psat){
+    v=1;
+    dv=0;
   }
+  else{
+    double t=pow(Psat/P,alpha);
+    v=t;
+    dv=-(alpha*t)/P;
+  }
+}
 
-  double
-  dkr_BC(double P){
-    if(P>=Psat) return 0;
-    return -alpha*pow(Psat/P,alpha)/P;
+void
+Physics::save(fstream& file){
+  file.write((char*)&g,sizeof(double));
+  file.write((char*)&rho,sizeof(double));
+  file.write((char*)&phi,sizeof(double));
+  file.write((char*)&k0,sizeof(double));
+  file.write((char*)&Psec,sizeof(double));
+  file.write((char*)&nivrivsat,sizeof(double));
+  file.write((char*)&model,sizeof(Model));
+  for(size_t i=0;i<max_model_parameters;++i){
+    file.write((char*)&model_data[i],sizeof(double));
   }
+}
 
-  void
-  kr_dkr_BC(double P,double& v,double& dv){
-    if(P>=Psat){
-      v=1;
-      dv=0;
-    }
-    else{
-      double t=pow(Psat/P,alpha);
-      v=t;
-    dv=-(alpha*t)/P;
-    }
+void
+Physics::load(fstream& file){
+  file.read((char*)&g,sizeof(double));
+  file.read((char*)&rho,sizeof(double));
+  file.read((char*)&phi,sizeof(double));
+  file.read((char*)&k0,sizeof(double));
+  file.read((char*)&Psec,sizeof(double));
+  file.read((char*)&nivrivsat,sizeof(double));
+  Model m;
+  file.read((char*)&m,sizeof(Model));
+  setModel(m);
+  for(size_t i=0;i<max_model_parameters;++i){
+    file.read((char*)&model_data[i],sizeof(double));
   }
 }
+
+double
+Physics::Psol(double hsoil,double hov){
+  double a=Psat-rho*g*nivrivsat;
+  return rho*g*(hov-hsoil)+a*nivrivsat/(hov-hsoil);
+}
+
+double
+Physics::dPsol(double hsoil,double hov){
+  double a=Psat-rho*g*nivrivsat;
+  double b=hov-hsoil;
+  return rho*g-a*nivrivsat/(b*b);
+}
+
+double
+Physics::invPsol(double hsoil,double Psol){
+  double a=Psat-rho*g*nivrivsat;
+  double b=a*nivrivsat;
+  return hsoil+(Psol+sqrt(Psol*Psol-4*rho*g*b))/(2*rho*g);
+}

src/physics.hpp

@@ -3,52 +3,62 @@
 
 #include <cmath>
 #include <cassert>
+#include <iostream>
+#include <fstream>
 
-//! The Physics namespace contains all physical parameters characterising the soil.
-namespace Physics{
+using namespace std;
 
+//! The Physics class contains all physical parameters characterising the soil.
+//! This class contains only static members.
+class Physics{
+public:
   enum Model{BrooksCorey};
 
   //! Set physics model
-  void setModel(Model model);
+  static void setModel(Model model);
 
   //! Gravity acceleration (m.s^-2)
-  extern double g;
+  static double g;
   //! Fluid density (g.l^(-1))
-  extern double rho;
+  static double rho;
 
   //! Porosity of the soil
-  extern double phi;
+  static double phi;
   //! Conductivity of the saturated soil
-  extern double k0;
+  static double k0;
 
-  extern double Psec;
+  static double Psec;
   //! Characterise the water pressure at the bottom of the overland water
-  extern double nivrivsat;
+  static double nivrivsat;
 
-   //! Return the saturation in function of the pressure
-  extern double (*s)(double);
+  //! Model used
+  static Model model;
+
+  //! Return the saturation in function of the pressure
+  static double (*s)(double);
   //! Return the derivtive of the saturation in function of the pressure
-  extern double (*ds)(double);
+  static double (*ds)(double);
   //! Set the saturation and its derivative in function of the pressure
-  extern void (*s_ds)(double,double&,double&);
+  static void (*s_ds)(double,double&,double&);
 
-  extern double (*s_inv)(double);
+  static double (*s_inv)(double);
 
   //! Return the relative conductivity in function of the pressure
-  extern double(*kr)(double);
+  static double (*kr)(double);
   //! Return the derivtive of the relative conductivity in function of the pressure
-  extern double (*dkr)(double);
+  static double (*dkr)(double);
   //! Set the relative conductivity and its derivative in function of the pressure
-  extern void (*kr_dkr)(double,double&,double&);
-
+  static void (*kr_dkr)(double,double&,double&);
 
   //---------------------
   // Models descriptions
   //---------------------
 
   //! Datas used to define the model
-  extern double model_data[6];
+  static double model_data[6];
+
+  //! Max model Datas
+  static const size_t max_model_parameters=6;
 
   //------------------------
   // Brooks and Corey model
@@ -61,19 +71,32 @@ namespace Physics{
   //model_data[3] -> alpha
 
   //! Brooks and Corey saturation map
-  double s_BC(double P);
+  static double s_BC(double P);
   //! Brooks and Corey derivative of the saturation map
-  double ds_BC(double P);
+  static double ds_BC(double P);
   //! Brooks and Corey saturation and its derivative setter
-  void (s_ds_BC)(double P,double& v,double& dv);
+  static void (s_ds_BC)(double P,double& v,double& dv);
 
-  double s_inv_BC(double s);
+  static double s_inv_BC(double s);
 
   //! Brooks and Corey relative conductivity map
-  double kr_BC(double P);
+  static double kr_BC(double P);
   //! Brooks and Corey derivative of the relative conductivity map
-  double dkr_BC(double P);
+  static double dkr_BC(double P);
   //! Brooks and Corey relative conductivity and its derivative setter
-  void (kr_dkr_BC)(double P,double& v,double& dv);
-}
+  static void (kr_dkr_BC)(double P,double& v,double& dv);
+
+  static double Psol(double hsoil,double hov);
+  static double dPsol(double hsoil,double hov);
+  static double invPsol(double hsoil,double Psol);
+  //--------------------
+  // File input/output
+  //--------------------
+
+  static void save(fstream& file);
+  static void load(fstream& file);
+
+};
+
+#define Psat Physics::model_data[0]
 #endif

src/piccard.hpp

@@ -0,0 +1,11 @@
+#ifndef PICCARD_HPP
+#define PICCARD_HPP
+
+class Piccard{
+private:
+public:
+  Piccard(Geometry& G);
+  bool eval(const Solution& S_in,Solution& S_out);
+};
+
+#endif

src/qt/.#input.hpp

@@ -1 +0,0 @@
-fromentin@atuan.36324:1624433263

src/qt/input_cloud.cpp

@@ -1,9 +1,9 @@
-#include "input_cloud.hpp"
+#include "qt/input/cloud.hpp"
 
-QtInputCloud::QtInputCloud(Source* _source){
-  source=_source;
-  cloud=source->addCloud();
+QtInputCloud::QtInputCloud(QtInputData* d,Cloud* c){
+  data=d;
 
+  cloud=(c==nullptr)?data->addCloud():c;
 
   init_groupbox=new QGroupBox("Initial values");
   amplitude_init_widget=new QWidget;
@@ -25,19 +25,19 @@ QtInputCloud::QtInputCloud(Source* _source){
   left_right_final_layout=new QHBoxLayout;
 
   amplitude_init_label=new QLabel("Amplitude: ");
-  amplitude_init_input=new QLineEdit(QString::number(cloud->amplitude_init));
+  amplitude_init_input=new QLineEdit();
   amplitude_init_layout->addWidget(amplitude_init_label);
   amplitude_init_layout->addWidget(amplitude_init_input);
   amplitude_init_widget->setLayout(amplitude_init_layout);
 
   left_init_label=new QLabel("Left: ");
   delta_left_init_label=new QLabel("  Left delta: ");
-  left_init_input=new QLineEdit(QString::number(cloud->left_init));
-  delta_left_init_input=new QLineEdit(QString::number(cloud->delta_left_init));
+  left_init_input=new QLineEdit();
+  delta_left_init_input=new QLineEdit();
   right_init_label=new QLabel("  Right: ");
   delta_right_init_label=new QLabel("  Right delta: ");
-  right_init_input=new QLineEdit(QString::number(cloud->right_init));
-  delta_right_init_input=new QLineEdit(QString::number(cloud->delta_right_init));
+  right_init_input=new QLineEdit();
+  delta_right_init_input=new QLineEdit();
 
   left_right_init_layout->addWidget(left_init_label);
   left_right_init_layout->addWidget(left_init_input,1);
@@ -53,21 +53,20 @@ QtInputCloud::QtInputCloud(Source* _source){
   init_layout->addWidget(left_right_init_widget);
   init_groupbox->setLayout(init_layout);
 
-
   amplitude_final_label=new QLabel("Amplitude: ");
-  amplitude_final_input=new QLineEdit(QString::number(cloud->amplitude_final));
+  amplitude_final_input=new QLineEdit();
   amplitude_final_layout->addWidget(amplitude_final_label);
   amplitude_final_layout->addWidget(amplitude_final_input);
   amplitude_final_widget->setLayout(amplitude_final_layout);
 
   left_final_label=new QLabel("Left: ");
   delta_left_final_label=new QLabel("  Left delta: ");
-  left_final_input=new QLineEdit(QString::number(cloud->left_final));
-  delta_left_final_input=new QLineEdit(QString::number(cloud->delta_left_final));
+  left_final_input=new QLineEdit();
+  delta_left_final_input=new QLineEdit();
   right_final_label=new QLabel("  Right: ");
   delta_right_final_label=new QLabel("  Right delta: ");
-  right_final_input=new QLineEdit(QString::number(cloud->right_final));
-  delta_right_final_input=new QLineEdit(QString::number(cloud->delta_right_final));
+  right_final_input=new QLineEdit();
+  delta_right_final_input=new QLineEdit();
 
   left_right_final_layout->addWidget(left_final_label);
   left_right_final_layout->addWidget(left_final_input,1);
@@ -113,6 +112,8 @@ QtInputCloud::QtInputCloud(Source* _source){
   setFrameShape(QFrame::Box);
 
   connect(remove_button,&QPushButton::clicked,this,&QtInputCloud::emitRemove);
+
+  getCloud();
 }
 
 
@@ -136,15 +137,36 @@ QtInputCloud::validate(){
   rf=right_final_input->text().toDouble();
   if(li>=ri) return right_init_input;
   if(lf>=rf) return right_final_input;
-  cloud->left_init=li;
-  cloud->right_init=ri;
-  cloud->left_final=lf;
-  cloud->right_final=rf;
+
+  setCloud();
+
+  return nullptr;
+}
+
+void
+QtInputCloud::setCloud(){
+  cloud->left_init=left_init_input->text().toDouble();
+  cloud->right_init=right_init_input->text().toDouble();
+  cloud->left_final=left_final_input->text().toDouble();
+  cloud->right_final=right_final_input->text().toDouble();
   cloud->amplitude_init=amplitude_init_input->text().toDouble();
   cloud->amplitude_final=amplitude_final_input->text().toDouble();
   cloud->delta_left_init=delta_left_init_input->text().toDouble();
   cloud->delta_right_init=delta_right_init_input->text().toDouble();
   cloud->delta_left_final=delta_left_final_input->text().toDouble();
   cloud->delta_right_final=delta_right_final_input->text().toDouble();
-  return nullptr;
+}
+
+void
+QtInputCloud::getCloud(){
+  amplitude_init_input->setText(QString::number(cloud->amplitude_init));
+  left_init_input->setText(QString::number(cloud->left_init));
+  right_init_input->setText(QString::number(cloud->right_init));
+  delta_left_init_input->setText(QString::number(cloud->delta_left_init));
+  delta_right_init_input->setText(QString::number(cloud->delta_right_init));
+  amplitude_final_input->setText(QString::number(cloud->amplitude_final));
+  left_final_input->setText(QString::number(cloud->left_final));
+  right_final_input->setText(QString::number(cloud->right_final));
+  delta_left_final_input->setText(QString::number(cloud->delta_left_final));
+  delta_right_final_input->setText(QString::number(cloud->delta_right_final));
 }

src/qt/input_cloud.hpp

@@ -11,13 +11,14 @@
 #include <QGroupBox>
 #include <QDoubleValidator>
 
-#include "source.hpp"
+#include "qt/input/data.hpp"
 
 class QtInputCloud:public QFrame{
   Q_OBJECT
 private:
-  Source* source;
+  QtInputData* data;
   Cloud* cloud;
+
   QGroupBox* init_groupbox;
   QGroupBox* final_groupbox;
   QWidget* init_final_widget;
@@ -64,9 +65,11 @@ private:
   QDoubleValidator* double_validator;
   QDoubleValidator* double_amplitude_validator;
 public:
-  QtInputCloud(Source*);
+  QtInputCloud(QtInputData* data,Cloud* cloud=nullptr);
   ~QtInputCloud();
   QWidget* validate();
+  void setCloud();
+  void getCloud();
 public slots:
   void emitRemove();
 signals:
@@ -75,7 +78,6 @@ signals:
 
 inline
 QtInputCloud::~QtInputCloud(){
-  source->removeCloud(cloud);
   delete double_validator;
   delete double_amplitude_validator;
 }

src/qt/input_clouds_tab.cpp

@@ -1,7 +1,7 @@
-#include "input_clouds_tab.hpp"
+#include "qt/input/clouds_tab.hpp"
 
-QtInputCloudsTab::QtInputCloudsTab(Source* _source):QWidget(){
-  source=_source;
+QtInputCloudsTab::QtInputCloudsTab(QtInputData* d):QWidget(){
+  data=d;
 
   main_layout=new QVBoxLayout;
   clouds_layout=new QVBoxLayout;
@@ -30,12 +30,18 @@ QtInputCloudsTab::QtInputCloudsTab(Source* _source):QWidget(){
 
 void
 QtInputCloudsTab::addCloud(){
-  QtInputCloud* input_cloud=new QtInputCloud(source);
+  QtInputCloud* input_cloud=new QtInputCloud(data);
   clouds_layout->addWidget(input_cloud);
   connect(input_cloud,&QtInputCloud::remove,this,&QtInputCloudsTab::removeCloud);
   emit sourcesChanged();
 }
 
+void
+QtInputCloudsTab::addExistingCloud(Cloud* cloud){
+  QtInputCloud* input_cloud=new QtInputCloud(data,cloud);
+  clouds_layout->addWidget(input_cloud);
+  connect(input_cloud,&QtInputCloud::remove,this,&QtInputCloudsTab::removeCloud);
+}
 void
 QtInputCloudsTab::removeCloud(QtInputCloud* input_cloud){
   disconnect(input_cloud,nullptr,nullptr,nullptr);
@@ -65,3 +71,10 @@ QtInputCloudsTab::updateSources(){
   }
   emit sourcesChanged();
 }
+
+void
+QtInputCloudsTab::getClouds(){
+  for(auto it=data->source.clouds.begin();it!=data->source.clouds.end();++it){
+   addExistingCloud(*it);
+  }
+}

src/qt/input_clouds_tab.hpp

@@ -12,14 +12,13 @@
 #include <QPalette>
 #include <QMessageBox>
 
-#include "qt/input_cloud.hpp"
-
-#include "source.hpp"
+#include "qt/input/cloud.hpp"
+#include "input_data.hpp"
 
 class QtInputCloudsTab:public QWidget{
   Q_OBJECT
 private:
-  Source* source;
+  QtInputData* data;
   QPushButton* add_button;
   QPushButton* refresh_button;
   QScrollArea* scroll_area;
@@ -29,11 +28,11 @@ private:
   QHBoxLayout* button_layout;
   QWidget* button_widget;
 public:
-  QtInputCloudsTab(Source* source);
+  QtInputCloudsTab(QtInputData* data);
   ~QtInputCloudsTab();
   QWidget* validate();
-  void save(fstream& file);
-  void load(fstream& file);
+  void getClouds();
+  void addExistingCloud(Cloud* cloud);
 public slots:
    void addCloud();
    void removeCloud(QtInputCloud*);

src/qt/input/data.cpp

@@ -0,0 +1,221 @@
+#include "qt/input/data.hpp"
+
+QtInputData::QtInputData(){
+  depth=1;
+  nZ_max=100;
+  initPhysics();
+  initSplines();
+  initGeometry();
+  initInitialState();
+}
+
+void QtInputData::initPhysics(){
+  Physics::g=9.81;
+  Physics::rho=1000;
+  Physics::phi=0.3;
+  Physics::k0=3e-5;
+  Physics::nivrivsat=0.02;
+  Physics::Psec=-10000;
+  Physics::model_data[0]=-2000;
+  Physics::model_data[1]=0;
+  Physics::model_data[2]=3;
+  Physics::model_data[3]=11;
+  Physics::model_data[4]=0;
+  Physics::model_data[5]=0;
+  Physics::setModel(Physics::BrooksCorey);
+}
+
+void
+QtInputData::initSplines() {
+  // Control points of splines
+  for(size_t i=0;i<np;++i){
+    double x=double(i)/(np-1);
+    point[i].x=x;
+    point[i].y=0.2;
+    point[np+i].x=x;
+    point[np+i].y=0.8;
+    point[2*np+i].x=x;
+    point[2*np+i].y=0.5;
+  }
+  // Control points of oveland points
+  for(size_t i=0;i<np_ov;++i){
+    double x=double(i+0.5)/np_ov;
+    point[3*np+i].x=x;
+    point[3*np+i].y=0.75;
+  }
+
+  // Splines
+  spline_bot.setPoints(point,np);
+  spline_bot.compute();
+  spline_soil.setPoints(&point[np],np);
+  spline_soil.compute();
+  spline_sat.setPoints(&point[2*np],np);
+  spline_sat.compute();
+}
+
+void
+QtInputData::initGeometry(){
+  geometry.lX=10;
+  geometry.nX=100;
+  geometry.dX=geometry.lX/(geometry.nX-1);
+  geometry.hsoil=new double[nmax_Qt];
+  geometry.dhsoil=new double[nmax_Qt];
+  geometry.hbot=new double[nmax_Qt];
+  geometry.dhbot=new double[nmax_Qt];
+  geometry.nZ=new size_t[nmax_Qt];
+  geometry.dZ=new double[nmax_Qt];
+  geometry.Z=new double*[nmax_Qt];
+  for(size_t i=0;i<nmax_Qt;++i){
+    geometry.Z[i]=new double[nmax_Qt];
+  }
+  updateGeometry();
+}
+
+void
+QtInputData::initInitialState(){
+  initial_state->hsat=new double[nmax_Qt];
+  initial_state->hov=new double[nmax_Qt];
+  initial_state->Pinit=new double*[nmax_Qt];
+  for(size_t i=0;i<nmax_Qt;++i){
+    initial_state->Pinit[i]=new double[nmax_Qt];
+  }
+  updateInitialState();
+}
+
+void
+QtInputData::updateGeometry(){
+  double prev_factor=factor;
+  double hs_max=-inf,hb_min=inf;
+  double v;
+  geometry.dX=geometry.lX/(geometry.nX-1);
+  for(size_t k=0;k<geometry.nX;++k){
+    double x=k*geometry.dX/geometry.lX;
+    v=geometry.hsoil[k]=spline_soil(x);
+    hs_max=max(v,hs_max);
+    geometry.dhsoil[k]=spline_soil.derivate(x);
+    v=geometry.hbot[k]=spline_bot(x);
+    hb_min=min(v,hb_min);
+    geometry.dhbot[k]=spline_bot.derivate(x);
+  }
+  double dZ_avg=depth/nZ_max;
+  factor=depth/(hs_max-hb_min);
+  for(size_t k=0;k<geometry.nX;++k){
+    geometry.hsoil[k]*=factor;
+    geometry.hbot[k]*=factor;
+  }
+  geometry.initZ(dZ_avg,false);
+  double r=factor/prev_factor;
+  for(auto it=initial_state->tanks.begin();it!=initial_state->tanks.end();++it){
+    Tank* tank=*it;
+    tank->top*=r;
+    tank->bottom*=r;
+    tank->delta_top*=r;
+    tank->delta_bottom*=r;
+  }
+  for(auto it=source.pumps.begin();it!=source.pumps.end();++it){
+    Pump* pump=*it;
+    pump->top_init*=r;
+    pump->bottom_init*=r;
+    pump->delta_top_init*=r;
+    pump->delta_bottom_init*=r;
+    pump->top_final*=r;
+    pump->bottom_final*=r;
+    pump->delta_top_final*=r;
+    pump->delta_bottom_final*=r;
+  }
+}
+
+void
+QtInputData::updateInitialState(){
+  for(size_t i=0;i<geometry.nX;++i){
+    double x=double(i)/(geometry.nX-1);
+    initial_state->hsat[i]=spline_sat(x)*factor;
+  }
+  initial_state->updatePressure();
+  updateOverland();
+}
+
+void
+QtInputData::save(fstream& file){
+  InputData::save(file);
+  file.write((char*)&depth,sizeof(double));
+  file.write((char*)&nZ_max,sizeof(size_t));
+  size_t temp=np;
+  file.write((char*)&temp,sizeof(size_t));
+  for(size_t i=0;i<3*np+np_ov;++i){
+    file.write((char*)&point[i].x,sizeof(double));
+    file.write((char*)&point[i].y,sizeof(double));
+  }
+}
+
+void
+QtInputData::load(fstream& file){
+  InputData::load(file);
+  file.read((char*)&depth,sizeof(double));
+  file.read((char*)&nZ_max,sizeof(size_t));
+  size_t temp;
+  file.read((char*)&temp,sizeof(size_t));
+  if(temp!=np){
+    cerr<<"Incorrect number of points"<<endl;
+    return;
+  }
+  for(size_t i=0;i<3*np+np_ov;++i){
+    file.read((char*)&point[i].x,sizeof(double));
+    file.read((char*)&point[i].y,sizeof(double));
+  }
+  updateSplineBot();
+  updateSplineSat();
+  updateSplineSoil();
+  updateGeometry();
+  updateInitialState();
+  updateOverland();
+}
+
+QtInputData::~QtInputData(){
+  delete[] geometry.hsoil;
+  delete[] geometry.dhsoil;
+  delete[] geometry.hbot;
+  delete[] geometry.dhbot;
+  delete[] geometry.nZ;
+  for(size_t i=0;i<nmax_Qt;++i){
+    delete[] geometry.Z[i];
+  }
+  delete[] geometry.Z;
+  geometry.hsoil=nullptr; //Bypass geometry destructor
+  delete[] initial_state->hsat;
+  for(size_t i=0;i<nmax_Qt;++i){
+    delete[] initial_state->Pinit[i];
+  }
+  delete[] initial_state->Pinit;
+  initial_state->hsat=nullptr; //Bypassfirst initial_state destructor
+}
+
+void
+QtInputData::updateOverland(){
+  for(size_t x=0;x<geometry.nX;++x){
+    initial_state->hov[x]=geometry.hsoil[x];
+  }
+  for(size_t i=0;i<np_ov;++i){
+    Point& P=point[3*np+i];
+    size_t x_ov=P.x*(geometry.nX-1);
+    double y=P.y*factor;
+    if(geometry.hsoil[x_ov]<y){
+      size_t x_left=findSoilOver(y,x_ov,-1);
+      size_t x_right=findSoilOver(y,x_ov,1);
+        for(size_t x=x_left;x<=x_right;++x){
+        initial_state->hov[x]=max(initial_state->hov[x],y);
+      }
+    }
+  }
+}
+
+size_t
+QtInputData::findSoilOver(double y,int x_ov,int x_step){
+  int x=x_ov;
+  while(x>=0 and x<geometry.nX and geometry.hsoil[x]<y){
+    x+=x_step;
+  }
+  x=min(x,(int)geometry.nX-1);
+  x=max(x,0);
+  return x;
+}

src/qt/input/data.hpp

@@ -0,0 +1,105 @@
+#ifndef QT_INPUT_DATA_HPP
+#define QT_INPUT_DATA_HPP
+
+#include "input_data.hpp"
+#include "math/spline.hpp"
+
+class QtInputData:public InputData{
+private:
+  size_t findSoilOver(double y,int x,int x_step);
+public:
+  static constexpr size_t const nmax_Qt=400;
+  static constexpr size_t const np=10;
+  static constexpr size_t const np_ov=5;
+  double depth;
+  size_t nZ_max;
+  Point point[3*np+np_ov];
+  Spline spline_soil;
+  Spline spline_bot;
+  Spline spline_sat;
+  QtInputData();
+  ~QtInputData();
+  void initSplines();
+  void initPhysics();
+  void initGeometry();
+  void initInitialState();
+  void updateSplineBot();
+  void updateSplineSat();
+  void updateSplineSoil();
+  void updateGeometry();
+  void updateInitialState();
+  void updateOverland();
+  double evalSplineBot(double x) const;
+  double evalSplineSat(double x) const;
+  double evalSplineSoil(double x) const;
+  Tank* addTank();
+  void removeTank(Tank* tank);
+  Pump* addPump();
+  void removePump(Pump* pump);
+  Cloud* addCloud();
+  void removeCloud(Cloud* cloud);
+  void save(fstream& file);
+  void load(fstream& file);
+};
+
+inline void
+QtInputData::updateSplineSoil(){
+  spline_soil.compute();
+}
+
+inline void
+QtInputData::updateSplineBot(){
+  spline_bot.compute();
+}
+
+inline void
+QtInputData::updateSplineSat(){
+  spline_sat.compute();
+}
+
+inline double
+QtInputData::evalSplineBot(double x) const{
+  return spline_bot(x);
+}
+
+inline double
+QtInputData::evalSplineSat(double x) const{
+  return spline_sat(x);
+}
+
+inline double
+QtInputData::evalSplineSoil(double x) const{
+  return spline_soil(x);
+}
+
+inline Tank*
+QtInputData::addTank(){
+  return initial_state->addTank();
+}
+
+inline void
+QtInputData::removeTank(Tank* tank){
+  initial_state->removeTank(tank);
+}
+
+inline Pump*
+QtInputData::addPump(){
+  return source.addPump();
+}
+
+inline void
+QtInputData::removePump(Pump* pump){
+  source.removePump(pump);
+}
+
+inline Cloud*
+QtInputData::addCloud(){
+  return source.addCloud();
+}
+
+inline void
+QtInputData::removeCloud(Cloud* cloud){
+  source.removeCloud(cloud);
+}
+
+#endif

src/qt/input_geometry.cpp

@@ -1,22 +1,24 @@
-#include "input_geometry.hpp"
+#include "qt/input/geometry.hpp"
 
-QtInputGeometry::QtInputGeometry(Geometry* _geometry):QWidget(){
-  geometry=_geometry;
+QtInputGeometry::QtInputGeometry(QtInputData* m):QWidget(){
+  data=m;
 
   main_layout=new QVBoxLayout;
   lX_label=new QLabel("Horizontal lenght of the domain (l<sub>X</sub>)");
   nX_label=new QLabel("Number of horizontal steps (n<sub>X</sub>)");
-  nZ_label=new QLabel("Maximal number of vertical steps (n<sub>Z</sub>)");
+  nZ_max_label=new QLabel("Maximal number of vertical steps (n<sub>Z</sub>)");
   depth_label=new QLabel("Maximal height between h<sub>soil</sub> and h<sub>bot</sub> (depth)");
-  lX_input=new QLineEdit(QString::number(geometry->lX));
-  nX_input=new QLineEdit(QString::number(geometry->nX));
-  nZ_input=new QLineEdit(QString::number(geometry->nZ_max));
-  depth_input=new QLineEdit(QString::number(geometry->depth));
 
+  lX_input=new QLineEdit();
+  nX_input=new QLineEdit();
+  nZ_max_input=new QLineEdit();
+  depth_input=new QLineEdit();
+
+//  loadGeometry();
   //Labels
   lX_label->setTextFormat(Qt::RichText);
   nX_label->setTextFormat(Qt::RichText);
-  nZ_label->setTextFormat(Qt::RichText);
+  nZ_max_label->setTextFormat(Qt::RichText);
   depth_label->setTextFormat(Qt::RichText);
 
   //Button
@@ -33,8 +35,8 @@ QtInputGeometry::QtInputGeometry(Geometry* _geometry):QWidget(){
   main_layout->addWidget(depth_label);
   main_layout->addWidget(depth_input);
   main_layout->addSpacing(vspace);
-  main_layout->addWidget(nZ_label);
-  main_layout->addWidget(nZ_input);
+  main_layout->addWidget(nZ_max_label);
+  main_layout->addWidget(nZ_max_input);
   main_layout->addSpacing(2*vspace);
   main_layout->addWidget(refresh_button);
   main_layout->addStretch();
@@ -46,13 +48,30 @@ QtInputGeometry::QtInputGeometry(Geometry* _geometry):QWidget(){
   positive_double_validator->setBottom(0);
   positive_int_validator=new QIntValidator;
   positive_int_validator->setBottom(1);
-  positive_int_validator->setTop(geometry->nmax_Qt);
+  positive_int_validator->setTop(QtInputData::nmax_Qt);
   lX_input->setValidator(positive_double_validator);
   nX_input->setValidator(positive_int_validator);
   depth_input->setValidator(positive_double_validator);
-  nZ_input->setValidator(positive_int_validator);
+  nZ_max_input->setValidator(positive_int_validator);
 
   connect(refresh_button,&QPushButton::clicked,this,&QtInputGeometry::emitGeometryChanged);
+
+  getGeometry();
+}
+
+void
+QtInputGeometry::emitGeometryChanged(){
+  QWidget* widget=validate();
+  if(widget!=nullptr){
+    QMessageBox msgBox;
+    msgBox.setText("Incorrect geometry entry");
+    msgBox.exec();
+    widget->setFocus();
+  }
+  else{
+    setGeometry();
+    emit geometryChanged();
+  }
 }
 
 QWidget*
@@ -60,34 +79,22 @@ QtInputGeometry::validate(){
   if(not lX_input->hasAcceptableInput()) return lX_input;
   if(not nX_input->hasAcceptableInput()) return nX_input;
   if(not depth_input->hasAcceptableInput()) return depth_input;
-  if(not nZ_input->hasAcceptableInput()) return nZ_input;
+  if(not nZ_max_input->hasAcceptableInput()) return nZ_max_input;
   return nullptr;
 }
 
 void
-QtInputGeometry::save(fstream& file){
-  double d;
-  size_t n;
-  d=lX_input->text().toDouble();
-  file.write((char*)&d,sizeof(double));
-  n=nX_input->text().toULong();
-  file.write((char*)&n,sizeof(size_t));
-  d=depth_input->text().toDouble();
-  file.write((char*)&d,sizeof(double));
-  n=nZ_input->text().toULong();
-  file.write((char*)&n,sizeof(size_t));
+QtInputGeometry::getGeometry(){
+  lX_input->setText(QString::number(data->geometry.lX));
+  nX_input->setText(QString::number(data->geometry.nX));
+  depth_input->setText(QString::number(data->depth));
+  nZ_max_input->setText(QString::number(data->nZ_max));
 }
 
 void
-QtInputGeometry::load(fstream& file){
-  double d;
-  size_t n;
-  file.read((char*)&d,sizeof(double));
-  lX_input->setText(QString::number(d));
-  file.read((char*)&n,sizeof(size_t));
-  nX_input->setText(QString::number(n));
-  file.read((char*)&d,sizeof(double));
-  depth_input->setText(QString::number(d));
-  file.read((char*)&n,sizeof(size_t));
-  nZ_input->setText(QString::number(n));
+QtInputGeometry::setGeometry(){
+  data->geometry.lX=lX_input->text().toDouble();
+  data->geometry.nX=nX_input->text().toULong();
+  data->depth=depth_input->text().toDouble();
+  data->nZ_max=nZ_max_input->text().toULong();
 }

src/qt/input_geometry.hpp

@@ -11,42 +11,35 @@
 #include <QDoubleValidator>
 #include <QIntValidator>
 #include <QPushButton>
-#include "geometry.hpp"
+#include <QMessageBox>
+#include "qt/input/data.hpp"
 
-static const double def_lX=10;
-static const size_t def_nX=100;
-static const size_t def_depth=1;
-static const size_t def_nZ=100;
 
 using namespace std;
 
 class QtInputGeometry:public QWidget{
   Q_OBJECT
 private:
-  Geometry* geometry;
+  QtInputData* data;
   QVBoxLayout* main_layout;
   QLabel* lX_label;
   QLabel* nX_label;
-  QLabel* nZ_label;
+  QLabel* nZ_max_label;
   QLabel* depth_label;
   QLineEdit* lX_input;
   QLineEdit* nX_input;
-  QLineEdit* nZ_input;
+  QLineEdit* nZ_max_input;
   QLineEdit* depth_input;
   QPushButton* refresh_button;
   QDoubleValidator* positive_double_validator;
   QIntValidator* positive_int_validator;
 
 public:
-  QtInputGeometry(Geometry* geometry);
+  QtInputGeometry(QtInputData* data);
   ~QtInputGeometry();
   QWidget* validate();
-  void save(fstream& file);
-  void load(fstream& file);
-  double get_lX();
-  size_t get_nX();
-  double get_depth();
-  size_t get_nZ_max();
+  void setGeometry();
+  void getGeometry();
 public slots:
   void emitGeometryChanged();
 signals:
@@ -58,29 +51,4 @@ QtInputGeometry::~QtInputGeometry(){
   delete positive_double_validator;
   delete positive_int_validator;
 }
-
-inline void
-QtInputGeometry::emitGeometryChanged(){
-  emit geometryChanged();
-}
-
-inline double
-QtInputGeometry::get_lX(){
-  return lX_input->text().toDouble();
-}
-
-inline size_t
-QtInputGeometry::get_nX(){
-  return nX_input->text().toULong();
-}
-
-inline double
-QtInputGeometry::get_depth(){
-  return depth_input->text().toDouble();
-}
-
-inline size_t
-QtInputGeometry::get_nZ_max(){
-  return nZ_input->text().toULong();
-}
 #endif

src/qt/input_initial_state.cpp

@@ -1,7 +1,7 @@
-#include "input_initial_state.hpp"
+#include "qt/input/initial_state.hpp"
 
-QtInputInitialState::QtInputInitialState(InitialState* _initial_state):QWidget(){
-  initial_state=_initial_state;
+QtInputInitialState::QtInputInitialState(QtInputData* d):QWidget(){
+  data=d;
 
   main_layout=new QVBoxLayout;
   tanks_layout=new QVBoxLayout;
@@ -30,12 +30,19 @@ QtInputInitialState::QtInputInitialState(InitialState* _initial_state):QWidget()
 
 void
 QtInputInitialState::addTank(){
-  QtInputTank* input_tank=new QtInputTank(initial_state);
+  QtInputTank* input_tank=new QtInputTank(data);
   tanks_layout->addWidget(input_tank);
   connect(input_tank,&QtInputTank::remove,this,&QtInputInitialState::removeTank);
   emit initialStateChanged();
 }
 
+void
+QtInputInitialState::addExistingTank(Tank* tank){
+  QtInputTank* input_tank=new QtInputTank(data,tank);
+  tanks_layout->addWidget(input_tank);
+  connect(input_tank,&QtInputTank::remove,this,&QtInputInitialState::removeTank);
+}
+
 void
 QtInputInitialState::removeTank(QtInputTank* input_tank){
   disconnect(input_tank,nullptr,nullptr,nullptr);
@@ -72,3 +79,12 @@ QtInputInitialState::save(fstream& file){
   //d=lX_input->text().toDouble();
   //file.write((char*)&d,sizeof(double));
 }
+
+
+void
+QtInputInitialState::getTanks(){
+  for(auto it=data->initial_state->tanks.begin();it!=data->initial_state->tanks.end();++it){
+   addExistingTank(*it);
+
+  }
+}

src/qt/input_initial_state.hpp

@@ -12,14 +12,13 @@
 #include <QPalette>
 #include <QMessageBox>
 
-#include "qt/input_tank.hpp"
-
-#include "initial_state.hpp"
+#include "qt/input/tank.hpp"
+#include "qt/input/data.hpp"
 
 class QtInputInitialState:public QWidget{
   Q_OBJECT
 private:
-  InitialState* initial_state;
+  QtInputData* data;
   QPushButton* add_button;
   QPushButton* refresh_button;
   QScrollArea* scroll_area;
@@ -29,15 +28,18 @@ private:
   QHBoxLayout* button_layout;
   QWidget* button_widget;
 public:
-  QtInputInitialState(InitialState* initial_state);
+  QtInputInitialState(QtInputData* data);
   ~QtInputInitialState();
   QWidget* validate();
   void save(fstream& file);
   void load(fstream& file);
+  void getTanks();
+  void addExistingTank(Tank* tank);
 public slots:
-   void addTank();
-   void removeTank(QtInputTank*);
-   void updateInitialState();
+  void addTank();
+  void removeTank(QtInputTank*);
+  void updateInitialState();
+
 signals:
    void initialStateChanged();
 };

src/qt/input.cpp

@@ -1,10 +1,8 @@
-#include "input.hpp"
+#include "qt/input/input.hpp"
 
 QtInput::QtInput():QWidget(){
-  geometry=new Geometry;
-  geometry->initDefault();
-  initial_state=new InitialState;
-  source=new Source;
+
+  data=new QtInputData;
 
   main_layout=new QVBoxLayout;
   button_layout=new QHBoxLayout;
@@ -12,19 +10,15 @@ QtInput::QtInput():QWidget(){
   button_widget=new QWidget;
   button_save=new QPushButton("Save");
   button_cancel=new QPushButton("Cancel");
-  button_run=new QPushButton("Run");
+  time_bar=new QScrollBar(Qt::Horizontal);
 
   input_physics=new QtInputPhysics;
   input_time=new QtInputTime;
-  input_geometry=new QtInputGeometry(geometry);
-  input_initial_state=new QtInputInitialState(initial_state);
-  input_pump_tab=new QtInputPumpTab(source);
-  input_clouds_tab=new QtInputCloudsTab(source);
-  input_view=new QtInputView(input_geometry,geometry,initial_state,source);
-
-  time_bar=new QScrollBar(Qt::Horizontal);
-//  input_view->setGeometry(geometry);
-//  input_view->updateGeometry();
+  input_geometry=new QtInputGeometry(data);
+  input_initial_state=new QtInputInitialState(data);
+  input_pump_tab=new QtInputPumpTab(data);
+  input_clouds_tab=new QtInputCloudsTab(data);
+  input_view=new QtInputView(data);
 
   //Tab
   tab_widget->addTab(input_physics,"Physics");
@@ -35,7 +29,6 @@ QtInput::QtInput():QWidget(){
   tab_widget->addTab(input_clouds_tab,"Clouds");
 
   //Buttons
-  button_layout->addWidget(button_run);
   button_layout->addWidget(button_save);
   button_layout->addWidget(button_cancel);
   button_widget->setLayout(button_layout);
@@ -50,22 +43,22 @@ QtInput::QtInput():QWidget(){
   //Conectors
   connect(button_save,&QPushButton::clicked,this,&QtInput::save);
   connect(button_cancel,&QPushButton::clicked,this,&QtInput::cancel);
-  connect(button_run,&QPushButton::clicked,this,&QtInput::run);
   connect(input_geometry,&QtInputGeometry::geometryChanged,this,&QtInput::updateGeometry);
+  connect(input_physics,&QtInputPhysics::physicsChanged,this,&QtInput::updateInitialState);
+  connect(input_view,&QtInputView::geometryChanged,this,&QtInput::updateGeometry);
   connect(tab_widget,&QTabWidget::currentChanged,this,&QtInput::changeTabIndex);
   connect(input_initial_state,&QtInputInitialState::initialStateChanged,this,&QtInput::updateInitialState);
   connect(time_bar,&QScrollBar::valueChanged,input_view,&QtInputView::setTime);
   connect(input_pump_tab,&QtInputPumpTab::sourcesChanged,this,&QtInput::updateSource);
   connect(input_clouds_tab,&QtInputCloudsTab::sourcesChanged,this,&QtInput::updateSource);
+  connect(input_view,&QtInputView::timeChanged,this,&QtInput::updateSource);
   previous_index=-1;
-
 }
 
 QtInput::QtInput(QString filename):QtInput(){
   load(filename.toStdString());
 }
 
-
 bool
 QtInput::validate(){
   QWidget* widget=input_physics->validate();
@@ -120,23 +113,15 @@ QtInput::save(){
   }
 }
 
-void
-QtInput::run(){
-  if(validate()){
-    emit run_signal();
-  }
-}
 
 void
 QtInput::save_input(string filename){
+  updateGeometry();
   fstream file;
   file.open(filename.c_str(),fstream::out|fstream::trunc|fstream::binary);
-  input_physics->save(file);
-  input_time->save(file);
-  input_geometry->save(file);
-  input_initial_state->save(file);
-  //input_pump_tab->save(file);
-  //input_clouds_tab->save(file);
+  input_physics->setPhysics();
+  input_time->setTime();
+  data->save(file);
   file.close();
 }
 
@@ -149,9 +134,13 @@ void
 QtInput::load(string filename){
   fstream file;
   file.open(filename.c_str(),fstream::in|fstream::binary);
-  input_physics->load(file);
-  input_time->load(file);
-  input_geometry->load(file);
+  data->load(file);
+  input_physics->getPhysics();
+  input_time->getTime();
+  input_geometry->getGeometry();
+  input_initial_state->getTanks();
+  input_pump_tab->getPumps();
+  input_clouds_tab->getClouds();
   file.close();
 }
 
@@ -181,35 +170,29 @@ QtInput::changeTabIndex(int index){
     break;
     default:
     break;
-
   }
   previous_index=index;
   input_view->update();
-
 }
 
 void
 QtInput::updateGeometry(){
-  QWidget* widget=input_geometry->validate();
-  if(widget!=nullptr){
-    QMessageBox msgBox;
-    msgBox.setText("Incorrect geometry entry");
-    msgBox.exec();
-    tab_widget->setCurrentWidget(input_geometry);
-    widget->setFocus();
-  }
-  else{
-    input_view->updateGeometry();
-  }
+  data->updateGeometry();
+  data->updateInitialState();
+  input_view->update();
 }
 
 void
 QtInput::updateInitialState(){
-  initial_state->updatePressure();
-  input_view->updateInitialState();
+  data->updateInitialState();
+  input_view->update();
 }
 
 void
 QtInput::updateSource(){
-  input_view->updateSource();
+  input_view->update();
+}
+
+QtInput::~QtInput(){
+  delete data;
 }

src/qt/input.hpp

@@ -10,14 +10,14 @@
 #include <QScrollBar>
 #include <fstream>
 
-#include "input_physics.hpp"
-#include "input_time.hpp"
-#include "input_geometry.hpp"
-#include "input_initial_state.hpp"
-#include "input_pump_tab.hpp"
-#include "input_clouds_tab.hpp"
-#include "input_view.hpp"
-#include "geometry.hpp"
+#include "qt/input/physics.hpp"
+#include "qt/input/time.hpp"
+#include "qt/input/geometry.hpp"
+#include "qt/input/initial_state.hpp"
+#include "qt/input/pump_tab.hpp"
+#include "qt/input/clouds_tab.hpp"
+#include "qt/input/view.hpp"
+#include "qt/input/data.hpp"
 
 using namespace std;
 
@@ -25,16 +25,14 @@ class QtInput:public QWidget{
   Q_OBJECT
 private:
   int previous_index;
-  Geometry* geometry;
-  InitialState* initial_state;
-  Source* source;
+  QtInputData* data;
   QVBoxLayout* main_layout;
   QHBoxLayout* button_layout;
   QTabWidget* tab_widget;
   QWidget* button_widget;
   QPushButton* button_save;
-  QPushButton* button_run;
   QPushButton* button_cancel;
+  QScrollBar* time_bar;
   QtInputPhysics* input_physics;
   QtInputTime* input_time;
   QtInputGeometry* input_geometry;
@@ -42,8 +40,7 @@ private:
   QtInputPumpTab* input_pump_tab;
   QtInputCloudsTab* input_clouds_tab;
   QtInputView* input_view;
-  QScrollBar* time_bar;
-
+  void initDefaultGeometry();
   bool validate();
   void load(string filename);
 public:
@@ -51,12 +48,11 @@ public:
   QtInput(QString filename);
   ~QtInput();
   void save_input(string filename);
-
+  const QtInputData* getData() const;
 signals:
   void run_signal();
   void exit_signal();
 private slots:
-  void run();
   void save();
   void cancel();
   void changeTabIndex(int index);
@@ -65,10 +61,9 @@ private slots:
   void updateSource();
 };
 
-inline QtInput::~QtInput(){
-  delete geometry;
-  delete initial_state;
-  delete source;
-}
 
+inline const QtInputData*
+QtInput::getData() const{
+  return data;
+}
 #endif

src/qt/mainwindow.cpp

@@ -1,10 +1,7 @@
-#include "mainwindow.hpp"
+#include "qt/input/main_window.hpp"
 
 QtMainWindow::QtMainWindow():QMainWindow(){
   input=nullptr;
-  //view_solution=nullptr;
-
-  //Actions
   new_act=new QAction("New input",this);
   load_act=new QAction("Load input",this);
   exit_act=new QAction("Exit",this);
@@ -20,32 +17,21 @@ QtMainWindow::QtMainWindow():QMainWindow(){
   connect(exit_act,&QAction::triggered,this,&QtMainWindow::exit);
 }
 
-QtMainWindow::QtMainWindow(string filename):QtMainWindow(){
-  input=new QtInput(QString::fromStdString(filename));
-  run_input();
-
-}
-
 void
 QtMainWindow::new_input(){
   input=new QtInput;
   setCentralWidget(input);
-  connect(input,&QtInput::run_signal,this,&QtMainWindow::run_input);
   connect(input,&QtInput::exit_signal,this,&QtMainWindow::exit_input);
-  new_act->setEnabled(false);
-  load_act->setEnabled(false);
 }
 
 void
-QtMainWindow::run_input(){
-  /*disconnect(input,nullptr,nullptr,nullptr);
+QtMainWindow::load_input(){
+  QString filename=QFileDialog::getOpenFileName(this,"Load input","inputs/","QT input file (*.input)");
+  input=new QtInput(filename);
+  setCentralWidget(input);
+  connect(input,&QtInput::exit_signal,this,&QtMainWindow::exit_input);
   new_act->setEnabled(false);
   load_act->setEnabled(false);
-  kernel=new Kernel(input->getPhysics(),input->getTime(),input->getGeometry());
-  delete input;
-  input=nullptr;
-  view_solution=new QtViewSolution(kernel);
-  setCentralWidget(view_solution);*/
 }
 
 void
@@ -58,18 +44,6 @@ QtMainWindow::exit_input(){
   load_act->setEnabled(true);
 }
 
-void
-QtMainWindow::load_input(){
-  QString filename=QFileDialog::getOpenFileName(this,"Load input","inputs/","QT input file (*.input)");
-  input=new QtInput(filename);
-  //input=new QtInput;
-  setCentralWidget(input);
-  connect(input,&QtInput::run_signal,this,&QtMainWindow::run_input);
-  connect(input,&QtInput::exit_signal,this,&QtMainWindow::exit_input);
-  new_act->setEnabled(false);
-  load_act->setEnabled(false);
-}
-
 void
 QtMainWindow::exit(){
   QApplication::quit();

src/qt/mainwindow.hpp

@@ -7,8 +7,8 @@
 #include <QAction>
 #include <QApplication>
 
-#include "qt/input.hpp"
-#include "kernel.hpp"
+#include "qt/input/input.hpp"
+
 
 class QtMainWindow:public QMainWindow{
   Q_OBJECT
@@ -17,17 +17,14 @@ private:
   QMenu* input_menu;
   QAction* new_act;
   QAction* load_act;
-  QAction* exit_act;
-  Kernel* kernel;
+    QAction* exit_act;
 public:
   QtMainWindow();
-  QtMainWindow(string filename);
   ~QtMainWindow();
 private slots:
   void new_input();
   void load_input();
   void exit();
-  void run_input();
   void exit_input();
 };
 

src/qt/input_physics.cpp

@@ -1,4 +1,4 @@
-#include "input_physics.hpp"
+#include "qt/input/physics.hpp"
 
 QtInputPhysics::QtInputPhysics():QWidget(){
   g_label=new QLabel("Gravity (g) : ");
@@ -11,21 +11,22 @@ QtInputPhysics::QtInputPhysics():QWidget(){
   phi_input=new QLineEdit();
   k0_input=new QLineEdit();
   nivrivsat_input=new QLineEdit();
-  for(size_t i=0;i<max_model_parameters;++i){
+  for(size_t i=0;i<Physics::max_model_parameters;++i){
     model_label[i]=new QLabel(this);
     model_input[i]=new QLineEdit(this);
   }
 
+
   //Text format
   rho_label->setTextFormat(Qt::RichText);
   phi_label->setTextFormat(Qt::RichText);
   k0_label->setTextFormat(Qt::RichText);
-  for(size_t i=0;i<max_model_parameters;++i){
+  for(size_t i=0;i<Physics::max_model_parameters;++i){
     model_label[i]->setTextFormat(Qt::RichText);
   }
 
   //Hide model parameters
-  for(size_t i=0;i<max_model_parameters;++i){
+  for(size_t i=0;i<Physics::max_model_parameters;++i){
     model_label[i]->hide();
     model_input[i]->hide();
   }
@@ -34,7 +35,6 @@ QtInputPhysics::QtInputPhysics():QWidget(){
   positive_double_validator=new QDoubleValidator;
   positive_double_validator->setBottom(0);
 
-
   g_input->setValidator(positive_double_validator);
   rho_input->setValidator(double_validator);
   phi_input->setValidator(double_validator);
@@ -47,9 +47,16 @@ QtInputPhysics::QtInputPhysics():QWidget(){
   phy_model_selection_box=new QComboBox;
   phy_model_selection_box->addItem("Brooks and Corey");
 
+  //Button
+  refresh_button=new QPushButton("Refresh");
+
+
+  //Widget
+  phy_base_widget=new QWidget;
 
   //Layouts
-  base_layout=new QHBoxLayout;
+  main_layout=new QVBoxLayout;
+  phy_base_layout=new QHBoxLayout;
   phy_common_layout=new QVBoxLayout;
   phy_model_layout=new QVBoxLayout;
 
@@ -73,7 +80,7 @@ QtInputPhysics::QtInputPhysics():QWidget(){
 
   //Phy model layout
   phy_model_layout->addWidget(phy_model_selection_box);
-  for(size_t i=0;i<max_model_parameters;++i){
+  for(size_t i=0;i<Physics::max_model_parameters;++i){
     phy_model_layout->addSpacing(vspace);
     phy_model_layout->addWidget(model_label[i]);
     phy_model_layout->addWidget(model_input[i]);
@@ -85,9 +92,14 @@ QtInputPhysics::QtInputPhysics():QWidget(){
   phy_model_box->setLayout(phy_model_layout);
 
   //Base widget
-  base_layout->addWidget(phy_common_box);
-  base_layout->addWidget(phy_model_box);
-  setLayout(base_layout);
+  phy_base_layout->addWidget(phy_common_box);
+  phy_base_layout->addWidget(phy_model_box);
+  phy_base_widget->setLayout(phy_base_layout);
+
+  //Main_lyaout
+  main_layout->addWidget(phy_base_widget);
+  main_layout->addWidget(refresh_button);
+  setLayout(main_layout);
 
   //Conections
   connect(phy_model_selection_box,QOverload<int>::of(&QComboBox::activated),this,&QtInputPhysics::modelChoosed);
@@ -96,14 +108,14 @@ QtInputPhysics::QtInputPhysics():QWidget(){
   nb_model_parameters=0;
   modelChoosed(phy_model_selection_box->currentIndex());
 
-  loadData();
-}
-
+  connect(refresh_button,&QPushButton::clicked,this,&QtInputPhysics::emitPhysicsChanged);
 
+  getPhysics();
+}
 
 void
 QtInputPhysics::modelChoosed(int index){
-  for(size_t i=0;i<nb_model_parameters;++i){
+  for(size_t i=0;i<Physics::max_model_parameters;++i){
     model_label[i]->hide();
     model_input[i]->hide();
   }
@@ -148,57 +160,42 @@ QtInputPhysics::validate(){
 }
 
 void
-QtInputPhysics::save(fstream& file){
-  double d;
-  size_t s;
-  d=g_input->text().toDouble();
-  file.write((char*)&d,sizeof(double));
-  d=rho_input->text().toDouble();
-  file.write((char*)&d,sizeof(double));
-  d=phi_input->text().toDouble();
-  file.write((char*)&d,sizeof(double));
-  d=k0_input->text().toDouble();
-  file.write((char*)&d,sizeof(double));
-  d=nivrivsat_input->text().toDouble();
-  file.write((char*)&d,sizeof(double));
-  s=phy_model_selection_box->currentIndex();
-  file.write((char*)&s,sizeof(size_t));
-  s=nb_model_parameters;
-  file.write((char*)&s,sizeof(size_t));
-  for(size_t i=0;i<nb_model_parameters;++i){
-    d=model_input[i]->text().toDouble();
-    file.write((char*)&d,sizeof(double));
+QtInputPhysics::getPhysics(){
+  g_input->setText(QString::number(Physics::g));
+  rho_input->setText(QString::number(Physics::rho));
+  phi_input->setText(QString::number(Physics::phi));
+  k0_input->setText(QString::number(Physics::k0));
+  nivrivsat_input->setText(QString::number(Physics::nivrivsat));
+  phy_model_selection_box->setCurrentIndex(Physics::model);
+  for(size_t i=0;i<Physics::max_model_parameters;++i){
+    model_input[i]->setText(QString::number(Physics::model_data[i]));
   }
 }
 
 void
-QtInputPhysics::load(fstream& file){
-  double d;
-  size_t s;
-  file.read((char*)&d,sizeof(double));
-  g_input->setText(QString::number(d));
-  file.read((char*)&d,sizeof(double));
-  rho_input->setText(QString::number(d));
-  file.read((char*)&d,sizeof(double));
-  phi_input->setText(QString::number(d));
-  file.read((char*)&d,sizeof(double));
-  k0_input->setText(QString::number(d));
-  file.read((char*)&d,sizeof(double));
-  nivrivsat_input->setText(QString::number(d));
-  file.read((char*)&s,sizeof(size_t));
-  phy_model_selection_box->setCurrentIndex(s);
-  file.read((char*)&s,sizeof(size_t));
-  for(size_t i=0;i<s;++i){
-    file.read((char*)&d,sizeof(double));
-    model_input[i]->setText(QString::number(d));
+QtInputPhysics::setPhysics(){
+  Physics::g=g_input->text().toDouble();
+  Physics::rho=rho_input->text().toDouble();
+  Physics::phi=phi_input->text().toDouble();
+  Physics::k0=k0_input->text().toDouble();
+  Physics::nivrivsat=nivrivsat_input->text().toDouble();
+  //Model must be already assigned by QtInputPhysics::modelChoosed
+  for(size_t i=0;i<Physics::max_model_parameters;++i){
+    Physics::model_data[i]=model_input[i]->text().toDouble();
   }
 }
 
 void
-QtInputPhysics::loadData(){
-  g_input->setText(QString::number(Physics::g));
-  rho_input->setText(QString::number(Physics::rho));
-  phi_input->setText(QString::number(Physics::phi));
-  k0_input->setText(QString::number(Physics::k0));
-  nivrivsat_input->setText(QString::number(Physics::nivrivsat));
+QtInputPhysics::emitPhysicsChanged(){
+  QWidget* widget=validate();
+  if(widget!=nullptr){
+    QMessageBox msgBox;
+    msgBox.setText("Incorrect geometry entry");
+    msgBox.exec();
+    widget->setFocus();
+  }
+  else{
+    setPhysics();
+    emit physicsChanged();
+  }
 }

src/qt/input_physics.hpp

@@ -1,7 +1,7 @@
 #ifndef QT_INPUT_PHYSICS_HPP
 #define QT_INPUT_PHYSICS_HPP
 
-#include <iostream>
+
 #include <QLabel>
 #include <QWidget>
 #include <QGroupBox>
@@ -10,13 +10,16 @@
 #include <QVBoxLayout>
 #include <QLineEdit>
 #include <QDoubleValidator>
+#include <QPushButton>
+#include <QMessageBox>
+
+#include <iostream>
 #include <fstream>
 
-#include "physics.hpp"
+#include "../../physics.hpp"
 
 using namespace std;
 
-static const size_t max_model_parameters=6;
 
 class QtInputPhysics:public QWidget{
   Q_OBJECT
@@ -24,22 +27,25 @@ private:
   QGroupBox* phy_common_box;
   QGroupBox* phy_model_box;
   QComboBox* phy_model_selection_box;
+  QWidget* phy_base_widget;
   QLabel* g_label;
   QLabel* rho_label;
   QLabel* phi_label;
   QLabel* k0_label;
   QLabel* nivrivsat_label;
-  QLabel* model_label[max_model_parameters];
+  QLabel* model_label[Physics::max_model_parameters];
   QLineEdit* g_input;
   QLineEdit* rho_input;
   QLineEdit* phi_input;
   QLineEdit* k0_input;
   QLineEdit* nivrivsat_input;
-  QLineEdit* model_input[max_model_parameters];
-  QHBoxLayout* base_layout;
+  QLineEdit* model_input[Physics::max_model_parameters];
+  QHBoxLayout* phy_base_layout;
   QVBoxLayout* phy_common_layout;
   QVBoxLayout* phy_model_layout;
+  QVBoxLayout* main_layout;
   QSpacerItem* phy_model_layout_spacer;
+  QPushButton* refresh_button;
   QDoubleValidator* double_validator;
   QDoubleValidator* positive_double_validator;
   size_t nb_model_parameters;
@@ -47,15 +53,18 @@ public:
   QtInputPhysics();
   ~QtInputPhysics();
   QWidget* validate();
-  void save(fstream& file);
-  void load(fstream& file);
-  void loadData();
+  void setPhysics();
+  void getPhysics();
 private slots:
   void modelChoosed(int index);
+  void emitPhysicsChanged();
+signals:
+  void physicsChanged();
 };
 
 inline QtInputPhysics::~QtInputPhysics(){
   delete double_validator;
   delete positive_double_validator;
 }
+
 #endif

src/qt/input_pump.cpp

@@ -1,9 +1,9 @@
-#include "input_pump.hpp"
+#include "qt/input/pump.hpp"
 
-QtInputPump::QtInputPump(Source* _source){
-  source=_source;
-  pump=source->addPump();
+QtInputPump::QtInputPump(QtInputData* d,Pump* p){
+  data=d;
 
+  pump=(p==nullptr)?data->addPump():p;
 
   init_groupbox=new QGroupBox("Initial values");
   amplitude_init_widget=new QWidget;
@@ -30,19 +30,19 @@ QtInputPump::QtInputPump(Source* _source){
 
 
   amplitude_init_label=new QLabel("Amplitude: ");
-  amplitude_init_input=new QLineEdit(QString::number(pump->amplitude_init));
+  amplitude_init_input=new QLineEdit();
   amplitude_init_layout->addWidget(amplitude_init_label);
   amplitude_init_layout->addWidget(amplitude_init_input);
   amplitude_init_widget->setLayout(amplitude_init_layout);
 
   left_init_label=new QLabel("Left: ");
   delta_left_init_label=new QLabel("  Left delta: ");
-  left_init_input=new QLineEdit(QString::number(pump->left_init));
-  delta_left_init_input=new QLineEdit(QString::number(pump->delta_left_init));
+  left_init_input=new QLineEdit();
+  delta_left_init_input=new QLineEdit();
   right_init_label=new QLabel("  Right: ");
   delta_right_init_label=new QLabel("  Right delta: ");
-  right_init_input=new QLineEdit(QString::number(pump->right_init));
-  delta_right_init_input=new QLineEdit(QString::number(pump->delta_right_init));
+  right_init_input=new QLineEdit();
+  delta_right_init_input=new QLineEdit();
 
   left_right_init_layout->addWidget(left_init_label);
   left_right_init_layout->addWidget(left_init_input,1);
@@ -56,12 +56,12 @@ QtInputPump::QtInputPump(Source* _source){
 
   bottom_init_label=new QLabel("Bottom: ");
   delta_bottom_init_label=new QLabel("  Bottom delta: ");
-  bottom_init_input=new QLineEdit(QString::number(pump->bottom_init));
-  delta_bottom_init_input=new QLineEdit(QString::number(pump->delta_bottom_init));
+  bottom_init_input=new QLineEdit();
+  delta_bottom_init_input=new QLineEdit();
   top_init_label=new QLabel("  Top: ");
   delta_top_init_label=new QLabel("  Top delta: ");
-  top_init_input=new QLineEdit(QString::number(pump->top_init));
-  delta_top_init_input=new QLineEdit(QString::number(pump->delta_top_init));
+  top_init_input=new QLineEdit();
+  delta_top_init_input=new QLineEdit();
 
   bottom_top_init_layout->addWidget(bottom_init_label);
   bottom_top_init_layout->addWidget(bottom_init_input,1);
@@ -80,19 +80,19 @@ QtInputPump::QtInputPump(Source* _source){
 
 
   amplitude_final_label=new QLabel("Amplitude: ");
-  amplitude_final_input=new QLineEdit(QString::number(pump->amplitude_final));
+  amplitude_final_input=new QLineEdit();
   amplitude_final_layout->addWidget(amplitude_final_label);
   amplitude_final_layout->addWidget(amplitude_final_input);
   amplitude_final_widget->setLayout(amplitude_final_layout);
 
   left_final_label=new QLabel("Left: ");
   delta_left_final_label=new QLabel("  Left delta: ");
-  left_final_input=new QLineEdit(QString::number(pump->left_final));
-  delta_left_final_input=new QLineEdit(QString::number(pump->delta_left_final));
+  left_final_input=new QLineEdit();
+  delta_left_final_input=new QLineEdit();
   right_final_label=new QLabel("  Right: ");
   delta_right_final_label=new QLabel("  Right delta: ");
-  right_final_input=new QLineEdit(QString::number(pump->right_final));
-  delta_right_final_input=new QLineEdit(QString::number(pump->delta_right_final));
+  right_final_input=new QLineEdit();
+  delta_right_final_input=new QLineEdit();
 
   left_right_final_layout->addWidget(left_final_label);
   left_right_final_layout->addWidget(left_final_input,1);
@@ -106,12 +106,12 @@ QtInputPump::QtInputPump(Source* _source){
 
   bottom_final_label=new QLabel("Bottom: ");
   delta_bottom_final_label=new QLabel("  Bottom delta: ");
-  bottom_final_input=new QLineEdit(QString::number(pump->bottom_final));
-  delta_bottom_final_input=new QLineEdit(QString::number(pump->delta_bottom_final));
+  bottom_final_input=new QLineEdit();
+  delta_bottom_final_input=new QLineEdit();
   top_final_label=new QLabel("  Top: ");
   delta_top_final_label=new QLabel("  Top delta: ");
-  top_final_input=new QLineEdit(QString::number(pump->top_final));
-  delta_top_final_input=new QLineEdit(QString::number(pump->delta_top_final));
+  top_final_input=new QLineEdit();
+  delta_top_final_input=new QLineEdit();
 
   bottom_top_final_layout->addWidget(bottom_final_label);
   bottom_top_final_layout->addWidget(bottom_final_input,1);
@@ -167,6 +167,18 @@ QtInputPump::QtInputPump(Source* _source){
   setFrameShape(QFrame::Box);
 
   connect(remove_button,&QPushButton::clicked,this,&QtInputPump::emitRemove);
+
+  if(p==nullptr){
+    pump->bottom_init*=data->factor;
+    pump->top_init*=data->factor;
+    pump->delta_bottom_init*=data->factor;
+    pump->delta_top_init*=data->factor;
+    pump->bottom_final*=data->factor;
+    pump->top_final*=data->factor;
+    pump->delta_bottom_final*=data->factor;
+    pump->delta_top_final*=data->factor;
+  }
+  getPump();
 }
 
 
@@ -204,23 +216,60 @@ QtInputPump::validate(){
   if(bi>=ti) return top_init_input;
   if(lf>=rf) return right_final_input;
   if(bf>=tf) return top_final_input;
+  setPump();
+  return nullptr;
+}
+
+void
+QtInputPump::setPump(){
+  double li,lf,ri,rf,ti,tf,bi,bf;
+  li=left_init_input->text().toDouble();
+  lf=left_final_input->text().toDouble();
+  ri=right_init_input->text().toDouble();
+  rf=right_final_input->text().toDouble();
+  ti=top_init_input->text().toDouble();
+  tf=top_final_input->text().toDouble();
+  bi=bottom_init_input->text().toDouble();
+  bf=bottom_final_input->text().toDouble();
   pump->left_init=li;
   pump->right_init=ri;
-  pump->top_init=ti;
-  pump->bottom_init=bi;
+  pump->top_init=ti*data->factor;
+  pump->bottom_init=bi*data->factor;
   pump->left_final=lf;
   pump->right_final=rf;
-  pump->top_final=tf;
-  pump->bottom_final=bf;
+  pump->top_final=tf*data->factor;
+  pump->bottom_final=bf*data->factor;
   pump->amplitude_init=amplitude_init_input->text().toDouble();
   pump->amplitude_final=amplitude_final_input->text().toDouble();
   pump->delta_left_init=delta_left_init_input->text().toDouble();
   pump->delta_right_init=delta_right_init_input->text().toDouble();
-  pump->delta_top_init=delta_top_init_input->text().toDouble();
-  pump->delta_bottom_init=delta_bottom_init_input->text().toDouble();
+  pump->delta_top_init=delta_top_init_input->text().toDouble()*data->factor;
+  pump->delta_bottom_init=delta_bottom_init_input->text().toDouble()*data->factor;
   pump->delta_left_final=delta_left_final_input->text().toDouble();
   pump->delta_right_final=delta_right_final_input->text().toDouble();
-  pump->delta_top_final=delta_top_final_input->text().toDouble();
-  pump->delta_bottom_final=delta_bottom_final_input->text().toDouble();
-  return nullptr;
+  pump->delta_top_final=delta_top_final_input->text().toDouble()*data->factor;
+  pump->delta_bottom_final=delta_bottom_final_input->text().toDouble()*data->factor;
+}
+
+void
+QtInputPump::getPump(){
+  double factor_inv=1/data->factor;
+  amplitude_init_input->setText(QString::number(pump->amplitude_init));
+  left_init_input->setText(QString::number(pump->left_init));
+  right_init_input->setText(QString::number(pump->right_init));
+  top_init_input->setText(QString::number(pump->top_init*factor_inv));
+  bottom_init_input->setText(QString::number(pump->bottom_init*factor_inv));
+  delta_left_init_input->setText(QString::number(pump->delta_left_init));
+  delta_right_init_input->setText(QString::number(pump->delta_right_init));
+  delta_top_init_input->setText(QString::number(pump->delta_top_init*factor_inv));
+  delta_bottom_init_input->setText(QString::number(pump->delta_bottom_init*factor_inv));
+  amplitude_final_input->setText(QString::number(pump->amplitude_final));
+  left_final_input->setText(QString::number(pump->left_final));
+  right_final_input->setText(QString::number(pump->right_final));
+  top_final_input->setText(QString::number(pump->top_final*factor_inv));
+  bottom_final_input->setText(QString::number(pump->bottom_final*factor_inv));
+  delta_left_final_input->setText(QString::number(pump->delta_left_final));
+  delta_right_final_input->setText(QString::number(pump->delta_right_final));
+  delta_top_final_input->setText(QString::number(pump->delta_top_final*factor_inv));
+  delta_bottom_final_input->setText(QString::number(pump->delta_bottom_final*factor_inv));
 }

src/qt/input_pump.hpp

@@ -10,13 +10,12 @@
 #include <QPushButton>
 #include <QGroupBox>
 #include <QDoubleValidator>
-
-#include "source.hpp"
+#include "qt/input/data.hpp"
 
 class QtInputPump:public QFrame{
   Q_OBJECT
 private:
-  Source* source;
+  QtInputData* data;
   Pump* pump;
   QGroupBox* init_groupbox;
   QGroupBox* final_groupbox;
@@ -84,9 +83,11 @@ private:
   QDoubleValidator* double_validator;
   QDoubleValidator* double_amplitude_validator;
 public:
-  QtInputPump(Source*);
+  QtInputPump(QtInputData* data,Pump* pump=nullptr);
   ~QtInputPump();
   QWidget* validate();
+  void setPump();
+  void getPump();
 public slots:
   void emitRemove();
 signals:
@@ -95,7 +96,6 @@ signals:
 
 inline
 QtInputPump::~QtInputPump(){
-  source->removePump(pump);
   delete double_validator;
   delete double_amplitude_validator;
 }

src/qt/input_pump_tab.cpp

@@ -1,7 +1,7 @@
-#include "input_pump_tab.hpp"
+#include "qt/input/pump_tab.hpp"
 
-QtInputPumpTab::QtInputPumpTab(Source* _source):QWidget(){
-  source=_source;
+QtInputPumpTab::QtInputPumpTab(QtInputData* d):QWidget(){
+  data=d;
 
   main_layout=new QVBoxLayout;
   pumps_layout=new QVBoxLayout;
@@ -30,12 +30,19 @@ QtInputPumpTab::QtInputPumpTab(Source* _source):QWidget(){
 
 void
 QtInputPumpTab::addPump(){
-  QtInputPump* input_pump=new QtInputPump(source);
+  QtInputPump* input_pump=new QtInputPump(data);
   pumps_layout->addWidget(input_pump);
   connect(input_pump,&QtInputPump::remove,this,&QtInputPumpTab::removePump);
   emit sourcesChanged();
 }
 
+void
+QtInputPumpTab::addExistingPump(Pump* pump){
+  QtInputPump* input_pump=new QtInputPump(data,pump);
+  pumps_layout->addWidget(input_pump);
+  connect(input_pump,&QtInputPump::remove,this,&QtInputPumpTab::removePump);
+}
+
 void
 QtInputPumpTab::removePump(QtInputPump* input_pump){
   disconnect(input_pump,nullptr,nullptr,nullptr);
@@ -65,3 +72,10 @@ QtInputPumpTab::updateSources(){
   }
   emit sourcesChanged();
 }
+
+void
+QtInputPumpTab::getPumps(){
+  for(auto it=data->source.pumps.begin();it!=data->source.pumps.end();++it){
+   addExistingPump(*it);
+  }
+}

src/qt/input_pump_tab.hpp

@@ -13,14 +13,14 @@
 #include <QPalette>
 #include <QMessageBox>
 
-#include "qt/input_pump.hpp"
+#include "qt/input/pump.hpp"
 
-#include "source.hpp"
+#include "qt/input/data.hpp"
 
 class QtInputPumpTab:public QWidget{
   Q_OBJECT
 private:
-  Source* source;
+  QtInputData* data;
   QPushButton* add_button;
   QPushButton* refresh_button;
   QScrollArea* scroll_area;
@@ -30,11 +30,11 @@ private:
   QHBoxLayout* button_layout;
   QWidget* button_widget;
 public:
-  QtInputPumpTab(Source* source);
+  QtInputPumpTab(QtInputData* data);
   ~QtInputPumpTab();
   QWidget* validate();
-  void save(fstream& file);
-  void load(fstream& file);
+  void addExistingPump(Pump* pump);
+  void getPumps();
 public slots:
    void addPump();
    void removePump(QtInputPump*);

src/qt/input_tank.cpp

@@ -1,8 +1,13 @@
-#include "input_tank.hpp"
+#include "qt/input/tank.hpp"
 
-QtInputTank::QtInputTank(InitialState* _initial_state):QFrame(){
-  initial_state=_initial_state;
-  tank=initial_state->addTank();
+QtInputTank::QtInputTank(QtInputData* d,Tank* t):QFrame(){
+  data=d;
+  if(t==nullptr){
+    tank=data->addTank();
+  }
+  else{
+    tank=t;
+  }
 
   saturation_widget=new QWidget;
   left_right_widget=new QWidget;
@@ -14,19 +19,19 @@ QtInputTank::QtInputTank(InitialState* _initial_state):QFrame(){
   bottom_top_layout=new QHBoxLayout;
 
   saturation_label=new QLabel("Saturation: ");
-  saturation_input=new QLineEdit(QString::number(tank->saturation));
+  saturation_input=new QLineEdit();
   saturation_layout->addWidget(saturation_label);
   saturation_layout->addWidget(saturation_input);
   saturation_widget->setLayout(saturation_layout);
 
   left_label=new QLabel("Left: ");
   delta_left_label=new QLabel("  Left delta: ");
-  left_input=new QLineEdit(QString::number(tank->left));
-  delta_left_input=new QLineEdit(QString::number(tank->delta_left));
+  left_input=new QLineEdit();
+  delta_left_input=new QLineEdit();
   right_label=new QLabel("  Right: ");
   delta_right_label=new QLabel("  Right delta: ");
-  right_input=new QLineEdit(QString::number(tank->right));
-  delta_right_input=new QLineEdit(QString::number(tank->delta_right));
+  right_input=new QLineEdit();
+  delta_right_input=new QLineEdit();
 
   left_right_layout->addWidget(left_label);
   left_right_layout->addWidget(left_input,1);
@@ -40,12 +45,12 @@ QtInputTank::QtInputTank(InitialState* _initial_state):QFrame(){
 
   bottom_label=new QLabel("Bottom: ");
   delta_bottom_label=new QLabel("  Bottom delta: ");
-  bottom_input=new QLineEdit(QString::number(tank->bottom));
-  delta_bottom_input=new QLineEdit(QString::number(tank->delta_bottom));
+  bottom_input=new QLineEdit();
+  delta_bottom_input=new QLineEdit();
   top_label=new QLabel("  Top: ");
   delta_top_label=new QLabel("  Top delta: ");
-  top_input=new QLineEdit(QString::number(tank->top));
-  delta_top_input=new QLineEdit(QString::number(tank->delta_top));
+  top_input=new QLineEdit();
+  delta_top_input=new QLineEdit();
 
   bottom_top_layout->addWidget(bottom_label);
   bottom_top_layout->addWidget(bottom_input,1);
@@ -81,10 +86,20 @@ QtInputTank::QtInputTank(InitialState* _initial_state):QFrame(){
   setFrameShape(QFrame::Box);
 
   connect(remove_button,&QPushButton::clicked,this,&QtInputTank::emitRemove);
+
+  if(t==nullptr){
+    tank->bottom*=data->factor;
+    tank->top*=data->factor;
+    tank->delta_bottom*=data->factor;
+    tank->delta_top*=data->factor;
+  }
+  getTank();
+
 }
 
 QWidget*
 QtInputTank::validate(){
+
   if(not saturation_input->hasAcceptableInput()) return saturation_input;
   if(not left_input->hasAcceptableInput()) return left_input;
   if(not right_input->hasAcceptableInput()) return right_input;
@@ -94,14 +109,25 @@ QtInputTank::validate(){
   if(not delta_right_input->hasAcceptableInput()) return delta_right_input;
   if(not delta_bottom_input->hasAcceptableInput()) return delta_bottom_input;
   if(not delta_top_input->hasAcceptableInput()) return delta_top_input;
-  double s,l,r,t,b,sl,sr,st,sb;
-  s=saturation_input->text().toDouble();
+  double l,r,t,b;
   l=left_input->text().toDouble();
   r=right_input->text().toDouble();
   t=top_input->text().toDouble();
   b=bottom_input->text().toDouble();
   if(l>=r) return right_input;
   if(b>=t) return top_input;
+  setTank();
+  return nullptr;
+}
+
+void
+QtInputTank::setTank(){
+  double s,l,r,t,b,sl,sr,st,sb;
+  s=saturation_input->text().toDouble();
+  l=left_input->text().toDouble();
+  r=right_input->text().toDouble();
+  t=top_input->text().toDouble();
+  b=bottom_input->text().toDouble();
   sl=delta_left_input->text().toDouble();
   sr=delta_right_input->text().toDouble();
   st=delta_top_input->text().toDouble();
@@ -109,12 +135,24 @@ QtInputTank::validate(){
   tank->saturation=s;
   tank->left=l;
   tank->right=r;
-  tank->bottom=b;
-  tank->top=t;
+  tank->bottom=b*data->factor;
+  tank->top=t*data->factor;
   tank->delta_left=sl;
   tank->delta_right=sr;
-  tank->delta_bottom=sb;
-  tank->delta_top=st;
-  return nullptr;
+  tank->delta_bottom=sb*data->factor;
+  tank->delta_top=st*data->factor;
+}
 
+void
+QtInputTank::getTank(){
+  double factor_inv=1/data->factor;
+  saturation_input->setText(QString::number(tank->saturation));
+  left_input->setText(QString::number(tank->left));
+  right_input->setText(QString::number(tank->right));
+  top_input->setText(QString::number(tank->top*factor_inv));
+  bottom_input->setText(QString::number(tank->bottom*factor_inv));
+  delta_left_input->setText(QString::number(tank->delta_left));
+  delta_right_input->setText(QString::number(tank->delta_right));
+  delta_top_input->setText(QString::number(tank->delta_top*factor_inv));
+  delta_bottom_input->setText(QString::number(tank->delta_bottom*factor_inv));
 }

src/qt/input_tank.hpp

@@ -9,12 +9,12 @@
 #include <QPushButton>
 #include <QDoubleValidator>
 
-#include "initial_state.hpp"
+#include "qt/input/data.hpp"
 
 class QtInputTank:public QFrame{
   Q_OBJECT
 private:
-  InitialState* initial_state;
+  QtInputData* data;
   Tank* tank;
   QWidget* saturation_widget;
   QWidget* left_right_widget;
@@ -44,9 +44,11 @@ private:
   QLineEdit* delta_bottom_input;
   QDoubleValidator* double_validator;
 public:
-  QtInputTank(InitialState*);
+  QtInputTank(QtInputData* data,Tank* tank=nullptr);
   ~QtInputTank();
   QWidget* validate();
+  void setTank();
+  void getTank();
 public slots:
   void emitRemove();
 signals:
@@ -55,7 +57,7 @@ signals:
 
 inline
 QtInputTank::~QtInputTank(){
-  initial_state->removeTank(tank);
+  //model->removeTank(tank);
   delete double_validator;
 }
 

src/qt/input_time.cpp

@@ -1,12 +1,12 @@
-#include "input_time.hpp"
+#include "qt/input/time.hpp"
 
 QtInputTime::QtInputTime():QWidget(){
   main_layout=new QVBoxLayout;
   T_label=new QLabel("Total duration of the simulation (T)");
   nT_label=new QLabel("Number of time steps (nT)");
 
-  T_input=new QLineEdit(QString::number(Time::T));
-  nT_input=new QLineEdit(QString::number(Time::nT));
+  T_input=new QLineEdit();
+  nT_input=new QLineEdit();
 
   positive_double_validator=new QDoubleValidator;
   positive_double_validator->setBottom(0);
@@ -14,7 +14,7 @@ QtInputTime::QtInputTime():QWidget(){
   positive_int_validator->setBottom(1);
   T_input->setValidator(positive_double_validator);
   nT_input->setValidator(positive_int_validator);
-  
+
   int vspace=20;
   main_layout->addWidget(T_label);
   main_layout->addWidget(T_input);
@@ -24,6 +24,8 @@ QtInputTime::QtInputTime():QWidget(){
   main_layout->addSpacing(vspace);
   main_layout->addStretch();
   setLayout(main_layout);
+
+  getTime();
 }
 
 QWidget*
@@ -34,20 +36,13 @@ QtInputTime::validate(){
 }
 
 void
-QtInputTime::save(fstream& file){
-  double d=T_input->text().toDouble();
-  file.write((char*)&d,sizeof(double));
-  size_t s=nT_input->text().toULong();
-  file.write((char*)&s,sizeof(size_t));
+QtInputTime::setTime(){
+  Time::T=T_input->text().toDouble();
+  Time::nT=nT_input->text().toULong();
 }
 
 void
-QtInputTime::load(fstream& file){
-  double d;
-  size_t s;
-  file.read((char*)&d,sizeof(double));
-  T_input->setText(QString::number(d));
-  file.read((char*)&s,sizeof(size_t));
-  nT_input->setText(QString::number(s));
+QtInputTime::getTime(){
+  T_input->setText(QString::number(Time::T));
+  nT_input->setText(QString::number(Time::nT));
 }
-

src/qt/input_time.hpp

@@ -10,7 +10,7 @@
 #include <QIntValidator>
 #include <fstream>
 
-#include "time.hpp"
+#include "../../time.hpp"
 
 using namespace std;
 
@@ -28,8 +28,8 @@ public:
   QtInputTime();
   ~QtInputTime();
   QWidget* validate();
-  void save(fstream& file);
-  void load(fstream& file);
+  void setTime();
+  void getTime();
 };
 
 inline
@@ -37,4 +37,5 @@ QtInputTime::~QtInputTime(){
   delete positive_double_validator;
   delete positive_int_validator;
 }
+
 #endif

src/qt/input/view.cpp

@@ -0,0 +1,321 @@
+#include "qt/input/view.hpp"
+
+QtInputView::QtInputView(QtInputData* m){
+  data=m;
+  selected=p_none;
+  status=Other;
+}
+
+void
+QtInputView::initializeGL(){
+  glClearColor(1,1,1,1);
+  glEnable(GL_DEPTH_TEST);
+  glEnable(GL_LIGHT0);
+  glEnable(GL_LIGHTING);
+  glColorMaterial(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE);
+  glEnable(GL_COLOR_MATERIAL);
+}
+
+void
+QtInputView::resizeGL(int w,int h){
+  glViewport(0,0,w,h);
+  glMatrixMode(GL_PROJECTION);
+  glLoadIdentity();
+  gluOrtho2D(0,1,0,1);
+  glMatrixMode(GL_MODELVIEW);
+  glLoadIdentity();
+}
+
+void
+QtInputView::paintGL(){
+  if(status==Geom or status==Init){
+    glBegin(GL_QUADS);
+    double xs=double(pointSize)/width();
+    double ys=double(pointSize)/height();
+    glColor3f(0,0,0);
+    size_t imin=(status==Geom)?0:2*QtInputData::np;
+    size_t imax=(status==Geom)?2*QtInputData::np:3*QtInputData::np;
+    for(size_t i=imin;i<imax;++i){
+      Point& P=data->point[i];
+      double x=P.x;
+      double y=P.y;
+      glVertex3f(x-xs,y-ys,1);
+      glVertex3f(x+xs,y-ys,1);
+      glVertex3f(x+xs,y+ys,1);
+      glVertex3f(x-xs,y+ys,1);
+    }
+    if(status==Init){
+      glColor3f(0.3,0.3,1);
+      for(size_t i=3*QtInputData::np;i<3*QtInputData::np+QtInputData::np_ov;++i){
+        Point& P=data->point[i];
+        double x=P.x;
+        double y=P.y;
+        glVertex3f(x-xs,y-ys,1);
+        glVertex3f(x+xs,y-ys,1);
+        glVertex3f(x+xs,y+ys,1);
+        glVertex3f(x-xs,y+ys,1);
+      }
+    }
+    glEnd();
+  }
+  for(auto it=data->initial_state->tanks.begin();it!=data->initial_state->tanks.end();++it){
+    paintTank(*it);
+  }
+  for(auto it=data->source.pumps.begin();it!=data->source.pumps.end();++it){
+    paintPump(*it);
+  }
+  for(auto it=data->source.clouds.begin();it!=data->source.clouds.end();++it){
+    paintCloud(*it);
+  }
+  glColor3f(0.6,0.6,0.6);
+  drawSpline(data->spline_bot);
+  glColor3f(0.6,0.3,0);
+  drawSpline(data->spline_soil);
+  glColor3f(0,1,1);
+  drawSpline(data->spline_sat);
+  Geometry& G=data->geometry;
+  glBegin(GL_TRIANGLES);
+  for(size_t i=0;i<G.nX-1;++i){
+    size_t max_left=G.nZ[i];
+    size_t max_right=G.nZ[i+1];
+    size_t left=0;
+    size_t right=0;
+      while(left<max_left-1 and right<max_right-1){
+      if(G.Z[i][left+1]<G.Z[i+1][right+1]){
+        drawTriangle(i,left,i+1,right,i,left+1);
+        ++left;
+      }
+      else{
+        drawTriangle(i,left,i+1,right,i+1,right+1);
+        ++right;
+      }
+    }
+    if(left==max_left-1){
+      while(right<max_right-1){
+        drawTriangle(i,left,i+1,right,i+1,right+1);
+        ++right;
+      }
+    }
+    if(right==max_right-1){
+      while(left<max_left-1){
+        drawTriangle(i,left,i+1,right,i,left+1);
+        ++left;
+      }
+    }
+  }
+  glEnd();
+  drawOverland();
+}
+
+void
+QtInputView::paintTank(Tank* tank){
+  double factor_inv=1.0/data->factor;
+  double l=tank->left;
+  double r=tank->right;
+  double b=tank->bottom*factor_inv;
+  double t=tank->top*factor_inv;
+  glLineWidth(2);
+  glColor3f(1,0,1);
+  glBegin(GL_LINE_LOOP);
+  glVertex2f(l,b);
+  glVertex2f(r,b);
+  glVertex2f(r,t);
+  glVertex2f(l,t);
+  glEnd();
+}
+
+void
+QtInputView::paintPump(Pump* pump){
+  double factor_inv=1.0/data->factor;
+  double a=pump->get_amplitude(time);
+  if(a==0) return;
+  double l=pump->get_left(time);
+  double r=pump->get_right(time);
+  double b=pump->get_bottom(time)*factor_inv;
+  double t=pump->get_top(time)*factor_inv;
+  double dl=pump->get_left_delta(time);
+  double dr=pump->get_right_delta(time);
+  double db=pump->get_bottom_delta(time)*factor_inv;
+  double dt=pump->get_top_delta(time)*factor_inv;
+  glLineWidth(2);
+  if(a>0) glColor3f(0,0.6,1);
+  else glColor3f(1,0.2,0);
+  glLineStipple(3, 0xAAAA);
+  glBegin(GL_LINE_LOOP);
+  glVertex2f(l,b);
+  glVertex2f(r,b);
+  glVertex2f(r,t);
+  glVertex2f(l,t);
+  glEnd();
+  glLineWidth(1);
+  glEnable(GL_LINE_STIPPLE);
+  glBegin(GL_LINE_LOOP);
+  glVertex2f(l-dl,b-db);
+  glVertex2f(r+dr,b-db);
+  glVertex2f(r+dr,t+dt);
+  glVertex2f(l-dl,t+dt);
+  glEnd();
+  glDisable(GL_LINE_STIPPLE);
+}
+
+void
+QtInputView::paintCloud(Cloud* cloud){
+  double a=cloud->get_amplitude(time);
+  double amax=cloud->get_amplitude_max();
+  if(a==0) return;
+  double l=cloud->get_left(time);
+  double r=cloud->get_right(time);
+  double dl=cloud->get_left_delta(time);
+  double dr=cloud->get_right_delta(time);
+  double middle=0.97;
+  double thick=0.03*a/amax;
+  double top=middle+thick;
+  double bottom=middle-thick;
+
+  glLineWidth(1);
+  glColor3f(0.5,0.5,0.5);
+  glLineWidth(1);
+  glBegin(GL_POLYGON);
+  glVertex2f(l-dl,middle);
+  glVertex2f(l,top);
+  glVertex2f(r,top);
+  glVertex2f(r+dr,middle);
+  glVertex2f(r,bottom);
+  glVertex2f(l,bottom);
+  glEnd();
+}
+
+void
+QtInputView::setColor(size_t ix,size_t iz){
+  double p=getP(ix,iz);
+  double s=Physics::s(p);
+  double r=(1-s);
+  double g=0.7*(1-s);
+  double b=0.8*s+0.4*(1-s);
+  glColor3f(r,g,b);
+}
+
+void
+QtInputView::drawTriangle(size_t ix1,size_t iz1,size_t ix2,size_t iz2,size_t ix3,size_t iz3){
+  Geometry& G=data->geometry;
+  double factor_inv=1/data->factor;
+  double dX=1/double(G.nX-1);
+  double x1=ix1*dX;
+  double x2=ix2*dX;
+  double x3=ix3*dX;
+  double y1=G.Z[ix1][iz1]*factor_inv;
+  double y2=G.Z[ix2][iz2]*factor_inv;
+  double y3=G.Z[ix3][iz3]*factor_inv;
+  setColor(ix1,iz1);
+  glVertex3f(x1,y1,0);
+  setColor(ix2,iz2);
+  glVertex3f(x2,y2,0);
+  setColor(ix3,iz3);
+  glVertex3f(x3,y3,0);
+}
+
+void
+QtInputView::drawOverland(){
+  Geometry& G=data->geometry;
+  double dX=1/double(G.nX-1);
+  double factor_inv=1/data->factor;
+  glColor3f(0.5,0.5,1);
+  glBegin(GL_QUADS);
+  for(size_t i=0;i<G.nX-1;++i){
+    glVertex3f(i*dX,G.hsoil[i]*factor_inv,0);
+    glVertex3f(i*dX,data->initial_state->hov[i]*factor_inv,0);
+    glVertex3f(i*dX+dX,data->initial_state->hov[i+1]*factor_inv,0);
+    glVertex3f(i*dX+dX,G.hsoil[i+1]*factor_inv,0);
+  }
+  glEnd();
+}
+
+void
+QtInputView::mousePressEvent(QMouseEvent* event){
+  if(status==Geom or status==Init){
+    size_t p=findPoint(event->x(),height()-event->y());
+    if(status==Geom){
+      selected=(p<2*QtInputData::np)?p:p_none;
+    }
+    else{
+      selected=(p>=2*QtInputData::np)?p:p_none;
+    }
+  }
+}
+
+void
+QtInputView::mouseMoveEvent(QMouseEvent* event){
+  if(selected!=p_none){
+    double mx=double(event->x())/width();
+    double my=1-double(event->y())/height();
+    moveSelected(mx,my);
+  }
+}
+
+size_t
+QtInputView::findPoint(int x,int y){
+  double w=width();
+  double h=height();
+  for(size_t i=0;i<3*QtInputData::np+QtInputData::np_ov;++i){
+    int px=data->point[i].x*w;
+    int py=data->point[i].y*h;
+    if(abs(x-px)<pointSize and abs(y-py)<pointSize) return i;
+  }
+  return p_none;
+}
+
+void
+QtInputView::drawSpline(Spline& S){
+  double w=width();
+  glLineWidth(3);
+  glBegin(GL_LINE_STRIP);
+  for(int i=0;i<=w;i+=2){
+    double x=i/w;
+    glVertex2f(x,S(x));
+  }
+  glEnd();
+}
+
+void
+QtInputView::moveSelected(double x,double y){
+  double xp=data->point[selected].x;
+  double yp=data->point[selected].y;
+  if((selected%QtInputData::np)!=0 and (selected%QtInputData::np)!=QtInputData::np-1){
+    float xmin=data->point[selected-1].x+min_d;
+    float xmax=data->point[selected+1].x-min_d;
+    if(xmin<=x and x<=xmax){
+      data->point[selected].x=x;
+    }
+  }
+  if(0<=y and y<=1){
+    data->point[selected].y=y;
+  }
+  if(selected<QtInputData::np) data->updateSplineBot();
+  else if(selected<2*QtInputData::np) data->updateSplineSoil();
+  else if(selected<3*QtInputData::np) data->updateSplineSat();
+  else data->updateOverland();
+
+  if(selected<2*QtInputData::np){
+    double w=width();
+    bool ok=true;
+    for(int i=0;i<=w;++i){
+      double x=i/w;
+      double ybot=data->evalSplineBot(x);
+      double ysoil=data->evalSplineSoil(x);
+      if(ybot<min_d or ybot>1-min_d or ysoil<min_d or ysoil>1-min_d or ybot>ysoil-min_d){
+        data->point[selected].x=xp;
+        data->point[selected].y=yp;
+        if(selected<QtInputData::np) data->updateSplineBot();
+        else data->updateSplineSoil();
+        update();
+        return;
+      }
+    }
+  }
+  update();
+}
+
+void
+QtInputView::mouseReleaseEvent(QMouseEvent* event){
+  if(selected<3*QtInputData::np) emit geometryChanged();
+}

src/qt/input_view.hpp

@@ -1,59 +1,62 @@
 #ifndef QT_INPUT_VIEW_HPP
 #define QT_INPUT_VIEW_HPP
 
-#include "initial_state.hpp"
-#include "input_geometry.hpp"
-#include "view.hpp"
+#include <QPainter>
+#include <QOpenGLWidget>
+#include <QMouseEvent>
+#include <GL/glut.h>
+
 #include "math/point.hpp"
 #include "math/spline.hpp"
-#include "geometry.hpp"
-#include "time.hpp"
-#include "source.hpp"
+#include "qt/input/geometry.hpp"
+#include "qt/input/initial_state.hpp"
+#include "qt/input/data.hpp"
 
 using namespace std;
 
-class QtInputView:public QtView{
+class QtInputView: public QOpenGLWidget{
   Q_OBJECT
 public:
   enum Status{Geom,Init,Sources,Other};
 private:
-  QtInputGeometry* input_geometry;
-  InitialState* initial_state;
-  Source* source;
   Status status;
-  static constexpr size_t const np=10;
+  QtInputData* data;
   static constexpr int const pointSize=6;
   static constexpr double const min_d=0.01;
+  static constexpr size_t const p_none=3*QtInputData::np+QtInputData::np_ov;
   size_t selected;
-  Point point[3*np];
-  Spline hsoil;
-  Spline hbot;
-  Spline hsat;
   double time;
-  void initPoints();
   size_t findPoint(int x,int y);
   void drawSpline(Spline& S);
+  void drawOverland();
   void moveSelected(double x,double y);
 public:
-  QtInputView(QtInputGeometry* input_geometry,Geometry* geometry,InitialState* initial_state,Source* source);
+  QtInputView(QtInputData* data);
   ~QtInputView();
   void paintGL() override;
   void mousePressEvent(QMouseEvent* event);
   void mouseMoveEvent(QMouseEvent* event);
   void mouseReleaseEvent(QMouseEvent* event);
   void setStatus(Status status);
+  void initializeGL();
+  void resizeGL(int x,int h);
+  void drawTriangle(size_t ix1,size_t iz1,size_t ix2,size_t iz2,size_t ix3,size_t iz3);
+  void setColor(size_t ix,size_t iz);
   double getP(size_t ix,size_t iz);
   void paintTank(Tank*);
   void paintPump(Pump*);
   void paintCloud(Cloud*);
+
 public slots:
-  void updateGeometry();
-  void updateInitialState();
-  void updateSource();
+  //void updateGeometry();
+  //void updateInitialState(double factor);
+  //void updateSource();
   void setTime(int v);
+signals:
+  void geometryChanged();
+  void timeChanged();
 };
 
-
 inline void
 QtInputView::setStatus(Status _status){
   status=_status;
@@ -61,8 +64,7 @@ QtInputView::setStatus(Status _status){
 
 inline double
 QtInputView::getP(size_t ix,size_t iz){
-  return initial_state->Pinit[ix][iz];
-
+  return data->initial_state->Pinit[ix][iz];
 }
 
 inline
@@ -72,6 +74,9 @@ QtInputView::~QtInputView(){
 inline void
 QtInputView::setTime(int v){
   time=double(v)/99;
-  updateSource();
+  emit timeChanged();
+
 }
+
+
 #endif

src/qt/input_view.cpp

@@ -1,255 +0,0 @@
-#include "qt/input_view.hpp"
-
-QtInputView::QtInputView(QtInputGeometry* _input_geometry,Geometry* geometry,InitialState* _initial_state,Source* _source):QtView(geometry){
-  input_geometry=_input_geometry;
-  initial_state=_initial_state;
-  source=_source;
-  selected=3*np;
-  status=Other;
-  initPoints();
-  hbot.setPoints(point,np);
-  hbot.compute();
-  hsoil.setPoints(&point[np],np);
-  hsoil.compute();
-  hsat.setPoints(&point[2*np],np);
-  hsat.compute();
-  initial_state->update(*geometry,hsat);
-}
-
-void
-QtInputView::initPoints(){
-  for(size_t i=0;i<np;++i){
-    double x=double(i)/(np-1);
-    point[i].x=x;
-    point[i].y=0.2;
-    point[np+i].x=x;
-    point[np+i].y=0.8;
-    point[2*np+i].x=x;
-    point[2*np+i].y=0.5;
-  }
-}
-
-void
-QtInputView::paintGL(){
-  if(status==Geom or status==Init){
-    glBegin(GL_QUADS);
-    double xs=double(pointSize)/width();
-    double ys=double(pointSize)/height();
-    glColor3f(0,0,0);
-    size_t imin=(status==Geom)?0:2*np;
-    size_t imax=(status==Geom)?2*np:3*np;
-    for(size_t i=imin;i<imax;++i){
-      Point& P=point[i];
-      double x=P.x;
-      double y=P.y;
-      glVertex3f(x-xs,y-ys,1);
-      glVertex3f(x+xs,y-ys,1);
-      glVertex3f(x+xs,y+ys,1);
-      glVertex3f(x-xs,y+ys,1);
-    }
-    glEnd();
-  }
-  for(auto it=initial_state->tanks.begin();it!=initial_state->tanks.end();++it){
-    paintTank(*it);
-  }
-  for(auto it=source->pumps.begin();it!=source->pumps.end();++it){
-    paintPump(*it);
-  }
-  for(auto it=source->clouds.begin();it!=source->clouds.end();++it){
-    paintCloud(*it);
-  }
-  glColor3f(0.6,0.6,0.6);
-  drawSpline(hbot);
-  glColor3f(0.6,0.3,0);
-  drawSpline(hsoil);
-  glColor3f(0,1,1);
-  drawSpline(hsat);
-  QtView::paintGL();
-}
-
-void
-QtInputView::paintTank(Tank* tank){
-  double l=tank->left;
-  double r=tank->right;
-  double b=tank->bottom;
-  double t=tank->top;
-  glLineWidth(2);
-  glColor3f(1,0,1);
-  glBegin(GL_LINE_LOOP);
-  glVertex2f(l,b);
-  glVertex2f(r,b);
-  glVertex2f(r,t);
-  glVertex2f(l,t);
-  glEnd();
-}
-
-void
-QtInputView::paintPump(Pump* pump){
-  double a=pump->get_amplitude(time);
-  if(a==0) return;
-  double l=pump->get_left(time);
-  double r=pump->get_right(time);
-  double b=pump->get_bottom(time);
-  double t=pump->get_top(time);
-  double dl=pump->get_left_delta(time);
-  double dr=pump->get_right_delta(time);
-  double db=pump->get_bottom_delta(time);
-  double dt=pump->get_top_delta(time);
-  glLineWidth(2);
-  if(a>0){
-    glColor3f(0,0.6,1);
-  }
-  else{
-    glColor3f(1,0.2,0);
-  }
-  glLineStipple(3, 0xAAAA);
-  glBegin(GL_LINE_LOOP);
-  glVertex2f(l,b);
-  glVertex2f(r,b);
-  glVertex2f(r,t);
-  glVertex2f(l,t);
-  glEnd();
-  glLineWidth(1);
-  glEnable(GL_LINE_STIPPLE);
-  glBegin(GL_LINE_LOOP);
-  glVertex2f(l-dl,b-db);
-  glVertex2f(r+dr,b-db);
-  glVertex2f(r+dr,t+dt);
-  glVertex2f(l-dl,t+dt);
-  glEnd();
-  glDisable(GL_LINE_STIPPLE);
-}
-
-void
-QtInputView::paintCloud(Cloud* cloud){
-  double a=cloud->get_amplitude(time);
-  double amax=cloud->get_amplitude_max();
-  if(a==0) return;
-  double l=cloud->get_left(time);
-  double r=cloud->get_right(time);
-  double dl=cloud->get_left_delta(time);
-  double dr=cloud->get_right_delta(time);
-  double middle=0.97;
-  double thick=0.03*a/amax;
-  double top=middle+thick;
-  double bottom=middle-thick;
-
-  glLineWidth(1);
-  glColor3f(0.5,0.5,0.5);
-  glLineWidth(1);
-  //glEnable(GL_LINE_STIPPLE);
-  glBegin(GL_POLYGON);
-  glVertex2f(l-dl,middle);
-  glVertex2f(l,top);
-  glVertex2f(r,top);
-  glVertex2f(r+dr,middle);
-  glVertex2f(r,bottom);
-  glVertex2f(l,bottom);
-  glEnd();
-  //glDisable(GL_LINE_STIPPLE);
-}
-
-void
-QtInputView::mousePressEvent(QMouseEvent* event){
-  if(status==Geom or status==Init){
-    size_t p=findPoint(event->x(),height()-event->y());
-    if(status==Geom){
-      selected=(p<2*np)?p:3*np;
-    }
-    else{
-      selected=(p>=2*np)?p:3*np;
-    }
-  }
-}
-
-void
-QtInputView::mouseMoveEvent(QMouseEvent* event){
-  if(selected<3*np){
-    double mx=double(event->x())/width();
-    double my=1-double(event->y())/height();
-    moveSelected(mx,my);
-  }
-}
-
-size_t
-QtInputView::findPoint(int x,int y){
-  double w=width();
-  double h=height();
-  for(size_t i=0;i<3*np;++i){
-    int px=point[i].x*w;
-    int py=point[i].y*h;
-    if(abs(x-px)<pointSize and abs(y-py)<pointSize) return i;
-  }
-  return 3*np;
-}
-
-void
-QtInputView::updateGeometry(){
-  geometry->update(input_geometry->get_lX(),input_geometry->get_nX(),input_geometry->get_depth(),input_geometry->get_nZ_max(),hsoil,hbot);
-  initial_state->update(*geometry,hsat);
-  update();
-}
-
-void
-QtInputView::updateInitialState(){
-  update();
-}
-
-void
-QtInputView::updateSource(){
-  update();
-}
-void
-QtInputView::drawSpline(Spline& S){
-  double w=width();
-  glLineWidth(3);
-  glBegin(GL_LINE_STRIP);
-  for(int i=0;i<=w;i+=2){
-    double x=i/w;
-    glVertex2f(x,S(x));
-  }
-  glEnd();
-}
-
-void
-QtInputView::moveSelected(double x,double y){
-  double xp=point[selected].x;
-  double yp=point[selected].y;
-  if((selected%np)!=0 and (selected%np)!=np-1){
-    float xmin=point[selected-1].x+min_d;
-    float xmax=point[selected+1].x-min_d;
-    if(xmin<=x and x<=xmax){
-      point[selected].x=x;
-    }
-  }
-  if(0<=y and y<=1){
-    point[selected].y=y;
-  }
-  if(selected<np) hbot.compute();
-  else if(selected<2*np) hsoil.compute();
-  else hsat.compute();
-
-  if(selected<2*np){
-    double w=width();
-    bool ok=true;
-    for(int i=0;i<=w;++i){
-      double x=i/w;
-      double ybot=hbot(x);
-      double ysoil=hsoil(x);
-      if(ybot<min_d or ybot>1-min_d or ysoil<min_d or ysoil>1-min_d or ybot>ysoil-min_d){
-        point[selected].x=xp;
-        point[selected].y=yp;
-        if(selected<np) hbot.compute();
-        else hsoil.compute();
-        update();
-        return;
-      }
-    }
-  }
-  update();
-}
-
-void
-QtInputView::mouseReleaseEvent(QMouseEvent* event){
-  if(selected<3*np) updateGeometry();
-}

src/qt/obsolete/geometry_spline.cpp

@@ -0,0 +1,72 @@
+#include "geometry_spline.hpp"
+
+void
+GeometrySpline::initDefault(){
+  lX=10;
+  nX=100;
+  nX_max=nmax_Qt;
+  dX=lX/(nX-1);
+  hsoil=new double[nmax_Qt];
+  dhsoil=new double[nmax_Qt];
+  hbot=new double[nmax_Qt];
+  dhbot=new double[nmax_Qt];
+  nZ=new size_t[nmax_Qt];
+  dZ=new double[nmax_Qt];
+  Z=new double*[nmax_Qt];
+  for(size_t i=0;i<nmax_Qt;++i){
+    Z[i]=new double[nmax_Qt];
+  }
+  for(size_t i=0;i<nX;++i){
+    hsoil[i]=0.8;
+    hbot[i]=0.2;
+    dhsoil[i]=0;
+    dhbot[i]=0;
+  }
+  depth=5;
+  nZ_max=100;
+  factor=5/0.6;
+  initZ(depth/nZ_max,false);
+}
+
+void
+GeometrySpline::update(double _lX,size_t _nX,double _depth,size_t _nZ_max,Spline& Ssoil,Spline& Sbot){
+  lX=_lX;
+  depth=_depth;
+  nX=_nX;
+  nZ_max=_nZ_max;
+  dX=lX/(nX-1);
+  double hs_max=-inf,hb_min=inf;
+  double v;
+  for(size_t k=0;k<nX;++k){
+    double x=k*dX/lX;
+    v=hsoil[k]=Ssoil(x);
+    hs_max=max(v,hs_max);
+    dhsoil[k]=Ssoil.derivate(x);
+    v=hbot[k]=Sbot(x);
+    hb_min=min(v,hb_min);
+    dhbot[k]=Sbot.derivate(x);
+  }
+  double dZ_avg=depth/_nZ_max;
+  factor=depth/(hs_max-hb_min);
+  for(size_t k=0;k<nX;++k){
+    hsoil[k]*=factor;
+    hbot[k]*=factor;
+  }
+  initZ(dZ_avg,false);
+}
+
+GeometrySpline::~GeometrySpline(){
+  if(hsoil!=nullptr){
+    delete[] hsoil;
+    delete[] dhsoil;
+    delete[] hbot;
+    delete[] dhbot;
+    delete[] dZ;
+    delete[] nZ;
+    for(size_t i=0;i<nX_max;++i){
+      delete[] Z[i];
+    }
+    delete[] Z;
+  }
+  hsoil=nullptr;
+}

src/qt/obsolete/geometry_spline.hpp

@@ -0,0 +1,23 @@
+#ifndef GEOMETRY_SPLINE_HPP
+#define GEOMETRY_SPLINE_HPP
+
+#include "geometry.hpp"
+
+class GeometrySpline:public Geometry{
+public:
+  size_t nX_max;
+  double factor;
+  static constexpr size_t const nmax_Qt=400;
+  GeometrySpline();
+  ~GeometrySpline();
+  void initDefault(Geometry* ge);
+  void update(double _lX,size_t _nX,double _depth,size_t _nZ_max,Spline& hbot,Spline& hsoil);
+};
+
+inline
+GeometrySpline::GeometrySpline():Geometry(){
+
+}
+
+
+#endif

src/qt/obsolete/initial_state_spline.cpp

@@ -0,0 +1,26 @@
+#include "initial_state_spline.hpp"
+
+void
+InitialStateSpline::update(Spline* _Ssat){
+  Ssat=_Ssat;
+  if(hsat!=nullptr){
+    delete[] hsat;
+    for(size_t i=0;i<geometry_spline->nX;++i){
+      delete[] Pinit[i];
+    }
+    delete[] Pinit;
+  }
+  size_t nX=geometry_spline->nX;
+  hsat=new double[nX];
+  Pinit=new double*[nX];
+
+  for(size_t i=0;i<nX;++i){
+    size_t nZ=geometry_spline->nZ[i];
+    Pinit[i]=new double[nZ];
+  }
+  for(size_t i=0;i<geometry_spline->nX;++i){
+    double x=double(i)/(geometry_spline->nX-1);
+    hsat[i]=(*Ssat)(x)*geometry_spline->factor;
+  }
+  updatePressure();
+}

src/qt/obsolete/initial_state_spline.hpp

@@ -0,0 +1,21 @@
+#ifndef INITIAL_STATE_SPLINE_HPP
+#define INITIAL_STATE_SPLINE_HPP
+
+#include "initial_state.hpp"
+#include "geometry_spline.hpp"
+
+class InitialStateSpline:public InitialState{
+private:
+  Spline* Ssat;
+  GeometrySpline* geometry_spline;
+public:
+  InitialStateSpline(GeometrySpline* geometry);
+  void update(Spline* Ssat);
+};
+
+inline
+InitialStateSpline::InitialStateSpline(GeometrySpline* g):InitialState((Geometry*)g){
+  geometry_spline=g;
+}
+
+#endif

src/qt/obsolete/input_overland.hpp

@@ -0,0 +1 @@
+

src/qt/view.cpp

@@ -1,7 +1,5 @@
 #include "view.hpp"
 
-
-
 void
 QtView::initializeGL(){
   glClearColor(1,1,1,1);
@@ -24,14 +22,15 @@ QtView::resizeGL(int w,int h){
 
 void
 QtView::paintGL(){
+  Geometry& G=data->geometry;
   glBegin(GL_TRIANGLES);
-  for(size_t i=0;i<geometry->nX-1;++i){
-    size_t max_left=geometry->nZ[i];
-    size_t max_right=geometry->nZ[i+1];
+  for(size_t i=0;i<G.nX-1;++i){
+    size_t max_left=G.nZ[i];
+    size_t max_right=G.nZ[i+1];
     size_t left=0;
     size_t right=0;
-    while(left<max_left-1 and right<max_right-1){
-      if(geometry->Z[i][left+1]<geometry->Z[i+1][right+1]){
+      while(left<max_left-1 and right<max_right-1){
+      if(G.Z[i][left+1]<G.Z[i+1][right+1]){
         drawTriangle(i,left,i+1,right,i,left+1);
         ++left;
       }
@@ -54,6 +53,7 @@ QtView::paintGL(){
     }
   }
   glEnd();
+  drawOverland();
 }
 
 void
@@ -68,18 +68,44 @@ QtView::setColor(size_t ix,size_t iz){
 
 void
 QtView::drawTriangle(size_t ix1,size_t iz1,size_t ix2,size_t iz2,size_t ix3,size_t iz3){
-  double dX=1/double(geometry->nX-1);
+  Geometry& G=data->geometry;
+  double factor_inv=1/data->factor;
+  double dX=1/double(G.nX-1);
   double x1=ix1*dX;
   double x2=ix2*dX;
   double x3=ix3*dX;
-  double y1=geometry->Z[ix1][iz1]/geometry->factor;
-  double y2=geometry->Z[ix2][iz2]/geometry->factor;
-  double y3=geometry->Z[ix3][iz3]/geometry->factor;
+  double y1=G.Z[ix1][iz1]*factor_inv;
+  double y2=G.Z[ix2][iz2]*factor_inv;
+  double y3=G.Z[ix3][iz3]*factor_inv;
   setColor(ix1,iz1);
   glVertex3f(x1,y1,0);
   setColor(ix2,iz2);
   glVertex3f(x2,y2,0);
   setColor(ix3,iz3);
   glVertex3f(x3,y3,0);
+}
+
+void
+QtView::drawOverland(){
+  Geometry& G=data->geometry;
+  double dX=1/double(G.nX-1);
+  double factor_inv=1/data->factor;
+  glColor3f(0.5,0.5,1);
+  glBegin(GL_QUADS);
+  for(size_t i=0;i<G.nX-1;++i){
+    glVertex3f(i*dX,G.hsoil[i]*factor_inv,0);
+    glVertex3f(i*dX,data->initial_state->hov[i]*factor_inv,0);
+    glVertex3f(i*dX+dX,data->initial_state->hov[i+1]*factor_inv,0);
+    glVertex3f(i*dX+dX,G.hsoil[i+1]*factor_inv,0);
+  }
+  glEnd();
+/*  glColor3f(0.2,0.2,1);
+
+  glLineWidth(2);
+  glBegin(GL_LINE_STRIP);
+  for(size_t i=0;i<G.nX;++i){
+    glVertex2f(i*dX,data->initial_state.hov[i]*factor_inv);
+  }
+  glEnd();*/
 
 }

src/qt/view.hpp

@@ -10,16 +10,16 @@
 
 #include <GL/glut.h>
 
-#include "geometry.hpp"
-#include "physics.hpp"
+#include "input_data.hpp"
 
 using namespace std;
 
 class QtView:public QOpenGLWidget{
 protected:
-  Geometry* geometry;
+  InputData* data;
+  void drawOverland();
 public:
-  QtView(Geometry* geometry);
+  QtView(InputData* data);
   void setGeometry(Geometry* geometry);
   void initializeGL();
   void paintGL();
@@ -31,8 +31,8 @@ public:
 };
 
 inline
-QtView::QtView(Geometry* _geometry):QOpenGLWidget(){
-  geometry=_geometry;
+QtView::QtView(InputData* d):QOpenGLWidget(){
+  data=d;
 }
 
 

src/qt/output/kernel.cpp

@@ -0,0 +1,18 @@
+#include "qt/output/kernel.hpp"
+
+QtKernel::QtKernel(string filename):Kernel(filename){
+}
+
+void
+QtKernel::run(){
+  QTime start=QTime::currentTime();
+  cout<<"QtKernel::run -> start"<<endl;
+  for(size_t i=1;i<Time::nT;++i){
+    next();
+  }
+  cout<<"QtKernel::run -> end"<<endl;
+  QTime end=QTime::currentTime();
+  float ms=start.msecsTo(end);
+  cout<<"Duration : "<<ms/1000<<endl;
+
+}

src/qt/output/kernel.hpp

@@ -0,0 +1,24 @@
+#ifndef QT_KERNEL_HPP
+#define QT_KERNEL_HPP
+
+#include <QThread>
+#include <QTime>
+#include "../../kernel.hpp"
+
+class QtKernel:public QObject,public Kernel{
+  Q_OBJECT
+public:
+  QtKernel(string filename);
+  ~QtKernel();
+public slots:
+  void run();
+signals:
+  void finished(int);
+};
+
+inline
+QtKernel::~QtKernel(){
+}
+
+
+#endif

src/qt/output/main_window.cpp

@@ -0,0 +1,7 @@
+#include "qt/output/main_window.hpp"
+
+QtMainWindow::QtMainWindow():QMainWindow(){
+  string filename="inputs/kernel-2.input";
+  output=new QtOutput(filename);
+  setCentralWidget(output);
+}

src/qt/output/main_window.hpp

@@ -0,0 +1,26 @@
+#ifndef QT_MAINWINDOW_HPP
+#define QT_MAINWINDOW_HPP
+
+#include <QMainWindow>
+#include <QMenu>
+#include <QMenuBar>
+#include <QAction>
+#include <QApplication>
+
+#include "qt/output/output.hpp"
+#include "qt/output/kernel.hpp"
+
+class QtMainWindow:public QMainWindow{
+  Q_OBJECT
+private:
+  QtOutput* output;
+  Kernel* kernel;
+public:
+  QtMainWindow();
+  ~QtMainWindow();
+};
+
+inline
+QtMainWindow::~QtMainWindow(){
+}
+#endif

src/qt/output/output.cpp

@@ -0,0 +1,49 @@
+#include "qt/output/output.hpp"
+
+QtOutput::QtOutput(string filename):QWidget(){
+
+  kernel=new QtKernel(filename);
+  kernel_thread=new QThread();
+  kernel->moveToThread(kernel_thread);
+
+  main_layout=new QVBoxLayout;
+  time_bar=new QScrollBar(Qt::Horizontal);
+  output_view=new QtOutputView(kernel);
+  info_widget=new QWidget;
+  info_layout=new QGridLayout;
+  time_label=new QLabel;
+
+  time_bar->setMinimum(0);
+  time_bar->setMaximum(Time::nT-1);
+  info_layout->addWidget(time_label);
+  info_widget->setLayout(info_layout);
+  main_layout->addWidget(output_view,1);
+  main_layout->addWidget(time_bar);
+  main_layout->addWidget(info_widget);
+  setLayout(main_layout);
+  step=0;
+  update();
+  connect(time_bar,&QScrollBar::valueChanged,this,&QtOutput::time_change);
+  connect(kernel_thread,&QThread::started,kernel,&QtKernel::run);
+  kernel_thread->start();
+}
+
+QtOutput::~QtOutput(){
+}
+
+void
+QtOutput::update(){
+  update_infos();
+  output_view->draw(step);
+}
+
+void
+QtOutput::update_infos(){
+  time_label->setText("Time : "+QString::number(step*Time::dT));
+}
+
+void
+QtOutput::time_change(int val){
+  step=val;
+  update();
+}

src/qt/view_solution.hpp

@@ -1,29 +1,32 @@
-#ifndef VIEW_SOLUTION_HPP
-#define VIEW_SOLUTION_HPP
+#ifndef QT_OUTPUT_HPP
+#define QT_OUTPUT_HPP
 
 #include <QWidget>
 #include <QVBoxLayout>
 #include <QScrollBar>
 #include <QGridLayout>
 #include <QLabel>
+#include <QThread>
 
-#include "qt/view_solution_geometry.hpp"
-#include "kernel.hpp"
+#include "qt/output/output_view.hpp"
+#include "qt/output/kernel.hpp"
 
-class QtViewSolution:public QWidget{
+class QtOutput:public QWidget{
   Q_OBJECT
 private:
+  QtKernel* kernel;
+  QThread* kernel_thread;
+
   QVBoxLayout* main_layout;
   QScrollBar* time_bar;
   QWidget* info_widget;
   QGridLayout* info_layout;
   QLabel* time_label;
-  Kernel* kernel;
-  QtViewSolutionGeometry* solution_geometry;
+  QtOutputView* output_view;
   size_t step;
 public:
-  QtViewSolution(Kernel* kernel);
-  ~QtViewSolution();
+  QtOutput(string filename);
+  ~QtOutput();
   void update();
   void update_infos();
 private slots:

src/qt/output/output_view.cpp

@@ -0,0 +1,226 @@
+#include "qt/output/output_view.hpp"
+
+QtOutputView::QtOutputView(QtKernel* k){
+  kernel=k;
+  nX=kernel->geometry.nX;
+  dX=kernel->geometry.dX;
+  lX=kernel->geometry.lX;
+  factor=kernel->factor;
+};
+
+void
+QtOutputView::initializeGL(){
+  glClearColor(1,1,1,1);
+  glEnable(GL_DEPTH_TEST);
+  glEnable(GL_LIGHT0);
+  glEnable(GL_LIGHTING);
+  glColorMaterial(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE);
+  glEnable(GL_COLOR_MATERIAL);
+}
+
+void
+QtOutputView::paintGL(){
+  glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
+  if(step>kernel->step) return;
+  double h=height();
+  double w=width();
+  glBegin(GL_TRIANGLES);
+  for(size_t i=0;i<nX-1;++i){
+    size_t max_left=kernel->nZ(i);
+    size_t max_right=kernel->nZ(i+1);
+    size_t left=0;
+    size_t right=0;
+    while(left<max_left-1 and right<max_right-1){
+      if(kernel->Z(i,left+1)<=kernel->Z(i+1,right+1)){
+        drawTriangle(i,left,i+1,right,i,left+1);
+        ++left;
+      }
+      else{
+        drawTriangle(i,left,i+1,right,i+1,right+1);
+        ++right;
+      }
+    }
+    if(left==max_left-1){
+      while(right<max_right-1){
+        drawTriangle(i,left,i+1,right,i+1,right+1);
+        ++right;
+      }
+    }
+    if(right==max_right-1){
+      while(left<max_left-1){
+        drawTriangle(i,left,i+1,right,i,left+1);
+        ++left;
+      }
+    }
+  }
+  glEnd();
+  //Draw area under hbot
+  glColor3f(0.6,0.6,0.6);
+  glLineWidth(3);
+  glBegin(GL_QUADS);
+  for(size_t ix=0;ix<nX-1;ix+=1){
+    double x=ix*dX;
+    glVertex3f(x,0,1);
+    glVertex3f(x,kernel->geometry.hbot[ix],0);
+    glVertex3f(x+dX,kernel->geometry.hbot[ix+1],0);
+    glVertex3f(x+dX,0,0);
+
+  }
+  glEnd();
+
+  //Draw hsoil
+  glColor3f(0.6,0.3,0);
+  glLineWidth(1);
+  glBegin(GL_LINE_STRIP);
+  for(size_t ix=0;ix<nX;ix+=1){
+    double x=ix*dX;
+    glVertex3f(x,kernel->geometry.hsoil[ix],0.8);
+  }
+  glEnd();
+
+  //Draw hydr
+  glColor3f(1,0,0);
+  glLineWidth(3);
+  glBegin(GL_LINE_STRIP);
+  for(size_t ix=0;ix<nX;ix+=1){
+    double x=ix*dX;
+    glVertex3f(x,kernel->solution[step]->hydr[ix]-Psat/(Physics::rho*Physics::g),1);
+  }
+  glEnd();
+
+  //Draw hov
+  glColor3f(0,0,0);
+  glLineWidth(1);
+  glBegin(GL_LINE_STRIP);
+  for(size_t ix=0;ix<nX;ix+=1){
+    double x=ix*dX;
+    double y=kernel->solution[step]->hov[ix];
+    glVertex3f(x,y,1);
+  }
+  glEnd();
+
+  drawVectors();
+}
+
+void
+QtOutputView::resizeGL(int w,int h){
+  glViewport(0,0,w,h);
+  glMatrixMode(GL_PROJECTION);
+  glLoadIdentity();
+  gluOrtho2D(0,lX,0,factor);
+  glMatrixMode(GL_MODELVIEW);
+  glLoadIdentity();
+}
+
+void
+QtOutputView::drawTriangle(size_t ix1,size_t iz1,size_t ix2,size_t iz2,size_t ix3,size_t iz3){
+  double factor_inv=1;//1/kernel->solution.factor;
+  double h=height();
+  double w=width();
+  double x1=ix1*dX;
+  double x2=ix2*dX;
+  double x3=ix3*dX;
+  double y1=kernel->Z(ix1,iz1)*factor_inv;;
+  double y2=kernel->Z(ix2,iz2)*factor_inv;;
+  double y3=kernel->Z(ix3,iz3)*factor_inv;;
+  setColor(ix1,iz1);
+  glVertex3f(x1,y1,0);
+  setColor(ix2,iz2);
+  glVertex3f(x2,y2,0);
+  setColor(ix3,iz3);
+  glVertex3f(x3,y3,0);
+}
+
+void
+QtOutputView::draw(size_t _step){
+  step=_step;
+  update();
+}
+
+void
+QtOutputView::mousePressEvent(QMouseEvent* event){
+  double x=double(event->x())/width()*lX;
+  double y=1-(double(event->y())/height());
+  displayInfos(x,y*factor);
+}
+
+void
+QtOutputView::setColor(size_t ix,size_t iz){
+  double p=getP(ix,iz);
+  double s=Physics::s(p);
+
+  double r=(1-s);
+  double g=0.7*(1-s);
+  double b=0.8*s+0.4*(1-s);
+  glColor3f(r,g,b);
+}
+
+void
+QtOutputView::drawVector(size_t ix,size_t iz,double u,double v){
+  double x1=ix*dX;
+  double y1=kernel->Z(ix,iz);
+  double x2=x1+u;
+  double y2=y1+v;
+  double x3=x2+(-0.866*u+0.5*v)*scale_arrow_head*lX;
+  double y3=y2+(-0.5*u-0.866*v)*scale_arrow_head*factor;
+  double x4=x2+(-0.866*u-0.5*v)*scale_arrow_head*lX;
+  double y4=y2+(0.5*u-0.866*v)*scale_arrow_head*factor;
+
+  glColor3f(1,1,0);
+  glBegin(GL_LINE_STRIP);
+  glVertex3f(x1,y1,0.8);
+  glVertex3f(x2,y2,0.8);
+  glEnd();
+  glBegin(GL_LINE_STRIP);
+  glVertex3f(x3,y3,0.8);
+  glVertex3f(x2,y2,0.8);
+  glVertex3f(x4,y4,0.8);
+  glEnd();
+
+}
+
+void
+QtOutputView::drawVectors(){
+  for(size_t ix=1;ix<nX-1;ix+=arrow_step_x){
+    //Case iz=0
+    double p1=getP(ix,0);
+    double p2=getP(ix,1);
+    double u=-Physics::k0*Physics::kr(getPl(ix)+Physics::rho*Physics::g*(getl(ix)-kernel->Z(ix,0)))*(getHydr(ix+1)-getHydr(ix-1))/(2*dX);
+    double v=-Physics::k0*Physics::kr(p1)*((p2-p1)/(Physics::rho*Physics::g*kernel->geometry.dZ[ix])+1);
+    drawVector(ix,0,scale_arrow*u,scale_arrow*v);
+    for(size_t iz=arrow_step_z;iz<kernel->nZ(ix)-1;iz+=arrow_step_z){
+      double p0=getP(ix,iz-1);
+      double p1=getP(ix,iz);
+      double p2=getP(ix,iz+1);
+      double u=-Physics::k0*Physics::kr(getPl(ix)+Physics::rho*Physics::g*(getl(ix)-kernel->Z(ix,iz)))*(getHydr(ix+1)-getHydr(ix-1))/(2*dX);
+      double v=-Physics::k0*Physics::kr(p1)*((p2-p0)/(2*Physics::rho*Physics::g*kernel->geometry.dZ[ix])+1);
+      drawVector(ix,iz,scale_arrow*u,scale_arrow*v);
+    }
+  }
+}
+
+void
+QtOutputView::displayInfos(double x,double z){
+  cout<<"---------------------------------"<<endl;
+  cout<<"x = "<<x<<endl;
+  cout<<"z = "<<z<<endl;
+  size_t ix=floor(x/dX+0.5);
+  if(ix==nX) --ix;
+  cout<<"ix = "<<ix<<endl;
+
+  double hbot=kernel->geometry.hbot[ix];
+  cout<<"hbot = "<<hbot<<endl;
+  double hsoil=kernel->geometry.hsoil[ix];
+  cout<<"hsoil = "<<hsoil<<endl;
+  cout<<"l = "<<getl(ix)<<endl;
+    if(z>hbot and z<hsoil){
+    double dZ=kernel->geometry.dZ[ix];
+    cout<<"dZ = "<<dZ<<endl;
+    size_t nz=kernel->geometry.nZ[ix];
+    cout<<"nZ = "<<nz<<endl;
+    size_t iz=floor((z-hbot)/dZ+0.5);
+    if(iz==nz) --iz;
+    cout<<"iz = "<<iz<<endl;
+    cout<<"P = "<<getP(ix,iz)<<endl;
+    }
+}